ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.26
Committed: Thu Jun 7 14:24:07 2012 UTC (12 years, 11 months ago) by mingyang
Content type: text/plain
Branch: MAIN
Changes since 1.25: +18 -5 lines
Log Message: 
id variable added

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 fabstoec 1.16 #include "MitPhysics/Utils/interface/PFMetCorrectionTools.h"
12     #include "MitAna/DataTree/interface/PFJetCol.h"
13     #include "MitAna/DataTree/interface/GenJetCol.h"
14 bendavid 1.1 #include "TDataMember.h"
15 bendavid 1.20 #include "TFile.h"
16 bendavid 1.1 #include <TNtuple.h>
17     #include <TRandom3.h>
18     #include <TSystem.h>
19    
20     using namespace mithep;
21    
22     ClassImp(mithep::PhotonTreeWriter)
23 mingyang 1.26 templateClassImp(mithep::PhotonTreeWriterPhoton)//ming: what's this?
24 bendavid 1.1 ClassImp(mithep::PhotonTreeWriterDiphotonEvent)
25    
26     //--------------------------------------------------------------------------------------------------
27     PhotonTreeWriter::PhotonTreeWriter(const char *name, const char *title) :
28     // Base Module...
29 paus 1.12 BaseMod (name,title),
30 bendavid 1.1 // define all the Branches to load
31 paus 1.12 fPhotonBranchName (Names::gkPhotonBrn),
32 bendavid 1.20 fPFPhotonName ("PFPhotons"),
33 paus 1.12 fElectronName (Names::gkElectronBrn),
34     fGoodElectronName (Names::gkElectronBrn),
35     fConversionName (Names::gkMvfConversionBrn),
36     fTrackBranchName (Names::gkTrackBrn),
37     fPileUpDenName (Names::gkPileupEnergyDensityBrn),
38     fPVName (Names::gkPVBeamSpotBrn),
39     fBeamspotName (Names::gkBeamSpotBrn),
40     fPFCandName (Names::gkPFCandidatesBrn),
41     fMCParticleName (Names::gkMCPartBrn),
42     fPileUpName (Names::gkPileupInfoBrn),
43     fSuperClusterName ("PFSuperClusters"),
44     fPFMetName ("PFMet"),
45     fPFJetName (Names::gkPFJetBrn),
46 fabstoec 1.19 funcorrPFJetName ("AKt5PFJets"),
47 fabstoec 1.17 fGenJetName ("AKT5GenJets"),
48 fabstoec 1.13 fLeptonTagElectronsName ("HggLeptonTagElectrons"),
49     fLeptonTagMuonsName ("HggLeptonTagMuons"),
50    
51 paus 1.12 fIsData (false),
52     fPhotonsFromBranch (kTRUE),
53     fPVFromBranch (kTRUE),
54     fGoodElectronsFromBranch(kTRUE),
55     fPFJetsFromBranch (kTRUE),
56 bendavid 1.1 // ----------------------------------------
57     // collections....
58 paus 1.12 fPhotons (0),
59 bendavid 1.20 fPFPhotons (0),
60 paus 1.12 fElectrons (0),
61     fConversions (0),
62     fTracks (0),
63     fPileUpDen (0),
64     fPV (0),
65     fBeamspot (0),
66     fPFCands (0),
67     fMCParticles (0),
68     fPileUp (0),
69     fSuperClusters (0),
70     fPFJets (0),
71 fabstoec 1.16 fGenJets (0),
72     funcorrPFJets (0),
73 fabstoec 1.13
74     fLeptonTagElectrons (0),
75     fLeptonTagMuons (0),
76    
77 paus 1.12 fLoopOnGoodElectrons (kFALSE),
78     fApplyElectronVeto (kTRUE),
79     fWriteDiphotonTree (kTRUE),
80     fWriteSingleTree (kTRUE),
81 bendavid 1.20 fEnablePFPhotons (kTRUE),
82 paus 1.12 fExcludeSinglePrompt (kFALSE),
83     fExcludeDoublePrompt (kFALSE),
84     fEnableJets (kFALSE),
85 bendavid 1.22 fApplyJetId (kFALSE),
86 fabstoec 1.13 fApplyLeptonTag (kFALSE),
87 fabstoec 1.15 fApplyBTag (kFALSE),
88 fabstoec 1.17 fApplyPFMetCorrections (kFALSE),
89 bendavid 1.20 fFillClusterArrays (kFALSE),
90 bendavid 1.22 fFillVertexTree (kFALSE),
91 paus 1.12 fPhFixDataFile (gSystem->Getenv("CMSSW_BASE") +
92     TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat")),
93     fTupleName ("hPhotonTree")
94 bendavid 1.1 {
95 paus 1.12 // Constructor
96 bendavid 1.1 }
97    
98 paus 1.12 PhotonTreeWriter::~PhotonTreeWriter()
99     {
100     // Destructor
101 bendavid 1.1 }
102    
103     //--------------------------------------------------------------------------------------------------
104     void PhotonTreeWriter::Process()
105     {
106     // ------------------------------------------------------------
107     // Process entries of the tree.
108     LoadEventObject(fPhotonBranchName, fPhotons);
109 paus 1.12 LoadEventObject(fGoodElectronName, fGoodElectrons);
110 fabstoec 1.13
111     // lepton tag collections
112     if( fApplyLeptonTag ) {
113     LoadEventObject(fLeptonTagElectronsName, fLeptonTagElectrons);
114     LoadEventObject(fLeptonTagMuonsName, fLeptonTagMuons);
115     }
116    
117 bendavid 1.2 const BaseCollection *egcol = 0;
118 paus 1.12 if (fLoopOnGoodElectrons)
119 bendavid 1.2 egcol = fGoodElectrons;
120 paus 1.12 else
121 bendavid 1.2 egcol = fPhotons;
122 paus 1.12 if (egcol->GetEntries()<1)
123     return;
124 bendavid 1.20
125     if (fEnablePFPhotons) LoadEventObject(fPFPhotonName, fPFPhotons);
126 bendavid 1.1 LoadEventObject(fElectronName, fElectrons);
127     LoadEventObject(fConversionName, fConversions);
128     LoadEventObject(fTrackBranchName, fTracks);
129     LoadEventObject(fPileUpDenName, fPileUpDen);
130     LoadEventObject(fPVName, fPV);
131     LoadEventObject(fBeamspotName, fBeamspot);
132     LoadEventObject(fPFCandName, fPFCands);
133 bendavid 1.2 LoadEventObject(fSuperClusterName, fSuperClusters);
134 paus 1.12 LoadEventObject(fPFMetName, fPFMet);
135 fabstoec 1.16 if (fEnableJets){
136 paus 1.12 LoadEventObject(fPFJetName, fPFJets);
137 fabstoec 1.19 //LoadEventObject(funcorrPFJetName, funcorrPFJets);
138 fabstoec 1.16 LoadBranch(funcorrPFJetName);
139 fabstoec 1.19 // if(!fIsData) LoadEventObject(fGenJetName, fGenJets);
140 fabstoec 1.16 }
141 bendavid 1.1 // ------------------------------------------------------------
142     // load event based information
143 paus 1.12 Int_t _numPU = -1.; // some sensible default values....
144 bendavid 1.1 Int_t _numPUminus = -1.; // some sensible default values....
145 paus 1.12 Int_t _numPUplus = -1.; // some sensible default values....
146 bendavid 1.1
147 paus 1.12 Float_t rho = -99.;
148 bendavid 1.1 if( fPileUpDen->GetEntries() > 0 )
149 paus 1.12 rho = (Double_t) fPileUpDen->At(0)->RhoRandomLowEta();
150 bendavid 1.1
151 mingyang 1.26 const BaseVertex *bsp = dynamic_cast<const BaseVertex*>(fBeamspot->At(0));//ming:what's this?
152 bendavid 1.20
153 bendavid 1.1 if( !fIsData ) {
154     LoadBranch(fMCParticleName);
155     LoadBranch(fPileUpName);
156 fabstoec 1.19 if (fEnableJets) LoadEventObject(fGenJetName, fGenJets);
157     } else fGenJets = NULL;
158 bendavid 1.1
159     if( !fIsData ) {
160     for (UInt_t i=0; i<fPileUp->GetEntries(); ++i) {
161     const PileupInfo *puinfo = fPileUp->At(i);
162     if (puinfo->GetBunchCrossing()==0) _numPU = puinfo->GetPU_NumInteractions();
163     else if (puinfo->GetBunchCrossing() == -1) _numPUminus = puinfo->GetPU_NumInteractions();
164     else if (puinfo->GetBunchCrossing() == 1) _numPUplus = puinfo->GetPU_NumInteractions();
165     }
166     }
167    
168     // in case of a MC event, try to find Higgs and Higgs decay Z poisition
169 paus 1.12 Float_t _pth = -100.;
170     Float_t _decayZ = -100.;
171 bendavid 1.1 Float_t _genmass = -100.;
172 paus 1.12 if (!fIsData)
173     FindHiggsPtAndZ(_pth, _decayZ, _genmass);
174 bendavid 1.1
175 fabstoec 1.19
176 mtouch 1.25 Double_t _evt = GetEventHeader()->EvtNum();
177 fabstoec 1.19 Double_t _spfMet = fPFMet->At(0)->SumEt();
178 mtouch 1.25
179 fabstoec 1.13 fDiphotonEvent->leptonTag = -1; // disabled
180    
181 bendavid 1.20 fDiphotonEvent->rho = fPileUpDen->At(0)->RhoKt6PFJets();
182     fDiphotonEvent->rho25 = fPileUpDen->At(0)->RhoRandomLowEta();
183     fDiphotonEvent->rhoold = fPileUpDen->At(0)->Rho();
184 bendavid 1.1 fDiphotonEvent->genHiggspt = _pth;
185     fDiphotonEvent->genHiggsZ = _decayZ;
186     fDiphotonEvent->genmass = _genmass;
187     fDiphotonEvent->gencostheta = -99.;
188     fDiphotonEvent->nVtx = fPV->GetEntries();
189 bendavid 1.2 fDiphotonEvent->bsX = fBeamspot->At(0)->X();
190     fDiphotonEvent->bsY = fBeamspot->At(0)->Y();
191     fDiphotonEvent->bsZ = fBeamspot->At(0)->Z();
192 bendavid 1.7 fDiphotonEvent->bsSigmaZ = fBeamspot->At(0)->SigmaZ();
193 bendavid 1.2 fDiphotonEvent->vtxX = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->X() : -99.;
194     fDiphotonEvent->vtxY = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Y() : -99.;
195 bendavid 1.1 fDiphotonEvent->vtxZ = (fDiphotonEvent->nVtx>0) ? fPV->At(0)->Z() : -99.;
196     fDiphotonEvent->numPU = _numPU;
197     fDiphotonEvent->numPUminus = _numPUminus;
198     fDiphotonEvent->numPUplus = _numPUplus;
199     fDiphotonEvent->mass = -99.;
200     fDiphotonEvent->ptgg = -99.;
201     fDiphotonEvent->costheta = -99.;
202 bendavid 1.2 fDiphotonEvent->mt = -99.;
203     fDiphotonEvent->cosphimet = -99.;
204     fDiphotonEvent->mtele = -99.;
205     fDiphotonEvent->cosphimetele = -99.;
206 bendavid 1.1 fDiphotonEvent->evt = GetEventHeader()->EvtNum();
207     fDiphotonEvent->run = GetEventHeader()->RunNum();
208     fDiphotonEvent->lumi = GetEventHeader()->LumiSec();
209     fDiphotonEvent->evtcat = -99;
210 bendavid 1.2 fDiphotonEvent->nobj = fPhotons->GetEntries();
211     fDiphotonEvent->pfmet = fPFMet->At(0)->Pt();
212     fDiphotonEvent->pfmetphi = fPFMet->At(0)->Phi();
213     fDiphotonEvent->pfmetx = fPFMet->At(0)->Px();
214     fDiphotonEvent->pfmety = fPFMet->At(0)->Py();
215 fabstoec 1.19 fDiphotonEvent->spfMet = _spfMet;
216 bendavid 1.1 fDiphotonEvent->masscor = -99.;
217     fDiphotonEvent->masscorerr = -99.;
218 bendavid 1.2 fDiphotonEvent->masscorele = -99.;
219     fDiphotonEvent->masscoreleerr = -99.;
220     fDiphotonEvent->ismc = GetEventHeader()->IsMC();
221    
222 bendavid 1.5 //jets
223     const Jet *jet1 = 0;
224     const Jet *jet2 = 0;
225     const Jet *jetcentral = 0;
226    
227     fDiphotonEvent->jet1pt = -99.;
228     fDiphotonEvent->jet1eta = -99.;
229     fDiphotonEvent->jet1phi = -99.;
230     fDiphotonEvent->jet1mass = -99.;
231     fDiphotonEvent->jet2pt = -99.;
232     fDiphotonEvent->jet2eta = -99.;
233     fDiphotonEvent->jet2phi = -99.;
234     fDiphotonEvent->jet2mass = -99.;
235     fDiphotonEvent->jetcentralpt = -99.;
236     fDiphotonEvent->jetcentraleta = -99.;
237     fDiphotonEvent->jetcentralphi = -99.;
238     fDiphotonEvent->jetcentralmass = -99.;
239     fDiphotonEvent->dijetpt = -99.;
240     fDiphotonEvent->dijeteta = -99.;
241     fDiphotonEvent->dijetphi = -99.;
242     fDiphotonEvent->dijetmass = -99.;
243     fDiphotonEvent->jetetaplus = -99.;
244     fDiphotonEvent->jetetaminus = -99.;
245     fDiphotonEvent->zeppenfeld = -99.;
246     fDiphotonEvent->dphidijetgg = -99.;
247    
248 fabstoec 1.19 //uncorrected jets
249     // const Jet *uncorrjet1 = 0;
250     // const Jet *uncorrjet2 = 0;
251     // const Jet *uncorrjetcentral = 0;
252    
253     // fDiphotonEvent->uncorrjet1pt = -99.;
254     // fDiphotonEvent->uncorrjet1eta = -99.;
255     // fDiphotonEvent->uncorrjet1phi = -99.;
256     // fDiphotonEvent->uncorrjet1mass = -99.;
257     // fDiphotonEvent->uncorrjet2pt = -99.;
258     // fDiphotonEvent->uncorrjet2eta = -99.;
259     // fDiphotonEvent->uncorrjet2phi = -99.;
260     // fDiphotonEvent->uncorrjet2mass = -99.;
261     // fDiphotonEvent->uncorrjetcentralpt = -99.;
262     // fDiphotonEvent->uncorrjetcentraleta = -99.;
263     // fDiphotonEvent->uncorrjetcentralphi = -99.;
264     // fDiphotonEvent->uncorrjetcentralmass = -99.;
265     // fDiphotonEvent->diuncorrjetpt = -99.;
266     // fDiphotonEvent->diuncorrjeteta = -99.;
267     // fDiphotonEvent->diuncorrjetphi = -99.;
268     // fDiphotonEvent->diuncorrjetmass = -99.;
269    
270    
271    
272 bendavid 1.1 Int_t nhitsbeforevtxmax = 1;
273 paus 1.12 if (!fApplyElectronVeto)
274     nhitsbeforevtxmax = 999;
275 bendavid 1.1
276 bendavid 1.2 if (egcol->GetEntries()>=2) {
277     const Particle *p1 = 0;
278     const Particle *p2 = 0;
279     const Photon *phHard = 0;
280     const Photon *phSoft = 0;
281     const Electron *ele1 = 0;
282     const Electron *ele2 = 0;
283     const SuperCluster *sc1 = 0;
284     const SuperCluster *sc2 = 0;
285    
286     if (fLoopOnGoodElectrons) {
287     ele1 = fGoodElectrons->At(0);
288     ele2 = fGoodElectrons->At(1);
289     p1 = ele1;
290     p2 = ele2;
291     sc1 = ele1->SCluster();
292     sc2 = ele2->SCluster();
293     phHard = PhotonTools::MatchedPhoton(ele1,fPhotons);
294     phSoft = PhotonTools::MatchedPhoton(ele2,fPhotons);
295     }
296     else {
297     phHard = fPhotons->At(0);
298     phSoft = fPhotons->At(1);
299     p1 = phHard;
300     p2 = phSoft;
301     sc1 = phHard->SCluster();
302     sc2 = phSoft->SCluster();
303     ele1 = PhotonTools::MatchedElectron(phHard,fGoodElectrons);
304     ele2 = PhotonTools::MatchedElectron(phSoft,fGoodElectrons);
305 bendavid 1.1 }
306    
307 paus 1.12 const DecayParticle *conv1 = PhotonTools::MatchedConversion(sc1,fConversions,bsp,
308     nhitsbeforevtxmax);
309     const DecayParticle *conv2 = PhotonTools::MatchedConversion(sc2,fConversions,bsp,
310     nhitsbeforevtxmax);
311 bendavid 1.2
312 bendavid 1.20 const SuperCluster *pfsc1 = PhotonTools::MatchedPFSC(sc1,fPFPhotons, fElectrons);
313     const SuperCluster *pfsc2 = PhotonTools::MatchedPFSC(sc2,fPFPhotons, fElectrons);
314 bendavid 1.1
315     const MCParticle *phgen1 = 0;
316     const MCParticle *phgen2 = 0;
317 paus 1.12 if (!fIsData) {
318 bendavid 1.6 phgen1 = PhotonTools::MatchMC(p1,fMCParticles,!fApplyElectronVeto);
319     phgen2 = PhotonTools::MatchMC(p2,fMCParticles,!fApplyElectronVeto);
320 bendavid 1.2 }
321    
322 paus 1.12 if (fExcludeSinglePrompt && (phgen1 || phgen2))
323     return;
324     if (fExcludeDoublePrompt && (phgen1 && phgen2))
325     return;
326 bendavid 1.2
327     if (!fLoopOnGoodElectrons && phHard->HasPV()) {
328     fDiphotonEvent->vtxX = phHard->PV()->X();
329     fDiphotonEvent->vtxY = phHard->PV()->Y();
330     fDiphotonEvent->vtxZ = phHard->PV()->Z();
331 bendavid 1.4 fDiphotonEvent->vtxprob = phHard->VtxProb();
332 bendavid 1.2 }
333    
334 fabstoec 1.15 // fill Btag information... set to true if One jet fullfills
335 fabstoec 1.19 fDiphotonEvent -> btagJet1 = -1.;
336     fDiphotonEvent -> btagJet1Pt = -99.;
337     fDiphotonEvent -> btagJet1Eta = -99.;
338    
339     fDiphotonEvent -> btagJet2 = -1.;
340     fDiphotonEvent -> btagJet2Pt = -99.;
341     fDiphotonEvent -> btagJet2Eta = -99.;
342 fabstoec 1.18
343 fabstoec 1.15 if( fApplyBTag && fEnableJets ) {
344     float highTag = 0.;
345     float highJetPt = 0.;
346     float highJetEta = 0.;
347 fabstoec 1.19
348     float trailTag = 0.;
349 fabstoec 1.18 float trailJetPt = 0.;
350     float trailJetEta = 0.;
351 fabstoec 1.19
352 fabstoec 1.15 for (UInt_t ijet=0; ijet<fPFJets->GetEntries();++ijet) {
353     const Jet *jet = fPFJets->At(ijet);
354 fabstoec 1.24 if ( jet->Pt() < 20. || jet->AbsEta() > 2.4 ) continue;
355 fabstoec 1.15 if ( jet->CombinedSecondaryVertexBJetTagsDisc() > highTag ) {
356 fabstoec 1.19
357 fabstoec 1.24 trailTag = highTag;
358     trailJetPt = highJetPt;
359 fabstoec 1.18 trailJetEta = highJetEta;
360 fabstoec 1.19
361 fabstoec 1.15 highTag = jet->CombinedSecondaryVertexBJetTagsDisc();
362     highJetPt = jet->Pt();
363     highJetEta = jet->Eta();
364 fabstoec 1.19
365 fabstoec 1.18 } else if ( jet->CombinedSecondaryVertexBJetTagsDisc() > trailTag ) {
366 fabstoec 1.19
367 fabstoec 1.18 trailTag = jet->CombinedSecondaryVertexBJetTagsDisc();
368     trailJetPt = jet->Pt();
369     trailJetEta = jet->Eta();
370 fabstoec 1.15 }
371     }
372 fabstoec 1.18 fDiphotonEvent -> btagJet1 = highTag;
373     fDiphotonEvent -> btagJet1Pt = highJetPt;
374 fabstoec 1.19 fDiphotonEvent -> btagJet1Eta = highJetEta;
375    
376 fabstoec 1.24 fDiphotonEvent -> btagJet2 = trailTag;
377     fDiphotonEvent -> btagJet2Pt = trailJetPt;
378     fDiphotonEvent -> btagJet2Eta = trailJetEta;
379 fabstoec 1.19
380 fabstoec 1.15 }
381    
382 bendavid 1.5 //fill jet variables
383 bendavid 1.22 const Vertex *selvtx = fPV->At(0);
384     if (!fLoopOnGoodElectrons && phHard->HasPV()) selvtx = phHard->PV();
385 bendavid 1.5 if (fEnableJets) {
386     for (UInt_t ijet=0; ijet<fPFJets->GetEntries();++ijet) {
387     const Jet *jet = fPFJets->At(ijet);
388 bendavid 1.8 if (jet->AbsEta()<4.7 && MathUtils::DeltaR(jet,p1)>0.5 && MathUtils::DeltaR(jet,p2)>0.5) {
389 bendavid 1.22 const PFJet *pfjet = dynamic_cast<const PFJet*>(jet);
390     if (!pfjet) continue;
391     if (fApplyJetId && !fJetId.passCut(pfjet,selvtx,fPV)) continue;
392 bendavid 1.5 if (!jet1) jet1 = jet;
393     else if (!jet2) jet2 = jet;
394     else if (!jetcentral && 0) jetcentral = jet;
395     }
396     if (jet1&&jet2&&jetcentral) break;
397 mtouch 1.25 }
398 bendavid 1.5 }
399    
400     if (jet1) {
401     fDiphotonEvent->jet1pt = jet1->Pt();
402     fDiphotonEvent->jet1eta = jet1->Eta();
403     fDiphotonEvent->jet1phi = jet1->Phi();
404     fDiphotonEvent->jet1mass = jet1->Mass();
405     }
406    
407     if (jet2) {
408     fDiphotonEvent->jet2pt = jet2->Pt();
409     fDiphotonEvent->jet2eta = jet2->Eta();
410     fDiphotonEvent->jet2phi = jet2->Phi();
411     fDiphotonEvent->jet2mass = jet2->Mass();
412     }
413    
414     if (jetcentral) {
415     fDiphotonEvent->jetcentralpt = jetcentral->Pt();
416     fDiphotonEvent->jetcentraleta = jetcentral->Eta();
417     fDiphotonEvent->jetcentralphi = jetcentral->Phi();
418     fDiphotonEvent->jetcentralmass = jetcentral->Mass();
419     }
420    
421     if (jet1&&jet2){
422     FourVectorM momjj = (jet1->Mom() + jet2->Mom());
423    
424     fDiphotonEvent->dijetpt = momjj.Pt();
425     fDiphotonEvent->dijeteta = momjj.Eta();
426     fDiphotonEvent->dijetphi = momjj.Phi();
427     fDiphotonEvent->dijetmass = momjj.M();
428    
429     if (jet1->Eta()>jet2->Eta()) {
430     fDiphotonEvent->jetetaplus = jet1->Eta();
431     fDiphotonEvent->jetetaminus = jet2->Eta();
432     }
433     else {
434     fDiphotonEvent->jetetaplus = jet2->Eta();
435     fDiphotonEvent->jetetaminus = jet1->Eta();
436     }
437     }
438    
439 mtouch 1.25
440    
441    
442    
443    
444 fabstoec 1.16 //added gen. info of whether a lep. or nutrino is from W or Z --Heng 02/14/2012 12:30 EST
445     Double_t _fromZ = -99;
446     Double_t _fromW = -99;
447     Float_t _zpt = -99;
448     Float_t _zEta = -99;
449     Float_t _allZpt = -99;
450     Float_t _allZEta = -99;
451    
452     if( !fIsData ){
453    
454     // loop over all GEN particles and look for nutrinos whoes mother is Z
455     for(UInt_t j=0; j<fMCParticles->GetEntries(); ++j) {
456     const MCParticle* p = fMCParticles->At(j);
457     if( p->AbsPdgId()==23 ||p->AbsPdgId()==32 || p->AbsPdgId()==33 ) {
458     _allZpt=p->Pt();
459     _allZEta=p->Eta();
460     if (p->HasDaughter(12,kFALSE) || p->HasDaughter(14,kFALSE) || p->HasDaughter(16,kFALSE) ||p->HasDaughter(18,kFALSE) ) {
461     _fromZ=1;
462     _zpt=p->Pt();
463     _zEta=p->Eta();
464     }
465     }
466     else _fromW=1;
467     }
468     }
469 mtouch 1.25
470 fabstoec 1.16 fDiphotonEvent->fromZ = _fromZ;
471     fDiphotonEvent->fromW = _fromW;
472     fDiphotonEvent->zpt = _zpt;
473     fDiphotonEvent->zEta = _zEta;
474     fDiphotonEvent->allZpt = _allZpt;
475     fDiphotonEvent->allZEta = _allZEta;
476    
477     Double_t _dphiMetgg = -99;
478     Double_t _cosdphiMetgg = -99;
479     Double_t _dphiPhPh = -99;
480 mtouch 1.25
481 bendavid 1.4 Double_t _mass = -99.;
482     Double_t _masserr = -99.;
483     Double_t _masserrsmeared = -99.;
484     Double_t _masserrwrongvtx = -99.;
485     Double_t _masserrsmearedwrongvtx = -99.;
486     Double_t _ptgg = -99.;
487 bendavid 1.5 Double_t _etagg = -99.;
488     Double_t _phigg = -99.;
489 bendavid 1.4 Double_t _costheta = -99.;
490 bendavid 1.2 PhotonTools::DiphotonR9EtaPtCats _evtcat = PhotonTools::kOctCat0;
491     if (phHard && phSoft) {
492 fabstoec 1.16 _dphiMetgg = MathUtils::DeltaPhi((phHard->Mom()+phSoft->Mom()).Phi(),fPFMet->At(0)->Phi());
493     _cosdphiMetgg = TMath::Cos(_dphiMetgg);
494     _dphiPhPh = MathUtils::DeltaPhi((phHard->Mom()).Phi(),(phSoft->Mom()).Phi());
495 bendavid 1.2 _mass = (phHard->Mom()+phSoft->Mom()).M();
496     _masserr = 0.5*_mass*TMath::Sqrt(phHard->EnergyErr()*phHard->EnergyErr()/phHard->E()/phHard->E() + phSoft->EnergyErr()*phSoft->EnergyErr()/phSoft->E()/phSoft->E());
497     _masserrsmeared = 0.5*_mass*TMath::Sqrt(phHard->EnergyErrSmeared()*phHard->EnergyErrSmeared()/phHard->E()/phHard->E() + phSoft->EnergyErrSmeared()*phSoft->EnergyErrSmeared()/phSoft->E()/phSoft->E());
498     _ptgg = (phHard->Mom()+phSoft->Mom()).Pt();
499 bendavid 1.5 _etagg = (phHard->Mom()+phSoft->Mom()).Eta();
500     _phigg = (phHard->Mom()+phSoft->Mom()).Phi();
501 bendavid 1.2 _costheta = ThreeVector(phHard->Mom()).Unit().Dot(ThreeVector(phSoft->Mom()).Unit());
502     _evtcat = PhotonTools::DiphotonR9EtaPtCat(phHard,phSoft);
503 bendavid 1.4
504 bendavid 1.5 const Double_t dz = sqrt(2.0)*5.8;
505 paus 1.12 Double_t deltamvtx = _mass*VertexTools::DeltaMassVtx(phHard->CaloPos().X(),
506     phHard->CaloPos().Y(),
507     phHard->CaloPos().Z(),
508     phSoft->CaloPos().X(),
509     phSoft->CaloPos().Y(),
510     phSoft->CaloPos().Z(),
511     fDiphotonEvent->vtxX,
512     fDiphotonEvent->vtxY,
513     fDiphotonEvent->vtxZ,
514     dz);
515 bendavid 1.5 fDiphotonEvent->deltamvtx = deltamvtx;
516    
517 paus 1.12 _masserrwrongvtx = TMath::Sqrt(_masserr*_masserr + deltamvtx*deltamvtx);
518 bendavid 1.4 _masserrsmearedwrongvtx = TMath::Sqrt(_masserrsmeared*_masserrsmeared + deltamvtx*deltamvtx);
519    
520 bendavid 1.5 if (jet1 && jet2) {
521     fDiphotonEvent->zeppenfeld = TMath::Abs(_etagg - 0.5*(jet1->Eta()+jet2->Eta()));
522     fDiphotonEvent->dphidijetgg = MathUtils::DeltaPhi( (jet1->Mom()+jet2->Mom()).Phi(), _phigg );
523     }
524 bendavid 1.4
525 bendavid 1.20 PFJetOArr pfjets;
526     for (UInt_t ijet=0; ijet<fPFJets->GetEntries(); ++ijet) {
527     const PFJet *pfjet = dynamic_cast<const PFJet*>(fPFJets->At(ijet));
528     if (pfjet && MathUtils::DeltaR(*pfjet,*phHard)>0.3 && MathUtils::DeltaR(*pfjet,*phSoft)>0.3) pfjets.Add(pfjet);
529     }
530    
531     PFCandidateOArr pfcands;
532     for (UInt_t icand=0; icand<fPFCands->GetEntries(); ++icand) {
533     const PFCandidate *pfcand = fPFCands->At(icand);
534     if (MathUtils::DeltaR(*pfcand,*phHard)>0.1 && MathUtils::DeltaR(*pfcand,*phSoft)>0.1) pfcands.Add(pfcand);
535     }
536    
537     const Vertex *firstvtx = fPV->At(0);
538     const Vertex *selvtx = fPV->At(0);
539    
540     if (!fLoopOnGoodElectrons && phHard->HasPV()) {
541     selvtx = phHard->PV();
542     }
543    
544 bendavid 1.21 if (0) //disable for now for performance reasons
545 bendavid 1.20 {
546     Met mmet = fMVAMet.GetMet( false,
547     0.,0.,0.,
548     fPFMet->At(0),
549     &pfcands,selvtx,fPV, fPileUpDen->At(0)->Rho(),
550     &pfjets,
551     int(fPV->GetEntries()),
552     kFALSE);
553    
554     TMatrixD *metcov = fMVAMet.GetMetCovariance();
555    
556     ThreeVector fullmet(mmet.Px() - phHard->Px() - phSoft->Px(),
557     mmet.Py() - phHard->Py() - phSoft->Py(),
558     0.);
559    
560     ROOT::Math::SMatrix<double,2,2,ROOT::Math::MatRepSym<double,2> > mcov;
561     mcov(0,0) = (*metcov)(0,0);
562     mcov(0,1) = (*metcov)(0,1);
563     mcov(1,0) = (*metcov)(1,0);
564     mcov(1,1) = (*metcov)(1,1);
565     ROOT::Math::SVector<double,2> vmet;
566     vmet(0) = fullmet.X();
567     vmet(1) = fullmet.Y();
568     mcov.Invert();
569     Double_t metsig = sqrt(ROOT::Math::Similarity(mcov,vmet));
570    
571     fDiphotonEvent->mvametsel = fullmet.Rho();
572     fDiphotonEvent->mvametselphi = fullmet.Phi();
573     fDiphotonEvent->mvametselx = fullmet.X();
574     fDiphotonEvent->mvametsely = fullmet.Y();
575     fDiphotonEvent->mvametselsig = metsig;
576     }
577    
578 bendavid 1.21 if (0) //disable for now for performance reasons
579 bendavid 1.20 {
580     Met mmet = fMVAMet.GetMet( false,
581     0.,0.,0.,
582     fPFMet->At(0),
583     &pfcands,firstvtx,fPV, fPileUpDen->At(0)->Rho(),
584     &pfjets,
585     int(fPV->GetEntries()),
586     kFALSE);
587    
588     TMatrixD *metcov = fMVAMet.GetMetCovariance();
589    
590     ThreeVector fullmet(mmet.Px() - phHard->Px() - phSoft->Px(),
591     mmet.Py() - phHard->Py() - phSoft->Py(),
592     0.);
593    
594     ROOT::Math::SMatrix<double,2,2,ROOT::Math::MatRepSym<double,2> > mcov;
595     mcov(0,0) = (*metcov)(0,0);
596     mcov(0,1) = (*metcov)(0,1);
597     mcov(1,0) = (*metcov)(1,0);
598     mcov(1,1) = (*metcov)(1,1);
599     ROOT::Math::SVector<double,2> vmet;
600     vmet(0) = fullmet.X();
601     vmet(1) = fullmet.Y();
602     mcov.Invert();
603     Double_t metsig = sqrt(ROOT::Math::Similarity(mcov,vmet));
604    
605     fDiphotonEvent->mvametfirst = fullmet.Rho();
606     fDiphotonEvent->mvametfirstphi = fullmet.Phi();
607     fDiphotonEvent->mvametfirstx = fullmet.X();
608     fDiphotonEvent->mvametfirsty = fullmet.Y();
609     fDiphotonEvent->mvametfirstsig = metsig;
610     }
611    
612 bendavid 1.2 }
613 fabstoec 1.19
614 mtouch 1.25
615 fabstoec 1.19 fDiphotonEvent->corrpfmet = -99.;
616     fDiphotonEvent->corrpfmetphi = -99.;
617     fDiphotonEvent->corrpfmetx = -99.;
618     fDiphotonEvent->corrpfmety = -99.;
619    
620     Met *corrMet =NULL;
621 bendavid 1.2
622 fabstoec 1.19 if (fApplyPFMetCorrections){
623     corrMet = new Met(fPFMet->At(0)->Px(),fPFMet->At(0)->Py());
624 fabstoec 1.17
625 fabstoec 1.19 if (!fIsData){
626 mtouch 1.25 PFMetCorrectionTools::correctMet(corrMet,phHard,phSoft,1,0,funcorrPFJets,fGenJets,fPFJets,_evt);
627 fabstoec 1.19 PFMetCorrectionTools::shiftMet(corrMet,fIsData,_spfMet);
628     }
629     else {
630     PFMetCorrectionTools::shiftMet(corrMet,fIsData,_spfMet);
631 mtouch 1.25 PFMetCorrectionTools::correctMet(corrMet,phHard,phSoft,0,1,funcorrPFJets,fGenJets,fPFJets,_evt);
632 fabstoec 1.19 }
633    
634     fDiphotonEvent->corrpfmet = corrMet->Pt();
635     fDiphotonEvent->corrpfmetphi = corrMet->Phi();
636     fDiphotonEvent->corrpfmetx = corrMet->Px();
637     fDiphotonEvent->corrpfmety = corrMet->Py();
638 fabstoec 1.17
639 fabstoec 1.19 delete corrMet;
640 fabstoec 1.17 }
641 fabstoec 1.16
642 mtouch 1.25
643 bendavid 1.2 Float_t _massele = -99.;
644     Float_t _ptee = -99.;
645     Float_t _costhetaele = -99.;
646     if (ele1 && ele2) {
647     _massele = (ele1->Mom()+ele2->Mom()).M();
648     _ptee = (ele1->Mom()+ele2->Mom()).Pt();
649     _costhetaele = ThreeVector(ele1->Mom()).Unit().Dot(ThreeVector(ele2->Mom()).Unit());
650 bendavid 1.1 }
651    
652     Float_t _gencostheta = -99.;
653     if (phgen1 && phgen2) {
654     _gencostheta = ThreeVector(phgen1->Mom()).Unit().Dot(ThreeVector(phgen2->Mom()).Unit());
655     }
656    
657     fDiphotonEvent->gencostheta = _gencostheta;
658 fabstoec 1.16 fDiphotonEvent->dphiMetgg = _dphiMetgg;
659     fDiphotonEvent->cosdphiMetgg = _cosdphiMetgg;
660     fDiphotonEvent->dphiPhPh = _dphiPhPh;
661 bendavid 1.1 fDiphotonEvent->mass = _mass;
662 bendavid 1.2 fDiphotonEvent->masserr = _masserr;
663     fDiphotonEvent->masserrsmeared = _masserrsmeared;
664 bendavid 1.4 fDiphotonEvent->masserrwrongvtx = _masserrwrongvtx;
665     fDiphotonEvent->masserrsmearedwrongvtx = _masserrsmearedwrongvtx;
666 bendavid 1.1 fDiphotonEvent->ptgg = _ptgg;
667 bendavid 1.5 fDiphotonEvent->etagg = _etagg;
668     fDiphotonEvent->phigg = _phigg;
669 bendavid 1.2 fDiphotonEvent->costheta = _costheta;;
670     fDiphotonEvent->massele = _massele;
671     fDiphotonEvent->ptee = _ptee;
672     fDiphotonEvent->costhetaele = _costhetaele;
673     fDiphotonEvent->evtcat = _evtcat;
674 bendavid 1.1
675 fabstoec 1.24 fDiphotonEvent->photons[0].SetVars(phHard,conv1,ele1,pfsc1,phgen1,fPhfixph,fPhfixele,fTracks,fPV,fPFCands,rho,fFillClusterArrays,fElectrons,fApplyElectronVeto);
676     fDiphotonEvent->photons[1].SetVars(phSoft,conv2,ele2,pfsc2,phgen2,fPhfixph,fPhfixele,fTracks,fPV,fPFCands,rho,fFillClusterArrays,fElectrons,fApplyElectronVeto);
677 bendavid 1.1
678     Float_t ph1ecor = fDiphotonEvent->photons[0].Ecor();
679     Float_t ph1ecorerr = fDiphotonEvent->photons[0].Ecorerr();
680     Float_t ph2ecor = fDiphotonEvent->photons[1].Ecor();
681     Float_t ph2ecorerr = fDiphotonEvent->photons[1].Ecorerr();
682 bendavid 1.2
683     Float_t ph1ecorele = fDiphotonEvent->photons[0].Ecorele();
684     Float_t ph1ecoreleerr = fDiphotonEvent->photons[0].Ecoreleerr();
685     Float_t ph2ecorele = fDiphotonEvent->photons[1].Ecorele();
686     Float_t ph2ecoreleerr = fDiphotonEvent->photons[1].Ecoreleerr();
687    
688 bendavid 1.1
689 paus 1.12 fDiphotonEvent->masscor = TMath::Sqrt(2.0*ph1ecor*ph2ecor*(1.0-fDiphotonEvent->costheta));
690     fDiphotonEvent->masscorerr = 0.5*fDiphotonEvent->masscor*
691     TMath::Sqrt(ph1ecorerr*ph1ecorerr/ph1ecor/ph1ecor + ph2ecorerr*ph2ecorerr/ph2ecor/ph2ecor);
692    
693     fDiphotonEvent->masscorele = TMath::Sqrt(2.0*ph1ecorele*ph2ecorele*
694     (1.0-fDiphotonEvent->costheta));
695     fDiphotonEvent->masscoreleerr = 0.5*fDiphotonEvent->masscorele*
696     TMath::Sqrt(ph1ecoreleerr*ph1ecoreleerr/ph1ecorele/ph1ecorele +
697     ph2ecoreleerr*ph2ecoreleerr/ph2ecorele/ph2ecorele);
698 bendavid 1.2
699 paus 1.12 //printf("r9 = %5f, photon sigieie = %5f, seed sigieie = %5f\n",phHard->R9(),
700     // phHard->CoviEtaiEta(),sqrt(phHard->SCluster()->Seed()->CoviEtaiEta()));
701 fabstoec 1.13
702 fabstoec 1.14 // MuonStuff
703     fDiphotonEvent-> muonPt = -99.;
704     fDiphotonEvent-> muonEta = -99.;
705     fDiphotonEvent-> muDR1 = -99.;
706     fDiphotonEvent-> muDR2 = -99.;
707     fDiphotonEvent-> muIso1 = -99.;
708     fDiphotonEvent-> muIso2 = -99.;
709     fDiphotonEvent-> muIso3 = -99.;
710     fDiphotonEvent-> muIso4 = -99.;
711     fDiphotonEvent-> muD0 = -99.;
712     fDiphotonEvent-> muDZ = -99.;
713     fDiphotonEvent-> muChi2 = -99.;
714     fDiphotonEvent-> muNhits = -99;
715     fDiphotonEvent-> muNpixhits = -99;
716     fDiphotonEvent-> muNegs = -99;
717     fDiphotonEvent-> muNMatch = -99;
718    
719     // Electron Stuff
720     fDiphotonEvent-> elePt = -99.;
721     fDiphotonEvent-> eleEta = -99.;
722     fDiphotonEvent-> eleSCEta = -99.;
723     fDiphotonEvent-> eleIso1 = -99.;
724     fDiphotonEvent-> eleIso2 = -99.;
725     fDiphotonEvent-> eleIso3 = -99.;
726     fDiphotonEvent-> eleIso4 = -99.;
727     fDiphotonEvent-> eleDist = -99.;
728     fDiphotonEvent-> eleDcot = -99.;
729     fDiphotonEvent-> eleCoviee = -99.;
730     fDiphotonEvent-> eleDphiin = -99.;
731     fDiphotonEvent-> eleDetain = -99.;
732     fDiphotonEvent-> eleDR1 = -99.;
733     fDiphotonEvent-> eleDR2 = -99.;
734     fDiphotonEvent-> eleMass1 = -99.;
735     fDiphotonEvent-> eleMass2 = -99.;
736     fDiphotonEvent-> eleNinnerHits = -99;
737    
738 fabstoec 1.13 if( fApplyLeptonTag ) {
739     // perform lepton tagging
740     // the diphoton event record will have one more entry; i.e. leptonTag
741     // leptonTag = -1 -> lepton-taggng was swicthed off
742     // = 0 -> event tagged as 'non-lepton-event'
743     // = +1 -> event tagged as muon-event
744     // = +2 -> event tagged as electron-event
745     fDiphotonEvent->leptonTag = 0;
746    
747     if ( fLeptonTagMuons->GetEntries() > 0 ) {
748 fabstoec 1.14
749 fabstoec 1.13 // need to have dR > 1 for with respect to both photons
750     for(UInt_t iMuon = 0; iMuon <fLeptonTagMuons->GetEntries(); ++iMuon) {
751     if(MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phHard) < 1.) continue;
752     if(MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phSoft) < 1.) continue;
753    
754     fDiphotonEvent->leptonTag = 1;
755 fabstoec 1.14
756     fDiphotonEvent-> muonPt = fLeptonTagMuons->At(iMuon)->Pt();
757     fDiphotonEvent-> muonEta = fLeptonTagMuons->At(iMuon)->Eta();
758     fDiphotonEvent-> muDR1 = MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phHard);
759     fDiphotonEvent-> muDR2 = MathUtils::DeltaR(fLeptonTagMuons->At(iMuon),phSoft);
760     fDiphotonEvent-> muIso1 = (fLeptonTagMuons->At(iMuon)->IsoR03SumPt() + fLeptonTagMuons->At(iMuon)->IsoR03EmEt() + fLeptonTagMuons->At(iMuon)->IsoR03HadEt() - fPileUpDen->At(0)->RhoRandomLowEta() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(iMuon)->Pt();
761     fDiphotonEvent-> muIso2 = (fLeptonTagMuons->At(iMuon)->IsoR03SumPt() + fLeptonTagMuons->At(iMuon)->IsoR03EmEt() + fLeptonTagMuons->At(iMuon)->IsoR03HadEt() - fPileUpDen->At(0)->RhoRandom() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(iMuon)->Pt();
762     fDiphotonEvent-> muIso3 = (fLeptonTagMuons->At(iMuon)->IsoR03SumPt() + fLeptonTagMuons->At(iMuon)->IsoR03EmEt() + fLeptonTagMuons->At(iMuon)->IsoR03HadEt() - fPileUpDen->At(0)->RhoLowEta() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(iMuon)->Pt();
763     fDiphotonEvent-> muIso4 = (fLeptonTagMuons->At(iMuon)->IsoR03SumPt() + fLeptonTagMuons->At(iMuon)->IsoR03EmEt() + fLeptonTagMuons->At(iMuon)->IsoR03HadEt() - fPileUpDen->At(0)->Rho() * TMath::Pi() * 0.3 * 0.3)/ fLeptonTagMuons->At(iMuon)->Pt();
764     fDiphotonEvent-> muD0 = TMath::Abs(fLeptonTagMuons->At(iMuon)->BestTrk()->D0Corrected(*fPV->At(0)));
765     fDiphotonEvent-> muDZ = TMath::Abs(fLeptonTagMuons->At(iMuon)->BestTrk()->DzCorrected(*fPV->At(0)));
766     fDiphotonEvent-> muChi2 = fLeptonTagMuons->At(iMuon)->GlobalTrk()->Chi2()/fLeptonTagMuons->At(iMuon)->GlobalTrk()->Ndof();
767    
768     fDiphotonEvent-> muNhits = fLeptonTagMuons->At(iMuon)->BestTrk()->NHits();
769     fDiphotonEvent-> muNpixhits = fLeptonTagMuons->At(iMuon)->BestTrk()->NPixelHits();
770     fDiphotonEvent-> muNegs = fLeptonTagMuons->At(iMuon)->NSegments();
771     fDiphotonEvent-> muNMatch = fLeptonTagMuons->At(iMuon)->NMatches();
772    
773 fabstoec 1.13 break;
774     }
775     }
776     if ( fDiphotonEvent->leptonTag < 1 && fLeptonTagElectrons->GetEntries() > 0 ) {
777     for(UInt_t iEle = 0; iEle < fLeptonTagElectrons->GetEntries(); ++iEle) {
778     if(MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phHard) < 1.) continue;
779     if(MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phSoft) < 1.) continue;
780    
781     // here we also check the mass ....
782     if ( TMath::Abs( (phHard->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M()-91.19 ) < 5. ) continue;
783     if ( TMath::Abs( (phSoft->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M()-91.19 ) < 5. ) continue;
784    
785     fDiphotonEvent->leptonTag = 2;
786 fabstoec 1.14
787     fDiphotonEvent-> elePt = fLeptonTagElectrons->At(iEle)->Pt();
788     fDiphotonEvent-> eleEta = fLeptonTagElectrons->At(iEle)->Eta();
789     fDiphotonEvent-> eleSCEta = fLeptonTagElectrons->At(iEle)->SCluster()->Eta();
790     fDiphotonEvent-> eleIso1 = (fLeptonTagElectrons->At(iEle)->TrackIsolationDr03() + fLeptonTagElectrons->At(iEle)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(iEle)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->RhoRandomLowEta() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
791     fDiphotonEvent-> eleIso2 = (fLeptonTagElectrons->At(iEle)->TrackIsolationDr03() + fLeptonTagElectrons->At(iEle)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(iEle)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->RhoRandom() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
792     fDiphotonEvent-> eleIso3 = (fLeptonTagElectrons->At(iEle)->TrackIsolationDr03() + fLeptonTagElectrons->At(iEle)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(iEle)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->RhoLowEta() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
793     fDiphotonEvent-> eleIso4 = (fLeptonTagElectrons->At(iEle)->TrackIsolationDr03() + fLeptonTagElectrons->At(iEle)->EcalRecHitIsoDr03() + fLeptonTagElectrons->At(iEle)->HcalTowerSumEtDr03() - fPileUpDen->At(0)->Rho() * TMath::Pi() * 0.3 * 0.3)/fDiphotonEvent-> elePt;
794     fDiphotonEvent-> eleDist = fLeptonTagElectrons->At(iEle)->ConvPartnerDist();
795     fDiphotonEvent-> eleDcot = fLeptonTagElectrons->At(iEle)->ConvPartnerDCotTheta();
796     fDiphotonEvent-> eleCoviee = fLeptonTagElectrons->At(iEle)->CoviEtaiEta();
797     fDiphotonEvent-> eleDphiin = TMath::Abs(fLeptonTagElectrons->At(iEle)->DeltaPhiSuperClusterTrackAtVtx());
798     fDiphotonEvent-> eleDetain = TMath::Abs(fLeptonTagElectrons->At(iEle)->DeltaEtaSuperClusterTrackAtVtx());
799     fDiphotonEvent-> eleDR1 = MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phHard);
800     fDiphotonEvent-> eleDR2 = MathUtils::DeltaR(fLeptonTagElectrons->At(iEle),phSoft);
801     fDiphotonEvent-> eleMass1 = (phHard->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M();
802     fDiphotonEvent-> eleMass2 = (phSoft->Mom()+fLeptonTagElectrons->At(iEle)->Mom()).M();
803     fDiphotonEvent-> eleNinnerHits = fLeptonTagElectrons->At(iEle)->Trk()->NExpectedHitsInner();
804    
805 fabstoec 1.13 break;
806     }
807     }
808     }
809    
810 paus 1.12 if (fWriteDiphotonTree)
811     hCiCTuple->Fill();
812 bendavid 1.22
813     if (fFillVertexTree && phHard && phSoft) {
814     for (UInt_t ivtx = 0; ivtx<fPV->GetEntries(); ++ivtx) {
815     const Vertex *v = fPV->At(ivtx);
816     fDiphotonVtx->SetVars(v,phHard,phSoft,fPFCands,ivtx,fPV->GetEntries(),_decayZ);
817     hVtxTree->Fill();
818     }
819     }
820    
821 bendavid 1.1 }
822    
823 paus 1.12 if (!fWriteSingleTree)
824     return;
825 bendavid 1.1
826 bendavid 1.2 for (UInt_t iph = 0; iph<egcol->GetEntries(); ++iph) {
827     const Particle *p = 0;
828     const Photon *ph = 0;
829     const Electron *ele = 0;
830     const SuperCluster *sc = 0;
831     if (fLoopOnGoodElectrons) {
832     ele = fGoodElectrons->At(iph);
833     p = ele;
834     sc = ele->SCluster();
835     ph = PhotonTools::MatchedPhoton(ele,fPhotons);
836     }
837     else {
838     ph = fPhotons->At(iph);
839     p = ph;
840     sc = ph->SCluster();
841     ele = PhotonTools::MatchedElectron(ph,fGoodElectrons);
842     }
843    
844 paus 1.12 const DecayParticle *conv = PhotonTools::MatchedConversion(sc,fConversions,bsp,
845     nhitsbeforevtxmax);
846 bendavid 1.20 const SuperCluster *pfsc = PhotonTools::MatchedPFSC(sc,fPFPhotons, fElectrons);
847    
848 bendavid 1.2 if (!fLoopOnGoodElectrons && ph->HasPV()) {
849 bendavid 1.1 fDiphotonEvent->vtxZ = ph->PV()->Z();
850     }
851    
852     const MCParticle *phgen = 0;
853     if( !fIsData ) {
854 bendavid 1.6 phgen = PhotonTools::MatchMC(p,fMCParticles,!fApplyElectronVeto);
855 bendavid 1.1 }
856    
857 bendavid 1.2 if (fExcludeSinglePrompt && phgen) return;
858    
859     fDiphotonEvent->mt = -99.;
860     fDiphotonEvent->cosphimet = -99.;
861     fDiphotonEvent->mtele = -99.;
862     fDiphotonEvent->cosphimetele = -99.;
863    
864     if (ph) {
865     fDiphotonEvent->cosphimet = TMath::Cos(ph->Phi()-fPFMet->At(0)->Phi());
866 paus 1.12 fDiphotonEvent->mt = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ph->Pt()*
867     (1.0-fDiphotonEvent->cosphimet));
868 bendavid 1.2 }
869 bendavid 1.1
870 bendavid 1.2 if (ele) {
871     fDiphotonEvent->cosphimetele = TMath::Cos(ele->Phi()-fPFMet->At(0)->Phi());
872 paus 1.12 fDiphotonEvent->mtele = TMath::Sqrt(2.0*fPFMet->At(0)->Pt()*ele->Pt()*
873     (1.0-fDiphotonEvent->cosphimetele));
874 bendavid 1.2 }
875    
876 fabstoec 1.24 fSinglePhoton->SetVars(ph,conv,ele,pfsc,phgen,fPhfixph,fPhfixele,fTracks,fPV,fPFCands,rho,fFillClusterArrays,
877 paus 1.12 fElectrons,fApplyElectronVeto);
878 bendavid 1.1 hCiCTupleSingle->Fill();
879     }
880 fabstoec 1.16
881    
882    
883 bendavid 1.1
884     return;
885     }
886    
887     //--------------------------------------------------------------------------------------------------
888     void PhotonTreeWriter::SlaveBegin()
889     {
890     // Run startup code on the computer (slave) doing the actual analysis. Here,
891     // we just request the photon collection branch.
892    
893 fabstoec 1.13 if( fApplyLeptonTag ) {
894     ReqEventObject(fLeptonTagElectronsName, fLeptonTagElectrons, false);
895     ReqEventObject(fLeptonTagMuonsName, fLeptonTagMuons, false);
896     }
897    
898 paus 1.12 ReqEventObject(fPhotonBranchName,fPhotons, fPhotonsFromBranch);
899 bendavid 1.20 if (fEnablePFPhotons) ReqEventObject(fPFPhotonName,fPFPhotons, true);
900 paus 1.12 ReqEventObject(fTrackBranchName, fTracks, true);
901     ReqEventObject(fElectronName, fElectrons, true);
902     ReqEventObject(fGoodElectronName,fGoodElectrons,fGoodElectronsFromBranch);
903     ReqEventObject(fPileUpDenName, fPileUpDen, true);
904     ReqEventObject(fPVName, fPV, fPVFromBranch);
905     ReqEventObject(fConversionName, fConversions, true);
906     ReqEventObject(fBeamspotName, fBeamspot, true);
907     ReqEventObject(fPFCandName, fPFCands, true);
908     ReqEventObject(fSuperClusterName,fSuperClusters,true);
909     ReqEventObject(fPFMetName, fPFMet, true);
910 fabstoec 1.16 if (fEnableJets){
911     ReqEventObject(fPFJetName, fPFJets, fPFJetsFromBranch);
912     ReqBranch(funcorrPFJetName, funcorrPFJets);
913 fabstoec 1.19 // if (!fIsData) ReqEventObject(fGenJetName, fGenJets, true);
914 fabstoec 1.16 }
915 bendavid 1.1 if (!fIsData) {
916 paus 1.12 ReqBranch(fPileUpName, fPileUp);
917     ReqBranch(fMCParticleName, fMCParticles);
918 fabstoec 1.19 if (fEnableJets) ReqEventObject(fGenJetName, fGenJets, true);
919 bendavid 1.1 }
920 bendavid 1.6 if (fIsData) {
921 paus 1.12 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
922     TString("/src/MitPhysics/data/PhotonFixGRPV22.dat");
923 bendavid 1.6 }
924     else {
925 paus 1.12 fPhFixDataFile = gSystem->Getenv("CMSSW_BASE") +
926     TString("/src/MitPhysics/data/PhotonFixSTART42V13.dat");
927 bendavid 1.6 }
928 bendavid 1.1
929 bendavid 1.2 fPhfixph.initialise("4_2",std::string(fPhFixDataFile));
930     fPhfixele.initialise("4_2e",std::string(fPhFixDataFile));
931    
932 bendavid 1.20 // fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
933     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
934     // TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
935     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_42.root"))),
936     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_42.root"))),
937     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1_42.root"))),
938     // TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2_42.root")))
939     // );
940    
941     fMVAMet.Initialize(TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
942     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
943     TString(getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")),
944     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmet_52.root"))),
945     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetphi_52.root"))),
946     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu1cov_52.root"))),
947     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/gbrmetu2cov_52.root")))
948     );
949 bendavid 1.22
950     fJetId.Initialize(JetIDMVA::kMedium,
951 bendavid 1.23 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_lowpt.weights.xml"))),
952     TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/mva_JetID_highpt.weights.xml"))),
953 bendavid 1.22 JetIDMVA::kCut,
954 bendavid 1.23 TString((getenv("CMSSW_BASE")+string("/src/MitPhysics/Utils/python/JetIdParams_cfi.py")))
955     );
956 bendavid 1.20
957 bendavid 1.1 fDiphotonEvent = new PhotonTreeWriterDiphotonEvent;
958 bendavid 1.20 fSinglePhoton = new PhotonTreeWriterPhoton<16>;
959 bendavid 1.1
960 bendavid 1.20 TFile *ftmp = TFile::Open(TString::Format("%s_tmp.root",GetName()),"RECREATE");
961    
962     if (fWriteDiphotonTree) {
963 paus 1.12 hCiCTuple = new TTree(fTupleName.Data(),fTupleName.Data());
964 bendavid 1.20 hCiCTuple->SetAutoSave(300e9);
965     }
966 bendavid 1.1 TString singlename = fTupleName + TString("Single");
967 bendavid 1.20 if (fWriteSingleTree) {
968 paus 1.12 hCiCTupleSingle = new TTree(singlename,singlename);
969 bendavid 1.20 hCiCTupleSingle->SetAutoSave(300e9);
970     }
971    
972 bendavid 1.1 //make flattish tree from classes so we don't have to rely on dictionaries for reading later
973     TClass *eclass = TClass::GetClass("mithep::PhotonTreeWriterDiphotonEvent");
974 bendavid 1.20 TClass *pclass = TClass::GetClass("mithep::PhotonTreeWriterPhoton<16>");
975 paus 1.12 TList *elist = eclass->GetListOfDataMembers();
976     TList *plist = pclass->GetListOfDataMembers();
977 bendavid 1.1
978     for (int i=0; i<elist->GetEntries(); ++i) {
979 mingyang 1.26 const TDataMember *tdm = static_cast<const TDataMember*>(elist->At(i));//ming
980 bendavid 1.20 if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
981 bendavid 1.1 TString typestring;
982 bendavid 1.20 if (TString(tdm->GetTypeName()).BeginsWith("Char_t")) typestring = "B";
983     else if (TString(tdm->GetTypeName()).BeginsWith("UChar_t")) typestring = "b";
984     else if (TString(tdm->GetTypeName()).BeginsWith("Short_t")) typestring = "S";
985     else if (TString(tdm->GetTypeName()).BeginsWith("UShort_t")) typestring = "s";
986     else if (TString(tdm->GetTypeName()).BeginsWith("Int_t")) typestring = "I";
987     else if (TString(tdm->GetTypeName()).BeginsWith("UInt_t")) typestring = "i";
988     else if (TString(tdm->GetTypeName()).BeginsWith("Float_t")) typestring = "F";
989     else if (TString(tdm->GetTypeName()).BeginsWith("Double_t")) typestring = "D";
990     else if (TString(tdm->GetTypeName()).BeginsWith("Long64_t")) typestring = "L";
991     else if (TString(tdm->GetTypeName()).BeginsWith("ULong64_t")) typestring = "l";
992     else if (TString(tdm->GetTypeName()).BeginsWith("Bool_t")) typestring = "O";
993     else continue;
994 bendavid 1.1 //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
995 mingyang 1.26 Char_t *addr = (Char_t*)fDiphotonEvent;//ming:?
996 bendavid 1.1 assert(sizeof(Char_t)==1);
997 bendavid 1.20 if (fWriteDiphotonTree) hCiCTuple->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
998     if (fWriteSingleTree) hCiCTupleSingle->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
999 bendavid 1.1 }
1000    
1001     for (int iph=0; iph<2; ++iph) {
1002     for (int i=0; i<plist->GetEntries(); ++i) {
1003     const TDataMember *tdm = static_cast<const TDataMember*>(plist->At(i));
1004 bendavid 1.20 if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
1005 bendavid 1.1 TString typestring;
1006 bendavid 1.20 if (TString(tdm->GetTypeName()).BeginsWith("Char_t")) typestring = "B";
1007     else if (TString(tdm->GetTypeName()).BeginsWith("UChar_t")) typestring = "b";
1008     else if (TString(tdm->GetTypeName()).BeginsWith("Short_t")) typestring = "S";
1009     else if (TString(tdm->GetTypeName()).BeginsWith("UShort_t")) typestring = "s";
1010     else if (TString(tdm->GetTypeName()).BeginsWith("Int_t")) typestring = "I";
1011     else if (TString(tdm->GetTypeName()).BeginsWith("UInt_t")) typestring = "i";
1012     else if (TString(tdm->GetTypeName()).BeginsWith("Float_t")) typestring = "F";
1013     else if (TString(tdm->GetTypeName()).BeginsWith("Double_t")) typestring = "D";
1014     else if (TString(tdm->GetTypeName()).BeginsWith("Long64_t")) typestring = "L";
1015     else if (TString(tdm->GetTypeName()).BeginsWith("ULong64_t")) typestring = "l";
1016     else if (TString(tdm->GetTypeName()).BeginsWith("Bool_t")) typestring = "O";
1017     else continue;
1018 bendavid 1.1 //printf("%s\n",tdm->GetTypeName());
1019     TString varname = TString::Format("ph%d.%s",iph+1,tdm->GetName());
1020 bendavid 1.20 if (tdm->GetArrayDim()==1) {
1021     varname = TString::Format("%s[%i]",varname.Data(),tdm->GetMaxIndex(0));
1022     }
1023    
1024     //printf("typename = %s, arraydim = %i, arraysize = %i,varname = %s\n", tdm->GetTypeName(), tdm->GetArrayDim(), tdm->GetMaxIndex(0), varname.Data());
1025 bendavid 1.1
1026     Char_t *addr = (Char_t*)&fDiphotonEvent->photons[iph];
1027     assert(sizeof(Char_t)==1);
1028 bendavid 1.20 if (fWriteDiphotonTree) hCiCTuple->Branch(varname,addr+tdm->GetOffset(),TString::Format("%s/%s",varname.Data(),typestring.Data()));
1029 bendavid 1.1
1030     if (iph==0) {
1031     TString singlename = TString::Format("ph.%s",tdm->GetName());
1032 bendavid 1.20 if (tdm->GetArrayDim()==1) {
1033     singlename = TString::Format("%s[%i]",singlename.Data(),tdm->GetMaxIndex(0));
1034     }
1035 bendavid 1.1 Char_t *addrsingle = (Char_t*)fSinglePhoton;
1036 bendavid 1.20 if (fWriteSingleTree) hCiCTupleSingle->Branch(singlename,addrsingle+tdm->GetOffset(),TString::Format("%s/%s",singlename.Data(),typestring.Data()));
1037 bendavid 1.1 }
1038     }
1039     }
1040 paus 1.12
1041     if (fWriteDiphotonTree)
1042     AddOutput(hCiCTuple);
1043     if (fWriteSingleTree)
1044     AddOutput(hCiCTupleSingle);
1045 bendavid 1.22
1046     if (fFillVertexTree) {
1047     fDiphotonVtx = new PhotonTreeWriterVtx;
1048     hVtxTree = new TTree("hVtxTree","hVtxTree");
1049     hVtxTree->SetAutoSave(300e9);
1050     AddOutput(hVtxTree);
1051    
1052     TClass *vclass = TClass::GetClass("mithep::PhotonTreeWriterVtx");
1053     TList *vlist = vclass->GetListOfDataMembers();
1054    
1055     for (int i=0; i<vlist->GetEntries(); ++i) {
1056     const TDataMember *tdm = static_cast<const TDataMember*>(vlist->At(i));
1057     if (!(tdm->IsBasic() && tdm->IsPersistent())) continue;
1058     TString typestring;
1059     if (TString(tdm->GetTypeName()).BeginsWith("Char_t")) typestring = "B";
1060     else if (TString(tdm->GetTypeName()).BeginsWith("UChar_t")) typestring = "b";
1061     else if (TString(tdm->GetTypeName()).BeginsWith("Short_t")) typestring = "S";
1062     else if (TString(tdm->GetTypeName()).BeginsWith("UShort_t")) typestring = "s";
1063     else if (TString(tdm->GetTypeName()).BeginsWith("Int_t")) typestring = "I";
1064     else if (TString(tdm->GetTypeName()).BeginsWith("UInt_t")) typestring = "i";
1065     else if (TString(tdm->GetTypeName()).BeginsWith("Float_t")) typestring = "F";
1066     else if (TString(tdm->GetTypeName()).BeginsWith("Double_t")) typestring = "D";
1067     else if (TString(tdm->GetTypeName()).BeginsWith("Long64_t")) typestring = "L";
1068     else if (TString(tdm->GetTypeName()).BeginsWith("ULong64_t")) typestring = "l";
1069     else if (TString(tdm->GetTypeName()).BeginsWith("Bool_t")) typestring = "O";
1070     else continue;
1071     //printf("%s %s: %i\n",tdm->GetTypeName(),tdm->GetName(),int(tdm->GetOffset()));
1072     Char_t *addr = (Char_t*)fDiphotonVtx;
1073     assert(sizeof(Char_t)==1);
1074     hVtxTree->Branch(tdm->GetName(),addr + tdm->GetOffset(),TString::Format("%s/%s",tdm->GetName(),typestring.Data()));
1075     }
1076     }
1077    
1078 bendavid 1.1 }
1079    
1080     // ----------------------------------------------------------------------------------------
1081     // some helpfer functions....
1082 paus 1.12 void PhotonTreeWriter::FindHiggsPtAndZ(Float_t& pt, Float_t& decayZ, Float_t& mass)
1083     {
1084     pt = -999.;
1085 bendavid 1.1 decayZ = -999.;
1086 paus 1.12 mass = -999.;
1087 bendavid 1.1
1088     // loop over all GEN particles and look for status 1 photons
1089     for(UInt_t i=0; i<fMCParticles->GetEntries(); ++i) {
1090     const MCParticle* p = fMCParticles->At(i);
1091 paus 1.12 if (p->Is(MCParticle::kH) || (!fApplyElectronVeto &&
1092     (p->AbsPdgId()==23 || p->AbsPdgId()==24))) {
1093     pt = p->Pt();
1094 bendavid 1.1 decayZ = p->DecayVertex().Z();
1095 paus 1.12 mass = p->Mass();
1096 bendavid 1.1 break;
1097     }
1098     }
1099    
1100     return;
1101     }
1102    
1103    
1104 paus 1.12 Float_t PhotonTreeWriter::GetEventCat(PhotonTools::CiCBaseLineCats cat1,
1105     PhotonTools::CiCBaseLineCats cat2) {
1106 bendavid 1.1
1107 paus 1.12 bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
1108     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
1109 bendavid 1.1
1110     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
1111     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
1112    
1113     if( ph1IsEB && ph2IsEB )
1114     return ( ph1IsHR9 && ph2IsHR9 ? 0. : 1.);
1115    
1116     return ( ph1IsHR9 && ph2IsHR9 ? 2. : 3.);
1117     }
1118    
1119 bendavid 1.20 template <int NClus>
1120     void PhotonTreeWriterPhoton<NClus>::SetVars(const Photon *p, const DecayParticle *c, const Electron *ele,
1121 fabstoec 1.24 const SuperCluster *pfsc, const MCParticle *m,
1122     PhotonFix &phfixph, PhotonFix &phfixele,
1123     const TrackCol* trackCol,const VertexCol* vtxCol,
1124     const PFCandidateCol* fPFCands,
1125     Double_t rho,
1126     Bool_t fillclusterarrays,
1127     const ElectronCol* els, Bool_t applyElectronVeto) {
1128 paus 1.12
1129     const SuperCluster *s = 0;
1130     if (p)
1131     s = p->SCluster();
1132     else
1133     s = ele->SCluster();
1134     const BasicCluster *b = s->Seed();
1135     const BasicCluster *b2 = 0;
1136     Double_t ebcmax = -99.;
1137     for (UInt_t i=0; i<s->ClusterSize(); ++i) {
1138     const BasicCluster *bc = s->Cluster(i);
1139     if (bc->Energy() > ebcmax && bc !=b) {
1140     b2 = bc;
1141     ebcmax = bc->Energy();
1142     }
1143     }
1144 bendavid 1.2
1145 paus 1.12 const BasicCluster *bclast = 0;
1146     Double_t ebcmin = 1e6;
1147     for (UInt_t i=0; i<s->ClusterSize(); ++i) {
1148     const BasicCluster *bc = s->Cluster(i);
1149     if (bc->Energy() < ebcmin && bc !=b) {
1150     bclast = bc;
1151     ebcmin = bc->Energy();
1152     }
1153     }
1154 bendavid 1.5
1155 paus 1.12 const BasicCluster *bclast2 = 0;
1156     ebcmin = 1e6;
1157     for (UInt_t i=0; i<s->ClusterSize(); ++i) {
1158     const BasicCluster *bc = s->Cluster(i);
1159     if (bc->Energy() < ebcmin && bc !=b && bc!=bclast) {
1160     bclast2 = bc;
1161     ebcmin = bc->Energy();
1162     }
1163     }
1164 bendavid 1.5
1165 paus 1.12 if (p) {
1166     hasphoton = kTRUE;
1167     e = p->E();
1168     pt = p->Pt();
1169     eta = p->Eta();
1170     phi = p->Phi();
1171     r9 = p->R9();
1172     e3x3 = p->E33();
1173     e5x5 = p->E55();
1174     hovere = p->HadOverEm();
1175 bendavid 1.20 hoveretower = p->HadOverEmTow();
1176 paus 1.12 sigietaieta = p->CoviEtaiEta();
1177     phcat = PhotonTools::CiCBaseLineCat(p);
1178     eerr = p->EnergyErr();
1179     eerrsmeared = p->EnergyErrSmeared();
1180     esmearing = p->EnergySmearing();
1181     idmva = p->IdMva();
1182     hcalisodr03 = p->HcalTowerSumEtDr03();
1183     ecalisodr03 = p->EcalRecHitIsoDr03();
1184     trkisohollowdr03 = p->HollowConeTrkIsoDr03();
1185     hcalisodr04 = p->HcalTowerSumEtDr04();
1186     ecalisodr04 = p->EcalRecHitIsoDr04();
1187     trkisohollowdr04 = p->HollowConeTrkIsoDr04();
1188    
1189    
1190     const Vertex *vtx = vtxCol->At(0);
1191     if (p->HasPV()) vtx = p->PV();
1192    
1193     UInt_t wVtxInd = 0;
1194    
1195     trackiso1 = IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0,
1196     trackCol, NULL, NULL,
1197     (!applyElectronVeto ? els : NULL) );
1198     //Question Ming:whyfPV->At(0) instead of selected vertex using ranking method?
1199 bendavid 1.2
1200 paus 1.12 // track iso worst vtx
1201     trackiso2 = IsolationTools::CiCTrackIsolation(p,vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0,
1202     trackCol, &wVtxInd,vtxCol,
1203     (!applyElectronVeto ? els : NULL) );
1204     combiso1 = ecalisodr03+hcalisodr04+trackiso1 - 0.17*rho;
1205 fabstoec 1.24 combiso2 = ecalisodr04+hcalisodr04+trackiso2 - 0.52*rho;
1206    
1207    
1208     // -----------------------------------------------------
1209     // PF-CiC4 Debug Stuff
1210     std::vector<double> debugVals;
1211     bool tmpPass = PhotonTools::PassCiCPFIsoSelection(p, vtx, fPFCands, vtxCol, rho, 20., &debugVals);
1212     if( debugVals.size() == 13 ) {
1213     pfcic4_tIso1 = debugVals[0];
1214     pfcic4_tIso2 = debugVals[1];
1215     pfcic4_tIso3 = debugVals[2];
1216     pfcic4_covIEtaIEta = debugVals[3];
1217     pfcic4_HoE = debugVals[4];
1218     pfcic4_R9 = debugVals[5];
1219     pfcic4_wVtxInd = debugVals[6];
1220     pfcic4_ecalIso3 = debugVals[7];
1221     pfcic4_ecalIso4 = debugVals[8];
1222     pfcic4_trackIsoSel03 = debugVals[9];
1223     pfcic4_trackIsoWorst04 = debugVals[10];
1224     pfcic4_combIso1 = debugVals[11];
1225     pfcic4_combIso2 = debugVals[12];
1226     }
1227     // -----------------------------------------------------
1228 mingyang 1.26 //id mva
1229     //2011
1230     idmva_tIso1abs=combiso1;
1231     idmva_tIso2abs=combiso2;
1232     idmva_tIso3abs=trackiso1;
1233     idmva_absIsoEcal=ecalisodr03;
1234     idmva_absIsoHcal=hcalisodr04;
1235     //2012
1236     idmva_CoviEtaiPhi=p->SCluster()->Seed()->CoviEtaiPhi();
1237     idmva_s4ratio=p->SCluster()->Seed()->E2x2()/p->SCluster()->Seed()->E5x5();
1238     idmva_GammaIso=IsolationTools::PFGammaIsolation(p,0.3,0,fPFCands);
1239     idmva_ChargedIso_selvtx=IsolationTools::PFChargedIsolation(p,vtx,0.3,0.,fPFCands);
1240     idmva_ChargedIso_worstvtx=IsolationTools::PFChargedIsolation(p,vtx,0.3,0.,fPFCands,&wVtxInd,vtxCol);
1241     idmva_PsEffWidthSigmaRR=sqrt(p->SCluster()->PsEffWidthSigmaXX()*p->SCluster()->PsEffWidthSigmaXX()+p->SCluster()->PsEffWidthSigmaYY()*p->SCluster()->PsEffWidthSigmaYY());
1242 paus 1.12 }
1243     else {
1244     hasphoton = kFALSE;
1245     e = -99.;
1246     pt = -99.;
1247     eta = -99.;
1248     phi = -99.;
1249     r9 = b->E3x3()/s->RawEnergy();
1250     e3x3 = b->E3x3();
1251     e5x5 = b->E5x5();
1252     hovere = ele->HadronicOverEm();
1253     sigietaieta = ele->CoviEtaiEta();
1254     phcat = -99;
1255     eerr = -99.;
1256     eerrsmeared = -99.;
1257     esmearing = 0.;
1258     idmva = -99.;
1259     hcalisodr03 = -99;
1260     ecalisodr03 = -99;
1261     trkisohollowdr03 = -99;
1262     hcalisodr04 = -99;
1263     ecalisodr04 = -99;
1264     trkisohollowdr04 = -99;
1265     trackiso1 = -99.;
1266     trackiso2 = -99.;
1267     combiso1 = -99.;
1268     combiso2 = -99.;
1269     }
1270    
1271     sce = s->Energy();
1272     scrawe = s->RawEnergy();
1273     scpse = s->PreshowerEnergy();
1274 bendavid 1.20 scpssigmaxx = s->PsEffWidthSigmaXX();
1275     scpssigmayy = s->PsEffWidthSigmaYY();
1276 paus 1.12 sceta = s->Eta();
1277     scphi = s->Phi();
1278     scnclusters = s->ClusterSize();
1279     scnhits = s->NHits();
1280     scetawidth = -99.;
1281     scphiwidth = -99.;
1282     if (p) {
1283     scetawidth = p->EtaWidth();
1284     scphiwidth = p->PhiWidth();
1285     }
1286     else {
1287     scetawidth = s->EtaWidth();
1288     scphiwidth = s->PhiWidth();
1289     }
1290     isbarrel = (s->AbsEta()<1.5);
1291     isr9reco = (isbarrel && r9>0.94) || (!isbarrel && r9>0.95);
1292     isr9cat = (r9>0.94);
1293    
1294     eseed = b->Energy();
1295     etaseed = b->Eta();
1296     phiseed = b->Phi();
1297     ietaseed = b->IEta();
1298     iphiseed = b->IPhi();
1299     ixseed = b->IX();
1300     iyseed = b->IY();
1301     etacryseed = b->EtaCry();
1302     phicryseed = b->PhiCry();
1303     xcryseed = b->XCry();
1304     ycryseed = b->YCry();
1305     thetaaxisseed = b->ThetaAxis();
1306     phiaxisseed = b->PhiAxis();
1307     sigietaietaseed = TMath::Sqrt(b->CoviEtaiEta());
1308     sigiphiphiseed = TMath::Sqrt(b->CoviPhiiPhi());
1309     if (isnan(sigiphiphiseed)) sigiphiphiseed = -99.;
1310     covietaiphiseed = b->CoviEtaiPhi();
1311     if (isnan(covietaiphiseed)) covietaiphiseed = -99.;
1312     e3x3seed = b->E3x3();
1313     e5x5seed = b->E5x5();
1314     emaxseed = b->EMax();
1315     e2ndseed = b->E2nd();
1316     etopseed = b->ETop();
1317     ebottomseed = b->EBottom();
1318     eleftseed = b->ELeft();
1319     erightseed = b->ERight();
1320     e1x3seed = b->E1x3();
1321     e3x1seed = b->E3x1();
1322     e1x5seed = b->E1x5();
1323     e2x2seed = b->E2x2();
1324     e4x4seed = b->E4x4();
1325     e2x5maxseed = b->E2x5Max();
1326     e2x5topseed = b->E2x5Top();
1327     e2x5bottomseed = b->E2x5Bottom();
1328     e2x5leftseed = b->E2x5Left();
1329     e2x5rightseed = b->E2x5Right();
1330     xseedseed = b->Pos().X();
1331     yseedseed = b->Pos().Y();
1332     zseedseed = b->Pos().Z();
1333     nhitsseed = b->NHits();
1334    
1335     if (b2) {
1336     ebc2 = b2->Energy();
1337     etabc2 = b2->Eta();
1338     phibc2 = b2->Phi();
1339     ietabc2 = b2->IEta();
1340     iphibc2 = b2->IPhi();
1341     ixbc2 = b2->IX();
1342     iybc2 = b2->IY();
1343     etacrybc2 = b2->EtaCry();
1344     phicrybc2 = b2->PhiCry();
1345     xcrybc2 = b2->XCry();
1346     ycrybc2 = b2->YCry();
1347     thetaaxisbc2 = b2->ThetaAxis();
1348     phiaxisbc2 = b2->PhiAxis();
1349     sigietaietabc2 = TMath::Sqrt(b2->CoviEtaiEta());
1350     sigiphiphibc2 = TMath::Sqrt(b2->CoviPhiiPhi());
1351     if (isnan(sigiphiphibc2)) sigiphiphibc2 = -99.;
1352     covietaiphibc2 = b2->CoviEtaiPhi();
1353     if (isnan(covietaiphibc2)) covietaiphibc2 = -99.;
1354     e3x3bc2 = b2->E3x3();
1355     e5x5bc2 = b2->E5x5();
1356     emaxbc2 = b2->EMax();
1357     e2ndbc2 = b2->E2nd();
1358     etopbc2 = b2->ETop();
1359     ebottombc2 = b2->EBottom();
1360     eleftbc2 = b2->ELeft();
1361     erightbc2 = b2->ERight();
1362     e1x3bc2 = b2->E1x3();
1363     e3x1bc2 = b2->E3x1();
1364     e1x5bc2 = b2->E1x5();
1365     e2x2bc2 = b2->E2x2();
1366     e4x4bc2 = b2->E4x4();
1367     e2x5maxbc2 = b2->E2x5Max();
1368     e2x5topbc2 = b2->E2x5Top();
1369     e2x5bottombc2 = b2->E2x5Bottom();
1370     e2x5leftbc2 = b2->E2x5Left();
1371     e2x5rightbc2 = b2->E2x5Right();
1372     xbc2bc2 = b2->Pos().X();
1373     ybc2bc2 = b2->Pos().Y();
1374     zbc2bc2 = b2->Pos().Z();
1375     nhitsbc2= b2->NHits();
1376     }
1377     else {
1378     ebc2 = 0;
1379     etabc2 = 0;
1380     phibc2 = 0;
1381     ietabc2 = 0;
1382     iphibc2 = 0;
1383     ixbc2 = 0;
1384     iybc2 = 0;
1385     etacrybc2 = 0;
1386     phicrybc2 = 0;
1387     xcrybc2 = 0;
1388     ycrybc2 = 0;
1389     thetaaxisbc2 = 0;
1390     phiaxisbc2 = 0;
1391     sigietaietabc2 = 0;
1392     sigiphiphibc2 = 0;
1393     covietaiphibc2 = 0;
1394     e3x3bc2 = 0;
1395     e5x5bc2 = 0;
1396     emaxbc2 = 0;
1397     e2ndbc2 = 0;
1398     etopbc2 = 0;
1399     ebottombc2 = 0;
1400     eleftbc2 = 0;
1401     erightbc2 = 0;
1402     e1x3bc2 = 0;
1403     e3x1bc2 = 0;
1404     e1x5bc2 = 0;
1405     e2x2bc2 = 0;
1406     e4x4bc2 = 0;
1407     e2x5maxbc2 = 0;
1408     e2x5topbc2 = 0;
1409     e2x5bottombc2 = 0;
1410     e2x5leftbc2 = 0;
1411     e2x5rightbc2 = 0;
1412     xbc2bc2 = 0;
1413     ybc2bc2 = 0;
1414     zbc2bc2 = 0;
1415     nhitsbc2 = 0;
1416     }
1417 bendavid 1.5
1418 paus 1.12 if (bclast) {
1419     ebclast = bclast->Energy();
1420     etabclast = bclast->Eta();
1421     phibclast = bclast->Phi();
1422     ietabclast = bclast->IEta();
1423     iphibclast = bclast->IPhi();
1424     ixbclast = bclast->IX();
1425     iybclast = bclast->IY();
1426     etacrybclast = bclast->EtaCry();
1427     phicrybclast = bclast->PhiCry();
1428     xcrybclast = bclast->XCry();
1429     ycrybclast = bclast->YCry();
1430     thetaaxisbclast = bclast->ThetaAxis();
1431     phiaxisbclast = bclast->PhiAxis();
1432     sigietaietabclast = TMath::Sqrt(bclast->CoviEtaiEta());
1433     sigiphiphibclast = TMath::Sqrt(bclast->CoviPhiiPhi());
1434     if (isnan(sigiphiphibclast)) sigiphiphibclast = -99.;
1435     covietaiphibclast = bclast->CoviEtaiPhi();
1436     if (isnan(covietaiphibclast)) covietaiphibclast = -99.;
1437     e3x3bclast = bclast->E3x3();
1438     nhitsbclast = bclast->NHits();
1439     }
1440     else {
1441     ebclast = 0;
1442     etabclast = 0;
1443     phibclast = 0;
1444     ietabclast = 0;
1445     iphibclast = 0;
1446     ixbclast = 0;
1447     iybclast = 0;
1448     etacrybclast = 0;
1449     phicrybclast = 0;
1450     xcrybclast = 0;
1451     ycrybclast = 0;
1452     thetaaxisbclast = 0;
1453     phiaxisbclast = 0;
1454     sigietaietabclast = 0;
1455     sigiphiphibclast = 0;
1456     covietaiphibclast = 0;
1457     e3x3bclast = 0;
1458     nhitsbclast = 0;
1459     }
1460 bendavid 1.5
1461 paus 1.12 if (bclast2) {
1462     ebclast2 = bclast2->Energy();
1463     etabclast2 = bclast2->Eta();
1464     phibclast2 = bclast2->Phi();
1465     ietabclast2 = bclast2->IEta();
1466     iphibclast2 = bclast2->IPhi();
1467     ixbclast2 = bclast2->IX();
1468     iybclast2 = bclast2->IY();
1469     etacrybclast2 = bclast2->EtaCry();
1470     phicrybclast2 = bclast2->PhiCry();
1471     xcrybclast2 = bclast2->XCry();
1472     ycrybclast2 = bclast2->YCry();
1473     thetaaxisbclast2 = bclast2->ThetaAxis();
1474     phiaxisbclast2 = bclast2->PhiAxis();
1475     sigietaietabclast2 = TMath::Sqrt(bclast2->CoviEtaiEta());
1476     sigiphiphibclast2 = TMath::Sqrt(bclast2->CoviPhiiPhi());
1477     if (isnan(sigiphiphibclast2)) sigiphiphibclast2 = -99.;
1478     covietaiphibclast2 = bclast2->CoviEtaiPhi();
1479     if (isnan(covietaiphibclast2)) covietaiphibclast2 = -99.;
1480     e3x3bclast2 = bclast2->E3x3();
1481     nhitsbclast2 = bclast2->NHits();
1482     }
1483     else {
1484     ebclast2 = 0;
1485     etabclast2 = 0;
1486     phibclast2 = 0;
1487     ietabclast2 = 0;
1488     iphibclast2 = 0;
1489     ixbclast2 = 0;
1490     iybclast2 = 0;
1491     etacrybclast2 = 0;
1492     phicrybclast2 = 0;
1493     xcrybclast2 = 0;
1494     ycrybclast2 = 0;
1495     thetaaxisbclast2 = 0;
1496     phiaxisbclast2 = 0;
1497     sigietaietabclast2 = 0;
1498     sigiphiphibclast2 = 0;
1499     covietaiphibclast2 = 0;
1500     e3x3bclast2 = 0;
1501     nhitsbclast2 = 0;
1502     }
1503 bendavid 1.5
1504 bendavid 1.20 for (UInt_t iclus=0; iclus<NClus; ++iclus) {
1505     if (fillclusterarrays && iclus < s->ClusterSize() ) {
1506     const BasicCluster *ib =s->Cluster(iclus);
1507    
1508     ebcs[iclus] = ib->Energy();
1509     etabcs[iclus] = ib->Eta();
1510     phibcs[iclus] = ib->Phi();
1511     ietabcs[iclus] = ib->IEta();
1512     iphibcs[iclus] = ib->IPhi();
1513     ixbcs[iclus] = ib->IX();
1514     iybcs[iclus] = ib->IY();
1515     etacrybcs[iclus] = ib->EtaCry();
1516     phicrybcs[iclus] = ib->PhiCry();
1517     xcrybcs[iclus] = ib->XCry();
1518     ycrybcs[iclus] = ib->YCry();
1519     sigietaietabcs[iclus] = TMath::Sqrt(ib->CoviEtaiEta());
1520     sigiphiphibcs[iclus] = TMath::Sqrt(ib->CoviPhiiPhi());
1521     covietaiphibcs[iclus] = ib->CoviEtaiPhi();
1522     sigetaetabcs[iclus] = TMath::Sqrt(ib->CovEtaEta());
1523     sigphiphibcs[iclus] = TMath::Sqrt(ib->CovPhiPhi());
1524     covetaphibcs[iclus] = ib->CovEtaPhi();
1525     e3x3bcs[iclus] = ib->E3x3();
1526     e5x5bcs[iclus] = ib->E5x5();
1527     emaxbcs[iclus] = ib->EMax();
1528     e2ndbcs[iclus] = ib->E2nd();
1529     etopbcs[iclus] = ib->ETop();
1530     ebottombcs[iclus] = ib->EBottom();
1531     eleftbcs[iclus] = ib->ELeft();
1532     erightbcs[iclus] = ib->ERight();
1533     e1x3bcs[iclus] = ib->E1x3();
1534     e3x1bcs[iclus] = ib->E3x1();
1535     e1x5bcs[iclus] = ib->E1x5();
1536     e2x2bcs[iclus] = ib->E2x2();
1537     e4x4bcs[iclus] = ib->E4x4();
1538     e2x5maxbcs[iclus] = ib->E2x5Max();
1539     e2x5topbcs[iclus] = ib->E2x5Top();
1540     e2x5bottombcs[iclus] = ib->E2x5Bottom();
1541     e2x5leftbcs[iclus] = ib->E2x5Left();
1542     e2x5rightbcs[iclus] = ib->E2x5Right();
1543     nhitsbcs[iclus]= ib->NHits();
1544     }
1545     else {
1546     ebcs[iclus] = -999;
1547     etabcs[iclus] = -999;
1548     phibcs[iclus] = -999;
1549     ietabcs[iclus] = -999;
1550     iphibcs[iclus] = -999;
1551     ixbcs[iclus] = -999;
1552     iybcs[iclus] = -999;
1553     etacrybcs[iclus] = -999;
1554     phicrybcs[iclus] = -999;
1555     xcrybcs[iclus] = -999;
1556     ycrybcs[iclus] = -999;
1557     sigietaietabcs[iclus] = -999;
1558     sigiphiphibcs[iclus] = -999;
1559     covietaiphibcs[iclus] = -999;
1560     sigetaetabcs[iclus] = -999;
1561     sigphiphibcs[iclus] = -999;
1562     covetaphibcs[iclus] = -999;
1563     e3x3bcs[iclus] = -999;
1564     e5x5bcs[iclus] = -999;
1565     emaxbcs[iclus] = -999;
1566     e2ndbcs[iclus] = -999;
1567     etopbcs[iclus] = -999;
1568     ebottombcs[iclus] = -999;
1569     eleftbcs[iclus] = -999;
1570     erightbcs[iclus] = -999;
1571     e1x3bcs[iclus] = -999;
1572     e3x1bcs[iclus] = -999;
1573     e1x5bcs[iclus] = -999;
1574     e2x2bcs[iclus] = -999;
1575     e4x4bcs[iclus] = -999;
1576     e2x5maxbcs[iclus] = -999;
1577     e2x5topbcs[iclus] = -999;
1578     e2x5bottombcs[iclus] = -999;
1579     e2x5leftbcs[iclus] = -999;
1580     e2x5rightbcs[iclus] = -999;
1581     nhitsbcs[iclus] = -999;
1582     }
1583     }
1584    
1585     for (UInt_t iclus=0; iclus<NClus; ++iclus) {
1586     if (fillclusterarrays && pfsc && iclus < pfsc->ClusterSize() ) {
1587     const BasicCluster *ib =pfsc->Cluster(iclus);
1588    
1589     epfbcs[iclus] = ib->Energy();
1590     etapfbcs[iclus] = ib->Eta();
1591     phipfbcs[iclus] = ib->Phi();
1592     ietapfbcs[iclus] = ib->IEta();
1593     iphipfbcs[iclus] = ib->IPhi();
1594     ixpfbcs[iclus] = ib->IX();
1595     iypfbcs[iclus] = ib->IY();
1596     etacrypfbcs[iclus] = ib->EtaCry();
1597     phicrypfbcs[iclus] = ib->PhiCry();
1598     xcrypfbcs[iclus] = ib->XCry();
1599     ycrypfbcs[iclus] = ib->YCry();
1600     sigietaietapfbcs[iclus] = TMath::Sqrt(ib->CoviEtaiEta());
1601     sigiphiphipfbcs[iclus] = TMath::Sqrt(ib->CoviPhiiPhi());
1602     covietaiphipfbcs[iclus] = ib->CoviEtaiPhi();
1603     sigetaetapfbcs[iclus] = TMath::Sqrt(ib->CovEtaEta());
1604     sigphiphipfbcs[iclus] = TMath::Sqrt(ib->CovPhiPhi());
1605     covetaphipfbcs[iclus] = ib->CovEtaPhi();
1606     e3x3pfbcs[iclus] = ib->E3x3();
1607     e5x5pfbcs[iclus] = ib->E5x5();
1608     emaxpfbcs[iclus] = ib->EMax();
1609     e2ndpfbcs[iclus] = ib->E2nd();
1610     etoppfbcs[iclus] = ib->ETop();
1611     ebottompfbcs[iclus] = ib->EBottom();
1612     eleftpfbcs[iclus] = ib->ELeft();
1613     erightpfbcs[iclus] = ib->ERight();
1614     e1x3pfbcs[iclus] = ib->E1x3();
1615     e3x1pfbcs[iclus] = ib->E3x1();
1616     e1x5pfbcs[iclus] = ib->E1x5();
1617     e2x2pfbcs[iclus] = ib->E2x2();
1618     e4x4pfbcs[iclus] = ib->E4x4();
1619     e2x5maxpfbcs[iclus] = ib->E2x5Max();
1620     e2x5toppfbcs[iclus] = ib->E2x5Top();
1621     e2x5bottompfbcs[iclus] = ib->E2x5Bottom();
1622     e2x5leftpfbcs[iclus] = ib->E2x5Left();
1623     e2x5rightpfbcs[iclus] = ib->E2x5Right();
1624     nhitspfbcs[iclus]= ib->NHits();
1625     }
1626     else {
1627     epfbcs[iclus] = -999;
1628     etapfbcs[iclus] = -999;
1629     phipfbcs[iclus] = -999;
1630     ietapfbcs[iclus] = -999;
1631     iphipfbcs[iclus] = -999;
1632     ixpfbcs[iclus] = -999;
1633     iypfbcs[iclus] = -999;
1634     etacrypfbcs[iclus] = -999;
1635     phicrypfbcs[iclus] = -999;
1636     xcrypfbcs[iclus] = -999;
1637     ycrypfbcs[iclus] = -999;
1638     sigietaietapfbcs[iclus] = -999;
1639     sigiphiphipfbcs[iclus] = -999;
1640     covietaiphipfbcs[iclus] = -999;
1641     sigetaetapfbcs[iclus] = -999;
1642     sigphiphipfbcs[iclus] = -999;
1643     covetaphipfbcs[iclus] = -999;
1644     e3x3pfbcs[iclus] = -999;
1645     e5x5pfbcs[iclus] = -999;
1646     emaxpfbcs[iclus] = -999;
1647     e2ndpfbcs[iclus] = -999;
1648     etoppfbcs[iclus] = -999;
1649     ebottompfbcs[iclus] = -999;
1650     eleftpfbcs[iclus] = -999;
1651     erightpfbcs[iclus] = -999;
1652     e1x3pfbcs[iclus] = -999;
1653     e3x1pfbcs[iclus] = -999;
1654     e1x5pfbcs[iclus] = -999;
1655     e2x2pfbcs[iclus] = -999;
1656     e4x4pfbcs[iclus] = -999;
1657     e2x5maxpfbcs[iclus] = -999;
1658     e2x5toppfbcs[iclus] = -999;
1659     e2x5bottompfbcs[iclus] = -999;
1660     e2x5leftpfbcs[iclus] = -999;
1661     e2x5rightpfbcs[iclus] = -999;
1662     nhitspfbcs[iclus] = -999;
1663     }
1664     }
1665    
1666     for (UInt_t iclus=0; iclus<100; ++iclus) {
1667     if (fillclusterarrays && pfsc && iclus < pfsc->NPsClusts() ) {
1668     const PsCluster *ib = pfsc->PsClust(iclus);
1669    
1670     epsc[iclus] = ib->Energy();
1671     etapsc[iclus] = ib->Eta();
1672     phipsc[iclus] = ib->Phi();
1673     planepsc[iclus] = ib->PsPlane();
1674     }
1675     else {
1676     epsc[iclus] = -999;
1677     etapsc[iclus] = -999;
1678     phipsc[iclus] = -999;
1679     planepsc[iclus] = 0;
1680     }
1681     }
1682    
1683 paus 1.12 //initialize photon energy corrections if needed
1684     /*if (!PhotonFix::initialised()) {
1685     PhotonFix::initialise("4_2",std::string((gSystem->Getenv("CMSSW_BASE") + TString("/src/MitPhysics/data/PhotonFix.dat")).Data()));
1686     }*/
1687    
1688     phfixph.setup(e,sceta,scphi,r9);
1689     phfixele.setup(e,sceta,scphi,r9);
1690    
1691     const Float_t dval = -99.;
1692     ecor = phfixph.fixedEnergy();
1693     ecorerr = phfixph.sigmaEnergy();
1694     ecorele = phfixele.fixedEnergy();
1695     ecoreleerr = phfixele.sigmaEnergy();
1696     if (phfixph.isbarrel()) {
1697     etac = phfixph.etaC();
1698     etas = phfixph.etaS();
1699     etam = phfixph.etaM();
1700     phic = phfixph.phiC();
1701     phis = phfixph.phiS();
1702     phim = phfixph.phiM();
1703     xz = dval;
1704     xc = dval;
1705     xs = dval;
1706     xm = dval;
1707     yz = dval;
1708     yc = dval;
1709     ys = dval;
1710     ym = dval;
1711     }
1712     else {
1713     etac = dval;
1714     etas = dval;
1715     etam = dval;
1716     phic = dval;
1717     phis = dval;
1718     phim = dval;
1719     xz = phfixph.xZ();
1720     xc = phfixph.xC();
1721     xs = phfixph.xS();
1722     xm = phfixph.xM();
1723     yz = phfixph.yZ();
1724     yc = phfixph.yC();
1725     ys = phfixph.yS();
1726     ym = phfixph.yM();
1727     }
1728    
1729     if (c) {
1730     hasconversion = kTRUE;
1731     convp = c->P();
1732     convpt = c->Pt();
1733     conveta = c->Eta();
1734     convphi = c->Phi();
1735     ThreeVector dirconvsc = ThreeVector(s->Point()) - c->Position();
1736     convdeta = c->Eta() - dirconvsc.Eta();
1737     convdphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
1738     convvtxrho = c->Position().Rho();
1739     convvtxz = c->Position().Z();
1740     convvtxphi = c->Position().Phi();
1741    
1742     const StableData *leadsd = dynamic_cast<const StableData*>(c->DaughterDat(0));
1743     const StableData *trailsd = dynamic_cast<const StableData*>(c->DaughterDat(1));
1744     if (leadsd->Pt()<trailsd->Pt()) {
1745     const StableData *sdtmp = leadsd;
1746     leadsd = trailsd;
1747     trailsd = sdtmp;
1748     }
1749    
1750     const Track *leadtrack = dynamic_cast<const StableParticle*>(leadsd->Original())->Trk();
1751     const Track *trailtrack = dynamic_cast<const StableParticle*>(trailsd->Original())->Trk();
1752    
1753     convleadpt = leadsd->Pt();
1754     convtrailpt = trailsd->Pt();
1755     convleadtrackpt = leadtrack->Pt();
1756     convleadtrackalgo = leadtrack->Algo();
1757     if (convleadtrackalgo==29) convleadtrackalgos=2; //gsf track
1758     else if (convleadtrackalgo==15 ||convleadtrackalgo==16) convleadtrackalgos=1; //ecal-seeded track
1759     else convleadtrackalgos = 0; //std iterative track
1760     convleadtrackcharge = leadtrack->Charge();
1761     convtrailtrackpt = trailtrack->Pt();
1762     convtrailtrackalgo = trailtrack->Algo();
1763     if (convtrailtrackalgo==29) convtrailtrackalgos=2; //gsf track
1764     else if (convtrailtrackalgo==15 ||convtrailtrackalgo==16) convtrailtrackalgos=1; //ecal-seeded track
1765     else convtrailtrackalgos = 0; //std iterative track
1766     trailtrackcharge = trailtrack->Charge();
1767     }
1768     else {
1769     hasconversion = kFALSE;
1770     convp = -99.;
1771     convpt = -99.;
1772     conveta = -99.;
1773     convphi = -99.;
1774     convdeta = -99.;
1775     convdphi = -99.;
1776     convvtxrho = -99.;
1777     convvtxz = -999.;
1778     convvtxphi = -99.;
1779     convleadpt = -99.;
1780     convtrailpt = -99.;
1781     convleadtrackpt = -99.;
1782     convleadtrackalgo = -99;
1783     convleadtrackalgos = -99;
1784     convleadtrackcharge = 0;
1785     convtrailtrackpt = -99.;
1786     convtrailtrackalgo = -99;
1787     convtrailtrackalgos = -99;
1788     trailtrackcharge = 0;
1789     }
1790    
1791     //electron quantities
1792     if (ele) {
1793     haselectron = kTRUE;
1794     eleisecaldriven = ele->IsEcalDriven();
1795     eleistrackerdriven = ele->IsTrackerDriven();
1796     elee = ele->E();
1797     elept = ele->Pt();
1798     eleeta = ele->Eta();
1799     elephi = ele->Phi();
1800     elecharge = ele->Charge();
1801     elefbrem = ele->FBrem();
1802     eledeta = ele->DeltaEtaSuperClusterTrackAtVtx();
1803     eledphi = ele->DeltaPhiSuperClusterTrackAtVtx();
1804     elep = s->Energy()/ele->ESuperClusterOverP();
1805     elepin = ele->PIn();
1806     elepout = ele->POut();
1807     }
1808     else {
1809     haselectron = kFALSE;
1810     eleisecaldriven = kFALSE;
1811     eleistrackerdriven = kFALSE;
1812     elee = -99.;
1813     elept = -99.;
1814     eleeta = -99.;
1815     elephi = -99.;
1816     elecharge = -99;
1817     elefbrem = -99.;
1818     eledeta = -99.;
1819     eledphi = -99.;
1820     elep = -99.;
1821     elepin = -99.;
1822     elepout = -99.;
1823     }
1824    
1825     //pf supercluster quantities
1826     if (pfsc) {
1827     haspfsc = kTRUE;
1828     pfsce = pfsc->Energy();
1829     pfscrawe = pfsc->RawEnergy();
1830     pfsceta = pfsc->Eta();
1831     pfscphi = pfsc->Phi();
1832 bendavid 1.20 pfscnclusters = pfsc->NClusters();
1833     pfscnhits = pfsc->NHits();
1834     pfscetawidth = pfsc->EtaWidth();
1835     pfscphiwidth = pfsc->PhiWidth();
1836     pfscnpsclusters = pfsc->NPsClusts();
1837 paus 1.12 }
1838     else {
1839     haspfsc = kFALSE;
1840     pfsce = -99.;
1841     pfscrawe = -99.;
1842     pfsceta = -99.;
1843     pfscphi = -99.;
1844 bendavid 1.20 pfscnclusters = 0;
1845     pfscnhits = 0;
1846     pfscetawidth = -99.;
1847     pfscphiwidth = -99.;
1848     pfscnpsclusters = 0;
1849 paus 1.12 }
1850    
1851     genz = -99.;
1852     if (m) {
1853     ispromptgen = kTRUE;
1854     gene = m->E();
1855     genpt = m->Pt();
1856     geneta = m->Eta();
1857     genphi = m->Phi();
1858     const MCParticle *mm = m->DistinctMother();
1859     if (mm) genz = mm->DecayVertex().Z();
1860     pdgid = m->PdgId();
1861     if (mm) motherpdgid = mm->PdgId();
1862     else motherpdgid = -99;
1863     }
1864     else {
1865     ispromptgen = kFALSE;
1866     gene = -99.;
1867     genpt = -99.;
1868     geneta = -99.;
1869     genphi = -99.;
1870     pdgid = -99;
1871     motherpdgid = -99;
1872     }
1873 bendavid 1.1 }
1874 bendavid 1.22
1875     void PhotonTreeWriterVtx::SetVars(const Vertex *v, const Photon *p1, const Photon *p2, const PFCandidateCol *pfcands, Int_t idx, Int_t numvtx, Float_t genvtxz)
1876     {
1877    
1878     //printf("start\n");
1879    
1880     n = idx;
1881     nvtx = numvtx;
1882     zgen = genvtxz;
1883    
1884     x = v->X();
1885     y = v->Y();
1886     z = v->Z();
1887    
1888     Double_t dsumpt = 0.;
1889     Double_t dsumptsq = 0.;
1890    
1891     nchalltoward = 0;
1892     nchalltransverse = 0;
1893     nchallaway = 0;
1894     nchcuttoward = 0;
1895     nchcuttransverse = 0;
1896     nchcutaway = 0;
1897    
1898     ThreeVector vtxmom;
1899    
1900     //printf("mom\n");
1901     FourVectorM diphoton = p1->MomVtx(v->Position()) + p2->MomVtx(v->Position());
1902     //printf("done mom\n");
1903     ptgg = diphoton.Pt();
1904     phigg = diphoton.Phi();
1905     etagg = diphoton.Eta();
1906     mgg = diphoton.M();
1907     pxgg = diphoton.Px();
1908     pygg = diphoton.Py();
1909     pzgg = diphoton.Pz();
1910    
1911     //printf("loop\n");
1912    
1913     for (UInt_t i = 0; i<pfcands->GetEntries(); ++i) {
1914     const PFCandidate *pfc = pfcands->At(i);
1915     if (pfc->PFType()!=PFCandidate::eHadron || !pfc->HasTrackerTrk()) continue;
1916     if (TMath::Abs( pfc->TrackerTrk()->DzCorrected(*v) ) > 0.2) continue;
1917     if (TMath::Abs( pfc->TrackerTrk()->D0Corrected(*v) ) > 0.1) continue;
1918    
1919     vtxmom += ThreeVector(pfc->Px(),pfc->Py(),pfc->Pz());
1920    
1921     dsumpt += pfc->Pt();
1922     dsumptsq += pfc->Pt()*pfc->Pt();
1923    
1924     Double_t dphi = TMath::Abs(MathUtils::DeltaPhi(*pfc,diphoton));
1925     if (dphi<(TMath::Pi()/3.0)) {
1926     ++nchalltoward;
1927     if (pfc->Pt()>0.5) ++nchcuttoward;
1928     }
1929     else if (dphi>(2.0*TMath::Pi()/3.0)) {
1930     ++nchallaway;
1931     if (pfc->Pt()>0.5) ++nchcutaway;
1932     }
1933     else {
1934     ++nchalltransverse;
1935     if (pfc->Pt()>0.5) ++nchcuttransverse;
1936     }
1937    
1938     }
1939    
1940     //printf("doneloop\n");
1941    
1942    
1943     sumpt = dsumpt;
1944     sumptsq = dsumptsq;
1945    
1946     pt = vtxmom.Rho();
1947     phi = vtxmom.Phi();
1948     eta = vtxmom.Eta();
1949     px = vtxmom.X();
1950     py = vtxmom.Y();
1951     pz = vtxmom.Z();
1952    
1953     //printf("done\n");
1954    
1955     return;
1956    
1957     }