ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.55
Committed: Mon Oct 7 12:47:49 2013 UTC (11 years, 7 months ago) by veverka
Content type: text/plain
Branch: MAIN
Changes since 1.54: +65 -18 lines
Log Message: 
Added photon indexes and renamed SuperClusters to PFSuperClusters

File Contents

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