ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.3
Committed: Fri Oct 7 09:56:36 2011 UTC (13 years, 7 months ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025c, Mit_025b, Mit_025a, Mit_025
Changes since 1.2: +2 -0 lines
Log Message:
misc updates

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 bendavid 1.3 esmearing = p->EnergySmearing();
542 bendavid 1.2 }
543     else {
544     hasphoton = kFALSE;
545     e = -99.;
546     pt = -99.;
547     eta = -99.;
548     phi = -99.;
549     r9 = b->E3x3()/s->RawEnergy();
550     e3x3 = b->E3x3();
551     e5x5 = b->E5x5();
552     hovere = ele->HadronicOverEm();
553     sigietaieta = ele->CoviEtaiEta();
554     phcat = -99;
555     eerr = -99.;
556     eerrsmeared = -99.;
557 bendavid 1.3 esmearing = 0.;
558 bendavid 1.2 }
559    
560    
561 bendavid 1.1 sce = s->Energy();
562     scrawe = s->RawEnergy();
563     scpse = s->PreshowerEnergy();
564     sceta = s->Eta();
565     scphi = s->Phi();
566     scnclusters = s->ClusterSize();
567     scnhits = s->NHits();
568 bendavid 1.2 scetawidth = s->EtaWidth();
569     scphiwidth = s->PhiWidth();
570 bendavid 1.1 isbarrel = (s->AbsEta()<1.5);
571     isr9reco = (isbarrel && r9>0.94) || (!isbarrel && r9>0.95);
572     isr9cat = (r9>0.94);
573    
574    
575    
576 bendavid 1.2
577 bendavid 1.1 sigiphiphi = TMath::Sqrt(b->CoviPhiiPhi());
578     if (isnan(sigiphiphi)) sigiphiphi = -99.;
579     covietaiphi = b->CoviEtaiPhi();
580     if (isnan(covietaiphi)) covietaiphi = -99.;
581     emax = b->EMax();
582     e2nd = b->E2nd();
583     etop = b->ETop();
584     ebottom = b->EBottom();
585     eleft = b->ELeft();
586     eright = b->ERight();
587     e1x3 = b->E1x3();
588     e3x1 = b->E3x1();
589     e1x5 = b->E1x5();
590     e2x2 = b->E2x2();
591     e4x4 = b->E4x4();
592     e2x5max = b->E2x5Max();
593     e2x5top = b->E2x5Top();
594     e2x5bottom = b->E2x5Bottom();
595     e2x5left = b->E2x5Left();
596     e2x5right = b->E2x5Right();
597 bendavid 1.2 eseed = b->Energy();
598    
599 bendavid 1.1 //initialize photon energy corrections if needed
600 bendavid 1.2 /*if (!PhotonFix::initialised()) {
601 bendavid 1.1 PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
602 bendavid 1.2 }*/
603    
604    
605     phfixph.setup(e,sceta,scphi,r9);
606     phfixele.setup(e,sceta,scphi,r9);
607 bendavid 1.1
608     const Float_t dval = -99.;
609 bendavid 1.2 ecor = phfixph.fixedEnergy();
610     ecorerr = phfixph.sigmaEnergy();
611     ecorele = phfixele.fixedEnergy();
612     ecoreleerr = phfixele.sigmaEnergy();
613     if (phfixph.isbarrel()) {
614     etac = phfixph.etaC();
615     etas = phfixph.etaS();
616     etam = phfixph.etaM();
617     phic = phfixph.phiC();
618     phis = phfixph.phiS();
619     phim = phfixph.phiM();
620 bendavid 1.1 xz = dval;
621     xc = dval;
622     xs = dval;
623     xm = dval;
624     yz = dval;
625     yc = dval;
626     ys = dval;
627     ym = dval;
628     }
629     else {
630     etac = dval;
631     etas = dval;
632     etam = dval;
633     phic = dval;
634     phis = dval;
635     phim = dval;
636 bendavid 1.2 xz = phfixph.xZ();
637     xc = phfixph.xC();
638     xs = phfixph.xS();
639     xm = phfixph.xM();
640     yz = phfixph.yZ();
641     yc = phfixph.yC();
642     ys = phfixph.yS();
643     ym = phfixph.yM();
644 bendavid 1.1 }
645    
646     if (c) {
647     hasconversion = kTRUE;
648     convp = c->P();
649     convpt = c->Pt();
650     conveta = c->Eta();
651     convphi = c->Phi();
652 bendavid 1.2 ThreeVector dirconvsc = ThreeVector(s->Point()) - c->Position();
653 bendavid 1.1 convdeta = c->Eta() - dirconvsc.Eta();
654     convdphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
655     convvtxrho = c->Position().Rho();
656     convvtxz = c->Position().Z();
657     convvtxphi = c->Position().Phi();
658    
659     const StableData *leadsd = dynamic_cast<const StableData*>(c->DaughterDat(0));
660     const StableData *trailsd = dynamic_cast<const StableData*>(c->DaughterDat(1));
661     if (leadsd->Pt()<trailsd->Pt()) {
662     const StableData *sdtmp = leadsd;
663     leadsd = trailsd;
664     trailsd = sdtmp;
665     }
666    
667     const Track *leadtrack = dynamic_cast<const StableParticle*>(leadsd->Original())->Trk();
668     const Track *trailtrack = dynamic_cast<const StableParticle*>(trailsd->Original())->Trk();
669    
670     convleadpt = leadsd->Pt();
671     convtrailpt = trailsd->Pt();
672     convleadtrackpt = leadtrack->Pt();
673     convleadtrackalgo = leadtrack->Algo();
674     if (convleadtrackalgo==29) convleadtrackalgos=2; //gsf track
675     else if (convleadtrackalgo==15 ||convleadtrackalgo==16) convleadtrackalgos=1; //ecal-seeded track
676     else convleadtrackalgos = 0; //std iterative track
677     convleadtrackcharge = leadtrack->Charge();
678     convtrailtrackpt = trailtrack->Pt();
679     convtrailtrackalgo = trailtrack->Algo();
680     if (convtrailtrackalgo==29) convtrailtrackalgos=2; //gsf track
681     else if (convtrailtrackalgo==15 ||convtrailtrackalgo==16) convtrailtrackalgos=1; //ecal-seeded track
682     else convtrailtrackalgos = 0; //std iterative track
683     trailtrackcharge = trailtrack->Charge();
684     }
685     else {
686     hasconversion = kFALSE;
687     convp = -99.;
688     convpt = -99.;
689     conveta = -99.;
690     convphi = -99.;
691     convdeta = -99.;
692     convdphi = -99.;
693     convvtxrho = -99.;
694     convvtxz = -999.;
695     convvtxphi = -99.;
696     convleadpt = -99.;
697     convtrailpt = -99.;
698     convleadtrackpt = -99.;
699     convleadtrackalgo = -99;
700     convleadtrackalgos = -99;
701     convleadtrackcharge = 0;
702     convtrailtrackpt = -99.;
703     convtrailtrackalgo = -99;
704     convtrailtrackalgos = -99;
705     trailtrackcharge = 0;
706     }
707    
708 bendavid 1.2 //electron quantities
709     if (ele) {
710     haselectron = kTRUE;
711     eleisecaldriven = ele->IsEcalDriven();
712     eleistrackerdriven = ele->IsTrackerDriven();
713     elee = ele->E();
714     elept = ele->Pt();
715     eleeta = ele->Eta();
716     elephi = ele->Phi();
717     elecharge = ele->Charge();
718     elefbrem = ele->FBrem();
719     eledeta = ele->DeltaEtaSuperClusterTrackAtVtx();
720     eledphi = ele->DeltaPhiSuperClusterTrackAtVtx();
721     elep = s->Energy()/ele->ESuperClusterOverP();
722     elepin = ele->PIn();
723     elepout = ele->POut();
724     }
725     else {
726     haselectron = kFALSE;
727     eleisecaldriven = kFALSE;
728     eleistrackerdriven = kFALSE;
729     elee = -99.;
730     elept = -99.;
731     eleeta = -99.;
732     elephi = -99.;
733     elecharge = -99;
734     elefbrem = -99.;
735     eledeta = -99.;
736     eledphi = -99.;
737     elep = -99.;
738     elepin = -99.;
739     elepout = -99.;
740     }
741    
742     //pf supercluster quantities
743     if (pfsc) {
744     haspfsc = kTRUE;
745     pfsce = pfsc->Energy();
746     pfscrawe = pfsc->RawEnergy();
747     pfsceta = pfsc->Eta();
748     pfscphi = pfsc->Phi();
749     }
750     else {
751     haspfsc = kFALSE;
752     pfsce = -99.;
753     pfscrawe = -99.;
754     pfsceta = -99.;
755     pfscphi = -99.;
756     }
757    
758 bendavid 1.1 genz = -99.;
759     if (m) {
760     ispromptgen = kTRUE;
761     gene = m->E();
762     genpt = m->Pt();
763     geneta = m->Eta();
764     genphi = m->Phi();
765     const MCParticle *mm = m->DistinctMother();
766     if (mm) genz = mm->DecayVertex().Z();
767     }
768     else {
769     ispromptgen = kFALSE;
770     gene = -99.;
771     genpt = -99.;
772     geneta = -99.;
773     genphi = -99.;
774     }
775 bendavid 1.2
776 bendavid 1.1 }
777 bendavid 1.2