ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/IsolationTools.cc
Revision: 1.29
Committed: Mon May 7 18:05:52 2012 UTC (12 years, 11 months ago) by ceballos
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_028, Mit_027, Mit_027a
Changes since 1.28: +15 -9 lines
Log Message: 
new rhos

File Contents

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