ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/PhotonTools.cc
Revision: 1.29
Committed: Mon May 28 14:24:49 2012 UTC (12 years, 11 months ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_028a
Changes since 1.28: +2 -2 lines
Log Message:
Add single leg conversion option for diphoton vertex selection

File Contents

# Content
1 // $Id: PhotonTools.cc,v 1.28 2012/05/27 16:08:50 bendavid Exp $
2
3 #include "MitPhysics/Utils/interface/PhotonTools.h"
4 #include "MitPhysics/Utils/interface/ElectronTools.h"
5 #include "MitPhysics/Utils/interface/IsolationTools.h"
6 #include "MitAna/DataTree/interface/StableData.h"
7 #include <TFile.h>
8 #include <TRandom3.h>
9
10 ClassImp(mithep::PhotonTools)
11
12 using namespace mithep;
13
14 //--------------------------------------------------------------------------------------------------
15 PhotonTools::PhotonTools()
16 {
17 // Constructor.
18 }
19
20 //--------------------------------------------------------------------------------------------------
21 PhotonTools::eScaleCats PhotonTools::EScaleCat(const Photon *p)
22 {
23 if (p->SCluster()->AbsEta()<1.0) {
24 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 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
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 rng->SetSeed(iSeed);
63 FourVectorM mom = p->Mom();
64 Double_t scale = rng->Gaus(1.,width);
65
66 if( scale > 0)
67 p->SetMom(scale*mom.X(), scale*mom.Y(), scale*mom.Z(), scale*mom.E());
68
69 return;
70 }
71
72 void PhotonTools::SmearPhotonError(Photon* p, Double_t width) {
73
74 if( !p ) return;
75 if( width < 0.) return;
76
77
78 Double_t smear = p->EnergySmearing();
79 p->SetEnergySmearing(TMath::Sqrt(smear*smear+width*width*p->E()*p->E()));
80
81 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 }
87
88 void PhotonTools::ScalePhotonR9(Photon* p, Double_t scale) {
89 p->SetR9(scale*p->R9());
90 }
91
92
93 void PhotonTools::ScalePhotonError(Photon* p, Double_t scale) {
94 p->SetEnergyErrSmeared(scale*p->EnergyErrSmeared());
95 }
96
97 //--------------------------------------------------------------------------------------------------
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 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 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 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 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 }
224
225 //--------------------------------------------------------------------------------------------------
226 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 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 const DecayParticle *match = 0;
273 Double_t rhosmallest = 999.;
274 for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
275 const DecayParticle *c = conversions->At(i);
276 ThreeVector dirconvsc = ThreeVector(sc->Point()) - c->Position();
277 Double_t dr = MathUtils::DeltaR(*c,dirconvsc);
278 Double_t rho = c->Position().Rho();
279 if (dr<dRMin && rho<rhosmallest && c->Prob()>probMin && c->LxyCorrected(vtx)>lxyMin) {
280 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 rhosmallest = rho;
284 match = c;
285 }
286 }
287
288 }
289
290 return match;
291
292 }
293
294 //--------------------------------------------------------------------------------------------------
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 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 }
349
350 PhotonTools::CiCBaseLineCats PhotonTools::CiCBaseLineCat(const Photon *p) {
351 if( p->SCluster()->AbsEta()<1.5 ) {
352 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
359 // shoud NEVER happen, but you never know...
360 return kCiCNoCat;
361 }
362
363 //--------------------------------------------------------------------------------------------------
364 const DecayParticle *PhotonTools::MatchedCiCConversion(const Photon *p, const DecayParticleCol *conversions,
365 Double_t dPhiMin,
366 Double_t dEtaMin,
367 Double_t dRMin,
368 bool print) {
369
370 // if there are no conversons, return
371 if ( !p || !conversions) return NULL;
372
373 const DecayParticle *match = NULL;
374
375 double minDeta = 999.;
376 double minDphi = 999.;
377 double minDR = 999.;
378
379 double phPhi = p->SCluster()->Phi();
380 double phEta = p->SCluster()->Eta();
381
382 if(print)
383 std::cout<<" --- conv match photon eta = "<<phEta<<" phi = "<<phPhi<<std::endl;
384
385
386 for (UInt_t i=0; i<conversions->GetEntries(); ++i) {
387 const DecayParticle *c = conversions->At(i);
388
389 if(print)
390 std::cout<< " c "<<i+1<<" pt = "<<c->Pt()<<std::endl;
391
392 if(c->Pt() < 1. ) continue; // is this refittedPirMomentum?
393
394 if (c->NDaughters()==2 && c->Prob() < 1e-6) continue;
395
396 //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
402 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
410 //---ETA correction
411 float Theta = 0.0 ;
412 float ZEcal = R_ECAL*sinh(convEta)+Zvertex;
413
414 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
418 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 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 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 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 }
453 }
454
455 //if(minDphi < dPhiMin && minDeta < dEtaMin)
456 if(minDR < dRMin)
457 return match;
458 else
459 return NULL;
460 }
461
462 bool PhotonTools::PassCiCSelection(const Photon* ph, const Vertex* vtx,
463 const TrackCol* trackCol,
464 const ElectronCol* eleCol,
465 const VertexCol* vtxCol,
466 double rho, double ptmin,
467 bool applyEleVeto,
468 bool print, float* kin) {
469
470
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
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 std::cout<<" photon Et = "<<ph->Et()<<std::endl;
531 std::cout<<" Eta = "<<ph->SCluster()->Eta()<<std::endl;
532 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 }
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 kin[18] = (float) ph->SCluster()->Eta();
564 kin[19] = (float) ph->SCluster()->Phi();
565 }
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 && ( dRTrack > cic4_allcuts_temp_sublead[_tCat-1+6*4] || !applyEleVeto ) ) ) passCuts = -1.;
581
582 if(print) std::cout<<" ---> "<<passCuts<<std::endl;
583
584 if(kin) {
585 kin[20] = passCuts;
586 kin[21] = (float) _tCat;
587 }
588
589 if(passCuts > 0.) return true;
590 return false;
591 }
592
593
594
595 bool PhotonTools::PassCiCPFIsoSelection(const Photon* ph,
596 const Vertex* vtx,
597 const PFCandidateCol* pfCol,
598 const VertexCol* vtxCol,
599 double rho, double ptmin)
600 {
601
602
603 // these values are taken from the H2GGlobe code... (actually from Marco/s mail)
604 float cic4_allcuts_temp_sublead[] = {
605 6.0, 4.7, 5.6, 3.6,
606 10.0, 6.5, 5.6, 4.4,
607 3.8, 2.5, 3.1, 2.2,
608 0.0108, 0.0102, 0.028, 0.028,
609 0.124, 0.092, 0.142, 0.063,
610 0.94, 0.28, 0.94, 0.24 }; // the last line is PixelmatchVeto and un-used
611
612 // cut on Et instead of Pt???
613 Bool_t isbarrel = ph->SCluster()->AbsEta()<1.5;
614
615
616
617 // compute all relevant observables first
618 double ecalIso3 = IsolationTools::PFGammaIsolation(ph, 0.3, 0.0, pfCol);
619 double ecalIso4 = IsolationTools::PFGammaIsolation(ph, 0.4, 0.0, pfCol);
620
621 unsigned int wVtxInd = 0;
622
623 double trackIsoSel03 = IsolationTools::PFChargedIsolation(ph, vtx, 0.3, 0.0, pfCol);
624
625 // track iso worst vtx
626 double trackIsoWorst04 = IsolationTools::PFChargedIsolation(ph, vtx, 0.4, 0.00, pfCol, &wVtxInd, vtxCol);
627
628 double combIso1 = ecalIso3+trackIsoSel03 - 0.09*rho;
629 double combIso2 = ecalIso4+trackIsoWorst04 - 0.23*rho;
630
631 double tIso1 = (combIso1) *50./ph->Et();
632 double tIso2 = (combIso2) *50./(ph->MomVtx(vtxCol->At(wVtxInd)->Position()).Pt());
633 //double tIso2 = (combIso2) *50./ph->Et();
634 double tIso3 = (trackIsoSel03)*50./ph->Et();
635
636 double covIEtaIEta =ph->CoviEtaiEta();
637 double HoE = ph->HadOverEm();
638
639 double R9 = ph->R9();
640
641 // check which category it is ...
642 int _tCat = 1;
643 if ( !isbarrel ) _tCat = 3;
644 if ( R9 < 0.94 ) _tCat++;
645
646 float passCuts = 1.;
647
648 if ( ph->Pt() <= ptmin ) passCuts = -1.;
649
650 // not needed anymore, do in pre-selection...
651 //if ( ph->SCluster()->AbsEta()>=2.5 || (ph->SCluster()->AbsEta()>=1.4442 && ph->SCluster()->AbsEta()<=1.566)) passCuts = -1.;
652
653 if( ! ( tIso1 < cic4_allcuts_temp_sublead[_tCat-1+0*4]
654 && tIso2 < cic4_allcuts_temp_sublead[_tCat-1+1*4]
655 && tIso3 < cic4_allcuts_temp_sublead[_tCat-1+2*4]
656 && covIEtaIEta < cic4_allcuts_temp_sublead[_tCat-1+3*4]
657 && HoE < cic4_allcuts_temp_sublead[_tCat-1+4*4]
658 && R9 > cic4_allcuts_temp_sublead[_tCat-1+5*4] ) ) passCuts = -1.;
659
660
661 if(passCuts > 0.) return true;
662 return false;
663 }
664
665
666
667
668
669 const MCParticle *PhotonTools::MatchMC(const Particle *ph, const MCParticleCol *c, Bool_t matchElectrons) {
670
671 // printf("Start loop\n");
672 for (UInt_t i=0; i<c->GetEntries(); ++i) {
673 const MCParticle *p = c->At(i);
674 // if (p->IsGenerated() && p->AbsPdgId()==11 && (p->DistinctMother()->AbsPdgId()==23|| p->DistinctMother()->AbsPdgId()==24)) {
675 // printf("pdgid = %i, status = %i, pt = %5f, mother pdgid = %i\n",p->PdgId(),p->Status(),p->Pt(),p->Mother()->PdgId());
676 // }
677 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)) {
678 return p;
679 }
680 if ( p->AbsPdgId()==22 && p->IsGenerated() && MathUtils::DeltaR(*ph,*p) < 0.3 && p->Mother() && (p->Mother()->AbsPdgId()==25 || p->Mother()->AbsPdgId()<=21) ) {
681 return p;
682 }
683 }
684 return 0;
685 }
686
687 PhotonTools::DiphotonR9EtaPtCats PhotonTools::DiphotonR9EtaPtCat(const Photon *p1, const Photon *p2) {
688
689 PhotonTools::CiCBaseLineCats cat1 = CiCBaseLineCat(p1);
690 PhotonTools::CiCBaseLineCats cat2 = CiCBaseLineCat(p2);
691
692 PhotonTools::DiphotonR9EtaPtCats evtcat = PhotonTools::kOctCat0;
693
694 bool ph1IsEB = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat2);
695 bool ph2IsEB = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat2);
696
697 bool ph1IsHR9 = (cat1 == PhotonTools::kCiCCat1 || cat1 == PhotonTools::kCiCCat3);
698 bool ph2IsHR9 = (cat2 == PhotonTools::kCiCCat1 || cat2 == PhotonTools::kCiCCat3);
699
700 if( ph1IsEB && ph2IsEB )
701 evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat0 : PhotonTools::kOctCat1);
702 else
703 evtcat = ( ph1IsHR9 && ph2IsHR9 ? PhotonTools::kOctCat2 : PhotonTools::kOctCat3);
704
705 float ptgg = (p1->Mom()+p2->Mom()).Pt();
706 if (ptgg<40.0) evtcat = PhotonTools::DiphotonR9EtaPtCats(evtcat + 4);
707
708
709 return evtcat;
710 }
711
712 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) {
713 float ScEta=p->SCluster()->Eta();
714 float Et=p->Et();
715 float R9=p->R9();
716 float HoE=p->HadOverEm();
717 float CovIEtaIEta=p->CoviEtaiEta();
718 float EcalIsoDr03=p->EcalRecHitIsoDr03();
719 float HcalIsoDr03=p->HcalTowerSumEtDr03();
720 float TrkIsoHollowDr03=p->HollowConeTrkIsoDr03();
721 float NewEcalIso=EcalIsoDr03-0.012*Et;
722 float NewHcalIso=HcalIsoDr03-0.005*Et;
723 float NewTrkIsoHollowDr03=TrkIsoHollowDr03-0.002*Et;
724 Bool_t IsBarrel=kFALSE;
725 Bool_t IsEndcap=kFALSE;
726 Bool_t PassEleVetoRaw = PhotonTools::PassElectronVetoConvRecovery(p, els, conversions, bs);
727 Bool_t PassEleVeto = (!applyElectronVeto && !invertElectronVeto) || (applyElectronVeto && !invertElectronVeto && PassEleVetoRaw) || (!applyElectronVeto && invertElectronVeto && !PassEleVetoRaw);
728 float AbsTrackIsoCIC=IsolationTools::CiCTrackIsolation(p,vtx, 0.3, 0.02, 0.0, 0.0, 0.1, 1.0,trackCol, NULL, NULL, (!applyElectronVeto ? els : NULL) );
729 float HcalEcalPUCorr=EcalIsoDr03+HcalIsoDr03-0.17*rho;
730
731 if(fabs(ScEta)<1.4442){IsBarrel=kTRUE;}
732 if(fabs(ScEta)>1.566 && fabs(ScEta)<2.5){IsEndcap=kTRUE;}
733 if((!IsBarrel) && (!IsEndcap)){
734 return kFALSE;
735 }
736 if(!PassEleVeto){
737 return kFALSE;
738 }
739 if(R9<=0.9){
740 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) {
741 return kTRUE;
742 }
743 }
744 if(R9>0.9){
745 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) {
746 return kTRUE;
747 }
748 }
749 return kFALSE;
750 }
751
752 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) {
753 float ScEta=p->SCluster()->Eta();
754 float Et=p->Et();
755 float R9=p->R9();
756 float HoE=p->HadOverEm();
757 float CovIEtaIEta=p->CoviEtaiEta();
758 float EcalIsoDr03=p->EcalRecHitIsoDr03();
759 float HcalIsoDr03=p->HcalTowerSumEtDr03();
760 float TrkIsoHollowDr03=p->HollowConeTrkIsoDr03();
761 float NewEcalIso=EcalIsoDr03-0.012*Et;
762 float NewHcalIso=HcalIsoDr03-0.005*Et;
763 float NewTrkIsoHollowDr03=TrkIsoHollowDr03-0.002*Et;
764 Bool_t IsBarrel=kFALSE;
765 Bool_t IsEndcap=kFALSE;
766 Bool_t PassEleVetoRaw = PhotonTools::PassElectronVetoConvRecovery(p, els, conversions, bs);
767 Bool_t PassEleVeto = (!applyElectronVeto && !invertElectronVeto) || (applyElectronVeto && !invertElectronVeto && PassEleVetoRaw) || (!applyElectronVeto && invertElectronVeto && !PassEleVetoRaw);
768
769 float ChargedIso_selvtx_DR002To0p02=IsolationTools::PFChargedIsolation(p,vtx, 0.2, 0.,fPFCands);
770
771 if(fabs(ScEta)<1.4442){IsBarrel=kTRUE;}
772 if(fabs(ScEta)>1.566 && fabs(ScEta)<2.5){IsEndcap=kTRUE;}
773 if((!IsBarrel) && (!IsEndcap)){
774 return kFALSE;
775 }
776 if(!PassEleVeto){
777 return kFALSE;
778 }
779 if(R9<=0.9){
780 if(HoE<0.075 && ((IsBarrel && CovIEtaIEta<0.014) || (IsEndcap && CovIEtaIEta<0.034)) && NewEcalIso<4 && NewHcalIso<4 && NewTrkIsoHollowDr03<4 && ChargedIso_selvtx_DR002To0p02<4) {
781 return kTRUE;
782 }
783 }
784 if(R9>0.9){
785 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) {
786 return kTRUE;
787 }
788 }
789 return kFALSE;
790 }