ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/PhotonTreeWriter.cc
Revision: 1.21
Committed: Sun May 27 17:08:11 2012 UTC (12 years, 11 months ago) by bendavid
Content type: text/plain
Branch: MAIN
Changes since 1.20: +2 -1 lines
Log Message:
disable mva met temporarily

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