ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/IsolationTools.cc
Revision: 1.21
Committed: Fri Sep 16 14:09:34 2011 UTC (13 years, 7 months ago) by ceballos
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025, Mit_025pre2
Changes since 1.20: +7 -4 lines
Log Message: 
fixing iso

File Contents

# User Rev Content
1 ceballos 1.21 // $Id: IsolationTools.cc,v 1.20 2011/06/27 12:32:21 fabstoec Exp $
2 loizides 1.1
3     #include "MitPhysics/Utils/interface/IsolationTools.h"
4 bendavid 1.11 #include "MitPhysics/Utils/interface/PhotonTools.h"
5 loizides 1.1 #include "MitCommon/MathTools/interface/MathUtils.h"
6    
7 loizides 1.4 ClassImp(mithep::IsolationTools)
8    
9 loizides 1.1 using namespace mithep;
10    
11     //--------------------------------------------------------------------------------------------------
12     Double_t IsolationTools::TrackIsolation(const Track *p, Double_t extRadius, Double_t intRadius,
13     Double_t ptLow, Double_t maxVtxZDist,
14 phedex 1.2 const Collection<Track> *tracks)
15 loizides 1.1 {
16     //Computes the Track Isolation: Summed Transverse Momentum of all tracks inside an
17     //annulus around the electron seed track.
18    
19     Double_t ptSum =0.;
20     for (UInt_t i=0; i<tracks->GetEntries();i++) {
21     Double_t tmpPt = tracks->At(i)->Pt();
22     Double_t deltaZ = fabs(p->Z0() - tracks->At(i)->Z0());
23    
24     //ignore the track if it is below the pt threshold
25     if (tmpPt < ptLow)
26     continue;
27     //ingore the track if it is too far away in Z
28     if (deltaZ > maxVtxZDist)
29     continue;
30    
31 phedex 1.2 Double_t dr = MathUtils::DeltaR(p->Phi(),p->Eta(),tracks->At(i)->Phi(), tracks->At(i)->Eta());
32 loizides 1.1 //add the track pt if it is inside the annulus
33     if ( dr < extRadius &&
34     dr >= intRadius ) {
35     ptSum += tmpPt;
36     }
37     }
38     return ptSum;
39     }
40    
41     //--------------------------------------------------------------------------------------------------
42     Double_t IsolationTools::EcalIsolation(const SuperCluster *sc, Double_t coneSize, Double_t etLow,
43 phedex 1.2 const Collection<BasicCluster> *basicClusters)
44 loizides 1.1 {
45     //Computes the Ecal Isolation: Summed Transverse Energy of all Basic Clusters inside a
46     //cone around the electron, excluding those that are inside the electron super cluster.
47    
48     Double_t ecalIsol=0.;
49     const BasicCluster *basicCluster= 0;
50     for (UInt_t i=0; i<basicClusters->GetEntries();i++) {
51     basicCluster = basicClusters->At(i);
52     Double_t basicClusterEnergy = basicCluster->Energy();
53     Double_t basicClusterEta = basicCluster->Eta();
54     Double_t basicClusterEt = basicClusterEnergy*sin(2*atan(exp(basicClusterEta)));
55    
56 bendavid 1.3 if (basicClusterEt > etLow) {
57 loizides 1.1 bool inSuperCluster = false;
58    
59     // loop over the basic clusters of the supercluster
60     // to make sure that the basic cluster is not inside
61     // the super cluster. We exclude those.
62     for (UInt_t j=0; j<sc->ClusterSize(); j++) {
63     const BasicCluster *tempBasicClusterInSuperCluster = sc->Cluster(j);
64     if (tempBasicClusterInSuperCluster == basicCluster) {
65     inSuperCluster = true;
66     }
67     }
68    
69     if (!inSuperCluster) {
70 phedex 1.2 Double_t dr = MathUtils::DeltaR(sc->Phi(), sc->Eta(),
71     basicCluster->Phi(),basicCluster->Eta());
72 loizides 1.1 if(dr < coneSize) {
73     ecalIsol += basicClusterEt;
74     }
75     }
76     }
77     }
78     return ecalIsol;
79     }
80    
81     //--------------------------------------------------------------------------------------------------
82     Double_t IsolationTools::CaloTowerHadIsolation(const ThreeVector *p, Double_t extRadius,
83     Double_t intRadius, Double_t etLow,
84 phedex 1.2 const Collection<CaloTower> *caloTowers)
85 loizides 1.1 {
86     //Computes the CaloTower Had Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
87     //inside an annulus around the electron super cluster position.
88    
89     Double_t sumEt = 0;
90     for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
91     Double_t caloTowerEt = caloTowers->At(i)->HadEt();
92     Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
93     p->Phi(), p->Eta());
94     if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
95     sumEt += caloTowerEt;
96     }
97     }
98     return sumEt;
99     }
100    
101     //--------------------------------------------------------------------------------------------------
102     Double_t IsolationTools::CaloTowerEmIsolation(const ThreeVector *p, Double_t extRadius,
103     Double_t intRadius, Double_t etLow,
104 phedex 1.2 const Collection<CaloTower> *caloTowers)
105 loizides 1.1 {
106     //Computes the CaloTower Em Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
107     //inside an annulus around the electron super cluster position.
108    
109     Double_t sumEt = 0;
110     for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
111     Double_t caloTowerEt = caloTowers->At(i)->EmEt();
112     Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
113     p->Phi(), p->Eta());
114     if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
115     sumEt += caloTowerEt;
116     }
117     }
118     return sumEt;
119     }
120 ceballos 1.5
121     //--------------------------------------------------------------------------------------------------
122     Double_t IsolationTools::PFMuonIsolation(const Muon *p, const Collection<PFCandidate> *PFCands,
123 ceballos 1.8 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
124 ceballos 1.21 Double_t extRadius, Double_t intRadiusGamma, Double_t intRadius)
125 mzanetti 1.15 {
126     //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
127     //annulus around the particle seed track.
128    
129     Double_t zLepton = 0.0;
130     if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
131    
132     Double_t ptSum =0.;
133     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
134     const PFCandidate *pf = PFCands->At(i);
135    
136     // exclude muon
137     if(pf->TrackerTrk() && p->TrackerTrk() &&
138     pf->TrackerTrk() == p->TrackerTrk()) continue;
139    
140 mzanetti 1.18 Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
141    
142     // pt cut applied to neutrals
143     if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
144    
145 mzanetti 1.15 // ignore the pf candidate if it is too far away in Z
146     Double_t deltaZ = 0.0;
147     if(pf->HasTrk()) {
148     deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
149     }
150 mzanetti 1.18 if (deltaZ >= delta_z) continue;
151    
152     // inner cone veto for gammas
153 ceballos 1.21 if (pf->PFType() == PFCandidate::eGamma && dr < intRadiusGamma) continue;
154 mzanetti 1.18
155 ceballos 1.21 // inner cone veto for tracks
156     if (dr < intRadius) continue;
157    
158 mzanetti 1.18 // add the pf pt if it is inside the extRadius
159 ceballos 1.21 if (dr < extRadius) ptSum += pf->Pt();
160 mzanetti 1.15
161     }
162     return ptSum;
163     }
164    
165     //--------------------------------------------------------------------------------------------------
166     Double_t IsolationTools::PFMuonIsolation(const Muon *p, const Collection<PFCandidate> *PFCands, const Vertex *vertex,
167     const MuonCol *goodMuons, const ElectronCol *goodElectrons,
168     Double_t delta_z, Double_t ptMin, Double_t extRadius,
169     Double_t intRadius, int isoType, Double_t beta)
170 ceballos 1.5 {
171     //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
172     //annulus around the particle seed track.
173    
174     Double_t zLepton = 0.0;
175 ceballos 1.8 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
176 ceballos 1.5
177     Double_t ptSum =0.;
178     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
179     const PFCandidate *pf = PFCands->At(i);
180    
181     Bool_t isGoodType = kFALSE;
182     // all particles
183 ceballos 1.6 if (isoType == 0) isGoodType = kTRUE;
184 ceballos 1.5 // charged particles only
185 ceballos 1.6 else if(isoType == 1 && pf->BestTrk()) isGoodType = kTRUE;
186 ceballos 1.5 // charged particles and gammas only
187 ceballos 1.6 else if(isoType == 2 &&
188     (pf->BestTrk() || pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
189 ceballos 1.8 // all particles, rejecting good leptons
190     else if(isoType == 3) isGoodType = kTRUE;
191 ceballos 1.5
192     if(isGoodType == kFALSE) continue;
193    
194 mzanetti 1.13 // pt cut applied to neutrals
195     if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
196 ceballos 1.6
197 ceballos 1.5 if(pf->TrackerTrk() && p->TrackerTrk() &&
198     pf->TrackerTrk() == p->TrackerTrk()) continue;
199    
200     Double_t deltaZ = 0.0;
201     if(pf->BestTrk()) {
202 ceballos 1.8 deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
203 ceballos 1.5 }
204    
205     // ignore the pf candidate if it is too far away in Z
206 ceballos 1.6 if (deltaZ >= delta_z)
207 ceballos 1.5 continue;
208    
209     Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
210 mzanetti 1.18
211     // inner cone veto for gammas
212     if (pf->PFType() == PFCandidate::eGamma && dr < intRadius) continue;
213    
214 ceballos 1.5 // add the pf pt if it is inside the extRadius and outside the intRadius
215 mzanetti 1.18 if (dr < extRadius ) {
216 ceballos 1.8 Bool_t isLepton = kFALSE;
217     if(goodMuons && isoType == 3){
218     for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
219     const Muon *m = goodMuons->At(nl);
220     if(pf->TrackerTrk() && m->TrackerTrk() &&
221     pf->TrackerTrk() == m->TrackerTrk()) {
222     isLepton = kTRUE;
223     break;
224     }
225     }
226     }
227     if(goodElectrons && isLepton == kFALSE && isoType == 3){
228     for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
229     const Electron *e = goodElectrons->At(nl);
230     if(pf->TrackerTrk() && e->TrackerTrk() &&
231     pf->TrackerTrk() == e->TrackerTrk()) {
232     isLepton = kTRUE;
233     break;
234     }
235     if(pf->GsfTrk() && e->GsfTrk() &&
236     pf->GsfTrk() == e->GsfTrk()) {
237     isLepton = kTRUE;
238     break;
239     }
240     }
241     }
242     if(isLepton == kFALSE){
243     if(pf->BestTrk()) ptSum += pf->Pt();
244     else ptSum += pf->Pt()*beta;
245     }
246 ceballos 1.5 }
247     }
248     return ptSum;
249     }
250 mzanetti 1.15
251 ceballos 1.5 //--------------------------------------------------------------------------------------------------
252     Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands,
253 ceballos 1.8 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
254 sixie 1.16 Double_t extRadius, Double_t intRadius)
255     {
256    
257     //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
258     //annulus around the particle seed track.
259    
260     Double_t zLepton = 0.0;
261     if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
262    
263     Double_t ptSum =0.;
264     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
265     const PFCandidate *pf = PFCands->At(i);
266    
267     // pt cut applied to neutrals
268     if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
269    
270     if(pf->TrackerTrk() && p->TrackerTrk() &&
271     pf->TrackerTrk() == p->TrackerTrk()) continue;
272    
273     if(pf->GsfTrk() && p->GsfTrk() &&
274     pf->GsfTrk() == p->GsfTrk()) continue;
275    
276     Double_t deltaZ = 0.0;
277     if(pf->BestTrk()) {
278     deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
279     }
280    
281     // ignore the pf candidate if it is too far away in Z
282     if (deltaZ >= delta_z)
283     continue;
284    
285     Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
286     // add the pf pt if it is inside the extRadius and outside the intRadius
287     if ( dr < extRadius &&
288     dr >= intRadius ) {
289    
290     //EtaStrip Veto for Gamma
291     if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;
292    
293     //InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
294     if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
295     && MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue;
296    
297     ptSum += pf->Pt();
298    
299     }
300     }
301     return ptSum;
302     }
303    
304    
305     //--------------------------------------------------------------------------------------------------
306     Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands,
307     const Vertex *vertex, const MuonCol *goodMuons,
308     const ElectronCol *goodElectrons,
309     Double_t delta_z, Double_t ptMin,
310 ceballos 1.8 Double_t extRadius, Double_t intRadius, int isoType,
311 sixie 1.16 Double_t beta)
312 ceballos 1.5 {
313     //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
314     //annulus around the particle seed track.
315    
316     Double_t zLepton = 0.0;
317 ceballos 1.8 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
318 ceballos 1.5
319     Double_t ptSum =0.;
320     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
321     const PFCandidate *pf = PFCands->At(i);
322    
323     Bool_t isGoodType = kFALSE;
324     // all particles
325 ceballos 1.6 if (isoType == 0) isGoodType = kTRUE;
326 ceballos 1.5 // charged particles only
327 ceballos 1.6 else if(isoType == 1 && pf->BestTrk()) isGoodType = kTRUE;
328 ceballos 1.5 // charged particles and gammas only
329 ceballos 1.6 else if(isoType == 2 &&
330     (pf->BestTrk() || pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
331 ceballos 1.8 // all particles, rejecting good leptons
332     else if(isoType == 3) isGoodType = kTRUE;
333 ceballos 1.5
334     if(isGoodType == kFALSE) continue;
335    
336 sixie 1.14 // pt cut applied to neutrals
337     if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
338 ceballos 1.6
339 ceballos 1.5 if(pf->TrackerTrk() && p->TrackerTrk() &&
340     pf->TrackerTrk() == p->TrackerTrk()) continue;
341    
342     if(pf->GsfTrk() && p->GsfTrk() &&
343     pf->GsfTrk() == p->GsfTrk()) continue;
344    
345     Double_t deltaZ = 0.0;
346     if(pf->BestTrk()) {
347 ceballos 1.8 deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
348 ceballos 1.5 }
349    
350     // ignore the pf candidate if it is too far away in Z
351 ceballos 1.6 if (deltaZ >= delta_z)
352 ceballos 1.5 continue;
353    
354     Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
355     // add the pf pt if it is inside the extRadius and outside the intRadius
356     if ( dr < extRadius &&
357     dr >= intRadius ) {
358 ceballos 1.8 Bool_t isLepton = kFALSE;
359     if(goodMuons && isoType == 3){
360     for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
361     const Muon *m = goodMuons->At(nl);
362     if(pf->TrackerTrk() && m->TrackerTrk() &&
363     pf->TrackerTrk() == m->TrackerTrk()) {
364     isLepton = kTRUE;
365     break;
366     }
367     }
368     }
369     if(goodElectrons && isLepton == kFALSE && isoType == 3){
370     for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
371     const Electron *e = goodElectrons->At(nl);
372     if(pf->TrackerTrk() && e->TrackerTrk() &&
373     pf->TrackerTrk() == e->TrackerTrk()) {
374     isLepton = kTRUE;
375     break;
376     }
377     if(pf->GsfTrk() && e->GsfTrk() &&
378     pf->GsfTrk() == e->GsfTrk()) {
379     isLepton = kTRUE;
380     break;
381     }
382     }
383     }
384 sixie 1.14
385     if (isLepton == kTRUE) continue;
386    
387     //EtaStrip Veto for Gamma
388     if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;
389    
390     //InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
391     if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
392     && MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue;
393    
394    
395     if(pf->BestTrk()) ptSum += pf->Pt();
396     else ptSum += pf->Pt()*beta;
397    
398    
399 ceballos 1.5 }
400     }
401     return ptSum;
402     }
403 ceballos 1.7 //--------------------------------------------------------------------------------------------------
404 ceballos 1.8 Double_t IsolationTools::BetaM(const TrackCol *tracks, const Muon *p, const Vertex *vertex,
405 ceballos 1.7 Double_t ptMin, Double_t delta_z, Double_t extRadius,
406     Double_t intRadius){
407    
408 ceballos 1.9 if(!tracks) return 1.0;
409 ceballos 1.7 if(tracks->GetEntries() <= 0) return 1.0;
410    
411     double Pt_jets_X = 0. ;
412     double Pt_jets_Y = 0. ;
413     double Pt_jets_X_tot = 0. ;
414     double Pt_jets_Y_tot = 0. ;
415    
416     for(int i=0;i<int(tracks->GetEntries());i++){
417     const Track *pTrack = tracks->At(i);
418    
419     if(pTrack && p->TrackerTrk() &&
420     pTrack == p->TrackerTrk()) continue;
421    
422     if(pTrack->Pt() <= ptMin) continue;
423    
424     Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
425     if ( dr < extRadius && dr >= intRadius ) {
426     Pt_jets_X_tot += pTrack->Px();
427     Pt_jets_Y_tot += pTrack->Py();
428 ceballos 1.8 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
429 ceballos 1.7 if(pDz < delta_z){
430     Pt_jets_X += pTrack->Px();
431     Pt_jets_Y += pTrack->Py();
432     }
433     }
434     }
435    
436     if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
437     return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);
438    
439     return 1.0;
440     }
441    
442     //--------------------------------------------------------------------------------------------------
443 ceballos 1.8 Double_t IsolationTools::BetaE(const TrackCol *tracks, const Electron *p, const Vertex *vertex,
444 ceballos 1.7 Double_t ptMin, Double_t delta_z, Double_t extRadius,
445     Double_t intRadius){
446    
447     if(!tracks) return 1.0;
448     if(tracks->GetEntries() <= 0) return 1.0;
449    
450     double Pt_jets_X = 0. ;
451     double Pt_jets_Y = 0. ;
452     double Pt_jets_X_tot = 0. ;
453     double Pt_jets_Y_tot = 0. ;
454    
455     for(int i=0;i<int(tracks->GetEntries());i++){
456     const Track *pTrack = tracks->At(i);
457    
458     if(pTrack && p->TrackerTrk() &&
459     pTrack == p->TrackerTrk()) continue;
460    
461     if(pTrack && p->GsfTrk() &&
462     pTrack == p->GsfTrk()) continue;
463    
464     if(pTrack->Pt() <= ptMin) continue;
465    
466     Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
467     if ( dr < extRadius && dr >= intRadius ) {
468     Pt_jets_X_tot += pTrack->Px();
469     Pt_jets_Y_tot += pTrack->Py();
470 ceballos 1.8 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
471 ceballos 1.7 if(pDz < delta_z){
472     Pt_jets_X += pTrack->Px();
473     Pt_jets_Y += pTrack->Py();
474     }
475     }
476     }
477    
478     if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
479     return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);
480    
481     return 1.0;
482     }
483 fabstoec 1.10
484    
485     // method added by F.Stoeckli: computes the track isolation with NO constrint on the OV-track compatibility
486     Double_t IsolationTools::TrackIsolationNoPV(const mithep::Particle* p, const BaseVertex* bsp,
487     Double_t extRadius,
488     Double_t intRadius,
489     Double_t ptLow,
490     Double_t etaStrip,
491     Double_t maxD0,
492     mithep::TrackQuality::EQuality quality,
493 bendavid 1.11 const mithep::Collection<mithep::Track> *tracks,
494     UInt_t maxNExpectedHitsInner,
495     const mithep::DecayParticleCol *conversions) {
496 fabstoec 1.10
497     // loop over all tracks
498     Double_t tPt = 0.;
499     for(UInt_t i=0; i<tracks->GetEntries(); ++i) {
500     const Track* t = tracks->At(i);
501     if ( t->Pt() < ptLow ) continue;
502     if ( ! t->Quality().Quality(quality) ) continue;
503     // only check for beamspot if available, otherwise ignore cut
504     if ( bsp && fabs(t->D0Corrected( *bsp) ) > maxD0) continue;
505 bendavid 1.11 if (t->NExpectedHitsInner()>maxNExpectedHitsInner) continue;
506     if (conversions && PhotonTools::MatchedConversion(t,conversions,bsp)) continue;
507 fabstoec 1.10 Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
508     Double_t dEta = fabs(t->Eta()-p->Eta());
509     if(dR < extRadius && dR > intRadius && dEta > etaStrip) tPt += t->Pt();
510     }
511     return tPt;
512     }
513    
514 fabstoec 1.19
515 fabstoec 1.20 Double_t IsolationTools::CiCTrackIsolation(const mithep::Photon* p,
516 fabstoec 1.19 const BaseVertex* theVtx,
517     Double_t extRadius,
518     Double_t intRadius,
519     Double_t ptLow,
520     Double_t etaStrip,
521     Double_t maxD0,
522     Double_t maxDZ,
523     const mithep::Collection<mithep::Track> *tracks,
524 fabstoec 1.20 unsigned int* worstVtxIndex,
525     const mithep::Collection<mithep::Vertex> *vtxs,
526     bool print) {
527 fabstoec 1.19
528     UInt_t numVtx = 1;
529     const BaseVertex* iVtx = theVtx;
530     if( vtxs ) {
531     numVtx = vtxs->GetEntries();
532     if (numVtx > 0)
533     iVtx = vtxs->At(0);
534     else
535     return 0.;
536     }
537    
538 fabstoec 1.20
539     if(print) {
540     std::cout<<" Testing photon with"<<std::endl;
541     std::cout<<" Et = "<<p->Et()<<std::endl;
542     std::cout<<" Eta = "<<p->Eta()<<std::endl;
543     std::cout<<" Phi = "<<p->Phi()<<std::endl;
544     }
545    
546 fabstoec 1.19 Double_t iIso = 0.;
547     Double_t maxIso = 0.;
548 fabstoec 1.20
549     if(worstVtxIndex)
550     *worstVtxIndex=0;
551 fabstoec 1.19
552     for(UInt_t i=0; i<numVtx; ++i) {
553 fabstoec 1.20
554     if(i>0) iVtx = vtxs->At(i);
555    
556    
557     if(print) {
558     std::cout<<" Vertex #"<<i<<std::endl;
559     std::cout<<" with X = "<<iVtx->X()<<std::endl;
560     std::cout<<" with Y = "<<iVtx->Y()<<std::endl;
561     std::cout<<" with Z = "<<iVtx->Z()<<std::endl;
562     }
563    
564     // compute the ph momentum with respect to this Vtx
565     FourVectorM phMom = p->MomVtx(iVtx->Position());
566    
567     if(print) {
568     std::cout<<" photon has changed to:"<<std::endl;
569     std::cout<<" Et = "<<phMom.Et()<<std::endl;
570     std::cout<<" eta = "<<phMom.Eta()<<std::endl;
571     std::cout<<" Phi = "<<phMom.Phi()<<std::endl;
572     }
573    
574 fabstoec 1.19 iIso = 0.;
575 fabstoec 1.20 for(UInt_t j=0; j<tracks->GetEntries(); ++j) {
576     const Track* t = tracks->At(j);
577    
578     //Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
579     //Double_t dEta = TMath::Abs(t->Eta()-p->Eta());
580    
581     Double_t dR = MathUtils::DeltaR(t->Mom(),phMom);
582     Double_t dEta = TMath::Abs(t->Eta()-phMom.Eta());
583    
584     if(print) {
585     std::cout<<" passing track #"<<j<<std::endl;
586     std::cout<<" pt = "<<t->Pt()<<std::endl;
587     std::cout<<" eta = "<<t->Eta()<<std::endl;
588     std::cout<<" phi = "<<t->Phi()<<std::endl;
589     std::cout<<" d0 = "<<fabs(t->D0Corrected( *iVtx ))<<std::endl;
590     std::cout<<" dZ = "<<fabs(t->DzCorrected( *iVtx ))<<std::endl;
591     std::cout<<" dR = "<<dR<<std::endl;
592     std::cout<<" dEta = "<<dEta<<std::endl;
593     std::cout<<" vx = "<<t->X0()<<std::endl;
594     std::cout<<" vy = "<<t->Y0()<<std::endl;
595     std::cout<<" vz = "<<t->Z0()<<std::endl;
596     }
597    
598 fabstoec 1.19 if ( t->Pt() < ptLow ) continue;
599     // only check for beamspot if available, otherwise ignore cut
600     if ( fabs(t->D0Corrected( *iVtx )) > maxD0) continue;
601     if ( fabs(t->DzCorrected( *iVtx )) > maxDZ) continue;
602 fabstoec 1.20
603    
604     if(dR < extRadius && dR > intRadius && dEta > etaStrip) {
605     iIso += t->Pt();
606    
607     if(print) {
608     std::cout<<" passing track #"<<j<<std::endl;
609     std::cout<<" pt = "<<t->Pt()<<std::endl;
610     std::cout<<" eta = "<<t->Eta()<<std::endl;
611     std::cout<<" phi = "<<t->Phi()<<std::endl;
612     std::cout<<" d0 = "<<fabs(t->D0Corrected( *iVtx ))<<std::endl;
613     std::cout<<" dZ = "<<fabs(t->DzCorrected( *iVtx ))<<std::endl;
614     std::cout<<" dR = "<<dR<<std::endl;
615     std::cout<<" dEta = "<<dEta<<std::endl;
616     std::cout<<" vx = "<<t->X0()<<std::endl;
617     std::cout<<" vy = "<<t->Y0()<<std::endl;
618     std::cout<<" vz = "<<t->Z0()<<std::endl;
619     std::cout<<" new tIso = "<<iIso<<std::endl;
620     }
621     }
622     }
623     if ( iIso > maxIso ) {
624     maxIso = iIso;
625     if(worstVtxIndex)
626     *worstVtxIndex=i;
627 fabstoec 1.19 }
628     }
629 fabstoec 1.20
630     if(print) {
631     if(worstVtxIndex)
632     std::cout<<" max TrkIso is given by Vtx #"<<*worstVtxIndex<<" with an amount of tIso = "<<maxIso<<std::endl;
633     else
634     std::cout<<" max TrkIso is given by Vtx #0 with an amount of tIso = "<<maxIso<<std::endl;
635     }
636     return maxIso;
637 fabstoec 1.19 }