ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.13
Committed: Tue Apr 24 11:45:53 2012 UTC (13 years ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.12: +58 -1 lines
Log Message:
added features for Hgg LeptonTag analysis

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 "TDataMember.h"
12 #include <TNtuple.h>
13 #include <TRandom3.h>
14 #include <TSystem.h>
15
16 using namespace mithep;
17
18 ClassImp(mithep::PhotonTreeWriter)
19 ClassImp(mithep::PhotonTreeWriterPhoton)
20 ClassImp(mithep::PhotonTreeWriterDiphotonEvent)
21
22 //--------------------------------------------------------------------------------------------------
23 PhotonTreeWriter::PhotonTreeWriter(const char *name, const char *title) :
24 // Base Module...
25 BaseMod (name,title),
26 // define all the Branches to load
27 fPhotonBranchName (Names::gkPhotonBrn),
28 fElectronName (Names::gkElectronBrn),
29 fGoodElectronName (Names::gkElectronBrn),
30 fConversionName (Names::gkMvfConversionBrn),
31 fTrackBranchName (Names::gkTrackBrn),
32 fPileUpDenName (Names::gkPileupEnergyDensityBrn),
33 fPVName (Names::gkPVBeamSpotBrn),
34 fBeamspotName (Names::gkBeamSpotBrn),
35 fPFCandName (Names::gkPFCandidatesBrn),
36 fMCParticleName (Names::gkMCPartBrn),
37 fPileUpName (Names::gkPileupInfoBrn),
38 fSuperClusterName ("PFSuperClusters"),
39 fPFMetName ("PFMet"),
40 fPFJetName (Names::gkPFJetBrn),
41
42 fLeptonTagElectronsName ("HggLeptonTagElectrons"),
43 fLeptonTagMuonsName ("HggLeptonTagMuons"),
44
45 fIsData (false),
46 fPhotonsFromBranch (kTRUE),
47 fPVFromBranch (kTRUE),
48 fGoodElectronsFromBranch(kTRUE),
49 fPFJetsFromBranch (kTRUE),
50 // ----------------------------------------
51 // collections....
52 fPhotons (0),
53 fElectrons (0),
54 fConversions (0),
55 fTracks (0),
56 fPileUpDen (0),
57 fPV (0),
58 fBeamspot (0),
59 fPFCands (0),
60 fMCParticles (0),
61 fPileUp (0),
62 fSuperClusters (0),
63 fPFJets (0),
64
65 fLeptonTagElectrons (0),
66 fLeptonTagMuons (0),
67
68 fLoopOnGoodElectrons (kFALSE),
69 fApplyElectronVeto (kTRUE),
70 fWriteDiphotonTree (kTRUE),
71 fWriteSingleTree (kTRUE),
72 fExcludeSinglePrompt (kFALSE),
73 fExcludeDoublePrompt (kFALSE),
74 fEnableJets (kFALSE),
75 fApplyLeptonTag (kFALSE),
76 fPhFixDataFile (gSystem->Getenv("CMSSW_BASE") +
77 TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat")),
78 fTupleName ("hPhotonTree")
79 {
80 // Constructor
81 }
82
83 PhotonTreeWriter::~PhotonTreeWriter()
84 {
85 // Destructor
86 }
87
88 //--------------------------------------------------------------------------------------------------
89 void PhotonTreeWriter::Process()
90 {
91 // ------------------------------------------------------------
92 // Process entries of the tree.
93 LoadEventObject(fPhotonBranchName, fPhotons);
94 LoadEventObject(fGoodElectronName, fGoodElectrons);
95
96 // lepton tag collections
97 if( fApplyLeptonTag ) {
98 LoadEventObject(fLeptonTagElectronsName, fLeptonTagElectrons);
99 LoadEventObject(fLeptonTagMuonsName, fLeptonTagMuons);
100 }
101
102 const BaseCollection *egcol = 0;
103 if (fLoopOnGoodElectrons)
104 egcol = fGoodElectrons;
105 else
106 egcol = fPhotons;
107 if (egcol->GetEntries()<1)
108 return;
109 LoadEventObject(fElectronName, fElectrons);
110 LoadEventObject(fConversionName, fConversions);
111 LoadEventObject(fTrackBranchName, fTracks);
112 LoadEventObject(fPileUpDenName, fPileUpDen);
113 LoadEventObject(fPVName, fPV);
114 LoadEventObject(fBeamspotName, fBeamspot);
115 LoadEventObject(fPFCandName, fPFCands);
116 LoadEventObject(fSuperClusterName, fSuperClusters);
117 LoadEventObject(fPFMetName, fPFMet);
118 if (fEnableJets)
119 LoadEventObject(fPFJetName, fPFJets);
120
121 // ------------------------------------------------------------
122 // load event based information
123 Int_t _numPU = -1.; // some sensible default values....
124 Int_t _numPUminus = -1.; // some sensible default values....
125 Int_t _numPUplus = -1.; // some sensible default values....
126
127 Float_t rho = -99.;
128 if( fPileUpDen->GetEntries() > 0 )
129 rho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
130
131 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
132
133 if( !fIsData ) {
134 LoadBranch(fMCParticleName);
135 LoadBranch(fPileUpName);
136 }
137
138 if( !fIsData ) {
139 for (UInt_t i=0; i<fPileUp->GetEntries(); ++i) {
140 const PileupInfo *puinfo = fPileUp->At(i);
141 if (puinfo->GetBunchCrossing()==0) _numPU = puinfo->GetPU_NumInteractions();
142 else if (puinfo->GetBunchCrossing() == -1) _numPUminus = puinfo->GetPU_NumInteractions();
143 else if (puinfo->GetBunchCrossing() == 1) _numPUplus = puinfo->GetPU_NumInteractions();
144 }
145 }
146
147 // in case of a MC event, try to find Higgs and Higgs decay Z poisition
148 Float_t _pth = -100.;
149 Float_t _decayZ = -100.;
150 Float_t _genmass = -100.;
151 if (!fIsData)
152 FindHiggsPtAndZ(_pth, _decayZ, _genmass);
153
154 fDiphotonEvent->leptonTag = -1; // disabled
155
156 fDiphotonEvent->rho = rho;
157 fDiphotonEvent->genHiggspt = _pth;
158 fDiphotonEvent->genHiggsZ = _decayZ;
159 fDiphotonEvent->genmass = _genmass;
160 fDiphotonEvent->gencostheta = -99.;
161 fDiphotonEvent->nVtx = fPV->GetEntries();
162 fDiphotonEvent->bsX = fBeamspot->At(0)->X();
163 fDiphotonEvent->bsY = fBeamspot->At(0)->Y();
164 fDiphotonEvent->bsZ = fBeamspot->At(0)->Z();
165 fDiphotonEvent->bsSigmaZ = fBeamspot->At(0)->SigmaZ();
166 fDiphotonEvent->vtxX = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->X() : -99.;
167 fDiphotonEvent->vtxY = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Y() : -99.;
168 fDiphotonEvent->vtxZ = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Z() : -99.;
169 fDiphotonEvent->numPU = _numPU;
170 fDiphotonEvent->numPUminus = _numPUminus;
171 fDiphotonEvent->numPUplus = _numPUplus;
172 fDiphotonEvent->mass = -99.;
173 fDiphotonEvent->ptgg = -99.;
174 fDiphotonEvent->costheta = -99.;
175 fDiphotonEvent->mt = -99.;
176 fDiphotonEvent->cosphimet = -99.;
177 fDiphotonEvent->mtele = -99.;
178 fDiphotonEvent->cosphimetele = -99.;
179 fDiphotonEvent->evt = GetEventHeader()->EvtNum();
180 fDiphotonEvent->run = GetEventHeader()->RunNum();
181 fDiphotonEvent->lumi = GetEventHeader()->LumiSec();
182 fDiphotonEvent->evtcat = -99;
183 fDiphotonEvent->nobj = fPhotons->GetEntries();
184 fDiphotonEvent->pfmet = fPFMet->At(0)->Pt();
185 fDiphotonEvent->pfmetphi = fPFMet->At(0)->Phi();
186 fDiphotonEvent->pfmetx = fPFMet->At(0)->Px();
187 fDiphotonEvent->pfmety = fPFMet->At(0)->Py();
188 fDiphotonEvent->masscor = -99.;
189 fDiphotonEvent->masscorerr = -99.;
190 fDiphotonEvent->masscorele = -99.;
191 fDiphotonEvent->masscoreleerr = -99.;
192 fDiphotonEvent->ismc = GetEventHeader()->IsMC();
193
194 //jets
195 const Jet *jet1 = 0;
196 const Jet *jet2 = 0;
197 const Jet *jetcentral = 0;
198
199 fDiphotonEvent->jet1pt = -99.;
200 fDiphotonEvent->jet1eta = -99.;
201 fDiphotonEvent->jet1phi = -99.;
202 fDiphotonEvent->jet1mass = -99.;
203 fDiphotonEvent->jet2pt = -99.;
204 fDiphotonEvent->jet2eta = -99.;
205 fDiphotonEvent->jet2phi = -99.;
206 fDiphotonEvent->jet2mass = -99.;
207 fDiphotonEvent->jetcentralpt = -99.;
208 fDiphotonEvent->jetcentraleta = -99.;
209 fDiphotonEvent->jetcentralphi = -99.;
210 fDiphotonEvent->jetcentralmass = -99.;
211 fDiphotonEvent->dijetpt = -99.;
212 fDiphotonEvent->dijeteta = -99.;
213 fDiphotonEvent->dijetphi = -99.;
214 fDiphotonEvent->dijetmass = -99.;
215 fDiphotonEvent->jetetaplus = -99.;
216 fDiphotonEvent->jetetaminus = -99.;
217 fDiphotonEvent->zeppenfeld = -99.;
218 fDiphotonEvent->dphidijetgg = -99.;
219
220 Int_t nhitsbeforevtxmax = 1;
221 if (!fApplyElectronVeto)
222 nhitsbeforevtxmax = 999;
223
224 if (egcol->GetEntries()>=2) {
225 const Particle *p1 = 0;
226 const Particle *p2 = 0;
227 const Photon *phHard = 0;
228 const Photon *phSoft = 0;
229 const Electron *ele1 = 0;
230 const Electron *ele2 = 0;
231 const SuperCluster *sc1 = 0;
232 const SuperCluster *sc2 = 0;
233
234 if (fLoopOnGoodElectrons) {
235 ele1 = fGoodElectrons->At(0);
236 ele2 = fGoodElectrons->At(1);
237 p1 = ele1;
238 p2 = ele2;
239 sc1 = ele1->SCluster();
240 sc2 = ele2->SCluster();
241 phHard = PhotonTools::MatchedPhoton(ele1,fPhotons);
242 phSoft = PhotonTools::MatchedPhoton(ele2,fPhotons);
243 }
244 else {
245 phHard = fPhotons->At(0);
246 phSoft = fPhotons->At(1);
247 p1 = phHard;
248 p2 = phSoft;
249 sc1 = phHard->SCluster();
250 sc2 = phSoft->SCluster();
251 ele1 = PhotonTools::MatchedElectron(phHard,fGoodElectrons);
252 ele2 = PhotonTools::MatchedElectron(phSoft,fGoodElectrons);
253 }
254
255 const DecayParticle *conv1 = PhotonTools::MatchedConversion(sc1,fConversions,bsp,
256 nhitsbeforevtxmax);
257 const DecayParticle *conv2 = PhotonTools::MatchedConversion(sc2,fConversions,bsp,
258 nhitsbeforevtxmax);
259
260 const SuperCluster *pfsc1 = PhotonTools::MatchedSC(sc1,fSuperClusters);
261 const SuperCluster *pfsc2 = PhotonTools::MatchedSC(sc2,fSuperClusters);
262
263 const MCParticle *phgen1 = 0;
264 const MCParticle *phgen2 = 0;
265 if (!fIsData) {
266 phgen1 = PhotonTools::MatchMC(p1,fMCParticles,!fApplyElectronVeto);
267 phgen2 = PhotonTools::MatchMC(p2,fMCParticles,!fApplyElectronVeto);
268 }
269
270 if (fExcludeSinglePrompt && (phgen1 || phgen2))
271 return;
272 if (fExcludeDoublePrompt && (phgen1 && phgen2))
273 return;
274
275 if (!fLoopOnGoodElectrons && phHard->HasPV()) {
276 fDiphotonEvent->vtxX = phHard->PV()->X();
277 fDiphotonEvent->vtxY = phHard->PV()->Y();
278 fDiphotonEvent->vtxZ = phHard->PV()->Z();
279 fDiphotonEvent->vtxprob = phHard->VtxProb();
280 }
281
282 //fill jet variables
283 if (fEnableJets) {
284 for (UInt_t ijet=0; ijet<fPFJets->GetEntries();++ijet) {
285 const Jet *jet = fPFJets->At(ijet);
286 if (jet->AbsEta()<4.7 && MathUtils::DeltaR(jet,p1)>0.5 && MathUtils::DeltaR(jet,p2)>0.5) {
287 if (!jet1) jet1 = jet;
288 else if (!jet2) jet2 = jet;
289 else if (!jetcentral && 0) jetcentral = jet;
290 }
291 if (jet1&&jet2&&jetcentral) break;
292 }
293 }
294
295 if (jet1) {
296 fDiphotonEvent->jet1pt = jet1->Pt();
297 fDiphotonEvent->jet1eta = jet1->Eta();
298 fDiphotonEvent->jet1phi = jet1->Phi();
299 fDiphotonEvent->jet1mass = jet1->Mass();
300 }
301
302 if (jet2) {
303 fDiphotonEvent->jet2pt = jet2->Pt();
304 fDiphotonEvent->jet2eta = jet2->Eta();
305 fDiphotonEvent->jet2phi = jet2->Phi();
306 fDiphotonEvent->jet2mass = jet2->Mass();
307 }
308
309 if (jetcentral) {
310 fDiphotonEvent->jetcentralpt = jetcentral->Pt();
311 fDiphotonEvent->jetcentraleta = jetcentral->Eta();
312 fDiphotonEvent->jetcentralphi = jetcentral->Phi();
313 fDiphotonEvent->jetcentralmass = jetcentral->Mass();
314 }
315
316 if (jet1&&jet2){
317 FourVectorM momjj = (jet1->Mom() + jet2->Mom());
318
319 fDiphotonEvent->dijetpt = momjj.Pt();
320 fDiphotonEvent->dijeteta = momjj.Eta();
321 fDiphotonEvent->dijetphi = momjj.Phi();
322 fDiphotonEvent->dijetmass = momjj.M();
323
324 if (jet1->Eta()>jet2->Eta()) {
325 fDiphotonEvent->jetetaplus = jet1->Eta();
326 fDiphotonEvent->jetetaminus = jet2->Eta();
327 }
328 else {
329 fDiphotonEvent->jetetaplus = jet2->Eta();
330 fDiphotonEvent->jetetaminus = jet1->Eta();
331 }
332 }
333
334 Double_t _mass = -99.;
335 Double_t _masserr = -99.;
336 Double_t _masserrsmeared = -99.;
337 Double_t _masserrwrongvtx = -99.;
338 Double_t _masserrsmearedwrongvtx = -99.;
339 Double_t _ptgg = -99.;
340 Double_t _etagg = -99.;
341 Double_t _phigg = -99.;
342 Double_t _costheta = -99.;
343 PhotonTools::DiphotonR9EtaPtCats _evtcat = PhotonTools::kOctCat0;
344 if (phHard && phSoft) {
345 _mass = (phHard->Mom()+phSoft->Mom()).M();
346 _masserr = 0.5*_mass*TMath::Sqrt(phHard->EnergyErr()*phHard->EnergyErr()/phHard->E()/phHard->E() + phSoft->EnergyErr()*phSoft->EnergyErr()/phSoft->E()/phSoft->E());
347 _masserrsmeared = 0.5*_mass*TMath::Sqrt(phHard->EnergyErrSmeared()*phHard->EnergyErrSmeared()/phHard->E()/phHard->E() + phSoft->EnergyErrSmeared()*phSoft->EnergyErrSmeared()/phSoft->E()/phSoft->E());
348 _ptgg = (phHard->Mom()+phSoft->Mom()).Pt();
349 _etagg = (phHard->Mom()+phSoft->Mom()).Eta();
350 _phigg = (phHard->Mom()+phSoft->Mom()).Phi();
351 _costheta = ThreeVector(phHard->Mom()).Unit().Dot(ThreeVector(phSoft->Mom()).Unit());
352 _evtcat = PhotonTools::DiphotonR9EtaPtCat(phHard,phSoft);
353
354 const Double_t dz = sqrt(2.0)*5.8;
355 Double_t deltamvtx = _mass*VertexTools::DeltaMassVtx(phHard->CaloPos().X(),
356 phHard->CaloPos().Y(),
357 phHard->CaloPos().Z(),
358 phSoft->CaloPos().X(),
359 phSoft->CaloPos().Y(),
360 phSoft->CaloPos().Z(),
361 fDiphotonEvent->vtxX,
362 fDiphotonEvent->vtxY,
363 fDiphotonEvent->vtxZ,
364 dz);
365 fDiphotonEvent->deltamvtx = deltamvtx;
366
367 _masserrwrongvtx = TMath::Sqrt(_masserr*_masserr + deltamvtx*deltamvtx);
368 _masserrsmearedwrongvtx = TMath::Sqrt(_masserrsmeared*_masserrsmeared + deltamvtx*deltamvtx);
369
370 if (jet1 && jet2) {
371 fDiphotonEvent->zeppenfeld = TMath::Abs(_etagg - 0.5*(jet1->Eta()+jet2->Eta()));
372 fDiphotonEvent->dphidijetgg = MathUtils::DeltaPhi( (jet1->Mom()+jet2->Mom()).Phi(), _phigg );
373 }
374
375 }
376
377
378 Float_t _massele = -99.;
379 Float_t _ptee = -99.;
380 Float_t _costhetaele = -99.;
381 if (ele1 && ele2) {
382 _massele = (ele1->Mom()+ele2->Mom()).M();
383 _ptee = (ele1->Mom()+ele2->Mom()).Pt();
384 _costhetaele = ThreeVector(ele1->Mom()).Unit().Dot(ThreeVector(ele2->Mom()).Unit());
385 }
386
387 Float_t _gencostheta = -99.;
388 if (phgen1 && phgen2) {
389 _gencostheta = ThreeVector(phgen1->Mom()).Unit().Dot(ThreeVector(phgen2->Mom()).Unit());
390 }
391
392 fDiphotonEvent->gencostheta = _gencostheta;
393 fDiphotonEvent->mass = _mass;
394 fDiphotonEvent->masserr = _masserr;
395 fDiphotonEvent->masserrsmeared = _masserrsmeared;
396 fDiphotonEvent->masserrwrongvtx = _masserrwrongvtx;
397 fDiphotonEvent->masserrsmearedwrongvtx = _masserrsmearedwrongvtx;
398 fDiphotonEvent->ptgg = _ptgg;
399 fDiphotonEvent->etagg = _etagg;
400 fDiphotonEvent->phigg = _phigg;
401 fDiphotonEvent->costheta = _costheta;;
402 fDiphotonEvent->massele = _massele;
403 fDiphotonEvent->ptee = _ptee;
404 fDiphotonEvent->costhetaele = _costhetaele;
405 fDiphotonEvent->evtcat = _evtcat;
406
407 fDiphotonEvent->photons[0].SetVars(phHard,conv1,ele1,pfsc1,phgen1,fPhfixph,fPhfixele,fTracks,fPV,rho,fElectrons,fApplyElectronVeto);
408 fDiphotonEvent->photons[1].SetVars(phSoft,conv2,ele2,pfsc2,phgen2,fPhfixph,fPhfixele,fTracks,fPV,rho,fElectrons,fApplyElectronVeto);
409
410 Float_t ph1ecor = fDiphotonEvent->photons[0].Ecor();
411 Float_t ph1ecorerr = fDiphotonEvent->photons[0].Ecorerr();
412 Float_t ph2ecor = fDiphotonEvent->photons[1].Ecor();
413 Float_t ph2ecorerr = fDiphotonEvent->photons[1].Ecorerr();
414
415 Float_t ph1ecorele = fDiphotonEvent->photons[0].Ecorele();
416 Float_t ph1ecoreleerr = fDiphotonEvent->photons[0].Ecoreleerr();
417 Float_t ph2ecorele = fDiphotonEvent->photons[1].Ecorele();
418 Float_t ph2ecoreleerr = fDiphotonEvent->photons[1].Ecoreleerr();
419
420
421 fDiphotonEvent->masscor = TMath::Sqrt(2.0*ph1ecor*ph2ecor*(1.0-fDiphotonEvent->costheta));
422 fDiphotonEvent->masscorerr = 0.5*fDiphotonEvent->masscor*
423 TMath::Sqrt(ph1ecorerr*ph1ecorerr/ph1ecor/ph1ecor + ph2ecorerr*ph2ecorerr/ph2ecor/ph2ecor);
424
425 fDiphotonEvent->masscorele = TMath::Sqrt(2.0*ph1ecorele*ph2ecorele*
426 (1.0-fDiphotonEvent->costheta));
427 fDiphotonEvent->masscoreleerr = 0.5*fDiphotonEvent->masscorele*
428 TMath::Sqrt(ph1ecoreleerr*ph1ecoreleerr/ph1ecorele/ph1ecorele +
429 ph2ecoreleerr*ph2ecoreleerr/ph2ecorele/ph2ecorele);
430
431 //printf("r9 = %5f, photon sigieie = %5f, seed sigieie = %5f\n",phHard->R9(),
432 // phHard->CoviEtaiEta(),sqrt(phHard->SCluster()->Seed()->CoviEtaiEta()));
433
434 if( fApplyLeptonTag ) {
435 // perform lepton tagging
436 // the diphoton event record will have one more entry; i.e. leptonTag
437 // leptonTag = -1 -> lepton-taggng was swicthed off
438 // = 0 -> event tagged as 'non-lepton-event'
439 // = +1 -> event tagged as muon-event
440 // = +2 -> event tagged as electron-event
441 fDiphotonEvent->leptonTag = 0;
442
443 if ( fLeptonTagMuons->GetEntries() > 0 ) {
444 // need to have dR > 1 for with respect to both photons
445 for(UInt_t iMuon = 0; iMuon <fLeptonTagMuons->GetEntries(); ++iMuon) {
446 if(MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phHard) < 1.) continue;
447 if(MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phSoft) < 1.) continue;
448
449 fDiphotonEvent->leptonTag = 1;
450 break;
451 }
452 }
453 if ( fDiphotonEvent->leptonTag < 1 && fLeptonTagElectrons->GetEntries() > 0 ) {
454 for(UInt_t iEle = 0; iEle < fLeptonTagElectrons->GetEntries(); ++iEle) {
455 if(MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phHard) < 1.) continue;
456 if(MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phSoft) < 1.) continue;
457
458 // here we also check the mass ....
459 if ( TMath::Abs( (phHard->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M()-91.19 ) < 5. ) continue;
460 if ( TMath::Abs( (phSoft->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M()-91.19 ) < 5. ) continue;
461
462 fDiphotonEvent->leptonTag = 2;
463 break;
464 }
465 }
466 }
467
468 if (fWriteDiphotonTree)
469 hCiCTuple->Fill();
470 }
471
472 if (!fWriteSingleTree)
473 return;
474
475 for (UInt_t iph = 0; iph<egcol->GetEntries(); ++iph) {
476 const Particle *p = 0;
477 const Photon *ph = 0;
478 const Electron *ele = 0;
479 const SuperCluster *sc = 0;
480 if (fLoopOnGoodElectrons) {
481 ele = fGoodElectrons->At(iph);
482 p = ele;
483 sc = ele->SCluster();
484 ph = PhotonTools::MatchedPhoton(ele,fPhotons);
485 }
486 else {
487 ph = fPhotons->At(iph);
488 p = ph;
489 sc = ph->SCluster();
490 ele = PhotonTools::MatchedElectron(ph,fGoodElectrons);
491 }
492
493 const DecayParticle *conv = PhotonTools::MatchedConversion(sc,fConversions,bsp,
494 nhitsbeforevtxmax);
495 const SuperCluster *pfsc = PhotonTools::MatchedSC(sc,fSuperClusters);
496
497 if (!fLoopOnGoodElectrons && ph->HasPV()) {
498 fDiphotonEvent->vtxZ = ph->PV()->Z();
499 }
500
501 const MCParticle *phgen = 0;
502 if( !fIsData ) {
503 phgen = PhotonTools::MatchMC(p,fMCParticles,!fApplyElectronVeto);
504 }
505
506 if (fExcludeSinglePrompt && phgen) return;
507
508 fDiphotonEvent->mt = -99.;
509 fDiphotonEvent->cosphimet = -99.;
510 fDiphotonEvent->mtele = -99.;
511 fDiphotonEvent->cosphimetele = -99.;
512
513 if (ph) {
514 fDiphotonEvent->cosphimet = TMath::Cos(ph->Phi()-fPFMet->At(0)->Phi());
515 fDiphotonEvent->mt = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ph->Pt()*
516 (1.0-fDiphotonEvent->cosphimet));
517 }
518
519 if (ele) {
520 fDiphotonEvent->cosphimetele = TMath::Cos(ele->Phi()-fPFMet->At(0)->Phi());
521 fDiphotonEvent->mtele = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ele->Pt()*
522 (1.0-fDiphotonEvent->cosphimetele));
523 }
524
525 fSinglePhoton->SetVars(ph,conv,ele,pfsc,phgen,fPhfixph,fPhfixele,fTracks,fPV,rho,
526 fElectrons,fApplyElectronVeto);
527 hCiCTupleSingle->Fill();
528 }
529
530 return;
531 }
532
533 //--------------------------------------------------------------------------------------------------
534 void PhotonTreeWriter::SlaveBegin()
535 {
536 // Run startup code on the computer (slave) doing the actual analysis. Here,
537 // we just request the photon collection branch.
538
539 if( fApplyLeptonTag ) {
540 ReqEventObject(fLeptonTagElectronsName, fLeptonTagElectrons, false);
541 ReqEventObject(fLeptonTagMuonsName, fLeptonTagMuons, false);
542 }
543
544 ReqEventObject(fPhotonBranchName,fPhotons, fPhotonsFromBranch);
545 ReqEventObject(fTrackBranchName, fTracks, true);
546 ReqEventObject(fElectronName, fElectrons, true);
547 ReqEventObject(fGoodElectronName,fGoodElectrons,fGoodElectronsFromBranch);
548 ReqEventObject(fPileUpDenName, fPileUpDen, true);
549 ReqEventObject(fPVName, fPV, fPVFromBranch);
550 ReqEventObject(fConversionName, fConversions, true);
551 ReqEventObject(fBeamspotName, fBeamspot, true);
552 ReqEventObject(fPFCandName, fPFCands, true);
553 ReqEventObject(fSuperClusterName,fSuperClusters,true);
554 ReqEventObject(fPFMetName, fPFMet, true);
555 if (fEnableJets)
556 ReqEventObject(fPFJetName, fPFJets, fPFJetsFromBranch);
557 if (!fIsData) {
558 ReqBranch(fPileUpName, fPileUp);
559 ReqBranch(fMCParticleName, fMCParticles);
560 }
561 if (fIsData) {
562 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
563 TString("/src/MitPhysics/data/PhotonFixGRPV22.dat");
564 }
565 else {
566 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
567 TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat");
568 }
569
570 fPhfixph.initialise("4_2",std::string(fPhFixDataFile));
571 fPhfixele.initialise("4_2e",std::string(fPhFixDataFile));
572
573 fDiphotonEvent = new PhotonTreeWriterDiphotonEvent;
574 fSinglePhoton = new PhotonTreeWriterPhoton;
575
576 if (fWriteDiphotonTree)
577 hCiCTuple = new TTree(fTupleName.Data(),fTupleName.Data());
578 TString singlename = fTupleName + TString("Single");
579 if (fWriteSingleTree)
580 hCiCTupleSingle = new TTree(singlename,singlename);
581
582 //make flattish tree from classes so we don't have to rely on dictionaries for reading later
583 TClass *eclass = TClass::GetClass("mithep::PhotonTreeWriterDiphotonEvent");
584 TClass *pclass = TClass::GetClass("mithep::PhotonTreeWriterPhoton");
585 TList *elist = eclass->GetListOfDataMembers();
586 TList *plist = pclass->GetListOfDataMembers();
587
588 for (int i=0; i<elist->GetEntries(); ++i) {
589 const TDataMember *tdm = static_cast<const TDataMember*>(elist->At(i));
590 if (!(tdm->IsBasic() && tdm->IsPersistent()))
591 continue;
592 TString typestring;
593 if (TString(tdm->GetTypeName())=="Char_t") typestring = "B";
594 else if (TString(tdm->GetTypeName())=="UChar_t") typestring = "b";
595 else if (TString(tdm->GetTypeName())=="Short_t") typestring = "S";
596 else if (TString(tdm->GetTypeName())=="UShort_t") typestring = "s";
597 else if (TString(tdm->GetTypeName())=="Int_t") typestring = "I";
598 else if (TString(tdm->GetTypeName())=="UInt_t") typestring = "i";
599 else if (TString(tdm->GetTypeName())=="Float_t") typestring = "F";
600 else if (TString(tdm->GetTypeName())=="Double_t") typestring = "D";
601 else if (TString(tdm->GetTypeName())=="Long64_t") typestring = "L";
602 else if (TString(tdm->GetTypeName())=="ULong64_t") typestring = "l";
603 else if (TString(tdm->GetTypeName())=="Bool_t") typestring = "O";
604 else
605 continue;
606 //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
607 Char_t *addr = (Char_t*)fDiphotonEvent;
608 assert(sizeof(Char_t)==1);
609
610 if (fWriteDiphotonTree)
611 hCiCTuple->Branch(tdm->GetName(),addr + tdm->GetOffset(),
612 TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
613 if (fWriteSingleTree)
614 hCiCTupleSingle->Branch(tdm->GetName(),addr + tdm->GetOffset(),
615 TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
616 }
617
618 for (int iph=0; iph<2; ++iph) {
619 for (int i=0; i<plist->GetEntries(); ++i) {
620 const TDataMember *tdm = static_cast<const TDataMember*>(plist->At(i));
621 if (!(tdm->IsBasic() && tdm->IsPersistent()))
622 continue;
623 TString typestring;
624 if (TString(tdm->GetTypeName())=="Char_t")
625 typestring = "B";
626 else if (TString(tdm->GetTypeName())=="UChar_t")
627 typestring = "b";
628 else if (TString(tdm->GetTypeName())=="Short_t")
629 typestring = "S";
630 else if (TString(tdm->GetTypeName())=="UShort_t")
631 typestring = "s";
632 else if (TString(tdm->GetTypeName())=="Int_t")
633 typestring = "I";
634 else if (TString(tdm->GetTypeName())=="UInt_t")
635 typestring = "i";
636 else if (TString(tdm->GetTypeName())=="Float_t")
637 typestring = "F";
638 else if (TString(tdm->GetTypeName())=="Double_t")
639 typestring = "D";
640 else if (TString(tdm->GetTypeName())=="Long64_t")
641 typestring = "L";
642 else if (TString(tdm->GetTypeName())=="ULong64_t")
643 typestring = "l";
644 else if (TString(tdm->GetTypeName())=="Bool_t")
645 typestring = "O";
646 else
647 continue;
648 //printf("%s\n",tdm->GetTypeName());
649 TString varname = TString::Format("ph%d.%s",iph+1,tdm->GetName());
650
651 Char_t *addr = (Char_t*)&fDiphotonEvent->photons[iph];
652 assert(sizeof(Char_t)==1);
653
654 if (fWriteDiphotonTree)
655 hCiCTuple->Branch(varname,addr+tdm->GetOffset(),
656 TString::Format("%s/%s",varname.Data(),typestring.Data()));
657
658 if (iph==0) {
659 TString singlename = TString::Format("ph.%s",tdm->GetName());
660 Char_t *addrsingle = (Char_t*)fSinglePhoton;
661 if (fWriteSingleTree)
662 hCiCTupleSingle->Branch(singlename,addrsingle+tdm->GetOffset(),
663 TString::Format("%s/%s",singlename.Data(),typestring.Data()));
664 }
665 }
666 }
667
668 if (fWriteDiphotonTree)
669 AddOutput(hCiCTuple);
670 if (fWriteSingleTree)
671 AddOutput(hCiCTupleSingle);
672 }
673
674 // ----------------------------------------------------------------------------------------
675 // some helpfer functions....
676 void PhotonTreeWriter::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass)
677 {
678 pt = -999.;
679 decayZ = -999.;
680 mass = -999.;
681
682 // loop over all GEN particles and look for status 1 photons
683 for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
684 const MCParticle* p = fMCParticles->At(i);
685 if (p->Is(MCParticle::kH) || (!fApplyElectronVeto &&
686 (p->AbsPdgId()==23 || p->AbsPdgId()==24))) {
687 pt = p->Pt();
688 decayZ = p->DecayVertex().Z();
689 mass = p->Mass();
690 break;
691 }
692 }
693
694 return;
695 }
696
697
698 Float_t PhotonTreeWriter::GetEventCat(PhotonTools::CiCBaseLineCats cat1,
699 PhotonTools::CiCBaseLineCats cat2) {
700
701 bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
702 bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
703
704 bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
705 bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
706
707 if( ph1IsEB && ph2IsEB )
708 return ( ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
709
710 return ( ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
711 }
712
713 void PhotonTreeWriterPhoton::SetVars(const Photon *p, const DecayParticle *c, const Electron *ele,
714 const SuperCluster *pfsc, const MCParticle *m,
715 PhotonFix &phfixph, PhotonFix &phfixele,
716 const TrackCol* trackCol,const VertexCol* vtxCol,Double_t rho,
717 const ElectronCol* els, Bool_t applyElectronVeto) {
718
719 const SuperCluster *s = 0;
720 if (p)
721 s = p->SCluster();
722 else
723 s = ele->SCluster();
724 const BasicCluster *b = s->Seed();
725 const BasicCluster *b2 = 0;
726 Double_t ebcmax = -99.;
727 for (UInt_t i=0; i<s->ClusterSize(); ++i) {
728 const BasicCluster *bc = s->Cluster(i);
729 if (bc->Energy() > ebcmax && bc !=b) {
730 b2 = bc;
731 ebcmax = bc->Energy();
732 }
733 }
734
735 const BasicCluster *bclast = 0;
736 Double_t ebcmin = 1e6;
737 for (UInt_t i=0; i<s->ClusterSize(); ++i) {
738 const BasicCluster *bc = s->Cluster(i);
739 if (bc->Energy() < ebcmin && bc !=b) {
740 bclast = bc;
741 ebcmin = bc->Energy();
742 }
743 }
744
745 const BasicCluster *bclast2 = 0;
746 ebcmin = 1e6;
747 for (UInt_t i=0; i<s->ClusterSize(); ++i) {
748 const BasicCluster *bc = s->Cluster(i);
749 if (bc->Energy() < ebcmin && bc !=b && bc!=bclast) {
750 bclast2 = bc;
751 ebcmin = bc->Energy();
752 }
753 }
754
755 if (p) {
756 hasphoton = kTRUE;
757 e = p->E();
758 pt = p->Pt();
759 eta = p->Eta();
760 phi = p->Phi();
761 r9 = p->R9();
762 e3x3 = p->E33();
763 e5x5 = p->E55();
764 hovere = p->HadOverEm();
765 sigietaieta = p->CoviEtaiEta();
766 phcat = PhotonTools::CiCBaseLineCat(p);
767 eerr = p->EnergyErr();
768 eerrsmeared = p->EnergyErrSmeared();
769 esmearing = p->EnergySmearing();
770 idmva = p->IdMva();
771 hcalisodr03 = p->HcalTowerSumEtDr03();
772 ecalisodr03 = p->EcalRecHitIsoDr03();
773 trkisohollowdr03 = p->HollowConeTrkIsoDr03();
774 hcalisodr04 = p->HcalTowerSumEtDr04();
775 ecalisodr04 = p->EcalRecHitIsoDr04();
776 trkisohollowdr04 = p->HollowConeTrkIsoDr04();
777
778
779 const Vertex *vtx = vtxCol->At(0);
780 if (p->HasPV()) vtx = p->PV();
781
782 UInt_t wVtxInd = 0;
783
784 trackiso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0,
785 trackCol, NULL, NULL,
786 (!applyElectronVeto ? els : NULL) );
787 //Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
788
789 // track iso worst vtx
790 trackiso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0,
791 trackCol, &wVtxInd,vtxCol,
792 (!applyElectronVeto ? els : NULL) );
793 combiso1 = ecalisodr03+hcalisodr04+trackiso1 - 0.17*rho;
794 combiso2 = ecalisodr04+hcalisodr04+trackiso2 - 0.52*rho;
795 }
796 else {
797 hasphoton = kFALSE;
798 e = -99.;
799 pt = -99.;
800 eta = -99.;
801 phi = -99.;
802 r9 = b->E3x3()/s->RawEnergy();
803 e3x3 = b->E3x3();
804 e5x5 = b->E5x5();
805 hovere = ele->HadronicOverEm();
806 sigietaieta = ele->CoviEtaiEta();
807 phcat = -99;
808 eerr = -99.;
809 eerrsmeared = -99.;
810 esmearing = 0.;
811 idmva = -99.;
812 hcalisodr03 = -99;
813 ecalisodr03 = -99;
814 trkisohollowdr03 = -99;
815 hcalisodr04 = -99;
816 ecalisodr04 = -99;
817 trkisohollowdr04 = -99;
818 trackiso1 = -99.;
819 trackiso2 = -99.;
820 combiso1 = -99.;
821 combiso2 = -99.;
822 }
823
824 sce = s->Energy();
825 scrawe = s->RawEnergy();
826 scpse = s->PreshowerEnergy();
827 sceta = s->Eta();
828 scphi = s->Phi();
829 scnclusters = s->ClusterSize();
830 scnhits = s->NHits();
831 scetawidth = -99.;
832 scphiwidth = -99.;
833 if (p) {
834 scetawidth = p->EtaWidth();
835 scphiwidth = p->PhiWidth();
836 }
837 else {
838 scetawidth = s->EtaWidth();
839 scphiwidth = s->PhiWidth();
840 }
841 isbarrel = (s->AbsEta()<1.5);
842 isr9reco = (isbarrel && r9>0.94) || (!isbarrel && r9>0.95);
843 isr9cat = (r9>0.94);
844
845 eseed = b->Energy();
846 etaseed = b->Eta();
847 phiseed = b->Phi();
848 ietaseed = b->IEta();
849 iphiseed = b->IPhi();
850 ixseed = b->IX();
851 iyseed = b->IY();
852 etacryseed = b->EtaCry();
853 phicryseed = b->PhiCry();
854 xcryseed = b->XCry();
855 ycryseed = b->YCry();
856 thetaaxisseed = b->ThetaAxis();
857 phiaxisseed = b->PhiAxis();
858 sigietaietaseed = TMath::Sqrt(b->CoviEtaiEta());
859 sigiphiphiseed = TMath::Sqrt(b->CoviPhiiPhi());
860 if (isnan(sigiphiphiseed)) sigiphiphiseed = -99.;
861 covietaiphiseed = b->CoviEtaiPhi();
862 if (isnan(covietaiphiseed)) covietaiphiseed = -99.;
863 e3x3seed = b->E3x3();
864 e5x5seed = b->E5x5();
865 emaxseed = b->EMax();
866 e2ndseed = b->E2nd();
867 etopseed = b->ETop();
868 ebottomseed = b->EBottom();
869 eleftseed = b->ELeft();
870 erightseed = b->ERight();
871 e1x3seed = b->E1x3();
872 e3x1seed = b->E3x1();
873 e1x5seed = b->E1x5();
874 e2x2seed = b->E2x2();
875 e4x4seed = b->E4x4();
876 e2x5maxseed = b->E2x5Max();
877 e2x5topseed = b->E2x5Top();
878 e2x5bottomseed = b->E2x5Bottom();
879 e2x5leftseed = b->E2x5Left();
880 e2x5rightseed = b->E2x5Right();
881 xseedseed = b->Pos().X();
882 yseedseed = b->Pos().Y();
883 zseedseed = b->Pos().Z();
884 nhitsseed = b->NHits();
885
886 if (b2) {
887 ebc2 = b2->Energy();
888 etabc2 = b2->Eta();
889 phibc2 = b2->Phi();
890 ietabc2 = b2->IEta();
891 iphibc2 = b2->IPhi();
892 ixbc2 = b2->IX();
893 iybc2 = b2->IY();
894 etacrybc2 = b2->EtaCry();
895 phicrybc2 = b2->PhiCry();
896 xcrybc2 = b2->XCry();
897 ycrybc2 = b2->YCry();
898 thetaaxisbc2 = b2->ThetaAxis();
899 phiaxisbc2 = b2->PhiAxis();
900 sigietaietabc2 = TMath::Sqrt(b2->CoviEtaiEta());
901 sigiphiphibc2 = TMath::Sqrt(b2->CoviPhiiPhi());
902 if (isnan(sigiphiphibc2)) sigiphiphibc2 = -99.;
903 covietaiphibc2 = b2->CoviEtaiPhi();
904 if (isnan(covietaiphibc2)) covietaiphibc2 = -99.;
905 e3x3bc2 = b2->E3x3();
906 e5x5bc2 = b2->E5x5();
907 emaxbc2 = b2->EMax();
908 e2ndbc2 = b2->E2nd();
909 etopbc2 = b2->ETop();
910 ebottombc2 = b2->EBottom();
911 eleftbc2 = b2->ELeft();
912 erightbc2 = b2->ERight();
913 e1x3bc2 = b2->E1x3();
914 e3x1bc2 = b2->E3x1();
915 e1x5bc2 = b2->E1x5();
916 e2x2bc2 = b2->E2x2();
917 e4x4bc2 = b2->E4x4();
918 e2x5maxbc2 = b2->E2x5Max();
919 e2x5topbc2 = b2->E2x5Top();
920 e2x5bottombc2 = b2->E2x5Bottom();
921 e2x5leftbc2 = b2->E2x5Left();
922 e2x5rightbc2 = b2->E2x5Right();
923 xbc2bc2 = b2->Pos().X();
924 ybc2bc2 = b2->Pos().Y();
925 zbc2bc2 = b2->Pos().Z();
926 nhitsbc2= b2->NHits();
927 }
928 else {
929 ebc2 = 0;
930 etabc2 = 0;
931 phibc2 = 0;
932 ietabc2 = 0;
933 iphibc2 = 0;
934 ixbc2 = 0;
935 iybc2 = 0;
936 etacrybc2 = 0;
937 phicrybc2 = 0;
938 xcrybc2 = 0;
939 ycrybc2 = 0;
940 thetaaxisbc2 = 0;
941 phiaxisbc2 = 0;
942 sigietaietabc2 = 0;
943 sigiphiphibc2 = 0;
944 covietaiphibc2 = 0;
945 e3x3bc2 = 0;
946 e5x5bc2 = 0;
947 emaxbc2 = 0;
948 e2ndbc2 = 0;
949 etopbc2 = 0;
950 ebottombc2 = 0;
951 eleftbc2 = 0;
952 erightbc2 = 0;
953 e1x3bc2 = 0;
954 e3x1bc2 = 0;
955 e1x5bc2 = 0;
956 e2x2bc2 = 0;
957 e4x4bc2 = 0;
958 e2x5maxbc2 = 0;
959 e2x5topbc2 = 0;
960 e2x5bottombc2 = 0;
961 e2x5leftbc2 = 0;
962 e2x5rightbc2 = 0;
963 xbc2bc2 = 0;
964 ybc2bc2 = 0;
965 zbc2bc2 = 0;
966 nhitsbc2 = 0;
967 }
968
969 if (bclast) {
970 ebclast = bclast->Energy();
971 etabclast = bclast->Eta();
972 phibclast = bclast->Phi();
973 ietabclast = bclast->IEta();
974 iphibclast = bclast->IPhi();
975 ixbclast = bclast->IX();
976 iybclast = bclast->IY();
977 etacrybclast = bclast->EtaCry();
978 phicrybclast = bclast->PhiCry();
979 xcrybclast = bclast->XCry();
980 ycrybclast = bclast->YCry();
981 thetaaxisbclast = bclast->ThetaAxis();
982 phiaxisbclast = bclast->PhiAxis();
983 sigietaietabclast = TMath::Sqrt(bclast->CoviEtaiEta());
984 sigiphiphibclast = TMath::Sqrt(bclast->CoviPhiiPhi());
985 if (isnan(sigiphiphibclast)) sigiphiphibclast = -99.;
986 covietaiphibclast = bclast->CoviEtaiPhi();
987 if (isnan(covietaiphibclast)) covietaiphibclast = -99.;
988 e3x3bclast = bclast->E3x3();
989 nhitsbclast = bclast->NHits();
990 }
991 else {
992 ebclast = 0;
993 etabclast = 0;
994 phibclast = 0;
995 ietabclast = 0;
996 iphibclast = 0;
997 ixbclast = 0;
998 iybclast = 0;
999 etacrybclast = 0;
1000 phicrybclast = 0;
1001 xcrybclast = 0;
1002 ycrybclast = 0;
1003 thetaaxisbclast = 0;
1004 phiaxisbclast = 0;
1005 sigietaietabclast = 0;
1006 sigiphiphibclast = 0;
1007 covietaiphibclast = 0;
1008 e3x3bclast = 0;
1009 nhitsbclast = 0;
1010 }
1011
1012 if (bclast2) {
1013 ebclast2 = bclast2->Energy();
1014 etabclast2 = bclast2->Eta();
1015 phibclast2 = bclast2->Phi();
1016 ietabclast2 = bclast2->IEta();
1017 iphibclast2 = bclast2->IPhi();
1018 ixbclast2 = bclast2->IX();
1019 iybclast2 = bclast2->IY();
1020 etacrybclast2 = bclast2->EtaCry();
1021 phicrybclast2 = bclast2->PhiCry();
1022 xcrybclast2 = bclast2->XCry();
1023 ycrybclast2 = bclast2->YCry();
1024 thetaaxisbclast2 = bclast2->ThetaAxis();
1025 phiaxisbclast2 = bclast2->PhiAxis();
1026 sigietaietabclast2 = TMath::Sqrt(bclast2->CoviEtaiEta());
1027 sigiphiphibclast2 = TMath::Sqrt(bclast2->CoviPhiiPhi());
1028 if (isnan(sigiphiphibclast2)) sigiphiphibclast2 = -99.;
1029 covietaiphibclast2 = bclast2->CoviEtaiPhi();
1030 if (isnan(covietaiphibclast2)) covietaiphibclast2 = -99.;
1031 e3x3bclast2 = bclast2->E3x3();
1032 nhitsbclast2 = bclast2->NHits();
1033 }
1034 else {
1035 ebclast2 = 0;
1036 etabclast2 = 0;
1037 phibclast2 = 0;
1038 ietabclast2 = 0;
1039 iphibclast2 = 0;
1040 ixbclast2 = 0;
1041 iybclast2 = 0;
1042 etacrybclast2 = 0;
1043 phicrybclast2 = 0;
1044 xcrybclast2 = 0;
1045 ycrybclast2 = 0;
1046 thetaaxisbclast2 = 0;
1047 phiaxisbclast2 = 0;
1048 sigietaietabclast2 = 0;
1049 sigiphiphibclast2 = 0;
1050 covietaiphibclast2 = 0;
1051 e3x3bclast2 = 0;
1052 nhitsbclast2 = 0;
1053 }
1054
1055 //initialize photon energy corrections if needed
1056 /*if (!PhotonFix::initialised()) {
1057 PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
1058 }*/
1059
1060 phfixph.setup(e,sceta,scphi,r9);
1061 phfixele.setup(e,sceta,scphi,r9);
1062
1063 const Float_t dval = -99.;
1064 ecor = phfixph.fixedEnergy();
1065 ecorerr = phfixph.sigmaEnergy();
1066 ecorele = phfixele.fixedEnergy();
1067 ecoreleerr = phfixele.sigmaEnergy();
1068 if (phfixph.isbarrel()) {
1069 etac = phfixph.etaC();
1070 etas = phfixph.etaS();
1071 etam = phfixph.etaM();
1072 phic = phfixph.phiC();
1073 phis = phfixph.phiS();
1074 phim = phfixph.phiM();
1075 xz = dval;
1076 xc = dval;
1077 xs = dval;
1078 xm = dval;
1079 yz = dval;
1080 yc = dval;
1081 ys = dval;
1082 ym = dval;
1083 }
1084 else {
1085 etac = dval;
1086 etas = dval;
1087 etam = dval;
1088 phic = dval;
1089 phis = dval;
1090 phim = dval;
1091 xz = phfixph.xZ();
1092 xc = phfixph.xC();
1093 xs = phfixph.xS();
1094 xm = phfixph.xM();
1095 yz = phfixph.yZ();
1096 yc = phfixph.yC();
1097 ys = phfixph.yS();
1098 ym = phfixph.yM();
1099 }
1100
1101 if (c) {
1102 hasconversion = kTRUE;
1103 convp = c->P();
1104 convpt = c->Pt();
1105 conveta = c->Eta();
1106 convphi = c->Phi();
1107 ThreeVector dirconvsc = ThreeVector(s->Point()) - c->Position();
1108 convdeta = c->Eta() - dirconvsc.Eta();
1109 convdphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
1110 convvtxrho = c->Position().Rho();
1111 convvtxz = c->Position().Z();
1112 convvtxphi = c->Position().Phi();
1113
1114 const StableData *leadsd = dynamic_cast<const StableData*>(c->DaughterDat(0));
1115 const StableData *trailsd = dynamic_cast<const StableData*>(c->DaughterDat(1));
1116 if (leadsd->Pt()<trailsd->Pt()) {
1117 const StableData *sdtmp = leadsd;
1118 leadsd = trailsd;
1119 trailsd = sdtmp;
1120 }
1121
1122 const Track *leadtrack = dynamic_cast<const StableParticle*>(leadsd->Original())->Trk();
1123 const Track *trailtrack = dynamic_cast<const StableParticle*>(trailsd->Original())->Trk();
1124
1125 convleadpt = leadsd->Pt();
1126 convtrailpt = trailsd->Pt();
1127 convleadtrackpt = leadtrack->Pt();
1128 convleadtrackalgo = leadtrack->Algo();
1129 if (convleadtrackalgo==29) convleadtrackalgos=2; //gsf track
1130 else if (convleadtrackalgo==15 ||convleadtrackalgo==16) convleadtrackalgos=1; //ecal-seeded track
1131 else convleadtrackalgos = 0; //std iterative track
1132 convleadtrackcharge = leadtrack->Charge();
1133 convtrailtrackpt = trailtrack->Pt();
1134 convtrailtrackalgo = trailtrack->Algo();
1135 if (convtrailtrackalgo==29) convtrailtrackalgos=2; //gsf track
1136 else if (convtrailtrackalgo==15 ||convtrailtrackalgo==16) convtrailtrackalgos=1; //ecal-seeded track
1137 else convtrailtrackalgos = 0; //std iterative track
1138 trailtrackcharge = trailtrack->Charge();
1139 }
1140 else {
1141 hasconversion = kFALSE;
1142 convp = -99.;
1143 convpt = -99.;
1144 conveta = -99.;
1145 convphi = -99.;
1146 convdeta = -99.;
1147 convdphi = -99.;
1148 convvtxrho = -99.;
1149 convvtxz = -999.;
1150 convvtxphi = -99.;
1151 convleadpt = -99.;
1152 convtrailpt = -99.;
1153 convleadtrackpt = -99.;
1154 convleadtrackalgo = -99;
1155 convleadtrackalgos = -99;
1156 convleadtrackcharge = 0;
1157 convtrailtrackpt = -99.;
1158 convtrailtrackalgo = -99;
1159 convtrailtrackalgos = -99;
1160 trailtrackcharge = 0;
1161 }
1162
1163 //electron quantities
1164 if (ele) {
1165 haselectron = kTRUE;
1166 eleisecaldriven = ele->IsEcalDriven();
1167 eleistrackerdriven = ele->IsTrackerDriven();
1168 elee = ele->E();
1169 elept = ele->Pt();
1170 eleeta = ele->Eta();
1171 elephi = ele->Phi();
1172 elecharge = ele->Charge();
1173 elefbrem = ele->FBrem();
1174 eledeta = ele->DeltaEtaSuperClusterTrackAtVtx();
1175 eledphi = ele->DeltaPhiSuperClusterTrackAtVtx();
1176 elep = s->Energy()/ele->ESuperClusterOverP();
1177 elepin = ele->PIn();
1178 elepout = ele->POut();
1179 }
1180 else {
1181 haselectron = kFALSE;
1182 eleisecaldriven = kFALSE;
1183 eleistrackerdriven = kFALSE;
1184 elee = -99.;
1185 elept = -99.;
1186 eleeta = -99.;
1187 elephi = -99.;
1188 elecharge = -99;
1189 elefbrem = -99.;
1190 eledeta = -99.;
1191 eledphi = -99.;
1192 elep = -99.;
1193 elepin = -99.;
1194 elepout = -99.;
1195 }
1196
1197 //pf supercluster quantities
1198 if (pfsc) {
1199 haspfsc = kTRUE;
1200 pfsce = pfsc->Energy();
1201 pfscrawe = pfsc->RawEnergy();
1202 pfsceta = pfsc->Eta();
1203 pfscphi = pfsc->Phi();
1204 }
1205 else {
1206 haspfsc = kFALSE;
1207 pfsce = -99.;
1208 pfscrawe = -99.;
1209 pfsceta = -99.;
1210 pfscphi = -99.;
1211 }
1212
1213 genz = -99.;
1214 if (m) {
1215 ispromptgen = kTRUE;
1216 gene = m->E();
1217 genpt = m->Pt();
1218 geneta = m->Eta();
1219 genphi = m->Phi();
1220 const MCParticle *mm = m->DistinctMother();
1221 if (mm) genz = mm->DecayVertex().Z();
1222 pdgid = m->PdgId();
1223 if (mm) motherpdgid = mm->PdgId();
1224 else motherpdgid = -99;
1225 }
1226 else {
1227 ispromptgen = kFALSE;
1228 gene = -99.;
1229 genpt = -99.;
1230 geneta = -99.;
1231 genphi = -99.;
1232 pdgid = -99;
1233 motherpdgid = -99;
1234 }
1235 }