ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/PhotonTools.cc
Revision: 1.6
Committed: Wed Jun 29 18:28:05 2011 UTC (13 years, 10 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
Changes since 1.5: +56 -13 lines
Log Message:
updated CiC module

File Contents

# User Rev Content
1 fabstoec 1.6 // $Id: PhotonTools.cc,v 1.5 2011/06/27 12:32:21 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    
9     ClassImp(mithep::PhotonTools)
10    
11     using namespace mithep;
12    
13     //--------------------------------------------------------------------------------------------------
14     PhotonTools::PhotonTools()
15     {
16     // Constructor.
17     }
18    
19     //--------------------------------------------------------------------------------------------------
20     Bool_t PhotonTools::PassConversionId(const Photon *p, const DecayParticle *c) {
21    
22     if (!c) return kTRUE;
23    
24     ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
25     Double_t deta = c->Eta()-dirconvsc.Eta();
26     Double_t dphi = MathUtils::DeltaPhi(c->Phi(),dirconvsc.Phi());
27     Double_t eoverp = p->SCluster()->Energy()/c->P();
28    
29     if (p->IsEB() && eoverp>2.0) return kFALSE;
30     if (p->IsEE() && eoverp>3.0) return kFALSE;
31    
32     if (p->IsEE() && TMath::Abs(deta)>0.01) return kFALSE;
33     if (p->IsEE() && TMath::Abs(dphi)>0.01) return kFALSE;
34    
35     return kTRUE;
36    
37     }
38    
39     //--------------------------------------------------------------------------------------------------
40     Bool_t PhotonTools::PassElectronVeto(const Photon *p, const ElectronCol *els) {
41    
42     Bool_t pass = kTRUE;
43     for (UInt_t i=0; i<els->GetEntries(); ++i) {
44     const Electron *e = els->At(i);
45     if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0) {
46     pass = kFALSE;
47     }
48     }
49    
50     return pass;
51     }
52    
53     //--------------------------------------------------------------------------------------------------
54 fabstoec 1.4 Double_t PhotonTools::ElectronVetoCiC(const Photon *p, const ElectronCol *els) {
55    
56     for (UInt_t i=0; i<els->GetEntries(); ++i) {
57     const Electron *e = els->At(i);
58 fabstoec 1.5 if (e->SCluster()==p->SCluster() ) {
59     //if( e->GsfTrk()->NExpectedHitsInner()==0 && e->GsfTrk()->Pt() > 2.5 ) {
60     if( e->GsfTrk()->NExpectedHitsInner()==0 ) {
61     double dEta = e->DeltaEtaSuperClusterTrackAtVtx();
62     double dPhi = e->DeltaPhiSuperClusterTrackAtVtx();
63     double dR = TMath::Sqrt(dEta*dEta+dPhi*dPhi);
64     return dR;
65     }
66     }
67     }
68     return 99.;
69 fabstoec 1.4 }
70    
71     //--------------------------------------------------------------------------------------------------
72 bendavid 1.1 Bool_t PhotonTools::PassElectronVetoConvRecovery(const Photon *p, const ElectronCol *els, const DecayParticleCol *conversions, const BaseVertex *v) {
73    
74     Bool_t pass = kTRUE;
75     for (UInt_t i=0; i<els->GetEntries(); ++i) {
76     const Electron *e = els->At(i);
77     if (e->SCluster()==p->SCluster() && e->GsfTrk()->NExpectedHitsInner()==0 && ElectronTools::PassConversionFilter(e, conversions,
78     v, 0, 1e-6, 2.0, kFALSE, kFALSE) ) {
79     pass = kFALSE;
80     }
81     }
82    
83     return pass;
84     }
85    
86     //--------------------------------------------------------------------------------------------------
87     Bool_t PhotonTools::PassTriggerMatching(const Photon *p, const TriggerObjectCol *trigobjs)
88     {
89    
90     for (UInt_t i=0; i<trigobjs->GetEntries(); ++i) {
91     const TriggerObject *trigobj = trigobjs->At(i);
92     if (trigobj->TriggerType()==TriggerObject::TriggerCluster || trigobj->TriggerType()==TriggerObject::TriggerElectron || trigobj->TriggerType()==TriggerObject::TriggerPhoton) {
93     if (MathUtils::DeltaR(p->SCluster(),trigobj)<0.3) {
94     return kTRUE;
95     }
96     }
97     }
98    
99     return kFALSE;
100    
101    
102     }
103    
104     //--------------------------------------------------------------------------------------------------
105     const DecayParticle *PhotonTools::MatchedConversion(const Photon *p, const DecayParticleCol *conversions,
106     const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
107     Double_t lxyMin, Double_t dRMin) {
108    
109     const DecayParticle *match = 0;
110 bendavid 1.3 Double_t rhosmallest = 999.;
111 bendavid 1.1 for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
112     const DecayParticle *c = conversions->At(i);
113     ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
114     Double_t dr = MathUtils::DeltaR(*c,dirconvsc);
115 bendavid 1.3 Double_t rho = c->Position().Rho();
116     if (dr<dRMin && rho<rhosmallest && c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
117 bendavid 1.1 Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
118     Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
119     if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
120 bendavid 1.3 rhosmallest = rho;
121 bendavid 1.1 match = c;
122     }
123     }
124    
125     }
126    
127     return match;
128    
129     }
130    
131 bendavid 1.2 //--------------------------------------------------------------------------------------------------
132     const DecayParticle *PhotonTools::MatchedConversion(const Track *t, const DecayParticleCol *conversions,
133     const BaseVertex *vtx, Int_t nWrongHitsMax, Double_t probMin,
134     Double_t lxyMin) {
135    
136     for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
137     const DecayParticle *c = conversions->At(i);
138     if (c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
139     Int_t nhb1 = dynamic_cast<const StableData*>(c->DaughterDat(0))->NHitsBeforeVtx();
140     Int_t nhb2 = dynamic_cast<const StableData*>(c->DaughterDat(1))->NHitsBeforeVtx();
141     if (TMath::Max(nhb1,nhb2)<=nWrongHitsMax) {
142     const Track *ct1 = dynamic_cast<const ChargedParticle*>(c->Daughter(0))->Trk();
143     const Track *ct2 = dynamic_cast<const ChargedParticle*>(c->Daughter(1))->Trk();
144     if (t==ct1 || t==ct2) return c;
145     }
146     }
147    
148     }
149    
150     return 0;
151    
152     }
153    
154 bendavid 1.1 PhotonTools::DiphotonR9EtaCats PhotonTools::DiphotonR9EtaCat(const Photon *p1, const Photon *p2) {
155    
156     if (p1->IsEB() && p2->IsEB()) {
157     if (p1->R9()>0.93 && p2->R9()>0.93) return kCat1;
158     else return kCat2;
159    
160     }
161     else {
162     if (p1->R9()>0.93 && p2->R9()>0.93) return kCat3;
163     else return kCat4;
164     }
165    
166     }
167    
168     PhotonTools::DiphotonR9EtaConversionCats PhotonTools::DiphotonR9EtaConversionCat(const Photon *p1, const Photon *p2, const DecayParticleCol *conversions, const BaseVertex *v) {
169    
170     const DecayParticle *conv1 = MatchedConversion(p1, conversions, v);
171     const DecayParticle *conv2 = MatchedConversion(p2, conversions, v);
172    
173     if (p1->IsEB() && p2->IsEB()) {
174     if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat1;
175     else if (conv1||conv2) return kNewCat2;
176     else return kNewCat3;
177    
178     }
179     else {
180     if (p1->R9()>0.93 && p2->R9()>0.93) return kNewCat4;
181     else if (conv1||conv2) return kNewCat5;
182     else return kNewCat6;
183     }
184    
185 fabstoec 1.4 }
186    
187     PhotonTools::CiCBaseLineCats PhotonTools::CiCBaseLineCat(const Photon *p) {
188     if( p->IsEB() ) {
189     if ( p->R9() > 0.94 ) return kCiCCat1;
190     else return kCiCCat2;
191     } else {
192     if ( p->R9() > 0.94 ) return kCiCCat3;
193     else return kCiCCat4;
194     }
195     }
196 fabstoec 1.5
197     //--------------------------------------------------------------------------------------------------
198     const DecayParticle *PhotonTools::MatchedCiCConversion(const Photon *p, const DecayParticleCol *conversions,
199     Double_t dPhiMin,
200     Double_t dEtaMin) {
201    
202    
203 fabstoec 1.6 const DecayParticle *match = NULL;
204    
205     double minDeta = 999.;
206     double minDphi = 999.;
207    
208     double phPhi = p->SCluster()->Phi();
209     double phEta = p->SCluster()->Eta();
210 fabstoec 1.5
211     for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
212     const DecayParticle *c = conversions->At(i);
213 fabstoec 1.6
214     if(c->Pt() < 1. ) continue; // is this refittedPirMomentum?
215    
216     //ThreeVector dirconvsc = ThreeVector(p->SCluster()->Point()) - c->Position();
217     //ThreeVector dirconvsc = p->CaloPos() - c->Position();
218    
219     double convPhi = c->Phi();
220     double convEta = c->Eta();
221    
222     const ThreeVector wrong(0,0,0);
223     double Zvertex = c->DzCorrected(wrong);
224     // ------------------ FROM GLOBE ----------------------
225     //---Definitions for ECAL
226     const float R_ECAL = 136.5;
227     const float Z_Endcap = 328.0;
228     const float etaBarrelEndcap = 1.479;
229    
230     //---ETA correction
231     float Theta = 0.0 ;
232     float ZEcal = R_ECAL*sinh(convEta)+Zvertex;
233    
234     if(ZEcal != 0.0) Theta = TMath::ATan(R_ECAL/ZEcal);
235     if(Theta<0.0) Theta = Theta+M_PI;
236     double fEta = - TMath::Log(TMath::Tan(0.5*Theta));
237 fabstoec 1.5
238 fabstoec 1.6 if( fabs(fEta) > etaBarrelEndcap ) {
239     float Zend = Z_Endcap ;
240     if(convEta<0.0 ) Zend = -Zend ;
241     float Zlen = Zend - Zvertex ;
242     float RR = Zlen/TMath::SinH(convEta);
243     Theta = TMath::ATan(RR/Zend);
244     if(Theta<0.0) Theta = Theta+M_PI;
245     fEta = -TMath::Log(TMath::Tan(0.5*Theta));
246     }
247     convEta = fEta;
248     // ---------------------------------------------------
249    
250     Double_t dphi = TMath::Abs(phPhi - convPhi);
251     if(dphi > M_PI) dphi = 2.*M_PI-dphi;
252     //Double_t deta = c->Eta()-dirconvsc.Eta();
253     Double_t deta = TMath::Abs(convEta-phEta);
254     if (dphi<dPhiMin && TMath::Abs(deta)<dEtaMin) {
255     if(dphi < minDphi && TMath::Abs(deta)<minDeta) {
256     minDphi = dphi;
257     minDeta = TMath::Abs(deta);
258     match = c;
259     }
260     }
261 fabstoec 1.5 }
262    
263 fabstoec 1.6 return match;
264 fabstoec 1.5 }
265    
266     bool PhotonTools::PassCiCSelection(Photon* ph, const Vertex* vtx,
267     const TrackCol* trackCol,
268     const ElectronCol* eleCol,
269     const VertexCol* vtxCol,
270     double rho, double ptmin,
271     bool print, float* kin) {
272    
273     // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
274     float cic4_allcuts_temp_sublead[] = {
275     3.8, 2.2, 1.77, 1.29,
276     11.7, 3.4, 3.9, 1.84,
277     3.5, 2.2, 2.3, 1.45,
278     0.0106, 0.0097, 0.028, 0.027,
279     0.082, 0.062, 0.065, 0.048,
280     0.94, 0.36, 0.94, 0.32,
281     1., 0.062, 0.97, 0.97,
282     1.5, 1.5, 1.5, 1.5 }; // the last line is PixelmatchVeto and un-used
283    
284     // cut on Et instead of Pt???
285     Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
286    
287     // compute all relevant observables first
288     double ecalIso3 = ph->EcalRecHitIsoDr03();
289     double ecalIso4 = ph->EcalRecHitIsoDr04();
290     double hcalIso4 = ph->HcalTowerSumEtDr04();
291    
292     unsigned int wVtxInd = 0;
293    
294     double trackIso1 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
295    
296     // track iso only
297     double trackIso3 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol);
298    
299     // track iso worst vtx
300     double trackIso2 = IsolationTools::CiCTrackIsolation(ph, vtx, 0.4, 0.02, 0.0, 0.0, 0.1, 1.0, trackCol, &wVtxInd, vtxCol);
301    
302     double combIso1 = ecalIso3+hcalIso4+trackIso1 - 0.17*rho;
303     double combIso2 = ecalIso4+hcalIso4+trackIso2 - 0.52*rho;
304    
305     double tIso1 = (combIso1) *50./ph->Et();
306     double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
307     //double tIso2 = (combIso2) *50./ph->Et();
308     double tIso3 = (trackIso3)*50./ph->Et();
309    
310     double covIEtaIEta =ph->CoviEtaiEta();
311     double HoE = ph->HadOverEm();
312    
313     double R9 = ph->R9();
314    
315     double dRTrack = PhotonTools::ElectronVetoCiC(ph, eleCol);
316    
317     // check which category it is ...
318     int _tCat = 1;
319     if ( !isbarrel ) _tCat = 3;
320     if ( R9 < 0.94 ) _tCat++;
321    
322     if(print) {
323     std::cout<<" -------------------------- "<<std::endl;
324     std::cout<<" photon Et = "<<ph->Et()<<std::endl;
325     std::cout<<" Eta = "<<ph->SCluster()->Eta()<<std::endl;
326     std::cout<<" HoE = "<<HoE<<std::endl;
327     std::cout<<" R9 = "<<R9<<std::endl;
328     std::cout<<" dR = "<<dRTrack<<std::endl;
329     std::cout<<" iso1 = "<<tIso1<<std::endl;
330     std::cout<<" iso2 = "<<tIso2<<std::endl;
331     std::cout<<" iso3 = "<<tIso3<<std::endl;
332     }
333    
334     if(kin) {
335     kin[0] = tIso1;
336     kin[1] = tIso2;
337     kin[2] = tIso3;
338     kin[3] = covIEtaIEta;
339     kin[4] = HoE;
340     kin[5] = R9;
341     kin[6] = dRTrack;
342     kin[7] = (float) ph->Pt();
343     kin[8] = (float) ph->Eta();
344     kin[9] = (float) ph->Phi();
345    
346     // iso quantities separate
347     kin[10] = ecalIso3;
348     kin[11] = ecalIso4;
349     kin[12] = hcalIso4;
350     kin[13] = trackIso1;
351     kin[14] = trackIso2;
352     kin[15] = trackIso3;
353    
354     kin[16] = (float) ph->Et();
355     kin[17] = (float) ph->E();
356    
357     }
358    
359     float passCuts = 1.;
360    
361     if ( ph->Pt() <= ptmin ) passCuts = -1.;
362    
363     // not needed anymore, do in pre-selection...
364     //if ( ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
365    
366     if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
367     && tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
368     && tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
369     && covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
370     && HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
371     && R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4]
372     && dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] )
373     ) {
374     passCuts = -1.;
375     }
376    
377     if(kin) {
378     kin[18] = passCuts;
379     kin[19] = (float) _tCat;
380     }
381    
382     if(passCuts > 0.) return true;
383     return false;
384     }