ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/PhotonTools.cc
Revision: 1.25
Committed: Fri May 18 16:20:33 2012 UTC (12 years, 11 months ago) by bendavid
Content type: text/plain
Branch: MAIN
Changes since 1.24: +34 -1 lines
Log Message:
add pf photon matching tools

File Contents

# User Rev Content
1 bendavid 1.25 // $Id: PhotonTools.cc,v 1.24 2012/03/22 16:25:19 bendavid 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     void PhotonTools::ScalePhoton(Photon* p, Double_t scale) {
50     if( !p ) return;
51     FourVectorM mom = p->Mom();
52     p->SetMom(scale*mom.X(), scale*mom.Y(), scale*mom.Z(), scale*mom.E());
53    
54     }
55    
56     void PhotonTools::SmearPhoton(Photon* p, TRandom3* rng, Double_t width, UInt_t iSeed) {
57    
58     if( !p ) return;
59     if( !rng) return;
60     if( width < 0.) return;
61    
62 fabstoec 1.10 rng->SetSeed(iSeed);
63 fabstoec 1.8 FourVectorM mom = p->Mom();
64     Double_t scale = rng->Gaus(1.,width);
65 fabstoec 1.10
66 fabstoec 1.8 if( scale > 0)
67     p->SetMom(scale*mom.X(), scale*mom.Y(), scale*mom.Z(), scale*mom.E());
68 bendavid 1.13
69     return;
70     }
71 fabstoec 1.8
72 bendavid 1.13 void PhotonTools::SmearPhotonError(Photon* p, Double_t width) {
73    
74     if( !p ) return;
75     if( width < 0.) return;
76    
77 bendavid 1.14
78     Double_t smear = p->EnergySmearing();
79     p->SetEnergySmearing(TMath::Sqrt(smear*smear+width*width*p->E()*p->E()));
80    
81 bendavid 1.13 Double_t err = p->EnergyErrSmeared();
82     if (err>=0.0) {
83     p->SetEnergyErrSmeared(TMath::Sqrt(err*err + width*width*p->E()*p->E()));
84     }
85    
86 fabstoec 1.8 }
87    
88 bendavid 1.16 void PhotonTools::ScalePhotonR9(Photon* p, Double_t scale) {
89     p->SetR9(scale*p->R9());
90     }
91    
92 bendavid 1.23
93 bendavid 1.16 void PhotonTools::ScalePhotonError(Photon* p, Double_t scale) {
94     p->SetEnergyErrSmeared(scale*p->EnergyErrSmeared());
95     }
96    
97 bendavid 1.1 //--------------------------------------------------------------------------------------------------
98     Bool_t PhotonTools::PassConversionId(const Photon *p, const DecayParticle *c) {
99    
100     if (!c) return kTRUE;
101    
102     ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
103     Double_t deta = c->Eta()-dirconvsc.Eta();
104     Double_t dphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
105     Double_t eoverp = p->SCluster()->Energy()/c->P();
106    
107     if (p->IsEB() && eoverp>2.0) return kFALSE;
108     if (p->IsEE() && eoverp>3.0) return kFALSE;
109    
110     if (p->IsEE() && TMath::Abs(deta)>0.01) return kFALSE;
111     if (p->IsEE() && TMath::Abs(dphi)>0.01) return kFALSE;
112    
113     return kTRUE;
114    
115     }
116    
117     //--------------------------------------------------------------------------------------------------
118     Bool_t PhotonTools::PassElectronVeto(const Photon *p, const ElectronCol *els) {
119    
120     Bool_t pass = kTRUE;
121     for (UInt_t i=0; i<els->GetEntries(); ++i) {
122     const Electron *e = els->At(i);
123     if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0) {
124     pass = kFALSE;
125     }
126     }
127    
128     return pass;
129     }
130    
131     //--------------------------------------------------------------------------------------------------
132 bendavid 1.13 const Electron *PhotonTools::MatchedElectron(const Photon *p, const ElectronCol *els) {
133    
134     for (UInt_t i=0; i<els->GetEntries(); ++i) {
135     const Electron *e = els->At(i);
136     if ( e->SCluster()==p->SCluster() ) {
137     return e;
138     }
139     }
140    
141     return 0;
142     }
143    
144     //--------------------------------------------------------------------------------------------------
145     const Photon *PhotonTools::MatchedPhoton(const Electron *e, const PhotonCol *phs) {
146    
147     for (UInt_t i=0; i<phs->GetEntries(); ++i) {
148     const Photon *p = phs->At(i);
149     if ( p->SCluster()==e->SCluster() ) {
150     return p;
151     }
152     }
153    
154     return 0;
155     }
156    
157     //--------------------------------------------------------------------------------------------------
158     const SuperCluster *PhotonTools::MatchedSC(const SuperCluster *psc, const SuperClusterCol *scs, Double_t drMin) {
159    
160     Double_t drsmallest = 999.;
161     const SuperCluster *match = 0;
162     for (UInt_t i=0; i<scs->GetEntries(); ++i) {
163     const SuperCluster *sc = scs->At(i);
164     Double_t dr = MathUtils::DeltaR(*sc,*psc);
165     if ( dr<drsmallest && dr<drMin ) {
166     drsmallest = dr;
167     match = sc;
168     }
169     }
170    
171     return match;
172     }
173    
174     //--------------------------------------------------------------------------------------------------
175 bendavid 1.25 const SuperCluster *PhotonTools::MatchedPFSC(const SuperCluster *psc, const PhotonCol *pfphos, const ElectronCol *eles, Double_t drMin) {
176    
177    
178     if (pfphos) {
179     for (UInt_t i=0; i<pfphos->GetEntries(); ++i) {
180     const Photon *pfpho = pfphos->At(i);
181     if (psc == pfpho->SCluster() && pfpho->PFSCluster()) return pfpho->PFSCluster();
182     }
183     }
184    
185     for (UInt_t i=0; i<eles->GetEntries(); ++i) {
186     const Electron *ele = eles->At(i);
187     if (psc == ele->SCluster() && ele->PFSCluster()) return ele->PFSCluster();
188     }
189    
190     Double_t drsmallest = 999.;
191     const SuperCluster *match = 0;
192     for (UInt_t i=0; i<eles->GetEntries(); ++i) {
193     const Electron *ele = eles->At(i);
194     //if (psc == ele->SCluster() && ele->HasPFSCluster()) return ele->PFSCluster();
195     if (!ele->PFSCluster()) continue;
196    
197     Double_t dr = MathUtils::DeltaR(*ele->PFSCluster(),*psc);
198     if ( dr<drsmallest && dr<drMin ) {
199     drsmallest = dr;
200     match = ele->PFSCluster();
201     }
202     }
203    
204     return match;
205     }
206    
207     //--------------------------------------------------------------------------------------------------
208 fabstoec 1.4 Double_t PhotonTools::ElectronVetoCiC(const Photon *p, const ElectronCol *els) {
209    
210     for (UInt_t i=0; i<els->GetEntries(); ++i) {
211     const Electron *e = els->At(i);
212 fabstoec 1.5 if (e->SCluster()==p->SCluster() ) {
213     //if( e->GsfTrk()->NExpectedHitsInner()==0 && e->GsfTrk()->Pt() > 2.5 ) {
214     if( e->GsfTrk()->NExpectedHitsInner()==0 ) {
215     double dEta = e->DeltaEtaSuperClusterTrackAtVtx();
216     double dPhi = e->DeltaPhiSuperClusterTrackAtVtx();
217     double dR = TMath::Sqrt(dEta*dEta+dPhi*dPhi);
218     return dR;
219     }
220     }
221     }
222     return 99.;
223 fabstoec 1.4 }
224    
225     //--------------------------------------------------------------------------------------------------
226 bendavid 1.1 Bool_t PhotonTools::PassElectronVetoConvRecovery(const Photon *p, const ElectronCol *els, const DecayParticleCol *conversions, const BaseVertex *v) {
227    
228     Bool_t pass = kTRUE;
229     for (UInt_t i=0; i<els->GetEntries(); ++i) {
230     const Electron *e = els->At(i);
231     if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0 && ElectronTools::PassConversionFilter(e, conversions,
232     v, 0, 1e-6, 2.0, kFALSE, kFALSE) ) {
233     pass = kFALSE;
234     }
235     }
236    
237     return pass;
238     }
239    
240     //--------------------------------------------------------------------------------------------------
241     Bool_t PhotonTools::PassTriggerMatching(const Photon *p, const TriggerObjectCol *trigobjs)
242     {
243    
244     for (UInt_t i=0; i<trigobjs->GetEntries(); ++i) {
245     const TriggerObject *trigobj = trigobjs->At(i);
246     if (trigobj->TriggerType()==TriggerObject::TriggerCluster || trigobj->TriggerType()==TriggerObject::TriggerElectron || trigobj->TriggerType()==TriggerObject::TriggerPhoton) {
247     if (MathUtils::DeltaR(p->SCluster(),trigobj)<0.3) {
248     return kTRUE;
249     }
250     }
251     }
252    
253     return kFALSE;
254    
255    
256     }
257    
258     //--------------------------------------------------------------------------------------------------
259     const DecayParticle *PhotonTools::MatchedConversion(const Photon *p, const DecayParticleCol *conversions,
260     const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
261     Double_t lxyMin, Double_t dRMin) {
262    
263 bendavid 1.13 return MatchedConversion(p->SCluster(), conversions, vtx, nWrongHitsMax, probMin, lxyMin, dRMin);
264    
265     }
266    
267     //--------------------------------------------------------------------------------------------------
268     const DecayParticle *PhotonTools::MatchedConversion(const SuperCluster *sc, const DecayParticleCol *conversions,
269     const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
270     Double_t lxyMin, Double_t dRMin) {
271    
272 bendavid 1.1 const DecayParticle *match = 0;
273 bendavid 1.3 Double_t rhosmallest = 999.;
274 bendavid 1.1 for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
275     const DecayParticle *c = conversions->At(i);
276 bendavid 1.13 ThreeVector dirconvsc = ThreeVector(sc->Point()) - c->Position();
277 bendavid 1.1 Double_t dr = MathUtils::DeltaR(*c,dirconvsc);
278 bendavid 1.3 Double_t rho = c->Position().Rho();
279     if (dr<dRMin && rho<rhosmallest && c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
280 bendavid 1.1 Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
281     Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
282     if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
283 bendavid 1.3 rhosmallest = rho;
284 bendavid 1.1 match = c;
285     }
286     }
287    
288     }
289    
290     return match;
291    
292     }
293    
294 bendavid 1.2 //--------------------------------------------------------------------------------------------------
295     const DecayParticle *PhotonTools::MatchedConversion(const Track *t, const DecayParticleCol *conversions,
296     const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
297     Double_t lxyMin) {
298    
299     for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
300     const DecayParticle *c = conversions->At(i);
301     if (c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
302     Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
303     Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
304     if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
305     const Track *ct1 = dynamic_cast<const ChargedParticle*>(c->Daughter(0))->Trk();
306     const Track *ct2 = dynamic_cast<const ChargedParticle*>(c->Daughter(1))->Trk();
307     if (t==ct1 || t==ct2) return c;
308     }
309     }
310    
311     }
312    
313     return 0;
314    
315     }
316    
317 bendavid 1.1 PhotonTools::DiphotonR9EtaCats PhotonTools::DiphotonR9EtaCat(const Photon *p1, const Photon *p2) {
318    
319     if (p1->IsEB() && p2->IsEB()) {
320     if (p1->R9()>0.93 && p2->R9()>0.93) return kCat1;
321     else return kCat2;
322    
323     }
324     else {
325     if (p1->R9()>0.93 && p2->R9()>0.93) return kCat3;
326     else return kCat4;
327     }
328    
329     }
330    
331     PhotonTools::DiphotonR9EtaConversionCats PhotonTools::DiphotonR9EtaConversionCat(const Photon *p1, const Photon *p2, const DecayParticleCol *conversions, const BaseVertex *v) {
332    
333     const DecayParticle *conv1 = MatchedConversion(p1, conversions, v);
334     const DecayParticle *conv2 = MatchedConversion(p2, conversions, v);
335    
336     if (p1->IsEB() && p2->IsEB()) {
337     if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat1;
338     else if (conv1||conv2) return kNewCat2;
339     else return kNewCat3;
340    
341     }
342     else {
343     if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat4;
344     else if (conv1||conv2) return kNewCat5;
345     else return kNewCat6;
346     }
347    
348 fabstoec 1.4 }
349    
350     PhotonTools::CiCBaseLineCats PhotonTools::CiCBaseLineCat(const Photon *p) {
351 bendavid 1.12 if( p->SCluster()->AbsEta()<1.5 ) {
352 fabstoec 1.4 if ( p->R9() > 0.94 ) return kCiCCat1;
353     else return kCiCCat2;
354     } else {
355     if ( p->R9() > 0.94 ) return kCiCCat3;
356     else return kCiCCat4;
357     }
358 fabstoec 1.8
359     // shoud NEVER happen, but you never know...
360     return kCiCNoCat;
361 fabstoec 1.4 }
362 fabstoec 1.5
363     //--------------------------------------------------------------------------------------------------
364     const DecayParticle *PhotonTools::MatchedCiCConversion(const Photon *p, const DecayParticleCol *conversions,
365     Double_t dPhiMin,
366 fabstoec 1.7 Double_t dEtaMin,
367     Double_t dRMin,
368     bool print) {
369 fabstoec 1.5
370 fabstoec 1.8 // if there are no conversons, return
371     if ( !p || !conversions) return NULL;
372    
373 fabstoec 1.6 const DecayParticle *match = NULL;
374    
375     double minDeta = 999.;
376     double minDphi = 999.;
377 fabstoec 1.7 double minDR = 999.;
378 fabstoec 1.6
379     double phPhi = p->SCluster()->Phi();
380     double phEta = p->SCluster()->Eta();
381 fabstoec 1.7
382     if(print)
383     std::cout<<" --- conv match photon eta = "<<phEta<<" phi = "<<phPhi<<std::endl;
384    
385    
386 fabstoec 1.5 for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
387     const DecayParticle *c = conversions->At(i);
388 fabstoec 1.6
389 fabstoec 1.7 if(print)
390     std::cout<< " c "<<i+1<<" pt = "<<c->Pt()<<std::endl;
391    
392 fabstoec 1.6 if(c->Pt() < 1. ) continue; // is this refittedPirMomentum?
393    
394 bendavid 1.24 if (c->Prob() < 1e-6) continue;
395    
396 fabstoec 1.6 //ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
397     //ThreeVector dirconvsc = p->CaloPos() - c->Position();
398    
399     double convPhi = c->Phi();
400     double convEta = c->Eta();
401 fabstoec 1.7
402 fabstoec 1.6 const ThreeVector wrong(0,0,0);
403     double Zvertex = c->DzCorrected(wrong);
404     // ------------------ FROM GLOBE ----------------------
405     //---Definitions for ECAL
406     const float R_ECAL = 136.5;
407     const float Z_Endcap = 328.0;
408     const float etaBarrelEndcap = 1.479;
409 fabstoec 1.7
410 fabstoec 1.6 //---ETA correction
411     float Theta = 0.0 ;
412     float ZEcal = R_ECAL*sinh(convEta)+Zvertex;
413 fabstoec 1.7
414 fabstoec 1.6 if(ZEcal != 0.0) Theta = TMath::ATan(R_ECAL/ZEcal);
415     if(Theta<0.0) Theta = Theta+M_PI;
416     double fEta = - TMath::Log(TMath::Tan(0.5*Theta));
417 fabstoec 1.5
418 fabstoec 1.6 if( fabs(fEta) > etaBarrelEndcap ) {
419     float Zend = Z_Endcap ;
420     if(convEta<0.0 ) Zend = -Zend ;
421     float Zlen = Zend - Zvertex ;
422     float RR = Zlen/TMath::SinH(convEta);
423     Theta = TMath::ATan(RR/Zend);
424     if(Theta<0.0) Theta = Theta+M_PI;
425     fEta = -TMath::Log(TMath::Tan(0.5*Theta));
426     }
427     // ---------------------------------------------------
428    
429 fabstoec 1.7 if(print) {
430     std::cout<<" eta bare = "<<convEta<<std::endl;
431     std::cout<<" eta new = "<<fEta<<std::endl;
432     std::cout<<" phi = "<<convPhi<<std::endl;
433     }
434    
435     convEta = fEta;
436    
437 fabstoec 1.6 Double_t dphi = TMath::Abs(phPhi - convPhi);
438     if(dphi > M_PI) dphi = 2.*M_PI-dphi;
439     //Double_t deta = c->Eta()-dirconvsc.Eta();
440     Double_t deta = TMath::Abs(convEta-phEta);
441 fabstoec 1.7 Double_t dR = TMath::Sqrt(dphi*dphi+deta*deta);
442     //if(dphi < minDphi && TMath::Abs(deta)<minDeta) {
443     if(dR < minDR) {
444     minDR = dR;
445     minDphi = dphi;
446     minDeta = TMath::Abs(deta);
447     match = c;
448    
449     if(print)
450     std::cout<<" conv "<<i+1<<" matches with dPhi = "<<minDphi<<" dEta = "<<minDeta<<std::endl;
451    
452 fabstoec 1.6 }
453 fabstoec 1.5 }
454    
455 fabstoec 1.7 //if(minDphi < dPhiMin && minDeta < dEtaMin)
456     if(minDR < dRMin)
457     return match;
458     else
459     return NULL;
460 fabstoec 1.5 }
461    
462 fabstoec 1.9 bool PhotonTools::PassCiCSelection(const Photon* ph, const Vertex* vtx,
463     const TrackCol* trackCol,
464     const ElectronCol* eleCol,
465     const VertexCol* vtxCol,
466 fabstoec 1.5 double rho, double ptmin,
467 fabstoec 1.9 bool applyEleVeto,
468 fabstoec 1.5 bool print, float* kin) {
469 fabstoec 1.9
470 fabstoec 1.5
471     // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
472     float cic4_allcuts_temp_sublead[] = {
473     3.8, 2.2, 1.77, 1.29,
474     11.7, 3.4, 3.9, 1.84,
475     3.5, 2.2, 2.3, 1.45,
476     0.0106, 0.0097, 0.028, 0.027,
477     0.082, 0.062, 0.065, 0.048,
478     0.94, 0.36, 0.94, 0.32,
479     1., 0.062, 0.97, 0.97,
480     1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
481    
482     // cut on Et instead of Pt???
483     Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
484    
485     // compute all relevant observables first
486     double ecalIso3 = ph->EcalRecHitIsoDr03();
487     double ecalIso4 = ph->EcalRecHitIsoDr04();
488     double hcalIso4 = ph->HcalTowerSumEtDr04();
489    
490     unsigned int wVtxInd = 0;
491    
492     double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
493    
494     // track iso only
495     double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
496    
497     // track iso worst vtx
498     double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
499    
500     double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
501     double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
502    
503     double tIso1 = (combIso1) *50./ph->Et();
504     double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
505     //double tIso2 = (combIso2) *50./ph->Et();
506     double tIso3 = (trackIso3)*50./ph->Et();
507    
508     double covIEtaIEta =ph->CoviEtaiEta();
509     double HoE = ph->HadOverEm();
510    
511     double R9 = ph->R9();
512    
513     double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);
514    
515     // check which category it is ...
516     int _tCat = 1;
517     if ( !isbarrel ) _tCat = 3;
518     if ( R9 < 0.94 ) _tCat++;
519    
520     if(print) {
521     std::cout<<" -------------------------- "<<std::endl;
522 fabstoec 1.15
523     std::cout<<" trackIso1 = "<<trackIso1<<std::endl;
524     std::cout<<" trackIso2 = "<<trackIso2<<std::endl;
525     std::cout<<" trackIso3 = "<<trackIso3<<std::endl;
526     std::cout<<" ecalIso3 = "<<ecalIso3<<std::endl;
527     std::cout<<" ecalIso4 = "<<ecalIso4<<std::endl;
528     std::cout<<" hcalIso4 = "<<hcalIso4<<std::endl;
529    
530 fabstoec 1.5 std::cout<<" photon Et = "<<ph->Et()<<std::endl;
531 mingyang 1.17 std::cout<<" Eta = "<<ph->SCluster()->Eta()<<std::endl;
532 fabstoec 1.15 std::cout<<" HoE = "<<HoE<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+4*4]<<")"<<std::endl;
533     std::cout<<" R9 = "<<R9<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+5*4]<<")"<<std::endl;
534     std::cout<<" dR = "<<dRTrack<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+6*4]<<")"<<std::endl;
535     std::cout<<" iso1 = "<<tIso1<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+0*4]<<")"<<std::endl;
536     std::cout<<" iso2 = "<<tIso2<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+1*4]<<")"<<std::endl;
537     std::cout<<" iso3 = "<<tIso3<<" ("<<cic4_allcuts_temp_sublead[_tCat-1+2*4]<<")"<<std::endl;
538 fabstoec 1.5 }
539    
540     if(kin) {
541     kin[0] = tIso1;
542     kin[1] = tIso2;
543     kin[2] = tIso3;
544     kin[3] = covIEtaIEta;
545     kin[4] = HoE;
546     kin[5] = R9;
547     kin[6] = dRTrack;
548     kin[7] = (float) ph->Pt();
549     kin[8] = (float) ph->Eta();
550     kin[9] = (float) ph->Phi();
551    
552     // iso quantities separate
553     kin[10] = ecalIso3;
554     kin[11] = ecalIso4;
555     kin[12] = hcalIso4;
556     kin[13] = trackIso1;
557     kin[14] = trackIso2;
558     kin[15] = trackIso3;
559    
560     kin[16] = (float) ph->Et();
561     kin[17] = (float) ph->E();
562    
563 bendavid 1.11 kin[18] = (float) ph->SCluster()->Eta();
564     kin[19] = (float) ph->SCluster()->Phi();
565 fabstoec 1.5 }
566    
567     float passCuts = 1.;
568    
569     if ( ph->Pt() <= ptmin ) passCuts = -1.;
570    
571     // not needed anymore, do in pre-selection...
572     //if ( ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
573    
574     if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
575     && tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
576     && tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
577     && covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
578     && HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
579     && R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4]
580 fabstoec 1.9 && ( dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] || !applyEleVeto ) ) ) passCuts = -1.;
581 fabstoec 1.5
582 fabstoec 1.15 if(print) std::cout<<" ---> "<<passCuts<<std::endl;
583    
584 fabstoec 1.5 if(kin) {
585 bendavid 1.11 kin[20] = passCuts;
586     kin[21] = (float) _tCat;
587 fabstoec 1.5 }
588    
589     if(passCuts > 0.) return true;
590     return false;
591     }
592 bendavid 1.11
593 bendavid 1.13 const MCParticle *PhotonTools::MatchMC(const Particle *ph, const MCParticleCol *c, Bool_t matchElectrons) {
594    
595     // printf("Start loop\n");
596 bendavid 1.11 for (UInt_t i=0; i<c->GetEntries(); ++i) {
597     const MCParticle *p = c->At(i);
598 bendavid 1.13 // if (p->IsGenerated() && p->AbsPdgId()==11 && (p->DistinctMother()->AbsPdgId()==23|| p->DistinctMother()->AbsPdgId()==24)) {
599     // printf("pdgid = %i, status = %i, pt = %5f, mother pdgid = %i\n",p->PdgId(),p->Status(),p->Pt(),p->Mother()->PdgId());
600 bendavid 1.12 // }
601     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)) {
602     return p;
603     }
604 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) ) {
605     return p;
606     }
607     }
608     return 0;
609 bendavid 1.12 }
610    
611     PhotonTools::DiphotonR9EtaPtCats PhotonTools::DiphotonR9EtaPtCat(const Photon *p1, const Photon *p2) {
612    
613     PhotonTools::CiCBaseLineCats cat1 = CiCBaseLineCat(p1);
614     PhotonTools::CiCBaseLineCats cat2 = CiCBaseLineCat(p2);
615    
616     PhotonTools::DiphotonR9EtaPtCats evtcat = PhotonTools::kOctCat0;
617    
618     bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
619     bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
620    
621     bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
622     bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
623    
624     if( ph1IsEB && ph2IsEB )
625     evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat0 : PhotonTools::kOctCat1);
626     else
627     evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat2 : PhotonTools::kOctCat3);
628    
629     float ptgg = (p1->Mom()+p2->Mom()).Pt();
630     if (ptgg<40.0) evtcat = PhotonTools::DiphotonR9EtaPtCats(evtcat + 4);
631    
632    
633     return evtcat;
634 fabstoec 1.15 }
635 mingyang 1.17
636 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) {
637 mingyang 1.17 float ScEta=p->SCluster()->Eta();
638     float Et=p->Et();
639     float R9=p->R9();
640     float HoE=p->HadOverEm();
641     float CovIEtaIEta=p->CoviEtaiEta();
642     float EcalIsoDr03=p->EcalRecHitIsoDr03();
643     float HcalIsoDr03=p->HcalTowerSumEtDr03();
644     float TrkIsoHollowDr03=p->HollowConeTrkIsoDr03();
645     float NewEcalIso=EcalIsoDr03-0.012*Et;
646     float NewHcalIso=HcalIsoDr03-0.005*Et;
647     float NewTrkIsoHollowDr03=TrkIsoHollowDr03-0.002*Et;
648     Bool_t IsBarrel=kFALSE;
649     Bool_t IsEndcap=kFALSE;
650 bendavid 1.22 Bool_t PassEleVetoRaw = PhotonTools::PassElectronVetoConvRecovery(p, els, conversions, bs);
651 bendavid 1.20 Bool_t PassEleVeto = (!applyElectronVeto && !invertElectronVeto) || (applyElectronVeto && !invertElectronVeto && PassEleVetoRaw) || (!applyElectronVeto && invertElectronVeto && !PassEleVetoRaw);
652 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) );
653 mingyang 1.18 float HcalEcalPUCorr=EcalIsoDr03+HcalIsoDr03-0.17*rho;
654    
655 mingyang 1.17 if(fabs(ScEta)<1.4442){IsBarrel=kTRUE;}
656     if(fabs(ScEta)>1.566 && fabs(ScEta)<2.5){IsEndcap=kTRUE;}
657     if((!IsBarrel) && (!IsEndcap)){
658     return kFALSE;
659     }
660     if(!PassEleVeto){
661     return kFALSE;
662     }
663     if(R9<=0.9){
664 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) {
665 mingyang 1.17 return kTRUE;
666     }
667     }
668     if(R9>0.9){
669 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) {
670 mingyang 1.17 return kTRUE;
671     }
672     }
673     return kFALSE;
674     }