ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/IsolationTools.cc
Revision: 1.27
Committed: Sat Apr 28 11:34:15 2012 UTC (13 years ago) by ceballos
Content type: text/plain
Branch: MAIN
Changes since 1.26: +103 -1 lines
Log Message:
new 2012 isolation

File Contents

# User Rev Content
1 ceballos 1.27 // $Id: IsolationTools.cc,v 1.26 2011/12/31 23:16:13 sixie 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 ceballos 1.23 Double_t IsolationTools::PFMuonIsolation(const Muon *p, const PFCandidateCol *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     if(pf->HasTrk()) {
147 sixie 1.24 Double_t deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
148     if (deltaZ >= delta_z) continue;
149 mzanetti 1.15 }
150 mzanetti 1.18
151     // inner cone veto for gammas
152 ceballos 1.21 if (pf->PFType() == PFCandidate::eGamma && dr < intRadiusGamma) continue;
153 mzanetti 1.18
154 ceballos 1.21 // inner cone veto for tracks
155     if (dr < intRadius) continue;
156    
157 mzanetti 1.18 // add the pf pt if it is inside the extRadius
158 ceballos 1.21 if (dr < extRadius) ptSum += pf->Pt();
159 mzanetti 1.15
160     }
161     return ptSum;
162     }
163    
164     //--------------------------------------------------------------------------------------------------
165 ceballos 1.23 Double_t IsolationTools::PFMuonIsolation(const Muon *p, const PFCandidateCol *PFCands,
166     const MuonCol *goodMuons, const ElectronCol *goodElectrons,
167     const Vertex *vertex, Double_t delta_z, Double_t ptMin,
168     Double_t extRadius, Double_t intRadiusGamma, Double_t intRadius)
169 ceballos 1.5 {
170     //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
171     //annulus around the particle seed track.
172    
173     Double_t zLepton = 0.0;
174 ceballos 1.8 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
175 ceballos 1.5
176     Double_t ptSum =0.;
177     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
178     const PFCandidate *pf = PFCands->At(i);
179    
180 ceballos 1.23 // exclude muon
181     if(pf->TrackerTrk() && p->TrackerTrk() &&
182     pf->TrackerTrk() == p->TrackerTrk()) continue;
183 ceballos 1.5
184 ceballos 1.23 Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
185    
186 mzanetti 1.13 // pt cut applied to neutrals
187     if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
188 ceballos 1.6
189 ceballos 1.23 // ignore the pf candidate if it is too far away in Z
190     if(pf->HasTrk()) {
191 sixie 1.24 Double_t deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
192     if (deltaZ >= delta_z) continue;
193 ceballos 1.5 }
194 ceballos 1.23
195 mzanetti 1.18 // inner cone veto for gammas
196 ceballos 1.23 if (pf->PFType() == PFCandidate::eGamma && dr < intRadiusGamma) continue;
197    
198     // inner cone veto for tracks
199     if (dr < intRadius) continue;
200 mzanetti 1.18
201 ceballos 1.5 // add the pf pt if it is inside the extRadius and outside the intRadius
202 mzanetti 1.18 if (dr < extRadius ) {
203 ceballos 1.8 Bool_t isLepton = kFALSE;
204 ceballos 1.23 if(goodMuons){
205 ceballos 1.8 for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
206     const Muon *m = goodMuons->At(nl);
207     if(pf->TrackerTrk() && m->TrackerTrk() &&
208     pf->TrackerTrk() == m->TrackerTrk()) {
209     isLepton = kTRUE;
210     break;
211     }
212     }
213     }
214 ceballos 1.23 if(goodElectrons && isLepton == kFALSE){
215 ceballos 1.8 for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
216     const Electron *e = goodElectrons->At(nl);
217     if(pf->TrackerTrk() && e->TrackerTrk() &&
218     pf->TrackerTrk() == e->TrackerTrk()) {
219     isLepton = kTRUE;
220     break;
221     }
222     if(pf->GsfTrk() && e->GsfTrk() &&
223     pf->GsfTrk() == e->GsfTrk()) {
224     isLepton = kTRUE;
225     break;
226     }
227     }
228     }
229 ceballos 1.23 if (isLepton == kTRUE) continue;
230     ptSum += pf->Pt();
231     }
232 ceballos 1.5 }
233     return ptSum;
234     }
235 mzanetti 1.15
236 ceballos 1.5 //--------------------------------------------------------------------------------------------------
237     Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands,
238 ceballos 1.8 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
239 sixie 1.26 Double_t extRadius, Double_t intRadius, Int_t PFCandidateType)
240 sixie 1.16 {
241    
242     //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
243     //annulus around the particle seed track.
244    
245     Double_t zLepton = 0.0;
246     if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
247    
248     Double_t ptSum =0.;
249     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
250     const PFCandidate *pf = PFCands->At(i);
251    
252 sixie 1.26 //select only specified PFCandidate types
253     if (PFCandidateType >= 0) {
254     if (pf->PFType() != PFCandidateType) continue;
255     }
256    
257 sixie 1.16 // pt cut applied to neutrals
258     if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
259    
260     if(pf->TrackerTrk() && p->TrackerTrk() &&
261     pf->TrackerTrk() == p->TrackerTrk()) continue;
262    
263     if(pf->GsfTrk() && p->GsfTrk() &&
264     pf->GsfTrk() == p->GsfTrk()) continue;
265    
266 sixie 1.25 // ignore the pf candidate if it is too far away in Z
267 sixie 1.16 if(pf->BestTrk()) {
268 sixie 1.25 Double_t deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
269     if (deltaZ >= delta_z) continue;
270 sixie 1.16 }
271    
272    
273     Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
274     // add the pf pt if it is inside the extRadius and outside the intRadius
275     if ( dr < extRadius &&
276     dr >= intRadius ) {
277    
278     //EtaStrip Veto for Gamma
279     if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;
280    
281     //InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
282     if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
283     && MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue;
284    
285     ptSum += pf->Pt();
286    
287     }
288     }
289     return ptSum;
290     }
291    
292    
293     //--------------------------------------------------------------------------------------------------
294     Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands,
295 ceballos 1.23 const MuonCol *goodMuons, const ElectronCol *goodElectrons,
296     const Vertex *vertex, Double_t delta_z, Double_t ptMin,
297 sixie 1.26 Double_t extRadius, Double_t intRadius, Int_t PFCandidateType)
298 ceballos 1.5 {
299 ceballos 1.23
300 ceballos 1.5 //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
301     //annulus around the particle seed track.
302    
303     Double_t zLepton = 0.0;
304 ceballos 1.8 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
305 ceballos 1.5
306     Double_t ptSum =0.;
307     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
308     const PFCandidate *pf = PFCands->At(i);
309    
310 sixie 1.26 //select only specified PFCandidate types
311     if (PFCandidateType >= 0) {
312     if (pf->PFType() != PFCandidateType) continue;
313     }
314    
315 sixie 1.14 // pt cut applied to neutrals
316     if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
317 ceballos 1.6
318 ceballos 1.5 if(pf->TrackerTrk() && p->TrackerTrk() &&
319     pf->TrackerTrk() == p->TrackerTrk()) continue;
320    
321     if(pf->GsfTrk() && p->GsfTrk() &&
322     pf->GsfTrk() == p->GsfTrk()) continue;
323    
324 sixie 1.25 // ignore the pf candidate if it is too far away in Z
325 ceballos 1.5 if(pf->BestTrk()) {
326 sixie 1.25 Double_t deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
327     if (deltaZ >= delta_z) continue;
328 ceballos 1.5 }
329    
330     Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
331     // add the pf pt if it is inside the extRadius and outside the intRadius
332     if ( dr < extRadius &&
333     dr >= intRadius ) {
334 ceballos 1.23
335     //EtaStrip Veto for Gamma
336     if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;
337    
338     //InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
339     if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
340     && MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue;
341    
342 ceballos 1.8 Bool_t isLepton = kFALSE;
343 ceballos 1.23 if(goodMuons){
344 ceballos 1.8 for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
345     const Muon *m = goodMuons->At(nl);
346     if(pf->TrackerTrk() && m->TrackerTrk() &&
347     pf->TrackerTrk() == m->TrackerTrk()) {
348     isLepton = kTRUE;
349     break;
350     }
351     }
352     }
353 ceballos 1.23 if(goodElectrons && isLepton == kFALSE){
354 ceballos 1.8 for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
355     const Electron *e = goodElectrons->At(nl);
356     if(pf->TrackerTrk() && e->TrackerTrk() &&
357     pf->TrackerTrk() == e->TrackerTrk()) {
358     isLepton = kTRUE;
359     break;
360     }
361     if(pf->GsfTrk() && e->GsfTrk() &&
362     pf->GsfTrk() == e->GsfTrk()) {
363     isLepton = kTRUE;
364     break;
365     }
366     }
367     }
368 sixie 1.14
369     if (isLepton == kTRUE) continue;
370 ceballos 1.23 ptSum += pf->Pt();
371 sixie 1.14
372 ceballos 1.5 }
373     }
374     return ptSum;
375     }
376 ceballos 1.7 //--------------------------------------------------------------------------------------------------
377 ceballos 1.27 Double_t IsolationTools::PFElectronIsolation2012(const Electron *ele, const Vertex *vertex,
378     const PFCandidateCol *PFCands,
379     const PileupEnergyDensityCol *PileupEnergyDensity,
380     ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaTarget,
381     const ElectronCol *goodElectrons,
382     const MuonCol *goodMuons, Double_t dRMax){
383    
384     Double_t tmpChargedIso_DR = 0;
385     Double_t tmpGammaIso_DR = 0;
386     Double_t tmpNeutralHadronIso_DR = 0;
387    
388     for (UInt_t p=0; p<PFCands->GetEntries();p++) {
389     const PFCandidate *pf = PFCands->At(p);
390    
391     //exclude the electron itself
392     if(pf->GsfTrk() && ele->GsfTrk() &&
393     pf->GsfTrk() == ele->GsfTrk()) continue;
394     if(pf->TrackerTrk() && ele->TrackerTrk() &&
395     pf->TrackerTrk() == ele->TrackerTrk()) continue;
396    
397     //************************************************************
398     // New Isolation Calculations
399     //************************************************************
400     Double_t dr = MathUtils::DeltaR(ele->Mom(), pf->Mom());
401    
402     if (dr < 1.0) {
403     Bool_t IsLeptonFootprint = kFALSE;
404     //************************************************************
405     // Lepton Footprint Removal
406     //************************************************************
407     for (UInt_t q=0; q < goodElectrons->GetEntries() ; ++q) {
408     //if pf candidate matches an electron passing ID cuts, then veto it
409     if(pf->GsfTrk() && goodElectrons->At(q)->GsfTrk() &&
410     pf->GsfTrk() == goodElectrons->At(q)->GsfTrk()) IsLeptonFootprint = kTRUE;
411     if(pf->TrackerTrk() && goodElectrons->At(q)->TrackerTrk() &&
412     pf->TrackerTrk() == goodElectrons->At(q)->TrackerTrk()) IsLeptonFootprint = kTRUE;
413     //if pf candidate lies in veto regions of electron passing ID cuts, then veto it
414     if(pf->BestTrk() && fabs(goodElectrons->At(q)->SCluster()->Eta()) >= 1.479
415     && MathUtils::DeltaR(goodElectrons->At(q)->Mom(), pf->Mom()) < 0.015) IsLeptonFootprint = kTRUE;
416     if(pf->PFType() == PFCandidate::eGamma && fabs(goodElectrons->At(q)->SCluster()->Eta()) >= 1.479 &&
417     MathUtils::DeltaR(goodElectrons->At(q)->Mom(), pf->Mom()) < 0.08) IsLeptonFootprint = kTRUE;
418     }
419     for (UInt_t q=0; q < goodMuons->GetEntries() ; ++q) {
420     //if pf candidate matches an muon passing ID cuts, then veto it
421     if(pf->TrackerTrk() && goodMuons->At(q)->TrackerTrk() &&
422     pf->TrackerTrk() == goodMuons->At(q)->TrackerTrk()) IsLeptonFootprint = kTRUE;
423     //if pf candidate lies in veto regions of muon passing ID cuts, then veto it
424     if(pf->BestTrk() && MathUtils::DeltaR(goodMuons->At(q)->Mom(), pf->Mom()) < 0.01) IsLeptonFootprint = kTRUE;
425     }
426    
427     if (!IsLeptonFootprint) {
428     Bool_t passVeto = kTRUE;
429     //Charged
430     if(pf->BestTrk()) {
431     // CMS DOESN"T WANT THIS
432     //if (!(fabs(pf->BestTrk()->DzCorrected(*vertex) - ele->BestTrk()->DzCorrected(*vertex)) < 0.2)) passVeto = kFALSE;
433     //************************************************************
434     // Veto any PFmuon, or PFEle
435     if (pf->PFType() == PFCandidate::eElectron || pf->PFType() == PFCandidate::eMuon) passVeto = kFALSE;
436     //************************************************************
437     //************************************************************
438     // Footprint Veto
439     if (fabs(ele->SCluster()->Eta()) >= 1.479 && dr < 0.015) passVeto = kFALSE;
440     //************************************************************
441     if (passVeto) {
442     if (dr < dRMax) tmpChargedIso_DR += pf->Pt();
443     } //pass veto
444    
445     }
446     //Gamma
447     else if (pf->PFType() == PFCandidate::eGamma) {
448     //************************************************************
449     // Footprint Veto
450     if (fabs(ele->SCluster()->Eta()) >= 1.479) {
451     if (dr < 0.08) passVeto = kFALSE;
452     }
453     //************************************************************
454    
455     if (passVeto) {
456     if (dr < dRMax) tmpGammaIso_DR += pf->Pt();
457     }
458     }
459     //NeutralHadron
460     else {
461     if (dr < dRMax) tmpNeutralHadronIso_DR += pf->Pt();
462     }
463     } //not lepton footprint
464     } //in 1.0 dr cone
465     } //loop over PF candidates
466    
467     Double_t Rho = 0;
468     if (!(TMath::IsNaN(PileupEnergyDensity->At(0)->Rho()) || isinf(PileupEnergyDensity->At(0)->Rho()))) Rho = PileupEnergyDensity->At(0)->Rho();
469    
470     Double_t IsoVar_ChargedIso_DR = tmpChargedIso_DR/ele->Pt();
471     Double_t IsoVar_NeutralIso_DR = tmpGammaIso_DR + tmpNeutralHadronIso_DR;
472     // Careful here, we have kEleNeutralIso04 only for now
473     if(dRMax != 0.4) assert(0);
474     IsoVar_NeutralIso_DR = TMath::Max((IsoVar_NeutralIso_DR - Rho*ElectronTools::ElectronEffectiveArea(ElectronTools::kEleNeutralIso04, ele->SCluster()->Eta(), EffectiveAreaTarget))/ele->Pt(), 0.0);
475    
476     return (IsoVar_ChargedIso_DR + IsoVar_NeutralIso_DR);
477     }
478     //--------------------------------------------------------------------------------------------------
479 ceballos 1.8 Double_t IsolationTools::BetaM(const TrackCol *tracks, const Muon *p, const Vertex *vertex,
480 ceballos 1.7 Double_t ptMin, Double_t delta_z, Double_t extRadius,
481     Double_t intRadius){
482    
483 ceballos 1.9 if(!tracks) return 1.0;
484 ceballos 1.7 if(tracks->GetEntries() <= 0) return 1.0;
485    
486     double Pt_jets_X = 0. ;
487     double Pt_jets_Y = 0. ;
488     double Pt_jets_X_tot = 0. ;
489     double Pt_jets_Y_tot = 0. ;
490    
491     for(int i=0;i<int(tracks->GetEntries());i++){
492     const Track *pTrack = tracks->At(i);
493    
494     if(pTrack && p->TrackerTrk() &&
495     pTrack == p->TrackerTrk()) continue;
496    
497     if(pTrack->Pt() <= ptMin) continue;
498    
499     Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
500     if ( dr < extRadius && dr >= intRadius ) {
501     Pt_jets_X_tot += pTrack->Px();
502     Pt_jets_Y_tot += pTrack->Py();
503 ceballos 1.8 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
504 ceballos 1.7 if(pDz < delta_z){
505     Pt_jets_X += pTrack->Px();
506     Pt_jets_Y += pTrack->Py();
507     }
508     }
509     }
510    
511     if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
512     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);
513    
514     return 1.0;
515     }
516    
517     //--------------------------------------------------------------------------------------------------
518 ceballos 1.8 Double_t IsolationTools::BetaE(const TrackCol *tracks, const Electron *p, const Vertex *vertex,
519 ceballos 1.7 Double_t ptMin, Double_t delta_z, Double_t extRadius,
520     Double_t intRadius){
521    
522     if(!tracks) return 1.0;
523     if(tracks->GetEntries() <= 0) return 1.0;
524    
525     double Pt_jets_X = 0. ;
526     double Pt_jets_Y = 0. ;
527     double Pt_jets_X_tot = 0. ;
528     double Pt_jets_Y_tot = 0. ;
529    
530     for(int i=0;i<int(tracks->GetEntries());i++){
531     const Track *pTrack = tracks->At(i);
532    
533     if(pTrack && p->TrackerTrk() &&
534     pTrack == p->TrackerTrk()) continue;
535    
536     if(pTrack && p->GsfTrk() &&
537     pTrack == p->GsfTrk()) continue;
538    
539     if(pTrack->Pt() <= ptMin) continue;
540    
541     Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
542     if ( dr < extRadius && dr >= intRadius ) {
543     Pt_jets_X_tot += pTrack->Px();
544     Pt_jets_Y_tot += pTrack->Py();
545 ceballos 1.8 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
546 ceballos 1.7 if(pDz < delta_z){
547     Pt_jets_X += pTrack->Px();
548     Pt_jets_Y += pTrack->Py();
549     }
550     }
551     }
552    
553     if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
554     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);
555    
556     return 1.0;
557     }
558 fabstoec 1.10
559    
560     // method added by F.Stoeckli: computes the track isolation with NO constrint on the OV-track compatibility
561     Double_t IsolationTools::TrackIsolationNoPV(const mithep::Particle* p, const BaseVertex* bsp,
562     Double_t extRadius,
563     Double_t intRadius,
564     Double_t ptLow,
565     Double_t etaStrip,
566     Double_t maxD0,
567     mithep::TrackQuality::EQuality quality,
568 bendavid 1.11 const mithep::Collection<mithep::Track> *tracks,
569     UInt_t maxNExpectedHitsInner,
570     const mithep::DecayParticleCol *conversions) {
571 fabstoec 1.10
572     // loop over all tracks
573     Double_t tPt = 0.;
574 fabstoec 1.22 //std::cout<<" *** TrackIso:"<<std::endl;
575 fabstoec 1.10 for(UInt_t i=0; i<tracks->GetEntries(); ++i) {
576     const Track* t = tracks->At(i);
577 fabstoec 1.22 if(t->Pt()>1. && false) {
578     std::cout<<" "<<i<<" pt = "<<t->Pt()<<" ("<<ptLow<<")"<<std::endl;
579     std::cout<<" d0 = "<<fabs(t->D0Corrected( *bsp) )<<" ("<<maxD0<<")"<<std::endl;
580     //std::cout<<" conv ? "<<PhotonTools::MatchedConversion(t,conversions,bsp)<<std::endl;
581     std::cout<<" dR = "<<MathUtils::DeltaR(t->Mom(),p->Mom())<<" ("<<extRadius<<","<<intRadius<<")"<<std::endl;
582     std::cout<<" dEta = "<<fabs(t->Eta()-p->Eta())<<" ("<<etaStrip<<")"<<std::endl;
583     }
584 fabstoec 1.10 if ( t->Pt() < ptLow ) continue;
585     if ( ! t->Quality().Quality(quality) ) continue;
586     // only check for beamspot if available, otherwise ignore cut
587     if ( bsp && fabs(t->D0Corrected( *bsp) ) > maxD0) continue;
588 bendavid 1.11 if (t->NExpectedHitsInner()>maxNExpectedHitsInner) continue;
589     if (conversions && PhotonTools::MatchedConversion(t,conversions,bsp)) continue;
590 fabstoec 1.10 Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
591     Double_t dEta = fabs(t->Eta()-p->Eta());
592     if(dR < extRadius && dR > intRadius && dEta > etaStrip) tPt += t->Pt();
593     }
594     return tPt;
595     }
596    
597 fabstoec 1.19
598 fabstoec 1.20 Double_t IsolationTools::CiCTrackIsolation(const mithep::Photon* p,
599 fabstoec 1.19 const BaseVertex* theVtx,
600     Double_t extRadius,
601     Double_t intRadius,
602     Double_t ptLow,
603     Double_t etaStrip,
604     Double_t maxD0,
605     Double_t maxDZ,
606     const mithep::Collection<mithep::Track> *tracks,
607 fabstoec 1.20 unsigned int* worstVtxIndex,
608     const mithep::Collection<mithep::Vertex> *vtxs,
609 fabstoec 1.22 const mithep::Collection<mithep::Electron> *eles,
610     bool print,
611     double* ptmax, double* dRmax) {
612 fabstoec 1.19
613     UInt_t numVtx = 1;
614     const BaseVertex* iVtx = theVtx;
615     if( vtxs ) {
616     numVtx = vtxs->GetEntries();
617     if (numVtx > 0)
618     iVtx = vtxs->At(0);
619     else
620     return 0.;
621     }
622    
623 fabstoec 1.22 // NEW for Electron T&P: need to remove the electron Gsf Track (applied if eles != NULL)
624     const Track* theGsfTrack = NULL;
625     if ( eles ) {
626     // find the electron that matches the Photon SC
627     for(UInt_t j=0; j<eles->GetEntries(); ++j) {
628     if ( eles->At(j)->SCluster() == p->SCluster() ) {
629     if( eles->At(j)->HasTrackerTrk() )
630     theGsfTrack = eles->At(j)->TrackerTrk();
631     break;
632     }
633     }
634     }
635 fabstoec 1.20
636     if(print) {
637     std::cout<<" Testing photon with"<<std::endl;
638     std::cout<<" Et = "<<p->Et()<<std::endl;
639     std::cout<<" Eta = "<<p->Eta()<<std::endl;
640     std::cout<<" Phi = "<<p->Phi()<<std::endl;
641     }
642    
643 fabstoec 1.19 Double_t iIso = 0.;
644     Double_t maxIso = 0.;
645 fabstoec 1.20
646     if(worstVtxIndex)
647     *worstVtxIndex=0;
648 fabstoec 1.19
649 fabstoec 1.22 double t_ptmax = 0.;
650     double t_dRmax = 0.;
651    
652 fabstoec 1.19 for(UInt_t i=0; i<numVtx; ++i) {
653 fabstoec 1.20
654     if(i>0) iVtx = vtxs->At(i);
655    
656    
657     if(print) {
658     std::cout<<" Vertex #"<<i<<std::endl;
659     std::cout<<" with X = "<<iVtx->X()<<std::endl;
660     std::cout<<" with Y = "<<iVtx->Y()<<std::endl;
661     std::cout<<" with Z = "<<iVtx->Z()<<std::endl;
662     }
663    
664 fabstoec 1.22 Photon* phTemp = new Photon(*p);
665    
666     // RESET CALO_POS!
667     phTemp->SetCaloPosXYZ(p->SCluster()->Point().X(),p->SCluster()->Point().Y(),p->SCluster()->Point().Z());
668    
669 fabstoec 1.20 // compute the ph momentum with respect to this Vtx
670 fabstoec 1.22 //FourVectorM phMom = p->MomVtx(iVtx->Position());
671     FourVectorM phMom = phTemp->MomVtx(iVtx->Position());
672    
673     delete phTemp;
674 fabstoec 1.20
675     if(print) {
676     std::cout<<" photon has changed to:"<<std::endl;
677     std::cout<<" Et = "<<phMom.Et()<<std::endl;
678     std::cout<<" eta = "<<phMom.Eta()<<std::endl;
679     std::cout<<" Phi = "<<phMom.Phi()<<std::endl;
680     }
681    
682 fabstoec 1.19 iIso = 0.;
683 fabstoec 1.20 for(UInt_t j=0; j<tracks->GetEntries(); ++j) {
684     const Track* t = tracks->At(j);
685 fabstoec 1.22 if( theGsfTrack && t == theGsfTrack ) continue;
686 fabstoec 1.20
687     //Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
688     //Double_t dEta = TMath::Abs(t->Eta()-p->Eta());
689    
690     Double_t dR = MathUtils::DeltaR(t->Mom(),phMom);
691     Double_t dEta = TMath::Abs(t->Eta()-phMom.Eta());
692    
693 fabstoec 1.22 if(print && t->Pt()>1. && false) {
694 fabstoec 1.20 std::cout<<" passing track #"<<j<<std::endl;
695     std::cout<<" pt = "<<t->Pt()<<std::endl;
696     std::cout<<" eta = "<<t->Eta()<<std::endl;
697     std::cout<<" phi = "<<t->Phi()<<std::endl;
698     std::cout<<" d0 = "<<fabs(t->D0Corrected( *iVtx ))<<std::endl;
699     std::cout<<" dZ = "<<fabs(t->DzCorrected( *iVtx ))<<std::endl;
700     std::cout<<" dR = "<<dR<<std::endl;
701     std::cout<<" dEta = "<<dEta<<std::endl;
702     std::cout<<" vx = "<<t->X0()<<std::endl;
703     std::cout<<" vy = "<<t->Y0()<<std::endl;
704     std::cout<<" vz = "<<t->Z0()<<std::endl;
705     }
706    
707 fabstoec 1.19 if ( t->Pt() < ptLow ) continue;
708     // only check for beamspot if available, otherwise ignore cut
709     if ( fabs(t->D0Corrected( *iVtx )) > maxD0) continue;
710     if ( fabs(t->DzCorrected( *iVtx )) > maxDZ) continue;
711 fabstoec 1.20
712    
713     if(dR < extRadius && dR > intRadius && dEta > etaStrip) {
714     iIso += t->Pt();
715    
716 fabstoec 1.22 if(t->Pt() > t_ptmax) {
717     t_ptmax=t->Pt();
718     t_dRmax=dR;
719     }
720    
721     if(print && t->Pt()>1.) {
722 fabstoec 1.20 std::cout<<" passing track #"<<j<<std::endl;
723     std::cout<<" pt = "<<t->Pt()<<std::endl;
724     std::cout<<" eta = "<<t->Eta()<<std::endl;
725     std::cout<<" phi = "<<t->Phi()<<std::endl;
726     std::cout<<" d0 = "<<fabs(t->D0Corrected( *iVtx ))<<std::endl;
727     std::cout<<" dZ = "<<fabs(t->DzCorrected( *iVtx ))<<std::endl;
728     std::cout<<" dR = "<<dR<<std::endl;
729     std::cout<<" dEta = "<<dEta<<std::endl;
730     std::cout<<" vx = "<<t->X0()<<std::endl;
731     std::cout<<" vy = "<<t->Y0()<<std::endl;
732     std::cout<<" vz = "<<t->Z0()<<std::endl;
733     std::cout<<" new tIso = "<<iIso<<std::endl;
734     }
735     }
736     }
737     if ( iIso > maxIso ) {
738     maxIso = iIso;
739     if(worstVtxIndex)
740     *worstVtxIndex=i;
741 fabstoec 1.19 }
742     }
743 fabstoec 1.20
744 fabstoec 1.22 if(ptmax) (*ptmax)=t_ptmax;
745     if(dRmax) (*dRmax)=t_dRmax;
746    
747 fabstoec 1.20 if(print) {
748     if(worstVtxIndex)
749     std::cout<<" max TrkIso is given by Vtx #"<<*worstVtxIndex<<" with an amount of tIso = "<<maxIso<<std::endl;
750     else
751     std::cout<<" max TrkIso is given by Vtx #0 with an amount of tIso = "<<maxIso<<std::endl;
752     }
753     return maxIso;
754 fabstoec 1.19 }