ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/PhotonTools.cc
Revision: 1.37
Committed: Mon Oct 8 17:37:40 2012 UTC (12 years, 7 months ago) by mingyang
Content type: text/plain
Branch: MAIN
Changes since 1.36: +27 -1 lines
Log Message: 
new photon cat added

File Contents

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