ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.2
Committed: Thu Sep 8 15:51:23 2011 UTC (13 years, 7 months ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025pre2, Mit_025pre1
Changes since 1.1: +307 -76 lines
Log Message:
Add tool and required changes for photon energy regression

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     }
236    
237     Float_t _mass = -99.;
238     Float_t _masserr = -99.;
239     Float_t _masserrsmeared = -99.;
240     Float_t _ptgg = -99.;
241     Float_t _costheta = -99.;
242     PhotonTools::DiphotonR9EtaPtCats _evtcat = PhotonTools::kOctCat0;
243     if (phHard && phSoft) {
244     _mass = (phHard->Mom()+phSoft->Mom()).M();
245     _masserr = 0.5*_mass*TMath::Sqrt(phHard->EnergyErr()*phHard->EnergyErr()/phHard->E()/phHard->E() + phSoft->EnergyErr()*phSoft->EnergyErr()/phSoft->E()/phSoft->E());
246     _masserrsmeared = 0.5*_mass*TMath::Sqrt(phHard->EnergyErrSmeared()*phHard->EnergyErrSmeared()/phHard->E()/phHard->E() + phSoft->EnergyErrSmeared()*phSoft->EnergyErrSmeared()/phSoft->E()/phSoft->E());
247     _ptgg = (phHard->Mom()+phSoft->Mom()).Pt();
248     _costheta = ThreeVector(phHard->Mom()).Unit().Dot(ThreeVector(phSoft->Mom()).Unit());
249     _evtcat = PhotonTools::DiphotonR9EtaPtCat(phHard,phSoft);
250     }
251    
252    
253     Float_t _massele = -99.;
254     Float_t _ptee = -99.;
255     Float_t _costhetaele = -99.;
256     if (ele1 && ele2) {
257     _massele = (ele1->Mom()+ele2->Mom()).M();
258     _ptee = (ele1->Mom()+ele2->Mom()).Pt();
259     _costhetaele = ThreeVector(ele1->Mom()).Unit().Dot(ThreeVector(ele2->Mom()).Unit());
260 bendavid 1.1 }
261    
262     Float_t _gencostheta = -99.;
263     if (phgen1 && phgen2) {
264     _gencostheta = ThreeVector(phgen1->Mom()).Unit().Dot(ThreeVector(phgen2->Mom()).Unit());
265     }
266    
267     fDiphotonEvent->gencostheta = _gencostheta;
268     fDiphotonEvent->mass = _mass;
269 bendavid 1.2 fDiphotonEvent->masserr = _masserr;
270     fDiphotonEvent->masserrsmeared = _masserrsmeared;
271 bendavid 1.1 fDiphotonEvent->ptgg = _ptgg;
272 bendavid 1.2 fDiphotonEvent->costheta = _costheta;;
273     fDiphotonEvent->massele = _massele;
274     fDiphotonEvent->ptee = _ptee;
275     fDiphotonEvent->costhetaele = _costhetaele;
276     fDiphotonEvent->evtcat = _evtcat;
277 bendavid 1.1
278 bendavid 1.2 fDiphotonEvent->photons[0].SetVars(phHard,conv1,ele1,pfsc1,phgen1,fPhfixph,fPhfixele);
279     fDiphotonEvent->photons[1].SetVars(phSoft,conv2,ele2,pfsc2,phgen2,fPhfixph,fPhfixele);
280 bendavid 1.1
281     Float_t ph1ecor = fDiphotonEvent->photons[0].Ecor();
282     Float_t ph1ecorerr = fDiphotonEvent->photons[0].Ecorerr();
283     Float_t ph2ecor = fDiphotonEvent->photons[1].Ecor();
284     Float_t ph2ecorerr = fDiphotonEvent->photons[1].Ecorerr();
285 bendavid 1.2
286     Float_t ph1ecorele = fDiphotonEvent->photons[0].Ecorele();
287     Float_t ph1ecoreleerr = fDiphotonEvent->photons[0].Ecoreleerr();
288     Float_t ph2ecorele = fDiphotonEvent->photons[1].Ecorele();
289     Float_t ph2ecoreleerr = fDiphotonEvent->photons[1].Ecoreleerr();
290    
291 bendavid 1.1
292     fDiphotonEvent->masscor = TMath::Sqrt(2.0*ph1ecor*ph2ecor*(1.0-fDiphotonEvent->costheta));
293     fDiphotonEvent->masscorerr = 0.5*fDiphotonEvent->masscor*TMath::Sqrt(ph1ecorerr*ph1ecorerr/ph1ecor/ph1ecor + ph2ecorerr*ph2ecorerr/ph2ecor/ph2ecor);
294    
295 bendavid 1.2 fDiphotonEvent->masscorele = TMath::Sqrt(2.0*ph1ecorele*ph2ecorele*(1.0-fDiphotonEvent->costheta));
296     fDiphotonEvent->masscoreleerr = 0.5*fDiphotonEvent->masscorele*TMath::Sqrt(ph1ecoreleerr*ph1ecoreleerr/ph1ecorele/ph1ecorele + ph2ecoreleerr*ph2ecoreleerr/ph2ecorele/ph2ecorele);
297    
298 bendavid 1.1 //printf("r9 = %5f, photon sigieie = %5f, seed sigieie = %5f\n",phHard->R9(),phHard->CoviEtaiEta(),sqrt(phHard->SCluster()->Seed()->CoviEtaiEta()));
299    
300     if (fWriteDiphotonTree) hCiCTuple->Fill();
301    
302     }
303    
304     if (!fWriteSingleTree) return;
305    
306    
307 bendavid 1.2 for (UInt_t iph = 0; iph<egcol->GetEntries(); ++iph) {
308 bendavid 1.1
309 bendavid 1.2 const Particle *p = 0;
310     const Photon *ph = 0;
311     const Electron *ele = 0;
312     const SuperCluster *sc = 0;
313     if (fLoopOnGoodElectrons) {
314     ele = fGoodElectrons->At(iph);
315     p = ele;
316     sc = ele->SCluster();
317     ph = PhotonTools::MatchedPhoton(ele,fPhotons);
318     }
319     else {
320     ph = fPhotons->At(iph);
321     p = ph;
322     sc = ph->SCluster();
323     ele = PhotonTools::MatchedElectron(ph,fGoodElectrons);
324     }
325    
326     const DecayParticle *conv = PhotonTools::MatchedConversion(sc,fConversions,bsp,nhitsbeforevtxmax);
327     const SuperCluster *pfsc = PhotonTools::MatchedSC(sc,fSuperClusters);
328    
329    
330    
331     if (!fLoopOnGoodElectrons && ph->HasPV()) {
332 bendavid 1.1 fDiphotonEvent->vtxZ = ph->PV()->Z();
333     }
334    
335     const MCParticle *phgen = 0;
336     if( !fIsData ) {
337 bendavid 1.2 phgen = PhotonTools::MatchMC(p,fMCParticles,fInvertElectronVeto);
338 bendavid 1.1 }
339    
340 bendavid 1.2 if (fExcludeSinglePrompt && phgen) return;
341    
342     fDiphotonEvent->mt = -99.;
343     fDiphotonEvent->cosphimet = -99.;
344     fDiphotonEvent->mtele = -99.;
345     fDiphotonEvent->cosphimetele = -99.;
346    
347     if (ph) {
348     fDiphotonEvent->cosphimet = TMath::Cos(ph->Phi()-fPFMet->At(0)->Phi());
349     fDiphotonEvent->mt = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ph->Pt()*(1.0-fDiphotonEvent->cosphimet));
350     }
351 bendavid 1.1
352 bendavid 1.2 if (ele) {
353     fDiphotonEvent->cosphimetele = TMath::Cos(ele->Phi()-fPFMet->At(0)->Phi());
354     fDiphotonEvent->mtele = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ele->Pt()*(1.0-fDiphotonEvent->cosphimetele));
355     }
356    
357     fSinglePhoton->SetVars(ph,conv,ele,pfsc,phgen,fPhfixph,fPhfixele);
358 bendavid 1.1 hCiCTupleSingle->Fill();
359    
360     }
361    
362    
363     return;
364    
365     }
366    
367     //--------------------------------------------------------------------------------------------------
368     void PhotonTreeWriter::SlaveBegin()
369     {
370     // Run startup code on the computer (slave) doing the actual analysis. Here,
371     // we just request the photon collection branch.
372    
373     ReqEventObject(fPhotonBranchName, fPhotons, fPhotonsFromBranch);
374     ReqEventObject(fTrackBranchName, fTracks, true);
375     ReqEventObject(fElectronName, fElectrons, true);
376     ReqEventObject(fGoodElectronName, fGoodElectrons, fGoodElectronsFromBranch);
377     ReqEventObject(fPileUpDenName, fPileUpDen, true);
378     ReqEventObject(fPVName, fPV, fPVFromBranch);
379     ReqEventObject(fConversionName, fConversions,true);
380     ReqEventObject(fBeamspotName, fBeamspot, true);
381     ReqEventObject(fPFCandName, fPFCands, true);
382 bendavid 1.2 ReqEventObject(fSuperClusterName, fSuperClusters, true);
383     ReqEventObject(fPFMetName, fPFMet, true);
384 bendavid 1.1
385     if (!fIsData) {
386     ReqBranch(fPileUpName, fPileUp);
387     ReqBranch(fMCParticleName, fMCParticles);
388     }
389    
390    
391    
392     //initialize photon energy corrections
393     //PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
394    
395 bendavid 1.2 fPhfixph.initialise("4_2",std::string(fPhFixDataFile));
396     fPhfixele.initialise("4_2e",std::string(fPhFixDataFile));
397    
398 bendavid 1.1 fDiphotonEvent = new PhotonTreeWriterDiphotonEvent;
399     fSinglePhoton = new PhotonTreeWriterPhoton;
400    
401    
402     if (fWriteDiphotonTree) hCiCTuple = new TTree(fTupleName.Data(),fTupleName.Data());
403     TString singlename = fTupleName + TString("Single");
404     if (fWriteSingleTree) hCiCTupleSingle = new TTree(singlename,singlename);
405    
406     //make flattish tree from classes so we don't have to rely on dictionaries for reading later
407     TClass *eclass = TClass::GetClass("mithep::PhotonTreeWriterDiphotonEvent");
408     TClass *pclass = TClass::GetClass("mithep::PhotonTreeWriterPhoton");
409     TList *elist = eclass->GetListOfDataMembers();
410     TList *plist = pclass->GetListOfDataMembers();
411    
412     for (int i=0; i<elist->GetEntries(); ++i) {
413     const TDataMember *tdm = static_cast<const TDataMember*>(elist->At(i));
414     if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
415     TString typestring;
416     if (TString(tdm->GetTypeName())=="Char_t") typestring = "B";
417     else if (TString(tdm->GetTypeName())=="UChar_t") typestring = "b";
418     else if (TString(tdm->GetTypeName())=="Short_t") typestring = "S";
419     else if (TString(tdm->GetTypeName())=="UShort_t") typestring = "s";
420     else if (TString(tdm->GetTypeName())=="Int_t") typestring = "I";
421     else if (TString(tdm->GetTypeName())=="UInt_t") typestring = "i";
422     else if (TString(tdm->GetTypeName())=="Float_t") typestring = "F";
423     else if (TString(tdm->GetTypeName())=="Double_t") typestring = "D";
424     else if (TString(tdm->GetTypeName())=="Long64_t") typestring = "L";
425     else if (TString(tdm->GetTypeName())=="ULong64_t") typestring = "l";
426     else if (TString(tdm->GetTypeName())=="Bool_t") typestring = "O";
427     else continue;
428     //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
429     Char_t *addr = (Char_t*)fDiphotonEvent;
430     assert(sizeof(Char_t)==1);
431     if (fWriteDiphotonTree) hCiCTuple->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
432     if (fWriteSingleTree) hCiCTupleSingle->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
433     }
434    
435     for (int iph=0; iph<2; ++iph) {
436     for (int i=0; i<plist->GetEntries(); ++i) {
437     const TDataMember *tdm = static_cast<const TDataMember*>(plist->At(i));
438     if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
439     TString typestring;
440     if (TString(tdm->GetTypeName())=="Char_t") typestring = "B";
441     else if (TString(tdm->GetTypeName())=="UChar_t") typestring = "b";
442     else if (TString(tdm->GetTypeName())=="Short_t") typestring = "S";
443     else if (TString(tdm->GetTypeName())=="UShort_t") typestring = "s";
444     else if (TString(tdm->GetTypeName())=="Int_t") typestring = "I";
445     else if (TString(tdm->GetTypeName())=="UInt_t") typestring = "i";
446     else if (TString(tdm->GetTypeName())=="Float_t") typestring = "F";
447     else if (TString(tdm->GetTypeName())=="Double_t") typestring = "D";
448     else if (TString(tdm->GetTypeName())=="Long64_t") typestring = "L";
449     else if (TString(tdm->GetTypeName())=="ULong64_t") typestring = "l";
450     else if (TString(tdm->GetTypeName())=="Bool_t") typestring = "O";
451     else continue;
452     //printf("%s\n",tdm->GetTypeName());
453     TString varname = TString::Format("ph%d.%s",iph+1,tdm->GetName());
454    
455     Char_t *addr = (Char_t*)&fDiphotonEvent->photons[iph];
456     assert(sizeof(Char_t)==1);
457     if (fWriteDiphotonTree) hCiCTuple->Branch(varname,addr+tdm->GetOffset(),TString::Format("%s/%s",varname.Data(),typestring.Data()));
458    
459     if (iph==0) {
460     TString singlename = TString::Format("ph.%s",tdm->GetName());
461     Char_t *addrsingle = (Char_t*)fSinglePhoton;
462     if (fWriteSingleTree) hCiCTupleSingle->Branch(singlename,addrsingle+tdm->GetOffset(),TString::Format("%s/%s",singlename.Data(),typestring.Data()));
463     }
464     }
465     }
466    
467    
468     if (fWriteDiphotonTree) AddOutput(hCiCTuple);
469     if (fWriteSingleTree) AddOutput(hCiCTupleSingle);
470    
471    
472     }
473    
474     // ----------------------------------------------------------------------------------------
475     // some helpfer functions....
476     void PhotonTreeWriter::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass) {
477    
478     pt = -999.;
479     decayZ = -999.;
480     mass = -999.;
481    
482     // loop over all GEN particles and look for status 1 photons
483     for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
484     const MCParticle* p = fMCParticles->At(i);
485 bendavid 1.2 if( p->Is(MCParticle::kH) || (fInvertElectronVeto && (p->AbsPdgId()==23||p->AbsPdgId()==24) ) ) {
486 bendavid 1.1 pt=p->Pt();
487     decayZ = p->DecayVertex().Z();
488     mass = p->Mass();
489     break;
490     }
491     }
492    
493     return;
494     }
495    
496    
497     Float_t PhotonTreeWriter::GetEventCat(PhotonTools::CiCBaseLineCats cat1, PhotonTools::CiCBaseLineCats cat2) {
498    
499     bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
500     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
501    
502     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
503     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
504    
505     if( ph1IsEB && ph2IsEB )
506     return ( ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
507    
508     return ( ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
509     }
510    
511 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) {
512    
513     const SuperCluster *s = 0;
514     if (p) {
515     s = p->SCluster();
516     }
517     else {
518     s = ele->SCluster();
519     }
520     const BasicCluster *b = s->Seed();
521    
522     // if (p && ele) {
523     // printf("p : r9 = %5f, e33 = %5f, e55 = %5f, hovere = %5f, sigieie = %5f\n",p->R9(),p->E33(),p->E55(),p->HadOverEm(),p->CoviEtaiEta());
524     // 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());
525     // }
526    
527     if (p) {
528     hasphoton = kTRUE;
529     e = p->E();
530     pt = p->Pt();
531     eta = p->Eta();
532     phi = p->Phi();
533     r9 = p->R9();
534     e3x3 = p->E33();
535     e5x5 = p->E55();
536     hovere = p->HadOverEm();
537     sigietaieta = p->CoviEtaiEta();
538     phcat = PhotonTools::CiCBaseLineCat(p);
539     eerr = p->EnergyErr();
540     eerrsmeared = p->EnergyErrSmeared();
541     }
542     else {
543     hasphoton = kFALSE;
544     e = -99.;
545     pt = -99.;
546     eta = -99.;
547     phi = -99.;
548     r9 = b->E3x3()/s->RawEnergy();
549     e3x3 = b->E3x3();
550     e5x5 = b->E5x5();
551     hovere = ele->HadronicOverEm();
552     sigietaieta = ele->CoviEtaiEta();
553     phcat = -99;
554     eerr = -99.;
555     eerrsmeared = -99.;
556     }
557    
558    
559 bendavid 1.1 sce = s->Energy();
560     scrawe = s->RawEnergy();
561     scpse = s->PreshowerEnergy();
562     sceta = s->Eta();
563     scphi = s->Phi();
564     scnclusters = s->ClusterSize();
565     scnhits = s->NHits();
566 bendavid 1.2 scetawidth = s->EtaWidth();
567     scphiwidth = s->PhiWidth();
568 bendavid 1.1 isbarrel = (s->AbsEta()<1.5);
569     isr9reco = (isbarrel && r9>0.94) || (!isbarrel && r9>0.95);
570     isr9cat = (r9>0.94);
571    
572    
573    
574 bendavid 1.2
575 bendavid 1.1 sigiphiphi = TMath::Sqrt(b->CoviPhiiPhi());
576     if (isnan(sigiphiphi)) sigiphiphi = -99.;
577     covietaiphi = b->CoviEtaiPhi();
578     if (isnan(covietaiphi)) covietaiphi = -99.;
579     emax = b->EMax();
580     e2nd = b->E2nd();
581     etop = b->ETop();
582     ebottom = b->EBottom();
583     eleft = b->ELeft();
584     eright = b->ERight();
585     e1x3 = b->E1x3();
586     e3x1 = b->E3x1();
587     e1x5 = b->E1x5();
588     e2x2 = b->E2x2();
589     e4x4 = b->E4x4();
590     e2x5max = b->E2x5Max();
591     e2x5top = b->E2x5Top();
592     e2x5bottom = b->E2x5Bottom();
593     e2x5left = b->E2x5Left();
594     e2x5right = b->E2x5Right();
595 bendavid 1.2 eseed = b->Energy();
596    
597 bendavid 1.1 //initialize photon energy corrections if needed
598 bendavid 1.2 /*if (!PhotonFix::initialised()) {
599 bendavid 1.1 PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
600 bendavid 1.2 }*/
601    
602    
603     phfixph.setup(e,sceta,scphi,r9);
604     phfixele.setup(e,sceta,scphi,r9);
605 bendavid 1.1
606     const Float_t dval = -99.;
607 bendavid 1.2 ecor = phfixph.fixedEnergy();
608     ecorerr = phfixph.sigmaEnergy();
609     ecorele = phfixele.fixedEnergy();
610     ecoreleerr = phfixele.sigmaEnergy();
611     if (phfixph.isbarrel()) {
612     etac = phfixph.etaC();
613     etas = phfixph.etaS();
614     etam = phfixph.etaM();
615     phic = phfixph.phiC();
616     phis = phfixph.phiS();
617     phim = phfixph.phiM();
618 bendavid 1.1 xz = dval;
619     xc = dval;
620     xs = dval;
621     xm = dval;
622     yz = dval;
623     yc = dval;
624     ys = dval;
625     ym = dval;
626     }
627     else {
628     etac = dval;
629     etas = dval;
630     etam = dval;
631     phic = dval;
632     phis = dval;
633     phim = dval;
634 bendavid 1.2 xz = phfixph.xZ();
635     xc = phfixph.xC();
636     xs = phfixph.xS();
637     xm = phfixph.xM();
638     yz = phfixph.yZ();
639     yc = phfixph.yC();
640     ys = phfixph.yS();
641     ym = phfixph.yM();
642 bendavid 1.1 }
643    
644     if (c) {
645     hasconversion = kTRUE;
646     convp = c->P();
647     convpt = c->Pt();
648     conveta = c->Eta();
649     convphi = c->Phi();
650 bendavid 1.2 ThreeVector dirconvsc = ThreeVector(s->Point()) - c->Position();
651 bendavid 1.1 convdeta = c->Eta() - dirconvsc.Eta();
652     convdphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
653     convvtxrho = c->Position().Rho();
654     convvtxz = c->Position().Z();
655     convvtxphi = c->Position().Phi();
656    
657     const StableData *leadsd = dynamic_cast<const StableData*>(c->DaughterDat(0));
658     const StableData *trailsd = dynamic_cast<const StableData*>(c->DaughterDat(1));
659     if (leadsd->Pt()<trailsd->Pt()) {
660     const StableData *sdtmp = leadsd;
661     leadsd = trailsd;
662     trailsd = sdtmp;
663     }
664    
665     const Track *leadtrack = dynamic_cast<const StableParticle*>(leadsd->Original())->Trk();
666     const Track *trailtrack = dynamic_cast<const StableParticle*>(trailsd->Original())->Trk();
667    
668     convleadpt = leadsd->Pt();
669     convtrailpt = trailsd->Pt();
670     convleadtrackpt = leadtrack->Pt();
671     convleadtrackalgo = leadtrack->Algo();
672     if (convleadtrackalgo==29) convleadtrackalgos=2; //gsf track
673     else if (convleadtrackalgo==15 ||convleadtrackalgo==16) convleadtrackalgos=1; //ecal-seeded track
674     else convleadtrackalgos = 0; //std iterative track
675     convleadtrackcharge = leadtrack->Charge();
676     convtrailtrackpt = trailtrack->Pt();
677     convtrailtrackalgo = trailtrack->Algo();
678     if (convtrailtrackalgo==29) convtrailtrackalgos=2; //gsf track
679     else if (convtrailtrackalgo==15 ||convtrailtrackalgo==16) convtrailtrackalgos=1; //ecal-seeded track
680     else convtrailtrackalgos = 0; //std iterative track
681     trailtrackcharge = trailtrack->Charge();
682     }
683     else {
684     hasconversion = kFALSE;
685     convp = -99.;
686     convpt = -99.;
687     conveta = -99.;
688     convphi = -99.;
689     convdeta = -99.;
690     convdphi = -99.;
691     convvtxrho = -99.;
692     convvtxz = -999.;
693     convvtxphi = -99.;
694     convleadpt = -99.;
695     convtrailpt = -99.;
696     convleadtrackpt = -99.;
697     convleadtrackalgo = -99;
698     convleadtrackalgos = -99;
699     convleadtrackcharge = 0;
700     convtrailtrackpt = -99.;
701     convtrailtrackalgo = -99;
702     convtrailtrackalgos = -99;
703     trailtrackcharge = 0;
704     }
705    
706 bendavid 1.2 //electron quantities
707     if (ele) {
708     haselectron = kTRUE;
709     eleisecaldriven = ele->IsEcalDriven();
710     eleistrackerdriven = ele->IsTrackerDriven();
711     elee = ele->E();
712     elept = ele->Pt();
713     eleeta = ele->Eta();
714     elephi = ele->Phi();
715     elecharge = ele->Charge();
716     elefbrem = ele->FBrem();
717     eledeta = ele->DeltaEtaSuperClusterTrackAtVtx();
718     eledphi = ele->DeltaPhiSuperClusterTrackAtVtx();
719     elep = s->Energy()/ele->ESuperClusterOverP();
720     elepin = ele->PIn();
721     elepout = ele->POut();
722     }
723     else {
724     haselectron = kFALSE;
725     eleisecaldriven = kFALSE;
726     eleistrackerdriven = kFALSE;
727     elee = -99.;
728     elept = -99.;
729     eleeta = -99.;
730     elephi = -99.;
731     elecharge = -99;
732     elefbrem = -99.;
733     eledeta = -99.;
734     eledphi = -99.;
735     elep = -99.;
736     elepin = -99.;
737     elepout = -99.;
738     }
739    
740     //pf supercluster quantities
741     if (pfsc) {
742     haspfsc = kTRUE;
743     pfsce = pfsc->Energy();
744     pfscrawe = pfsc->RawEnergy();
745     pfsceta = pfsc->Eta();
746     pfscphi = pfsc->Phi();
747     }
748     else {
749     haspfsc = kFALSE;
750     pfsce = -99.;
751     pfscrawe = -99.;
752     pfsceta = -99.;
753     pfscphi = -99.;
754     }
755    
756 bendavid 1.1 genz = -99.;
757     if (m) {
758     ispromptgen = kTRUE;
759     gene = m->E();
760     genpt = m->Pt();
761     geneta = m->Eta();
762     genphi = m->Phi();
763     const MCParticle *mm = m->DistinctMother();
764     if (mm) genz = mm->DecayVertex().Z();
765     }
766     else {
767     ispromptgen = kFALSE;
768     gene = -99.;
769     genpt = -99.;
770     geneta = -99.;
771     genphi = -99.;
772     }
773 bendavid 1.2
774 bendavid 1.1 }
775 bendavid 1.2