ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.5
Committed: Sun Dec 11 00:03:04 2011 UTC (13 years, 5 months ago) by bendavid
Content type: text/plain
Branch: MAIN
Changes since 1.4: +320 -47 lines
Log Message:
more photon updates for mva plus id analysis

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