ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.35
Committed: Sat Oct 13 21:09:57 2012 UTC (12 years, 7 months ago) by mingyang
Content type: text/plain
Branch: MAIN
Changes since 1.34: +5 -5 lines
Log Message:
30/120 --> 30./120.

File Contents

# Content
1 #include "MitPhysics/Mods/interface/PhotonTreeWriter.h"
2 #include "MitAna/DataTree/interface/PhotonCol.h"
3 #include "MitAna/DataTree/interface/PFCandidateCol.h"
4 #include "MitAna/DataTree/interface/StableData.h"
5 #include "MitAna/DataTree/interface/StableParticle.h"
6 #include "MitAna/DataTree/interface/PFMet.h"
7 #include "MitPhysics/Init/interface/ModNames.h"
8 #include "MitPhysics/Utils/interface/IsolationTools.h"
9 #include "MitPhysics/Utils/interface/PhotonTools.h"
10 #include "MitPhysics/Utils/interface/VertexTools.h"
11 #include "MitPhysics/Utils/interface/PFMetCorrectionTools.h"
12 #include "MitAna/DataTree/interface/PFJetCol.h"
13 #include "MitAna/DataTree/interface/GenJetCol.h"
14 #include "TDataMember.h"
15 #include "TFile.h"
16 #include <TNtuple.h>
17 #include <TRandom3.h>
18 #include <TSystem.h>
19
20 using namespace mithep;
21
22 ClassImp(mithep::PhotonTreeWriter)
23 templateClassImp(mithep::PhotonTreeWriterPhoton)//ming: what's this?
24 ClassImp(mithep::PhotonTreeWriterDiphotonEvent)
25
26 //--------------------------------------------------------------------------------------------------
27 PhotonTreeWriter::PhotonTreeWriter(const char *name, const char *title) :
28 // Base Module...
29 BaseMod (name,title),
30 // define all the Branches to load
31 fPhotonBranchName (Names::gkPhotonBrn),
32 fPFPhotonName ("PFPhotons"),
33 fElectronName (Names::gkElectronBrn),
34 fGoodElectronName (Names::gkElectronBrn),
35 fConversionName (Names::gkMvfConversionBrn),
36 fTrackBranchName (Names::gkTrackBrn),
37 fPileUpDenName (Names::gkPileupEnergyDensityBrn),
38 fPVName (Names::gkPVBeamSpotBrn),
39 fBeamspotName (Names::gkBeamSpotBrn),
40 fPFCandName (Names::gkPFCandidatesBrn),
41 fPFNoPileUpName ("PFNoPileUp"),
42 fPFPileUpName ("PFPileUp"),
43 fMCParticleName (Names::gkMCPartBrn),
44 fMCEventInfoName (Names::gkMCEvtInfoBrn),
45 fPileUpName (Names::gkPileupInfoBrn),
46 fSuperClusterName ("PFSuperClusters"),
47 fPFMetName ("PFMet"),
48 fPFJetName (Names::gkPFJetBrn),
49 funcorrPFJetName ("AKt5PFJets"),
50 fGenJetName ("AKT5GenJets"),
51 fLeptonTagElectronsName ("HggLeptonTagElectrons"),
52 fLeptonTagMuonsName ("HggLeptonTagMuons"),
53
54 fIsData (false),
55 fPhotonsFromBranch (kTRUE),
56 fPVFromBranch (kTRUE),
57 fGoodElectronsFromBranch(kTRUE),
58 fPFJetsFromBranch (kTRUE),
59 // ----------------------------------------
60 // collections....
61 fPhotons (0),
62 fPFPhotons (0),
63 fElectrons (0),
64 fConversions (0),
65 fTracks (0),
66 fPileUpDen (0),
67 fPV (0),
68 fBeamspot (0),
69 fPFCands (0),
70 fMCParticles (0),
71 fMCEventInfo (0),
72 fPileUp (0),
73 fSuperClusters (0),
74 fPFJets (0),
75 fGenJets (0),
76 funcorrPFJets (0),
77
78 fLeptonTagElectrons (0),
79 fLeptonTagMuons (0),
80 fPFNoPileUpCands (0),
81 fPFPileUpCands (0),
82
83 fLoopOnGoodElectrons (kFALSE),
84 fApplyElectronVeto (kTRUE),
85 fWriteDiphotonTree (kTRUE),
86 fWriteSingleTree (kTRUE),
87 fEnablePFPhotons (kTRUE),
88 fExcludeSinglePrompt (kFALSE),
89 fExcludeDoublePrompt (kFALSE),
90 fEnableJets (kFALSE),
91 fApplyJetId (kFALSE),
92 fApplyLeptonTag (kFALSE),
93 fApplyVBFTag (kFALSE),
94 fApplyBTag (kFALSE),
95 fApplyPFMetCorrections (kFALSE),
96 fFillClusterArrays (kFALSE),
97 fFillVertexTree (kFALSE),
98 fDo2012LepTag (kFALSE),
99 fPhFixDataFile (gSystem->Getenv("CMSSW_BASE") +
100 TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat")),
101 fBeamspotWidth (5.8),
102 fTupleName ("hPhotonTree")
103 {
104 // Constructor
105 }
106
107 PhotonTreeWriter::~PhotonTreeWriter()
108 {
109 // Destructor
110 }
111
112 //--------------------------------------------------------------------------------------------------
113 void PhotonTreeWriter::Process()
114 {
115 // ------------------------------------------------------------
116 // Process entries of the tree.
117 LoadEventObject(fPhotonBranchName, fPhotons);
118 LoadEventObject(fGoodElectronName, fGoodElectrons);
119
120 // lepton tag collections
121 if( fApplyLeptonTag ) {
122 LoadEventObject(fLeptonTagElectronsName, fLeptonTagElectrons);
123 LoadEventObject(fLeptonTagMuonsName, fLeptonTagMuons);
124 }
125
126 const BaseCollection *egcol = 0;
127 if (fLoopOnGoodElectrons)
128 egcol = fGoodElectrons;
129 else
130 egcol = fPhotons;
131 if (egcol->GetEntries()<1)
132 return;
133
134 if (fEnablePFPhotons) LoadEventObject(fPFPhotonName, fPFPhotons);
135 LoadEventObject(fElectronName, fElectrons);
136 LoadEventObject(fConversionName, fConversions);
137 LoadEventObject(fTrackBranchName, fTracks);
138 LoadEventObject(fPileUpDenName, fPileUpDen);
139 LoadEventObject(fPVName, fPV);
140 LoadEventObject(fBeamspotName, fBeamspot);
141 LoadEventObject(fPFCandName, fPFCands);
142 LoadEventObject(fSuperClusterName, fSuperClusters);
143 LoadEventObject(fPFMetName, fPFMet);
144 // LoadEventObject(fPFNoPileUpName, fPFNoPileUpCands);
145 // LoadEventObject(fPFPileUpName, fPFPileUpCands);
146
147 if (fEnableJets){
148 LoadEventObject(fPFJetName, fPFJets);
149 //LoadEventObject(funcorrPFJetName, funcorrPFJets);
150 LoadBranch(funcorrPFJetName);
151 // if(!fIsData) LoadEventObject(fGenJetName, fGenJets);
152 }
153
154 // ------------------------------------------------------------
155 // load event based information
156 Int_t _numPU = -1.; // some sensible default values....
157 Int_t _numPUminus = -1.; // some sensible default values....
158 Int_t _numPUplus = -1.; // some sensible default values....
159
160 Float_t rho = -99.;
161 if( fPileUpDen->GetEntries() > 0 )
162 rho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
163
164 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));//ming:what's this?
165
166 if( !fIsData ) {
167 LoadBranch(fMCParticleName);
168 LoadBranch(fMCEventInfoName);
169 LoadBranch(fPileUpName);
170 if (fEnableJets) LoadEventObject(fGenJetName, fGenJets);
171 } else fGenJets = NULL;
172
173 if( !fIsData ) {
174 for (UInt_t i=0; i<fPileUp->GetEntries(); ++i) {
175 const PileupInfo *puinfo = fPileUp->At(i);
176 if (puinfo->GetBunchCrossing()==0) _numPU = puinfo->GetPU_NumInteractions();
177 else if (puinfo->GetBunchCrossing() == -1) _numPUminus = puinfo->GetPU_NumInteractions();
178 else if (puinfo->GetBunchCrossing() == 1) _numPUplus = puinfo->GetPU_NumInteractions();
179 }
180 }
181
182 // in case of a MC event, try to find Higgs and Higgs decay Z poisition
183 Float_t _pth = -100.;
184 Float_t _decayZ = -100.;
185 Float_t _genmass = -100.;
186 if (!fIsData)
187 FindHiggsPtAndZ(_pth, _decayZ, _genmass);
188
189
190 fDiphotonEvent->genz = -999.;
191 if (!fIsData) {
192 fDiphotonEvent->genz = fMCParticles->At(0)->DecayVertex().Z();
193 }
194
195 fDiphotonEvent->mcprocid = -999;
196 if (!fIsData) {
197 fDiphotonEvent->mcprocid = fMCEventInfo->ProcessId();
198 }
199
200 Double_t _evt = GetEventHeader()->EvtNum();
201
202 Double_t _spfMet = fPFMet->At(0)->SumEt();
203
204 fDiphotonEvent->leptonTag = -1; // disabled
205
206 fDiphotonEvent->rho = fPileUpDen->At(0)->RhoKt6PFJets();
207 fDiphotonEvent->rho25 = fPileUpDen->At(0)->RhoRandomLowEta();
208 fDiphotonEvent->rhoold = fPileUpDen->At(0)->Rho();
209 fDiphotonEvent->genHiggspt = _pth;
210 fDiphotonEvent->genHiggsZ = _decayZ;
211 fDiphotonEvent->genmass = _genmass;
212 fDiphotonEvent->gencostheta = -99.;
213 fDiphotonEvent->nVtx = fPV->GetEntries();
214 fDiphotonEvent->bsX = fBeamspot->At(0)->X();
215 fDiphotonEvent->bsY = fBeamspot->At(0)->Y();
216 fDiphotonEvent->bsZ = fBeamspot->At(0)->Z();
217 fDiphotonEvent->bsSigmaZ = fBeamspot->At(0)->SigmaZ();
218 fDiphotonEvent->vtxX = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->X() : -99.;
219 fDiphotonEvent->vtxY = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Y() : -99.;
220 fDiphotonEvent->vtxZ = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Z() : -99.;
221 fDiphotonEvent->numPU = _numPU;
222 fDiphotonEvent->numPUminus = _numPUminus;
223 fDiphotonEvent->numPUplus = _numPUplus;
224 fDiphotonEvent->mass = -99.;
225 fDiphotonEvent->ptgg = -99.;
226 fDiphotonEvent->costheta = -99.;
227 fDiphotonEvent->mt = -99.;
228 fDiphotonEvent->cosphimet = -99.;
229 fDiphotonEvent->mtele = -99.;
230 fDiphotonEvent->cosphimetele = -99.;
231 fDiphotonEvent->evt = GetEventHeader()->EvtNum();
232 fDiphotonEvent->run = GetEventHeader()->RunNum();
233 fDiphotonEvent->lumi = GetEventHeader()->LumiSec();
234 fDiphotonEvent->evtcat = -99;
235 fDiphotonEvent->nobj = fPhotons->GetEntries();
236 fDiphotonEvent->pfmet = fPFMet->At(0)->Pt();
237 fDiphotonEvent->pfmetphi = fPFMet->At(0)->Phi();
238 fDiphotonEvent->pfmetx = fPFMet->At(0)->Px();
239 fDiphotonEvent->pfmety = fPFMet->At(0)->Py();
240 fDiphotonEvent->spfMet = _spfMet;
241 fDiphotonEvent->masscor = -99.;
242 fDiphotonEvent->masscorerr = -99.;
243 fDiphotonEvent->masscorele = -99.;
244 fDiphotonEvent->masscoreleerr = -99.;
245 fDiphotonEvent->ismc = GetEventHeader()->IsMC();
246
247 //jets
248 const Jet *jet1 = 0;
249 const Jet *jet2 = 0;
250 const Jet *jetcentral = 0;
251
252 fDiphotonEvent->jet1pt = -99.;
253 fDiphotonEvent->jet1eta = -99.;
254 fDiphotonEvent->jet1phi = -99.;
255 fDiphotonEvent->jet1mass = -99.;
256 fDiphotonEvent->jet2pt = -99.;
257 fDiphotonEvent->jet2eta = -99.;
258 fDiphotonEvent->jet2phi = -99.;
259 fDiphotonEvent->jet2mass = -99.;
260 fDiphotonEvent->jetcentralpt = -99.;
261 fDiphotonEvent->jetcentraleta = -99.;
262 fDiphotonEvent->jetcentralphi = -99.;
263 fDiphotonEvent->jetcentralmass = -99.;
264 fDiphotonEvent->dijetpt = -99.;
265 fDiphotonEvent->dijeteta = -99.;
266 fDiphotonEvent->dijetphi = -99.;
267 fDiphotonEvent->dijetmass = -99.;
268 fDiphotonEvent->jetetaplus = -99.;
269 fDiphotonEvent->jetetaminus = -99.;
270 fDiphotonEvent->zeppenfeld = -99.;
271 fDiphotonEvent->dphidijetgg = -99.;
272
273 //uncorrected jets
274 // const Jet *uncorrjet1 = 0;
275 // const Jet *uncorrjet2 = 0;
276 // const Jet *uncorrjetcentral = 0;
277
278 // fDiphotonEvent->uncorrjet1pt = -99.;
279 // fDiphotonEvent->uncorrjet1eta = -99.;
280 // fDiphotonEvent->uncorrjet1phi = -99.;
281 // fDiphotonEvent->uncorrjet1mass = -99.;
282 // fDiphotonEvent->uncorrjet2pt = -99.;
283 // fDiphotonEvent->uncorrjet2eta = -99.;
284 // fDiphotonEvent->uncorrjet2phi = -99.;
285 // fDiphotonEvent->uncorrjet2mass = -99.;
286 // fDiphotonEvent->uncorrjetcentralpt = -99.;
287 // fDiphotonEvent->uncorrjetcentraleta = -99.;
288 // fDiphotonEvent->uncorrjetcentralphi = -99.;
289 // fDiphotonEvent->uncorrjetcentralmass = -99.;
290 // fDiphotonEvent->diuncorrjetpt = -99.;
291 // fDiphotonEvent->diuncorrjeteta = -99.;
292 // fDiphotonEvent->diuncorrjetphi = -99.;
293 // fDiphotonEvent->diuncorrjetmass = -99.;
294
295
296
297 Int_t nhitsbeforevtxmax = 1;
298 if (!fApplyElectronVeto)
299 nhitsbeforevtxmax = 999;
300
301 if (egcol->GetEntries()>=2) {
302 const Particle *p1 = 0;
303 const Particle *p2 = 0;
304 const Photon *phHard = 0;
305 const Photon *phSoft = 0;
306 const Electron *ele1 = 0;
307 const Electron *ele2 = 0;
308 const SuperCluster *sc1 = 0;
309 const SuperCluster *sc2 = 0;
310
311 if (fLoopOnGoodElectrons) {
312 ele1 = fGoodElectrons->At(0);
313 ele2 = fGoodElectrons->At(1);
314 p1 = ele1;
315 p2 = ele2;
316 sc1 = ele1->SCluster();
317 sc2 = ele2->SCluster();
318 phHard = PhotonTools::MatchedPhoton(ele1,fPhotons);
319 phSoft = PhotonTools::MatchedPhoton(ele2,fPhotons);
320 }
321 else {
322 phHard = fPhotons->At(0);
323 phSoft = fPhotons->At(1);
324 p1 = phHard;
325 p2 = phSoft;
326 sc1 = phHard->SCluster();
327 sc2 = phSoft->SCluster();
328 ele1 = PhotonTools::MatchedElectron(phHard,fGoodElectrons);
329 ele2 = PhotonTools::MatchedElectron(phSoft,fGoodElectrons);
330 }
331
332 const DecayParticle *conv1 = PhotonTools::MatchedConversion(sc1,fConversions,bsp,
333 nhitsbeforevtxmax);
334 const DecayParticle *conv2 = PhotonTools::MatchedConversion(sc2,fConversions,bsp,
335 nhitsbeforevtxmax);
336
337 const SuperCluster *pfsc1 = PhotonTools::MatchedPFSC(sc1,fPFPhotons, fElectrons);
338 const SuperCluster *pfsc2 = PhotonTools::MatchedPFSC(sc2,fPFPhotons, fElectrons);
339
340 const MCParticle *phgen1 = 0;
341 const MCParticle *phgen2 = 0;
342 if (!fIsData) {
343 phgen1 = PhotonTools::MatchMC(p1,fMCParticles,!fApplyElectronVeto);
344 phgen2 = PhotonTools::MatchMC(p2,fMCParticles,!fApplyElectronVeto);
345 }
346
347 if (fExcludeSinglePrompt && (phgen1 || phgen2))
348 return;
349 if (fExcludeDoublePrompt && (phgen1 && phgen2))
350 return;
351
352 if (!fLoopOnGoodElectrons && phHard->HasPV()) {
353 fDiphotonEvent->vtxX = phHard->PV()->X();
354 fDiphotonEvent->vtxY = phHard->PV()->Y();
355 fDiphotonEvent->vtxZ = phHard->PV()->Z();
356 fDiphotonEvent->vtxprob = phHard->VtxProb();
357 }
358
359 // fill Btag information... set to true if One jet fullfills
360 fDiphotonEvent -> btagJet1 = -1.;
361 fDiphotonEvent -> btagJet1Pt = -99.;
362 fDiphotonEvent -> btagJet1Eta = -99.;
363
364 fDiphotonEvent -> btagJet2 = -1.;
365 fDiphotonEvent -> btagJet2Pt = -99.;
366 fDiphotonEvent -> btagJet2Eta = -99.;
367
368 if( fApplyBTag && fEnableJets ) {
369 float highTag = 0.;
370 float highJetPt = 0.;
371 float highJetEta = 0.;
372
373 float trailTag = 0.;
374 float trailJetPt = 0.;
375 float trailJetEta = 0.;
376
377 for (UInt_t ijet=0; ijet<fPFJets->GetEntries();++ijet) {
378 const Jet *jet = fPFJets->At(ijet);
379 if ( jet->Pt() < 20. || jet->AbsEta() > 2.4 ) continue;
380 if ( jet->CombinedSecondaryVertexBJetTagsDisc() > highTag ) {
381
382 trailTag = highTag;
383 trailJetPt = highJetPt;
384 trailJetEta = highJetEta;
385
386 highTag = jet->CombinedSecondaryVertexBJetTagsDisc();
387 highJetPt = jet->Pt();
388 highJetEta = jet->Eta();
389
390 } else if ( jet->CombinedSecondaryVertexBJetTagsDisc() > trailTag ) {
391
392 trailTag = jet->CombinedSecondaryVertexBJetTagsDisc();
393 trailJetPt = jet->Pt();
394 trailJetEta = jet->Eta();
395 }
396 }
397 fDiphotonEvent -> btagJet1 = highTag;
398 fDiphotonEvent -> btagJet1Pt = highJetPt;
399 fDiphotonEvent -> btagJet1Eta = highJetEta;
400
401 fDiphotonEvent -> btagJet2 = trailTag;
402 fDiphotonEvent -> btagJet2Pt = trailJetPt;
403 fDiphotonEvent -> btagJet2Eta = trailJetEta;
404
405 }
406
407 //fill jet variables
408 const Vertex *selvtx = fPV->At(0);
409 if (!fLoopOnGoodElectrons && phHard->HasPV()) selvtx = phHard->PV();
410 if (fEnableJets) {
411 for (UInt_t ijet=0; ijet<fPFJets->GetEntries();++ijet) {
412 const Jet *jet = fPFJets->At(ijet);
413 if (jet->AbsEta()<4.7 && MathUtils::DeltaR(jet,p1)>0.5 && MathUtils::DeltaR(jet,p2)>0.5) {
414 const PFJet *pfjet = dynamic_cast<const PFJet*>(jet);
415 if (!pfjet) continue;
416 if (fApplyJetId && !fJetId.passCut(pfjet,selvtx,fPV)) continue;
417 if (!jet1) jet1 = jet;
418 else if (!jet2) jet2 = jet;
419 else if (!jetcentral && 0) jetcentral = jet;
420 }
421 if (jet1&&jet2&&jetcentral) break;
422 }
423 }
424
425 if (jet1) {
426 fDiphotonEvent->jet1pt = jet1->Pt();
427 fDiphotonEvent->jet1eta = jet1->Eta();
428 fDiphotonEvent->jet1phi = jet1->Phi();
429 fDiphotonEvent->jet1mass = jet1->Mass();
430 }
431
432 if (jet2) {
433 fDiphotonEvent->jet2pt = jet2->Pt();
434 fDiphotonEvent->jet2eta = jet2->Eta();
435 fDiphotonEvent->jet2phi = jet2->Phi();
436 fDiphotonEvent->jet2mass = jet2->Mass();
437 }
438
439 if (jetcentral) {
440 fDiphotonEvent->jetcentralpt = jetcentral->Pt();
441 fDiphotonEvent->jetcentraleta = jetcentral->Eta();
442 fDiphotonEvent->jetcentralphi = jetcentral->Phi();
443 fDiphotonEvent->jetcentralmass = jetcentral->Mass();
444 }
445
446 if (jet1&&jet2){
447 FourVectorM momjj = (jet1->Mom() + jet2->Mom());
448
449 fDiphotonEvent->dijetpt = momjj.Pt();
450 fDiphotonEvent->dijeteta = momjj.Eta();
451 fDiphotonEvent->dijetphi = momjj.Phi();
452 fDiphotonEvent->dijetmass = momjj.M();
453
454 if (jet1->Eta()>jet2->Eta()) {
455 fDiphotonEvent->jetetaplus = jet1->Eta();
456 fDiphotonEvent->jetetaminus = jet2->Eta();
457 }
458 else {
459 fDiphotonEvent->jetetaplus = jet2->Eta();
460 fDiphotonEvent->jetetaminus = jet1->Eta();
461 }
462 }
463
464 //added gen. info of whether a lep. or nutrino is from W or Z --Heng 02/14/2012 12:30 EST
465 Double_t _fromZ = -99;
466 Double_t _fromW = -99;
467 Float_t _zpt = -99;
468 Float_t _zEta = -99;
469 Float_t _allZpt = -99;
470 Float_t _allZEta = -99;
471
472 if( !fIsData ){
473
474 // loop over all GEN particles and look for nutrinos whoes mother is Z
475 for(UInt_t j=0; j<fMCParticles->GetEntries(); ++j) {
476 const MCParticle* p = fMCParticles->At(j);
477 if( p->AbsPdgId()==23 ||p->AbsPdgId()==32 || p->AbsPdgId()==33 ) {
478 _allZpt=p->Pt();
479 _allZEta=p->Eta();
480 if (p->HasDaughter(12,kFALSE) || p->HasDaughter(14,kFALSE) || p->HasDaughter(16,kFALSE) ||p->HasDaughter(18,kFALSE) ) {
481 _fromZ=1;
482 _zpt=p->Pt();
483 _zEta=p->Eta();
484 }
485 }
486 else _fromW=1;
487 }
488 }
489
490 fDiphotonEvent->fromZ = _fromZ;
491 fDiphotonEvent->fromW = _fromW;
492 fDiphotonEvent->zpt = _zpt;
493 fDiphotonEvent->zEta = _zEta;
494 fDiphotonEvent->allZpt = _allZpt;
495 fDiphotonEvent->allZEta = _allZEta;
496
497 Double_t _dphiMetgg = -99;
498 Double_t _cosdphiMetgg = -99;
499 Double_t _dphiPhPh = -99;
500
501 Double_t _mass = -99.;
502 Double_t _masserr = -99.;
503 Double_t _masserrsmeared = -99.;
504 Double_t _masserrwrongvtx = -99.;
505 Double_t _masserrsmearedwrongvtx = -99.;
506 Double_t _ptgg = -99.;
507 Double_t _etagg = -99.;
508 Double_t _phigg = -99.;
509 Double_t _costheta = -99.;
510 PhotonTools::DiphotonR9EtaPtCats _evtcat = PhotonTools::kOctCat0;
511 if (phHard && phSoft) {
512 _dphiMetgg = MathUtils::DeltaPhi((phHard->Mom()+phSoft->Mom()).Phi(),fPFMet->At(0)->Phi());
513 _cosdphiMetgg = TMath::Cos(_dphiMetgg);
514 _dphiPhPh = MathUtils::DeltaPhi((phHard->Mom()).Phi(),(phSoft->Mom()).Phi());
515 _mass = (phHard->Mom()+phSoft->Mom()).M();
516 _masserr = 0.5*_mass*TMath::Sqrt(phHard->EnergyErr()*phHard->EnergyErr()/phHard->E()/phHard->E() + phSoft->EnergyErr()*phSoft->EnergyErr()/phSoft->E()/phSoft->E());
517 _masserrsmeared = 0.5*_mass*TMath::Sqrt(phHard->EnergyErrSmeared()*phHard->EnergyErrSmeared()/phHard->E()/phHard->E() + phSoft->EnergyErrSmeared()*phSoft->EnergyErrSmeared()/phSoft->E()/phSoft->E());
518 _ptgg = (phHard->Mom()+phSoft->Mom()).Pt();
519 _etagg = (phHard->Mom()+phSoft->Mom()).Eta();
520 _phigg = (phHard->Mom()+phSoft->Mom()).Phi();
521 _costheta = ThreeVector(phHard->Mom()).Unit().Dot(ThreeVector(phSoft->Mom()).Unit());
522 _evtcat = PhotonTools::DiphotonR9EtaPtCat(phHard,phSoft);
523
524 const Double_t dz = sqrt(2.0)*fBeamspotWidth;
525 Double_t deltamvtx = _mass*VertexTools::DeltaMassVtx(phHard->CaloPos().X(),
526 phHard->CaloPos().Y(),
527 phHard->CaloPos().Z(),
528 phSoft->CaloPos().X(),
529 phSoft->CaloPos().Y(),
530 phSoft->CaloPos().Z(),
531 fDiphotonEvent->vtxX,
532 fDiphotonEvent->vtxY,
533 fDiphotonEvent->vtxZ,
534 dz);
535 fDiphotonEvent->deltamvtx = deltamvtx;
536
537 _masserrwrongvtx = TMath::Sqrt(_masserr*_masserr + deltamvtx*deltamvtx);
538 _masserrsmearedwrongvtx = TMath::Sqrt(_masserrsmeared*_masserrsmeared + deltamvtx*deltamvtx);
539
540 if (jet1 && jet2) {
541 fDiphotonEvent->zeppenfeld = TMath::Abs(_etagg - 0.5*(jet1->Eta()+jet2->Eta()));
542 fDiphotonEvent->dphidijetgg = MathUtils::DeltaPhi( (jet1->Mom()+jet2->Mom()).Phi(), _phigg );
543 }
544
545 PFJetOArr pfjets;
546 for (UInt_t ijet=0; ijet<fPFJets->GetEntries(); ++ijet) {
547 const PFJet *pfjet = dynamic_cast<const PFJet*>(fPFJets->At(ijet));
548 if (pfjet && MathUtils::DeltaR(*pfjet,*phHard)>0.3 && MathUtils::DeltaR(*pfjet,*phSoft)>0.3) pfjets.Add(pfjet);
549 }
550
551 PFCandidateOArr pfcands;
552 for (UInt_t icand=0; icand<fPFCands->GetEntries(); ++icand) {
553 const PFCandidate *pfcand = fPFCands->At(icand);
554 if (MathUtils::DeltaR(*pfcand,*phHard)>0.1 && MathUtils::DeltaR(*pfcand,*phSoft)>0.1) pfcands.Add(pfcand);
555 }
556
557 const Vertex *firstvtx = fPV->At(0);
558 const Vertex *selvtx = fPV->At(0);
559
560 if (!fLoopOnGoodElectrons && phHard->HasPV()) {
561 selvtx = phHard->PV();
562 }
563
564 if (0) //disable for now for performance reasons
565 {
566 Met mmet = fMVAMet.GetMet( false,
567 0.,0.,0.,
568 0.,0.,0.,
569 fPFMet->At(0),
570 &pfcands,selvtx,fPV, fPileUpDen->At(0)->Rho(),
571 &pfjets,
572 int(fPV->GetEntries()),
573 kFALSE);
574
575 TMatrixD *metcov = fMVAMet.GetMetCovariance();
576
577 ThreeVector fullmet(mmet.Px() - phHard->Px() - phSoft->Px(),
578 mmet.Py() - phHard->Py() - phSoft->Py(),
579 0.);
580
581 ROOT::Math::SMatrix<double,2,2,ROOT::Math::MatRepSym<double,2> > mcov;
582 mcov(0,0) = (*metcov)(0,0);
583 mcov(0,1) = (*metcov)(0,1);
584 mcov(1,0) = (*metcov)(1,0);
585 mcov(1,1) = (*metcov)(1,1);
586 ROOT::Math::SVector<double,2> vmet;
587 vmet(0) = fullmet.X();
588 vmet(1) = fullmet.Y();
589 mcov.Invert();
590 Double_t metsig = sqrt(ROOT::Math::Similarity(mcov,vmet));
591
592 fDiphotonEvent->mvametsel = fullmet.Rho();
593 fDiphotonEvent->mvametselphi = fullmet.Phi();
594 fDiphotonEvent->mvametselx = fullmet.X();
595 fDiphotonEvent->mvametsely = fullmet.Y();
596 fDiphotonEvent->mvametselsig = metsig;
597 }
598
599 if (0) //disable for now for performance reasons
600 {
601 Met mmet = fMVAMet.GetMet( false,
602 0.,0.,0.,
603 0.,0.,0.,
604 fPFMet->At(0),
605 &pfcands,firstvtx,fPV, fPileUpDen->At(0)->Rho(),
606 &pfjets,
607 int(fPV->GetEntries()),
608 kFALSE);
609
610 TMatrixD *metcov = fMVAMet.GetMetCovariance();
611
612 ThreeVector fullmet(mmet.Px() - phHard->Px() - phSoft->Px(),
613 mmet.Py() - phHard->Py() - phSoft->Py(),
614 0.);
615
616 ROOT::Math::SMatrix<double,2,2,ROOT::Math::MatRepSym<double,2> > mcov;
617 mcov(0,0) = (*metcov)(0,0);
618 mcov(0,1) = (*metcov)(0,1);
619 mcov(1,0) = (*metcov)(1,0);
620 mcov(1,1) = (*metcov)(1,1);
621 ROOT::Math::SVector<double,2> vmet;
622 vmet(0) = fullmet.X();
623 vmet(1) = fullmet.Y();
624 mcov.Invert();
625 Double_t metsig = sqrt(ROOT::Math::Similarity(mcov,vmet));
626
627 fDiphotonEvent->mvametfirst = fullmet.Rho();
628 fDiphotonEvent->mvametfirstphi = fullmet.Phi();
629 fDiphotonEvent->mvametfirstx = fullmet.X();
630 fDiphotonEvent->mvametfirsty = fullmet.Y();
631 fDiphotonEvent->mvametfirstsig = metsig;
632 }
633
634 }
635
636
637 fDiphotonEvent->corrpfmet = -99.;
638 fDiphotonEvent->corrpfmetphi = -99.;
639 fDiphotonEvent->corrpfmetx = -99.;
640 fDiphotonEvent->corrpfmety = -99.;
641
642 Met *corrMet =NULL;
643
644 if (fApplyPFMetCorrections){
645 corrMet = new Met(fPFMet->At(0)->Px(),fPFMet->At(0)->Py());
646
647 if (!fIsData){
648 PFMetCorrectionTools::correctMet(corrMet,phHard,phSoft,1,0,funcorrPFJets,fGenJets,fPFJets,_evt);
649 PFMetCorrectionTools::shiftMet(corrMet,fIsData,_spfMet);
650 }
651 else {
652 PFMetCorrectionTools::shiftMet(corrMet,fIsData,_spfMet);
653 PFMetCorrectionTools::correctMet(corrMet,phHard,phSoft,0,1,funcorrPFJets,fGenJets,fPFJets,_evt);
654 }
655
656 fDiphotonEvent->corrpfmet = corrMet->Pt();
657 fDiphotonEvent->corrpfmetphi = corrMet->Phi();
658 fDiphotonEvent->corrpfmetx = corrMet->Px();
659 fDiphotonEvent->corrpfmety = corrMet->Py();
660
661 delete corrMet;
662 }
663
664
665 Float_t _massele = -99.;
666 Float_t _ptee = -99.;
667 Float_t _costhetaele = -99.;
668 if (ele1 && ele2) {
669 _massele = (ele1->Mom()+ele2->Mom()).M();
670 _ptee = (ele1->Mom()+ele2->Mom()).Pt();
671 _costhetaele = ThreeVector(ele1->Mom()).Unit().Dot(ThreeVector(ele2->Mom()).Unit());
672 }
673
674 Float_t _gencostheta = -99.;
675 if (phgen1 && phgen2) {
676 _gencostheta = ThreeVector(phgen1->Mom()).Unit().Dot(ThreeVector(phgen2->Mom()).Unit());
677 }
678
679 fDiphotonEvent->gencostheta = _gencostheta;
680 fDiphotonEvent->dphiMetgg = _dphiMetgg;
681 fDiphotonEvent->cosdphiMetgg = _cosdphiMetgg;
682 fDiphotonEvent->dphiPhPh = _dphiPhPh;
683 fDiphotonEvent->mass = _mass;
684 fDiphotonEvent->masserr = _masserr;
685 fDiphotonEvent->masserrsmeared = _masserrsmeared;
686 fDiphotonEvent->masserrwrongvtx = _masserrwrongvtx;
687 fDiphotonEvent->masserrsmearedwrongvtx = _masserrsmearedwrongvtx;
688 fDiphotonEvent->ptgg = _ptgg;
689 fDiphotonEvent->etagg = _etagg;
690 fDiphotonEvent->phigg = _phigg;
691 fDiphotonEvent->costheta = _costheta;;
692 fDiphotonEvent->massele = _massele;
693 fDiphotonEvent->ptee = _ptee;
694 fDiphotonEvent->costhetaele = _costhetaele;
695 fDiphotonEvent->evtcat = _evtcat;
696
697 fDiphotonEvent->photons[0].SetVars(phHard,conv1,ele1,pfsc1,phgen1,fPhfixph,fPhfixele,fTracks,fPV,fPFCands,rho,fFillClusterArrays,fElectrons,fConversions,bsp,fApplyElectronVeto);
698 fDiphotonEvent->photons[1].SetVars(phSoft,conv2,ele2,pfsc2,phgen2,fPhfixph,fPhfixele,fTracks,fPV,fPFCands,rho,fFillClusterArrays,fElectrons,fConversions,bsp,fApplyElectronVeto);
699
700 Float_t ph1ecor = fDiphotonEvent->photons[0].Ecor();
701 Float_t ph1ecorerr = fDiphotonEvent->photons[0].Ecorerr();
702 Float_t ph2ecor = fDiphotonEvent->photons[1].Ecor();
703 Float_t ph2ecorerr = fDiphotonEvent->photons[1].Ecorerr();
704
705 Float_t ph1ecorele = fDiphotonEvent->photons[0].Ecorele();
706 Float_t ph1ecoreleerr = fDiphotonEvent->photons[0].Ecoreleerr();
707 Float_t ph2ecorele = fDiphotonEvent->photons[1].Ecorele();
708 Float_t ph2ecoreleerr = fDiphotonEvent->photons[1].Ecoreleerr();
709
710
711 fDiphotonEvent->masscor = TMath::Sqrt(2.0*ph1ecor*ph2ecor*(1.0-fDiphotonEvent->costheta));
712 fDiphotonEvent->masscorerr = 0.5*fDiphotonEvent->masscor*
713 TMath::Sqrt(ph1ecorerr*ph1ecorerr/ph1ecor/ph1ecor + ph2ecorerr*ph2ecorerr/ph2ecor/ph2ecor);
714
715 fDiphotonEvent->masscorele = TMath::Sqrt(2.0*ph1ecorele*ph2ecorele*
716 (1.0-fDiphotonEvent->costheta));
717 fDiphotonEvent->masscoreleerr = 0.5*fDiphotonEvent->masscorele*
718 TMath::Sqrt(ph1ecoreleerr*ph1ecoreleerr/ph1ecorele/ph1ecorele +
719 ph2ecoreleerr*ph2ecoreleerr/ph2ecorele/ph2ecorele);
720
721 //printf("r9 = %5f, photon sigieie = %5f, seed sigieie = %5f\n",phHard->R9(),
722 // phHard->CoviEtaiEta(),sqrt(phHard->SCluster()->Seed()->CoviEtaiEta()));
723
724 // MuonStuff
725 fDiphotonEvent-> muonPt = -99.;
726 fDiphotonEvent-> muonEta = -99.;
727 fDiphotonEvent-> muDR1 = -99.;
728 fDiphotonEvent-> muDR2 = -99.;
729 fDiphotonEvent-> muIso1 = -99.;
730 fDiphotonEvent-> muIso2 = -99.;
731 fDiphotonEvent-> muIso3 = -99.;
732 fDiphotonEvent-> muIso4 = -99.;
733 fDiphotonEvent-> muD0 = -99.;
734 fDiphotonEvent-> muDZ = -99.;
735 fDiphotonEvent-> muChi2 = -99.;
736 fDiphotonEvent-> muNhits = -99;
737 fDiphotonEvent-> muNpixhits = -99;
738 fDiphotonEvent-> muNegs = -99;
739 fDiphotonEvent-> muNMatch = -99;
740
741 // Electron Stuff
742 fDiphotonEvent-> elePt = -99.;
743 fDiphotonEvent-> eleEta = -99.;
744 fDiphotonEvent-> eleSCEta = -99.;
745 fDiphotonEvent-> eleIso1 = -99.;
746 fDiphotonEvent-> eleIso2 = -99.;
747 fDiphotonEvent-> eleIso3 = -99.;
748 fDiphotonEvent-> eleIso4 = -99.;
749 fDiphotonEvent-> eleDist = -99.;
750 fDiphotonEvent-> eleDcot = -99.;
751 fDiphotonEvent-> eleCoviee = -99.;
752 fDiphotonEvent-> eleDphiin = -99.;
753 fDiphotonEvent-> eleDetain = -99.;
754 fDiphotonEvent-> eleDR1 = -99.;
755 fDiphotonEvent-> eleDR2 = -99.;
756 fDiphotonEvent-> eleMass1 = -99.;
757 fDiphotonEvent-> eleMass2 = -99.;
758 fDiphotonEvent-> eleNinnerHits = -99;
759
760 if( fApplyLeptonTag ) {
761 // perform lepton tagging
762 // the diphoton event record will have one more entry; i.e. leptonTag
763 // leptonTag = -1 -> lepton-taggng was swicthed off
764 // = 0 -> event tagged as 'non-lepton-event'
765 // = +1 -> event tagged as muon-event
766 // = +2 -> event tagged as electron-event
767 fDiphotonEvent->leptonTag = 0;
768
769 if ( fLeptonTagMuons->GetEntries() > 0 ) {
770 // need to have dR > 1 for with respect to both photons ***changed to 0.7 for 2012
771 if( (MathUtils::DeltaR(fLeptonTagMuons->At(0),phHard) >= 1.0) &&
772 (MathUtils::DeltaR(fLeptonTagMuons->At(0),phSoft) >= 1.0) &&
773 ((phHard->Pt()/(phHard->Mom() + phSoft->Mom()).M())>(45./120.)) &&
774 ((phSoft->Pt()/(phHard->Mom() + phSoft->Mom()).M())>(30./120.)) ){
775
776 fDiphotonEvent->leptonTag = 1;
777
778 fDiphotonEvent-> muonPt = fLeptonTagMuons->At(0)->Pt();
779 fDiphotonEvent-> muonEta = fLeptonTagMuons->At(0)->Eta();
780 fDiphotonEvent-> muDR1 = MathUtils::DeltaR(fLeptonTagMuons->At(0),phHard);
781 fDiphotonEvent-> muDR2 = MathUtils::DeltaR(fLeptonTagMuons->At(0),phSoft);
782
783 fDiphotonEvent-> muIso1 = (fLeptonTagMuons->At(0)->IsoR03SumPt() + fLeptonTagMuons->At(0)->IsoR03EmEt() + fLeptonTagMuons->At(0)->IsoR03HadEt() - fPileUpDen->At(0)->RhoRandomLowEta() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(0)->Pt();
784 fDiphotonEvent-> muIso2 = (fLeptonTagMuons->At(0)->IsoR03SumPt() + fLeptonTagMuons->At(0)->IsoR03EmEt() + fLeptonTagMuons->At(0)->IsoR03HadEt() - fPileUpDen->At(0)->RhoRandom() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(0)->Pt();
785 fDiphotonEvent-> muIso3 = (fLeptonTagMuons->At(0)->IsoR03SumPt() + fLeptonTagMuons->At(0)->IsoR03EmEt() + fLeptonTagMuons->At(0)->IsoR03HadEt() - fPileUpDen->At(0)->RhoLowEta() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(0)->Pt();
786 fDiphotonEvent-> muIso4 = (fLeptonTagMuons->At(0)->IsoR03SumPt() + fLeptonTagMuons->At(0)->IsoR03EmEt() + fLeptonTagMuons->At(0)->IsoR03HadEt() - fPileUpDen->At(0)->Rho() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(0)->Pt();
787 fDiphotonEvent-> muD0 = TMath::Abs(fLeptonTagMuons->At(0)->BestTrk()->D0Corrected(*fPV->At(0)));
788 fDiphotonEvent-> muDZ = TMath::Abs(fLeptonTagMuons->At(0)->BestTrk()->DzCorrected(*fPV->At(0)));
789 fDiphotonEvent-> muChi2 = fLeptonTagMuons->At(0)->GlobalTrk()->Chi2()/fLeptonTagMuons->At(0)->GlobalTrk()->Ndof();
790
791 fDiphotonEvent-> muNhits = fLeptonTagMuons->At(0)->BestTrk()->NHits();
792 fDiphotonEvent-> muNpixhits = fLeptonTagMuons->At(0)->BestTrk()->NPixelHits();
793 fDiphotonEvent-> muNegs = fLeptonTagMuons->At(0)->NSegments();
794 fDiphotonEvent-> muNMatch = fLeptonTagMuons->At(0)->NMatches();
795 }
796 }
797 if ( fDiphotonEvent->leptonTag < 1 && fLeptonTagElectrons->GetEntries() > 0 ) {
798 if( (MathUtils::DeltaR(fLeptonTagElectrons->At(0),phHard) >= 1) &&
799 (MathUtils::DeltaR(fLeptonTagElectrons->At(0),phSoft) >= 1) &&
800 (PhotonTools::ElectronVetoCiC(phHard,fLeptonTagElectrons) >= 1) &&
801 (PhotonTools::ElectronVetoCiC(phSoft,fLeptonTagElectrons) >= 1) &&
802 (TMath::Abs( (phHard->Mom()+fLeptonTagElectrons->At(0)->Mom()).M()-91.19 ) >= 10) &&
803 (TMath::Abs( (phSoft->Mom()+fLeptonTagElectrons->At(0)->Mom()).M()-91.19 ) >= 10) &&
804 ((phHard->Pt()/(phHard->Mom() + phSoft->Mom()).M())>(45./120.)) &&
805 ((phSoft->Pt()/(phHard->Mom() + phSoft->Mom()).M())>(30./120.)) ){
806
807 fDiphotonEvent->leptonTag = 2;
808
809 fDiphotonEvent-> elePt = fLeptonTagElectrons->At(0)->Pt();
810 fDiphotonEvent-> eleEta = fLeptonTagElectrons->At(0)->Eta();
811 fDiphotonEvent-> eleSCEta = fLeptonTagElectrons->At(0)->SCluster()->Eta();
812 fDiphotonEvent-> eleIso1 = (fLeptonTagElectrons->At(0)->TrackIsolationDr03() + fLeptonTagElectrons->At(0)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(0)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->RhoRandomLowEta() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
813
814 fDiphotonEvent-> eleIso2 = -99.;
815 // fDiphotonEvent-> eleIso2 = ElectronTools::ElectronEffectiveArea(ElectronTools::kEleGammaIso03,fLeptonTagElectrons->At(0)->SCluster()->Eta(), ElectronTools::kEleEAData2012) + ElectronTools::ElectronEffectiveArea(ElectronTools::kEleNeutralHadronIso03, fLeptonTagElectrons->At(0)->SCluster()->Eta(), ElectronTools::kEleEAData2012) ;
816
817 fDiphotonEvent-> eleIso3 = (fLeptonTagElectrons->At(0)->TrackIsolationDr03() + fLeptonTagElectrons->At(0)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(0)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->RhoLowEta() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
818 fDiphotonEvent-> eleIso4 = (fLeptonTagElectrons->At(0)->TrackIsolationDr03() + fLeptonTagElectrons->At(0)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(0)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->Rho() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
819 fDiphotonEvent-> eleDist = fLeptonTagElectrons->At(0)->ConvPartnerDist();
820 fDiphotonEvent-> eleDcot = fLeptonTagElectrons->At(0)->ConvPartnerDCotTheta();
821 fDiphotonEvent-> eleCoviee = fLeptonTagElectrons->At(0)->CoviEtaiEta();
822 fDiphotonEvent-> eleDphiin = TMath::Abs(fLeptonTagElectrons->At(0)->DeltaPhiSuperClusterTrackAtVtx());
823 fDiphotonEvent-> eleDetain = TMath::Abs(fLeptonTagElectrons->At(0)->DeltaEtaSuperClusterTrackAtVtx());
824 fDiphotonEvent-> eleDR1 = MathUtils::DeltaR(fLeptonTagElectrons->At(0),phHard);
825 fDiphotonEvent-> eleDR2 = MathUtils::DeltaR(fLeptonTagElectrons->At(0),phSoft);
826 fDiphotonEvent-> eleMass1 = (phHard->Mom()+fLeptonTagElectrons->At(0)->Mom()).M();
827 fDiphotonEvent-> eleMass2 = (phSoft->Mom()+fLeptonTagElectrons->At(0)->Mom()).M();
828 fDiphotonEvent-> eleNinnerHits = fLeptonTagElectrons->At(0)->Trk()->NExpectedHitsInner();
829 }
830 }
831 }
832
833 //vbf tag
834 fDiphotonEvent->vbfTag = -1;
835 fDiphotonEvent->vbfbdt = -99;
836 if(fApplyVBFTag){
837 fDiphotonEvent->vbfTag = 0;
838 jet1pt_vbf = fDiphotonEvent->jet1pt;
839 jet2pt_vbf = fDiphotonEvent->jet2pt;
840 deltajeteta_vbf = TMath::Abs(fDiphotonEvent->jet1eta - fDiphotonEvent->jet2eta);
841 dijetmass_vbf = fDiphotonEvent->dijetmass;
842 zeppenfeld_vbf = TMath::Abs(fDiphotonEvent->zeppenfeld);
843 dphidijetgg_vbf = TMath::Abs(fDiphotonEvent->dphidijetgg);
844 diphoptOverdiphomass_vbf = fDiphotonEvent->ptgg/fDiphotonEvent->mass;
845 pho1ptOverdiphomass_vbf = phHard->Pt()/fDiphotonEvent->mass;
846 pho2ptOverdiphomass_vbf = phSoft->Pt()/fDiphotonEvent->mass;
847
848 if(jet1 && jet2 && (pho1ptOverdiphomass_vbf > 40./120.) && (pho2ptOverdiphomass_vbf > 30./120.) && (jet1pt_vbf > 30) && (jet2pt_vbf > 20) && (dijetmass_vbf > 250)){
849 fDiphotonEvent->vbfbdt = fMVAVBF.GetMVAbdtValue(jet1pt_vbf,jet2pt_vbf,deltajeteta_vbf,dijetmass_vbf,zeppenfeld_vbf,dphidijetgg_vbf,diphoptOverdiphomass_vbf,pho1ptOverdiphomass_vbf,pho2ptOverdiphomass_vbf);
850
851 if((fDiphotonEvent->vbfbdt>=0.985) && (fDiphotonEvent->vbfbdt<=1)){
852 fDiphotonEvent->vbfTag = 2;
853 }
854 if((fDiphotonEvent->vbfbdt>=0.93) && (fDiphotonEvent->vbfbdt<0.985)){
855 fDiphotonEvent->vbfTag = 1;
856 }
857 }
858 }
859 //printf("vbfbdt:%f\n",fDiphotonEvent->vbfbdt);
860 if (fWriteDiphotonTree)
861 hCiCTuple->Fill();
862
863 if (fFillVertexTree && phHard && phSoft) {
864 for (UInt_t ivtx = 0; ivtx<fPV->GetEntries(); ++ivtx) {
865 const Vertex *v = fPV->At(ivtx);
866 fDiphotonVtx->SetVars(v,phHard,phSoft,fPFCands,ivtx,fPV->GetEntries(),_decayZ);
867 hVtxTree->Fill();
868 }
869 }
870
871 }
872
873 if (!fWriteSingleTree)
874 return;
875
876 for (UInt_t iph = 0; iph<egcol->GetEntries(); ++iph) {
877 const Particle *p = 0;
878 const Photon *ph = 0;
879 const Electron *ele = 0;
880 const SuperCluster *sc = 0;
881 if (fLoopOnGoodElectrons) {
882 ele = fGoodElectrons->At(iph);
883 p = ele;
884 sc = ele->SCluster();
885 ph = PhotonTools::MatchedPhoton(ele,fPhotons);
886 }
887 else {
888 ph = fPhotons->At(iph);
889 p = ph;
890 sc = ph->SCluster();
891 ele = PhotonTools::MatchedElectron(ph,fGoodElectrons);
892 }
893
894 const DecayParticle *conv = PhotonTools::MatchedConversion(sc,fConversions,bsp,
895 nhitsbeforevtxmax);
896 const SuperCluster *pfsc = PhotonTools::MatchedPFSC(sc,fPFPhotons, fElectrons);
897
898 if (!fLoopOnGoodElectrons && ph->HasPV()) {
899 fDiphotonEvent->vtxZ = ph->PV()->Z();
900 }
901
902 const MCParticle *phgen = 0;
903 if( !fIsData ) {
904 phgen = PhotonTools::MatchMC(p,fMCParticles,!fApplyElectronVeto);
905 }
906
907 if (fExcludeSinglePrompt && phgen) return;
908
909 fDiphotonEvent->mt = -99.;
910 fDiphotonEvent->cosphimet = -99.;
911 fDiphotonEvent->mtele = -99.;
912 fDiphotonEvent->cosphimetele = -99.;
913
914 if (ph) {
915 fDiphotonEvent->cosphimet = TMath::Cos(ph->Phi()-fPFMet->At(0)->Phi());
916 fDiphotonEvent->mt = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ph->Pt()*
917 (1.0-fDiphotonEvent->cosphimet));
918 }
919
920 if (ele) {
921 fDiphotonEvent->cosphimetele = TMath::Cos(ele->Phi()-fPFMet->At(0)->Phi());
922 fDiphotonEvent->mtele = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ele->Pt()*
923 (1.0-fDiphotonEvent->cosphimetele));
924 }
925
926 fSinglePhoton->SetVars(ph,conv,ele,pfsc,phgen,fPhfixph,fPhfixele,fTracks,fPV,fPFCands,rho,fFillClusterArrays,
927 fElectrons,fConversions,bsp,fApplyElectronVeto);
928 hCiCTupleSingle->Fill();
929 }
930
931
932
933 return;
934 }
935
936 //--------------------------------------------------------------------------------------------------
937 void PhotonTreeWriter::SlaveBegin()
938 {
939 // Run startup code on the computer (slave) doing the actual analysis. Here,
940 // we just request the photon collection branch.
941
942 if( fApplyLeptonTag ) {
943 ReqEventObject(fLeptonTagElectronsName, fLeptonTagElectrons, false);
944 ReqEventObject(fLeptonTagMuonsName, fLeptonTagMuons, false);
945 }
946
947 // ReqEventObject(fPFNoPileUpName, fPFNoPileUpCands, false);
948 // ReqEventObject(fPFPileUpName, fPFPileUpCands, false);
949
950
951
952
953
954 ReqEventObject(fPhotonBranchName,fPhotons, fPhotonsFromBranch);
955 if (fEnablePFPhotons) ReqEventObject(fPFPhotonName,fPFPhotons, true);
956 ReqEventObject(fTrackBranchName, fTracks, true);
957 ReqEventObject(fElectronName, fElectrons, true);
958 ReqEventObject(fGoodElectronName,fGoodElectrons,fGoodElectronsFromBranch);
959 ReqEventObject(fPileUpDenName, fPileUpDen, true);
960 ReqEventObject(fPVName, fPV, fPVFromBranch);
961 ReqEventObject(fConversionName, fConversions, true);
962 ReqEventObject(fBeamspotName, fBeamspot, true);
963 ReqEventObject(fPFCandName, fPFCands, true);
964 ReqEventObject(fSuperClusterName,fSuperClusters,true);
965 ReqEventObject(fPFMetName, fPFMet, true);
966 if (fEnableJets){
967 ReqEventObject(fPFJetName, fPFJets, fPFJetsFromBranch);
968 ReqBranch(funcorrPFJetName, funcorrPFJets);
969 // if (!fIsData) ReqEventObject(fGenJetName, fGenJets, true);
970 }
971 if (!fIsData) {
972 ReqBranch(fPileUpName, fPileUp);
973 ReqBranch(fMCParticleName, fMCParticles);
974 ReqBranch(fMCEventInfoName, fMCEventInfo);
975 if (fEnableJets) ReqEventObject(fGenJetName, fGenJets, true);
976 }
977 if (fIsData) {
978 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
979 TString("/src/MitPhysics/data/PhotonFixGRPV22.dat");
980 }
981 else {
982 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
983 TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat");
984 }
985
986 fPhfixph.initialise("4_2",std::string(fPhFixDataFile));
987 fPhfixele.initialise("4_2e",std::string(fPhFixDataFile));
988
989 // fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
990 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
991 // TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
992 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_42.root"))),
993 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_42.root"))),
994 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1_42.root"))),
995 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2_42.root")))
996 // );
997
998 fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
999 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
1000 TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
1001 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_52.root"))),
1002 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_52.root"))),
1003 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1cov_52.root"))),
1004 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2cov_52.root")))
1005 );
1006
1007 fJetId.Initialize(JetIDMVA::kMedium,
1008 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
1009 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
1010 JetIDMVA::kCut,
1011 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")))
1012
1013 );
1014
1015 fMVAVBF.InitializeMVA();
1016
1017 fDiphotonEvent = new PhotonTreeWriterDiphotonEvent;
1018 fSinglePhoton = new PhotonTreeWriterPhoton<16>;
1019
1020 TFile *ftmp = TFile::Open(TString::Format("%s_tmp.root",GetName()),"RECREATE");
1021
1022 if (fWriteDiphotonTree) {
1023 hCiCTuple = new TTree(fTupleName.Data(),fTupleName.Data());
1024 hCiCTuple->SetAutoSave(300e9);
1025 }
1026 TString singlename = fTupleName + TString("Single");
1027 if (fWriteSingleTree) {
1028 hCiCTupleSingle = new TTree(singlename,singlename);
1029 hCiCTupleSingle->SetAutoSave(300e9);
1030 }
1031
1032 //make flattish tree from classes so we don't have to rely on dictionaries for reading later
1033 TClass *eclass = TClass::GetClass("mithep::PhotonTreeWriterDiphotonEvent");
1034 TClass *pclass = TClass::GetClass("mithep::PhotonTreeWriterPhoton<16>");
1035 TList *elist = eclass->GetListOfDataMembers();
1036 TList *plist = pclass->GetListOfDataMembers();
1037
1038 for (int i=0; i<elist->GetEntries(); ++i) {
1039 const TDataMember *tdm = static_cast<const TDataMember*>(elist->At(i));//ming
1040 if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
1041 TString typestring;
1042 if (TString(tdm->GetTypeName()).BeginsWith("Char_t")) typestring = "B";
1043 else if (TString(tdm->GetTypeName()).BeginsWith("UChar_t")) typestring = "b";
1044 else if (TString(tdm->GetTypeName()).BeginsWith("Short_t")) typestring = "S";
1045 else if (TString(tdm->GetTypeName()).BeginsWith("UShort_t")) typestring = "s";
1046 else if (TString(tdm->GetTypeName()).BeginsWith("Int_t")) typestring = "I";
1047 else if (TString(tdm->GetTypeName()).BeginsWith("UInt_t")) typestring = "i";
1048 else if (TString(tdm->GetTypeName()).BeginsWith("Float_t")) typestring = "F";
1049 else if (TString(tdm->GetTypeName()).BeginsWith("Double_t")) typestring = "D";
1050 else if (TString(tdm->GetTypeName()).BeginsWith("Long64_t")) typestring = "L";
1051 else if (TString(tdm->GetTypeName()).BeginsWith("ULong64_t")) typestring = "l";
1052 else if (TString(tdm->GetTypeName()).BeginsWith("Bool_t")) typestring = "O";
1053 else continue;
1054 //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
1055 Char_t *addr = (Char_t*)fDiphotonEvent;//ming:?
1056 assert(sizeof(Char_t)==1);
1057 if (fWriteDiphotonTree) hCiCTuple->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
1058 if (fWriteSingleTree) hCiCTupleSingle->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
1059 }
1060
1061 for (int iph=0; iph<2; ++iph) {
1062 for (int i=0; i<plist->GetEntries(); ++i) {
1063 const TDataMember *tdm = static_cast<const TDataMember*>(plist->At(i));
1064 if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
1065 TString typestring;
1066 if (TString(tdm->GetTypeName()).BeginsWith("Char_t")) typestring = "B";
1067 else if (TString(tdm->GetTypeName()).BeginsWith("UChar_t")) typestring = "b";
1068 else if (TString(tdm->GetTypeName()).BeginsWith("Short_t")) typestring = "S";
1069 else if (TString(tdm->GetTypeName()).BeginsWith("UShort_t")) typestring = "s";
1070 else if (TString(tdm->GetTypeName()).BeginsWith("Int_t")) typestring = "I";
1071 else if (TString(tdm->GetTypeName()).BeginsWith("UInt_t")) typestring = "i";
1072 else if (TString(tdm->GetTypeName()).BeginsWith("Float_t")) typestring = "F";
1073 else if (TString(tdm->GetTypeName()).BeginsWith("Double_t")) typestring = "D";
1074 else if (TString(tdm->GetTypeName()).BeginsWith("Long64_t")) typestring = "L";
1075 else if (TString(tdm->GetTypeName()).BeginsWith("ULong64_t")) typestring = "l";
1076 else if (TString(tdm->GetTypeName()).BeginsWith("Bool_t")) typestring = "O";
1077 else continue;
1078 //printf("%s\n",tdm->GetTypeName());
1079 TString varname = TString::Format("ph%d.%s",iph+1,tdm->GetName());
1080 if (tdm->GetArrayDim()==1) {
1081 varname = TString::Format("%s[%i]",varname.Data(),tdm->GetMaxIndex(0));
1082 }
1083
1084 //printf("typename = %s, arraydim = %i, arraysize = %i,varname = %s\n", tdm->GetTypeName(), tdm->GetArrayDim(), tdm->GetMaxIndex(0), varname.Data());
1085
1086 Char_t *addr = (Char_t*)&fDiphotonEvent->photons[iph];
1087 assert(sizeof(Char_t)==1);
1088 if (fWriteDiphotonTree) hCiCTuple->Branch(varname,addr+tdm->GetOffset(),TString::Format("%s/%s",varname.Data(),typestring.Data()));
1089
1090 if (iph==0) {
1091 TString singlename = TString::Format("ph.%s",tdm->GetName());
1092 if (tdm->GetArrayDim()==1) {
1093 singlename = TString::Format("%s[%i]",singlename.Data(),tdm->GetMaxIndex(0));
1094 }
1095 Char_t *addrsingle = (Char_t*)fSinglePhoton;
1096 if (fWriteSingleTree) hCiCTupleSingle->Branch(singlename,addrsingle+tdm->GetOffset(),TString::Format("%s/%s",singlename.Data(),typestring.Data()));
1097 }
1098 }
1099 }
1100
1101 if (fWriteDiphotonTree)
1102 AddOutput(hCiCTuple);
1103 if (fWriteSingleTree)
1104 AddOutput(hCiCTupleSingle);
1105
1106 if (fFillVertexTree) {
1107 fDiphotonVtx = new PhotonTreeWriterVtx;
1108 hVtxTree = new TTree("hVtxTree","hVtxTree");
1109 hVtxTree->SetAutoSave(300e9);
1110 AddOutput(hVtxTree);
1111
1112 TClass *vclass = TClass::GetClass("mithep::PhotonTreeWriterVtx");
1113 TList *vlist = vclass->GetListOfDataMembers();
1114
1115 for (int i=0; i<vlist->GetEntries(); ++i) {
1116 const TDataMember *tdm = static_cast<const TDataMember*>(vlist->At(i));
1117 if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
1118 TString typestring;
1119 if (TString(tdm->GetTypeName()).BeginsWith("Char_t")) typestring = "B";
1120 else if (TString(tdm->GetTypeName()).BeginsWith("UChar_t")) typestring = "b";
1121 else if (TString(tdm->GetTypeName()).BeginsWith("Short_t")) typestring = "S";
1122 else if (TString(tdm->GetTypeName()).BeginsWith("UShort_t")) typestring = "s";
1123 else if (TString(tdm->GetTypeName()).BeginsWith("Int_t")) typestring = "I";
1124 else if (TString(tdm->GetTypeName()).BeginsWith("UInt_t")) typestring = "i";
1125 else if (TString(tdm->GetTypeName()).BeginsWith("Float_t")) typestring = "F";
1126 else if (TString(tdm->GetTypeName()).BeginsWith("Double_t")) typestring = "D";
1127 else if (TString(tdm->GetTypeName()).BeginsWith("Long64_t")) typestring = "L";
1128 else if (TString(tdm->GetTypeName()).BeginsWith("ULong64_t")) typestring = "l";
1129 else if (TString(tdm->GetTypeName()).BeginsWith("Bool_t")) typestring = "O";
1130 else continue;
1131 //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
1132 Char_t *addr = (Char_t*)fDiphotonVtx;
1133 assert(sizeof(Char_t)==1);
1134 hVtxTree->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
1135 }
1136 }
1137
1138 }
1139
1140 // ----------------------------------------------------------------------------------------
1141 // some helpfer functions....
1142 void PhotonTreeWriter::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass)
1143 {
1144 pt = -999.;
1145 decayZ = -999.;
1146 mass = -999.;
1147
1148 // loop over all GEN particles and look for status 1 photons
1149 for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
1150 const MCParticle* p = fMCParticles->At(i);
1151 if (p->Is(MCParticle::kH) || (!fApplyElectronVeto &&
1152 (p->AbsPdgId()==23 || p->AbsPdgId()==24))) {
1153 pt = p->Pt();
1154 decayZ = p->DecayVertex().Z();
1155 mass = p->Mass();
1156 break;
1157 }
1158 }
1159
1160 return;
1161 }
1162
1163
1164 Float_t PhotonTreeWriter::GetEventCat(PhotonTools::CiCBaseLineCats cat1,
1165 PhotonTools::CiCBaseLineCats cat2) {
1166
1167 bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
1168 bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
1169
1170 bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
1171 bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
1172
1173 if( ph1IsEB && ph2IsEB )
1174 return ( ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
1175
1176 return ( ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
1177 }
1178
1179 template <int NClus>
1180 void PhotonTreeWriterPhoton<NClus>::SetVars(const Photon *p, const DecayParticle *c, const Electron *ele,
1181 const SuperCluster *pfsc, const MCParticle *m,
1182 PhotonFix &phfixph, PhotonFix &phfixele,
1183 const TrackCol* trackCol,const VertexCol* vtxCol,
1184 const PFCandidateCol* fPFCands,
1185 Double_t rho,
1186 Bool_t fillclusterarrays,
1187 const ElectronCol* els, const DecayParticleCol *convs, const BaseVertex *bs, Bool_t applyElectronVeto) {
1188
1189 const SuperCluster *s = 0;
1190 if (p)
1191 s = p->SCluster();
1192 else
1193 s = ele->SCluster();
1194 const BasicCluster *b = s->Seed();
1195 const BasicCluster *b2 = 0;
1196 Double_t ebcmax = -99.;
1197 for (UInt_t i=0; i<s->ClusterSize(); ++i) {
1198 const BasicCluster *bc = s->Cluster(i);
1199 if (bc->Energy() > ebcmax && bc !=b) {
1200 b2 = bc;
1201 ebcmax = bc->Energy();
1202 }
1203 }
1204
1205 const BasicCluster *bclast = 0;
1206 Double_t ebcmin = 1e6;
1207 for (UInt_t i=0; i<s->ClusterSize(); ++i) {
1208 const BasicCluster *bc = s->Cluster(i);
1209 if (bc->Energy() < ebcmin && bc !=b) {
1210 bclast = bc;
1211 ebcmin = bc->Energy();
1212 }
1213 }
1214
1215 const BasicCluster *bclast2 = 0;
1216 ebcmin = 1e6;
1217 for (UInt_t i=0; i<s->ClusterSize(); ++i) {
1218 const BasicCluster *bc = s->Cluster(i);
1219 if (bc->Energy() < ebcmin && bc !=b && bc!=bclast) {
1220 bclast2 = bc;
1221 ebcmin = bc->Energy();
1222 }
1223 }
1224
1225 if (p) {
1226 hasphoton = kTRUE;
1227 e = p->E();
1228 pt = p->Pt();
1229 eta = p->Eta();
1230 phi = p->Phi();
1231 r9 = p->R9();
1232 e3x3 = p->E33();
1233 e5x5 = p->E55();
1234 hovere = p->HadOverEm();
1235 hoveretower = p->HadOverEmTow();
1236 sigietaieta = p->CoviEtaiEta();
1237 phcat = PhotonTools::CiCBaseLineCat(p);
1238 eerr = p->EnergyErr();
1239 eerrsmeared = p->EnergyErrSmeared();
1240 esmearing = p->EnergySmearing();
1241 idmva = p->IdMva();
1242 hcalisodr03 = p->HcalTowerSumEtDr03();
1243 ecalisodr03 = p->EcalRecHitIsoDr03();
1244 trkisohollowdr03 = p->HollowConeTrkIsoDr03();
1245 hcalisodr04 = p->HcalTowerSumEtDr04();
1246 ecalisodr04 = p->EcalRecHitIsoDr04();
1247 trkisohollowdr04 = p->HollowConeTrkIsoDr04();
1248
1249 passeleveto = PhotonTools::PassElectronVetoConvRecovery(p, els, convs, bs);
1250
1251 const Vertex *vtx = vtxCol->At(0);
1252 if (p->HasPV()) vtx = p->PV();
1253
1254 UInt_t wVtxInd = 0;
1255
1256 trackiso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0,
1257 trackCol, NULL, NULL,
1258 (!applyElectronVeto ? els : NULL) );
1259 //Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
1260
1261 // track iso worst vtx
1262 trackiso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0,
1263 trackCol, &wVtxInd,vtxCol,
1264 (!applyElectronVeto ? els : NULL) );
1265 combiso1 = ecalisodr03+hcalisodr04+trackiso1 - 0.17*rho;
1266 combiso2 = ecalisodr04+hcalisodr04+trackiso2 - 0.52*rho;
1267
1268
1269 // -----------------------------------------------------
1270 // PF-CiC4 Debug Stuff
1271 std::vector<double> debugVals;
1272 bool tmpPass = PhotonTools::PassCiCPFIsoSelection(p, vtx, fPFCands, vtxCol, rho, 20., &debugVals);
1273 if( debugVals.size() == 13 ) {
1274 pfcic4_tIso1 = debugVals[0];
1275 pfcic4_tIso2 = debugVals[1];
1276 pfcic4_tIso3 = debugVals[2];
1277 pfcic4_covIEtaIEta = debugVals[3];
1278 pfcic4_HoE = debugVals[4];
1279 pfcic4_R9 = debugVals[5];
1280 pfcic4_wVtxInd = debugVals[6];
1281 pfcic4_ecalIso3 = debugVals[7];
1282 pfcic4_ecalIso4 = debugVals[8];
1283 pfcic4_trackIsoSel03 = debugVals[9];
1284 pfcic4_trackIsoWorst04 = debugVals[10];
1285 pfcic4_combIso1 = debugVals[11];
1286 pfcic4_combIso2 = debugVals[12];
1287 }
1288 // -----------------------------------------------------
1289 //id mva
1290 //2011
1291 idmva_tIso1abs=combiso1;
1292 idmva_tIso2abs=combiso2;
1293 idmva_tIso3abs=trackiso1;
1294 idmva_absIsoEcal=ecalisodr03;
1295 idmva_absIsoHcal=hcalisodr04;
1296 //2012
1297 idmva_CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
1298 idmva_s4ratio=p->S4Ratio();
1299 idmva_GammaIso=IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
1300 idmva_ChargedIso_selvtx=IsolationTools::PFChargedIsolation(p,vtx,0.3,0.,fPFCands);
1301 idmva_ChargedIso_worstvtx=IsolationTools::PFChargedIsolation(p,vtx,0.3,0.,fPFCands,&wVtxInd,vtxCol);
1302 idmva_PsEffWidthSigmaRR=p->EffSigmaRR();
1303 }
1304 else {
1305 hasphoton = kFALSE;
1306 e = -99.;
1307 pt = -99.;
1308 eta = -99.;
1309 phi = -99.;
1310 r9 = b->E3x3()/s->RawEnergy();
1311 e3x3 = b->E3x3();
1312 e5x5 = b->E5x5();
1313 hovere = ele->HadronicOverEm();
1314 sigietaieta = ele->CoviEtaiEta();
1315 phcat = -99;
1316 eerr = -99.;
1317 eerrsmeared = -99.;
1318 esmearing = 0.;
1319 idmva = -99.;
1320 hcalisodr03 = -99;
1321 ecalisodr03 = -99;
1322 trkisohollowdr03 = -99;
1323 hcalisodr04 = -99;
1324 ecalisodr04 = -99;
1325 trkisohollowdr04 = -99;
1326 trackiso1 = -99.;
1327 trackiso2 = -99.;
1328 combiso1 = -99.;
1329 combiso2 = -99.;
1330 }
1331
1332 sce = s->Energy();
1333 scrawe = s->RawEnergy();
1334 scpse = s->PreshowerEnergy();
1335 scpssigmaxx = s->PsEffWidthSigmaXX();
1336 scpssigmayy = s->PsEffWidthSigmaYY();
1337 sceta = s->Eta();
1338 scphi = s->Phi();
1339 scnclusters = s->ClusterSize();
1340 scnhits = s->NHits();
1341 scetawidth = -99.;
1342 scphiwidth = -99.;
1343 if (p) {
1344 scetawidth = p->EtaWidth();
1345 scphiwidth = p->PhiWidth();
1346 }
1347 else {
1348 scetawidth = s->EtaWidth();
1349 scphiwidth = s->PhiWidth();
1350 }
1351 isbarrel = (s->AbsEta()<1.5);
1352 isr9reco = (isbarrel && r9>0.94) || (!isbarrel && r9>0.95);
1353 isr9cat = (r9>0.94);
1354
1355 eseed = b->Energy();
1356 etaseed = b->Eta();
1357 phiseed = b->Phi();
1358 ietaseed = b->IEta();
1359 iphiseed = b->IPhi();
1360 ixseed = b->IX();
1361 iyseed = b->IY();
1362 etacryseed = b->EtaCry();
1363 phicryseed = b->PhiCry();
1364 xcryseed = b->XCry();
1365 ycryseed = b->YCry();
1366 thetaaxisseed = b->ThetaAxis();
1367 phiaxisseed = b->PhiAxis();
1368 sigietaietaseed = TMath::Sqrt(b->CoviEtaiEta());
1369 sigiphiphiseed = TMath::Sqrt(b->CoviPhiiPhi());
1370 if (isnan(sigiphiphiseed)) sigiphiphiseed = -99.;
1371 covietaiphiseed = b->CoviEtaiPhi();
1372 if (isnan(covietaiphiseed)) covietaiphiseed = -99.;
1373 e3x3seed = b->E3x3();
1374 e5x5seed = b->E5x5();
1375 emaxseed = b->EMax();
1376 e2ndseed = b->E2nd();
1377 etopseed = b->ETop();
1378 ebottomseed = b->EBottom();
1379 eleftseed = b->ELeft();
1380 erightseed = b->ERight();
1381 e1x3seed = b->E1x3();
1382 e3x1seed = b->E3x1();
1383 e1x5seed = b->E1x5();
1384 e2x2seed = b->E2x2();
1385 e4x4seed = b->E4x4();
1386 e2x5maxseed = b->E2x5Max();
1387 e2x5topseed = b->E2x5Top();
1388 e2x5bottomseed = b->E2x5Bottom();
1389 e2x5leftseed = b->E2x5Left();
1390 e2x5rightseed = b->E2x5Right();
1391 xseedseed = b->Pos().X();
1392 yseedseed = b->Pos().Y();
1393 zseedseed = b->Pos().Z();
1394 nhitsseed = b->NHits();
1395
1396 if (b2) {
1397 ebc2 = b2->Energy();
1398 etabc2 = b2->Eta();
1399 phibc2 = b2->Phi();
1400 ietabc2 = b2->IEta();
1401 iphibc2 = b2->IPhi();
1402 ixbc2 = b2->IX();
1403 iybc2 = b2->IY();
1404 etacrybc2 = b2->EtaCry();
1405 phicrybc2 = b2->PhiCry();
1406 xcrybc2 = b2->XCry();
1407 ycrybc2 = b2->YCry();
1408 thetaaxisbc2 = b2->ThetaAxis();
1409 phiaxisbc2 = b2->PhiAxis();
1410 sigietaietabc2 = TMath::Sqrt(b2->CoviEtaiEta());
1411 sigiphiphibc2 = TMath::Sqrt(b2->CoviPhiiPhi());
1412 if (isnan(sigiphiphibc2)) sigiphiphibc2 = -99.;
1413 covietaiphibc2 = b2->CoviEtaiPhi();
1414 if (isnan(covietaiphibc2)) covietaiphibc2 = -99.;
1415 e3x3bc2 = b2->E3x3();
1416 e5x5bc2 = b2->E5x5();
1417 emaxbc2 = b2->EMax();
1418 e2ndbc2 = b2->E2nd();
1419 etopbc2 = b2->ETop();
1420 ebottombc2 = b2->EBottom();
1421 eleftbc2 = b2->ELeft();
1422 erightbc2 = b2->ERight();
1423 e1x3bc2 = b2->E1x3();
1424 e3x1bc2 = b2->E3x1();
1425 e1x5bc2 = b2->E1x5();
1426 e2x2bc2 = b2->E2x2();
1427 e4x4bc2 = b2->E4x4();
1428 e2x5maxbc2 = b2->E2x5Max();
1429 e2x5topbc2 = b2->E2x5Top();
1430 e2x5bottombc2 = b2->E2x5Bottom();
1431 e2x5leftbc2 = b2->E2x5Left();
1432 e2x5rightbc2 = b2->E2x5Right();
1433 xbc2bc2 = b2->Pos().X();
1434 ybc2bc2 = b2->Pos().Y();
1435 zbc2bc2 = b2->Pos().Z();
1436 nhitsbc2= b2->NHits();
1437 }
1438 else {
1439 ebc2 = 0;
1440 etabc2 = 0;
1441 phibc2 = 0;
1442 ietabc2 = 0;
1443 iphibc2 = 0;
1444 ixbc2 = 0;
1445 iybc2 = 0;
1446 etacrybc2 = 0;
1447 phicrybc2 = 0;
1448 xcrybc2 = 0;
1449 ycrybc2 = 0;
1450 thetaaxisbc2 = 0;
1451 phiaxisbc2 = 0;
1452 sigietaietabc2 = 0;
1453 sigiphiphibc2 = 0;
1454 covietaiphibc2 = 0;
1455 e3x3bc2 = 0;
1456 e5x5bc2 = 0;
1457 emaxbc2 = 0;
1458 e2ndbc2 = 0;
1459 etopbc2 = 0;
1460 ebottombc2 = 0;
1461 eleftbc2 = 0;
1462 erightbc2 = 0;
1463 e1x3bc2 = 0;
1464 e3x1bc2 = 0;
1465 e1x5bc2 = 0;
1466 e2x2bc2 = 0;
1467 e4x4bc2 = 0;
1468 e2x5maxbc2 = 0;
1469 e2x5topbc2 = 0;
1470 e2x5bottombc2 = 0;
1471 e2x5leftbc2 = 0;
1472 e2x5rightbc2 = 0;
1473 xbc2bc2 = 0;
1474 ybc2bc2 = 0;
1475 zbc2bc2 = 0;
1476 nhitsbc2 = 0;
1477 }
1478
1479 if (bclast) {
1480 ebclast = bclast->Energy();
1481 etabclast = bclast->Eta();
1482 phibclast = bclast->Phi();
1483 ietabclast = bclast->IEta();
1484 iphibclast = bclast->IPhi();
1485 ixbclast = bclast->IX();
1486 iybclast = bclast->IY();
1487 etacrybclast = bclast->EtaCry();
1488 phicrybclast = bclast->PhiCry();
1489 xcrybclast = bclast->XCry();
1490 ycrybclast = bclast->YCry();
1491 thetaaxisbclast = bclast->ThetaAxis();
1492 phiaxisbclast = bclast->PhiAxis();
1493 sigietaietabclast = TMath::Sqrt(bclast->CoviEtaiEta());
1494 sigiphiphibclast = TMath::Sqrt(bclast->CoviPhiiPhi());
1495 if (isnan(sigiphiphibclast)) sigiphiphibclast = -99.;
1496 covietaiphibclast = bclast->CoviEtaiPhi();
1497 if (isnan(covietaiphibclast)) covietaiphibclast = -99.;
1498 e3x3bclast = bclast->E3x3();
1499 nhitsbclast = bclast->NHits();
1500 }
1501 else {
1502 ebclast = 0;
1503 etabclast = 0;
1504 phibclast = 0;
1505 ietabclast = 0;
1506 iphibclast = 0;
1507 ixbclast = 0;
1508 iybclast = 0;
1509 etacrybclast = 0;
1510 phicrybclast = 0;
1511 xcrybclast = 0;
1512 ycrybclast = 0;
1513 thetaaxisbclast = 0;
1514 phiaxisbclast = 0;
1515 sigietaietabclast = 0;
1516 sigiphiphibclast = 0;
1517 covietaiphibclast = 0;
1518 e3x3bclast = 0;
1519 nhitsbclast = 0;
1520 }
1521
1522 if (bclast2) {
1523 ebclast2 = bclast2->Energy();
1524 etabclast2 = bclast2->Eta();
1525 phibclast2 = bclast2->Phi();
1526 ietabclast2 = bclast2->IEta();
1527 iphibclast2 = bclast2->IPhi();
1528 ixbclast2 = bclast2->IX();
1529 iybclast2 = bclast2->IY();
1530 etacrybclast2 = bclast2->EtaCry();
1531 phicrybclast2 = bclast2->PhiCry();
1532 xcrybclast2 = bclast2->XCry();
1533 ycrybclast2 = bclast2->YCry();
1534 thetaaxisbclast2 = bclast2->ThetaAxis();
1535 phiaxisbclast2 = bclast2->PhiAxis();
1536 sigietaietabclast2 = TMath::Sqrt(bclast2->CoviEtaiEta());
1537 sigiphiphibclast2 = TMath::Sqrt(bclast2->CoviPhiiPhi());
1538 if (isnan(sigiphiphibclast2)) sigiphiphibclast2 = -99.;
1539 covietaiphibclast2 = bclast2->CoviEtaiPhi();
1540 if (isnan(covietaiphibclast2)) covietaiphibclast2 = -99.;
1541 e3x3bclast2 = bclast2->E3x3();
1542 nhitsbclast2 = bclast2->NHits();
1543 }
1544 else {
1545 ebclast2 = 0;
1546 etabclast2 = 0;
1547 phibclast2 = 0;
1548 ietabclast2 = 0;
1549 iphibclast2 = 0;
1550 ixbclast2 = 0;
1551 iybclast2 = 0;
1552 etacrybclast2 = 0;
1553 phicrybclast2 = 0;
1554 xcrybclast2 = 0;
1555 ycrybclast2 = 0;
1556 thetaaxisbclast2 = 0;
1557 phiaxisbclast2 = 0;
1558 sigietaietabclast2 = 0;
1559 sigiphiphibclast2 = 0;
1560 covietaiphibclast2 = 0;
1561 e3x3bclast2 = 0;
1562 nhitsbclast2 = 0;
1563 }
1564
1565 for (UInt_t iclus=0; iclus<NClus; ++iclus) {
1566 if (fillclusterarrays && iclus < s->ClusterSize() ) {
1567 const BasicCluster *ib =s->Cluster(iclus);
1568
1569 ebcs[iclus] = ib->Energy();
1570 etabcs[iclus] = ib->Eta();
1571 phibcs[iclus] = ib->Phi();
1572 ietabcs[iclus] = ib->IEta();
1573 iphibcs[iclus] = ib->IPhi();
1574 ixbcs[iclus] = ib->IX();
1575 iybcs[iclus] = ib->IY();
1576 etacrybcs[iclus] = ib->EtaCry();
1577 phicrybcs[iclus] = ib->PhiCry();
1578 xcrybcs[iclus] = ib->XCry();
1579 ycrybcs[iclus] = ib->YCry();
1580 sigietaietabcs[iclus] = TMath::Sqrt(ib->CoviEtaiEta());
1581 sigiphiphibcs[iclus] = TMath::Sqrt(ib->CoviPhiiPhi());
1582 covietaiphibcs[iclus] = ib->CoviEtaiPhi();
1583 sigetaetabcs[iclus] = TMath::Sqrt(ib->CovEtaEta());
1584 sigphiphibcs[iclus] = TMath::Sqrt(ib->CovPhiPhi());
1585 covetaphibcs[iclus] = ib->CovEtaPhi();
1586 e3x3bcs[iclus] = ib->E3x3();
1587 e5x5bcs[iclus] = ib->E5x5();
1588 emaxbcs[iclus] = ib->EMax();
1589 e2ndbcs[iclus] = ib->E2nd();
1590 etopbcs[iclus] = ib->ETop();
1591 ebottombcs[iclus] = ib->EBottom();
1592 eleftbcs[iclus] = ib->ELeft();
1593 erightbcs[iclus] = ib->ERight();
1594 e1x3bcs[iclus] = ib->E1x3();
1595 e3x1bcs[iclus] = ib->E3x1();
1596 e1x5bcs[iclus] = ib->E1x5();
1597 e2x2bcs[iclus] = ib->E2x2();
1598 e4x4bcs[iclus] = ib->E4x4();
1599 e2x5maxbcs[iclus] = ib->E2x5Max();
1600 e2x5topbcs[iclus] = ib->E2x5Top();
1601 e2x5bottombcs[iclus] = ib->E2x5Bottom();
1602 e2x5leftbcs[iclus] = ib->E2x5Left();
1603 e2x5rightbcs[iclus] = ib->E2x5Right();
1604 nhitsbcs[iclus]= ib->NHits();
1605 }
1606 else {
1607 ebcs[iclus] = -999;
1608 etabcs[iclus] = -999;
1609 phibcs[iclus] = -999;
1610 ietabcs[iclus] = -999;
1611 iphibcs[iclus] = -999;
1612 ixbcs[iclus] = -999;
1613 iybcs[iclus] = -999;
1614 etacrybcs[iclus] = -999;
1615 phicrybcs[iclus] = -999;
1616 xcrybcs[iclus] = -999;
1617 ycrybcs[iclus] = -999;
1618 sigietaietabcs[iclus] = -999;
1619 sigiphiphibcs[iclus] = -999;
1620 covietaiphibcs[iclus] = -999;
1621 sigetaetabcs[iclus] = -999;
1622 sigphiphibcs[iclus] = -999;
1623 covetaphibcs[iclus] = -999;
1624 e3x3bcs[iclus] = -999;
1625 e5x5bcs[iclus] = -999;
1626 emaxbcs[iclus] = -999;
1627 e2ndbcs[iclus] = -999;
1628 etopbcs[iclus] = -999;
1629 ebottombcs[iclus] = -999;
1630 eleftbcs[iclus] = -999;
1631 erightbcs[iclus] = -999;
1632 e1x3bcs[iclus] = -999;
1633 e3x1bcs[iclus] = -999;
1634 e1x5bcs[iclus] = -999;
1635 e2x2bcs[iclus] = -999;
1636 e4x4bcs[iclus] = -999;
1637 e2x5maxbcs[iclus] = -999;
1638 e2x5topbcs[iclus] = -999;
1639 e2x5bottombcs[iclus] = -999;
1640 e2x5leftbcs[iclus] = -999;
1641 e2x5rightbcs[iclus] = -999;
1642 nhitsbcs[iclus] = -999;
1643 }
1644 }
1645
1646 for (UInt_t iclus=0; iclus<NClus; ++iclus) {
1647 if (fillclusterarrays && pfsc && iclus < pfsc->ClusterSize() ) {
1648 const BasicCluster *ib =pfsc->Cluster(iclus);
1649
1650 epfbcs[iclus] = ib->Energy();
1651 etapfbcs[iclus] = ib->Eta();
1652 phipfbcs[iclus] = ib->Phi();
1653 ietapfbcs[iclus] = ib->IEta();
1654 iphipfbcs[iclus] = ib->IPhi();
1655 ixpfbcs[iclus] = ib->IX();
1656 iypfbcs[iclus] = ib->IY();
1657 etacrypfbcs[iclus] = ib->EtaCry();
1658 phicrypfbcs[iclus] = ib->PhiCry();
1659 xcrypfbcs[iclus] = ib->XCry();
1660 ycrypfbcs[iclus] = ib->YCry();
1661 sigietaietapfbcs[iclus] = TMath::Sqrt(ib->CoviEtaiEta());
1662 sigiphiphipfbcs[iclus] = TMath::Sqrt(ib->CoviPhiiPhi());
1663 covietaiphipfbcs[iclus] = ib->CoviEtaiPhi();
1664 sigetaetapfbcs[iclus] = TMath::Sqrt(ib->CovEtaEta());
1665 sigphiphipfbcs[iclus] = TMath::Sqrt(ib->CovPhiPhi());
1666 covetaphipfbcs[iclus] = ib->CovEtaPhi();
1667 e3x3pfbcs[iclus] = ib->E3x3();
1668 e5x5pfbcs[iclus] = ib->E5x5();
1669 emaxpfbcs[iclus] = ib->EMax();
1670 e2ndpfbcs[iclus] = ib->E2nd();
1671 etoppfbcs[iclus] = ib->ETop();
1672 ebottompfbcs[iclus] = ib->EBottom();
1673 eleftpfbcs[iclus] = ib->ELeft();
1674 erightpfbcs[iclus] = ib->ERight();
1675 e1x3pfbcs[iclus] = ib->E1x3();
1676 e3x1pfbcs[iclus] = ib->E3x1();
1677 e1x5pfbcs[iclus] = ib->E1x5();
1678 e2x2pfbcs[iclus] = ib->E2x2();
1679 e4x4pfbcs[iclus] = ib->E4x4();
1680 e2x5maxpfbcs[iclus] = ib->E2x5Max();
1681 e2x5toppfbcs[iclus] = ib->E2x5Top();
1682 e2x5bottompfbcs[iclus] = ib->E2x5Bottom();
1683 e2x5leftpfbcs[iclus] = ib->E2x5Left();
1684 e2x5rightpfbcs[iclus] = ib->E2x5Right();
1685 nhitspfbcs[iclus]= ib->NHits();
1686 }
1687 else {
1688 epfbcs[iclus] = -999;
1689 etapfbcs[iclus] = -999;
1690 phipfbcs[iclus] = -999;
1691 ietapfbcs[iclus] = -999;
1692 iphipfbcs[iclus] = -999;
1693 ixpfbcs[iclus] = -999;
1694 iypfbcs[iclus] = -999;
1695 etacrypfbcs[iclus] = -999;
1696 phicrypfbcs[iclus] = -999;
1697 xcrypfbcs[iclus] = -999;
1698 ycrypfbcs[iclus] = -999;
1699 sigietaietapfbcs[iclus] = -999;
1700 sigiphiphipfbcs[iclus] = -999;
1701 covietaiphipfbcs[iclus] = -999;
1702 sigetaetapfbcs[iclus] = -999;
1703 sigphiphipfbcs[iclus] = -999;
1704 covetaphipfbcs[iclus] = -999;
1705 e3x3pfbcs[iclus] = -999;
1706 e5x5pfbcs[iclus] = -999;
1707 emaxpfbcs[iclus] = -999;
1708 e2ndpfbcs[iclus] = -999;
1709 etoppfbcs[iclus] = -999;
1710 ebottompfbcs[iclus] = -999;
1711 eleftpfbcs[iclus] = -999;
1712 erightpfbcs[iclus] = -999;
1713 e1x3pfbcs[iclus] = -999;
1714 e3x1pfbcs[iclus] = -999;
1715 e1x5pfbcs[iclus] = -999;
1716 e2x2pfbcs[iclus] = -999;
1717 e4x4pfbcs[iclus] = -999;
1718 e2x5maxpfbcs[iclus] = -999;
1719 e2x5toppfbcs[iclus] = -999;
1720 e2x5bottompfbcs[iclus] = -999;
1721 e2x5leftpfbcs[iclus] = -999;
1722 e2x5rightpfbcs[iclus] = -999;
1723 nhitspfbcs[iclus] = -999;
1724 }
1725 }
1726
1727 for (UInt_t iclus=0; iclus<100; ++iclus) {
1728 if (fillclusterarrays && pfsc && iclus < pfsc->NPsClusts() ) {
1729 const PsCluster *ib = pfsc->PsClust(iclus);
1730
1731 epsc[iclus] = ib->Energy();
1732 etapsc[iclus] = ib->Eta();
1733 phipsc[iclus] = ib->Phi();
1734 planepsc[iclus] = ib->PsPlane();
1735 }
1736 else {
1737 epsc[iclus] = -999;
1738 etapsc[iclus] = -999;
1739 phipsc[iclus] = -999;
1740 planepsc[iclus] = 0;
1741 }
1742 }
1743
1744 //initialize photon energy corrections if needed
1745 /*if (!PhotonFix::initialised()) {
1746 PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
1747 }*/
1748
1749 phfixph.setup(e,sceta,scphi,r9);
1750 phfixele.setup(e,sceta,scphi,r9);
1751
1752 const Float_t dval = -99.;
1753 ecor = phfixph.fixedEnergy();
1754 ecorerr = phfixph.sigmaEnergy();
1755 ecorele = phfixele.fixedEnergy();
1756 ecoreleerr = phfixele.sigmaEnergy();
1757 if (phfixph.isbarrel()) {
1758 etac = phfixph.etaC();
1759 etas = phfixph.etaS();
1760 etam = phfixph.etaM();
1761 phic = phfixph.phiC();
1762 phis = phfixph.phiS();
1763 phim = phfixph.phiM();
1764 xz = dval;
1765 xc = dval;
1766 xs = dval;
1767 xm = dval;
1768 yz = dval;
1769 yc = dval;
1770 ys = dval;
1771 ym = dval;
1772 }
1773 else {
1774 etac = dval;
1775 etas = dval;
1776 etam = dval;
1777 phic = dval;
1778 phis = dval;
1779 phim = dval;
1780 xz = phfixph.xZ();
1781 xc = phfixph.xC();
1782 xs = phfixph.xS();
1783 xm = phfixph.xM();
1784 yz = phfixph.yZ();
1785 yc = phfixph.yC();
1786 ys = phfixph.yS();
1787 ym = phfixph.yM();
1788 }
1789
1790 if (c) {
1791 hasconversion = kTRUE;
1792 convp = c->P();
1793 convpt = c->Pt();
1794 conveta = c->Eta();
1795 convphi = c->Phi();
1796 ThreeVector dirconvsc = ThreeVector(s->Point()) - c->Position();
1797 convdeta = c->Eta() - dirconvsc.Eta();
1798 convdphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
1799 convvtxrho = c->Position().Rho();
1800 convvtxz = c->Position().Z();
1801 convvtxphi = c->Position().Phi();
1802
1803 const StableData *leadsd = dynamic_cast<const StableData*>(c->DaughterDat(0));
1804 const StableData *trailsd = dynamic_cast<const StableData*>(c->DaughterDat(1));
1805 if (leadsd->Pt()<trailsd->Pt()) {
1806 const StableData *sdtmp = leadsd;
1807 leadsd = trailsd;
1808 trailsd = sdtmp;
1809 }
1810
1811 const Track *leadtrack = dynamic_cast<const StableParticle*>(leadsd->Original())->Trk();
1812 const Track *trailtrack = dynamic_cast<const StableParticle*>(trailsd->Original())->Trk();
1813
1814 convleadpt = leadsd->Pt();
1815 convtrailpt = trailsd->Pt();
1816 convleadtrackpt = leadtrack->Pt();
1817 convleadtrackalgo = leadtrack->Algo();
1818 if (convleadtrackalgo==29) convleadtrackalgos=2; //gsf track
1819 else if (convleadtrackalgo==15 ||convleadtrackalgo==16) convleadtrackalgos=1; //ecal-seeded track
1820 else convleadtrackalgos = 0; //std iterative track
1821 convleadtrackcharge = leadtrack->Charge();
1822 convtrailtrackpt = trailtrack->Pt();
1823 convtrailtrackalgo = trailtrack->Algo();
1824 if (convtrailtrackalgo==29) convtrailtrackalgos=2; //gsf track
1825 else if (convtrailtrackalgo==15 ||convtrailtrackalgo==16) convtrailtrackalgos=1; //ecal-seeded track
1826 else convtrailtrackalgos = 0; //std iterative track
1827 trailtrackcharge = trailtrack->Charge();
1828 }
1829 else {
1830 hasconversion = kFALSE;
1831 convp = -99.;
1832 convpt = -99.;
1833 conveta = -99.;
1834 convphi = -99.;
1835 convdeta = -99.;
1836 convdphi = -99.;
1837 convvtxrho = -99.;
1838 convvtxz = -999.;
1839 convvtxphi = -99.;
1840 convleadpt = -99.;
1841 convtrailpt = -99.;
1842 convleadtrackpt = -99.;
1843 convleadtrackalgo = -99;
1844 convleadtrackalgos = -99;
1845 convleadtrackcharge = 0;
1846 convtrailtrackpt = -99.;
1847 convtrailtrackalgo = -99;
1848 convtrailtrackalgos = -99;
1849 trailtrackcharge = 0;
1850 }
1851
1852 //electron quantities
1853 if (ele) {
1854 haselectron = kTRUE;
1855 eleisecaldriven = ele->IsEcalDriven();
1856 eleistrackerdriven = ele->IsTrackerDriven();
1857 elee = ele->E();
1858 elept = ele->Pt();
1859 eleeta = ele->Eta();
1860 elephi = ele->Phi();
1861 elecharge = ele->Charge();
1862 elefbrem = ele->FBrem();
1863 eledeta = ele->DeltaEtaSuperClusterTrackAtVtx();
1864 eledphi = ele->DeltaPhiSuperClusterTrackAtVtx();
1865 elep = s->Energy()/ele->ESuperClusterOverP();
1866 elepin = ele->PIn();
1867 elepout = ele->POut();
1868 }
1869 else {
1870 haselectron = kFALSE;
1871 eleisecaldriven = kFALSE;
1872 eleistrackerdriven = kFALSE;
1873 elee = -99.;
1874 elept = -99.;
1875 eleeta = -99.;
1876 elephi = -99.;
1877 elecharge = -99;
1878 elefbrem = -99.;
1879 eledeta = -99.;
1880 eledphi = -99.;
1881 elep = -99.;
1882 elepin = -99.;
1883 elepout = -99.;
1884 }
1885
1886 //pf supercluster quantities
1887 if (pfsc) {
1888 haspfsc = kTRUE;
1889 pfsce = pfsc->Energy();
1890 pfscrawe = pfsc->RawEnergy();
1891 pfsceta = pfsc->Eta();
1892 pfscphi = pfsc->Phi();
1893 pfscnclusters = pfsc->NClusters();
1894 pfscnhits = pfsc->NHits();
1895 pfscetawidth = pfsc->EtaWidth();
1896 pfscphiwidth = pfsc->PhiWidth();
1897 pfscnpsclusters = pfsc->NPsClusts();
1898 }
1899 else {
1900 haspfsc = kFALSE;
1901 pfsce = -99.;
1902 pfscrawe = -99.;
1903 pfsceta = -99.;
1904 pfscphi = -99.;
1905 pfscnclusters = 0;
1906 pfscnhits = 0;
1907 pfscetawidth = -99.;
1908 pfscphiwidth = -99.;
1909 pfscnpsclusters = 0;
1910 }
1911
1912 genz = -99.;
1913 if (m) {
1914 ispromptgen = kTRUE;
1915 gene = m->E();
1916 genpt = m->Pt();
1917 geneta = m->Eta();
1918 genphi = m->Phi();
1919 const MCParticle *mm = m->DistinctMother();
1920 if (mm) genz = mm->DecayVertex().Z();
1921 pdgid = m->PdgId();
1922 if (mm) motherpdgid = mm->PdgId();
1923 else motherpdgid = -99;
1924 }
1925 else {
1926 ispromptgen = kFALSE;
1927 gene = -99.;
1928 genpt = -99.;
1929 geneta = -99.;
1930 genphi = -99.;
1931 pdgid = -99;
1932 motherpdgid = -99;
1933 }
1934 }
1935
1936 void PhotonTreeWriterVtx::SetVars(const Vertex *v, const Photon *p1, const Photon *p2, const PFCandidateCol *pfcands, Int_t idx, Int_t numvtx, Float_t genvtxz)
1937 {
1938
1939 //printf("start\n");
1940
1941 n = idx;
1942 nvtx = numvtx;
1943 zgen = genvtxz;
1944
1945 x = v->X();
1946 y = v->Y();
1947 z = v->Z();
1948
1949 Double_t dsumpt = 0.;
1950 Double_t dsumptsq = 0.;
1951
1952 nchalltoward = 0;
1953 nchalltransverse = 0;
1954 nchallaway = 0;
1955 nchcuttoward = 0;
1956 nchcuttransverse = 0;
1957 nchcutaway = 0;
1958
1959 ThreeVector vtxmom;
1960
1961 //printf("mom\n");
1962 FourVectorM diphoton = p1->MomVtx(v->Position()) + p2->MomVtx(v->Position());
1963 //printf("done mom\n");
1964 ptgg = diphoton.Pt();
1965 phigg = diphoton.Phi();
1966 etagg = diphoton.Eta();
1967 mgg = diphoton.M();
1968 pxgg = diphoton.Px();
1969 pygg = diphoton.Py();
1970 pzgg = diphoton.Pz();
1971
1972 //printf("loop\n");
1973
1974 for (UInt_t i = 0; i<pfcands->GetEntries(); ++i) {
1975 const PFCandidate *pfc = pfcands->At(i);
1976 if (pfc->PFType()!=PFCandidate::eHadron || !pfc->HasTrackerTrk()) continue;
1977 if (TMath::Abs( pfc->TrackerTrk()->DzCorrected(*v) ) > 0.2) continue;
1978 if (TMath::Abs( pfc->TrackerTrk()->D0Corrected(*v) ) > 0.1) continue;
1979
1980 vtxmom += ThreeVector(pfc->Px(),pfc->Py(),pfc->Pz());
1981
1982 dsumpt += pfc->Pt();
1983 dsumptsq += pfc->Pt()*pfc->Pt();
1984
1985 Double_t dphi = TMath::Abs(MathUtils::DeltaPhi(*pfc,diphoton));
1986 if (dphi<(TMath::Pi()/3.0)) {
1987 ++nchalltoward;
1988 if (pfc->Pt()>0.5) ++nchcuttoward;
1989 }
1990 else if (dphi>(2.0*TMath::Pi()/3.0)) {
1991 ++nchallaway;
1992 if (pfc->Pt()>0.5) ++nchcutaway;
1993 }
1994 else {
1995 ++nchalltransverse;
1996 if (pfc->Pt()>0.5) ++nchcuttransverse;
1997 }
1998
1999 }
2000
2001 //printf("doneloop\n");
2002
2003
2004 sumpt = dsumpt;
2005 sumptsq = dsumptsq;
2006
2007 pt = vtxmom.Rho();
2008 phi = vtxmom.Phi();
2009 eta = vtxmom.Eta();
2010 px = vtxmom.X();
2011 py = vtxmom.Y();
2012 pz = vtxmom.Z();
2013
2014 //printf("done\n");
2015
2016 return;
2017
2018 }