ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/PhotonTools.cc
Revision: 1.41
Committed: Mon Feb 11 13:26:28 2013 UTC (12 years, 2 months ago) by mingyang
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_029c, Mit_029b, Mit_029a
Changes since 1.40: +44 -1 lines
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 mingyang 1.41 // $Id: PhotonTools.cc,v 1.40 2013/01/08 13:24:01 fabstoec Exp $
2 bendavid 1.1
3     #include "MitPhysics/Utils/interface/PhotonTools.h"
4     #include "MitPhysics/Utils/interface/ElectronTools.h"
5 fabstoec 1.5 #include "MitPhysics/Utils/interface/IsolationTools.h"
6 bendavid 1.1 #include "MitAna/DataTree/interface/StableData.h"
7     #include <TFile.h>
8 fabstoec 1.8 #include <TRandom3.h>
9 bendavid 1.1
10     ClassImp(mithep::PhotonTools)
11    
12     using namespace mithep;
13    
14     //--------------------------------------------------------------------------------------------------
15     PhotonTools::PhotonTools()
16     {
17     // Constructor.
18     }
19    
20 bendavid 1.16 //--------------------------------------------------------------------------------------------------
21     PhotonTools::eScaleCats PhotonTools::EScaleCat(const Photon *p)
22     {
23     if (p->SCluster()->AbsEta()<1.0) {
24 bendavid 1.19 if (p->R9()>0.94) {
25     Int_t ieta = p->SCluster()->Seed()->IEta();
26     Int_t iphi = p->SCluster()->Seed()->IPhi();
27     Bool_t central = ( ((std::abs(ieta)-1)/5+1 >= 2 && (std::abs(ieta)-1)/5+1 < 5) || ((std::abs(ieta)-1)/5+1 >= 7 && (std::abs(ieta)-1)/5+1 < 9 ) || ((std::abs(ieta)-1)/5+1 >= 11 && (std::abs(ieta)-1)/5+1 < 13) || ((std::abs(ieta)-1)/5+1 >= 15 && (std::abs(ieta)-1)/5+1 < 17) ) && (iphi %20) > 5 && (iphi%20) <16;
28    
29     if (central) return kEBlowEtaGoldCenter;
30     else return kEBlowEtaGoldGap;
31     }
32 bendavid 1.16 else return kEBlowEtaBad;
33     }
34     else if (p->SCluster()->AbsEta()<1.5) {
35     if (p->R9()>0.94) return kEBhighEtaGold;
36     else return kEBhighEtaBad;
37     }
38     else if (p->SCluster()->AbsEta()<2.0) {
39     if (p->R9()>0.94) return kEElowEtaGold;
40     else return kEElowEtaBad;
41     }
42     else {
43     if (p->R9()>0.94) return kEEhighEtaGold;
44     else return kEEhighEtaBad;
45     }
46    
47     }
48 fabstoec 1.8
49 mingyang 1.37 PhotonTools::eScaleCats PhotonTools::EScaleCatHCP(const Photon *p)
50     {
51     if (p->SCluster()->AbsEta()<1.0) {
52     Int_t ieta = p->SCluster()->Seed()->IEta();
53     Int_t iphi = p->SCluster()->Seed()->IPhi();
54     Bool_t central = ( ((std::abs(ieta)-1)/5+1 >= 2 && (std::abs(ieta)-1)/5+1 < 5) || ((std::abs(ieta)-1)/5+1 >= 7 && (std::abs(ieta)-1)/5+1 < 9 ) || ((std::abs(ieta)-1)/5+1 >= 11 && (std::abs(ieta)-1)/5+1 < 13) || ((std::abs(ieta)-1)/5+1 >= 15 && (std::abs(ieta)-1)/5+1 < 17) ) && (iphi %20) > 5 && (iphi%20) <16;
55     if(p->R9()>0.94 && central) return kEBlowEtaGoldCenter;
56     else if(p->R9()>0.94 && (!central)) return kEBlowEtaGoldGap;
57     else if(p->R9()<=0.94 && central) return kEBlowEtaBadCenter;
58     else return kEBlowEtaBadGap;
59     }
60     else if (p->SCluster()->AbsEta()<1.5) {
61     if (p->R9()>0.94) return kEBhighEtaGold;
62     else return kEBhighEtaBad;
63     }
64     else if (p->SCluster()->AbsEta()<2.0) {
65     if (p->R9()>0.94) return kEElowEtaGold;
66     else return kEElowEtaBad;
67     }
68     else {
69     if (p->R9()>0.94) return kEEhighEtaGold;
70     else return kEEhighEtaBad;
71     }
72    
73     }
74    
75 fabstoec 1.8 void PhotonTools::ScalePhoton(Photon* p, Double_t scale) {
76     if( !p ) return;
77     FourVectorM mom = p->Mom();
78     p->SetMom(scale*mom.X(), scale*mom.Y(), scale*mom.Z(), scale*mom.E());
79    
80     }
81    
82     void PhotonTools::SmearPhoton(Photon* p, TRandom3* rng, Double_t width, UInt_t iSeed) {
83    
84     if( !p ) return;
85     if( !rng) return;
86     if( width < 0.) return;
87    
88 fabstoec 1.10 rng->SetSeed(iSeed);
89 fabstoec 1.8 FourVectorM mom = p->Mom();
90     Double_t scale = rng->Gaus(1.,width);
91 fabstoec 1.10
92 fabstoec 1.40 if( scale > 0) {
93 fabstoec 1.8 p->SetMom(scale*mom.X(), scale*mom.Y(), scale*mom.Z(), scale*mom.E());
94 fabstoec 1.40
95     // moved here from SmearPhotonError, in order to being able to
96     // use different smearing for Error and actual Momentum (any reason not to? Josh?)
97     Double_t smear = p->EnergySmearing();
98     p->SetEnergySmearing(TMath::Sqrt(smear*smear+width*width*p->E()*p->E()));
99    
100     }
101    
102 bendavid 1.13 return;
103     }
104 fabstoec 1.8
105 bendavid 1.13 void PhotonTools::SmearPhotonError(Photon* p, Double_t width) {
106    
107     if( !p ) return;
108     if( width < 0.) return;
109    
110 fabstoec 1.40 // WARNING: This is not set correctly if different smearing is used for error and actual momentum...
111     // therefore move this into 'SmearPhoton'.... any reason not to ? Josh ?
112     //Double_t smear = p->EnergySmearing();
113     //p->SetEnergySmearing(TMath::Sqrt(smear*smear+width*width*p->E()*p->E()));
114 bendavid 1.14
115 bendavid 1.13 Double_t err = p->EnergyErrSmeared();
116     if (err>=0.0) {
117     p->SetEnergyErrSmeared(TMath::Sqrt(err*err + width*width*p->E()*p->E()));
118     }
119    
120 fabstoec 1.8 }
121    
122 bendavid 1.16 void PhotonTools::ScalePhotonR9(Photon* p, Double_t scale) {
123     p->SetR9(scale*p->R9());
124     }
125    
126 bendavid 1.23
127 bendavid 1.16 void PhotonTools::ScalePhotonError(Photon* p, Double_t scale) {
128     p->SetEnergyErrSmeared(scale*p->EnergyErrSmeared());
129     }
130    
131 bendavid 1.1 //--------------------------------------------------------------------------------------------------
132     Bool_t PhotonTools::PassConversionId(const Photon *p, const DecayParticle *c) {
133    
134     if (!c) return kTRUE;
135    
136     ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
137     Double_t deta = c->Eta()-dirconvsc.Eta();
138     Double_t dphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
139     Double_t eoverp = p->SCluster()->Energy()/c->P();
140    
141     if (p->IsEB() && eoverp>2.0) return kFALSE;
142     if (p->IsEE() && eoverp>3.0) return kFALSE;
143    
144     if (p->IsEE() && TMath::Abs(deta)>0.01) return kFALSE;
145     if (p->IsEE() && TMath::Abs(dphi)>0.01) return kFALSE;
146    
147     return kTRUE;
148    
149     }
150    
151     //--------------------------------------------------------------------------------------------------
152     Bool_t PhotonTools::PassElectronVeto(const Photon *p, const ElectronCol *els) {
153    
154     Bool_t pass = kTRUE;
155     for (UInt_t i=0; i<els->GetEntries(); ++i) {
156     const Electron *e = els->At(i);
157     if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0) {
158     pass = kFALSE;
159     }
160     }
161    
162     return pass;
163     }
164    
165     //--------------------------------------------------------------------------------------------------
166 bendavid 1.13 const Electron *PhotonTools::MatchedElectron(const Photon *p, const ElectronCol *els) {
167    
168     for (UInt_t i=0; i<els->GetEntries(); ++i) {
169     const Electron *e = els->At(i);
170     if ( e->SCluster()==p->SCluster() ) {
171     return e;
172     }
173     }
174    
175     return 0;
176     }
177    
178     //--------------------------------------------------------------------------------------------------
179     const Photon *PhotonTools::MatchedPhoton(const Electron *e, const PhotonCol *phs) {
180    
181     for (UInt_t i=0; i<phs->GetEntries(); ++i) {
182     const Photon *p = phs->At(i);
183     if ( p->SCluster()==e->SCluster() ) {
184     return p;
185     }
186     }
187    
188     return 0;
189     }
190    
191     //--------------------------------------------------------------------------------------------------
192     const SuperCluster *PhotonTools::MatchedSC(const SuperCluster *psc, const SuperClusterCol *scs, Double_t drMin) {
193    
194     Double_t drsmallest = 999.;
195     const SuperCluster *match = 0;
196     for (UInt_t i=0; i<scs->GetEntries(); ++i) {
197     const SuperCluster *sc = scs->At(i);
198     Double_t dr = MathUtils::DeltaR(*sc,*psc);
199     if ( dr<drsmallest && dr<drMin ) {
200     drsmallest = dr;
201     match = sc;
202     }
203     }
204    
205     return match;
206     }
207    
208     //--------------------------------------------------------------------------------------------------
209 bendavid 1.25 const SuperCluster *PhotonTools::MatchedPFSC(const SuperCluster *psc, const PhotonCol *pfphos, const ElectronCol *eles, Double_t drMin) {
210    
211    
212     if (pfphos) {
213     for (UInt_t i=0; i<pfphos->GetEntries(); ++i) {
214     const Photon *pfpho = pfphos->At(i);
215     if (psc == pfpho->SCluster() && pfpho->PFSCluster()) return pfpho->PFSCluster();
216     }
217     }
218    
219     for (UInt_t i=0; i<eles->GetEntries(); ++i) {
220     const Electron *ele = eles->At(i);
221     if (psc == ele->SCluster() && ele->PFSCluster()) return ele->PFSCluster();
222     }
223    
224     Double_t drsmallest = 999.;
225     const SuperCluster *match = 0;
226     for (UInt_t i=0; i<eles->GetEntries(); ++i) {
227     const Electron *ele = eles->At(i);
228     //if (psc == ele->SCluster() && ele->HasPFSCluster()) return ele->PFSCluster();
229     if (!ele->PFSCluster()) continue;
230    
231     Double_t dr = MathUtils::DeltaR(*ele->PFSCluster(),*psc);
232     if ( dr<drsmallest && dr<drMin ) {
233     drsmallest = dr;
234     match = ele->PFSCluster();
235     }
236     }
237    
238     return match;
239     }
240    
241     //--------------------------------------------------------------------------------------------------
242 fabstoec 1.4 Double_t PhotonTools::ElectronVetoCiC(const Photon *p, const ElectronCol *els) {
243    
244     for (UInt_t i=0; i<els->GetEntries(); ++i) {
245     const Electron *e = els->At(i);
246 mingyang 1.41 //if (MathUtils::DeltaR(*e->SCluster(),*p->SCluster())<1) {
247 fabstoec 1.5 if (e->SCluster()==p->SCluster() ) {
248     //if( e->GsfTrk()->NExpectedHitsInner()==0 && e->GsfTrk()->Pt() > 2.5 ) {
249     if( e->GsfTrk()->NExpectedHitsInner()==0 ) {
250     double dEta = e->DeltaEtaSuperClusterTrackAtVtx();
251     double dPhi = e->DeltaPhiSuperClusterTrackAtVtx();
252     double dR = TMath::Sqrt(dEta*dEta+dPhi*dPhi);
253     return dR;
254     }
255     }
256     }
257     return 99.;
258 fabstoec 1.4 }
259    
260     //--------------------------------------------------------------------------------------------------
261 bendavid 1.1 Bool_t PhotonTools::PassElectronVetoConvRecovery(const Photon *p, const ElectronCol *els, const DecayParticleCol *conversions, const BaseVertex *v) {
262    
263     Bool_t pass = kTRUE;
264     for (UInt_t i=0; i<els->GetEntries(); ++i) {
265     const Electron *e = els->At(i);
266 fabstoec 1.34
267     // HACVK to match CMSSW bug...
268 bendavid 1.1 if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0 && ElectronTools::PassConversionFilter(e, conversions,
269 fabstoec 1.31 v, 0, 0., 1e-6, kTRUE, kFALSE) ) {
270 fabstoec 1.34 // v, 0, 1e-6, 2.0, kFALSE, kFALSE) ) {
271 bendavid 1.1 pass = kFALSE;
272     }
273     }
274    
275     return pass;
276     }
277    
278     //--------------------------------------------------------------------------------------------------
279     Bool_t PhotonTools::PassTriggerMatching(const Photon *p, const TriggerObjectCol *trigobjs)
280     {
281    
282     for (UInt_t i=0; i<trigobjs->GetEntries(); ++i) {
283     const TriggerObject *trigobj = trigobjs->At(i);
284     if (trigobj->TriggerType()==TriggerObject::TriggerCluster || trigobj->TriggerType()==TriggerObject::TriggerElectron || trigobj->TriggerType()==TriggerObject::TriggerPhoton) {
285     if (MathUtils::DeltaR(p->SCluster(),trigobj)<0.3) {
286     return kTRUE;
287     }
288     }
289     }
290    
291     return kFALSE;
292    
293    
294     }
295    
296     //--------------------------------------------------------------------------------------------------
297     const DecayParticle *PhotonTools::MatchedConversion(const Photon *p, const DecayParticleCol *conversions,
298     const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
299     Double_t lxyMin, Double_t dRMin) {
300    
301 bendavid 1.13 return MatchedConversion(p->SCluster(), conversions, vtx, nWrongHitsMax, probMin, lxyMin, dRMin);
302    
303     }
304    
305     //--------------------------------------------------------------------------------------------------
306     const DecayParticle *PhotonTools::MatchedConversion(const SuperCluster *sc, const DecayParticleCol *conversions,
307     const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
308     Double_t lxyMin, Double_t dRMin) {
309    
310 bendavid 1.1 const DecayParticle *match = 0;
311 bendavid 1.3 Double_t rhosmallest = 999.;
312 bendavid 1.1 for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
313     const DecayParticle *c = conversions->At(i);
314 bendavid 1.13 ThreeVector dirconvsc = ThreeVector(sc->Point()) - c->Position();
315 bendavid 1.1 Double_t dr = MathUtils::DeltaR(*c,dirconvsc);
316 bendavid 1.3 Double_t rho = c->Position().Rho();
317     if (dr<dRMin && rho<rhosmallest && c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
318 bendavid 1.1 Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
319     Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
320     if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
321 bendavid 1.3 rhosmallest = rho;
322 bendavid 1.1 match = c;
323     }
324     }
325    
326     }
327    
328     return match;
329    
330     }
331    
332 bendavid 1.2 //--------------------------------------------------------------------------------------------------
333     const DecayParticle *PhotonTools::MatchedConversion(const Track *t, const DecayParticleCol *conversions,
334     const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
335     Double_t lxyMin) {
336    
337     for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
338     const DecayParticle *c = conversions->At(i);
339     if (c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
340     Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
341     Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
342     if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
343     const Track *ct1 = dynamic_cast<const ChargedParticle*>(c->Daughter(0))->Trk();
344     const Track *ct2 = dynamic_cast<const ChargedParticle*>(c->Daughter(1))->Trk();
345     if (t==ct1 || t==ct2) return c;
346     }
347     }
348    
349     }
350    
351     return 0;
352    
353     }
354    
355 bendavid 1.1 PhotonTools::DiphotonR9EtaCats PhotonTools::DiphotonR9EtaCat(const Photon *p1, const Photon *p2) {
356    
357     if (p1->IsEB() && p2->IsEB()) {
358     if (p1->R9()>0.93 && p2->R9()>0.93) return kCat1;
359     else return kCat2;
360    
361     }
362     else {
363     if (p1->R9()>0.93 && p2->R9()>0.93) return kCat3;
364     else return kCat4;
365     }
366    
367     }
368    
369     PhotonTools::DiphotonR9EtaConversionCats PhotonTools::DiphotonR9EtaConversionCat(const Photon *p1, const Photon *p2, const DecayParticleCol *conversions, const BaseVertex *v) {
370    
371     const DecayParticle *conv1 = MatchedConversion(p1, conversions, v);
372     const DecayParticle *conv2 = MatchedConversion(p2, conversions, v);
373    
374     if (p1->IsEB() && p2->IsEB()) {
375     if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat1;
376     else if (conv1||conv2) return kNewCat2;
377     else return kNewCat3;
378    
379     }
380     else {
381     if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat4;
382     else if (conv1||conv2) return kNewCat5;
383     else return kNewCat6;
384     }
385    
386 fabstoec 1.4 }
387    
388     PhotonTools::CiCBaseLineCats PhotonTools::CiCBaseLineCat(const Photon *p) {
389 bendavid 1.12 if( p->SCluster()->AbsEta()<1.5 ) {
390 fabstoec 1.4 if ( p->R9() > 0.94 ) return kCiCCat1;
391     else return kCiCCat2;
392     } else {
393     if ( p->R9() > 0.94 ) return kCiCCat3;
394     else return kCiCCat4;
395     }
396 fabstoec 1.8
397     // shoud NEVER happen, but you never know...
398     return kCiCNoCat;
399 fabstoec 1.4 }
400 fabstoec 1.5
401     //--------------------------------------------------------------------------------------------------
402     const DecayParticle *PhotonTools::MatchedCiCConversion(const Photon *p, const DecayParticleCol *conversions,
403     Double_t dPhiMin,
404 fabstoec 1.7 Double_t dEtaMin,
405     Double_t dRMin,
406 fabstoec 1.39 bool print,
407     int* numLegs,
408     int* convIdx) {
409 fabstoec 1.5
410 fabstoec 1.8 // if there are no conversons, return
411     if ( !p || !conversions) return NULL;
412    
413 fabstoec 1.6 const DecayParticle *match = NULL;
414 fabstoec 1.39 int matchIdx = -1;
415 fabstoec 1.6
416     double minDeta = 999.;
417     double minDphi = 999.;
418 fabstoec 1.7 double minDR = 999.;
419 fabstoec 1.6
420     double phPhi = p->SCluster()->Phi();
421     double phEta = p->SCluster()->Eta();
422 fabstoec 1.39
423 fabstoec 1.7 if(print)
424     std::cout<<" --- conv match photon eta = "<<phEta<<" phi = "<<phPhi<<std::endl;
425 fabstoec 1.39
426 fabstoec 1.7
427 fabstoec 1.5 for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
428     const DecayParticle *c = conversions->At(i);
429 fabstoec 1.39
430 fabstoec 1.7 if(print)
431 fabstoec 1.39 std::cout<< " c "<<i+1<<" pt = "<<c->Pt()<<" with Nlegs = "<<c->NDaughters()<<" "<<c->Position().X()<<" "<<c->Position().Y()<<" "<<c->Position().Z()<<" "<<std::endl;
432 fabstoec 1.34
433 bendavid 1.38 if (c->Pt() < 10. ) continue; // is this refittedPairMomentum?
434 bendavid 1.29 if (c->NDaughters()==2 && c->Prob() < 1e-6) continue;
435 bendavid 1.24
436 bendavid 1.38 ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
437     Double_t dR = MathUtils::DeltaR(dirconvsc,c->Mom());
438 fabstoec 1.7
439     if(dR < minDR) {
440     minDR = dR;
441 fabstoec 1.39 // minDphi = dphi;
442     // minDeta = TMath::Abs(deta);
443 fabstoec 1.7 match = c;
444 fabstoec 1.39 matchIdx = (int) i;
445 fabstoec 1.7
446 fabstoec 1.39 if(print) {
447     //std::cout<<" conv "<<i+1<<" matches with dPhi = "<<minDphi<<" dEta = "<<minDeta<<std::endl;
448     std::cout<<" conv "<<i+1<<" matches with dR = "<<minDR<<std::endl;
449     }
450 fabstoec 1.6 }
451 fabstoec 1.5 }
452    
453 fabstoec 1.7 //if(minDphi < dPhiMin && minDeta < dEtaMin)
454 fabstoec 1.39 if( minDR < dRMin && match ) {
455     if ( numLegs ) (*numLegs) = match->NDaughters();
456     if ( convIdx ) (*convIdx) = matchIdx;
457     if(print)
458     std::cout<<" best conversion is chosen"<<std::endl;
459 fabstoec 1.7 return match;
460 fabstoec 1.39 } else {
461     if(print)
462     std::cout<<" NO conversion is chosen"<<std::endl;
463 fabstoec 1.7 return NULL;
464 fabstoec 1.39 }
465 fabstoec 1.5 }
466    
467 fabstoec 1.9 bool PhotonTools::PassCiCSelection(const Photon* ph, const Vertex* vtx,
468     const TrackCol* trackCol,
469     const ElectronCol* eleCol,
470     const VertexCol* vtxCol,
471 fabstoec 1.5 double rho, double ptmin,
472 fabstoec 1.9 bool applyEleVeto,
473 fabstoec 1.5 bool print, float* kin) {
474 fabstoec 1.9
475 fabstoec 1.5
476     // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
477     float cic4_allcuts_temp_sublead[] = {
478     3.8, 2.2, 1.77, 1.29,
479     11.7, 3.4, 3.9, 1.84,
480     3.5, 2.2, 2.3, 1.45,
481     0.0106, 0.0097, 0.028, 0.027,
482     0.082, 0.062, 0.065, 0.048,
483     0.94, 0.36, 0.94, 0.32,
484     1., 0.062, 0.97, 0.97,
485     1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
486    
487     // cut on Et instead of Pt???
488     Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
489    
490     // compute all relevant observables first
491     double ecalIso3 = ph->EcalRecHitIsoDr03();
492     double ecalIso4 = ph->EcalRecHitIsoDr04();
493     double hcalIso4 = ph->HcalTowerSumEtDr04();
494    
495     unsigned int wVtxInd = 0;
496    
497     double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
498    
499     // track iso only
500     double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
501    
502     // track iso worst vtx
503     double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
504    
505     double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
506     double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
507    
508     double tIso1 = (combIso1) *50./ph->Et();
509     double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
510     //double tIso2 = (combIso2) *50./ph->Et();
511     double tIso3 = (trackIso3)*50./ph->Et();
512    
513     double covIEtaIEta =ph->CoviEtaiEta();
514     double HoE = ph->HadOverEm();
515    
516     double R9 = ph->R9();
517    
518     double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);
519    
520     // check which category it is ...
521     int _tCat = 1;
522     if ( !isbarrel ) _tCat = 3;
523     if ( R9 < 0.94 ) _tCat++;
524    
525     if(print) {
526     std::cout<<" -------------------------- "<<std::endl;
527 fabstoec 1.15
528     std::cout<<" trackIso1 = "<<trackIso1<<std::endl;
529     std::cout<<" trackIso2 = "<<trackIso2<<std::endl;
530     std::cout<<" trackIso3 = "<<trackIso3<<std::endl;
531     std::cout<<" ecalIso3 = "<<ecalIso3<<std::endl;
532     std::cout<<" ecalIso4 = "<<ecalIso4<<std::endl;
533     std::cout<<" hcalIso4 = "<<hcalIso4<<std::endl;
534    
535 fabstoec 1.5 std::cout<<" photon Et = "<<ph->Et()<<std::endl;
536 fabstoec 1.30 std::cout<<" Eta = "<<ph->SCluster()->Eta()<<std::endl;
537 fabstoec 1.15 std::cout<<" HoE = "<<HoE<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+4*4]<<")"<<std::endl;
538     std::cout<<" R9 = "<<R9<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+5*4]<<")"<<std::endl;
539     std::cout<<" dR = "<<dRTrack<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+6*4]<<")"<<std::endl;
540     std::cout<<" iso1 = "<<tIso1<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+0*4]<<")"<<std::endl;
541     std::cout<<" iso2 = "<<tIso2<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+1*4]<<")"<<std::endl;
542     std::cout<<" iso3 = "<<tIso3<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+2*4]<<")"<<std::endl;
543 fabstoec 1.5 }
544    
545     if(kin) {
546     kin[0] = tIso1;
547     kin[1] = tIso2;
548     kin[2] = tIso3;
549     kin[3] = covIEtaIEta;
550     kin[4] = HoE;
551     kin[5] = R9;
552     kin[6] = dRTrack;
553     kin[7] = (float) ph->Pt();
554     kin[8] = (float) ph->Eta();
555     kin[9] = (float) ph->Phi();
556    
557     // iso quantities separate
558     kin[10] = ecalIso3;
559     kin[11] = ecalIso4;
560     kin[12] = hcalIso4;
561     kin[13] = trackIso1;
562     kin[14] = trackIso2;
563     kin[15] = trackIso3;
564    
565     kin[16] = (float) ph->Et();
566     kin[17] = (float) ph->E();
567    
568 bendavid 1.11 kin[18] = (float) ph->SCluster()->Eta();
569     kin[19] = (float) ph->SCluster()->Phi();
570 fabstoec 1.5 }
571    
572     float passCuts = 1.;
573    
574     if ( ph->Pt() <= ptmin ) passCuts = -1.;
575    
576     // not needed anymore, do in pre-selection...
577     //if ( ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
578    
579     if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
580     && tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
581     && tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
582     && covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
583     && HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
584     && R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4]
585 fabstoec 1.9 && ( dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] || !applyEleVeto ) ) ) passCuts = -1.;
586 fabstoec 1.5
587 fabstoec 1.15 if(print) std::cout<<" ---> "<<passCuts<<std::endl;
588    
589 fabstoec 1.5 if(kin) {
590 bendavid 1.11 kin[20] = passCuts;
591     kin[21] = (float) _tCat;
592 fabstoec 1.5 }
593    
594     if(passCuts > 0.) return true;
595     return false;
596     }
597 bendavid 1.11
598 bendavid 1.27 bool PhotonTools::PassCiCPFIsoSelection(const Photon* ph,
599 fabstoec 1.30 const Vertex* vtx,
600     const PFCandidateCol* pfCol,
601     const VertexCol* vtxCol,
602     double rho, double ptmin,
603     std::vector<double>* kin // store variables for debugging...
604     ) {
605 bendavid 1.27
606     // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
607     float cic4_allcuts_temp_sublead[] = {
608     6.0, 4.7, 5.6, 3.6,
609     10.0, 6.5, 5.6, 4.4,
610     3.8, 2.5, 3.1, 2.2,
611     0.0108, 0.0102, 0.028, 0.028,
612     0.124, 0.092, 0.142, 0.063,
613 fabstoec 1.30 0.94, 0.28, 0.94, 0.24 };
614 bendavid 1.27
615     // cut on Et instead of Pt???
616     Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
617 fabstoec 1.39
618 bendavid 1.27
619     // compute all relevant observables first
620     double ecalIso3 = IsolationTools::PFGammaIsolation(ph, 0.3, 0.0, pfCol);
621     double ecalIso4 = IsolationTools::PFGammaIsolation(ph, 0.4, 0.0, pfCol);
622    
623     unsigned int wVtxInd = 0;
624    
625     double trackIsoSel03 = IsolationTools::PFChargedIsolation(ph, vtx, 0.3, 0.0, pfCol);
626    
627     // track iso worst vtx
628     double trackIsoWorst04 = IsolationTools::PFChargedIsolation(ph, vtx, 0.4, 0.00, pfCol, &wVtxInd, vtxCol);
629    
630 fabstoec 1.30 double combIso1 = ecalIso3+trackIsoSel03 + 2.5 - 0.09*rho;
631     double combIso2 = ecalIso4+trackIsoWorst04 + 2.5 - 0.23*rho;
632 bendavid 1.27
633     double tIso1 = (combIso1) *50./ph->Et();
634     double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
635     //double tIso2 = (combIso2) *50./ph->Et();
636     double tIso3 = (trackIsoSel03)*50./ph->Et();
637    
638     double covIEtaIEta =ph->CoviEtaiEta();
639     double HoE = ph->HadOverEm();
640    
641     double R9 = ph->R9();
642    
643     // check which category it is ...
644     int _tCat = 1;
645     if ( !isbarrel ) _tCat = 3;
646     if ( R9 < 0.94 ) _tCat++;
647    
648     float passCuts = 1.;
649    
650 fabstoec 1.30 if( kin ) {
651     kin->resize(0);
652    
653     kin->push_back(tIso1);
654     kin->push_back(tIso2);
655     kin->push_back(tIso3);
656     kin->push_back(covIEtaIEta);
657     kin->push_back(HoE);
658     kin->push_back(R9);
659    
660     kin->push_back( (double) wVtxInd );
661     kin->push_back( ecalIso3 );
662     kin->push_back( ecalIso4 );
663    
664     kin->push_back( trackIsoSel03 );
665     kin->push_back( trackIsoWorst04 );
666    
667     kin->push_back( combIso1 );
668     kin->push_back( combIso2 );
669     }
670    
671    
672 bendavid 1.27 if ( ph->Pt() <= ptmin ) passCuts = -1.;
673    
674     // not needed anymore, do in pre-selection...
675     //if ( ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
676    
677     if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
678     && tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
679     && tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
680     && covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
681     && HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
682     && R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4] ) ) passCuts = -1.;
683    
684    
685     if(passCuts > 0.) return true;
686     return false;
687     }
688    
689    
690    
691    
692    
693 bendavid 1.13 const MCParticle *PhotonTools::MatchMC(const Particle *ph, const MCParticleCol *c, Bool_t matchElectrons) {
694    
695     // printf("Start loop\n");
696 bendavid 1.11 for (UInt_t i=0; i<c->GetEntries(); ++i) {
697     const MCParticle *p = c->At(i);
698 bendavid 1.13 // if (p->IsGenerated() && p->AbsPdgId()==11 && (p->DistinctMother()->AbsPdgId()==23|| p->DistinctMother()->AbsPdgId()==24)) {
699     // printf("pdgid = %i, status = %i, pt = %5f, mother pdgid = %i\n",p->PdgId(),p->Status(),p->Pt(),p->Mother()->PdgId());
700 bendavid 1.12 // }
701     if (matchElectrons && p->AbsPdgId()==11 && p->IsGenerated() && p->Status()==3 && MathUtils::DeltaR(*ph,*p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==23 || p->Mother()->AbsPdgId()==24 || p->Mother()->AbsPdgId()==22)) {
702     return p;
703     }
704 bendavid 1.11 if ( p->AbsPdgId()==22 && p->IsGenerated() && MathUtils::DeltaR(*ph,*p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==25 || p->Mother()->AbsPdgId()<=21) ) {
705     return p;
706     }
707     }
708     return 0;
709 bendavid 1.12 }
710    
711     PhotonTools::DiphotonR9EtaPtCats PhotonTools::DiphotonR9EtaPtCat(const Photon *p1, const Photon *p2) {
712    
713     PhotonTools::CiCBaseLineCats cat1 = CiCBaseLineCat(p1);
714     PhotonTools::CiCBaseLineCats cat2 = CiCBaseLineCat(p2);
715    
716     PhotonTools::DiphotonR9EtaPtCats evtcat = PhotonTools::kOctCat0;
717    
718     bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
719     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
720    
721     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
722     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
723    
724     if( ph1IsEB && ph2IsEB )
725     evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat0 : PhotonTools::kOctCat1);
726     else
727     evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat2 : PhotonTools::kOctCat3);
728    
729     float ptgg = (p1->Mom()+p2->Mom()).Pt();
730     if (ptgg<40.0) evtcat = PhotonTools::DiphotonR9EtaPtCats(evtcat + 4);
731    
732    
733     return evtcat;
734 fabstoec 1.15 }
735 mingyang 1.17
736 bendavid 1.22 Bool_t PhotonTools::PassSinglePhotonPresel(const Photon *p,const ElectronCol *els, const DecayParticleCol *conversions, const BaseVertex *bs, const TrackCol* trackCol,const Vertex *vtx, double rho, Bool_t applyElectronVeto, Bool_t invertElectronVeto) {
737 mingyang 1.17 float ScEta=p->SCluster()->Eta();
738     float Et=p->Et();
739     float R9=p->R9();
740     float HoE=p->HadOverEm();
741     float CovIEtaIEta=p->CoviEtaiEta();
742     float EcalIsoDr03=p->EcalRecHitIsoDr03();
743     float HcalIsoDr03=p->HcalTowerSumEtDr03();
744     float TrkIsoHollowDr03=p->HollowConeTrkIsoDr03();
745     float NewEcalIso=EcalIsoDr03-0.012*Et;
746     float NewHcalIso=HcalIsoDr03-0.005*Et;
747     float NewTrkIsoHollowDr03=TrkIsoHollowDr03-0.002*Et;
748     Bool_t IsBarrel=kFALSE;
749     Bool_t IsEndcap=kFALSE;
750 bendavid 1.22 Bool_t PassEleVetoRaw = PhotonTools::PassElectronVetoConvRecovery(p, els, conversions, bs);
751 bendavid 1.20 Bool_t PassEleVeto = (!applyElectronVeto && !invertElectronVeto) || (applyElectronVeto && !invertElectronVeto && PassEleVetoRaw) || (!applyElectronVeto && invertElectronVeto && !PassEleVetoRaw);
752 bendavid 1.22 float AbsTrackIsoCIC=IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0,trackCol, NULL, NULL, (!applyElectronVeto ? els : NULL) );
753 mingyang 1.18 float HcalEcalPUCorr=EcalIsoDr03+HcalIsoDr03-0.17*rho;
754    
755 mingyang 1.17 if(fabs(ScEta)<1.4442){IsBarrel=kTRUE;}
756     if(fabs(ScEta)>1.566 && fabs(ScEta)<2.5){IsEndcap=kTRUE;}
757     if((!IsBarrel) && (!IsEndcap)){
758     return kFALSE;
759     }
760     if(!PassEleVeto){
761     return kFALSE;
762     }
763     if(R9<=0.9){
764 bendavid 1.21 if(HoE<0.075 && ((IsBarrel && CovIEtaIEta<0.014) || (IsEndcap && CovIEtaIEta<0.034)) && NewEcalIso<4 && NewHcalIso<4 && NewTrkIsoHollowDr03<4 && HcalEcalPUCorr<3 && AbsTrackIsoCIC<2.8 && TrkIsoHollowDr03<4.0) {
765 mingyang 1.17 return kTRUE;
766     }
767     }
768     if(R9>0.9){
769 bendavid 1.21 if(((IsBarrel && HoE<0.082 && CovIEtaIEta<0.014) || (IsEndcap && HoE <0.075 && CovIEtaIEta<0.034)) && NewEcalIso<50 && NewHcalIso<50 && NewTrkIsoHollowDr03<50 && HcalEcalPUCorr<3 && AbsTrackIsoCIC<2.8 && TrkIsoHollowDr03<4.0) {
770 mingyang 1.17 return kTRUE;
771     }
772     }
773     return kFALSE;
774     }
775 mingyang 1.26
776 fabstoec 1.36
777     Bool_t PhotonTools::PassSinglePhotonPreselPFISO_NoTrigger(const Photon *p,const ElectronCol *els, const DecayParticleCol *conversions, const BaseVertex *bs, const TrackCol* trackCol,const Vertex *vtx, double rho, const PFCandidateCol *fPFCands, Bool_t applyElectronVeto, Bool_t invertElectronVeto) {
778    
779     float ScEta=p->SCluster()->Eta();
780     float Et=p->Et();
781     float R9=p->R9();
782     float HoE=p->HadOverEm();
783     float CovIEtaIEta=p->CoviEtaiEta();
784     float EcalIsoDr03=p->EcalRecHitIsoDr03();
785     float HcalIsoDr03=p->HcalTowerSumEtDr03();
786     float TrkIsoHollowDr03=p->HollowConeTrkIsoDr03();
787    
788     float NewEcalIso=EcalIsoDr03-0.012*Et;
789     float NewHcalIso=HcalIsoDr03-0.005*Et;
790     float NewTrkIsoHollowDr03=TrkIsoHollowDr03-0.002*Et;
791    
792    
793     Bool_t IsBarrel=kFALSE;
794     Bool_t IsEndcap=kFALSE;
795     Bool_t PassEleVetoRaw = PhotonTools::PassElectronVetoConvRecovery(p, els, conversions, bs);
796     Bool_t PassEleVeto = (!applyElectronVeto && !invertElectronVeto) || (applyElectronVeto && !invertElectronVeto && PassEleVetoRaw) || (!applyElectronVeto && invertElectronVeto && !PassEleVetoRaw);
797     float ChargedIso_selvtx_DR002To0p02=IsolationTools::PFChargedIsolation(p,vtx, 0.2, 0.,fPFCands);
798     if(fabs(ScEta)<1.4442){IsBarrel=kTRUE;}
799     if(fabs(ScEta)>1.566 && fabs(ScEta)<2.5){IsEndcap=kTRUE;}
800     if((!IsBarrel) && (!IsEndcap)){
801     return kFALSE;
802     }
803     if(!PassEleVeto){
804     return kFALSE;
805     }
806     if( ((IsBarrel && CovIEtaIEta<0.02) || (IsEndcap && CovIEtaIEta<0.06)) && ChargedIso_selvtx_DR002To0p02<4) {
807     return kTRUE;
808     }
809    
810     return kFALSE;
811     }
812    
813    
814    
815 mingyang 1.26 Bool_t PhotonTools::PassSinglePhotonPreselPFISO(const Photon *p,const ElectronCol *els, const DecayParticleCol *conversions, const BaseVertex *bs, const TrackCol* trackCol,const Vertex *vtx, double rho, const PFCandidateCol *fPFCands, Bool_t applyElectronVeto, Bool_t invertElectronVeto) {
816 fabstoec 1.34
817    
818 mingyang 1.26 float ScEta=p->SCluster()->Eta();
819     float Et=p->Et();
820     float R9=p->R9();
821     float HoE=p->HadOverEm();
822     float CovIEtaIEta=p->CoviEtaiEta();
823     float EcalIsoDr03=p->EcalRecHitIsoDr03();
824     float HcalIsoDr03=p->HcalTowerSumEtDr03();
825     float TrkIsoHollowDr03=p->HollowConeTrkIsoDr03();
826 fabstoec 1.34
827 mingyang 1.26 float NewEcalIso=EcalIsoDr03-0.012*Et;
828     float NewHcalIso=HcalIsoDr03-0.005*Et;
829     float NewTrkIsoHollowDr03=TrkIsoHollowDr03-0.002*Et;
830 fabstoec 1.34
831    
832 mingyang 1.26 Bool_t IsBarrel=kFALSE;
833     Bool_t IsEndcap=kFALSE;
834     Bool_t PassEleVetoRaw = PhotonTools::PassElectronVetoConvRecovery(p, els, conversions, bs);
835     Bool_t PassEleVeto = (!applyElectronVeto && !invertElectronVeto) || (applyElectronVeto && !invertElectronVeto && PassEleVetoRaw) || (!applyElectronVeto && invertElectronVeto && !PassEleVetoRaw);
836 bendavid 1.27 float ChargedIso_selvtx_DR002To0p02=IsolationTools::PFChargedIsolation(p,vtx, 0.2, 0.,fPFCands);
837 mingyang 1.26 if(fabs(ScEta)<1.4442){IsBarrel=kTRUE;}
838     if(fabs(ScEta)>1.566 && fabs(ScEta)<2.5){IsEndcap=kTRUE;}
839     if((!IsBarrel) && (!IsEndcap)){
840     return kFALSE;
841     }
842     if(!PassEleVeto){
843     return kFALSE;
844     }
845     if(R9<=0.9){
846     if(HoE<0.075 && ((IsBarrel && CovIEtaIEta<0.014) || (IsEndcap && CovIEtaIEta<0.034)) && NewEcalIso<4 && NewHcalIso<4 && NewTrkIsoHollowDr03<4 && ChargedIso_selvtx_DR002To0p02<4) {
847     return kTRUE;
848     }
849     }
850     if(R9>0.9){
851     if(((IsBarrel && HoE<0.082 && CovIEtaIEta<0.014) || (IsEndcap && HoE <0.075 && CovIEtaIEta<0.034)) && NewEcalIso<50 && NewHcalIso<50 && NewTrkIsoHollowDr03<50 && ChargedIso_selvtx_DR002To0p02<4) {
852     return kTRUE;
853     }
854     }
855     return kFALSE;
856     }
857 fabstoec 1.32
858 fabstoec 1.34 bool PhotonTools::PassVgamma2011Selection(const Photon* ph, double rho) {
859    
860     bool isEB = (ph->SCluster()->AbsEta()<1.5);
861    
862     if (ph->HasPixelSeed()) return false; // ? is this what we want ?
863     if (ph->HadOverEm() > 0.05) return false;
864     if (ph->CoviEtaiEta() > (isEB ? 0.011 : 0.03) ) return false;
865     if (ph->HollowConeTrkIsoDr04() > ( 2.0 + 0.001 *ph->Pt() +
866     (isEB ? 0.0167 : 0.032)*rho )) return false;
867     if (ph->EcalRecHitIsoDr04() > ( 4.2 + 0.006 *ph->Pt() +
868     (isEB ? 0.183 : 0.090)*rho )) return false;
869     if (ph->HcalTowerSumEtDr04() > ( 2.2 + 0.0025*ph->Pt() +
870     (isEB ? 0.062 : 0.180)*rho )) return false;
871    
872     // spike cleaning...
873 fabstoec 1.35 if ( ph->CoviEtaiEta() < 0.001 && TMath::Sqrt(TMath::Abs(ph->SCluster()->Seed()->CoviPhiiPhi())) < 0.001)
874     return false;
875 fabstoec 1.34
876     return true;
877     }
878 fabstoec 1.32
879 mingyang 1.41
880     Bool_t PhotonTools::PassSinglePhotonPreselPFISONoEcal(const Photon *p,const ElectronCol *els, const DecayParticleCol *conversions, const BaseVertex *bs, const TrackCol* trackCol,const Vertex *vtx, double rho, const PFCandidateCol *fPFCands, Bool_t applyElectronVeto, Bool_t invertElectronVeto) {
881    
882    
883     float ScEta=p->SCluster()->Eta();
884     float Et=p->Et();
885     float R9=p->R9();
886     float HoE=p->HadOverEm();
887     float CovIEtaIEta=p->CoviEtaiEta();
888     float HcalIsoDr03=p->HcalTowerSumEtDr03();
889     float TrkIsoHollowDr03=p->HollowConeTrkIsoDr03();
890    
891     float NewHcalIso=HcalIsoDr03-0.005*Et;
892     float NewTrkIsoHollowDr03=TrkIsoHollowDr03-0.002*Et;
893    
894    
895     Bool_t IsBarrel=kFALSE;
896     Bool_t IsEndcap=kFALSE;
897     Bool_t PassEleVetoRaw = PhotonTools::PassElectronVetoConvRecovery(p, els, conversions, bs);
898     Bool_t PassEleVeto = (!applyElectronVeto && !invertElectronVeto) || (applyElectronVeto && !invertElectronVeto && PassEleVetoRaw) || (!applyElectronVeto && invertElectronVeto && !PassEleVetoRaw);
899     float ChargedIso_selvtx_DR002To0p02=IsolationTools::PFChargedIsolation(p,vtx, 0.2, 0.,fPFCands);
900     if(fabs(ScEta)<1.4442){IsBarrel=kTRUE;}
901     if(fabs(ScEta)>1.566 && fabs(ScEta)<2.5){IsEndcap=kTRUE;}
902     if((!IsBarrel) && (!IsEndcap)){
903     return kFALSE;
904     }
905     if(!PassEleVeto){
906     return kFALSE;
907     }
908     if(R9<=0.9){
909     if(HoE<0.075 && ((IsBarrel && CovIEtaIEta<0.014) || (IsEndcap && CovIEtaIEta<0.034)) && NewHcalIso<4 && NewTrkIsoHollowDr03<4 && ChargedIso_selvtx_DR002To0p02<4) {
910     return kTRUE;
911     }
912     }
913     if(R9>0.9){
914     if(((IsBarrel && HoE<0.082 && CovIEtaIEta<0.014) || (IsEndcap && HoE <0.075 && CovIEtaIEta<0.034)) && NewHcalIso<50 && NewTrkIsoHollowDr03<50 && ChargedIso_selvtx_DR002To0p02<4) {
915     return kTRUE;
916     }
917     }
918     return kFALSE;
919     }
920    
921 fabstoec 1.32 void PhotonTools::ScalePhotonShowerShapes(Photon* p, PhotonTools::ShowerShapeScales scale) {
922    
923     switch(scale) {
924    
925     case kNoShowerShapeScaling:
926     break;
927    
928     case k2011ShowerShape:
929     //regression sigmaE
930     if (p->SCluster()->AbsEta()<1.5)
931     PhotonTools::ScalePhotonError(p,1.07);
932     else
933     PhotonTools::ScalePhotonError(p,1.045);
934    
935     //R9
936     if (p->SCluster()->AbsEta()<1.5) p->SetR9(1.0035*p->R9());
937     else p->SetR9(1.0035*p->R9());
938     //CoviEtaiEta(SigiEtaiEta)
939     if (p->SCluster()->AbsEta()<1.5) p->SetCoviEtaiEta(0.87*p->CoviEtaiEta()+0.0011);
940     else p->SetCoviEtaiEta(0.99*p->CoviEtaiEta());
941     //EtaWidth
942     if (p->SCluster()->AbsEta()<1.5) p->SetEtaWidth(0.99*p->EtaWidth());
943     else p->SetEtaWidth(0.99*p->EtaWidth());
944     //PhiWidth
945     if (p->SCluster()->AbsEta()<1.5) p->SetPhiWidth(0.99*p->PhiWidth());
946     else p->SetPhiWidth(0.99*p->PhiWidth());
947     break;
948    
949 bendavid 1.33 case k2012ShowerShape:
950 fabstoec 1.32 //R9
951     if (p->SCluster()->AbsEta()<1.5) p->SetR9(1.0045*p->R9()+0.001);
952     else p->SetR9(1.0086*p->R9()-0.0007);
953     //CoviEtaiEta(SigiEtaiEta)
954     if (p->SCluster()->AbsEta()<1.5) p->SetCoviEtaiEta(0.891832*p->CoviEtaiEta()+0.0009133);
955     else p->SetCoviEtaiEta(0.9947*p->CoviEtaiEta()+0.00003);
956     //EtaWidth
957     if (p->SCluster()->AbsEta()<1.5) p->SetEtaWidth(1.04302*p->EtaWidth()- 0.000618);
958     else p->SetEtaWidth(0.903254*p->EtaWidth()+ 0.001346);
959     //PhiWidth
960     if (p->SCluster()->AbsEta()<1.5) p->SetPhiWidth(1.00002*p->PhiWidth()- 0.000371);
961     else p->SetPhiWidth(0.99992*p->PhiWidth()+ 4.8e-07);
962     //s4Ratio
963     if (p->SCluster()->AbsEta()<1.5) p->SetS4Ratio(1.01894*p->S4Ratio()-0.01034);
964     else p->SetS4Ratio(1.04969*p->S4Ratio()-0.03642);
965     break;
966     default:
967     break;
968    
969     }
970     return;
971     }
972