ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.73
Committed: Tue Dec 10 15:31:21 2013 UTC (11 years, 5 months ago) by bendavid
Content type: text/plain
Branch: MAIN
Changes since 1.72: +63 -35 lines
Log Message: 
various lepton tagging fixes

File Contents

# Content
1 #include <time.h>
2 #include "MitPhysics/Mods/interface/PhotonTreeWriter.h"
3 #include "MitAna/DataTree/interface/PhotonCol.h"
4 #include "MitAna/DataTree/interface/PFCandidateCol.h"
5 #include "MitAna/DataTree/interface/StableData.h"
6 #include "MitAna/DataTree/interface/StableParticle.h"
7 #include "MitAna/DataTree/interface/PFMet.h"
8 #include "MitPhysics/Init/interface/ModNames.h"
9 #include "MitPhysics/Utils/interface/IsolationTools.h"
10 #include "MitPhysics/Utils/interface/PhotonTools.h"
11 #include "MitPhysics/Utils/interface/VertexTools.h"
12 #include "MitPhysics/Utils/interface/PFMetCorrectionTools.h"
13 #include "MitPhysics/Utils/interface/JetTools.h"
14 #include "MitAna/DataTree/interface/PFJetCol.h"
15 #include "MitAna/DataTree/interface/GenJetCol.h"
16 #include "TDataMember.h"
17 #include "TLorentzVector.h"
18 #include "TFile.h"
19 #include <TNtuple.h>
20 #include <TRandom3.h>
21 #include <TSystem.h>
22
23 using namespace mithep;
24
25 ClassImp(mithep::PhotonTreeWriter)
26 templateClassImp(mithep::PhotonTreeWriterPhoton)//ming: what's this?
27 ClassImp(mithep::PhotonTreeWriterDiphotonEvent)
28
29 //--------------------------------------------------------------------------------------------------
30 PhotonTreeWriter::PhotonTreeWriter(const char *name, const char *title) :
31 // Base Module...
32 BaseMod (name,title),
33 // define all the Branches to load
34 fPhotonBranchName (Names::gkPhotonBrn),
35 fPFPhotonName ("PFPhotons"),
36 fElectronName (Names::gkElectronBrn),
37 fGoodElectronName (Names::gkElectronBrn),
38 fConversionName (Names::gkMvfConversionBrn),
39 fPFConversionName ("PFPhotonConversions"),
40 fTrackBranchName (Names::gkTrackBrn),
41 fPileUpDenName (Names::gkPileupEnergyDensityBrn),
42 fPVName (Names::gkPVBeamSpotBrn),
43 fBeamspotName (Names::gkBeamSpotBrn),
44 fPFCandName (Names::gkPFCandidatesBrn),
45 fPFNoPileUpName ("PFNoPileUp"),
46 fPFPileUpName ("PFPileUp"),
47 fMCParticleName (Names::gkMCPartBrn),
48 fMCEventInfoName (Names::gkMCEvtInfoBrn),
49 fPileUpName (Names::gkPileupInfoBrn),
50 fPFSuperClusterName ("PFSuperClusters"),
51 fPFMetName ("PFMet"),
52 fPFJetName (Names::gkPFJetBrn),
53 funcorrPFJetName ("AKt5PFJets"),
54 fGenJetName ("AKT5GenJets"),
55 fLeptonTagElectronsName ("HggLeptonTagElectrons"),
56 fLeptonTagMuonsName ("HggLeptonTagMuons"),
57 fLeptonTagSoftElectronsName ("HggLeptonTagSoftElectrons"),
58 fLeptonTagSoftMuonsName ("HggLeptonTagSoftMuons"),
59
60 fIsData (false),
61 fIsCutBased (false),
62 fPhotonsFromBranch (kTRUE),
63 fPVFromBranch (kTRUE),
64 fGoodElectronsFromBranch(kTRUE),
65 fPFJetsFromBranch (kTRUE),
66 // ----------------------------------------
67 // flag for synchronization, adds vertex variables
68 // should be on for synching trees
69 fDoSynching (kFALSE),
70
71 // ----------------------------------------
72 // collections....
73 fPhotons (0),
74 fPFPhotons (0),
75 fElectrons (0),
76 fConversions (0),
77 fPFConversions (0),
78 fTracks (0),
79 fPileUpDen (0),
80 fPV (0),
81 fBeamspot (0),
82 fPFCands (0),
83 fMCParticles (0),
84 fMCEventInfo (0),
85 fPileUp (0),
86 fPFSuperClusters (0),
87 fPFJets (0),
88 fGenJets (0),
89 funcorrPFJets (0),
90
91 fLeptonTagElectrons (0),
92 fLeptonTagMuons (0),
93 fPFNoPileUpCands (0),
94 fPFPileUpCands (0),
95
96 fLoopOnGoodElectrons (kFALSE),
97 fApplyElectronVeto (kTRUE),
98 fWriteDiphotonTree (kTRUE),
99 fWriteSingleTree (kTRUE),
100 fEnablePFPhotons (kTRUE),
101 fExcludeSinglePrompt (kFALSE),
102 fExcludeDoublePrompt (kFALSE),
103 fEnableJets (kFALSE),
104 fEnableGenJets (kFALSE),
105 fApplyJetId (kFALSE),
106 fApplyLeptonTag (kFALSE),
107 fApplyVHLepTag (kFALSE),
108 fApplyVHHadTag (kFALSE),
109 fApplyVBFTag (kFALSE),
110 fApplyTTHTag (kFALSE),
111 fApplyBTag (kFALSE),
112 fApplyPFMetCorrections (kFALSE),
113 fFillClusterArrays (kFALSE),
114 fFillVertexTree (kFALSE),
115 fDo2012LepTag (kFALSE),
116 fVerbosityLevel (0),
117 fPhFixDataFile (gSystem->Getenv("CMSSW_BASE") +
118 TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat")),
119 fBeamspotWidth (5.8),
120 fTmpFile (0),
121
122 // JV: moved up the initializtion of fTupleName to avoid compilation warning
123 fTupleName ("hPhotonTree"),
124
125 fElectronIDMVA(0),
126 fElectronMVAWeights_Subdet0Pt10To20(""),
127 fElectronMVAWeights_Subdet1Pt10To20(""),
128 fElectronMVAWeights_Subdet2Pt10To20(""),
129 fElectronMVAWeights_Subdet0Pt20ToInf(""),
130 fElectronMVAWeights_Subdet1Pt20ToInf(""),
131 fElectronMVAWeights_Subdet2Pt20ToInf(""),
132
133 fTheRhoType(RhoUtilities::DEFAULT),
134 fProcessedEvents(0)
135
136 {
137 // Constructor
138 }
139
140 PhotonTreeWriter::~PhotonTreeWriter()
141 {
142 // Destructor
143 // Deal with the temporary file here?
144 // fTmpFile->Write();
145 fTmpFile->Close();
146 // TString shellcmd = TString("rm ") + TString(fTmpFile->GetName());
147 // delete fTmpFile;
148 // cout << shellcmd.Data() << endl;
149 // gSystem->Exec(shellcmd.Data());
150
151 }
152
153 //--------------------------------------------------------------------------------------------------
154 void PhotonTreeWriter::SlaveTerminate()
155 {
156
157 if (hCiCTuple) fTmpFile->WriteTObject(hCiCTuple,hCiCTuple->GetName());
158 if (hCiCTupleSingle) fTmpFile->WriteTObject(hCiCTuple,hCiCTupleSingle->GetName());
159
160 }
161
162
163 //--------------------------------------------------------------------------------------------------
164 void PhotonTreeWriter::Process()
165 {
166
167 fProcessedEvents++;
168
169 if (fVerbosityLevel > 0) {
170 PhotonTreeWriter::LogEventInfo();
171 }
172
173 //if(GetEventHeader()->EvtNum()==9008 || GetEventHeader()->EvtNum()==9008 || GetEventHeader()->EvtNum()==9010){
174 // printf("check photontreewriter 0\n");
175 // }
176 // ------------------------------------------------------------
177 // Process entries of the tree.
178 LoadEventObject(fPhotonBranchName, fPhotons);
179 LoadEventObject(fGoodElectronName, fGoodElectrons);
180
181 // lepton tag collections
182 if( fApplyLeptonTag || fApplyVHLepTag ) {
183 LoadEventObject(fLeptonTagElectronsName, fLeptonTagElectrons);
184 LoadEventObject(fLeptonTagMuonsName, fLeptonTagMuons);
185 }
186
187 if( fApplyVHLepTag ) {
188 LoadEventObject(fLeptonTagSoftElectronsName, fLeptonTagSoftElectrons);
189 LoadEventObject(fLeptonTagSoftMuonsName, fLeptonTagSoftMuons);
190 }
191
192 const BaseCollection *egcol = 0;
193 if (fLoopOnGoodElectrons)
194 egcol = fGoodElectrons;
195 else
196 egcol = fPhotons;
197 if (egcol->GetEntries()<1)
198 return;
199
200 if (fEnablePFPhotons) LoadEventObject(fPFPhotonName, fPFPhotons);
201 LoadEventObject(fElectronName, fElectrons);
202 LoadEventObject(fConversionName, fConversions);
203 if ( fDoSynching ) LoadEventObject(fPFConversionName, fPFConversions);
204 LoadEventObject(fTrackBranchName, fTracks);
205 LoadEventObject(fPileUpDenName, fPileUpDen);
206 LoadEventObject(fPVName, fPV);
207 LoadEventObject(fBeamspotName, fBeamspot);
208 LoadEventObject(fPFCandName, fPFCands);
209 LoadEventObject(fPFSuperClusterName, fPFSuperClusters);
210 LoadEventObject(fPFMetName, fPFMet);
211 // LoadEventObject(fPFNoPileUpName, fPFNoPileUpCands);
212 // LoadEventObject(fPFPileUpName, fPFPileUpCands);
213
214 if (fEnableJets){
215 LoadEventObject(fPFJetName, fPFJets);
216 //LoadEventObject(funcorrPFJetName, funcorrPFJets);
217 LoadBranch(funcorrPFJetName);
218 // if(!fIsData) LoadEventObject(fGenJetName, fGenJets);
219 }
220
221
222 // ------------------------------------------------------------
223 // load event based information
224 Int_t _numPU = -1.; // some sensible default values....
225 Int_t _numPUminus = -1.; // some sensible default values....
226 Int_t _numPUplus = -1.; // some sensible default values....
227
228 Float_t rho = -99.;
229 if( fPileUpDen->GetEntries() > 0 )
230 rho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();//m
231
232 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));//ming:what's this?
233
234 if( !fIsData ) {
235 LoadBranch(fMCParticleName);
236 LoadBranch(fMCEventInfoName);
237 LoadBranch(fPileUpName);
238 if (fEnableGenJets) LoadEventObject(fGenJetName, fGenJets);
239 } else fGenJets = NULL;
240
241 if( !fIsData ) {
242 for (UInt_t i=0; i<fPileUp->GetEntries(); ++i) {
243 const PileupInfo *puinfo = fPileUp->At(i);
244 if (puinfo->GetBunchCrossing()==0) _numPU = puinfo->GetPU_NumInteractions();
245 else if (puinfo->GetBunchCrossing() == -1) _numPUminus = puinfo->GetPU_NumInteractions();
246 else if (puinfo->GetBunchCrossing() == 1) _numPUplus = puinfo->GetPU_NumInteractions();
247 }
248 }
249
250 // in case of a MC event, try to find Higgs and Higgs decay Z poisition
251 Float_t _pth = -100.;
252 Float_t _decayZ = -100.;
253 Float_t _genmass = -100.;
254 if (!fIsData)
255 FindHiggsPtAndZ(_pth, _decayZ, _genmass);
256
257
258 fDiphotonEvent->genz = -999.;
259 if (!fIsData) {
260 fDiphotonEvent->genz = fMCParticles->At(0)->DecayVertex().Z();
261 }
262
263 fDiphotonEvent->mcprocid = -999;
264 if (!fIsData) {
265 fDiphotonEvent->mcprocid = fMCEventInfo->ProcessId();
266 fDiphotonEvent->mcweight = fMCEventInfo->Weight();
267 }
268
269 Double_t _evt = GetEventHeader()->EvtNum();
270
271 Double_t _spfMet = fPFMet->At(0)->SumEt();
272
273 fDiphotonEvent->leptonTag = -1; // disabled
274 fDiphotonEvent->VHLepTag = -1; // disabled
275 fDiphotonEvent->VHHadTag = -1; // disabled
276
277 // ====================================================
278 // Vtx synching stuff...
279 fDiphotonEvent->vtxInd1 = -1;
280 fDiphotonEvent->vtxInd2 = -1;
281 fDiphotonEvent->vtxInd3 = -1;
282
283 fDiphotonEvent->vtxBestPtbal = -1.;
284 fDiphotonEvent->vtxBestPtasym = -1.;
285 fDiphotonEvent->vtxBestSumpt2 = -1.;
286 fDiphotonEvent->vtxBestP2Conv = -1.;
287
288 fDiphotonEvent->vtxMva1 = -1.;
289 fDiphotonEvent->vtxMva2 = -1.;
290 fDiphotonEvent->vtxMva3 = -1.;
291
292 fDiphotonEvent->vtxNleg1 = -1;
293 fDiphotonEvent->vtxNleg2 = -1;
294 fDiphotonEvent->vtxNconv = -1;
295 // ====================================================
296
297
298 fDiphotonEvent->rho = fPileUpDen->At(0)->RhoKt6PFJets();
299 fDiphotonEvent->rho25 = fPileUpDen->At(0)->RhoRandomLowEta();
300 fDiphotonEvent->rhoold = fPileUpDen->At(0)->Rho();
301 fDiphotonEvent->genHiggspt = _pth;
302 fDiphotonEvent->genHiggsZ = _decayZ;
303 fDiphotonEvent->genmass = _genmass;
304 fDiphotonEvent->gencostheta = -99.;
305 fDiphotonEvent->nVtx = fPV->GetEntries();
306 fDiphotonEvent->bsX = fBeamspot->At(0)->X();
307 fDiphotonEvent->bsY = fBeamspot->At(0)->Y();
308 fDiphotonEvent->bsZ = fBeamspot->At(0)->Z();
309 fDiphotonEvent->bsSigmaZ = fBeamspot->At(0)->SigmaZ();
310 fDiphotonEvent->vtxX = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->X() : -99.;
311 fDiphotonEvent->vtxY = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Y() : -99.;
312 fDiphotonEvent->vtxZ = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Z() : -99.;
313 fDiphotonEvent->numPU = _numPU;
314 fDiphotonEvent->numPUminus = _numPUminus;
315 fDiphotonEvent->numPUplus = _numPUplus;
316 fDiphotonEvent->mass = -99.;
317 fDiphotonEvent->ptgg = -99.;
318 fDiphotonEvent->costheta = -99.;
319 fDiphotonEvent->mt = -99.;
320 fDiphotonEvent->cosphimet = -99.;
321 fDiphotonEvent->mtele = -99.;
322 fDiphotonEvent->cosphimetele = -99.;
323 fDiphotonEvent->evt = GetEventHeader()->EvtNum();
324 fDiphotonEvent->run = GetEventHeader()->RunNum();
325 fDiphotonEvent->lumi = GetEventHeader()->LumiSec();
326 fDiphotonEvent->evtcat = -99;
327 fDiphotonEvent->nobj = fPhotons->GetEntries();
328 fDiphotonEvent->pfmet = fPFMet->At(0)->Pt();
329 fDiphotonEvent->pfmetphi = fPFMet->At(0)->Phi();
330 fDiphotonEvent->pfmetx = fPFMet->At(0)->Px();
331 fDiphotonEvent->pfmety = fPFMet->At(0)->Py();
332 fDiphotonEvent->spfMet = _spfMet;
333 fDiphotonEvent->masscor = -99.;
334 fDiphotonEvent->masscorerr = -99.;
335 fDiphotonEvent->masscorele = -99.;
336 fDiphotonEvent->masscoreleerr = -99.;
337 fDiphotonEvent->ismc = GetEventHeader()->IsMC();
338
339 //jets
340 const Jet *jet1 = 0;
341 const Jet *jet2 = 0;
342 const Jet *jetcentral = 0;
343
344 fDiphotonEvent->jetleadNoIDpt = -99.;
345 fDiphotonEvent->jetleadNoIDeta = -99.;
346 fDiphotonEvent->jetleadNoIDphi = -99.;
347 fDiphotonEvent->jetleadNoIDmass = -99.;
348 fDiphotonEvent->jet1pt = -99.;
349 fDiphotonEvent->jet1eta = -99.;
350 fDiphotonEvent->jet1phi = -99.;
351 fDiphotonEvent->jet1mass = -99.;
352 fDiphotonEvent->jet2pt = -99.;
353 fDiphotonEvent->jet2eta = -99.;
354 fDiphotonEvent->jet2phi = -99.;
355 fDiphotonEvent->jet2mass = -99.;
356 fDiphotonEvent->jetcentralpt = -99.;
357 fDiphotonEvent->jetcentraleta = -99.;
358 fDiphotonEvent->jetcentralphi = -99.;
359 fDiphotonEvent->jetcentralmass = -99.;
360 fDiphotonEvent->dijetpt = -99.;
361 fDiphotonEvent->dijeteta = -99.;
362 fDiphotonEvent->dijetphi = -99.;
363 fDiphotonEvent->dijetmass = -99.;
364 fDiphotonEvent->jetetaplus = -99.;
365 fDiphotonEvent->jetetaminus = -99.;
366 fDiphotonEvent->zeppenfeld = -99.;
367 fDiphotonEvent->dphidijetgg = -99.;
368
369 //uncorrected jets
370 // const Jet *uncorrjet1 = 0;
371 // const Jet *uncorrjet2 = 0;
372 // const Jet *uncorrjetcentral = 0;
373
374 // fDiphotonEvent->uncorrjet1pt = -99.;
375 // fDiphotonEvent->uncorrjet1eta = -99.;
376 // fDiphotonEvent->uncorrjet1phi = -99.;
377 // fDiphotonEvent->uncorrjet1mass = -99.;
378 // fDiphotonEvent->uncorrjet2pt = -99.;
379 // fDiphotonEvent->uncorrjet2eta = -99.;
380 // fDiphotonEvent->uncorrjet2phi = -99.;
381 // fDiphotonEvent->uncorrjet2mass = -99.;
382 // fDiphotonEvent->uncorrjetcentralpt = -99.;
383 // fDiphotonEvent->uncorrjetcentraleta = -99.;
384 // fDiphotonEvent->uncorrjetcentralphi = -99.;
385 // fDiphotonEvent->uncorrjetcentralmass = -99.;
386 // fDiphotonEvent->diuncorrjetpt = -99.;
387 // fDiphotonEvent->diuncorrjeteta = -99.;
388 // fDiphotonEvent->diuncorrjetphi = -99.;
389 // fDiphotonEvent->diuncorrjetmass = -99.;
390
391
392
393 Int_t nhitsbeforevtxmax = 1;
394 if (!fApplyElectronVeto)
395 nhitsbeforevtxmax = 999;
396
397 if (egcol->GetEntries()>=2) {
398 const Particle *p1 = 0;
399 const Particle *p2 = 0;
400 const Photon *phHard = 0;
401 const Photon *phSoft = 0;
402 const Electron *ele1 = 0;
403 const Electron *ele2 = 0;
404 const SuperCluster *sc1 = 0;
405 const SuperCluster *sc2 = 0;
406
407 if (fLoopOnGoodElectrons) {
408 ele1 = fGoodElectrons->At(0);
409 ele2 = fGoodElectrons->At(1);
410 p1 = ele1;
411 p2 = ele2;
412 sc1 = ele1->SCluster();
413 sc2 = ele2->SCluster();
414 phHard = PhotonTools::MatchedPhoton(ele1,fPhotons);
415 phSoft = PhotonTools::MatchedPhoton(ele2,fPhotons);
416 }
417 else {
418 phHard = fPhotons->At(0);
419 phSoft = fPhotons->At(1);
420 p1 = phHard;
421 p2 = phSoft;
422 sc1 = phHard->SCluster();
423 sc2 = phSoft->SCluster();
424 ele1 = PhotonTools::MatchedElectron(phHard,fGoodElectrons);
425 ele2 = PhotonTools::MatchedElectron(phSoft,fGoodElectrons);
426 }
427
428 const DecayParticle *conv1 = PhotonTools::MatchedConversion(sc1,fConversions,bsp,
429 nhitsbeforevtxmax);
430 const DecayParticle *conv2 = PhotonTools::MatchedConversion(sc2,fConversions,bsp,
431 nhitsbeforevtxmax);
432
433 const SuperCluster *pfsc1 = PhotonTools::MatchedPFSC(sc1,fPFPhotons, fElectrons);
434 const SuperCluster *pfsc2 = PhotonTools::MatchedPFSC(sc2,fPFPhotons, fElectrons);
435
436 const MCParticle *phgen1 = 0;
437 const MCParticle *phgen2 = 0;
438 if (!fIsData) {
439 phgen1 = PhotonTools::MatchMC(p1,fMCParticles,!fApplyElectronVeto);
440 phgen2 = PhotonTools::MatchMC(p2,fMCParticles,!fApplyElectronVeto);
441 }
442
443 if (fExcludeSinglePrompt && (phgen1 || phgen2))
444 return;
445 if (fExcludeDoublePrompt && (phgen1 && phgen2))
446 return;
447
448 if (!fLoopOnGoodElectrons && phHard->HasPV()) {
449 fDiphotonEvent->vtxX = phHard->PV()->X();
450 fDiphotonEvent->vtxY = phHard->PV()->Y();
451 fDiphotonEvent->vtxZ = phHard->PV()->Z();
452 fDiphotonEvent->vtxprob = phHard->VtxProb();
453 }
454
455 // fill Btag information... set to true if One jet fullfills
456 fDiphotonEvent -> btagJet1 = -1.;
457 fDiphotonEvent -> btagJet1Pt = -99.;
458 fDiphotonEvent -> btagJet1Eta = -99.;
459
460 fDiphotonEvent -> btagJet2 = -1.;
461 fDiphotonEvent -> btagJet2Pt = -99.;
462 fDiphotonEvent -> btagJet2Eta = -99.;
463
464 if( fApplyBTag && fEnableJets ) {
465 float highTag = 0.;
466 float highJetPt = 0.;
467 float highJetEta = 0.;
468
469 float trailTag = 0.;
470 float trailJetPt = 0.;
471 float trailJetEta = 0.;
472
473 for (UInt_t ijet=0; ijet<fPFJets->GetEntries();++ijet) {
474 const Jet *jet = fPFJets->At(ijet);
475 if ( jet->Pt() < 20. || jet->AbsEta() > 2.4 ) continue;
476 if ( jet->CombinedSecondaryVertexBJetTagsDisc() > highTag ) {
477
478 trailTag = highTag;
479 trailJetPt = highJetPt;
480 trailJetEta = highJetEta;
481
482 highTag = jet->CombinedSecondaryVertexBJetTagsDisc();
483 highJetPt = jet->Pt();
484 highJetEta = jet->Eta();
485
486 } else if ( jet->CombinedSecondaryVertexBJetTagsDisc() > trailTag ) {
487
488 trailTag = jet->CombinedSecondaryVertexBJetTagsDisc();
489 trailJetPt = jet->Pt();
490 trailJetEta = jet->Eta();
491 }
492 }
493 fDiphotonEvent -> btagJet1 = highTag;
494 fDiphotonEvent -> btagJet1Pt = highJetPt;
495 fDiphotonEvent -> btagJet1Eta = highJetEta;
496
497 fDiphotonEvent -> btagJet2 = trailTag;
498 fDiphotonEvent -> btagJet2Pt = trailJetPt;
499 fDiphotonEvent -> btagJet2Eta = trailJetEta;
500
501 }
502
503
504 //fill jet variables
505 const Vertex *selvtx = fPV->At(0);
506 if (!fLoopOnGoodElectrons && phHard->HasPV()) selvtx = phHard->PV();
507 if (fEnableJets) {
508 for (UInt_t ijet=0; ijet<fPFJets->GetEntries();++ijet) {
509 const Jet *jet = fPFJets->At(ijet);
510 if (jet->AbsEta()<4.7 && MathUtils::DeltaR(jet,p1)>0.5 && MathUtils::DeltaR(jet,p2)>0.5) {
511 const PFJet *pfjet = dynamic_cast<const PFJet*>(jet);
512 if (!pfjet) continue;
513 if (fApplyJetId && !fJetId.passCut(pfjet,selvtx,fPV)) continue;
514 if (!jet1) jet1 = jet;
515 else if (!jet2) jet2 = jet;
516 else if (!jetcentral && 0) jetcentral = jet;
517 }
518 if (jet1&&jet2&&jetcentral) break;
519 }
520
521 for(UInt_t ijet=0; ijet<fPFJets->GetEntries();++ijet){
522
523 const Jet *jet = fPFJets->At(ijet);
524 if (MathUtils::DeltaR(jet,p1)>0.5 && MathUtils::DeltaR(jet,p2)>0.5) {
525 fDiphotonEvent->jetleadNoIDpt = jet->Pt();
526 fDiphotonEvent->jetleadNoIDeta = jet->Eta();
527 fDiphotonEvent->jetleadNoIDphi = jet->Phi();
528 fDiphotonEvent->jetleadNoIDmass = jet->Mass();
529 break;
530 }
531 }
532
533 if (jet1) {
534 fDiphotonEvent->jet1pt = jet1->Pt();
535 fDiphotonEvent->jet1eta = jet1->Eta();
536 fDiphotonEvent->jet1phi = jet1->Phi();
537 fDiphotonEvent->jet1mass = jet1->Mass();
538 }
539
540 if (jet2) {
541 fDiphotonEvent->jet2pt = jet2->Pt();
542 fDiphotonEvent->jet2eta = jet2->Eta();
543 fDiphotonEvent->jet2phi = jet2->Phi();
544 fDiphotonEvent->jet2mass = jet2->Mass();
545 }
546
547 if (jetcentral) {
548 fDiphotonEvent->jetcentralpt = jetcentral->Pt();
549 fDiphotonEvent->jetcentraleta = jetcentral->Eta();
550 fDiphotonEvent->jetcentralphi = jetcentral->Phi();
551 fDiphotonEvent->jetcentralmass = jetcentral->Mass();
552 }
553
554 if (jet1&&jet2){
555 FourVectorM momjj = (jet1->Mom() + jet2->Mom());
556
557 fDiphotonEvent->dijetpt = momjj.Pt();
558 fDiphotonEvent->dijeteta = momjj.Eta();
559 fDiphotonEvent->dijetphi = momjj.Phi();
560 fDiphotonEvent->dijetmass = momjj.M();
561
562 if (jet1->Eta()>jet2->Eta()) {
563 fDiphotonEvent->jetetaplus = jet1->Eta();
564 fDiphotonEvent->jetetaminus = jet2->Eta();
565 }
566 else {
567 fDiphotonEvent->jetetaplus = jet2->Eta();
568 fDiphotonEvent->jetetaminus = jet1->Eta();
569 }
570 }
571 }
572
573 //added gen. info of whether a lep. or nutrino is from W or Z --Heng 02/14/2012 12:30 EST
574 Double_t _fromZ = -99;
575 Double_t _fromW = -99;
576 Float_t _zpt = -99;
577 Float_t _zEta = -99;
578 Float_t _allZpt = -99;
579 Float_t _allZEta = -99;
580
581 if( !fIsData ){
582
583 // loop over all GEN particles and look for nutrinos whoes mother is Z
584 for(UInt_t j=0; j<fMCParticles->GetEntries(); ++j) {
585 const MCParticle* p = fMCParticles->At(j);
586 if( p->AbsPdgId()==23 ||p->AbsPdgId()==32 || p->AbsPdgId()==33 ) {
587 _allZpt=p->Pt();
588 _allZEta=p->Eta();
589 if (p->HasDaughter(12,kFALSE) || p->HasDaughter(14,kFALSE) || p->HasDaughter(16,kFALSE) ||p->HasDaughter(18,kFALSE) ) {
590 _fromZ=1;
591 _zpt=p->Pt();
592 _zEta=p->Eta();
593 }
594 }
595 else _fromW=1;
596 }
597 }
598
599 fDiphotonEvent->fromZ = _fromZ;
600 fDiphotonEvent->fromW = _fromW;
601 fDiphotonEvent->zpt = _zpt;
602 fDiphotonEvent->zEta = _zEta;
603 fDiphotonEvent->allZpt = _allZpt;
604 fDiphotonEvent->allZEta = _allZEta;
605
606 Double_t _dphiMetgg = -99;
607 Double_t _cosdphiMetgg = -99;
608 Double_t _dphiPhPh = -99;
609
610 Double_t _mass = -99.;
611 Double_t _masserr = -99.;
612 Double_t _masserrsmeared = -99.;
613 Double_t _masserrwrongvtx = -99.;
614 Double_t _masserrsmearedwrongvtx = -99.;
615 Double_t _ptgg = -99.;
616 Double_t _etagg = -99.;
617 Double_t _phigg = -99.;
618 Double_t _costheta = -99.;
619 PhotonTools::DiphotonR9EtaPtCats _evtcat = PhotonTools::kOctCat0;
620
621 const Vertex* realVtx = NULL;
622
623 if (phHard && phSoft) {
624 _dphiMetgg = MathUtils::DeltaPhi((phHard->Mom()+phSoft->Mom()).Phi(),fPFMet->At(0)->Phi());
625 _cosdphiMetgg = TMath::Cos(_dphiMetgg);
626 _dphiPhPh = MathUtils::DeltaPhi((phHard->Mom()).Phi(),(phSoft->Mom()).Phi());
627 _mass = (phHard->Mom()+phSoft->Mom()).M();
628 _masserr = 0.5*_mass*TMath::Sqrt(phHard->EnergyErr()*phHard->EnergyErr()/phHard->E()/phHard->E() + phSoft->EnergyErr()*phSoft->EnergyErr()/phSoft->E()/phSoft->E());
629 _masserrsmeared = 0.5*_mass*TMath::Sqrt(phHard->EnergyErrSmeared()*phHard->EnergyErrSmeared()/phHard->E()/phHard->E() + phSoft->EnergyErrSmeared()*phSoft->EnergyErrSmeared()/phSoft->E()/phSoft->E());
630 _ptgg = (phHard->Mom()+phSoft->Mom()).Pt();
631 _etagg = (phHard->Mom()+phSoft->Mom()).Eta();
632 _phigg = (phHard->Mom()+phSoft->Mom()).Phi();
633 _costheta = ThreeVector(phHard->Mom()).Unit().Dot(ThreeVector(phSoft->Mom()).Unit());
634 _evtcat = PhotonTools::DiphotonR9EtaPtCat(phHard,phSoft);
635
636 const Double_t dz = sqrt(2.0)*fBeamspotWidth;
637 Double_t deltamvtx = _mass*VertexTools::DeltaMassVtx(phHard->CaloPos().X(),
638 phHard->CaloPos().Y(),
639 phHard->CaloPos().Z(),
640 phSoft->CaloPos().X(),
641 phSoft->CaloPos().Y(),
642 phSoft->CaloPos().Z(),
643 fDiphotonEvent->vtxX,
644 fDiphotonEvent->vtxY,
645 fDiphotonEvent->vtxZ,
646 dz);
647 fDiphotonEvent->deltamvtx = deltamvtx;
648
649 _masserrwrongvtx = TMath::Sqrt(_masserr*_masserr + deltamvtx*deltamvtx);
650 _masserrsmearedwrongvtx = TMath::Sqrt(_masserrsmeared*_masserrsmeared + deltamvtx*deltamvtx);
651
652
653
654 // =================================================================================
655 // this is for synching the Vtx stuff
656 if ( fDoSynching ) {
657 //fill conversion collection for vertex selection, adding single leg conversions if needed
658 //note that momentum of single leg conversions needs to be recomputed from the track
659 //as it is not filled properly
660 DecayParticleOArr vtxconversions;
661 if ( true ) {
662 vtxconversions.SetOwner(kTRUE);
663 for (UInt_t iconv=0; iconv<fConversions->GetEntries(); ++iconv) {
664 DecayParticle *conv = new DecayParticle(*fConversions->At(iconv));
665 vtxconversions.AddOwned(conv);
666 }
667
668 for (UInt_t iconv=0; iconv<fPFConversions->GetEntries(); ++iconv) {
669 const DecayParticle *c = fPFConversions->At(iconv);
670 if (c->NDaughters()!=1) continue;
671
672 DecayParticle *conv = new DecayParticle(*c);
673 const Track *trk = static_cast<const StableParticle*>(conv->Daughter(0))->Trk();
674 conv->SetMom(trk->Px(), trk->Py(), trk->Pz(), trk->P());
675 vtxconversions.AddOwned(conv);
676 }
677 }
678 else {
679 for (UInt_t iconv=0; iconv<fConversions->GetEntries(); ++iconv) {
680 const DecayParticle *c = fConversions->At(iconv);
681 vtxconversions.Add(c);
682 }
683 }
684
685
686 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
687 double tmpVtxProb = 0.;
688
689 std::vector<int> debugInds; // this hold the Vtx indices for the first three guys
690 std::vector<double> debugVals; // this holds hte mva input/output for the best ranked
691 std::vector<int> debugConv; // this holds number of legs for frist and second photon (0 if no conversion) and total number of conversions
692
693
694 realVtx = fVtxTools.findVtxBasicRanking(phHard, phSoft, bsp, fPV,
695 &vtxconversions,kTRUE,tmpVtxProb,
696 &debugInds, &debugVals, &debugConv
697 );
698
699 fDiphotonEvent->vtxInd1 = debugInds[0];
700 fDiphotonEvent->vtxInd2 = debugInds[1];
701 fDiphotonEvent->vtxInd3 = debugInds[2];
702
703 fDiphotonEvent->vtxConv1Z = debugVals[0];
704 fDiphotonEvent->vtxConv1DZ = debugVals[1];
705 fDiphotonEvent->vtxConv1Prob = debugVals[2];
706
707 fDiphotonEvent->vtxConv2Z = debugVals[3];
708 fDiphotonEvent->vtxConv2DZ = debugVals[4];
709 fDiphotonEvent->vtxConv2Prob = debugVals[5];
710
711
712 fDiphotonEvent->vtxBestPtbal = debugVals[6];
713 fDiphotonEvent->vtxBestPtasym = debugVals[7];
714 fDiphotonEvent->vtxBestSumpt2 = debugVals[8];
715 fDiphotonEvent->vtxBestP2Conv = debugVals[9];
716
717 fDiphotonEvent->vtxMva1Z = debugVals[10];
718 fDiphotonEvent->vtxMva2Z = debugVals[11];
719 fDiphotonEvent->vtxMva3Z = debugVals[12];
720
721 fDiphotonEvent->vtxMva1 = debugVals[13];
722 fDiphotonEvent->vtxMva2 = debugVals[14];
723 fDiphotonEvent->vtxMva3 = debugVals[15];
724
725 fDiphotonEvent->vtxNleg1 = debugConv[0];
726 fDiphotonEvent->vtxNleg2 = debugConv[1];
727 fDiphotonEvent->vtxConvIdx1 = debugConv[2];
728 fDiphotonEvent->vtxConvIdx2 = debugConv[3];
729
730 fDiphotonEvent->vtxNconv = debugConv[4];
731
732 if( false ) {
733 printf("---------------------------------------------------------------------\n");
734 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);
735 printf("---------------------------------------------------------------------\n");
736 }
737 }
738
739 // end of Vtx synching stuff ...
740 // =================================================================================
741
742 PFJetOArr pfjets;
743 if (fEnableJets){
744 if (jet1 && jet2) {
745 fDiphotonEvent->zeppenfeld = TMath::Abs(_etagg - 0.5*(jet1->Eta()+jet2->Eta()));
746 fDiphotonEvent->dphidijetgg = MathUtils::DeltaPhi( (jet1->Mom()+jet2->Mom()).Phi(), _phigg );
747 }
748
749 //PFJetOArr pfjets;
750 for (UInt_t ijet=0; ijet<fPFJets->GetEntries(); ++ijet) {
751 const PFJet *pfjet = dynamic_cast<const PFJet*>(fPFJets->At(ijet));
752 if (pfjet && MathUtils::DeltaR(*pfjet,*phHard)>0.3 && MathUtils::DeltaR(*pfjet,*phSoft)>0.3) pfjets.Add(pfjet);
753 }
754 }
755
756 PFCandidateOArr pfcands;
757 for (UInt_t icand=0; icand<fPFCands->GetEntries(); ++icand) {
758 const PFCandidate *pfcand = fPFCands->At(icand);
759 if (MathUtils::DeltaR(*pfcand,*phHard)>0.1 && MathUtils::DeltaR(*pfcand,*phSoft)>0.1) pfcands.Add(pfcand);
760 }
761
762 const Vertex *firstvtx = fPV->At(0);
763 const Vertex *selvtx = fPV->At(0);
764
765 if (!fLoopOnGoodElectrons && phHard->HasPV()) {
766 selvtx = phHard->PV();
767 }
768
769 if (0) //disable for now for performance reasons
770 {
771 Met mmet = fMVAMet.GetMet( false,
772 0.,0.,0.,
773 0.,0.,0.,
774 fPFMet->At(0),
775 &pfcands,selvtx,fPV, fPileUpDen->At(0)->Rho(),
776 &pfjets,
777 int(fPV->GetEntries()),
778 kFALSE);
779
780 TMatrixD *metcov = fMVAMet.GetMetCovariance();
781
782 ThreeVector fullmet(mmet.Px() - phHard->Px() - phSoft->Px(),
783 mmet.Py() - phHard->Py() - phSoft->Py(),
784 0.);
785
786 ROOT::Math::SMatrix<double,2,2,ROOT::Math::MatRepSym<double,2> > mcov;
787 mcov(0,0) = (*metcov)(0,0);
788 mcov(0,1) = (*metcov)(0,1);
789 mcov(1,0) = (*metcov)(1,0);
790 mcov(1,1) = (*metcov)(1,1);
791 ROOT::Math::SVector<double,2> vmet;
792 vmet(0) = fullmet.X();
793 vmet(1) = fullmet.Y();
794 mcov.Invert();
795 Double_t metsig = sqrt(ROOT::Math::Similarity(mcov,vmet));
796
797 fDiphotonEvent->mvametsel = fullmet.Rho();
798 fDiphotonEvent->mvametselphi = fullmet.Phi();
799 fDiphotonEvent->mvametselx = fullmet.X();
800 fDiphotonEvent->mvametsely = fullmet.Y();
801 fDiphotonEvent->mvametselsig = metsig;
802 }
803
804 if (0) //disable for now for performance reasons
805 {
806 Met mmet = fMVAMet.GetMet( false,
807 0.,0.,0.,
808 0.,0.,0.,
809 fPFMet->At(0),
810 &pfcands,firstvtx,fPV, fPileUpDen->At(0)->Rho(),
811 &pfjets,
812 int(fPV->GetEntries()),
813 kFALSE);
814
815 TMatrixD *metcov = fMVAMet.GetMetCovariance();
816
817 ThreeVector fullmet(mmet.Px() - phHard->Px() - phSoft->Px(),
818 mmet.Py() - phHard->Py() - phSoft->Py(),
819 0.);
820
821 ROOT::Math::SMatrix<double,2,2,ROOT::Math::MatRepSym<double,2> > mcov;
822 mcov(0,0) = (*metcov)(0,0);
823 mcov(0,1) = (*metcov)(0,1);
824 mcov(1,0) = (*metcov)(1,0);
825 mcov(1,1) = (*metcov)(1,1);
826 ROOT::Math::SVector<double,2> vmet;
827 vmet(0) = fullmet.X();
828 vmet(1) = fullmet.Y();
829 mcov.Invert();
830 Double_t metsig = sqrt(ROOT::Math::Similarity(mcov,vmet));
831
832 fDiphotonEvent->mvametfirst = fullmet.Rho();
833 fDiphotonEvent->mvametfirstphi = fullmet.Phi();
834 fDiphotonEvent->mvametfirstx = fullmet.X();
835 fDiphotonEvent->mvametfirsty = fullmet.Y();
836 fDiphotonEvent->mvametfirstsig = metsig;
837 }
838
839 }
840
841 fDiphotonEvent->corrpfmet = -99.;
842 fDiphotonEvent->corrpfmetphi = -99.;
843 fDiphotonEvent->corrpfmetx = -99.;
844 fDiphotonEvent->corrpfmety = -99.;
845
846 Met *corrMet =NULL;
847
848 if (fApplyPFMetCorrections){
849 corrMet = new Met(fPFMet->At(0)->Px(),fPFMet->At(0)->Py());
850
851 if (!fIsData){
852 PFMetCorrectionTools::correctMet(corrMet,phHard,phSoft,1,0,funcorrPFJets,fGenJets,fPFJets,_evt);
853 PFMetCorrectionTools::shiftMet(corrMet,fIsData,_spfMet);
854 }
855 else {
856 PFMetCorrectionTools::shiftMet(corrMet,fIsData,_spfMet);
857 PFMetCorrectionTools::correctMet(corrMet,phHard,phSoft,0,1,funcorrPFJets,fGenJets,fPFJets,_evt);
858 }
859
860 fDiphotonEvent->corrpfmet = corrMet->Pt();
861 fDiphotonEvent->corrpfmetphi = corrMet->Phi();
862 fDiphotonEvent->corrpfmetx = corrMet->Px();
863 fDiphotonEvent->corrpfmety = corrMet->Py();
864
865 delete corrMet;
866 }
867
868
869 Float_t _massele = -99.;
870 Float_t _ptee = -99.;
871 Float_t _costhetaele = -99.;
872 if (ele1 && ele2) {
873 _massele = (ele1->Mom()+ele2->Mom()).M();
874 _ptee = (ele1->Mom()+ele2->Mom()).Pt();
875 _costhetaele = ThreeVector(ele1->Mom()).Unit().Dot(ThreeVector(ele2->Mom()).Unit());
876 }
877
878 Float_t _gencostheta = -99.;
879 if (phgen1 && phgen2) {
880 _gencostheta = ThreeVector(phgen1->Mom()).Unit().Dot(ThreeVector(phgen2->Mom()).Unit());
881 }
882
883 fDiphotonEvent->gencostheta = _gencostheta;
884 fDiphotonEvent->dphiMetgg = _dphiMetgg;
885 fDiphotonEvent->cosdphiMetgg = _cosdphiMetgg;
886 fDiphotonEvent->dphiPhPh = _dphiPhPh;
887 fDiphotonEvent->mass = _mass;
888 fDiphotonEvent->masserr = _masserr;
889 fDiphotonEvent->masserrsmeared = _masserrsmeared;
890 fDiphotonEvent->masserrwrongvtx = _masserrwrongvtx;
891 fDiphotonEvent->masserrsmearedwrongvtx = _masserrsmearedwrongvtx;
892 fDiphotonEvent->ptgg = _ptgg;
893 fDiphotonEvent->etagg = _etagg;
894 fDiphotonEvent->phigg = _phigg;
895 fDiphotonEvent->costheta = _costheta;;
896 fDiphotonEvent->massele = _massele;
897 fDiphotonEvent->ptee = _ptee;
898 fDiphotonEvent->costhetaele = _costhetaele;
899 fDiphotonEvent->evtcat = _evtcat;
900
901 fDiphotonEvent->photons[0].SetVars(phHard, conv1, ele1, pfsc1, phgen1,
902 fPhfixph, fPhfixele, fTracks, fPV,
903 fPFCands, rho, fFillClusterArrays,
904 fPhotons, fPFSuperClusters,
905 fElectrons, fConversions, bsp,
906 fApplyElectronVeto, realVtx);
907 fDiphotonEvent->photons[1].SetVars(phSoft, conv2, ele2, pfsc2, phgen2,
908 fPhfixph, fPhfixele, fTracks, fPV,
909 fPFCands, rho, fFillClusterArrays,
910 fPhotons, fPFSuperClusters,
911 fElectrons, fConversions,
912 bsp, fApplyElectronVeto, realVtx);
913
914 Float_t ph1ecor = fDiphotonEvent->photons[0].Ecor();
915 Float_t ph1ecorerr = fDiphotonEvent->photons[0].Ecorerr();
916 Float_t ph2ecor = fDiphotonEvent->photons[1].Ecor();
917 Float_t ph2ecorerr = fDiphotonEvent->photons[1].Ecorerr();
918
919 Float_t ph1ecorele = fDiphotonEvent->photons[0].Ecorele();
920 Float_t ph1ecoreleerr = fDiphotonEvent->photons[0].Ecoreleerr();
921 Float_t ph2ecorele = fDiphotonEvent->photons[1].Ecorele();
922 Float_t ph2ecoreleerr = fDiphotonEvent->photons[1].Ecoreleerr();
923
924 fDiphotonEvent->masscor = TMath::Sqrt(2.0*ph1ecor*ph2ecor*(1.0-fDiphotonEvent->costheta));
925 fDiphotonEvent->masscorerr = 0.5*fDiphotonEvent->masscor*
926 TMath::Sqrt(ph1ecorerr*ph1ecorerr/ph1ecor/ph1ecor + ph2ecorerr*ph2ecorerr/ph2ecor/ph2ecor);
927
928 fDiphotonEvent->masscorele = TMath::Sqrt(2.0*ph1ecorele*ph2ecorele*
929 (1.0-fDiphotonEvent->costheta));
930 fDiphotonEvent->masscoreleerr = 0.5*fDiphotonEvent->masscorele*
931 TMath::Sqrt(ph1ecoreleerr*ph1ecoreleerr/ph1ecorele/ph1ecorele +
932 ph2ecoreleerr*ph2ecoreleerr/ph2ecorele/ph2ecorele);
933
934 //printf("r9 = %5f, photon sigieie = %5f, seed sigieie = %5f\n",phHard->R9(),
935 // phHard->CoviEtaiEta(),sqrt(phHard->SCluster()->Seed()->CoviEtaiEta()));
936
937 // MuonStuff
938 fDiphotonEvent-> muonPt = -99.;
939 fDiphotonEvent-> muonEta = -99.;
940 fDiphotonEvent-> muonPhi = -99.;
941 fDiphotonEvent-> muDR1 = -99.;
942 fDiphotonEvent-> muDR2 = -99.;
943 fDiphotonEvent-> muIso1 = -99.;
944 fDiphotonEvent-> muIso2 = -99.;
945 fDiphotonEvent-> muIso3 = -99.;
946 fDiphotonEvent-> muIso4 = -99.;
947 fDiphotonEvent-> muD0 = -99.;
948 fDiphotonEvent-> muDZ = -99.;
949 fDiphotonEvent-> muChi2 = -99.;
950 fDiphotonEvent-> muNhits = -99;
951 fDiphotonEvent-> muNpixhits = -99;
952 fDiphotonEvent-> muNegs = -99;
953 fDiphotonEvent-> muNMatch = -99;
954
955 // Dimuon Stuff
956 fDiphotonEvent-> mu1Pt = -99.;
957 fDiphotonEvent-> mu1Eta = -99.;
958 fDiphotonEvent-> mu1Phi = -99.;
959 fDiphotonEvent-> mu2Pt = -99.;
960 fDiphotonEvent-> mu2Eta = -99.;
961 fDiphotonEvent-> mu2Phi = -99.;
962
963 // Electron Stuff
964 fDiphotonEvent-> elePt = -99.;
965 fDiphotonEvent-> eleEta = -99.;
966 fDiphotonEvent-> elePhi = -99.;
967 fDiphotonEvent-> eleSCEta = -99.;
968 fDiphotonEvent-> eleIso1 = -99.;
969 fDiphotonEvent-> eleIso2 = -99.;
970 fDiphotonEvent-> eleIso3 = -99.;
971 fDiphotonEvent-> eleIso4 = -99.;
972 fDiphotonEvent-> eleDist = -99.;
973 fDiphotonEvent-> eleDcot = -99.;
974 fDiphotonEvent-> eleCoviee = -99.;
975 fDiphotonEvent-> eleDphiin = -99.;
976 fDiphotonEvent-> eleDetain = -99.;
977 fDiphotonEvent-> eleDR1 = -99.;
978 fDiphotonEvent-> eleDR2 = -99.;
979 fDiphotonEvent-> eleMass1 = -99.;
980 fDiphotonEvent-> eleMass2 = -99.;
981 fDiphotonEvent-> eleNinnerHits = -99;
982 fDiphotonEvent-> eleIdMva = -99.;
983
984 // Dielectron Stuff
985 fDiphotonEvent-> ele1Pt = -99.;
986 fDiphotonEvent-> ele1Eta = -99.;
987 fDiphotonEvent-> ele1Phi = -99.;
988 fDiphotonEvent-> ele2Pt = -99.;
989 fDiphotonEvent-> ele2Eta = -99.;
990 fDiphotonEvent-> ele2Phi = -99.;
991
992
993 if( fApplyLeptonTag ) {
994 ApplyLeptonTag(phHard, phSoft, selvtx);
995 }
996
997 if( fApplyVHLepTag ) {
998 ApplyVHLepTag(phHard, phSoft, selvtx);
999 }
1000
1001 fDiphotonEvent->costhetastar = -99.;
1002
1003 if( fApplyVHHadTag && jet1 && jet2) {
1004 ApplyVHHadTag(phHard, phSoft, selvtx, jet1, jet2);
1005 }
1006
1007 //vbf tag
1008 fDiphotonEvent->vbfTag = -1;
1009 fDiphotonEvent->vbfbdt = -99;
1010 if(fApplyVBFTag){
1011 fDiphotonEvent->vbfTag = 0;
1012 jet1pt_vbf = fDiphotonEvent->jet1pt;
1013 jet2pt_vbf = fDiphotonEvent->jet2pt;
1014 deltajeteta_vbf = TMath::Abs(fDiphotonEvent->jet1eta - fDiphotonEvent->jet2eta);
1015 dijetmass_vbf = fDiphotonEvent->dijetmass;
1016 zeppenfeld_vbf = TMath::Abs(fDiphotonEvent->zeppenfeld);
1017 dphidijetgg_vbf = TMath::Abs(fDiphotonEvent->dphidijetgg);
1018 diphoptOverdiphomass_vbf = fDiphotonEvent->ptgg/fDiphotonEvent->mass;
1019 pho1ptOverdiphomass_vbf = phHard->Pt()/fDiphotonEvent->mass;
1020 pho2ptOverdiphomass_vbf = phSoft->Pt()/fDiphotonEvent->mass;
1021
1022 if(jet1 && jet2 && (pho1ptOverdiphomass_vbf > 40./120.) && (pho2ptOverdiphomass_vbf > 30./120.) && (jet1pt_vbf > 30) && (jet2pt_vbf > 20) && (dijetmass_vbf > 250)){
1023 fDiphotonEvent->vbfbdt = fMVAVBF.GetMVAbdtValue(jet1pt_vbf,jet2pt_vbf,deltajeteta_vbf,dijetmass_vbf,zeppenfeld_vbf,dphidijetgg_vbf,diphoptOverdiphomass_vbf,pho1ptOverdiphomass_vbf,pho2ptOverdiphomass_vbf);
1024
1025 if((fDiphotonEvent->vbfbdt>=0.985) && (fDiphotonEvent->vbfbdt<=1)){
1026 fDiphotonEvent->vbfTag = 2;
1027 }
1028 if((fDiphotonEvent->vbfbdt>=0.93) && (fDiphotonEvent->vbfbdt<0.985)){
1029 fDiphotonEvent->vbfTag = 1;
1030 }
1031 }
1032 } // End of VBF tag stuff
1033
1034 // ttH tag stuff
1035 fDiphotonEvent->tthTag = -1;
1036 if (fApplyTTHTag && phHard && phSoft && selvtx) {
1037 ApplyTTHTag(phHard, phSoft, selvtx);
1038 }
1039
1040 fDiphotonEvent->numJets = -99;
1041 fDiphotonEvent->numBJets = -99;
1042 fDiphotonEvent->bjetcsv = -99;
1043
1044 if (fDoSynching) {
1045 double minJetPt = 20.;
1046 double maxJetAbsEta = 4.7;
1047 fDiphotonEvent->numJets = NumberOfJets (phHard, phSoft, selvtx, minJetPt,
1048 maxJetAbsEta);
1049 fDiphotonEvent->numBJets = NumberOfBJets(phHard, phSoft, selvtx, minJetPt,
1050 maxJetAbsEta);
1051 /// The lines below cause a segmentation violation in the call to
1052 /// GetSelectedPFJets
1053 // DeltaRVetoVector vetos(2);
1054 // vetos.push_back(DeltaRVeto(static_cast<const Particle*>(phHard), 0.5));
1055 // vetos.push_back(DeltaRVeto(static_cast<const Particle*>(phSoft), 0.5));
1056 // // for (unsigned
1057 // PFJetVector *jets = GetSelectedPFJets(vetos, *selvtx, minJetPt,
1058 // maxJetAbsEta);
1059 // fDiphotonEvent->numJets = jets->size();
1060 // delete jets;
1061
1062 }
1063
1064
1065 //printf("vbfbdt:%f\n",fDiphotonEvent->vbfbdt);
1066 if (fWriteDiphotonTree)
1067 hCiCTuple->Fill();
1068
1069 if (fFillVertexTree && phHard && phSoft) {
1070 for (UInt_t ivtx = 0; ivtx<fPV->GetEntries(); ++ivtx) {
1071 const Vertex *v = fPV->At(ivtx);
1072 fDiphotonVtx->SetVars(v,phHard,phSoft,fPFCands,ivtx,fPV->GetEntries(),_decayZ);
1073 hVtxTree->Fill();
1074 }
1075 }
1076 }
1077
1078 if (!fWriteSingleTree)
1079 return;
1080
1081 for (UInt_t iph = 0; iph<egcol->GetEntries(); ++iph) {
1082 const Particle *p = 0;
1083 const Photon *ph = 0;
1084 const Electron *ele = 0;
1085 const SuperCluster *sc = 0;
1086 if (fLoopOnGoodElectrons) {
1087 ele = fGoodElectrons->At(iph);
1088 p = ele;
1089 sc = ele->SCluster();
1090 ph = PhotonTools::MatchedPhoton(ele,fPhotons);
1091 }
1092 else {
1093 ph = fPhotons->At(iph);
1094 p = ph;
1095 sc = ph->SCluster();
1096 ele = PhotonTools::MatchedElectron(ph,fGoodElectrons);
1097 }
1098
1099 const DecayParticle *conv = PhotonTools::MatchedConversion(sc,fConversions,bsp,
1100 nhitsbeforevtxmax);
1101 const SuperCluster *pfsc = PhotonTools::MatchedPFSC(sc,fPFPhotons, fElectrons);
1102
1103 if (!fLoopOnGoodElectrons && ph->HasPV()) {
1104 fDiphotonEvent->vtxZ = ph->PV()->Z();
1105 }
1106
1107 const MCParticle *phgen = 0;
1108 if( !fIsData ) {
1109 phgen = PhotonTools::MatchMC(p,fMCParticles,!fApplyElectronVeto);
1110 }
1111
1112 if (fExcludeSinglePrompt && phgen) return;
1113
1114 fDiphotonEvent->mt = -99.;
1115 fDiphotonEvent->cosphimet = -99.;
1116 fDiphotonEvent->mtele = -99.;
1117 fDiphotonEvent->cosphimetele = -99.;
1118
1119 if (ph) {
1120 fDiphotonEvent->cosphimet = TMath::Cos(ph->Phi()-fPFMet->At(0)->Phi());
1121 fDiphotonEvent->mt = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ph->Pt()*
1122 (1.0-fDiphotonEvent->cosphimet));
1123 }
1124
1125 if (ele) {
1126 fDiphotonEvent->cosphimetele = TMath::Cos(ele->Phi()-fPFMet->At(0)->Phi());
1127 fDiphotonEvent->mtele = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ele->Pt()*
1128 (1.0-fDiphotonEvent->cosphimetele));
1129 }
1130
1131 fSinglePhoton->SetVars(ph, conv, ele, pfsc, phgen, fPhfixph, fPhfixele,
1132 fTracks, fPV, fPFCands, rho, fFillClusterArrays,
1133 fPhotons, fPFSuperClusters, fElectrons, fConversions,
1134 bsp, fApplyElectronVeto);
1135 hCiCTupleSingle->Fill();
1136 }
1137
1138
1139
1140 return;
1141 }
1142
1143 //--------------------------------------------------------------------------------------------------
1144 void PhotonTreeWriter::SlaveBegin()
1145 {
1146 // Run startup code on the computer (slave) doing the actual analysis. Here,
1147 // we just request the photon collection branch.
1148
1149 if( fApplyLeptonTag || fApplyVHLepTag ) {
1150 ReqEventObject(fLeptonTagElectronsName, fLeptonTagElectrons, false);
1151 ReqEventObject(fLeptonTagMuonsName, fLeptonTagMuons, false);
1152 }
1153
1154 if( fApplyVHLepTag ) {
1155 ReqEventObject(fLeptonTagSoftElectronsName, fLeptonTagSoftElectrons, false);
1156 ReqEventObject(fLeptonTagSoftMuonsName, fLeptonTagSoftMuons, false);
1157 }
1158
1159 // ReqEventObject(fPFNoPileUpName, fPFNoPileUpCands, false);
1160 // ReqEventObject(fPFPileUpName, fPFPileUpCands, false);
1161
1162 ReqEventObject(fPhotonBranchName,fPhotons, fPhotonsFromBranch);
1163 if (fEnablePFPhotons) ReqEventObject(fPFPhotonName,fPFPhotons, true);
1164 ReqEventObject(fTrackBranchName, fTracks, true);
1165 ReqEventObject(fElectronName, fElectrons, true);
1166 ReqEventObject(fGoodElectronName,fGoodElectrons,fGoodElectronsFromBranch);
1167 ReqEventObject(fPileUpDenName, fPileUpDen, true);
1168 ReqEventObject(fPVName, fPV, fPVFromBranch);
1169 ReqEventObject(fConversionName, fConversions, true);
1170 if ( fDoSynching ) ReqEventObject(fPFConversionName, fPFConversions, true);
1171 ReqEventObject(fBeamspotName, fBeamspot, true);
1172 ReqEventObject(fPFCandName, fPFCands, true);
1173 ReqEventObject(fPFSuperClusterName,fPFSuperClusters,true);
1174 ReqEventObject(fPFMetName, fPFMet, true);
1175 if (fEnableJets){
1176 ReqEventObject(fPFJetName, fPFJets, fPFJetsFromBranch);
1177 ReqBranch(funcorrPFJetName, funcorrPFJets);
1178 // if (!fIsData) ReqEventObject(fGenJetName, fGenJets, true);
1179 }
1180 if (!fIsData) {
1181 ReqBranch(fPileUpName, fPileUp);
1182 ReqBranch(fMCParticleName, fMCParticles);
1183 ReqBranch(fMCEventInfoName, fMCEventInfo);
1184 if (fEnableGenJets) ReqEventObject(fGenJetName, fGenJets, true);
1185 }
1186 if (fIsData) {
1187 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
1188 TString("/src/MitPhysics/data/PhotonFixGRPV22.dat");
1189 }
1190 else {
1191 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
1192 TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat");
1193 }
1194
1195 fPhfixph.initialise("4_2",std::string(fPhFixDataFile));
1196 fPhfixele.initialise("4_2e",std::string(fPhFixDataFile));
1197
1198 // fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
1199 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
1200 // TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
1201 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_42.root"))),
1202 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_42.root"))),
1203 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1_42.root"))),
1204 // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2_42.root")))
1205 // );
1206
1207 fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
1208 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
1209 TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
1210 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_52.root"))),
1211 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_52.root"))),
1212 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1cov_52.root"))),
1213 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2cov_52.root")))
1214 );
1215
1216 fJetId.Initialize(JetIDMVA::kMedium,
1217 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
1218 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
1219 JetIDMVA::kCut,
1220 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")))
1221
1222 );
1223
1224 fMVAVBF.InitializeMVA();
1225
1226
1227 if( fDoSynching ) {
1228 fVtxTools.InitP(2);
1229 fElectronIDMVA = new ElectronIDMVA();
1230 fElectronIDMVA->Initialize("BDTG method",
1231 fElectronMVAWeights_Subdet0Pt10To20,
1232 fElectronMVAWeights_Subdet1Pt10To20,
1233 fElectronMVAWeights_Subdet2Pt10To20,
1234 fElectronMVAWeights_Subdet0Pt20ToInf,
1235 fElectronMVAWeights_Subdet1Pt20ToInf,
1236 fElectronMVAWeights_Subdet2Pt20ToInf,
1237 ElectronIDMVA::kIDEGamma2012NonTrigV1,
1238 fTheRhoType);
1239 }
1240
1241
1242 fDiphotonEvent = new PhotonTreeWriterDiphotonEvent;
1243 fSinglePhoton = new PhotonTreeWriterPhoton<16>;
1244
1245 fTmpFile = TFile::Open(TString::Format("%s_tmp.root",GetName()),"RECREATE");
1246
1247 if (fWriteDiphotonTree) {
1248 hCiCTuple = new TTree(fTupleName.Data(),fTupleName.Data());
1249 hCiCTuple->SetAutoSave(300e9);
1250 hCiCTuple->SetDirectory(fTmpFile);
1251 }
1252 TString singlename = fTupleName + TString("Single");
1253 if (fWriteSingleTree) {
1254 hCiCTupleSingle = new TTree(singlename,singlename);
1255 hCiCTupleSingle->SetAutoSave(300e9);
1256 hCiCTupleSingle->SetDirectory(fTmpFile);
1257 }
1258
1259 //make flattish tree from classes so we don't have to rely on dictionaries for reading later
1260 TClass *eclass = TClass::GetClass("mithep::PhotonTreeWriterDiphotonEvent");
1261 TClass *pclass = TClass::GetClass("mithep::PhotonTreeWriterPhoton<16>");
1262 TList *elist = eclass->GetListOfDataMembers();
1263 TList *plist = pclass->GetListOfDataMembers();
1264
1265 for (int i=0; i<elist->GetEntries(); ++i) {
1266 const TDataMember *tdm = static_cast<const TDataMember*>(elist->At(i));//ming
1267 if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
1268 TString typestring;
1269 if (TString(tdm->GetTypeName()).BeginsWith("Char_t")) typestring = "B";
1270 else if (TString(tdm->GetTypeName()).BeginsWith("UChar_t")) typestring = "b";
1271 else if (TString(tdm->GetTypeName()).BeginsWith("Short_t")) typestring = "S";
1272 else if (TString(tdm->GetTypeName()).BeginsWith("UShort_t")) typestring = "s";
1273 else if (TString(tdm->GetTypeName()).BeginsWith("Int_t")) typestring = "I";
1274 else if (TString(tdm->GetTypeName()).BeginsWith("UInt_t")) typestring = "i";
1275 else if (TString(tdm->GetTypeName()).BeginsWith("Float_t")) typestring = "F";
1276 else if (TString(tdm->GetTypeName()).BeginsWith("Double_t")) typestring = "D";
1277 else if (TString(tdm->GetTypeName()).BeginsWith("Long64_t")) typestring = "L";
1278 else if (TString(tdm->GetTypeName()).BeginsWith("ULong64_t")) typestring = "l";
1279 else if (TString(tdm->GetTypeName()).BeginsWith("Bool_t")) typestring = "O";
1280 else continue;
1281 //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
1282 Char_t *addr = (Char_t*)fDiphotonEvent;//ming:?
1283 assert(sizeof(Char_t)==1);
1284 if (fWriteDiphotonTree) hCiCTuple->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
1285 if (fWriteSingleTree) hCiCTupleSingle->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
1286 }
1287
1288 for (int iph=0; iph<2; ++iph) {
1289 for (int i=0; i<plist->GetEntries(); ++i) {
1290 const TDataMember *tdm = static_cast<const TDataMember*>(plist->At(i));
1291 if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
1292 TString typestring;
1293 if (TString(tdm->GetTypeName()).BeginsWith("Char_t")) typestring = "B";
1294 else if (TString(tdm->GetTypeName()).BeginsWith("UChar_t")) typestring = "b";
1295 else if (TString(tdm->GetTypeName()).BeginsWith("Short_t")) typestring = "S";
1296 else if (TString(tdm->GetTypeName()).BeginsWith("UShort_t")) typestring = "s";
1297 else if (TString(tdm->GetTypeName()).BeginsWith("Int_t")) typestring = "I";
1298 else if (TString(tdm->GetTypeName()).BeginsWith("UInt_t")) typestring = "i";
1299 else if (TString(tdm->GetTypeName()).BeginsWith("Float_t")) typestring = "F";
1300 else if (TString(tdm->GetTypeName()).BeginsWith("Double_t")) typestring = "D";
1301 else if (TString(tdm->GetTypeName()).BeginsWith("Long64_t")) typestring = "L";
1302 else if (TString(tdm->GetTypeName()).BeginsWith("ULong64_t")) typestring = "l";
1303 else if (TString(tdm->GetTypeName()).BeginsWith("Bool_t")) typestring = "O";
1304 else continue;
1305 //printf("%s\n",tdm->GetTypeName());
1306 TString varname = TString::Format("ph%d.%s",iph+1,tdm->GetName());
1307 if (tdm->GetArrayDim()==1) {
1308 varname = TString::Format("%s[%i]",varname.Data(),tdm->GetMaxIndex(0));
1309 }
1310
1311 //printf("typename = %s, arraydim = %i, arraysize = %i,varname = %s\n", tdm->GetTypeName(), tdm->GetArrayDim(), tdm->GetMaxIndex(0), varname.Data());
1312
1313 Char_t *addr = (Char_t*)&fDiphotonEvent->photons[iph];
1314 assert(sizeof(Char_t)==1);
1315 if (fWriteDiphotonTree) hCiCTuple->Branch(varname,addr+tdm->GetOffset(),TString::Format("%s/%s",varname.Data(),typestring.Data()));
1316
1317 if (iph==0) {
1318 TString singlename = TString::Format("ph.%s",tdm->GetName());
1319 if (tdm->GetArrayDim()==1) {
1320 singlename = TString::Format("%s[%i]",singlename.Data(),tdm->GetMaxIndex(0));
1321 }
1322 Char_t *addrsingle = (Char_t*)fSinglePhoton;
1323 if (fWriteSingleTree) hCiCTupleSingle->Branch(singlename,addrsingle+tdm->GetOffset(),TString::Format("%s/%s",singlename.Data(),typestring.Data()));
1324 }
1325 }
1326 }
1327
1328 if (fWriteDiphotonTree)
1329 AddOutput(hCiCTuple);
1330 if (fWriteSingleTree)
1331 AddOutput(hCiCTupleSingle);
1332
1333 if (fFillVertexTree) {
1334 fDiphotonVtx = new PhotonTreeWriterVtx;
1335 hVtxTree = new TTree("hVtxTree","hVtxTree");
1336 hVtxTree->SetAutoSave(300e9);
1337 AddOutput(hVtxTree);
1338
1339 TClass *vclass = TClass::GetClass("mithep::PhotonTreeWriterVtx");
1340 TList *vlist = vclass->GetListOfDataMembers();
1341
1342 for (int i=0; i<vlist->GetEntries(); ++i) {
1343 const TDataMember *tdm = static_cast<const TDataMember*>(vlist->At(i));
1344 if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
1345 TString typestring;
1346 if (TString(tdm->GetTypeName()).BeginsWith("Char_t")) typestring = "B";
1347 else if (TString(tdm->GetTypeName()).BeginsWith("UChar_t")) typestring = "b";
1348 else if (TString(tdm->GetTypeName()).BeginsWith("Short_t")) typestring = "S";
1349 else if (TString(tdm->GetTypeName()).BeginsWith("UShort_t")) typestring = "s";
1350 else if (TString(tdm->GetTypeName()).BeginsWith("Int_t")) typestring = "I";
1351 else if (TString(tdm->GetTypeName()).BeginsWith("UInt_t")) typestring = "i";
1352 else if (TString(tdm->GetTypeName()).BeginsWith("Float_t")) typestring = "F";
1353 else if (TString(tdm->GetTypeName()).BeginsWith("Double_t")) typestring = "D";
1354 else if (TString(tdm->GetTypeName()).BeginsWith("Long64_t")) typestring = "L";
1355 else if (TString(tdm->GetTypeName()).BeginsWith("ULong64_t")) typestring = "l";
1356 else if (TString(tdm->GetTypeName()).BeginsWith("Bool_t")) typestring = "O";
1357 else continue;
1358 //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
1359 Char_t *addr = (Char_t*)fDiphotonVtx;
1360 assert(sizeof(Char_t)==1);
1361 hVtxTree->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
1362 }
1363 }
1364
1365 }
1366
1367 // ----------------------------------------------------------------------------------------
1368 // some helpfer functions....
1369 void PhotonTreeWriter::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass)
1370 {
1371 pt = -999.;
1372 decayZ = -999.;
1373 mass = -999.;
1374
1375 // loop over all GEN particles and look for status 1 photons
1376 for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
1377 const MCParticle* p = fMCParticles->At(i);
1378 if (p->Is(MCParticle::kH) || (!fApplyElectronVeto &&
1379 (p->AbsPdgId()==23 || p->AbsPdgId()==24))) {
1380 pt = p->Pt();
1381 decayZ = p->DecayVertex().Z();
1382 mass = p->Mass();
1383 break;
1384 }
1385 }
1386
1387 return;
1388 }
1389
1390
1391 Float_t PhotonTreeWriter::GetEventCat(PhotonTools::CiCBaseLineCats cat1,
1392 PhotonTools::CiCBaseLineCats cat2) {
1393
1394 bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
1395 bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
1396
1397 bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
1398 bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
1399
1400 if( ph1IsEB && ph2IsEB )
1401 return ( ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
1402
1403 return ( ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
1404 }
1405
1406 template <int NClus>
1407 void
1408 PhotonTreeWriterPhoton<NClus>::SetVars(const Photon *p,
1409 const DecayParticle *c,
1410 const Electron *ele,
1411 const SuperCluster *pfsc,
1412 const MCParticle *m,
1413 PhotonFix &phfixph,
1414 PhotonFix &phfixele,
1415 const TrackCol* trackCol,
1416 const VertexCol* vtxCol,
1417 const PFCandidateCol* fPFCands,
1418 Double_t rho,
1419 Bool_t fillclusterarrays,
1420 const PhotonCol* ps,
1421 const SuperClusterCol* scs,
1422 const ElectronCol* els,
1423 const DecayParticleCol *convs,
1424 const BaseVertex *bs,
1425 Bool_t applyElectronVeto,
1426 const Vertex* realVtx)
1427 {
1428
1429 const SuperCluster *s = 0;
1430 if (p)
1431 s = p->SCluster();
1432 else
1433 s = ele->SCluster();
1434 const BasicCluster *b = s->Seed();
1435 const BasicCluster *b2 = 0;
1436 Double_t ebcmax = -99.;
1437 for (UInt_t i=0; i<s->ClusterSize(); ++i) {
1438 const BasicCluster *bc = s->Cluster(i);
1439 if (bc->Energy() > ebcmax && bc !=b) {
1440 b2 = bc;
1441 ebcmax = bc->Energy();
1442 }
1443 }
1444
1445 const BasicCluster *bclast = 0;
1446 Double_t ebcmin = 1e6;
1447 for (UInt_t i=0; i<s->ClusterSize(); ++i) {
1448 const BasicCluster *bc = s->Cluster(i);
1449 if (bc->Energy() < ebcmin && bc !=b) {
1450 bclast = bc;
1451 ebcmin = bc->Energy();
1452 }
1453 }
1454
1455 const BasicCluster *bclast2 = 0;
1456 ebcmin = 1e6;
1457 for (UInt_t i=0; i<s->ClusterSize(); ++i) {
1458 const BasicCluster *bc = s->Cluster(i);
1459 if (bc->Energy() < ebcmin && bc !=b && bc!=bclast) {
1460 bclast2 = bc;
1461 ebcmin = bc->Energy();
1462 }
1463 }
1464
1465 if (p) {
1466 hasphoton = kTRUE;
1467 index = PhotonTreeWriter::IndexOfElementInCollection(p, ps);
1468
1469 e = p->E();
1470 pt = p->Pt();
1471 eta = p->Eta();
1472 phi = p->Phi();
1473 r9 = p->R9();
1474 e3x3 = p->E33();
1475 e5x5 = p->E55();
1476 hovere = p->HadOverEm();
1477 hoveretower = p->HadOverEmTow();
1478 sigietaieta = p->CoviEtaiEta();
1479 phcat = PhotonTools::CiCBaseLineCat(p);
1480 eerr = p->EnergyErr();
1481 eerrsmeared = p->EnergyErrSmeared();
1482 esmearing = p->EnergySmearing();
1483 eerrsmearing = p->EnergyErrSmearing();
1484 idmva = p->IdMva();
1485 hcalisodr03 = p->HcalTowerSumEtDr03();
1486 ecalisodr03 = p->EcalRecHitIsoDr03();
1487 trkisohollowdr03 = p->HollowConeTrkIsoDr03();
1488 hcalisodr04 = p->HcalTowerSumEtDr04();
1489 ecalisodr04 = p->EcalRecHitIsoDr04();
1490 trkisohollowdr04 = p->HollowConeTrkIsoDr04();
1491
1492 passeleveto = PhotonTools::PassElectronVetoConvRecovery(p, els, convs, bs);
1493
1494 const Vertex *vtx = vtxCol->At(0);
1495 if (p->HasPV()) vtx = p->PV();
1496 if ( realVtx ) vtx = realVtx;
1497
1498 UInt_t wVtxInd = 0;
1499
1500 trackiso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0,
1501 trackCol, NULL, NULL,
1502 (!applyElectronVeto ? els : NULL) );
1503 //Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
1504
1505 // track iso worst vtx
1506 trackiso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0,
1507 trackCol, &wVtxInd,vtxCol,
1508 (!applyElectronVeto ? els : NULL) );
1509 combiso1 = ecalisodr03+hcalisodr04+trackiso1 - 0.17*rho;
1510 combiso2 = ecalisodr04+hcalisodr04+trackiso2 - 0.52*rho;
1511
1512
1513 // -----------------------------------------------------
1514 // PF-CiC4 Debug Stuff
1515 std::vector<double> debugVals;
1516 PhotonTools::PassCiCPFIsoSelection(p, vtx, fPFCands, vtxCol, rho, 20., &debugVals);
1517 if( debugVals.size() == 13 ) {
1518 pfcic4_tIso1 = debugVals[0];
1519 pfcic4_tIso2 = debugVals[1];
1520 pfcic4_tIso3 = debugVals[2];
1521 pfcic4_covIEtaIEta = debugVals[3];
1522 pfcic4_HoE = debugVals[4];
1523 pfcic4_R9 = debugVals[5];
1524 pfcic4_wVtxInd = debugVals[6];
1525 pfcic4_ecalIso3 = debugVals[7];
1526 pfcic4_ecalIso4 = debugVals[8];
1527 pfcic4_trackIsoSel03 = debugVals[9];
1528 pfcic4_trackIsoWorst04 = debugVals[10];
1529 pfcic4_combIso1 = debugVals[11];
1530 pfcic4_combIso2 = debugVals[12];
1531 }
1532 // -----------------------------------------------------
1533 //id mva
1534 //2011
1535 idmva_tIso1abs=combiso1;
1536 idmva_tIso2abs=combiso2;
1537 idmva_tIso3abs=trackiso1;
1538 idmva_absIsoEcal=ecalisodr03;
1539 idmva_absIsoHcal=hcalisodr04;
1540 //2012
1541 idmva_CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
1542 idmva_s4ratio=p->S4Ratio();
1543 idmva_GammaIso=IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
1544 idmva_ChargedIso_selvtx=IsolationTools::PFChargedIsolation(p,vtx,0.3,0.,fPFCands);
1545 idmva_ChargedIso_0p2_selvtx=IsolationTools::PFChargedIsolation(p,vtx,0.2,0.,fPFCands);
1546 idmva_ChargedIso_worstvtx=IsolationTools::PFChargedIsolation(p,vtx,0.3,0.,fPFCands,&wVtxInd,vtxCol);
1547 idmva_PsEffWidthSigmaRR=p->EffSigmaRR();
1548 }
1549 else {
1550 hasphoton = kFALSE;
1551 index = (UInt_t) -1;
1552 e = -99.;
1553 pt = -99.;
1554 eta = -99.;
1555 phi = -99.;
1556 r9 = b->E3x3()/s->RawEnergy();
1557 e3x3 = b->E3x3();
1558 e5x5 = b->E5x5();
1559 hovere = ele->HadronicOverEm();
1560 sigietaieta = ele->CoviEtaiEta();
1561 phcat = -99;
1562 eerr = -99.;
1563 eerrsmeared = -99.;
1564 esmearing = 0.;
1565 idmva = -99.;
1566 hcalisodr03 = -99;
1567 ecalisodr03 = -99;
1568 trkisohollowdr03 = -99;
1569 hcalisodr04 = -99;
1570 ecalisodr04 = -99;
1571 trkisohollowdr04 = -99;
1572 trackiso1 = -99.;
1573 trackiso2 = -99.;
1574 combiso1 = -99.;
1575 combiso2 = -99.;
1576 }
1577
1578 // TODO: fix the bug with supercluster index
1579 scindex = PhotonTreeWriter::IndexOfNearestSuperClusterInCollection(s, scs);
1580 /// DEBUG
1581 // cout << "JV: s, scs, index: " << s << ", " << scs << ", " << scindex << endl;
1582 sce = s->Energy();
1583 scrawe = s->RawEnergy();
1584 scpse = s->PreshowerEnergy();
1585 scpssigmaxx = s->PsEffWidthSigmaXX();
1586 scpssigmayy = s->PsEffWidthSigmaYY();
1587 sceta = s->Eta();
1588 scphi = s->Phi();
1589 scnclusters = s->ClusterSize();
1590 scnhits = s->NHits();
1591 scetawidth = -99.;
1592 scphiwidth = -99.;
1593 if (p) {
1594 scetawidth = p->EtaWidth();
1595 scphiwidth = p->PhiWidth();
1596 }
1597 else {
1598 scetawidth = s->EtaWidth();
1599 scphiwidth = s->PhiWidth();
1600 }
1601 isbarrel = (s->AbsEta()<1.5);
1602 isr9reco = (isbarrel && r9>0.94) || (!isbarrel && r9>0.95);
1603 isr9cat = (r9>0.94);
1604
1605 eseed = b->Energy();
1606 etaseed = b->Eta();
1607 phiseed = b->Phi();
1608 ietaseed = b->IEta();
1609 iphiseed = b->IPhi();
1610 ixseed = b->IX();
1611 iyseed = b->IY();
1612 etacryseed = b->EtaCry();
1613 phicryseed = b->PhiCry();
1614 xcryseed = b->XCry();
1615 ycryseed = b->YCry();
1616 thetaaxisseed = b->ThetaAxis();
1617 phiaxisseed = b->PhiAxis();
1618 sigietaietaseed = TMath::Sqrt(b->CoviEtaiEta());
1619 sigiphiphiseed = TMath::Sqrt(b->CoviPhiiPhi());
1620 if (isnan(sigiphiphiseed)) sigiphiphiseed = -99.;
1621 covietaiphiseed = b->CoviEtaiPhi();
1622 if (isnan(covietaiphiseed)) covietaiphiseed = -99.;
1623 e3x3seed = b->E3x3();
1624 e5x5seed = b->E5x5();
1625 emaxseed = b->EMax();
1626 e2ndseed = b->E2nd();
1627 etopseed = b->ETop();
1628 ebottomseed = b->EBottom();
1629 eleftseed = b->ELeft();
1630 erightseed = b->ERight();
1631 e1x3seed = b->E1x3();
1632 e3x1seed = b->E3x1();
1633 e1x5seed = b->E1x5();
1634 e2x2seed = b->E2x2();
1635 e4x4seed = b->E4x4();
1636 e2x5maxseed = b->E2x5Max();
1637 e2x5topseed = b->E2x5Top();
1638 e2x5bottomseed = b->E2x5Bottom();
1639 e2x5leftseed = b->E2x5Left();
1640 e2x5rightseed = b->E2x5Right();
1641 xseedseed = b->Pos().X();
1642 yseedseed = b->Pos().Y();
1643 zseedseed = b->Pos().Z();
1644 nhitsseed = b->NHits();
1645
1646 if (b2) {
1647 ebc2 = b2->Energy();
1648 etabc2 = b2->Eta();
1649 phibc2 = b2->Phi();
1650 ietabc2 = b2->IEta();
1651 iphibc2 = b2->IPhi();
1652 ixbc2 = b2->IX();
1653 iybc2 = b2->IY();
1654 etacrybc2 = b2->EtaCry();
1655 phicrybc2 = b2->PhiCry();
1656 xcrybc2 = b2->XCry();
1657 ycrybc2 = b2->YCry();
1658 thetaaxisbc2 = b2->ThetaAxis();
1659 phiaxisbc2 = b2->PhiAxis();
1660 sigietaietabc2 = TMath::Sqrt(b2->CoviEtaiEta());
1661 sigiphiphibc2 = TMath::Sqrt(b2->CoviPhiiPhi());
1662 if (isnan(sigiphiphibc2)) sigiphiphibc2 = -99.;
1663 covietaiphibc2 = b2->CoviEtaiPhi();
1664 if (isnan(covietaiphibc2)) covietaiphibc2 = -99.;
1665 e3x3bc2 = b2->E3x3();
1666 e5x5bc2 = b2->E5x5();
1667 emaxbc2 = b2->EMax();
1668 e2ndbc2 = b2->E2nd();
1669 etopbc2 = b2->ETop();
1670 ebottombc2 = b2->EBottom();
1671 eleftbc2 = b2->ELeft();
1672 erightbc2 = b2->ERight();
1673 e1x3bc2 = b2->E1x3();
1674 e3x1bc2 = b2->E3x1();
1675 e1x5bc2 = b2->E1x5();
1676 e2x2bc2 = b2->E2x2();
1677 e4x4bc2 = b2->E4x4();
1678 e2x5maxbc2 = b2->E2x5Max();
1679 e2x5topbc2 = b2->E2x5Top();
1680 e2x5bottombc2 = b2->E2x5Bottom();
1681 e2x5leftbc2 = b2->E2x5Left();
1682 e2x5rightbc2 = b2->E2x5Right();
1683 xbc2bc2 = b2->Pos().X();
1684 ybc2bc2 = b2->Pos().Y();
1685 zbc2bc2 = b2->Pos().Z();
1686 nhitsbc2= b2->NHits();
1687 }
1688 else {
1689 ebc2 = 0;
1690 etabc2 = 0;
1691 phibc2 = 0;
1692 ietabc2 = 0;
1693 iphibc2 = 0;
1694 ixbc2 = 0;
1695 iybc2 = 0;
1696 etacrybc2 = 0;
1697 phicrybc2 = 0;
1698 xcrybc2 = 0;
1699 ycrybc2 = 0;
1700 thetaaxisbc2 = 0;
1701 phiaxisbc2 = 0;
1702 sigietaietabc2 = 0;
1703 sigiphiphibc2 = 0;
1704 covietaiphibc2 = 0;
1705 e3x3bc2 = 0;
1706 e5x5bc2 = 0;
1707 emaxbc2 = 0;
1708 e2ndbc2 = 0;
1709 etopbc2 = 0;
1710 ebottombc2 = 0;
1711 eleftbc2 = 0;
1712 erightbc2 = 0;
1713 e1x3bc2 = 0;
1714 e3x1bc2 = 0;
1715 e1x5bc2 = 0;
1716 e2x2bc2 = 0;
1717 e4x4bc2 = 0;
1718 e2x5maxbc2 = 0;
1719 e2x5topbc2 = 0;
1720 e2x5bottombc2 = 0;
1721 e2x5leftbc2 = 0;
1722 e2x5rightbc2 = 0;
1723 xbc2bc2 = 0;
1724 ybc2bc2 = 0;
1725 zbc2bc2 = 0;
1726 nhitsbc2 = 0;
1727 }
1728
1729 if (bclast) {
1730 ebclast = bclast->Energy();
1731 etabclast = bclast->Eta();
1732 phibclast = bclast->Phi();
1733 ietabclast = bclast->IEta();
1734 iphibclast = bclast->IPhi();
1735 ixbclast = bclast->IX();
1736 iybclast = bclast->IY();
1737 etacrybclast = bclast->EtaCry();
1738 phicrybclast = bclast->PhiCry();
1739 xcrybclast = bclast->XCry();
1740 ycrybclast = bclast->YCry();
1741 thetaaxisbclast = bclast->ThetaAxis();
1742 phiaxisbclast = bclast->PhiAxis();
1743 sigietaietabclast = TMath::Sqrt(bclast->CoviEtaiEta());
1744 sigiphiphibclast = TMath::Sqrt(bclast->CoviPhiiPhi());
1745 if (isnan(sigiphiphibclast)) sigiphiphibclast = -99.;
1746 covietaiphibclast = bclast->CoviEtaiPhi();
1747 if (isnan(covietaiphibclast)) covietaiphibclast = -99.;
1748 e3x3bclast = bclast->E3x3();
1749 nhitsbclast = bclast->NHits();
1750 }
1751 else {
1752 ebclast = 0;
1753 etabclast = 0;
1754 phibclast = 0;
1755 ietabclast = 0;
1756 iphibclast = 0;
1757 ixbclast = 0;
1758 iybclast = 0;
1759 etacrybclast = 0;
1760 phicrybclast = 0;
1761 xcrybclast = 0;
1762 ycrybclast = 0;
1763 thetaaxisbclast = 0;
1764 phiaxisbclast = 0;
1765 sigietaietabclast = 0;
1766 sigiphiphibclast = 0;
1767 covietaiphibclast = 0;
1768 e3x3bclast = 0;
1769 nhitsbclast = 0;
1770 }
1771
1772 if (bclast2) {
1773 ebclast2 = bclast2->Energy();
1774 etabclast2 = bclast2->Eta();
1775 phibclast2 = bclast2->Phi();
1776 ietabclast2 = bclast2->IEta();
1777 iphibclast2 = bclast2->IPhi();
1778 ixbclast2 = bclast2->IX();
1779 iybclast2 = bclast2->IY();
1780 etacrybclast2 = bclast2->EtaCry();
1781 phicrybclast2 = bclast2->PhiCry();
1782 xcrybclast2 = bclast2->XCry();
1783 ycrybclast2 = bclast2->YCry();
1784 thetaaxisbclast2 = bclast2->ThetaAxis();
1785 phiaxisbclast2 = bclast2->PhiAxis();
1786 sigietaietabclast2 = TMath::Sqrt(bclast2->CoviEtaiEta());
1787 sigiphiphibclast2 = TMath::Sqrt(bclast2->CoviPhiiPhi());
1788 if (isnan(sigiphiphibclast2)) sigiphiphibclast2 = -99.;
1789 covietaiphibclast2 = bclast2->CoviEtaiPhi();
1790 if (isnan(covietaiphibclast2)) covietaiphibclast2 = -99.;
1791 e3x3bclast2 = bclast2->E3x3();
1792 nhitsbclast2 = bclast2->NHits();
1793 }
1794 else {
1795 ebclast2 = 0;
1796 etabclast2 = 0;
1797 phibclast2 = 0;
1798 ietabclast2 = 0;
1799 iphibclast2 = 0;
1800 ixbclast2 = 0;
1801 iybclast2 = 0;
1802 etacrybclast2 = 0;
1803 phicrybclast2 = 0;
1804 xcrybclast2 = 0;
1805 ycrybclast2 = 0;
1806 thetaaxisbclast2 = 0;
1807 phiaxisbclast2 = 0;
1808 sigietaietabclast2 = 0;
1809 sigiphiphibclast2 = 0;
1810 covietaiphibclast2 = 0;
1811 e3x3bclast2 = 0;
1812 nhitsbclast2 = 0;
1813 }
1814
1815 for (UInt_t iclus=0; iclus<NClus; ++iclus) {
1816 if (fillclusterarrays && iclus < s->ClusterSize() ) {
1817 const BasicCluster *ib =s->Cluster(iclus);
1818
1819 ebcs[iclus] = ib->Energy();
1820 etabcs[iclus] = ib->Eta();
1821 phibcs[iclus] = ib->Phi();
1822 ietabcs[iclus] = ib->IEta();
1823 iphibcs[iclus] = ib->IPhi();
1824 ixbcs[iclus] = ib->IX();
1825 iybcs[iclus] = ib->IY();
1826 etacrybcs[iclus] = ib->EtaCry();
1827 phicrybcs[iclus] = ib->PhiCry();
1828 xcrybcs[iclus] = ib->XCry();
1829 ycrybcs[iclus] = ib->YCry();
1830 sigietaietabcs[iclus] = TMath::Sqrt(ib->CoviEtaiEta());
1831 sigiphiphibcs[iclus] = TMath::Sqrt(ib->CoviPhiiPhi());
1832 covietaiphibcs[iclus] = ib->CoviEtaiPhi();
1833 sigetaetabcs[iclus] = TMath::Sqrt(ib->CovEtaEta());
1834 sigphiphibcs[iclus] = TMath::Sqrt(ib->CovPhiPhi());
1835 covetaphibcs[iclus] = ib->CovEtaPhi();
1836 e3x3bcs[iclus] = ib->E3x3();
1837 e5x5bcs[iclus] = ib->E5x5();
1838 emaxbcs[iclus] = ib->EMax();
1839 e2ndbcs[iclus] = ib->E2nd();
1840 etopbcs[iclus] = ib->ETop();
1841 ebottombcs[iclus] = ib->EBottom();
1842 eleftbcs[iclus] = ib->ELeft();
1843 erightbcs[iclus] = ib->ERight();
1844 e1x3bcs[iclus] = ib->E1x3();
1845 e3x1bcs[iclus] = ib->E3x1();
1846 e1x5bcs[iclus] = ib->E1x5();
1847 e2x2bcs[iclus] = ib->E2x2();
1848 e4x4bcs[iclus] = ib->E4x4();
1849 e2x5maxbcs[iclus] = ib->E2x5Max();
1850 e2x5topbcs[iclus] = ib->E2x5Top();
1851 e2x5bottombcs[iclus] = ib->E2x5Bottom();
1852 e2x5leftbcs[iclus] = ib->E2x5Left();
1853 e2x5rightbcs[iclus] = ib->E2x5Right();
1854 nhitsbcs[iclus]= ib->NHits();
1855 }
1856 else {
1857 ebcs[iclus] = -999;
1858 etabcs[iclus] = -999;
1859 phibcs[iclus] = -999;
1860 ietabcs[iclus] = -999;
1861 iphibcs[iclus] = -999;
1862 ixbcs[iclus] = -999;
1863 iybcs[iclus] = -999;
1864 etacrybcs[iclus] = -999;
1865 phicrybcs[iclus] = -999;
1866 xcrybcs[iclus] = -999;
1867 ycrybcs[iclus] = -999;
1868 sigietaietabcs[iclus] = -999;
1869 sigiphiphibcs[iclus] = -999;
1870 covietaiphibcs[iclus] = -999;
1871 sigetaetabcs[iclus] = -999;
1872 sigphiphibcs[iclus] = -999;
1873 covetaphibcs[iclus] = -999;
1874 e3x3bcs[iclus] = -999;
1875 e5x5bcs[iclus] = -999;
1876 emaxbcs[iclus] = -999;
1877 e2ndbcs[iclus] = -999;
1878 etopbcs[iclus] = -999;
1879 ebottombcs[iclus] = -999;
1880 eleftbcs[iclus] = -999;
1881 erightbcs[iclus] = -999;
1882 e1x3bcs[iclus] = -999;
1883 e3x1bcs[iclus] = -999;
1884 e1x5bcs[iclus] = -999;
1885 e2x2bcs[iclus] = -999;
1886 e4x4bcs[iclus] = -999;
1887 e2x5maxbcs[iclus] = -999;
1888 e2x5topbcs[iclus] = -999;
1889 e2x5bottombcs[iclus] = -999;
1890 e2x5leftbcs[iclus] = -999;
1891 e2x5rightbcs[iclus] = -999;
1892 nhitsbcs[iclus] = -999;
1893 }
1894 }
1895
1896 for (UInt_t iclus=0; iclus<NClus; ++iclus) {
1897 if (fillclusterarrays && pfsc && iclus < pfsc->ClusterSize() ) {
1898 const BasicCluster *ib =pfsc->Cluster(iclus);
1899
1900 epfbcs[iclus] = ib->Energy();
1901 etapfbcs[iclus] = ib->Eta();
1902 phipfbcs[iclus] = ib->Phi();
1903 ietapfbcs[iclus] = ib->IEta();
1904 iphipfbcs[iclus] = ib->IPhi();
1905 ixpfbcs[iclus] = ib->IX();
1906 iypfbcs[iclus] = ib->IY();
1907 etacrypfbcs[iclus] = ib->EtaCry();
1908 phicrypfbcs[iclus] = ib->PhiCry();
1909 xcrypfbcs[iclus] = ib->XCry();
1910 ycrypfbcs[iclus] = ib->YCry();
1911 sigietaietapfbcs[iclus] = TMath::Sqrt(ib->CoviEtaiEta());
1912 sigiphiphipfbcs[iclus] = TMath::Sqrt(ib->CoviPhiiPhi());
1913 covietaiphipfbcs[iclus] = ib->CoviEtaiPhi();
1914 sigetaetapfbcs[iclus] = TMath::Sqrt(ib->CovEtaEta());
1915 sigphiphipfbcs[iclus] = TMath::Sqrt(ib->CovPhiPhi());
1916 covetaphipfbcs[iclus] = ib->CovEtaPhi();
1917 e3x3pfbcs[iclus] = ib->E3x3();
1918 e5x5pfbcs[iclus] = ib->E5x5();
1919 emaxpfbcs[iclus] = ib->EMax();
1920 e2ndpfbcs[iclus] = ib->E2nd();
1921 etoppfbcs[iclus] = ib->ETop();
1922 ebottompfbcs[iclus] = ib->EBottom();
1923 eleftpfbcs[iclus] = ib->ELeft();
1924 erightpfbcs[iclus] = ib->ERight();
1925 e1x3pfbcs[iclus] = ib->E1x3();
1926 e3x1pfbcs[iclus] = ib->E3x1();
1927 e1x5pfbcs[iclus] = ib->E1x5();
1928 e2x2pfbcs[iclus] = ib->E2x2();
1929 e4x4pfbcs[iclus] = ib->E4x4();
1930 e2x5maxpfbcs[iclus] = ib->E2x5Max();
1931 e2x5toppfbcs[iclus] = ib->E2x5Top();
1932 e2x5bottompfbcs[iclus] = ib->E2x5Bottom();
1933 e2x5leftpfbcs[iclus] = ib->E2x5Left();
1934 e2x5rightpfbcs[iclus] = ib->E2x5Right();
1935 nhitspfbcs[iclus]= ib->NHits();
1936 }
1937 else {
1938 epfbcs[iclus] = -999;
1939 etapfbcs[iclus] = -999;
1940 phipfbcs[iclus] = -999;
1941 ietapfbcs[iclus] = -999;
1942 iphipfbcs[iclus] = -999;
1943 ixpfbcs[iclus] = -999;
1944 iypfbcs[iclus] = -999;
1945 etacrypfbcs[iclus] = -999;
1946 phicrypfbcs[iclus] = -999;
1947 xcrypfbcs[iclus] = -999;
1948 ycrypfbcs[iclus] = -999;
1949 sigietaietapfbcs[iclus] = -999;
1950 sigiphiphipfbcs[iclus] = -999;
1951 covietaiphipfbcs[iclus] = -999;
1952 sigetaetapfbcs[iclus] = -999;
1953 sigphiphipfbcs[iclus] = -999;
1954 covetaphipfbcs[iclus] = -999;
1955 e3x3pfbcs[iclus] = -999;
1956 e5x5pfbcs[iclus] = -999;
1957 emaxpfbcs[iclus] = -999;
1958 e2ndpfbcs[iclus] = -999;
1959 etoppfbcs[iclus] = -999;
1960 ebottompfbcs[iclus] = -999;
1961 eleftpfbcs[iclus] = -999;
1962 erightpfbcs[iclus] = -999;
1963 e1x3pfbcs[iclus] = -999;
1964 e3x1pfbcs[iclus] = -999;
1965 e1x5pfbcs[iclus] = -999;
1966 e2x2pfbcs[iclus] = -999;
1967 e4x4pfbcs[iclus] = -999;
1968 e2x5maxpfbcs[iclus] = -999;
1969 e2x5toppfbcs[iclus] = -999;
1970 e2x5bottompfbcs[iclus] = -999;
1971 e2x5leftpfbcs[iclus] = -999;
1972 e2x5rightpfbcs[iclus] = -999;
1973 nhitspfbcs[iclus] = -999;
1974 }
1975 }
1976
1977 for (UInt_t iclus=0; iclus<100; ++iclus) {
1978 if (fillclusterarrays && pfsc && iclus < pfsc->NPsClusts() ) {
1979 const PsCluster *ib = pfsc->PsClust(iclus);
1980
1981 epsc[iclus] = ib->Energy();
1982 etapsc[iclus] = ib->Eta();
1983 phipsc[iclus] = ib->Phi();
1984 planepsc[iclus] = ib->PsPlane();
1985 }
1986 else {
1987 epsc[iclus] = -999;
1988 etapsc[iclus] = -999;
1989 phipsc[iclus] = -999;
1990 planepsc[iclus] = 0;
1991 }
1992 }
1993
1994 //initialize photon energy corrections if needed
1995 /*if (!PhotonFix::initialised()) {
1996 PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
1997 }*/
1998
1999 phfixph.setup(e,sceta,scphi,r9);
2000 phfixele.setup(e,sceta,scphi,r9);
2001
2002 const Float_t dval = -99.;
2003 ecor = phfixph.fixedEnergy();
2004 ecorerr = phfixph.sigmaEnergy();
2005 ecorele = phfixele.fixedEnergy();
2006 ecoreleerr = phfixele.sigmaEnergy();
2007 if (phfixph.isbarrel()) {
2008 etac = phfixph.etaC();
2009 etas = phfixph.etaS();
2010 etam = phfixph.etaM();
2011 phic = phfixph.phiC();
2012 phis = phfixph.phiS();
2013 phim = phfixph.phiM();
2014 xz = dval;
2015 xc = dval;
2016 xs = dval;
2017 xm = dval;
2018 yz = dval;
2019 yc = dval;
2020 ys = dval;
2021 ym = dval;
2022 }
2023 else {
2024 etac = dval;
2025 etas = dval;
2026 etam = dval;
2027 phic = dval;
2028 phis = dval;
2029 phim = dval;
2030 xz = phfixph.xZ();
2031 xc = phfixph.xC();
2032 xs = phfixph.xS();
2033 xm = phfixph.xM();
2034 yz = phfixph.yZ();
2035 yc = phfixph.yC();
2036 ys = phfixph.yS();
2037 ym = phfixph.yM();
2038 }
2039
2040 if (c) {
2041 hasconversion = kTRUE;
2042 convp = c->P();
2043 convpt = c->Pt();
2044 conveta = c->Eta();
2045 convphi = c->Phi();
2046 ThreeVector dirconvsc = ThreeVector(s->Point()) - c->Position();
2047 convdeta = c->Eta() - dirconvsc.Eta();
2048 convdphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
2049 convvtxrho = c->Position().Rho();
2050 convvtxz = c->Position().Z();
2051 convvtxphi = c->Position().Phi();
2052
2053 const StableData *leadsd = dynamic_cast<const StableData*>(c->DaughterDat(0));
2054 const StableData *trailsd = dynamic_cast<const StableData*>(c->DaughterDat(1));
2055 if (leadsd->Pt()<trailsd->Pt()) {
2056 const StableData *sdtmp = leadsd;
2057 leadsd = trailsd;
2058 trailsd = sdtmp;
2059 }
2060
2061 const Track *leadtrack = dynamic_cast<const StableParticle*>(leadsd->Original())->Trk();
2062 const Track *trailtrack = dynamic_cast<const StableParticle*>(trailsd->Original())->Trk();
2063
2064 convleadpt = leadsd->Pt();
2065 convtrailpt = trailsd->Pt();
2066 convleadtrackpt = leadtrack->Pt();
2067 convleadtrackalgo = leadtrack->Algo();
2068 if (convleadtrackalgo==29) convleadtrackalgos=2; //gsf track
2069 else if (convleadtrackalgo==15 ||convleadtrackalgo==16) convleadtrackalgos=1; //ecal-seeded track
2070 else convleadtrackalgos = 0; //std iterative track
2071 convleadtrackcharge = leadtrack->Charge();
2072 convtrailtrackpt = trailtrack->Pt();
2073 convtrailtrackalgo = trailtrack->Algo();
2074 if (convtrailtrackalgo==29) convtrailtrackalgos=2; //gsf track
2075 else if (convtrailtrackalgo==15 ||convtrailtrackalgo==16) convtrailtrackalgos=1; //ecal-seeded track
2076 else convtrailtrackalgos = 0; //std iterative track
2077 trailtrackcharge = trailtrack->Charge();
2078 }
2079 else {
2080 hasconversion = kFALSE;
2081 convp = -99.;
2082 convpt = -99.;
2083 conveta = -99.;
2084 convphi = -99.;
2085 convdeta = -99.;
2086 convdphi = -99.;
2087 convvtxrho = -99.;
2088 convvtxz = -999.;
2089 convvtxphi = -99.;
2090 convleadpt = -99.;
2091 convtrailpt = -99.;
2092 convleadtrackpt = -99.;
2093 convleadtrackalgo = -99;
2094 convleadtrackalgos = -99;
2095 convleadtrackcharge = 0;
2096 convtrailtrackpt = -99.;
2097 convtrailtrackalgo = -99;
2098 convtrailtrackalgos = -99;
2099 trailtrackcharge = 0;
2100 }
2101
2102 //electron quantities
2103 if (ele) {
2104 haselectron = kTRUE;
2105 eleisecaldriven = ele->IsEcalDriven();
2106 eleistrackerdriven = ele->IsTrackerDriven();
2107 elee = ele->E();
2108 elept = ele->Pt();
2109 eleeta = ele->Eta();
2110 elephi = ele->Phi();
2111 elecharge = ele->Charge();
2112 elefbrem = ele->FBrem();
2113 eledeta = ele->DeltaEtaSuperClusterTrackAtVtx();
2114 eledphi = ele->DeltaPhiSuperClusterTrackAtVtx();
2115 elep = s->Energy()/ele->ESuperClusterOverP();
2116 elepin = ele->PIn();
2117 elepout = ele->POut();
2118 }
2119 else {
2120 haselectron = kFALSE;
2121 eleisecaldriven = kFALSE;
2122 eleistrackerdriven = kFALSE;
2123 elee = -99.;
2124 elept = -99.;
2125 eleeta = -99.;
2126 elephi = -99.;
2127 elecharge = -99;
2128 elefbrem = -99.;
2129 eledeta = -99.;
2130 eledphi = -99.;
2131 elep = -99.;
2132 elepin = -99.;
2133 elepout = -99.;
2134 }
2135
2136 //pf supercluster quantities
2137 if (pfsc) {
2138 haspfsc = kTRUE;
2139 pfsce = pfsc->Energy();
2140 pfscrawe = pfsc->RawEnergy();
2141 pfsceta = pfsc->Eta();
2142 pfscphi = pfsc->Phi();
2143 pfscnclusters = pfsc->NClusters();
2144 pfscnhits = pfsc->NHits();
2145 pfscetawidth = pfsc->EtaWidth();
2146 pfscphiwidth = pfsc->PhiWidth();
2147 pfscnpsclusters = pfsc->NPsClusts();
2148 }
2149 else {
2150 haspfsc = kFALSE;
2151 pfsce = -99.;
2152 pfscrawe = -99.;
2153 pfsceta = -99.;
2154 pfscphi = -99.;
2155 pfscnclusters = 0;
2156 pfscnhits = 0;
2157 pfscetawidth = -99.;
2158 pfscphiwidth = -99.;
2159 pfscnpsclusters = 0;
2160 }
2161
2162 genz = -99.;
2163 if (m) {
2164 ispromptgen = kTRUE;
2165 gene = m->E();
2166 genpt = m->Pt();
2167 geneta = m->Eta();
2168 genphi = m->Phi();
2169 const MCParticle *mm = m->DistinctMother();
2170 if (mm) genz = mm->DecayVertex().Z();
2171 pdgid = m->PdgId();
2172 if (mm) motherpdgid = mm->PdgId();
2173 else motherpdgid = -99;
2174 }
2175 else {
2176 ispromptgen = kFALSE;
2177 gene = -99.;
2178 genpt = -99.;
2179 geneta = -99.;
2180 genphi = -99.;
2181 pdgid = -99;
2182 motherpdgid = -99;
2183 }
2184 }
2185
2186 void PhotonTreeWriterVtx::SetVars(const Vertex *v, const Photon *p1, const Photon *p2, const PFCandidateCol *pfcands, Int_t idx, Int_t numvtx, Float_t genvtxz)
2187 {
2188
2189 //printf("start\n");
2190
2191 n = idx;
2192 nvtx = numvtx;
2193 zgen = genvtxz;
2194
2195 x = v->X();
2196 y = v->Y();
2197 z = v->Z();
2198
2199 Double_t dsumpt = 0.;
2200 Double_t dsumptsq = 0.;
2201
2202 nchalltoward = 0;
2203 nchalltransverse = 0;
2204 nchallaway = 0;
2205 nchcuttoward = 0;
2206 nchcuttransverse = 0;
2207 nchcutaway = 0;
2208
2209 ThreeVector vtxmom;
2210
2211 //printf("mom\n");
2212 FourVectorM diphoton = p1->MomVtx(v->Position()) + p2->MomVtx(v->Position());
2213 //printf("done mom\n");
2214 ptgg = diphoton.Pt();
2215 phigg = diphoton.Phi();
2216 etagg = diphoton.Eta();
2217 mgg = diphoton.M();
2218 pxgg = diphoton.Px();
2219 pygg = diphoton.Py();
2220 pzgg = diphoton.Pz();
2221
2222 //printf("loop\n");
2223
2224 for (UInt_t i = 0; i<pfcands->GetEntries(); ++i) {
2225 const PFCandidate *pfc = pfcands->At(i);
2226 if (pfc->PFType()!=PFCandidate::eHadron || !pfc->HasTrackerTrk()) continue;
2227 if (TMath::Abs( pfc->TrackerTrk()->DzCorrected(*v) ) > 0.2) continue;
2228 if (TMath::Abs( pfc->TrackerTrk()->D0Corrected(*v) ) > 0.1) continue;
2229
2230 vtxmom += ThreeVector(pfc->Px(),pfc->Py(),pfc->Pz());
2231
2232 dsumpt += pfc->Pt();
2233 dsumptsq += pfc->Pt()*pfc->Pt();
2234
2235 Double_t dphi = TMath::Abs(MathUtils::DeltaPhi(*pfc,diphoton));
2236 if (dphi<(TMath::Pi()/3.0)) {
2237 ++nchalltoward;
2238 if (pfc->Pt()>0.5) ++nchcuttoward;
2239 }
2240 else if (dphi>(2.0*TMath::Pi()/3.0)) {
2241 ++nchallaway;
2242 if (pfc->Pt()>0.5) ++nchcutaway;
2243 }
2244 else {
2245 ++nchalltransverse;
2246 if (pfc->Pt()>0.5) ++nchcuttransverse;
2247 }
2248
2249 }
2250
2251 //printf("doneloop\n");
2252
2253
2254 sumpt = dsumpt;
2255 sumptsq = dsumptsq;
2256
2257 pt = vtxmom.Rho();
2258 phi = vtxmom.Phi();
2259 eta = vtxmom.Eta();
2260 px = vtxmom.X();
2261 py = vtxmom.Y();
2262 pz = vtxmom.Z();
2263
2264 //printf("done\n");
2265
2266 return;
2267
2268 }
2269
2270
2271 //_____________________________________________________________________________
2272 void PhotonTreeWriter::ApplyLeptonTag(const Photon *phHard,
2273 const Photon *phSoft,
2274 const Vertex *selvtx)
2275 {
2276 // perform flavor-based lepton tagging (used before the legacy paper of 2013)
2277 // the diphoton event record will have one more entry; i.e. leptonTag
2278 // leptonTag = -1 -> lepton-taggng was swicthed off
2279 // = 0 -> event tagged as 'non-lepton-event'
2280 // = +1 -> event tagged as muon-event
2281 // = +2 -> event tagged as electron-event
2282 fDiphotonEvent->leptonTag = 0;
2283 Int_t closestVtx = 0;
2284 const Muon *muon = GetLeptonTagMuon(phHard, phSoft);
2285 if (muon) {
2286 // muon tagged
2287 fDiphotonEvent->leptonTag = 2;
2288 SetLeptonTagMuonVars(phHard, phSoft, muon);
2289 }
2290
2291 const Electron *electron = 0;
2292 if (fDiphotonEvent->leptonTag < 1) {
2293 electron = GetLeptonTagElectron(phHard, phSoft);
2294 }
2295
2296 if (electron &&
2297 MassOfPairIsWithinWindowAroundMZ(phHard, electron, 10) == false &&
2298 MassOfPairIsWithinWindowAroundMZ(phSoft, electron, 10) == false) {
2299 // electron tagged
2300 fDiphotonEvent->leptonTag = 1;
2301
2302 fDiphotonEvent-> elePt = fLeptonTagElectrons->At(0)->Pt();
2303 fDiphotonEvent-> eleEta = fLeptonTagElectrons->At(0)->Eta();
2304 fDiphotonEvent-> elePhi = fLeptonTagElectrons->At(0)->Phi();
2305 fDiphotonEvent-> eleSCEta = fLeptonTagElectrons->At(0)->SCluster()->Eta();
2306 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;
2307
2308 fDiphotonEvent-> eleIso2 = -99.;
2309
2310 if ( fDoSynching ) {
2311 Double_t distVtx = 999.0;
2312 for(UInt_t nv=0; nv<fPV->GetEntries(); nv++){
2313 double dz = TMath::Abs(fLeptonTagElectrons->At(0)->GsfTrk()->DzCorrected(*fPV->At(nv)));
2314 if(dz < distVtx) {
2315 distVtx = dz;
2316 closestVtx = nv;
2317 }
2318 }
2319 fDiphotonEvent-> eleIdMva = fElectronIDMVA->MVAValue(fLeptonTagElectrons->At(0), fPV->At(closestVtx));
2320 }
2321
2322
2323 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;
2324 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;
2325 fDiphotonEvent-> eleDist = fLeptonTagElectrons->At(0)->ConvPartnerDist();
2326 fDiphotonEvent-> eleDcot = fLeptonTagElectrons->At(0)->ConvPartnerDCotTheta();
2327 fDiphotonEvent-> eleCoviee = fLeptonTagElectrons->At(0)->CoviEtaiEta();
2328 fDiphotonEvent-> eleDphiin = TMath::Abs(fLeptonTagElectrons->At(0)->DeltaPhiSuperClusterTrackAtVtx());
2329 fDiphotonEvent-> eleDetain = TMath::Abs(fLeptonTagElectrons->At(0)->DeltaEtaSuperClusterTrackAtVtx());
2330 fDiphotonEvent-> eleDR1 = MathUtils::DeltaR(fLeptonTagElectrons->At(0),phHard);
2331 fDiphotonEvent-> eleDR2 = MathUtils::DeltaR(fLeptonTagElectrons->At(0),phSoft);
2332 fDiphotonEvent-> eleMass1 = (phHard->Mom()+fLeptonTagElectrons->At(0)->Mom()).M();
2333 fDiphotonEvent-> eleMass2 = (phSoft->Mom()+fLeptonTagElectrons->At(0)->Mom()).M();
2334 fDiphotonEvent-> eleNinnerHits = fLeptonTagElectrons->At(0)->Trk()->NExpectedHitsInner();
2335 } // electron tagged
2336
2337 if(false){
2338 if(fDiphotonEvent->evt==79737729 || fDiphotonEvent->evt== 871378986 || fDiphotonEvent->evt==528937923 || fDiphotonEvent->evt== 261543921){
2339 printf("ming sync check ele: 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,fDiphotonEvent->mass,fDiphotonEvent->elePt,fDiphotonEvent->eleEta,fDiphotonEvent->eleSCEta,closestVtx);
2340 //return;
2341 }
2342 if(fDiphotonEvent->evt==333643114 || fDiphotonEvent->evt==89022540 || fDiphotonEvent->evt==8983064 || fDiphotonEvent->evt==876316897 || fDiphotonEvent->evt==541603559 || fDiphotonEvent->evt==223740859) {
2343 printf("ming sync check muon: run:%d evt:%d lumi:%d leptonTag:%d numMuons:%d mass:%f\n muonPt:%f muonEta:%f\n\n",fDiphotonEvent->run,fDiphotonEvent->evt,fDiphotonEvent->lumi,fDiphotonEvent->leptonTag,fLeptonTagMuons->GetEntries(),fDiphotonEvent->mass,fDiphotonEvent->muonPt,fDiphotonEvent->muonEta);
2344 //return;
2345 }
2346 }
2347 } // void PhotonTreeWriter::ApplyLeptonTag(..)
2348
2349
2350 //_____________________________________________________________________________
2351 const Muon* PhotonTreeWriter::GetLeptonTagMuon(const Photon *phHard,
2352 const Photon *phSoft)
2353 {
2354 // need to have dR > 1 for with respect to both photons ***changed to 0.7 for 2012
2355 if (fVerbosityLevel > 90) {
2356 cout << "++ Entering GetLeptonTagMuon ..." << endl
2357 << "++ fLeptonTagMuons->GetEntries(): "
2358 << fLeptonTagMuons->GetEntries() << endl;
2359 const Muon *muon = 0;
2360 if (fLeptonTagMuons->GetEntries() > 0) {
2361 muon = fLeptonTagMuons->At(0);
2362 }
2363 if (muon) {
2364 cout << "++ fLeptonTagMuons->At(0): "
2365 << fLeptonTagMuons->At(0) << endl
2366 << "++ MathUtils::DeltaR(fLeptonTagMuons->At(0), phHard): "
2367 << MathUtils::DeltaR(fLeptonTagMuons->At(0), phHard) << endl
2368 << "++ MathUtils::DeltaR(fLeptonTagMuons->At(0), phSoft): "
2369 << MathUtils::DeltaR(fLeptonTagMuons->At(0), phSoft) << endl
2370 << "++ muon pt :" << muon->Pt () << endl
2371 << "++ muon eta:" << muon->Eta() << endl
2372 << "++ muon phi:" << muon->Phi() << endl;
2373 }
2374 cout << "++ Exiting GetLeptonTagMuon ..." << endl;
2375 }
2376
2377 for (unsigned int imu=0; imu<fLeptonTagMuons->GetEntries(); ++imu) {
2378 const Muon *mu = fLeptonTagMuons->At(imu);
2379 if (fLeptonTagMuons->GetEntries() > 0 &&
2380 fLeptonTagMuons->At(0) != 0 &&
2381 MathUtils::DeltaR(mu, phHard) >= 1.0 &&
2382 MathUtils::DeltaR(mu, phSoft) >= 1.0) {
2383 return mu;
2384 }
2385 }
2386
2387 return 0;
2388
2389 } // const Muon* PhotonTreeWriter::GetLeptonTagMuon(..)
2390
2391
2392 //_____________________________________________________________________________
2393 void PhotonTreeWriter::SetLeptonTagMuonVars(const Photon *phHard,
2394 const Photon *phSoft,
2395 const Muon *muon)
2396 {
2397 fDiphotonEvent-> muonPt = muon->Pt();
2398 fDiphotonEvent-> muonEta = muon->Eta();
2399 fDiphotonEvent-> muonPhi = muon->Phi();
2400 fDiphotonEvent-> muDR1 = MathUtils::DeltaR(muon, phHard);
2401 fDiphotonEvent-> muDR2 = MathUtils::DeltaR(muon, phSoft);
2402
2403 Float_t combinedIso = (muon->IsoR03SumPt() +
2404 muon->IsoR03EmEt() +
2405 muon->IsoR03HadEt());
2406
2407 Float_t coneArea = TMath::Pi() * 0.3 * 0.3;
2408
2409 Float_t rho1 = fPileUpDen->At(0)->RhoRandomLowEta();
2410 Float_t rho2 = fPileUpDen->At(0)->RhoRandom();
2411 Float_t rho3 = fPileUpDen->At(0)->RhoLowEta();
2412 Float_t rho4 = fPileUpDen->At(0)->Rho();
2413
2414 fDiphotonEvent-> muIso1 = (combinedIso - rho1 * coneArea) / muon->Pt();
2415 fDiphotonEvent-> muIso2 = (combinedIso - rho2 * coneArea) / muon->Pt();
2416 fDiphotonEvent-> muIso3 = (combinedIso - rho3 * coneArea) / muon->Pt();
2417 fDiphotonEvent-> muIso4 = (combinedIso - rho4 * coneArea) / muon->Pt();
2418
2419 fDiphotonEvent-> muD0 = TMath::Abs(muon->BestTrk()->D0Corrected(*fPV->At(0)));
2420 fDiphotonEvent-> muDZ = TMath::Abs(muon->BestTrk()->DzCorrected(*fPV->At(0)));
2421 fDiphotonEvent-> muChi2 = muon->GlobalTrk()->Chi2()/muon->GlobalTrk()->Ndof();
2422
2423 fDiphotonEvent-> muNhits = muon->BestTrk()->NHits();
2424 fDiphotonEvent-> muNpixhits = muon->BestTrk()->NPixelHits();
2425 fDiphotonEvent-> muNegs = muon->NSegments();
2426 fDiphotonEvent-> muNMatch = muon->NMatches();
2427 } // void PhotonTreeWriter::SetLeptonTagMuonVars(..)
2428
2429
2430 //_____________________________________________________________________________
2431 const Electron* PhotonTreeWriter::GetLeptonTagElectron(const Photon *phHard,
2432 const Photon *phSoft)
2433 {
2434 // need to have dR > 1 for with respect to both photons ***changed to 0.7 for 2012
2435 if (fVerbosityLevel > 90) {
2436 cout << "++ Entering GetLeptonTagElectron ..." << endl
2437 << "++ fLeptonTagElectrons->GetEntries(): "
2438 << fLeptonTagElectrons->GetEntries() << endl
2439 << "++ PhotonTools::ElectronVetoCiC(phHard, fLeptonTagElectrons): "
2440 << PhotonTools::ElectronVetoCiC(phHard, fLeptonTagElectrons) << endl
2441 << "++ PhotonTools::ElectronVetoCiC(phSoft, fLeptonTagElectrons): "
2442 << PhotonTools::ElectronVetoCiC(phSoft, fLeptonTagElectrons) << endl
2443 << "++ PhotonTools::ElectronVetoCiC(phHard, fElectrons): "
2444 << PhotonTools::ElectronVetoCiC(phHard, fElectrons) << endl
2445 << "++ PhotonTools::ElectronVetoCiC(phSoft, fElectrons): "
2446 << PhotonTools::ElectronVetoCiC(phSoft, fElectrons) << endl;
2447 if (fLeptonTagElectrons->GetEntries() > 0 && fLeptonTagElectrons->At(0) != 0) {
2448 cout << "++ fLeptonTagElectrons->At(0): "
2449 << fLeptonTagElectrons->At(0) << endl
2450 << "++ MathUtils::DeltaR(fLeptonTagElectrons->At(0), phHard): "
2451 << MathUtils::DeltaR(fLeptonTagElectrons->At(0), phHard) << endl
2452 << "++ MathUtils::DeltaR(fLeptonTagElectrons->At(0), phSoft): "
2453 << MathUtils::DeltaR(fLeptonTagElectrons->At(0), phSoft) << endl
2454 << "++ electron pt : " << fLeptonTagElectrons->At(0)->Pt () << endl
2455 << "++ electron eta: " << fLeptonTagElectrons->At(0)->Eta() << endl
2456 << "++ electron phi: " << fLeptonTagElectrons->At(0)->Phi() << endl;
2457 }
2458 cout << "++ Exiting GetLeptonTagElectron ..." << endl;
2459 }
2460
2461 //loop over all electrons
2462 for (unsigned int iele=0; iele<fLeptonTagElectrons->GetEntries(); ++iele) {
2463 const Electron *ele = fLeptonTagElectrons->At(iele);
2464 if (fLeptonTagElectrons->GetEntries() > 0 &&
2465 fLeptonTagElectrons->At(0) != 0 &&
2466 PhotonTools::ElectronVetoCiC(phHard, fLeptonTagElectrons) >= 1 &&
2467 PhotonTools::ElectronVetoCiC(phSoft, fLeptonTagElectrons) >= 1 &&
2468 PhotonTools::ElectronVetoCiC(phHard, fElectrons) >= 1 &&
2469 PhotonTools::ElectronVetoCiC(phSoft, fElectrons) >= 1 &&
2470 MathUtils::DeltaR(ele, phHard) >= 1 &&
2471 MathUtils::DeltaR(ele, phSoft) >= 1){
2472 return ele;
2473 }
2474 }
2475
2476 return 0;
2477
2478 } // const Electron* PhotonTreeWriter::GetLeptonTagElectron(..)
2479
2480
2481 //_____________________________________________________________________________
2482 bool PhotonTreeWriter::MassOfPairIsWithinWindowAroundMZ(
2483 const Particle * particle1,
2484 const Particle * particle2,
2485 Float_t halfWindowSize,
2486 Float_t MZ
2487 ) {
2488 Float_t mass = (particle1->Mom() + particle2->Mom()).M();
2489 return TMath::Abs(mass - MZ) < halfWindowSize;
2490 } // bool PhotonTreeWriter::MassOfPairIsWithinWindowAroundMZ(..)
2491
2492
2493 //_____________________________________________________________________________
2494 void PhotonTreeWriter::ApplyVHLepTag(const Photon *phHard,
2495 const Photon *phSoft,
2496 const Vertex *selvtx)
2497 {
2498
2499 // Perform flavor-based lepton tagging (used since the legacy paper of 2013)
2500 // the diphoton event record will have one more entry; i.e. leptonTag
2501 // VHLepTag = -1 -> lepton-tagging was swicthed off
2502 // = 0 -> event tagged as 'non-lepton event'
2503 // = +1 -> event tagged as a low-MET low-S/sqrt(B) event
2504 // = +2 -> event tagged as a high-MET/dilepton high-S/sqrt(B) event
2505 // = +3 -> event tagged as both a loose and a tight VH leptonic event
2506 // TODO: Set the selected vertex to the lepton vertex if tagged as a
2507 // VH(lep) event.
2508
2509 //printf("ApplyVHLepTag: electrons = %i, softelectrons = %i, muons = %i, softmuons = %i\n",fLeptonTagElectrons->GetEntries(), fLeptonTagSoftElectrons->GetEntries(), fLeptonTagMuons->GetEntries(), fLeptonTagSoftMuons->GetEntries());
2510
2511 if (fVerbosityLevel > 90) {
2512 cout << "+ Entering ApplyVHLepTag ..." << endl
2513 << "+ fElectrons->GetEntries(): " << fElectrons->GetEntries() << endl
2514 << "+ fLeptonTagElectrons->GetEntries(): "
2515 << fLeptonTagElectrons->GetEntries() << endl
2516 << "+ fLeptonTagMuons->GetEntries(): "
2517 << fLeptonTagMuons->GetEntries() << endl
2518 << "+ corrpfmet: " << fDiphotonEvent->corrpfmet << endl;
2519 }
2520
2521 fDiphotonEvent->VHLepTag = 0; // non-lepton event
2522 bool isVHLepLoose = false;
2523 bool isVHLepTight = false;
2524
2525 if (VHLepHasDielectron(phHard, phSoft) ||
2526 VHLepHasDimuon (phHard, phSoft)) isVHLepTight = true;
2527
2528 if (fDiphotonEvent->leptonTag < 0) {
2529 ApplyLeptonTag(phHard, phSoft, selvtx);
2530 }
2531
2532 //printf("check dilepton, tight = %i\n",int(isVHLepTight));
2533
2534 const Muon *muon = GetLeptonTagMuon(phHard, phSoft);
2535 if (muon && VHLepNumberOfJets(phHard, phSoft, selvtx, muon) <= 2) {
2536 if (fDiphotonEvent->corrpfmet > 45.) isVHLepTight = true; // high MET event
2537 else isVHLepLoose = true; // low MET event
2538 } // Found a good VH(lep) tag muon.
2539
2540 //printf("check muon, tight = %i, loose = %i\n",int(isVHLepTight),int(isVHLepLoose));
2541
2542 const Electron *electron = GetLeptonTagElectron(phHard, phSoft);
2543 if (electron && VHLepNumberOfJets(phHard, phSoft, selvtx, electron) <= 2) {
2544 if (fDiphotonEvent->corrpfmet > 45.) {
2545 isVHLepTight = true; // High MET event.
2546 } else if (!MassOfPairIsWithinWindowAroundMZ(phHard, electron, 10) &&
2547 !MassOfPairIsWithinWindowAroundMZ(phSoft, electron, 10)) {
2548 isVHLepLoose = true; // low MET event
2549 } // Low MET event.
2550 } // Found a good VH(lep) tag electron.
2551
2552 //printf("check electron, tight = %i, loose = %i\n",int(isVHLepTight),int(isVHLepLoose));
2553
2554 if (isVHLepLoose) fDiphotonEvent->VHLepTag += 1;
2555 if (isVHLepTight) fDiphotonEvent->VHLepTag += 2;
2556
2557 if (fVerbosityLevel > 90) {
2558 cout << "+ muon: " << muon << endl;
2559 if (muon) {
2560 cout << "+ njets(mu): "
2561 << VHLepNumberOfJets(phHard, phSoft, selvtx, muon) << endl;
2562 }
2563 cout << "+ electron: " << electron << endl;
2564 if (electron) {
2565 cout << "+ njets(ele): "
2566 << VHLepNumberOfJets(phHard, phSoft, selvtx, electron) << endl
2567 << "+ MassOfPairIsWithinWindowAroundMZ(phHard, electron, 10)"
2568 << MassOfPairIsWithinWindowAroundMZ(phHard, electron, 10) << endl
2569 << "+ MassOfPairIsWithinWindowAroundMZ(phSoft, electron, 10)"
2570 << MassOfPairIsWithinWindowAroundMZ(phSoft, electron, 10) << endl;
2571 }
2572 cout << "+ isVHLepTight: " << isVHLepTight << endl
2573 << "+ isVHLepLoose: " << isVHLepLoose << endl
2574 << "+ VHLepTag: " << fDiphotonEvent->VHLepTag << endl
2575 << "+ Exiting ApplyVHLepTag ... " << endl;
2576 } /// fVerbosityLevel > 90
2577
2578 } // void PhotonTreeWriter::ApplyVHLepTag(..)
2579
2580
2581 //_____________________________________________________________________________
2582 void PhotonTreeWriter::ApplyVHHadTag(const Photon *phHard,
2583 const Photon *phSoft,
2584 const Vertex *selvtx,
2585 const Jet *jet1,
2586 const Jet *jet2)
2587 {
2588
2589 // Perform VH hadronic tagging (used since the legacy paper of 2013)
2590 // the diphoton event record will have one more entry; i.e. VHHadTag
2591 // VHHadTag = -1 -> VH(had)-tagging was swicthed off
2592 // = 0 -> event tagged as 'non-VH(had) event'
2593 // = +1 -> event tagged as a VH(had) event
2594 // TODO: Should the selected vertex be updated using the vertex of the jets?
2595
2596 fDiphotonEvent->VHHadTag = 0; // non-VH(had) event
2597
2598 // Calculate |cos(theta*)|
2599 // Inspired by Globe:
2600 // https://github.com/h2gglobe/h2gglobe/blob/ae4356ac0d18e6f77da6e0420ab6f5168e353315/PhotonAnalysis/src/PhotonAnalysis.cc#L4369-L4377
2601 FourVectorM pgg = phHard->Mom() + phSoft->Mom(); // diphoton 4-vector
2602 FourVectorM pjj = jet1 ->Mom() + jet2 ->Mom(); // dijet 4-vector
2603 FourVectorM p4 = pjj + pgg; // 4-body 4-vector
2604 TLorentzVector Vstar(p4.X(), p4.Y(), p4.Z(), p4.T());
2605 TLorentzVector H_Vstar(pgg.X(), pgg.Y(), pgg.Z(), pgg.T());
2606 H_Vstar.Boost(-Vstar.BoostVector());
2607 double cosThetaStar = -H_Vstar.CosTheta();
2608 double absCosThetaStar = TMath::Abs(cosThetaStar);
2609
2610 if (fDoSynching) {
2611 fDiphotonEvent->costhetastar = cosThetaStar;
2612 }
2613
2614 Float_t &mass = fDiphotonEvent->mass;
2615 Float_t reducedMass = mass / 120.;
2616 UInt_t nJets = NumberOfJets(phHard, phSoft, selvtx, 40., 2.4);
2617 Float_t mjj = fDiphotonEvent->dijetmass;
2618
2619 /// See L2007-2013 of the Hgg AN 2013/253 v3
2620 if (fDiphotonEvent->ptgg > 130. * reducedMass &&
2621 nJets >= 2 &&
2622 60. < mjj && mjj < 120. &&
2623 absCosThetaStar < 0.5) {
2624 // Tag this event as a VH-hadronic one!
2625 fDiphotonEvent->VHHadTag = 1;
2626 }
2627 } // void PhotonTreeWriter::ApplyVHHadTag(..)
2628
2629
2630 //_____________________________________________________________________________
2631 bool PhotonTreeWriter::VHLepHasDielectron(const Photon *phHard,
2632 const Photon *phSoft) {
2633 if (fLeptonTagSoftElectrons->GetEntries() < 2) return false;
2634
2635 if (PhotonTools::ElectronVetoCiC(phHard, fLeptonTagSoftElectrons) < 1 ||
2636 PhotonTools::ElectronVetoCiC(phSoft, fLeptonTagSoftElectrons) < 1 ||
2637 PhotonTools::ElectronVetoCiC(phHard, fElectrons) < 1 ||
2638 PhotonTools::ElectronVetoCiC(phSoft, fElectrons) < 1) {
2639 return false;
2640 }
2641
2642 vector<UInt_t> goodElectrons;
2643
2644 // Loop over electrons.
2645 for (UInt_t iele=0; iele < fLeptonTagSoftElectrons->GetEntries(); ++iele){
2646 const Electron *ele = fLeptonTagSoftElectrons->At(iele);
2647 if (MathUtils::DeltaR(ele, phHard) < 0.5) continue;
2648 if (MathUtils::DeltaR(ele, phSoft) < 0.5) continue;
2649 goodElectrons.push_back(iele);
2650 }
2651
2652 if (goodElectrons.size() < 2) return false;
2653
2654 // Loop over pairs of selected electrons.
2655 for (UInt_t iiele1 = 0; iiele1 < goodElectrons.size() - 1; ++iiele1) {
2656 UInt_t iele1 = goodElectrons[iiele1];
2657 const Electron *ele1 = fLeptonTagSoftElectrons->At(iele1);
2658 for (UInt_t iiele2 = iiele1 + 1; iiele2 < goodElectrons.size(); ++iiele2) {
2659 UInt_t iele2 = goodElectrons[iiele2];
2660 const Electron *ele2 = fLeptonTagSoftElectrons->At(iele2);
2661 Double_t mass12 = (ele1->Mom() + ele2->Mom()).M();
2662 if (mass12 < 70. || 110. < mass12) continue;
2663 if (fVerbosityLevel > 0) {
2664 cout << " Found a tagging dielectron!" << endl << flush;
2665 }
2666 fDiphotonEvent->ele1Pt = ele1->Pt ();
2667 fDiphotonEvent->ele1Eta = ele1->Eta();
2668 fDiphotonEvent->ele1Phi = ele1->Phi();
2669 fDiphotonEvent->ele2Pt = ele2->Pt ();
2670 fDiphotonEvent->ele2Eta = ele2->Eta();
2671 fDiphotonEvent->ele2Phi = ele2->Phi();
2672 return true;
2673 }
2674 }
2675
2676 return false;
2677
2678 } // bool PhotonTreeWriter::VHLepHasDielectron(..)
2679
2680
2681 //_____________________________________________________________________________
2682 bool PhotonTreeWriter::VHLepHasDimuon(const Photon *phHard,
2683 const Photon *phSoft) {
2684 if (fLeptonTagSoftMuons->GetEntries() < 2) return false;
2685
2686 if (fVerbosityLevel > 0) {
2687 cout << "PhotonTreeWriter::VHLepHasDimuon: Found "
2688 << fLeptonTagSoftMuons->GetEntries() << " muons!" << endl;
2689 }
2690
2691 vector<UInt_t> goodMuons;
2692
2693 // Loop over muons and apply the cut Delta R (mu, pho) > 0.5 .
2694 for (UInt_t imu=0; imu < fLeptonTagSoftMuons->GetEntries(); ++imu){
2695 const Muon *mu = fLeptonTagSoftMuons->At(imu);
2696 if (MathUtils::DeltaR(mu, phHard) < 0.5) continue;
2697 if (MathUtils::DeltaR(mu, phSoft) < 0.5) continue;
2698 goodMuons.push_back(imu);
2699 }
2700
2701 if (goodMuons.size() < 2) return false;
2702
2703 // Loop over muon pairs of selected muons and apply the cut 70 < mass(mu1, mu2) < 110.
2704 for (UInt_t iimu1 = 0; iimu1 < goodMuons.size() - 1; ++iimu1) {
2705 UInt_t imu1 = goodMuons[iimu1];
2706 const Muon *mu1 = fLeptonTagSoftMuons->At(imu1);
2707 for (UInt_t iimu2 = iimu1 + 1; iimu2 < goodMuons.size(); ++iimu2) {
2708 UInt_t imu2 = goodMuons[iimu2];
2709 const Muon *mu2 = fLeptonTagSoftMuons->At(imu2);
2710 Double_t mass12 = (mu1->Mom() + mu2->Mom()).M();
2711 if (mass12 < 70. || 110. < mass12) continue;
2712 if (fVerbosityLevel > 0) {
2713 cout << " Found a tagging dimoun!" << endl << flush;
2714 }
2715 fDiphotonEvent->mu1Pt = mu1->Pt ();
2716 fDiphotonEvent->mu1Eta = mu1->Eta();
2717 fDiphotonEvent->mu1Phi = mu1->Phi();
2718 fDiphotonEvent->mu2Pt = mu2->Pt ();
2719 fDiphotonEvent->mu2Eta = mu2->Eta();
2720 fDiphotonEvent->mu2Phi = mu2->Phi();
2721 return true;
2722 }
2723 }
2724
2725 return false;
2726
2727 } // PhotonTreeWriter::VHLepHasDimuon(..)
2728
2729
2730 //_____________________________________________________________________________
2731 UInt_t PhotonTreeWriter::VHLepNumberOfJets(const Photon *phHard,
2732 const Photon *phSoft,
2733 const Vertex *selvtx,
2734 const Particle *lepton) {
2735
2736 UInt_t nJets = 0;
2737
2738 // Loop over jets, count those passing selection
2739 // Use same ID as for the tth tag
2740 for(UInt_t ijet=0; ijet < fPFJets->GetEntries(); ++ijet){
2741 const Jet *jet = fPFJets->At(ijet);
2742 // Apply jet selection, see L116 and L125 of the AN
2743 if (jet->Pt() < 20. || jet->AbsEta() > 2.4) continue;
2744 // Apply the cut Delta R(photon, jet) < 0.5.
2745 if (MathUtils::DeltaR(jet, phHard) < 0.5) continue;
2746 if (MathUtils::DeltaR(jet, phSoft) < 0.5) continue;
2747 // Apply the cut Delta R(photon, lepton) < 0.5.
2748 if (MathUtils::DeltaR(jet, lepton) < 0.5) continue;
2749 // Make sure we have a PF jet
2750 const PFJet *pfjet = dynamic_cast<const PFJet*>(jet);
2751 if (!pfjet) continue;
2752 if (!JetTools::passPFLooseId(pfjet)) continue;
2753 // Apply the jet ID / pileup removal as given in Table 4
2754 Double_t betaStar = JetTools::betaStarClassic(pfjet, selvtx, fPV);
2755 if (betaStar > 0.2 * log(fPV->GetEntries() - 0.64)) continue;
2756 if (JetTools::dR2Mean(pfjet, -1) > 0.065) continue;
2757 // this jet passes, count it in
2758 ++nJets;
2759 } // End of loop over jets
2760
2761 return nJets;
2762
2763 } // PhotonTreeWriter::VHLepNumberOfJets(..)
2764
2765
2766 //_____________________________________________________________________________
2767 // Applies the ttH tag given precelected leading and trailing photons
2768 // phHard and phSoft and the corresponding (pre?) selected vertex selvtx.
2769 // The result is stored as an integer value of the tthTag variable
2770 // entering the diphoton event record.
2771 void PhotonTreeWriter::ApplyTTHTag(const Photon *phHard,
2772 const Photon *phSoft,
2773 const Vertex *selvtx)
2774 {
2775 // ttH tag = -1 .. not applied
2776 // 0 .. not tagged
2777 // 1 .. tagged as a hadronic ttH event
2778 // 2 .. tagged as a leptonic ttH event
2779 // 3 .. tagged as both a hadronic and a leptonic ttH event
2780
2781 if (fVerbosityLevel > 90) {
2782 cout << "+ Entering ApplyTTHTagTag ..." << endl;
2783 cout << "+ fLeptonTagElectrons->GetEntries(): " << fLeptonTagElectrons->GetEntries() << endl;
2784 cout << "+ fLeptonTagMuons->GetEntries(): " << fLeptonTagMuons->GetEntries() << endl;
2785 }
2786
2787
2788 //printf("ApplyTTHLepTag: electrons = %i, softelectrons = %i, muons = %i, softmuons = %i\n",fLeptonTagElectrons->GetEntries(), fLeptonTagSoftElectrons->GetEntries(), fLeptonTagMuons->GetEntries(), fLeptonTagSoftMuons->GetEntries());
2789
2790
2791 fDiphotonEvent->tthTag = 0;
2792
2793 // Selection taken from the AN2012_480_V6 of 24 April 2013 further
2794 // refferred to as "the AN"
2795 // And from "the Hgg AN"
2796 // http://cms.cern.ch/iCMS/jsp/openfile.jsp?tp=draft&files=AN2013_253_v3.pdf
2797
2798 const Particle *lepton = TTHSelectLepton(phHard, phSoft, selvtx);
2799
2800 //printf("lepton = %p\n",(void*)lepton);
2801
2802 // Init jet object counters
2803 UInt_t nJets = 0;
2804 UInt_t nBJets = 0;
2805
2806 // Loop over jets, count those passing selection.
2807 for(UInt_t ijet=0; ijet < fPFJets->GetEntries(); ++ijet){
2808 const Jet *jet = fPFJets->At(ijet);
2809 // Apply jet selection, see L116 and L125 of the AN
2810 if (jet->Pt() < 25. || jet->AbsEta() > 2.4) continue;
2811 // Apply Delta R(photon, jet), see email from Francesco Micheli
2812 // sent 31 Aug 2013
2813 if (MathUtils::DeltaR(jet, phHard) < 0.5) continue;
2814 if (MathUtils::DeltaR(jet, phSoft) < 0.5) continue;
2815 // Apply Delta R(lepton, jet), see emails from Franceso Micheli
2816 // from 24 and 26 Nov 2013
2817 if (lepton && MathUtils::DeltaR(jet, lepton) < 0.5) continue;
2818 // Make sure we have a PF jet
2819 const PFJet *pfjet = dynamic_cast<const PFJet*>(jet);
2820 if (!pfjet) continue;
2821 if (!JetTools::passPFLooseId(pfjet)) continue;
2822 // Apply the jet ID as given in Table 4
2823 Double_t betaStar = JetTools::betaStarClassic(pfjet, selvtx, fPV);
2824 if (betaStar > 0.2 * log(fPV->GetEntries() - 0.64)) continue;
2825 if (JetTools::dR2Mean(pfjet, -1) > 0.065) continue;
2826 // this jet passes, count it in
2827 ++nJets;
2828 // Select b-jets that pass the CSV medium working point, see L128 of the AN
2829 // and https://twiki.cern.ch/twiki/bin/viewauth/CMS/BTagPerformanceOP
2830 if (jet->CombinedSecondaryVertexBJetTagsDisc() < 0.679) continue;
2831 ++nBJets;
2832 } // End of loop over jets
2833
2834 // Check the hadron tag, see Table 7 near L196 of the AN
2835 if (nJets >= 5 && nBJets >= 1) fDiphotonEvent->tthTag += 1;
2836
2837 // Check the lepton tag, see Table 7 near L196 of the AN
2838 // It has a precedence if both the leptonic and hadronic tags pass.
2839 // (private e-mail from Francesco Micheli on 15 July 2013).
2840 if (nJets >= 2 && nBJets >= 1 && lepton) fDiphotonEvent->tthTag += 2;
2841
2842 if (fVerbosityLevel > 90) {
2843 cout << "+ nJets: " << nJets << endl
2844 << "+ nBJets: " << nBJets << endl
2845 << "+ lepton: " << lepton << endl
2846 << "+ tthTag: " << fDiphotonEvent->tthTag << endl;
2847 }
2848
2849
2850 } // void PhotonTreeWriter::ApplyTTHTag()
2851
2852
2853 //_____________________________________________________________________________
2854 UInt_t PhotonTreeWriter::NumberOfJets(const Photon *phHard,
2855 const Photon *phSoft,
2856 const Vertex *selvtx,
2857 const double minJetPt,
2858 const double maxAbsEta) {
2859
2860 UInt_t nJets = 0;
2861
2862 // Loop over jets, count those passing selection
2863 // Use same ID as for the tth tag
2864 for(UInt_t ijet=0; ijet < fPFJets->GetEntries(); ++ijet){
2865 const Jet *jet = fPFJets->At(ijet);
2866 // Apply jet selection, see L116 and L125 of the AN
2867 if (jet->Pt() < minJetPt || jet->AbsEta() > maxAbsEta) continue;
2868 // Apply the cut Delta R(photon, jet) < 0.5.
2869 if (MathUtils::DeltaR(jet, phHard) < 0.5) continue;
2870 if (MathUtils::DeltaR(jet, phSoft) < 0.5) continue;
2871 // Make sure we have a PF jet
2872 const PFJet *pfjet = dynamic_cast<const PFJet*>(jet);
2873 if (!pfjet) continue;
2874 if (!JetTools::passPFLooseId(pfjet)) continue;
2875 // Apply the jet ID / pileup removal as given in Table 4
2876 Double_t betaStar = JetTools::betaStarClassic(pfjet, selvtx, fPV);
2877 if (betaStar > 0.2 * log(fPV->GetEntries() - 0.64)) continue;
2878 if (JetTools::dR2Mean(pfjet, -1) > 0.065) continue;
2879 // this jet passes, count it in
2880 ++nJets;
2881 } // End of loop over jets
2882
2883 return nJets;
2884
2885 } // NumberOfJets
2886
2887
2888 //_____________________________________________________________________________
2889 // PhotonTreeWriter::PFJetVector *
2890 // PhotonTreeWriter::GetSelectedPFJets(const DeltaRVetoVector &drVetos,
2891 // const Vertex &vertex,
2892 // const double minPt,
2893 // const double maxAbsEta) {
2894 //
2895 // PFJetVector *pfjets = new PFJetVector();
2896 //
2897 // // Loop over jets, count those passing selection
2898 // // Use same ID as for the tth tag
2899 // for(UInt_t ijet=0; ijet < fPFJets->GetEntries(); ++ijet){
2900 // const Jet *jet = fPFJets->At(ijet);
2901 // // Apply jet selection, see L116 and L125 of the AN
2902 // if (jet->Pt() < minPt || jet->AbsEta() > maxAbsEta) continue;
2903 // // Apply the vetos Delta R(jet, particle) > minDeltaR.
2904 // DeltaRVetoVector::const_iterator drVeto = drVetos.begin();
2905 // for (; drVeto < drVetos.end(); ++drVeto) {
2906 // const Particle *particle = drVeto->first;
2907 // std::cout << "run: " << fDiphotonEvent->run
2908 // << ", lumi: " << fDiphotonEvent->lumi
2909 // << ", events: " << fDiphotonEvent->evt
2910 // /// The line below causes a segmentation fault!
2911 // << ", veto particle pt: " << particle->Mom().Pt()
2912 // << ", eta: " << particle->Mom().Eta()
2913 // << ", phi: " << particle->Mom().Phi()
2914 // << "\n";
2915 // double minDeltaR = drVeto->second;
2916 // if (MathUtils::DeltaR(particle, jet) < minDeltaR) break;
2917 // } /// Loop over Delta R vetos.
2918 // if (drVeto < drVetos.end()) continue; /// failed a Delta R veto
2919 // // Make sure we have a PF jet
2920 // const PFJet *pfjet = dynamic_cast<const PFJet*>(jet);
2921 // if (!pfjet) continue;
2922 // if (!JetTools::passPFLooseId(pfjet)) continue;
2923 // // Apply the jet ID / pileup removal as given in Table 4
2924 // Double_t betaStar = JetTools::betaStarClassic(pfjet, &vertex, fPV);
2925 // if (betaStar > 0.2 * log(fPV->GetEntries() - 0.64)) continue;
2926 // if (JetTools::dR2Mean(pfjet, -1) > 0.065) continue;
2927 // // this jet passes, count it in
2928 // pfjets->push_back(pfjet);
2929 // } // End of loop over jets
2930 //
2931 // return pfjets;
2932 //
2933 // } // GetSelectedPFJets
2934
2935
2936 //_____________________________________________________________________________
2937 UInt_t PhotonTreeWriter::NumberOfBJets(const Photon *phHard,
2938 const Photon *phSoft,
2939 const Vertex *selvtx,
2940 const double minJetPt,
2941 const double maxAbsEta) {
2942
2943 UInt_t nBJets = 0;
2944
2945 // Loop over jets, count those passing selection
2946 // Use same ID as for the tth tag
2947 for(UInt_t ijet=0; ijet < fPFJets->GetEntries(); ++ijet){
2948 const Jet *jet = fPFJets->At(ijet);
2949 // Apply jet selection, see L116 and L125 of the AN
2950 if (jet->Pt() < minJetPt || jet->AbsEta() > maxAbsEta) continue;
2951 // Apply the cut Delta R(photon, jet) < 0.5.
2952 if (MathUtils::DeltaR(jet, phHard) < 0.5) continue;
2953 if (MathUtils::DeltaR(jet, phSoft) < 0.5) continue;
2954 // Make sure we have a PF jet
2955 const PFJet *pfjet = dynamic_cast<const PFJet*>(jet);
2956 if (!pfjet) continue;
2957 if (!JetTools::passPFLooseId(pfjet)) continue;
2958 // Apply the jet ID / pileup removal as given in Table 4
2959 Double_t betaStar = JetTools::betaStarClassic(pfjet, selvtx, fPV);
2960 if (betaStar > 0.2 * log(fPV->GetEntries() - 0.64)) continue;
2961 if (JetTools::dR2Mean(pfjet, -1) > 0.065) continue;
2962 // Select b-jets that pass the CSV medium working point, see L128 of the AN
2963 // and https://twiki.cern.ch/twiki/bin/viewauth/CMS/BTagPerformanceOP
2964 if (jet->CombinedSecondaryVertexBJetTagsDisc() < 0.679) continue;
2965 // this jet passes, count it in
2966 ++nBJets;
2967 if (fDoSynching && fDiphotonEvent->bjetcsv <= -99) {
2968 fDiphotonEvent->bjetcsv = jet->CombinedSecondaryVertexBJetTagsDisc();
2969 }
2970 } // End of loop over jets
2971
2972 return nBJets;
2973 } // NumberOfBJets
2974
2975
2976 //_____________________________________________________________________________
2977 const Particle *
2978 PhotonTreeWriter::TTHSelectLepton(const Photon *phHard,
2979 const Photon *phSoft,
2980 const Vertex *selvtx)
2981 {
2982 if (fVerbosityLevel > 90) {
2983 cout << "++ Entering TTHSelectLepton ..." << endl;
2984 }
2985 const Electron * electron = TTHSelectElectron(phHard, phSoft, selvtx);
2986 const Muon * muon = TTHSelectMuon (phHard, phSoft, selvtx);
2987 const Particle * lepton = 0;
2988 /// Select the lepton with the higher Pt (private discussion
2989 /// with Francesco Micheli on 26 Nov 2013)
2990 if (electron && muon) {
2991 if (electron->Pt() > muon->Pt()) {lepton = electron;}
2992 else {lepton = muon ;}
2993 } else if (electron && !muon) {
2994 lepton = electron;
2995 } else if (!electron && muon) {
2996 lepton = muon ;
2997 }
2998 if (fVerbosityLevel > 90) {
2999 cout << "++ electron: " << electron << endl
3000 << "++ muon: " << muon << endl
3001 << "++ lepton: " << lepton << endl
3002 << "++ Exiting TTHSelectLepton ..." << endl;
3003 }
3004 return lepton;
3005 } // TTHSelectLepton
3006
3007
3008 //_____________________________________________________________________________
3009 // Returns the highest-MVA-ID electron passing all cuts, 0 if none found.
3010 const Electron *
3011 PhotonTreeWriter::TTHSelectElectron(const Photon *phHard,
3012 const Photon *phSoft,
3013 const Vertex *selvtx)
3014 {
3015 if (fVerbosityLevel > 90) {
3016 cout << "+++ Entering TTHSelectElectron ..." << endl
3017 << "+++ PhotonTools::ElectronVetoCiC(phHard, fLeptonTagElectrons): "
3018 << PhotonTools::ElectronVetoCiC(phHard, fLeptonTagElectrons) << endl
3019 << "+++ PhotonTools::ElectronVetoCiC(phSoft, fLeptonTagElectrons): "
3020 << PhotonTools::ElectronVetoCiC(phSoft, fLeptonTagElectrons) << endl
3021 << "+++ PhotonTools::ElectronVetoCiC(phHard, fElectrons): "
3022 << PhotonTools::ElectronVetoCiC(phHard, fElectrons) << endl
3023 << "+++ PhotonTools::ElectronVetoCiC(phSoft, fElectrons): "
3024 << PhotonTools::ElectronVetoCiC(phSoft, fElectrons) << endl;
3025 }
3026
3027 if (PhotonTools::ElectronVetoCiC(phHard, fLeptonTagElectrons) >= 1 &&
3028 PhotonTools::ElectronVetoCiC(phSoft, fLeptonTagElectrons) >= 1 &&
3029 PhotonTools::ElectronVetoCiC(phHard, fElectrons) >= 1 &&
3030 PhotonTools::ElectronVetoCiC(phSoft, fElectrons) >= 1){
3031 // Loop over electrons, apply all cuts, find the one wiht hightes ID MVA
3032 for (UInt_t iele=0; iele < fLeptonTagElectrons->GetEntries(); ++iele) {
3033 const Electron *ele = fLeptonTagElectrons->At(iele);
3034
3035 if (fVerbosityLevel > 90) {
3036 cout << "+++ iele: " << iele << endl
3037 << "+++ ele: " << ele << endl
3038 << "+++ ele->Pt(): " << ele->Pt() << endl
3039 << "+++ ele->AbsEta(): " << ele->AbsEta() << endl
3040 << "+++ MathUtils::DeltaR(ele, phHard): "
3041 << MathUtils::DeltaR(ele, phHard) << endl
3042 << "+++ MathUtils::DeltaR(ele, phSoft): "
3043 << MathUtils::DeltaR(ele, phSoft) << endl
3044 << "+++ MassOfPairIsWithinWindowAroundMZ(ele, phHard, 10): "
3045 << MassOfPairIsWithinWindowAroundMZ(ele, phHard, 10) << endl
3046 << "+++ MassOfPairIsWithinWindowAroundMZ(ele, phSoft, 10): "
3047 << MassOfPairIsWithinWindowAroundMZ(ele, phSoft, 10) << endl;
3048 //<< "+++ GetElectronIdMva(ele): " << GetElectronIdMva(ele) << endl;
3049 }
3050
3051 // Apply kinematic cuts, see L133 and L134 of the AN
3052 //(JOSH: In fact these are already applied in the ElectronIDMod as we have configured it)
3053 //if (ele->Pt() < 20. || ele->AbsEta() > 2.5) continue;
3054 //keep Pt cut here in any case we reorganize logic later
3055 if (ele->Pt()<20.) continue;
3056
3057 // Require separation between this electron and both photons,
3058 // see the slide 7, bullet 2
3059 if (MathUtils::DeltaR(ele, phHard) < 0.5) continue;
3060 if (MathUtils::DeltaR(ele, phSoft) < 0.5) continue;
3061 // Require electron-photon mass outside of a 20 GeV window around MZ
3062 //JOSH: Shouldn't be applied for ttH category
3063 // if (MassOfPairIsWithinWindowAroundMZ(ele, phHard, 10)) continue;
3064 // if (MassOfPairIsWithinWindowAroundMZ(ele, phSoft, 10)) continue;
3065 // Electron ID MVA arbitration (private discussion with
3066 // Francesco Micheli on 26 Nov 2013)
3067 // double idMva = GetElectronIdMva(ele);
3068 // if (idMva > maxIdMva) {
3069 // maxIdMva = idMva;
3070 // selectedElectron = ele;
3071 // }
3072 //JOSH: electrons here are already sorted by mva output, so can just take the first one and break
3073 //(in fact the call to GetElectronIdMva was segfaulting for a reason which I didn't care to debug)
3074 return ele;
3075 } // Loop over electrons
3076 }
3077
3078 return 0;
3079
3080 } // TTHSelectElectron
3081
3082
3083 //_____________________________________________________________________________
3084 // Returns the highest-pt muon passing all the cuts, 0 if none found.
3085 const Muon *
3086 PhotonTreeWriter::TTHSelectMuon(const Photon *phHard,
3087 const Photon *phSoft,
3088 const Vertex *selvtx)
3089 {
3090 if (fVerbosityLevel > 90) {
3091 cout << "+++ Entering TTHSelectMuon ..." << endl;
3092 }
3093 const Muon *selectedMuon = 0;
3094 // Loop over muons
3095 for (UInt_t imu=0; imu < fLeptonTagMuons->GetEntries(); ++imu) {
3096 const Muon *mu = fLeptonTagMuons->At(imu);
3097 // Apply kinematic cuts, see L132 and L134 of the AN
3098 if (mu->Pt() < 20. || mu->AbsEta() > 2.4) continue;
3099 // Same as for electrons, require separation between this muon and both
3100 // photons.
3101 // Also confirmed by Francesco Micheli in an e-mail from 15 July 2013
3102 if (MathUtils::DeltaR(mu, phHard) < 0.5) continue;
3103 if (MathUtils::DeltaR(mu, phSoft) < 0.5) continue;
3104 // The first muon found is the highest-pt one because muons are pt-ordered.
3105 selectedMuon = mu;
3106 break;
3107 }
3108 if (fVerbosityLevel > 90) {
3109 cout << "+++ fLeptonTagMuons->GetEntries(): "
3110 << fLeptonTagMuons->GetEntries() << endl
3111 << "+++ Looping over muons ..." << endl;
3112 for (UInt_t imu=0; imu < fLeptonTagMuons->GetEntries(); ++imu) {
3113 const Muon *mu = fLeptonTagMuons->At(imu);
3114 cout << "++++ imu: " << imu << endl
3115 << "++++ mu: " << mu << endl
3116 << "++++ mu pt: " << mu->Pt() << endl
3117 << "++++ mu |eta|: " << mu->AbsEta() << endl
3118 << "++++ MathUtils::DeltaR(mu, phHard): "
3119 << MathUtils::DeltaR(mu, phHard) << endl
3120 << "++++ MathUtils::DeltaR(mu, phSoft): "
3121 << MathUtils::DeltaR(mu, phSoft) << endl;
3122 }
3123 cout << "+++ Finished loop over muons ..." << endl
3124 << "+++ selectedMuon: " << selectedMuon << endl
3125 << "+++ Exiting TTHSelectMuon ..." << endl << flush;
3126 }
3127 return selectedMuon;
3128 } // TTHSelectMuon
3129
3130
3131 //_____________________________________________________________________________
3132 template<typename Element, typename Collection>
3133 UInt_t PhotonTreeWriter::IndexOfElementInCollection(
3134 const Element * element,
3135 const Collection * collection
3136 )
3137 {
3138 UInt_t index = 0;
3139 for (; index < collection->GetEntries() && index < (UInt_t) -1; ++index) {
3140 if (element == collection->At(index)) {
3141 break;
3142 }
3143 }
3144 return index;
3145 }
3146
3147 //_____________________________________________________________________________
3148 UInt_t PhotonTreeWriter::IndexOfNearestSuperClusterInCollection(
3149 const SuperCluster *element,
3150 const SuperClusterCol *collection
3151 )
3152 {
3153 double minMass = 999.;
3154 UInt_t minIndex = 0;
3155 if (collection->GetEntries() > 0) {
3156 minMass = SuperClusterPairMass(element, collection->At(0));
3157 }
3158 for (UInt_t index = 1;
3159 index < collection->GetEntries() && index < (UInt_t) -1; ++index) {
3160 double mass = SuperClusterPairMass(element, collection->At(index));
3161 if (mass < minMass) {
3162 minMass = mass;
3163 minIndex = index;
3164 }
3165 }
3166 return minIndex;
3167 }
3168
3169
3170 //_____________________________________________________________________________
3171 double PhotonTreeWriter::SuperClusterPairMass(const SuperCluster *sc1,
3172 const SuperCluster *sc2)
3173 {
3174 FourVectorM p1 = SuperClusterFourVectorM(sc1);
3175 FourVectorM p2 = SuperClusterFourVectorM(sc2);
3176 return (p1 + p2).M();
3177 }
3178
3179
3180 //_____________________________________________________________________________
3181 FourVectorM PhotonTreeWriter::SuperClusterFourVectorM(const SuperCluster *sc)
3182 {
3183 double e = sc->Energy();
3184 double eta = sc->Eta();
3185 double phi = sc->Phi();
3186 double pt = e / TMath::CosH(eta);
3187 return FourVectorM(pt, eta, phi, 0.);
3188 }
3189
3190
3191 //_____________________________________________________________________________
3192 void PhotonTreeWriter::Terminate()
3193 {
3194 // Run finishing code on the computer (slave) that did the analysis
3195 }
3196
3197
3198 //_____________________________________________________________________________
3199 void PhotonTreeWriter::LogEventInfo()
3200 {
3201 time_t rawtime;
3202 struct tm * timeinfo;
3203 char buffer[80];
3204
3205 time(&rawtime);
3206 timeinfo = localtime(&rawtime);
3207 strftime(buffer, 80, "%c", timeinfo);
3208
3209 const EventHeader *event = GetEventHeader();
3210 cout << buffer << ": Processing " << fProcessedEvents << ". record"
3211 << ", run " << event->RunNum()
3212 << ", lumi " << event->LumiSec()
3213 << ", event " << event->EvtNum() << endl;
3214 }
3215
3216
3217 //______________________________________________________________________________
3218 double
3219 PhotonTreeWriter::GetElectronIdMva(const Electron *electron)
3220 {
3221 Int_t closestVtx = 0;
3222 Double_t distVtx = 999.0;
3223 /// Find the closest vertex in dz.
3224 for(UInt_t nv=0; nv<fPV->GetEntries(); nv++){
3225 double dz = TMath::Abs(electron->GsfTrk()->DzCorrected(*fPV->At(nv)));
3226 if(dz < distVtx) {
3227 distVtx = dz;
3228 closestVtx = nv;
3229 }
3230 }
3231 return fElectronIDMVA->MVAValue(electron, fPV->At(closestVtx));
3232 } // GetElectronIdMva