ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.41
Committed: Wed Jan 16 23:27:21 2013 UTC (12 years, 3 months ago) by mingyang
Content type: text/plain
Branch: MAIN
Changes since 1.40: +51 -34 lines
Log Message: 
addtional eleveto add for photon for electron tag

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