ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.14
Committed: Wed May 2 16:57:19 2012 UTC (13 years ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.13: +74 -0 lines
Log Message:
fixed/updated LeptonTag

File Contents

# User Rev Content
1 bendavid 1.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 bendavid 1.2 #include "MitAna/DataTree/interface/PFMet.h"
7 bendavid 1.1 #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 paus 1.12 BaseMod (name,title),
26 bendavid 1.1 // define all the Branches to load
27 paus 1.12 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 fabstoec 1.13
42     fLeptonTagElectronsName ("HggLeptonTagElectrons"),
43     fLeptonTagMuonsName ("HggLeptonTagMuons"),
44    
45 paus 1.12 fIsData (false),
46     fPhotonsFromBranch (kTRUE),
47     fPVFromBranch (kTRUE),
48     fGoodElectronsFromBranch(kTRUE),
49     fPFJetsFromBranch (kTRUE),
50 bendavid 1.1 // ----------------------------------------
51     // collections....
52 paus 1.12 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 fabstoec 1.13
65     fLeptonTagElectrons (0),
66     fLeptonTagMuons (0),
67    
68 paus 1.12 fLoopOnGoodElectrons (kFALSE),
69     fApplyElectronVeto (kTRUE),
70     fWriteDiphotonTree (kTRUE),
71     fWriteSingleTree (kTRUE),
72     fExcludeSinglePrompt (kFALSE),
73     fExcludeDoublePrompt (kFALSE),
74     fEnableJets (kFALSE),
75 fabstoec 1.13 fApplyLeptonTag (kFALSE),
76 paus 1.12 fPhFixDataFile (gSystem->Getenv("CMSSW_BASE") +
77     TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat")),
78     fTupleName ("hPhotonTree")
79 bendavid 1.1 {
80 paus 1.12 // Constructor
81 bendavid 1.1 }
82    
83 paus 1.12 PhotonTreeWriter::~PhotonTreeWriter()
84     {
85     // Destructor
86 bendavid 1.1 }
87    
88     //--------------------------------------------------------------------------------------------------
89     void PhotonTreeWriter::Process()
90     {
91     // ------------------------------------------------------------
92     // Process entries of the tree.
93     LoadEventObject(fPhotonBranchName, fPhotons);
94 paus 1.12 LoadEventObject(fGoodElectronName, fGoodElectrons);
95 fabstoec 1.13
96     // lepton tag collections
97     if( fApplyLeptonTag ) {
98     LoadEventObject(fLeptonTagElectronsName, fLeptonTagElectrons);
99     LoadEventObject(fLeptonTagMuonsName, fLeptonTagMuons);
100     }
101    
102 bendavid 1.2 const BaseCollection *egcol = 0;
103 paus 1.12 if (fLoopOnGoodElectrons)
104 bendavid 1.2 egcol = fGoodElectrons;
105 paus 1.12 else
106 bendavid 1.2 egcol = fPhotons;
107 paus 1.12 if (egcol->GetEntries()<1)
108     return;
109 bendavid 1.1 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 bendavid 1.2 LoadEventObject(fSuperClusterName, fSuperClusters);
117 paus 1.12 LoadEventObject(fPFMetName, fPFMet);
118     if (fEnableJets)
119     LoadEventObject(fPFJetName, fPFJets);
120 bendavid 1.1
121     // ------------------------------------------------------------
122     // load event based information
123 paus 1.12 Int_t _numPU = -1.; // some sensible default values....
124 bendavid 1.1 Int_t _numPUminus = -1.; // some sensible default values....
125 paus 1.12 Int_t _numPUplus = -1.; // some sensible default values....
126 bendavid 1.1
127 paus 1.12 Float_t rho = -99.;
128 bendavid 1.1 if( fPileUpDen->GetEntries() > 0 )
129 paus 1.12 rho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
130 bendavid 1.1
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 paus 1.12 Float_t _pth = -100.;
149     Float_t _decayZ = -100.;
150 bendavid 1.1 Float_t _genmass = -100.;
151 paus 1.12 if (!fIsData)
152     FindHiggsPtAndZ(_pth, _decayZ, _genmass);
153 bendavid 1.1
154 fabstoec 1.13 fDiphotonEvent->leptonTag = -1; // disabled
155    
156 paus 1.12 fDiphotonEvent->rho = rho;
157 bendavid 1.1 fDiphotonEvent->genHiggspt = _pth;
158     fDiphotonEvent->genHiggsZ = _decayZ;
159     fDiphotonEvent->genmass = _genmass;
160     fDiphotonEvent->gencostheta = -99.;
161     fDiphotonEvent->nVtx = fPV->GetEntries();
162 bendavid 1.2 fDiphotonEvent->bsX = fBeamspot->At(0)->X();
163     fDiphotonEvent->bsY = fBeamspot->At(0)->Y();
164     fDiphotonEvent->bsZ = fBeamspot->At(0)->Z();
165 bendavid 1.7 fDiphotonEvent->bsSigmaZ = fBeamspot->At(0)->SigmaZ();
166 bendavid 1.2 fDiphotonEvent->vtxX = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->X() : -99.;
167     fDiphotonEvent->vtxY = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Y() : -99.;
168 bendavid 1.1 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 bendavid 1.2 fDiphotonEvent->mt = -99.;
176     fDiphotonEvent->cosphimet = -99.;
177     fDiphotonEvent->mtele = -99.;
178     fDiphotonEvent->cosphimetele = -99.;
179 bendavid 1.1 fDiphotonEvent->evt = GetEventHeader()->EvtNum();
180     fDiphotonEvent->run = GetEventHeader()->RunNum();
181     fDiphotonEvent->lumi = GetEventHeader()->LumiSec();
182     fDiphotonEvent->evtcat = -99;
183 bendavid 1.2 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 bendavid 1.1 fDiphotonEvent->masscor = -99.;
189     fDiphotonEvent->masscorerr = -99.;
190 bendavid 1.2 fDiphotonEvent->masscorele = -99.;
191     fDiphotonEvent->masscoreleerr = -99.;
192     fDiphotonEvent->ismc = GetEventHeader()->IsMC();
193    
194 bendavid 1.5 //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 bendavid 1.1 Int_t nhitsbeforevtxmax = 1;
221 paus 1.12 if (!fApplyElectronVeto)
222     nhitsbeforevtxmax = 999;
223 bendavid 1.1
224 bendavid 1.2 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 bendavid 1.1 }
254    
255 paus 1.12 const DecayParticle *conv1 = PhotonTools::MatchedConversion(sc1,fConversions,bsp,
256     nhitsbeforevtxmax);
257     const DecayParticle *conv2 = PhotonTools::MatchedConversion(sc2,fConversions,bsp,
258     nhitsbeforevtxmax);
259 bendavid 1.2
260     const SuperCluster *pfsc1 = PhotonTools::MatchedSC(sc1,fSuperClusters);
261     const SuperCluster *pfsc2 = PhotonTools::MatchedSC(sc2,fSuperClusters);
262 bendavid 1.1
263     const MCParticle *phgen1 = 0;
264     const MCParticle *phgen2 = 0;
265 paus 1.12 if (!fIsData) {
266 bendavid 1.6 phgen1 = PhotonTools::MatchMC(p1,fMCParticles,!fApplyElectronVeto);
267     phgen2 = PhotonTools::MatchMC(p2,fMCParticles,!fApplyElectronVeto);
268 bendavid 1.2 }
269    
270 paus 1.12 if (fExcludeSinglePrompt && (phgen1 || phgen2))
271     return;
272     if (fExcludeDoublePrompt && (phgen1 && phgen2))
273     return;
274 bendavid 1.2
275     if (!fLoopOnGoodElectrons && phHard->HasPV()) {
276     fDiphotonEvent->vtxX = phHard->PV()->X();
277     fDiphotonEvent->vtxY = phHard->PV()->Y();
278     fDiphotonEvent->vtxZ = phHard->PV()->Z();
279 bendavid 1.4 fDiphotonEvent->vtxprob = phHard->VtxProb();
280 bendavid 1.2 }
281    
282 bendavid 1.5 //fill jet variables
283     if (fEnableJets) {
284     for (UInt_t ijet=0; ijet<fPFJets->GetEntries();++ijet) {
285     const Jet *jet = fPFJets->At(ijet);
286 bendavid 1.8 if (jet->AbsEta()<4.7 && MathUtils::DeltaR(jet,p1)>0.5 && MathUtils::DeltaR(jet,p2)>0.5) {
287 bendavid 1.5 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 bendavid 1.4 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 bendavid 1.5 Double_t _etagg = -99.;
341     Double_t _phigg = -99.;
342 bendavid 1.4 Double_t _costheta = -99.;
343 bendavid 1.2 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 bendavid 1.5 _etagg = (phHard->Mom()+phSoft->Mom()).Eta();
350     _phigg = (phHard->Mom()+phSoft->Mom()).Phi();
351 bendavid 1.2 _costheta = ThreeVector(phHard->Mom()).Unit().Dot(ThreeVector(phSoft->Mom()).Unit());
352     _evtcat = PhotonTools::DiphotonR9EtaPtCat(phHard,phSoft);
353 bendavid 1.4
354 bendavid 1.5 const Double_t dz = sqrt(2.0)*5.8;
355 paus 1.12 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 bendavid 1.5 fDiphotonEvent->deltamvtx = deltamvtx;
366    
367 paus 1.12 _masserrwrongvtx = TMath::Sqrt(_masserr*_masserr + deltamvtx*deltamvtx);
368 bendavid 1.4 _masserrsmearedwrongvtx = TMath::Sqrt(_masserrsmeared*_masserrsmeared + deltamvtx*deltamvtx);
369    
370 bendavid 1.5 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 bendavid 1.4
375 bendavid 1.2 }
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 bendavid 1.1 }
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 bendavid 1.2 fDiphotonEvent->masserr = _masserr;
395     fDiphotonEvent->masserrsmeared = _masserrsmeared;
396 bendavid 1.4 fDiphotonEvent->masserrwrongvtx = _masserrwrongvtx;
397     fDiphotonEvent->masserrsmearedwrongvtx = _masserrsmearedwrongvtx;
398 bendavid 1.1 fDiphotonEvent->ptgg = _ptgg;
399 bendavid 1.5 fDiphotonEvent->etagg = _etagg;
400     fDiphotonEvent->phigg = _phigg;
401 bendavid 1.2 fDiphotonEvent->costheta = _costheta;;
402     fDiphotonEvent->massele = _massele;
403     fDiphotonEvent->ptee = _ptee;
404     fDiphotonEvent->costhetaele = _costhetaele;
405     fDiphotonEvent->evtcat = _evtcat;
406 bendavid 1.1
407 paus 1.12 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 bendavid 1.1
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 bendavid 1.2
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 bendavid 1.1
421 paus 1.12 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 bendavid 1.2
431 paus 1.12 //printf("r9 = %5f, photon sigieie = %5f, seed sigieie = %5f\n",phHard->R9(),
432     // phHard->CoviEtaiEta(),sqrt(phHard->SCluster()->Seed()->CoviEtaiEta()));
433 fabstoec 1.13
434 fabstoec 1.14 // MuonStuff
435     fDiphotonEvent-> muonPt = -99.;
436     fDiphotonEvent-> muonEta = -99.;
437     fDiphotonEvent-> muDR1 = -99.;
438     fDiphotonEvent-> muDR2 = -99.;
439     fDiphotonEvent-> muIso1 = -99.;
440     fDiphotonEvent-> muIso2 = -99.;
441     fDiphotonEvent-> muIso3 = -99.;
442     fDiphotonEvent-> muIso4 = -99.;
443     fDiphotonEvent-> muD0 = -99.;
444     fDiphotonEvent-> muDZ = -99.;
445     fDiphotonEvent-> muChi2 = -99.;
446     fDiphotonEvent-> muNhits = -99;
447     fDiphotonEvent-> muNpixhits = -99;
448     fDiphotonEvent-> muNegs = -99;
449     fDiphotonEvent-> muNMatch = -99;
450    
451     // Electron Stuff
452     fDiphotonEvent-> elePt = -99.;
453     fDiphotonEvent-> eleEta = -99.;
454     fDiphotonEvent-> eleSCEta = -99.;
455     fDiphotonEvent-> eleIso1 = -99.;
456     fDiphotonEvent-> eleIso2 = -99.;
457     fDiphotonEvent-> eleIso3 = -99.;
458     fDiphotonEvent-> eleIso4 = -99.;
459     fDiphotonEvent-> eleDist = -99.;
460     fDiphotonEvent-> eleDcot = -99.;
461     fDiphotonEvent-> eleCoviee = -99.;
462     fDiphotonEvent-> eleDphiin = -99.;
463     fDiphotonEvent-> eleDetain = -99.;
464     fDiphotonEvent-> eleDR1 = -99.;
465     fDiphotonEvent-> eleDR2 = -99.;
466     fDiphotonEvent-> eleMass1 = -99.;
467     fDiphotonEvent-> eleMass2 = -99.;
468     fDiphotonEvent-> eleNinnerHits = -99;
469    
470 fabstoec 1.13 if( fApplyLeptonTag ) {
471     // perform lepton tagging
472     // the diphoton event record will have one more entry; i.e. leptonTag
473     // leptonTag = -1 -> lepton-taggng was swicthed off
474     // = 0 -> event tagged as 'non-lepton-event'
475     // = +1 -> event tagged as muon-event
476     // = +2 -> event tagged as electron-event
477     fDiphotonEvent->leptonTag = 0;
478    
479     if ( fLeptonTagMuons->GetEntries() > 0 ) {
480 fabstoec 1.14
481 fabstoec 1.13 // need to have dR > 1 for with respect to both photons
482     for(UInt_t iMuon = 0; iMuon <fLeptonTagMuons->GetEntries(); ++iMuon) {
483     if(MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phHard) < 1.) continue;
484     if(MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phSoft) < 1.) continue;
485    
486     fDiphotonEvent->leptonTag = 1;
487 fabstoec 1.14
488     fDiphotonEvent-> muonPt = fLeptonTagMuons->At(iMuon)->Pt();
489     fDiphotonEvent-> muonEta = fLeptonTagMuons->At(iMuon)->Eta();
490     fDiphotonEvent-> muDR1 = MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phHard);
491     fDiphotonEvent-> muDR2 = MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phSoft);
492     fDiphotonEvent-> muIso1 = (fLeptonTagMuons->At(iMuon)->IsoR03SumPt() + fLeptonTagMuons->At(iMuon)->IsoR03EmEt() + fLeptonTagMuons->At(iMuon)->IsoR03HadEt() - fPileUpDen->At(0)->RhoRandomLowEta() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(iMuon)->Pt();
493     fDiphotonEvent-> muIso2 = (fLeptonTagMuons->At(iMuon)->IsoR03SumPt() + fLeptonTagMuons->At(iMuon)->IsoR03EmEt() + fLeptonTagMuons->At(iMuon)->IsoR03HadEt() - fPileUpDen->At(0)->RhoRandom() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(iMuon)->Pt();
494     fDiphotonEvent-> muIso3 = (fLeptonTagMuons->At(iMuon)->IsoR03SumPt() + fLeptonTagMuons->At(iMuon)->IsoR03EmEt() + fLeptonTagMuons->At(iMuon)->IsoR03HadEt() - fPileUpDen->At(0)->RhoLowEta() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(iMuon)->Pt();
495     fDiphotonEvent-> muIso4 = (fLeptonTagMuons->At(iMuon)->IsoR03SumPt() + fLeptonTagMuons->At(iMuon)->IsoR03EmEt() + fLeptonTagMuons->At(iMuon)->IsoR03HadEt() - fPileUpDen->At(0)->Rho() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(iMuon)->Pt();
496     fDiphotonEvent-> muD0 = TMath::Abs(fLeptonTagMuons->At(iMuon)->BestTrk()->D0Corrected(*fPV->At(0)));
497     fDiphotonEvent-> muDZ = TMath::Abs(fLeptonTagMuons->At(iMuon)->BestTrk()->DzCorrected(*fPV->At(0)));
498     fDiphotonEvent-> muChi2 = fLeptonTagMuons->At(iMuon)->GlobalTrk()->Chi2()/fLeptonTagMuons->At(iMuon)->GlobalTrk()->Ndof();
499    
500     fDiphotonEvent-> muNhits = fLeptonTagMuons->At(iMuon)->BestTrk()->NHits();
501     fDiphotonEvent-> muNpixhits = fLeptonTagMuons->At(iMuon)->BestTrk()->NPixelHits();
502     fDiphotonEvent-> muNegs = fLeptonTagMuons->At(iMuon)->NSegments();
503     fDiphotonEvent-> muNMatch = fLeptonTagMuons->At(iMuon)->NMatches();
504    
505 fabstoec 1.13 break;
506     }
507     }
508     if ( fDiphotonEvent->leptonTag < 1 && fLeptonTagElectrons->GetEntries() > 0 ) {
509     for(UInt_t iEle = 0; iEle < fLeptonTagElectrons->GetEntries(); ++iEle) {
510     if(MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phHard) < 1.) continue;
511     if(MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phSoft) < 1.) continue;
512    
513     // here we also check the mass ....
514     if ( TMath::Abs( (phHard->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M()-91.19 ) < 5. ) continue;
515     if ( TMath::Abs( (phSoft->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M()-91.19 ) < 5. ) continue;
516    
517     fDiphotonEvent->leptonTag = 2;
518 fabstoec 1.14
519     fDiphotonEvent-> elePt = fLeptonTagElectrons->At(iEle)->Pt();
520     fDiphotonEvent-> eleEta = fLeptonTagElectrons->At(iEle)->Eta();
521     fDiphotonEvent-> eleSCEta = fLeptonTagElectrons->At(iEle)->SCluster()->Eta();
522     fDiphotonEvent-> eleIso1 = (fLeptonTagElectrons->At(iEle)->TrackIsolationDr03() + fLeptonTagElectrons->At(iEle)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(iEle)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->RhoRandomLowEta() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
523     fDiphotonEvent-> eleIso2 = (fLeptonTagElectrons->At(iEle)->TrackIsolationDr03() + fLeptonTagElectrons->At(iEle)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(iEle)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->RhoRandom() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
524     fDiphotonEvent-> eleIso3 = (fLeptonTagElectrons->At(iEle)->TrackIsolationDr03() + fLeptonTagElectrons->At(iEle)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(iEle)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->RhoLowEta() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
525     fDiphotonEvent-> eleIso4 = (fLeptonTagElectrons->At(iEle)->TrackIsolationDr03() + fLeptonTagElectrons->At(iEle)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(iEle)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->Rho() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
526     fDiphotonEvent-> eleDist = fLeptonTagElectrons->At(iEle)->ConvPartnerDist();
527     fDiphotonEvent-> eleDcot = fLeptonTagElectrons->At(iEle)->ConvPartnerDCotTheta();
528     fDiphotonEvent-> eleCoviee = fLeptonTagElectrons->At(iEle)->CoviEtaiEta();
529     fDiphotonEvent-> eleDphiin = TMath::Abs(fLeptonTagElectrons->At(iEle)->DeltaPhiSuperClusterTrackAtVtx());
530     fDiphotonEvent-> eleDetain = TMath::Abs(fLeptonTagElectrons->At(iEle)->DeltaEtaSuperClusterTrackAtVtx());
531     fDiphotonEvent-> eleDR1 = MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phHard);
532     fDiphotonEvent-> eleDR2 = MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phSoft);
533     fDiphotonEvent-> eleMass1 = (phHard->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M();
534     fDiphotonEvent-> eleMass2 = (phSoft->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M();
535     fDiphotonEvent-> eleNinnerHits = fLeptonTagElectrons->At(iEle)->Trk()->NExpectedHitsInner();
536    
537 fabstoec 1.13 break;
538     }
539     }
540     }
541    
542 paus 1.12 if (fWriteDiphotonTree)
543     hCiCTuple->Fill();
544 bendavid 1.1 }
545    
546 paus 1.12 if (!fWriteSingleTree)
547     return;
548 bendavid 1.1
549 bendavid 1.2 for (UInt_t iph = 0; iph<egcol->GetEntries(); ++iph) {
550     const Particle *p = 0;
551     const Photon *ph = 0;
552     const Electron *ele = 0;
553     const SuperCluster *sc = 0;
554     if (fLoopOnGoodElectrons) {
555     ele = fGoodElectrons->At(iph);
556     p = ele;
557     sc = ele->SCluster();
558     ph = PhotonTools::MatchedPhoton(ele,fPhotons);
559     }
560     else {
561     ph = fPhotons->At(iph);
562     p = ph;
563     sc = ph->SCluster();
564     ele = PhotonTools::MatchedElectron(ph,fGoodElectrons);
565     }
566    
567 paus 1.12 const DecayParticle *conv = PhotonTools::MatchedConversion(sc,fConversions,bsp,
568     nhitsbeforevtxmax);
569     const SuperCluster *pfsc = PhotonTools::MatchedSC(sc,fSuperClusters);
570 bendavid 1.2
571     if (!fLoopOnGoodElectrons && ph->HasPV()) {
572 bendavid 1.1 fDiphotonEvent->vtxZ = ph->PV()->Z();
573     }
574    
575     const MCParticle *phgen = 0;
576     if( !fIsData ) {
577 bendavid 1.6 phgen = PhotonTools::MatchMC(p,fMCParticles,!fApplyElectronVeto);
578 bendavid 1.1 }
579    
580 bendavid 1.2 if (fExcludeSinglePrompt && phgen) return;
581    
582     fDiphotonEvent->mt = -99.;
583     fDiphotonEvent->cosphimet = -99.;
584     fDiphotonEvent->mtele = -99.;
585     fDiphotonEvent->cosphimetele = -99.;
586    
587     if (ph) {
588     fDiphotonEvent->cosphimet = TMath::Cos(ph->Phi()-fPFMet->At(0)->Phi());
589 paus 1.12 fDiphotonEvent->mt = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ph->Pt()*
590     (1.0-fDiphotonEvent->cosphimet));
591 bendavid 1.2 }
592 bendavid 1.1
593 bendavid 1.2 if (ele) {
594     fDiphotonEvent->cosphimetele = TMath::Cos(ele->Phi()-fPFMet->At(0)->Phi());
595 paus 1.12 fDiphotonEvent->mtele = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ele->Pt()*
596     (1.0-fDiphotonEvent->cosphimetele));
597 bendavid 1.2 }
598    
599 paus 1.12 fSinglePhoton->SetVars(ph,conv,ele,pfsc,phgen,fPhfixph,fPhfixele,fTracks,fPV,rho,
600     fElectrons,fApplyElectronVeto);
601 bendavid 1.1 hCiCTupleSingle->Fill();
602     }
603    
604     return;
605     }
606    
607     //--------------------------------------------------------------------------------------------------
608     void PhotonTreeWriter::SlaveBegin()
609     {
610     // Run startup code on the computer (slave) doing the actual analysis. Here,
611     // we just request the photon collection branch.
612    
613 fabstoec 1.13 if( fApplyLeptonTag ) {
614     ReqEventObject(fLeptonTagElectronsName, fLeptonTagElectrons, false);
615     ReqEventObject(fLeptonTagMuonsName, fLeptonTagMuons, false);
616     }
617    
618 paus 1.12 ReqEventObject(fPhotonBranchName,fPhotons, fPhotonsFromBranch);
619     ReqEventObject(fTrackBranchName, fTracks, true);
620     ReqEventObject(fElectronName, fElectrons, true);
621     ReqEventObject(fGoodElectronName,fGoodElectrons,fGoodElectronsFromBranch);
622     ReqEventObject(fPileUpDenName, fPileUpDen, true);
623     ReqEventObject(fPVName, fPV, fPVFromBranch);
624     ReqEventObject(fConversionName, fConversions, true);
625     ReqEventObject(fBeamspotName, fBeamspot, true);
626     ReqEventObject(fPFCandName, fPFCands, true);
627     ReqEventObject(fSuperClusterName,fSuperClusters,true);
628     ReqEventObject(fPFMetName, fPFMet, true);
629     if (fEnableJets)
630     ReqEventObject(fPFJetName, fPFJets, fPFJetsFromBranch);
631 bendavid 1.1 if (!fIsData) {
632 paus 1.12 ReqBranch(fPileUpName, fPileUp);
633     ReqBranch(fMCParticleName, fMCParticles);
634 bendavid 1.1 }
635 bendavid 1.6 if (fIsData) {
636 paus 1.12 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
637     TString("/src/MitPhysics/data/PhotonFixGRPV22.dat");
638 bendavid 1.6 }
639     else {
640 paus 1.12 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
641     TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat");
642 bendavid 1.6 }
643 bendavid 1.1
644 bendavid 1.2 fPhfixph.initialise("4_2",std::string(fPhFixDataFile));
645     fPhfixele.initialise("4_2e",std::string(fPhFixDataFile));
646    
647 bendavid 1.1 fDiphotonEvent = new PhotonTreeWriterDiphotonEvent;
648 paus 1.12 fSinglePhoton = new PhotonTreeWriterPhoton;
649 bendavid 1.1
650 paus 1.12 if (fWriteDiphotonTree)
651     hCiCTuple = new TTree(fTupleName.Data(),fTupleName.Data());
652 bendavid 1.1 TString singlename = fTupleName + TString("Single");
653 paus 1.12 if (fWriteSingleTree)
654     hCiCTupleSingle = new TTree(singlename,singlename);
655 bendavid 1.1
656     //make flattish tree from classes so we don't have to rely on dictionaries for reading later
657     TClass *eclass = TClass::GetClass("mithep::PhotonTreeWriterDiphotonEvent");
658     TClass *pclass = TClass::GetClass("mithep::PhotonTreeWriterPhoton");
659 paus 1.12 TList *elist = eclass->GetListOfDataMembers();
660     TList *plist = pclass->GetListOfDataMembers();
661 bendavid 1.1
662     for (int i=0; i<elist->GetEntries(); ++i) {
663     const TDataMember *tdm = static_cast<const TDataMember*>(elist->At(i));
664 paus 1.12 if (!(tdm->IsBasic() && tdm->IsPersistent()))
665     continue;
666 bendavid 1.1 TString typestring;
667 paus 1.12 if (TString(tdm->GetTypeName())=="Char_t") typestring = "B";
668 bendavid 1.1 else if (TString(tdm->GetTypeName())=="UChar_t") typestring = "b";
669     else if (TString(tdm->GetTypeName())=="Short_t") typestring = "S";
670     else if (TString(tdm->GetTypeName())=="UShort_t") typestring = "s";
671     else if (TString(tdm->GetTypeName())=="Int_t") typestring = "I";
672     else if (TString(tdm->GetTypeName())=="UInt_t") typestring = "i";
673     else if (TString(tdm->GetTypeName())=="Float_t") typestring = "F";
674     else if (TString(tdm->GetTypeName())=="Double_t") typestring = "D";
675     else if (TString(tdm->GetTypeName())=="Long64_t") typestring = "L";
676     else if (TString(tdm->GetTypeName())=="ULong64_t") typestring = "l";
677     else if (TString(tdm->GetTypeName())=="Bool_t") typestring = "O";
678 paus 1.12 else
679     continue;
680 bendavid 1.1 //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
681     Char_t *addr = (Char_t*)fDiphotonEvent;
682     assert(sizeof(Char_t)==1);
683 paus 1.12
684     if (fWriteDiphotonTree)
685     hCiCTuple->Branch(tdm->GetName(),addr + tdm->GetOffset(),
686     TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
687     if (fWriteSingleTree)
688     hCiCTupleSingle->Branch(tdm->GetName(),addr + tdm->GetOffset(),
689     TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
690 bendavid 1.1 }
691    
692     for (int iph=0; iph<2; ++iph) {
693     for (int i=0; i<plist->GetEntries(); ++i) {
694     const TDataMember *tdm = static_cast<const TDataMember*>(plist->At(i));
695 paus 1.12 if (!(tdm->IsBasic() && tdm->IsPersistent()))
696     continue;
697 bendavid 1.1 TString typestring;
698 paus 1.12 if (TString(tdm->GetTypeName())=="Char_t")
699     typestring = "B";
700     else if (TString(tdm->GetTypeName())=="UChar_t")
701     typestring = "b";
702     else if (TString(tdm->GetTypeName())=="Short_t")
703     typestring = "S";
704     else if (TString(tdm->GetTypeName())=="UShort_t")
705     typestring = "s";
706     else if (TString(tdm->GetTypeName())=="Int_t")
707     typestring = "I";
708     else if (TString(tdm->GetTypeName())=="UInt_t")
709     typestring = "i";
710     else if (TString(tdm->GetTypeName())=="Float_t")
711     typestring = "F";
712     else if (TString(tdm->GetTypeName())=="Double_t")
713     typestring = "D";
714     else if (TString(tdm->GetTypeName())=="Long64_t")
715     typestring = "L";
716     else if (TString(tdm->GetTypeName())=="ULong64_t")
717     typestring = "l";
718     else if (TString(tdm->GetTypeName())=="Bool_t")
719     typestring = "O";
720     else
721     continue;
722 bendavid 1.1 //printf("%s\n",tdm->GetTypeName());
723     TString varname = TString::Format("ph%d.%s",iph+1,tdm->GetName());
724    
725     Char_t *addr = (Char_t*)&fDiphotonEvent->photons[iph];
726     assert(sizeof(Char_t)==1);
727 paus 1.12
728     if (fWriteDiphotonTree)
729     hCiCTuple->Branch(varname,addr+tdm->GetOffset(),
730     TString::Format("%s/%s",varname.Data(),typestring.Data()));
731 bendavid 1.1
732     if (iph==0) {
733     TString singlename = TString::Format("ph.%s",tdm->GetName());
734     Char_t *addrsingle = (Char_t*)fSinglePhoton;
735 paus 1.12 if (fWriteSingleTree)
736     hCiCTupleSingle->Branch(singlename,addrsingle+tdm->GetOffset(),
737     TString::Format("%s/%s",singlename.Data(),typestring.Data()));
738 bendavid 1.1 }
739     }
740     }
741 paus 1.12
742     if (fWriteDiphotonTree)
743     AddOutput(hCiCTuple);
744     if (fWriteSingleTree)
745     AddOutput(hCiCTupleSingle);
746 bendavid 1.1 }
747    
748     // ----------------------------------------------------------------------------------------
749     // some helpfer functions....
750 paus 1.12 void PhotonTreeWriter::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass)
751     {
752     pt = -999.;
753 bendavid 1.1 decayZ = -999.;
754 paus 1.12 mass = -999.;
755 bendavid 1.1
756     // loop over all GEN particles and look for status 1 photons
757     for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
758     const MCParticle* p = fMCParticles->At(i);
759 paus 1.12 if (p->Is(MCParticle::kH) || (!fApplyElectronVeto &&
760     (p->AbsPdgId()==23 || p->AbsPdgId()==24))) {
761     pt = p->Pt();
762 bendavid 1.1 decayZ = p->DecayVertex().Z();
763 paus 1.12 mass = p->Mass();
764 bendavid 1.1 break;
765     }
766     }
767    
768     return;
769     }
770    
771    
772 paus 1.12 Float_t PhotonTreeWriter::GetEventCat(PhotonTools::CiCBaseLineCats cat1,
773     PhotonTools::CiCBaseLineCats cat2) {
774 bendavid 1.1
775 paus 1.12 bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
776     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
777 bendavid 1.1
778     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
779     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
780    
781     if( ph1IsEB && ph2IsEB )
782     return ( ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
783    
784     return ( ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
785     }
786    
787 paus 1.12 void PhotonTreeWriterPhoton::SetVars(const Photon *p, const DecayParticle *c, const Electron *ele,
788     const SuperCluster *pfsc, const MCParticle *m,
789     PhotonFix &phfixph, PhotonFix &phfixele,
790     const TrackCol* trackCol,const VertexCol* vtxCol,Double_t rho,
791     const ElectronCol* els, Bool_t applyElectronVeto) {
792    
793     const SuperCluster *s = 0;
794     if (p)
795     s = p->SCluster();
796     else
797     s = ele->SCluster();
798     const BasicCluster *b = s->Seed();
799     const BasicCluster *b2 = 0;
800     Double_t ebcmax = -99.;
801     for (UInt_t i=0; i<s->ClusterSize(); ++i) {
802     const BasicCluster *bc = s->Cluster(i);
803     if (bc->Energy() > ebcmax && bc !=b) {
804     b2 = bc;
805     ebcmax = bc->Energy();
806     }
807     }
808 bendavid 1.2
809 paus 1.12 const BasicCluster *bclast = 0;
810     Double_t ebcmin = 1e6;
811     for (UInt_t i=0; i<s->ClusterSize(); ++i) {
812     const BasicCluster *bc = s->Cluster(i);
813     if (bc->Energy() < ebcmin && bc !=b) {
814     bclast = bc;
815     ebcmin = bc->Energy();
816     }
817     }
818 bendavid 1.5
819 paus 1.12 const BasicCluster *bclast2 = 0;
820     ebcmin = 1e6;
821     for (UInt_t i=0; i<s->ClusterSize(); ++i) {
822     const BasicCluster *bc = s->Cluster(i);
823     if (bc->Energy() < ebcmin && bc !=b && bc!=bclast) {
824     bclast2 = bc;
825     ebcmin = bc->Energy();
826     }
827     }
828 bendavid 1.5
829 paus 1.12 if (p) {
830     hasphoton = kTRUE;
831     e = p->E();
832     pt = p->Pt();
833     eta = p->Eta();
834     phi = p->Phi();
835     r9 = p->R9();
836     e3x3 = p->E33();
837     e5x5 = p->E55();
838     hovere = p->HadOverEm();
839     sigietaieta = p->CoviEtaiEta();
840     phcat = PhotonTools::CiCBaseLineCat(p);
841     eerr = p->EnergyErr();
842     eerrsmeared = p->EnergyErrSmeared();
843     esmearing = p->EnergySmearing();
844     idmva = p->IdMva();
845     hcalisodr03 = p->HcalTowerSumEtDr03();
846     ecalisodr03 = p->EcalRecHitIsoDr03();
847     trkisohollowdr03 = p->HollowConeTrkIsoDr03();
848     hcalisodr04 = p->HcalTowerSumEtDr04();
849     ecalisodr04 = p->EcalRecHitIsoDr04();
850     trkisohollowdr04 = p->HollowConeTrkIsoDr04();
851    
852    
853     const Vertex *vtx = vtxCol->At(0);
854     if (p->HasPV()) vtx = p->PV();
855    
856     UInt_t wVtxInd = 0;
857    
858     trackiso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0,
859     trackCol, NULL, NULL,
860     (!applyElectronVeto ? els : NULL) );
861     //Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
862 bendavid 1.2
863 paus 1.12 // track iso worst vtx
864     trackiso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0,
865     trackCol, &wVtxInd,vtxCol,
866     (!applyElectronVeto ? els : NULL) );
867     combiso1 = ecalisodr03+hcalisodr04+trackiso1 - 0.17*rho;
868     combiso2 = ecalisodr04+hcalisodr04+trackiso2 - 0.52*rho;
869     }
870     else {
871     hasphoton = kFALSE;
872     e = -99.;
873     pt = -99.;
874     eta = -99.;
875     phi = -99.;
876     r9 = b->E3x3()/s->RawEnergy();
877     e3x3 = b->E3x3();
878     e5x5 = b->E5x5();
879     hovere = ele->HadronicOverEm();
880     sigietaieta = ele->CoviEtaiEta();
881     phcat = -99;
882     eerr = -99.;
883     eerrsmeared = -99.;
884     esmearing = 0.;
885     idmva = -99.;
886     hcalisodr03 = -99;
887     ecalisodr03 = -99;
888     trkisohollowdr03 = -99;
889     hcalisodr04 = -99;
890     ecalisodr04 = -99;
891     trkisohollowdr04 = -99;
892     trackiso1 = -99.;
893     trackiso2 = -99.;
894     combiso1 = -99.;
895     combiso2 = -99.;
896     }
897    
898     sce = s->Energy();
899     scrawe = s->RawEnergy();
900     scpse = s->PreshowerEnergy();
901     sceta = s->Eta();
902     scphi = s->Phi();
903     scnclusters = s->ClusterSize();
904     scnhits = s->NHits();
905     scetawidth = -99.;
906     scphiwidth = -99.;
907     if (p) {
908     scetawidth = p->EtaWidth();
909     scphiwidth = p->PhiWidth();
910     }
911     else {
912     scetawidth = s->EtaWidth();
913     scphiwidth = s->PhiWidth();
914     }
915     isbarrel = (s->AbsEta()<1.5);
916     isr9reco = (isbarrel && r9>0.94) || (!isbarrel && r9>0.95);
917     isr9cat = (r9>0.94);
918    
919     eseed = b->Energy();
920     etaseed = b->Eta();
921     phiseed = b->Phi();
922     ietaseed = b->IEta();
923     iphiseed = b->IPhi();
924     ixseed = b->IX();
925     iyseed = b->IY();
926     etacryseed = b->EtaCry();
927     phicryseed = b->PhiCry();
928     xcryseed = b->XCry();
929     ycryseed = b->YCry();
930     thetaaxisseed = b->ThetaAxis();
931     phiaxisseed = b->PhiAxis();
932     sigietaietaseed = TMath::Sqrt(b->CoviEtaiEta());
933     sigiphiphiseed = TMath::Sqrt(b->CoviPhiiPhi());
934     if (isnan(sigiphiphiseed)) sigiphiphiseed = -99.;
935     covietaiphiseed = b->CoviEtaiPhi();
936     if (isnan(covietaiphiseed)) covietaiphiseed = -99.;
937     e3x3seed = b->E3x3();
938     e5x5seed = b->E5x5();
939     emaxseed = b->EMax();
940     e2ndseed = b->E2nd();
941     etopseed = b->ETop();
942     ebottomseed = b->EBottom();
943     eleftseed = b->ELeft();
944     erightseed = b->ERight();
945     e1x3seed = b->E1x3();
946     e3x1seed = b->E3x1();
947     e1x5seed = b->E1x5();
948     e2x2seed = b->E2x2();
949     e4x4seed = b->E4x4();
950     e2x5maxseed = b->E2x5Max();
951     e2x5topseed = b->E2x5Top();
952     e2x5bottomseed = b->E2x5Bottom();
953     e2x5leftseed = b->E2x5Left();
954     e2x5rightseed = b->E2x5Right();
955     xseedseed = b->Pos().X();
956     yseedseed = b->Pos().Y();
957     zseedseed = b->Pos().Z();
958     nhitsseed = b->NHits();
959    
960     if (b2) {
961     ebc2 = b2->Energy();
962     etabc2 = b2->Eta();
963     phibc2 = b2->Phi();
964     ietabc2 = b2->IEta();
965     iphibc2 = b2->IPhi();
966     ixbc2 = b2->IX();
967     iybc2 = b2->IY();
968     etacrybc2 = b2->EtaCry();
969     phicrybc2 = b2->PhiCry();
970     xcrybc2 = b2->XCry();
971     ycrybc2 = b2->YCry();
972     thetaaxisbc2 = b2->ThetaAxis();
973     phiaxisbc2 = b2->PhiAxis();
974     sigietaietabc2 = TMath::Sqrt(b2->CoviEtaiEta());
975     sigiphiphibc2 = TMath::Sqrt(b2->CoviPhiiPhi());
976     if (isnan(sigiphiphibc2)) sigiphiphibc2 = -99.;
977     covietaiphibc2 = b2->CoviEtaiPhi();
978     if (isnan(covietaiphibc2)) covietaiphibc2 = -99.;
979     e3x3bc2 = b2->E3x3();
980     e5x5bc2 = b2->E5x5();
981     emaxbc2 = b2->EMax();
982     e2ndbc2 = b2->E2nd();
983     etopbc2 = b2->ETop();
984     ebottombc2 = b2->EBottom();
985     eleftbc2 = b2->ELeft();
986     erightbc2 = b2->ERight();
987     e1x3bc2 = b2->E1x3();
988     e3x1bc2 = b2->E3x1();
989     e1x5bc2 = b2->E1x5();
990     e2x2bc2 = b2->E2x2();
991     e4x4bc2 = b2->E4x4();
992     e2x5maxbc2 = b2->E2x5Max();
993     e2x5topbc2 = b2->E2x5Top();
994     e2x5bottombc2 = b2->E2x5Bottom();
995     e2x5leftbc2 = b2->E2x5Left();
996     e2x5rightbc2 = b2->E2x5Right();
997     xbc2bc2 = b2->Pos().X();
998     ybc2bc2 = b2->Pos().Y();
999     zbc2bc2 = b2->Pos().Z();
1000     nhitsbc2= b2->NHits();
1001     }
1002     else {
1003     ebc2 = 0;
1004     etabc2 = 0;
1005     phibc2 = 0;
1006     ietabc2 = 0;
1007     iphibc2 = 0;
1008     ixbc2 = 0;
1009     iybc2 = 0;
1010     etacrybc2 = 0;
1011     phicrybc2 = 0;
1012     xcrybc2 = 0;
1013     ycrybc2 = 0;
1014     thetaaxisbc2 = 0;
1015     phiaxisbc2 = 0;
1016     sigietaietabc2 = 0;
1017     sigiphiphibc2 = 0;
1018     covietaiphibc2 = 0;
1019     e3x3bc2 = 0;
1020     e5x5bc2 = 0;
1021     emaxbc2 = 0;
1022     e2ndbc2 = 0;
1023     etopbc2 = 0;
1024     ebottombc2 = 0;
1025     eleftbc2 = 0;
1026     erightbc2 = 0;
1027     e1x3bc2 = 0;
1028     e3x1bc2 = 0;
1029     e1x5bc2 = 0;
1030     e2x2bc2 = 0;
1031     e4x4bc2 = 0;
1032     e2x5maxbc2 = 0;
1033     e2x5topbc2 = 0;
1034     e2x5bottombc2 = 0;
1035     e2x5leftbc2 = 0;
1036     e2x5rightbc2 = 0;
1037     xbc2bc2 = 0;
1038     ybc2bc2 = 0;
1039     zbc2bc2 = 0;
1040     nhitsbc2 = 0;
1041     }
1042 bendavid 1.5
1043 paus 1.12 if (bclast) {
1044     ebclast = bclast->Energy();
1045     etabclast = bclast->Eta();
1046     phibclast = bclast->Phi();
1047     ietabclast = bclast->IEta();
1048     iphibclast = bclast->IPhi();
1049     ixbclast = bclast->IX();
1050     iybclast = bclast->IY();
1051     etacrybclast = bclast->EtaCry();
1052     phicrybclast = bclast->PhiCry();
1053     xcrybclast = bclast->XCry();
1054     ycrybclast = bclast->YCry();
1055     thetaaxisbclast = bclast->ThetaAxis();
1056     phiaxisbclast = bclast->PhiAxis();
1057     sigietaietabclast = TMath::Sqrt(bclast->CoviEtaiEta());
1058     sigiphiphibclast = TMath::Sqrt(bclast->CoviPhiiPhi());
1059     if (isnan(sigiphiphibclast)) sigiphiphibclast = -99.;
1060     covietaiphibclast = bclast->CoviEtaiPhi();
1061     if (isnan(covietaiphibclast)) covietaiphibclast = -99.;
1062     e3x3bclast = bclast->E3x3();
1063     nhitsbclast = bclast->NHits();
1064     }
1065     else {
1066     ebclast = 0;
1067     etabclast = 0;
1068     phibclast = 0;
1069     ietabclast = 0;
1070     iphibclast = 0;
1071     ixbclast = 0;
1072     iybclast = 0;
1073     etacrybclast = 0;
1074     phicrybclast = 0;
1075     xcrybclast = 0;
1076     ycrybclast = 0;
1077     thetaaxisbclast = 0;
1078     phiaxisbclast = 0;
1079     sigietaietabclast = 0;
1080     sigiphiphibclast = 0;
1081     covietaiphibclast = 0;
1082     e3x3bclast = 0;
1083     nhitsbclast = 0;
1084     }
1085 bendavid 1.5
1086 paus 1.12 if (bclast2) {
1087     ebclast2 = bclast2->Energy();
1088     etabclast2 = bclast2->Eta();
1089     phibclast2 = bclast2->Phi();
1090     ietabclast2 = bclast2->IEta();
1091     iphibclast2 = bclast2->IPhi();
1092     ixbclast2 = bclast2->IX();
1093     iybclast2 = bclast2->IY();
1094     etacrybclast2 = bclast2->EtaCry();
1095     phicrybclast2 = bclast2->PhiCry();
1096     xcrybclast2 = bclast2->XCry();
1097     ycrybclast2 = bclast2->YCry();
1098     thetaaxisbclast2 = bclast2->ThetaAxis();
1099     phiaxisbclast2 = bclast2->PhiAxis();
1100     sigietaietabclast2 = TMath::Sqrt(bclast2->CoviEtaiEta());
1101     sigiphiphibclast2 = TMath::Sqrt(bclast2->CoviPhiiPhi());
1102     if (isnan(sigiphiphibclast2)) sigiphiphibclast2 = -99.;
1103     covietaiphibclast2 = bclast2->CoviEtaiPhi();
1104     if (isnan(covietaiphibclast2)) covietaiphibclast2 = -99.;
1105     e3x3bclast2 = bclast2->E3x3();
1106     nhitsbclast2 = bclast2->NHits();
1107     }
1108     else {
1109     ebclast2 = 0;
1110     etabclast2 = 0;
1111     phibclast2 = 0;
1112     ietabclast2 = 0;
1113     iphibclast2 = 0;
1114     ixbclast2 = 0;
1115     iybclast2 = 0;
1116     etacrybclast2 = 0;
1117     phicrybclast2 = 0;
1118     xcrybclast2 = 0;
1119     ycrybclast2 = 0;
1120     thetaaxisbclast2 = 0;
1121     phiaxisbclast2 = 0;
1122     sigietaietabclast2 = 0;
1123     sigiphiphibclast2 = 0;
1124     covietaiphibclast2 = 0;
1125     e3x3bclast2 = 0;
1126     nhitsbclast2 = 0;
1127     }
1128 bendavid 1.5
1129 paus 1.12 //initialize photon energy corrections if needed
1130     /*if (!PhotonFix::initialised()) {
1131     PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
1132     }*/
1133    
1134     phfixph.setup(e,sceta,scphi,r9);
1135     phfixele.setup(e,sceta,scphi,r9);
1136    
1137     const Float_t dval = -99.;
1138     ecor = phfixph.fixedEnergy();
1139     ecorerr = phfixph.sigmaEnergy();
1140     ecorele = phfixele.fixedEnergy();
1141     ecoreleerr = phfixele.sigmaEnergy();
1142     if (phfixph.isbarrel()) {
1143     etac = phfixph.etaC();
1144     etas = phfixph.etaS();
1145     etam = phfixph.etaM();
1146     phic = phfixph.phiC();
1147     phis = phfixph.phiS();
1148     phim = phfixph.phiM();
1149     xz = dval;
1150     xc = dval;
1151     xs = dval;
1152     xm = dval;
1153     yz = dval;
1154     yc = dval;
1155     ys = dval;
1156     ym = dval;
1157     }
1158     else {
1159     etac = dval;
1160     etas = dval;
1161     etam = dval;
1162     phic = dval;
1163     phis = dval;
1164     phim = dval;
1165     xz = phfixph.xZ();
1166     xc = phfixph.xC();
1167     xs = phfixph.xS();
1168     xm = phfixph.xM();
1169     yz = phfixph.yZ();
1170     yc = phfixph.yC();
1171     ys = phfixph.yS();
1172     ym = phfixph.yM();
1173     }
1174    
1175     if (c) {
1176     hasconversion = kTRUE;
1177     convp = c->P();
1178     convpt = c->Pt();
1179     conveta = c->Eta();
1180     convphi = c->Phi();
1181     ThreeVector dirconvsc = ThreeVector(s->Point()) - c->Position();
1182     convdeta = c->Eta() - dirconvsc.Eta();
1183     convdphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
1184     convvtxrho = c->Position().Rho();
1185     convvtxz = c->Position().Z();
1186     convvtxphi = c->Position().Phi();
1187    
1188     const StableData *leadsd = dynamic_cast<const StableData*>(c->DaughterDat(0));
1189     const StableData *trailsd = dynamic_cast<const StableData*>(c->DaughterDat(1));
1190     if (leadsd->Pt()<trailsd->Pt()) {
1191     const StableData *sdtmp = leadsd;
1192     leadsd = trailsd;
1193     trailsd = sdtmp;
1194     }
1195    
1196     const Track *leadtrack = dynamic_cast<const StableParticle*>(leadsd->Original())->Trk();
1197     const Track *trailtrack = dynamic_cast<const StableParticle*>(trailsd->Original())->Trk();
1198    
1199     convleadpt = leadsd->Pt();
1200     convtrailpt = trailsd->Pt();
1201     convleadtrackpt = leadtrack->Pt();
1202     convleadtrackalgo = leadtrack->Algo();
1203     if (convleadtrackalgo==29) convleadtrackalgos=2; //gsf track
1204     else if (convleadtrackalgo==15 ||convleadtrackalgo==16) convleadtrackalgos=1; //ecal-seeded track
1205     else convleadtrackalgos = 0; //std iterative track
1206     convleadtrackcharge = leadtrack->Charge();
1207     convtrailtrackpt = trailtrack->Pt();
1208     convtrailtrackalgo = trailtrack->Algo();
1209     if (convtrailtrackalgo==29) convtrailtrackalgos=2; //gsf track
1210     else if (convtrailtrackalgo==15 ||convtrailtrackalgo==16) convtrailtrackalgos=1; //ecal-seeded track
1211     else convtrailtrackalgos = 0; //std iterative track
1212     trailtrackcharge = trailtrack->Charge();
1213     }
1214     else {
1215     hasconversion = kFALSE;
1216     convp = -99.;
1217     convpt = -99.;
1218     conveta = -99.;
1219     convphi = -99.;
1220     convdeta = -99.;
1221     convdphi = -99.;
1222     convvtxrho = -99.;
1223     convvtxz = -999.;
1224     convvtxphi = -99.;
1225     convleadpt = -99.;
1226     convtrailpt = -99.;
1227     convleadtrackpt = -99.;
1228     convleadtrackalgo = -99;
1229     convleadtrackalgos = -99;
1230     convleadtrackcharge = 0;
1231     convtrailtrackpt = -99.;
1232     convtrailtrackalgo = -99;
1233     convtrailtrackalgos = -99;
1234     trailtrackcharge = 0;
1235     }
1236    
1237     //electron quantities
1238     if (ele) {
1239     haselectron = kTRUE;
1240     eleisecaldriven = ele->IsEcalDriven();
1241     eleistrackerdriven = ele->IsTrackerDriven();
1242     elee = ele->E();
1243     elept = ele->Pt();
1244     eleeta = ele->Eta();
1245     elephi = ele->Phi();
1246     elecharge = ele->Charge();
1247     elefbrem = ele->FBrem();
1248     eledeta = ele->DeltaEtaSuperClusterTrackAtVtx();
1249     eledphi = ele->DeltaPhiSuperClusterTrackAtVtx();
1250     elep = s->Energy()/ele->ESuperClusterOverP();
1251     elepin = ele->PIn();
1252     elepout = ele->POut();
1253     }
1254     else {
1255     haselectron = kFALSE;
1256     eleisecaldriven = kFALSE;
1257     eleistrackerdriven = kFALSE;
1258     elee = -99.;
1259     elept = -99.;
1260     eleeta = -99.;
1261     elephi = -99.;
1262     elecharge = -99;
1263     elefbrem = -99.;
1264     eledeta = -99.;
1265     eledphi = -99.;
1266     elep = -99.;
1267     elepin = -99.;
1268     elepout = -99.;
1269     }
1270    
1271     //pf supercluster quantities
1272     if (pfsc) {
1273     haspfsc = kTRUE;
1274     pfsce = pfsc->Energy();
1275     pfscrawe = pfsc->RawEnergy();
1276     pfsceta = pfsc->Eta();
1277     pfscphi = pfsc->Phi();
1278     }
1279     else {
1280     haspfsc = kFALSE;
1281     pfsce = -99.;
1282     pfscrawe = -99.;
1283     pfsceta = -99.;
1284     pfscphi = -99.;
1285     }
1286    
1287     genz = -99.;
1288     if (m) {
1289     ispromptgen = kTRUE;
1290     gene = m->E();
1291     genpt = m->Pt();
1292     geneta = m->Eta();
1293     genphi = m->Phi();
1294     const MCParticle *mm = m->DistinctMother();
1295     if (mm) genz = mm->DecayVertex().Z();
1296     pdgid = m->PdgId();
1297     if (mm) motherpdgid = mm->PdgId();
1298     else motherpdgid = -99;
1299     }
1300     else {
1301     ispromptgen = kFALSE;
1302     gene = -99.;
1303     genpt = -99.;
1304     geneta = -99.;
1305     genphi = -99.;
1306     pdgid = -99;
1307     motherpdgid = -99;
1308     }
1309 bendavid 1.1 }