ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.79
Committed: Mon Dec 16 16:52:36 2013 UTC (11 years, 4 months ago) by mingyang
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.78: +21 -5 lines
Log Message: 
r9 rescale added

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