ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.4
Committed: Fri Nov 18 00:07:16 2011 UTC (13 years, 5 months ago) by bendavid
Content type: text/plain
Branch: MAIN
Changes since 1.3: +77 -6 lines
Log Message:
8 cat energy smearing and scaling, bdt vertex selection and probability, scaled pt cuts for photons

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.1
43    
44     fIsData (false),
45     fPhotonsFromBranch (true),
46     fPVFromBranch (true),
47     fGoodElectronsFromBranch (kTRUE),
48    
49     // ----------------------------------------
50     // collections....
51     fPhotons (0),
52     fElectrons (0),
53     fConversions (0),
54     fTracks (0),
55     fPileUpDen (0),
56     fPV (0),
57     fBeamspot (0),
58     fPFCands (0),
59     fMCParticles (0),
60     fPileUp (0),
61 bendavid 1.2 fSuperClusters (0),
62 bendavid 1.1
63 bendavid 1.2 fLoopOnGoodElectrons(kFALSE),
64 bendavid 1.1 fInvertElectronVeto(kFALSE),
65    
66     fWriteDiphotonTree(kTRUE),
67     fWriteSingleTree(kTRUE),
68 bendavid 1.2 fExcludeSinglePrompt(kFALSE),
69     fExcludeDoublePrompt(kFALSE),
70     fPhFixDataFile(gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat")),
71 bendavid 1.1 fTupleName ("hPhotonTree")
72    
73    
74     {
75     // Constructor.
76     }
77    
78     PhotonTreeWriter::~PhotonTreeWriter(){
79     ;
80     }
81    
82     //--------------------------------------------------------------------------------------------------
83     void PhotonTreeWriter::Process()
84     {
85     // ------------------------------------------------------------
86     // Process entries of the tree.
87 bendavid 1.2
88 bendavid 1.1 LoadEventObject(fPhotonBranchName, fPhotons);
89 bendavid 1.2 LoadEventObject(fGoodElectronName, fGoodElectrons);
90    
91     const BaseCollection *egcol = 0;
92     if (fLoopOnGoodElectrons) {
93     egcol = fGoodElectrons;
94     }
95     else {
96     egcol = fPhotons;
97     }
98 bendavid 1.1
99 bendavid 1.2 if (egcol->GetEntries()<1) return;
100 bendavid 1.1
101     LoadEventObject(fElectronName, fElectrons);
102     LoadEventObject(fConversionName, fConversions);
103     LoadEventObject(fTrackBranchName, fTracks);
104     LoadEventObject(fPileUpDenName, fPileUpDen);
105     LoadEventObject(fPVName, fPV);
106     LoadEventObject(fBeamspotName, fBeamspot);
107     LoadEventObject(fPFCandName, fPFCands);
108 bendavid 1.2 LoadEventObject(fSuperClusterName, fSuperClusters);
109     LoadEventObject(fPFMetName, fPFMet);
110 bendavid 1.1
111     // ------------------------------------------------------------
112     // load event based information
113     Int_t _numPU = -1.; // some sensible default values....
114     Int_t _numPUminus = -1.; // some sensible default values....
115     Int_t _numPUplus = -1.; // some sensible default values....
116    
117     Float_t _tRho = -99.;
118     if( fPileUpDen->GetEntries() > 0 )
119     _tRho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
120    
121     const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));
122    
123     if( !fIsData ) {
124     LoadBranch(fMCParticleName);
125     LoadBranch(fPileUpName);
126     }
127    
128     if( !fIsData ) {
129     for (UInt_t i=0; i<fPileUp->GetEntries(); ++i) {
130     const PileupInfo *puinfo = fPileUp->At(i);
131     if (puinfo->GetBunchCrossing()==0) _numPU = puinfo->GetPU_NumInteractions();
132     else if (puinfo->GetBunchCrossing() == -1) _numPUminus = puinfo->GetPU_NumInteractions();
133     else if (puinfo->GetBunchCrossing() == 1) _numPUplus = puinfo->GetPU_NumInteractions();
134     }
135     }
136    
137     // in case of a MC event, try to find Higgs and Higgs decay Z poisition
138     Float_t _pth = -100.;
139     Float_t _decayZ = -100.;
140     Float_t _genmass = -100.;
141     if( !fIsData ) FindHiggsPtAndZ(_pth, _decayZ, _genmass);
142    
143     fDiphotonEvent->rho = _tRho;
144     fDiphotonEvent->genHiggspt = _pth;
145     fDiphotonEvent->genHiggsZ = _decayZ;
146     fDiphotonEvent->genmass = _genmass;
147     fDiphotonEvent->gencostheta = -99.;
148     fDiphotonEvent->nVtx = fPV->GetEntries();
149 bendavid 1.2 fDiphotonEvent->bsX = fBeamspot->At(0)->X();
150     fDiphotonEvent->bsY = fBeamspot->At(0)->Y();
151     fDiphotonEvent->bsZ = fBeamspot->At(0)->Z();
152     fDiphotonEvent->vtxX = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->X() : -99.;
153     fDiphotonEvent->vtxY = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Y() : -99.;
154 bendavid 1.1 fDiphotonEvent->vtxZ = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Z() : -99.;
155     fDiphotonEvent->numPU = _numPU;
156     fDiphotonEvent->numPUminus = _numPUminus;
157     fDiphotonEvent->numPUplus = _numPUplus;
158     fDiphotonEvent->mass = -99.;
159     fDiphotonEvent->ptgg = -99.;
160     fDiphotonEvent->costheta = -99.;
161 bendavid 1.2 fDiphotonEvent->mt = -99.;
162     fDiphotonEvent->cosphimet = -99.;
163     fDiphotonEvent->mtele = -99.;
164     fDiphotonEvent->cosphimetele = -99.;
165 bendavid 1.1 fDiphotonEvent->evt = GetEventHeader()->EvtNum();
166     fDiphotonEvent->run = GetEventHeader()->RunNum();
167     fDiphotonEvent->lumi = GetEventHeader()->LumiSec();
168     fDiphotonEvent->evtcat = -99;
169 bendavid 1.2 fDiphotonEvent->nobj = fPhotons->GetEntries();
170     fDiphotonEvent->pfmet = fPFMet->At(0)->Pt();
171     fDiphotonEvent->pfmetphi = fPFMet->At(0)->Phi();
172     fDiphotonEvent->pfmetx = fPFMet->At(0)->Px();
173     fDiphotonEvent->pfmety = fPFMet->At(0)->Py();
174 bendavid 1.1 fDiphotonEvent->masscor = -99.;
175     fDiphotonEvent->masscorerr = -99.;
176 bendavid 1.2 fDiphotonEvent->masscorele = -99.;
177     fDiphotonEvent->masscoreleerr = -99.;
178     fDiphotonEvent->ismc = GetEventHeader()->IsMC();
179    
180 bendavid 1.1 Int_t nhitsbeforevtxmax = 1;
181     if (fInvertElectronVeto) nhitsbeforevtxmax = 999;
182    
183 bendavid 1.2 if (egcol->GetEntries()>=2) {
184 bendavid 1.1
185 bendavid 1.2 const Particle *p1 = 0;
186     const Particle *p2 = 0;
187     const Photon *phHard = 0;
188     const Photon *phSoft = 0;
189     const Electron *ele1 = 0;
190     const Electron *ele2 = 0;
191     const SuperCluster *sc1 = 0;
192     const SuperCluster *sc2 = 0;
193    
194     if (fLoopOnGoodElectrons) {
195     ele1 = fGoodElectrons->At(0);
196     ele2 = fGoodElectrons->At(1);
197     p1 = ele1;
198     p2 = ele2;
199     sc1 = ele1->SCluster();
200     sc2 = ele2->SCluster();
201     phHard = PhotonTools::MatchedPhoton(ele1,fPhotons);
202     phSoft = PhotonTools::MatchedPhoton(ele2,fPhotons);
203     }
204     else {
205     phHard = fPhotons->At(0);
206     phSoft = fPhotons->At(1);
207     p1 = phHard;
208     p2 = phSoft;
209     sc1 = phHard->SCluster();
210     sc2 = phSoft->SCluster();
211     ele1 = PhotonTools::MatchedElectron(phHard,fGoodElectrons);
212     ele2 = PhotonTools::MatchedElectron(phSoft,fGoodElectrons);
213 bendavid 1.1 }
214    
215 bendavid 1.2 const DecayParticle *conv1 = PhotonTools::MatchedConversion(sc1,fConversions,bsp,nhitsbeforevtxmax);
216     const DecayParticle *conv2 = PhotonTools::MatchedConversion(sc2,fConversions,bsp,nhitsbeforevtxmax);
217    
218     const SuperCluster *pfsc1 = PhotonTools::MatchedSC(sc1,fSuperClusters);
219     const SuperCluster *pfsc2 = PhotonTools::MatchedSC(sc2,fSuperClusters);
220 bendavid 1.1
221     const MCParticle *phgen1 = 0;
222     const MCParticle *phgen2 = 0;
223     if( !fIsData ) {
224 bendavid 1.2 phgen1 = PhotonTools::MatchMC(p1,fMCParticles,fInvertElectronVeto);
225     phgen2 = PhotonTools::MatchMC(p2,fMCParticles,fInvertElectronVeto);
226     }
227    
228     if (fExcludeSinglePrompt && (phgen1 || phgen2) ) return;
229     if (fExcludeDoublePrompt && (phgen1 && phgen2) ) return;
230    
231     if (!fLoopOnGoodElectrons && phHard->HasPV()) {
232     fDiphotonEvent->vtxX = phHard->PV()->X();
233     fDiphotonEvent->vtxY = phHard->PV()->Y();
234     fDiphotonEvent->vtxZ = phHard->PV()->Z();
235 bendavid 1.4 fDiphotonEvent->vtxprob = phHard->VtxProb();
236 bendavid 1.2 }
237    
238 bendavid 1.4 Double_t _mass = -99.;
239     Double_t _masserr = -99.;
240     Double_t _masserrsmeared = -99.;
241     Double_t _masserrwrongvtx = -99.;
242     Double_t _masserrsmearedwrongvtx = -99.;
243     Double_t _ptgg = -99.;
244     Double_t _costheta = -99.;
245 bendavid 1.2 PhotonTools::DiphotonR9EtaPtCats _evtcat = PhotonTools::kOctCat0;
246     if (phHard && phSoft) {
247     _mass = (phHard->Mom()+phSoft->Mom()).M();
248     _masserr = 0.5*_mass*TMath::Sqrt(phHard->EnergyErr()*phHard->EnergyErr()/phHard->E()/phHard->E() + phSoft->EnergyErr()*phSoft->EnergyErr()/phSoft->E()/phSoft->E());
249     _masserrsmeared = 0.5*_mass*TMath::Sqrt(phHard->EnergyErrSmeared()*phHard->EnergyErrSmeared()/phHard->E()/phHard->E() + phSoft->EnergyErrSmeared()*phSoft->EnergyErrSmeared()/phSoft->E()/phSoft->E());
250     _ptgg = (phHard->Mom()+phSoft->Mom()).Pt();
251     _costheta = ThreeVector(phHard->Mom()).Unit().Dot(ThreeVector(phSoft->Mom()).Unit());
252     _evtcat = PhotonTools::DiphotonR9EtaPtCat(phHard,phSoft);
253 bendavid 1.4
254     const Double_t dz = 5.8;
255     Double_t deltamvtx = _mass*VertexTools::DeltaMassVtx(phHard->CaloPos().X(), phHard->CaloPos().Y(), phHard->CaloPos().Z(),
256     phSoft->CaloPos().X(), phSoft->CaloPos().Y(), phSoft->CaloPos().Z(),
257     dz);
258    
259     _masserrwrongvtx = TMath::Sqrt(_masserr*_masserr + deltamvtx*deltamvtx);
260     _masserrsmearedwrongvtx = TMath::Sqrt(_masserrsmeared*_masserrsmeared + deltamvtx*deltamvtx);
261    
262    
263 bendavid 1.2 }
264    
265    
266     Float_t _massele = -99.;
267     Float_t _ptee = -99.;
268     Float_t _costhetaele = -99.;
269     if (ele1 && ele2) {
270     _massele = (ele1->Mom()+ele2->Mom()).M();
271     _ptee = (ele1->Mom()+ele2->Mom()).Pt();
272     _costhetaele = ThreeVector(ele1->Mom()).Unit().Dot(ThreeVector(ele2->Mom()).Unit());
273 bendavid 1.1 }
274    
275     Float_t _gencostheta = -99.;
276     if (phgen1 && phgen2) {
277     _gencostheta = ThreeVector(phgen1->Mom()).Unit().Dot(ThreeVector(phgen2->Mom()).Unit());
278     }
279    
280     fDiphotonEvent->gencostheta = _gencostheta;
281     fDiphotonEvent->mass = _mass;
282 bendavid 1.2 fDiphotonEvent->masserr = _masserr;
283     fDiphotonEvent->masserrsmeared = _masserrsmeared;
284 bendavid 1.4 fDiphotonEvent->masserrwrongvtx = _masserrwrongvtx;
285     fDiphotonEvent->masserrsmearedwrongvtx = _masserrsmearedwrongvtx;
286 bendavid 1.1 fDiphotonEvent->ptgg = _ptgg;
287 bendavid 1.2 fDiphotonEvent->costheta = _costheta;;
288     fDiphotonEvent->massele = _massele;
289     fDiphotonEvent->ptee = _ptee;
290     fDiphotonEvent->costhetaele = _costhetaele;
291     fDiphotonEvent->evtcat = _evtcat;
292 bendavid 1.1
293 bendavid 1.2 fDiphotonEvent->photons[0].SetVars(phHard,conv1,ele1,pfsc1,phgen1,fPhfixph,fPhfixele);
294     fDiphotonEvent->photons[1].SetVars(phSoft,conv2,ele2,pfsc2,phgen2,fPhfixph,fPhfixele);
295 bendavid 1.1
296     Float_t ph1ecor = fDiphotonEvent->photons[0].Ecor();
297     Float_t ph1ecorerr = fDiphotonEvent->photons[0].Ecorerr();
298     Float_t ph2ecor = fDiphotonEvent->photons[1].Ecor();
299     Float_t ph2ecorerr = fDiphotonEvent->photons[1].Ecorerr();
300 bendavid 1.2
301     Float_t ph1ecorele = fDiphotonEvent->photons[0].Ecorele();
302     Float_t ph1ecoreleerr = fDiphotonEvent->photons[0].Ecoreleerr();
303     Float_t ph2ecorele = fDiphotonEvent->photons[1].Ecorele();
304     Float_t ph2ecoreleerr = fDiphotonEvent->photons[1].Ecoreleerr();
305    
306 bendavid 1.1
307     fDiphotonEvent->masscor = TMath::Sqrt(2.0*ph1ecor*ph2ecor*(1.0-fDiphotonEvent->costheta));
308     fDiphotonEvent->masscorerr = 0.5*fDiphotonEvent->masscor*TMath::Sqrt(ph1ecorerr*ph1ecorerr/ph1ecor/ph1ecor + ph2ecorerr*ph2ecorerr/ph2ecor/ph2ecor);
309    
310 bendavid 1.2 fDiphotonEvent->masscorele = TMath::Sqrt(2.0*ph1ecorele*ph2ecorele*(1.0-fDiphotonEvent->costheta));
311     fDiphotonEvent->masscoreleerr = 0.5*fDiphotonEvent->masscorele*TMath::Sqrt(ph1ecoreleerr*ph1ecoreleerr/ph1ecorele/ph1ecorele + ph2ecoreleerr*ph2ecoreleerr/ph2ecorele/ph2ecorele);
312    
313 bendavid 1.1 //printf("r9 = %5f, photon sigieie = %5f, seed sigieie = %5f\n",phHard->R9(),phHard->CoviEtaiEta(),sqrt(phHard->SCluster()->Seed()->CoviEtaiEta()));
314    
315     if (fWriteDiphotonTree) hCiCTuple->Fill();
316    
317     }
318    
319     if (!fWriteSingleTree) return;
320    
321    
322 bendavid 1.2 for (UInt_t iph = 0; iph<egcol->GetEntries(); ++iph) {
323 bendavid 1.1
324 bendavid 1.2 const Particle *p = 0;
325     const Photon *ph = 0;
326     const Electron *ele = 0;
327     const SuperCluster *sc = 0;
328     if (fLoopOnGoodElectrons) {
329     ele = fGoodElectrons->At(iph);
330     p = ele;
331     sc = ele->SCluster();
332     ph = PhotonTools::MatchedPhoton(ele,fPhotons);
333     }
334     else {
335     ph = fPhotons->At(iph);
336     p = ph;
337     sc = ph->SCluster();
338     ele = PhotonTools::MatchedElectron(ph,fGoodElectrons);
339     }
340    
341     const DecayParticle *conv = PhotonTools::MatchedConversion(sc,fConversions,bsp,nhitsbeforevtxmax);
342     const SuperCluster *pfsc = PhotonTools::MatchedSC(sc,fSuperClusters);
343    
344    
345    
346     if (!fLoopOnGoodElectrons && ph->HasPV()) {
347 bendavid 1.1 fDiphotonEvent->vtxZ = ph->PV()->Z();
348     }
349    
350     const MCParticle *phgen = 0;
351     if( !fIsData ) {
352 bendavid 1.2 phgen = PhotonTools::MatchMC(p,fMCParticles,fInvertElectronVeto);
353 bendavid 1.1 }
354    
355 bendavid 1.2 if (fExcludeSinglePrompt && phgen) return;
356    
357     fDiphotonEvent->mt = -99.;
358     fDiphotonEvent->cosphimet = -99.;
359     fDiphotonEvent->mtele = -99.;
360     fDiphotonEvent->cosphimetele = -99.;
361    
362     if (ph) {
363     fDiphotonEvent->cosphimet = TMath::Cos(ph->Phi()-fPFMet->At(0)->Phi());
364     fDiphotonEvent->mt = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ph->Pt()*(1.0-fDiphotonEvent->cosphimet));
365     }
366 bendavid 1.1
367 bendavid 1.2 if (ele) {
368     fDiphotonEvent->cosphimetele = TMath::Cos(ele->Phi()-fPFMet->At(0)->Phi());
369     fDiphotonEvent->mtele = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ele->Pt()*(1.0-fDiphotonEvent->cosphimetele));
370     }
371    
372     fSinglePhoton->SetVars(ph,conv,ele,pfsc,phgen,fPhfixph,fPhfixele);
373 bendavid 1.1 hCiCTupleSingle->Fill();
374    
375     }
376    
377    
378     return;
379    
380     }
381    
382     //--------------------------------------------------------------------------------------------------
383     void PhotonTreeWriter::SlaveBegin()
384     {
385     // Run startup code on the computer (slave) doing the actual analysis. Here,
386     // we just request the photon collection branch.
387    
388     ReqEventObject(fPhotonBranchName, fPhotons, fPhotonsFromBranch);
389     ReqEventObject(fTrackBranchName, fTracks, true);
390     ReqEventObject(fElectronName, fElectrons, true);
391     ReqEventObject(fGoodElectronName, fGoodElectrons, fGoodElectronsFromBranch);
392     ReqEventObject(fPileUpDenName, fPileUpDen, true);
393     ReqEventObject(fPVName, fPV, fPVFromBranch);
394     ReqEventObject(fConversionName, fConversions,true);
395     ReqEventObject(fBeamspotName, fBeamspot, true);
396     ReqEventObject(fPFCandName, fPFCands, true);
397 bendavid 1.2 ReqEventObject(fSuperClusterName, fSuperClusters, true);
398     ReqEventObject(fPFMetName, fPFMet, true);
399 bendavid 1.1
400     if (!fIsData) {
401     ReqBranch(fPileUpName, fPileUp);
402     ReqBranch(fMCParticleName, fMCParticles);
403     }
404    
405    
406    
407     //initialize photon energy corrections
408     //PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
409    
410 bendavid 1.2 fPhfixph.initialise("4_2",std::string(fPhFixDataFile));
411     fPhfixele.initialise("4_2e",std::string(fPhFixDataFile));
412    
413 bendavid 1.1 fDiphotonEvent = new PhotonTreeWriterDiphotonEvent;
414     fSinglePhoton = new PhotonTreeWriterPhoton;
415    
416    
417     if (fWriteDiphotonTree) hCiCTuple = new TTree(fTupleName.Data(),fTupleName.Data());
418     TString singlename = fTupleName + TString("Single");
419     if (fWriteSingleTree) hCiCTupleSingle = new TTree(singlename,singlename);
420    
421     //make flattish tree from classes so we don't have to rely on dictionaries for reading later
422     TClass *eclass = TClass::GetClass("mithep::PhotonTreeWriterDiphotonEvent");
423     TClass *pclass = TClass::GetClass("mithep::PhotonTreeWriterPhoton");
424     TList *elist = eclass->GetListOfDataMembers();
425     TList *plist = pclass->GetListOfDataMembers();
426    
427     for (int i=0; i<elist->GetEntries(); ++i) {
428     const TDataMember *tdm = static_cast<const TDataMember*>(elist->At(i));
429     if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
430     TString typestring;
431     if (TString(tdm->GetTypeName())=="Char_t") typestring = "B";
432     else if (TString(tdm->GetTypeName())=="UChar_t") typestring = "b";
433     else if (TString(tdm->GetTypeName())=="Short_t") typestring = "S";
434     else if (TString(tdm->GetTypeName())=="UShort_t") typestring = "s";
435     else if (TString(tdm->GetTypeName())=="Int_t") typestring = "I";
436     else if (TString(tdm->GetTypeName())=="UInt_t") typestring = "i";
437     else if (TString(tdm->GetTypeName())=="Float_t") typestring = "F";
438     else if (TString(tdm->GetTypeName())=="Double_t") typestring = "D";
439     else if (TString(tdm->GetTypeName())=="Long64_t") typestring = "L";
440     else if (TString(tdm->GetTypeName())=="ULong64_t") typestring = "l";
441     else if (TString(tdm->GetTypeName())=="Bool_t") typestring = "O";
442     else continue;
443     //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
444     Char_t *addr = (Char_t*)fDiphotonEvent;
445     assert(sizeof(Char_t)==1);
446     if (fWriteDiphotonTree) hCiCTuple->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
447     if (fWriteSingleTree) hCiCTupleSingle->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
448     }
449    
450     for (int iph=0; iph<2; ++iph) {
451     for (int i=0; i<plist->GetEntries(); ++i) {
452     const TDataMember *tdm = static_cast<const TDataMember*>(plist->At(i));
453     if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
454     TString typestring;
455     if (TString(tdm->GetTypeName())=="Char_t") typestring = "B";
456     else if (TString(tdm->GetTypeName())=="UChar_t") typestring = "b";
457     else if (TString(tdm->GetTypeName())=="Short_t") typestring = "S";
458     else if (TString(tdm->GetTypeName())=="UShort_t") typestring = "s";
459     else if (TString(tdm->GetTypeName())=="Int_t") typestring = "I";
460     else if (TString(tdm->GetTypeName())=="UInt_t") typestring = "i";
461     else if (TString(tdm->GetTypeName())=="Float_t") typestring = "F";
462     else if (TString(tdm->GetTypeName())=="Double_t") typestring = "D";
463     else if (TString(tdm->GetTypeName())=="Long64_t") typestring = "L";
464     else if (TString(tdm->GetTypeName())=="ULong64_t") typestring = "l";
465     else if (TString(tdm->GetTypeName())=="Bool_t") typestring = "O";
466     else continue;
467     //printf("%s\n",tdm->GetTypeName());
468     TString varname = TString::Format("ph%d.%s",iph+1,tdm->GetName());
469    
470     Char_t *addr = (Char_t*)&fDiphotonEvent->photons[iph];
471     assert(sizeof(Char_t)==1);
472     if (fWriteDiphotonTree) hCiCTuple->Branch(varname,addr+tdm->GetOffset(),TString::Format("%s/%s",varname.Data(),typestring.Data()));
473    
474     if (iph==0) {
475     TString singlename = TString::Format("ph.%s",tdm->GetName());
476     Char_t *addrsingle = (Char_t*)fSinglePhoton;
477     if (fWriteSingleTree) hCiCTupleSingle->Branch(singlename,addrsingle+tdm->GetOffset(),TString::Format("%s/%s",singlename.Data(),typestring.Data()));
478     }
479     }
480     }
481    
482    
483     if (fWriteDiphotonTree) AddOutput(hCiCTuple);
484     if (fWriteSingleTree) AddOutput(hCiCTupleSingle);
485    
486    
487     }
488    
489     // ----------------------------------------------------------------------------------------
490     // some helpfer functions....
491     void PhotonTreeWriter::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass) {
492    
493     pt = -999.;
494     decayZ = -999.;
495     mass = -999.;
496    
497     // loop over all GEN particles and look for status 1 photons
498     for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
499     const MCParticle* p = fMCParticles->At(i);
500 bendavid 1.2 if( p->Is(MCParticle::kH) || (fInvertElectronVeto && (p->AbsPdgId()==23||p->AbsPdgId()==24) ) ) {
501 bendavid 1.1 pt=p->Pt();
502     decayZ = p->DecayVertex().Z();
503     mass = p->Mass();
504     break;
505     }
506     }
507    
508     return;
509     }
510    
511    
512     Float_t PhotonTreeWriter::GetEventCat(PhotonTools::CiCBaseLineCats cat1, PhotonTools::CiCBaseLineCats cat2) {
513    
514     bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
515     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
516    
517     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
518     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
519    
520     if( ph1IsEB && ph2IsEB )
521     return ( ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
522    
523     return ( ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
524     }
525    
526 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) {
527    
528     const SuperCluster *s = 0;
529     if (p) {
530     s = p->SCluster();
531     }
532     else {
533     s = ele->SCluster();
534     }
535     const BasicCluster *b = s->Seed();
536 bendavid 1.4 const BasicCluster *b2 = 0;
537     Double_t ebcmax = -99.;
538     for (UInt_t i=0; i<s->ClusterSize(); ++i) {
539     const BasicCluster *bc = s->Cluster(i);
540     if (bc->Energy() > ebcmax && bc !=b) {
541     b2 = bc;
542     ebcmax = bc->Energy();
543     }
544     }
545 bendavid 1.2
546     // if (p && ele) {
547     // printf("p : r9 = %5f, e33 = %5f, e55 = %5f, hovere = %5f, sigieie = %5f\n",p->R9(),p->E33(),p->E55(),p->HadOverEm(),p->CoviEtaiEta());
548     // printf("ele/b: r9 = %5f, e33 = %5f, e55 = %5f, hovere = %5f, sigieie = %5f\n",b->E3x3()/s->RawEnergy(),b->E3x3(),b->E5x5(),ele->HadronicOverEm(),ele->CoviEtaiEta());
549     // }
550    
551     if (p) {
552     hasphoton = kTRUE;
553     e = p->E();
554     pt = p->Pt();
555     eta = p->Eta();
556     phi = p->Phi();
557     r9 = p->R9();
558     e3x3 = p->E33();
559     e5x5 = p->E55();
560     hovere = p->HadOverEm();
561     sigietaieta = p->CoviEtaiEta();
562     phcat = PhotonTools::CiCBaseLineCat(p);
563     eerr = p->EnergyErr();
564     eerrsmeared = p->EnergyErrSmeared();
565 bendavid 1.3 esmearing = p->EnergySmearing();
566 bendavid 1.2 }
567     else {
568     hasphoton = kFALSE;
569     e = -99.;
570     pt = -99.;
571     eta = -99.;
572     phi = -99.;
573     r9 = b->E3x3()/s->RawEnergy();
574     e3x3 = b->E3x3();
575     e5x5 = b->E5x5();
576     hovere = ele->HadronicOverEm();
577     sigietaieta = ele->CoviEtaiEta();
578     phcat = -99;
579     eerr = -99.;
580     eerrsmeared = -99.;
581 bendavid 1.3 esmearing = 0.;
582 bendavid 1.2 }
583    
584    
585 bendavid 1.1 sce = s->Energy();
586     scrawe = s->RawEnergy();
587     scpse = s->PreshowerEnergy();
588     sceta = s->Eta();
589     scphi = s->Phi();
590     scnclusters = s->ClusterSize();
591     scnhits = s->NHits();
592 bendavid 1.2 scetawidth = s->EtaWidth();
593     scphiwidth = s->PhiWidth();
594 bendavid 1.1 isbarrel = (s->AbsEta()<1.5);
595     isr9reco = (isbarrel && r9>0.94) || (!isbarrel && r9>0.95);
596     isr9cat = (r9>0.94);
597    
598    
599    
600 bendavid 1.2
601 bendavid 1.1 sigiphiphi = TMath::Sqrt(b->CoviPhiiPhi());
602     if (isnan(sigiphiphi)) sigiphiphi = -99.;
603     covietaiphi = b->CoviEtaiPhi();
604     if (isnan(covietaiphi)) covietaiphi = -99.;
605     emax = b->EMax();
606     e2nd = b->E2nd();
607     etop = b->ETop();
608     ebottom = b->EBottom();
609     eleft = b->ELeft();
610     eright = b->ERight();
611     e1x3 = b->E1x3();
612     e3x1 = b->E3x1();
613     e1x5 = b->E1x5();
614     e2x2 = b->E2x2();
615     e4x4 = b->E4x4();
616     e2x5max = b->E2x5Max();
617     e2x5top = b->E2x5Top();
618     e2x5bottom = b->E2x5Bottom();
619     e2x5left = b->E2x5Left();
620     e2x5right = b->E2x5Right();
621 bendavid 1.4 xseed = b->Pos().X();
622     yseed = b->Pos().Y();
623     zseed = b->Pos().Z();
624    
625     eseed = b->Energy();
626     etaseed = b->Eta();
627     phiseed = b->Phi();
628     ietaseed = b->IEta();
629     iphiseed = b->IPhi();
630     ixseed = b->IX();
631     iyseed = b->IY();
632     etacryseed = b->EtaCry();
633     phicryseed = b->PhiCry();
634     xcryseed = b->XCry();
635     ycryseed = b->YCry();
636     thetaaxisseed = b->ThetaAxis();
637     phiaxisseed = b->PhiAxis();
638    
639     if (b2) {
640     ebc2 = b2->Energy();
641     etabc2 = b2->Eta();
642     phibc2 = b2->Phi();
643     ietabc2 = b2->IEta();
644     iphibc2 = b2->IPhi();
645     ixbc2 = b2->IX();
646     iybc2 = b2->IY();
647     etacrybc2 = b2->EtaCry();
648     phicrybc2 = b2->PhiCry();
649     xcrybc2 = b2->XCry();
650     ycrybc2 = b2->YCry();
651     thetaaxisbc2 = b2->ThetaAxis();
652     phiaxisbc2 = b2->PhiAxis();
653     }
654     else {
655     ebc2 = 0.;
656     etabc2 = 0.;
657     phibc2 = 0.;
658     ietabc2 = 0.;
659     iphibc2 = 0.;
660     ixbc2 = 0.;
661     iybc2 = 0.;
662     etacrybc2 = 0.;
663     phicrybc2 = 0.;
664     xcrybc2 = 0.;
665     ycrybc2 = 0.;
666     thetaaxisbc2 = 0.;
667     phiaxisbc2 = 0.;
668     }
669 bendavid 1.2
670 bendavid 1.1 //initialize photon energy corrections if needed
671 bendavid 1.2 /*if (!PhotonFix::initialised()) {
672 bendavid 1.1 PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
673 bendavid 1.2 }*/
674    
675    
676     phfixph.setup(e,sceta,scphi,r9);
677     phfixele.setup(e,sceta,scphi,r9);
678 bendavid 1.1
679     const Float_t dval = -99.;
680 bendavid 1.2 ecor = phfixph.fixedEnergy();
681     ecorerr = phfixph.sigmaEnergy();
682     ecorele = phfixele.fixedEnergy();
683     ecoreleerr = phfixele.sigmaEnergy();
684     if (phfixph.isbarrel()) {
685     etac = phfixph.etaC();
686     etas = phfixph.etaS();
687     etam = phfixph.etaM();
688     phic = phfixph.phiC();
689     phis = phfixph.phiS();
690     phim = phfixph.phiM();
691 bendavid 1.1 xz = dval;
692     xc = dval;
693     xs = dval;
694     xm = dval;
695     yz = dval;
696     yc = dval;
697     ys = dval;
698     ym = dval;
699     }
700     else {
701     etac = dval;
702     etas = dval;
703     etam = dval;
704     phic = dval;
705     phis = dval;
706     phim = dval;
707 bendavid 1.2 xz = phfixph.xZ();
708     xc = phfixph.xC();
709     xs = phfixph.xS();
710     xm = phfixph.xM();
711     yz = phfixph.yZ();
712     yc = phfixph.yC();
713     ys = phfixph.yS();
714     ym = phfixph.yM();
715 bendavid 1.1 }
716    
717     if (c) {
718     hasconversion = kTRUE;
719     convp = c->P();
720     convpt = c->Pt();
721     conveta = c->Eta();
722     convphi = c->Phi();
723 bendavid 1.2 ThreeVector dirconvsc = ThreeVector(s->Point()) - c->Position();
724 bendavid 1.1 convdeta = c->Eta() - dirconvsc.Eta();
725     convdphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
726     convvtxrho = c->Position().Rho();
727     convvtxz = c->Position().Z();
728     convvtxphi = c->Position().Phi();
729    
730     const StableData *leadsd = dynamic_cast<const StableData*>(c->DaughterDat(0));
731     const StableData *trailsd = dynamic_cast<const StableData*>(c->DaughterDat(1));
732     if (leadsd->Pt()<trailsd->Pt()) {
733     const StableData *sdtmp = leadsd;
734     leadsd = trailsd;
735     trailsd = sdtmp;
736     }
737    
738     const Track *leadtrack = dynamic_cast<const StableParticle*>(leadsd->Original())->Trk();
739     const Track *trailtrack = dynamic_cast<const StableParticle*>(trailsd->Original())->Trk();
740    
741     convleadpt = leadsd->Pt();
742     convtrailpt = trailsd->Pt();
743     convleadtrackpt = leadtrack->Pt();
744     convleadtrackalgo = leadtrack->Algo();
745     if (convleadtrackalgo==29) convleadtrackalgos=2; //gsf track
746     else if (convleadtrackalgo==15 ||convleadtrackalgo==16) convleadtrackalgos=1; //ecal-seeded track
747     else convleadtrackalgos = 0; //std iterative track
748     convleadtrackcharge = leadtrack->Charge();
749     convtrailtrackpt = trailtrack->Pt();
750     convtrailtrackalgo = trailtrack->Algo();
751     if (convtrailtrackalgo==29) convtrailtrackalgos=2; //gsf track
752     else if (convtrailtrackalgo==15 ||convtrailtrackalgo==16) convtrailtrackalgos=1; //ecal-seeded track
753     else convtrailtrackalgos = 0; //std iterative track
754     trailtrackcharge = trailtrack->Charge();
755     }
756     else {
757     hasconversion = kFALSE;
758     convp = -99.;
759     convpt = -99.;
760     conveta = -99.;
761     convphi = -99.;
762     convdeta = -99.;
763     convdphi = -99.;
764     convvtxrho = -99.;
765     convvtxz = -999.;
766     convvtxphi = -99.;
767     convleadpt = -99.;
768     convtrailpt = -99.;
769     convleadtrackpt = -99.;
770     convleadtrackalgo = -99;
771     convleadtrackalgos = -99;
772     convleadtrackcharge = 0;
773     convtrailtrackpt = -99.;
774     convtrailtrackalgo = -99;
775     convtrailtrackalgos = -99;
776     trailtrackcharge = 0;
777     }
778    
779 bendavid 1.2 //electron quantities
780     if (ele) {
781     haselectron = kTRUE;
782     eleisecaldriven = ele->IsEcalDriven();
783     eleistrackerdriven = ele->IsTrackerDriven();
784     elee = ele->E();
785     elept = ele->Pt();
786     eleeta = ele->Eta();
787     elephi = ele->Phi();
788     elecharge = ele->Charge();
789     elefbrem = ele->FBrem();
790     eledeta = ele->DeltaEtaSuperClusterTrackAtVtx();
791     eledphi = ele->DeltaPhiSuperClusterTrackAtVtx();
792     elep = s->Energy()/ele->ESuperClusterOverP();
793     elepin = ele->PIn();
794     elepout = ele->POut();
795     }
796     else {
797     haselectron = kFALSE;
798     eleisecaldriven = kFALSE;
799     eleistrackerdriven = kFALSE;
800     elee = -99.;
801     elept = -99.;
802     eleeta = -99.;
803     elephi = -99.;
804     elecharge = -99;
805     elefbrem = -99.;
806     eledeta = -99.;
807     eledphi = -99.;
808     elep = -99.;
809     elepin = -99.;
810     elepout = -99.;
811     }
812    
813     //pf supercluster quantities
814     if (pfsc) {
815     haspfsc = kTRUE;
816     pfsce = pfsc->Energy();
817     pfscrawe = pfsc->RawEnergy();
818     pfsceta = pfsc->Eta();
819     pfscphi = pfsc->Phi();
820     }
821     else {
822     haspfsc = kFALSE;
823     pfsce = -99.;
824     pfscrawe = -99.;
825     pfsceta = -99.;
826     pfscphi = -99.;
827     }
828    
829 bendavid 1.1 genz = -99.;
830     if (m) {
831     ispromptgen = kTRUE;
832     gene = m->E();
833     genpt = m->Pt();
834     geneta = m->Eta();
835     genphi = m->Phi();
836     const MCParticle *mm = m->DistinctMother();
837     if (mm) genz = mm->DecayVertex().Z();
838     }
839     else {
840     ispromptgen = kFALSE;
841     gene = -99.;
842     genpt = -99.;
843     geneta = -99.;
844     genphi = -99.;
845     }
846 bendavid 1.2
847 bendavid 1.1 }
848 bendavid 1.2