ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/IsolationTools.cc
Revision: 1.22
Committed: Tue Oct 18 11:27:19 2011 UTC (13 years, 6 months ago) by fabstoec
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025c, Mit_025b, Mit_025a
Changes since 1.21: +47 -5 lines
Log Message: 
Added electron Veto inversion to CiCTrackIso & MVA Tools

File Contents

# User Rev Content
1 fabstoec 1.22 // $Id: IsolationTools.cc,v 1.21 2011/09/16 14:09:34 ceballos 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 fabstoec 1.22 //std::cout<<" *** TrackIso:"<<std::endl;
500 fabstoec 1.10 for(UInt_t i=0; i<tracks->GetEntries(); ++i) {
501     const Track* t = tracks->At(i);
502 fabstoec 1.22 if(t->Pt()>1. && false) {
503     std::cout<<" "<<i<<" pt = "<<t->Pt()<<" ("<<ptLow<<")"<<std::endl;
504     std::cout<<" d0 = "<<fabs(t->D0Corrected( *bsp) )<<" ("<<maxD0<<")"<<std::endl;
505     //std::cout<<" conv ? "<<PhotonTools::MatchedConversion(t,conversions,bsp)<<std::endl;
506     std::cout<<" dR = "<<MathUtils::DeltaR(t->Mom(),p->Mom())<<" ("<<extRadius<<","<<intRadius<<")"<<std::endl;
507     std::cout<<" dEta = "<<fabs(t->Eta()-p->Eta())<<" ("<<etaStrip<<")"<<std::endl;
508     }
509 fabstoec 1.10 if ( t->Pt() < ptLow ) continue;
510     if ( ! t->Quality().Quality(quality) ) continue;
511     // only check for beamspot if available, otherwise ignore cut
512     if ( bsp && fabs(t->D0Corrected( *bsp) ) > maxD0) continue;
513 bendavid 1.11 if (t->NExpectedHitsInner()>maxNExpectedHitsInner) continue;
514     if (conversions && PhotonTools::MatchedConversion(t,conversions,bsp)) continue;
515 fabstoec 1.10 Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
516     Double_t dEta = fabs(t->Eta()-p->Eta());
517     if(dR < extRadius && dR > intRadius && dEta > etaStrip) tPt += t->Pt();
518     }
519     return tPt;
520     }
521    
522 fabstoec 1.19
523 fabstoec 1.20 Double_t IsolationTools::CiCTrackIsolation(const mithep::Photon* p,
524 fabstoec 1.19 const BaseVertex* theVtx,
525     Double_t extRadius,
526     Double_t intRadius,
527     Double_t ptLow,
528     Double_t etaStrip,
529     Double_t maxD0,
530     Double_t maxDZ,
531     const mithep::Collection<mithep::Track> *tracks,
532 fabstoec 1.20 unsigned int* worstVtxIndex,
533     const mithep::Collection<mithep::Vertex> *vtxs,
534 fabstoec 1.22 const mithep::Collection<mithep::Electron> *eles,
535     bool print,
536     double* ptmax, double* dRmax) {
537 fabstoec 1.19
538     UInt_t numVtx = 1;
539     const BaseVertex* iVtx = theVtx;
540     if( vtxs ) {
541     numVtx = vtxs->GetEntries();
542     if (numVtx > 0)
543     iVtx = vtxs->At(0);
544     else
545     return 0.;
546     }
547    
548 fabstoec 1.22 // NEW for Electron T&P: need to remove the electron Gsf Track (applied if eles != NULL)
549     const Track* theGsfTrack = NULL;
550     if ( eles ) {
551     // find the electron that matches the Photon SC
552     for(UInt_t j=0; j<eles->GetEntries(); ++j) {
553     if ( eles->At(j)->SCluster() == p->SCluster() ) {
554     if( eles->At(j)->HasTrackerTrk() )
555     theGsfTrack = eles->At(j)->TrackerTrk();
556     break;
557     }
558     }
559     }
560 fabstoec 1.20
561     if(print) {
562     std::cout<<" Testing photon with"<<std::endl;
563     std::cout<<" Et = "<<p->Et()<<std::endl;
564     std::cout<<" Eta = "<<p->Eta()<<std::endl;
565     std::cout<<" Phi = "<<p->Phi()<<std::endl;
566     }
567    
568 fabstoec 1.19 Double_t iIso = 0.;
569     Double_t maxIso = 0.;
570 fabstoec 1.20
571     if(worstVtxIndex)
572     *worstVtxIndex=0;
573 fabstoec 1.19
574 fabstoec 1.22 double t_ptmax = 0.;
575     double t_dRmax = 0.;
576    
577 fabstoec 1.19 for(UInt_t i=0; i<numVtx; ++i) {
578 fabstoec 1.20
579     if(i>0) iVtx = vtxs->At(i);
580    
581    
582     if(print) {
583     std::cout<<" Vertex #"<<i<<std::endl;
584     std::cout<<" with X = "<<iVtx->X()<<std::endl;
585     std::cout<<" with Y = "<<iVtx->Y()<<std::endl;
586     std::cout<<" with Z = "<<iVtx->Z()<<std::endl;
587     }
588    
589 fabstoec 1.22 Photon* phTemp = new Photon(*p);
590    
591     // RESET CALO_POS!
592     phTemp->SetCaloPosXYZ(p->SCluster()->Point().X(),p->SCluster()->Point().Y(),p->SCluster()->Point().Z());
593    
594 fabstoec 1.20 // compute the ph momentum with respect to this Vtx
595 fabstoec 1.22 //FourVectorM phMom = p->MomVtx(iVtx->Position());
596     FourVectorM phMom = phTemp->MomVtx(iVtx->Position());
597    
598     delete phTemp;
599 fabstoec 1.20
600     if(print) {
601     std::cout<<" photon has changed to:"<<std::endl;
602     std::cout<<" Et = "<<phMom.Et()<<std::endl;
603     std::cout<<" eta = "<<phMom.Eta()<<std::endl;
604     std::cout<<" Phi = "<<phMom.Phi()<<std::endl;
605     }
606    
607 fabstoec 1.19 iIso = 0.;
608 fabstoec 1.20 for(UInt_t j=0; j<tracks->GetEntries(); ++j) {
609     const Track* t = tracks->At(j);
610 fabstoec 1.22 if( theGsfTrack && t == theGsfTrack ) continue;
611 fabstoec 1.20
612     //Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
613     //Double_t dEta = TMath::Abs(t->Eta()-p->Eta());
614    
615     Double_t dR = MathUtils::DeltaR(t->Mom(),phMom);
616     Double_t dEta = TMath::Abs(t->Eta()-phMom.Eta());
617    
618 fabstoec 1.22 if(print && t->Pt()>1. && false) {
619 fabstoec 1.20 std::cout<<" passing track #"<<j<<std::endl;
620     std::cout<<" pt = "<<t->Pt()<<std::endl;
621     std::cout<<" eta = "<<t->Eta()<<std::endl;
622     std::cout<<" phi = "<<t->Phi()<<std::endl;
623     std::cout<<" d0 = "<<fabs(t->D0Corrected( *iVtx ))<<std::endl;
624     std::cout<<" dZ = "<<fabs(t->DzCorrected( *iVtx ))<<std::endl;
625     std::cout<<" dR = "<<dR<<std::endl;
626     std::cout<<" dEta = "<<dEta<<std::endl;
627     std::cout<<" vx = "<<t->X0()<<std::endl;
628     std::cout<<" vy = "<<t->Y0()<<std::endl;
629     std::cout<<" vz = "<<t->Z0()<<std::endl;
630     }
631    
632 fabstoec 1.19 if ( t->Pt() < ptLow ) continue;
633     // only check for beamspot if available, otherwise ignore cut
634     if ( fabs(t->D0Corrected( *iVtx )) > maxD0) continue;
635     if ( fabs(t->DzCorrected( *iVtx )) > maxDZ) continue;
636 fabstoec 1.20
637    
638     if(dR < extRadius && dR > intRadius && dEta > etaStrip) {
639     iIso += t->Pt();
640    
641 fabstoec 1.22 if(t->Pt() > t_ptmax) {
642     t_ptmax=t->Pt();
643     t_dRmax=dR;
644     }
645    
646     if(print && t->Pt()>1.) {
647 fabstoec 1.20 std::cout<<" passing track #"<<j<<std::endl;
648     std::cout<<" pt = "<<t->Pt()<<std::endl;
649     std::cout<<" eta = "<<t->Eta()<<std::endl;
650     std::cout<<" phi = "<<t->Phi()<<std::endl;
651     std::cout<<" d0 = "<<fabs(t->D0Corrected( *iVtx ))<<std::endl;
652     std::cout<<" dZ = "<<fabs(t->DzCorrected( *iVtx ))<<std::endl;
653     std::cout<<" dR = "<<dR<<std::endl;
654     std::cout<<" dEta = "<<dEta<<std::endl;
655     std::cout<<" vx = "<<t->X0()<<std::endl;
656     std::cout<<" vy = "<<t->Y0()<<std::endl;
657     std::cout<<" vz = "<<t->Z0()<<std::endl;
658     std::cout<<" new tIso = "<<iIso<<std::endl;
659     }
660     }
661     }
662     if ( iIso > maxIso ) {
663     maxIso = iIso;
664     if(worstVtxIndex)
665     *worstVtxIndex=i;
666 fabstoec 1.19 }
667     }
668 fabstoec 1.20
669 fabstoec 1.22 if(ptmax) (*ptmax)=t_ptmax;
670     if(dRmax) (*dRmax)=t_dRmax;
671    
672 fabstoec 1.20 if(print) {
673     if(worstVtxIndex)
674     std::cout<<" max TrkIso is given by Vtx #"<<*worstVtxIndex<<" with an amount of tIso = "<<maxIso<<std::endl;
675     else
676     std::cout<<" max TrkIso is given by Vtx #0 with an amount of tIso = "<<maxIso<<std::endl;
677     }
678     return maxIso;
679 fabstoec 1.19 }