ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/IsolationTools.cc
Revision: 1.14
Committed: Tue May 10 13:04:15 2011 UTC (13 years, 11 months ago) by sixie
Content type: text/plain
Branch: MAIN
Changes since 1.13: +20 -6 lines
Log Message:
update PFIso for electrons, include strip veto for gammas, one tower veto for neutral hadrons

File Contents

# User Rev Content
1 sixie 1.14 // $Id: IsolationTools.cc,v 1.13 2011/05/10 10:39:14 mzanetti 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     Double_t extRadius, Double_t intRadius, int isoType,
125     Double_t beta, const MuonCol *goodMuons,
126     const ElectronCol *goodElectrons)
127 ceballos 1.5 {
128     //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
129     //annulus around the particle seed track.
130    
131     Double_t zLepton = 0.0;
132 ceballos 1.8 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
133 ceballos 1.5
134     Double_t ptSum =0.;
135     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
136     const PFCandidate *pf = PFCands->At(i);
137    
138     Bool_t isGoodType = kFALSE;
139     // all particles
140 ceballos 1.6 if (isoType == 0) isGoodType = kTRUE;
141 ceballos 1.5 // charged particles only
142 ceballos 1.6 else if(isoType == 1 && pf->BestTrk()) isGoodType = kTRUE;
143 ceballos 1.5 // charged particles and gammas only
144 ceballos 1.6 else if(isoType == 2 &&
145     (pf->BestTrk() || pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
146 ceballos 1.8 // all particles, rejecting good leptons
147     else if(isoType == 3) isGoodType = kTRUE;
148 ceballos 1.5
149     if(isGoodType == kFALSE) continue;
150    
151 mzanetti 1.13 // 0.1 pt cut applied to charged
152 ceballos 1.12 if( pf->BestTrk() && pf->Pt() <= 0.1) continue;
153 mzanetti 1.13 // pt cut applied to neutrals
154     if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
155 ceballos 1.6
156 ceballos 1.5 if(pf->TrackerTrk() && p->TrackerTrk() &&
157     pf->TrackerTrk() == p->TrackerTrk()) continue;
158    
159     Double_t deltaZ = 0.0;
160     if(pf->BestTrk()) {
161 ceballos 1.8 deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
162 ceballos 1.5 }
163    
164     // ignore the pf candidate if it is too far away in Z
165 ceballos 1.6 if (deltaZ >= delta_z)
166 ceballos 1.5 continue;
167    
168     Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
169     // add the pf pt if it is inside the extRadius and outside the intRadius
170     if ( dr < extRadius &&
171     dr >= intRadius ) {
172 ceballos 1.8 Bool_t isLepton = kFALSE;
173     if(goodMuons && isoType == 3){
174     for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
175     const Muon *m = goodMuons->At(nl);
176     if(pf->TrackerTrk() && m->TrackerTrk() &&
177     pf->TrackerTrk() == m->TrackerTrk()) {
178     isLepton = kTRUE;
179     break;
180     }
181     }
182     }
183     if(goodElectrons && isLepton == kFALSE && isoType == 3){
184     for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
185     const Electron *e = goodElectrons->At(nl);
186     if(pf->TrackerTrk() && e->TrackerTrk() &&
187     pf->TrackerTrk() == e->TrackerTrk()) {
188     isLepton = kTRUE;
189     break;
190     }
191     if(pf->GsfTrk() && e->GsfTrk() &&
192     pf->GsfTrk() == e->GsfTrk()) {
193     isLepton = kTRUE;
194     break;
195     }
196     }
197     }
198     if(isLepton == kFALSE){
199     if(pf->BestTrk()) ptSum += pf->Pt();
200     else ptSum += pf->Pt()*beta;
201     }
202 ceballos 1.5 }
203     }
204     return ptSum;
205     }
206     //--------------------------------------------------------------------------------------------------
207     Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands,
208 ceballos 1.8 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
209     Double_t extRadius, Double_t intRadius, int isoType,
210     Double_t beta, const MuonCol *goodMuons,
211     const ElectronCol *goodElectrons)
212 ceballos 1.5 {
213     //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
214     //annulus around the particle seed track.
215    
216     Double_t zLepton = 0.0;
217 ceballos 1.8 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
218 ceballos 1.5
219     Double_t ptSum =0.;
220     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
221     const PFCandidate *pf = PFCands->At(i);
222    
223     Bool_t isGoodType = kFALSE;
224     // all particles
225 ceballos 1.6 if (isoType == 0) isGoodType = kTRUE;
226 ceballos 1.5 // charged particles only
227 ceballos 1.6 else if(isoType == 1 && pf->BestTrk()) isGoodType = kTRUE;
228 ceballos 1.5 // charged particles and gammas only
229 ceballos 1.6 else if(isoType == 2 &&
230     (pf->BestTrk() || pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
231 ceballos 1.8 // all particles, rejecting good leptons
232     else if(isoType == 3) isGoodType = kTRUE;
233 ceballos 1.5
234     if(isGoodType == kFALSE) continue;
235    
236 sixie 1.14
237     // 0.1 pt cut applied to charged
238 ceballos 1.12 if( pf->BestTrk() && pf->Pt() <= 0.1) continue;
239    
240 sixie 1.14 // pt cut applied to neutrals
241     if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
242 ceballos 1.6
243 ceballos 1.5 if(pf->TrackerTrk() && p->TrackerTrk() &&
244     pf->TrackerTrk() == p->TrackerTrk()) continue;
245    
246     if(pf->GsfTrk() && p->GsfTrk() &&
247     pf->GsfTrk() == p->GsfTrk()) continue;
248    
249     Double_t deltaZ = 0.0;
250     if(pf->BestTrk()) {
251 ceballos 1.8 deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
252 ceballos 1.5 }
253    
254     // ignore the pf candidate if it is too far away in Z
255 ceballos 1.6 if (deltaZ >= delta_z)
256 ceballos 1.5 continue;
257    
258     Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
259     // add the pf pt if it is inside the extRadius and outside the intRadius
260     if ( dr < extRadius &&
261     dr >= intRadius ) {
262 ceballos 1.8 Bool_t isLepton = kFALSE;
263     if(goodMuons && isoType == 3){
264     for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
265     const Muon *m = goodMuons->At(nl);
266     if(pf->TrackerTrk() && m->TrackerTrk() &&
267     pf->TrackerTrk() == m->TrackerTrk()) {
268     isLepton = kTRUE;
269     break;
270     }
271     }
272     }
273     if(goodElectrons && isLepton == kFALSE && isoType == 3){
274     for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
275     const Electron *e = goodElectrons->At(nl);
276     if(pf->TrackerTrk() && e->TrackerTrk() &&
277     pf->TrackerTrk() == e->TrackerTrk()) {
278     isLepton = kTRUE;
279     break;
280     }
281     if(pf->GsfTrk() && e->GsfTrk() &&
282     pf->GsfTrk() == e->GsfTrk()) {
283     isLepton = kTRUE;
284     break;
285     }
286     }
287     }
288 sixie 1.14
289     if (isLepton == kTRUE) continue;
290    
291     //EtaStrip Veto for Gamma
292     if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;
293    
294     //InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
295     if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
296     && MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue;
297    
298    
299     if(pf->BestTrk()) ptSum += pf->Pt();
300     else ptSum += pf->Pt()*beta;
301    
302    
303 ceballos 1.5 }
304     }
305     return ptSum;
306     }
307 ceballos 1.7 //--------------------------------------------------------------------------------------------------
308 ceballos 1.8 Double_t IsolationTools::BetaM(const TrackCol *tracks, const Muon *p, const Vertex *vertex,
309 ceballos 1.7 Double_t ptMin, Double_t delta_z, Double_t extRadius,
310     Double_t intRadius){
311    
312 ceballos 1.9 if(!tracks) return 1.0;
313 ceballos 1.7 if(tracks->GetEntries() <= 0) return 1.0;
314    
315     double Pt_jets_X = 0. ;
316     double Pt_jets_Y = 0. ;
317     double Pt_jets_X_tot = 0. ;
318     double Pt_jets_Y_tot = 0. ;
319    
320     for(int i=0;i<int(tracks->GetEntries());i++){
321     const Track *pTrack = tracks->At(i);
322    
323     if(pTrack && p->TrackerTrk() &&
324     pTrack == p->TrackerTrk()) continue;
325    
326     if(pTrack->Pt() <= ptMin) continue;
327    
328     Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
329     if ( dr < extRadius && dr >= intRadius ) {
330     Pt_jets_X_tot += pTrack->Px();
331     Pt_jets_Y_tot += pTrack->Py();
332 ceballos 1.8 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
333 ceballos 1.7 if(pDz < delta_z){
334     Pt_jets_X += pTrack->Px();
335     Pt_jets_Y += pTrack->Py();
336     }
337     }
338     }
339    
340     if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
341     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);
342    
343     return 1.0;
344     }
345    
346     //--------------------------------------------------------------------------------------------------
347 ceballos 1.8 Double_t IsolationTools::BetaE(const TrackCol *tracks, const Electron *p, const Vertex *vertex,
348 ceballos 1.7 Double_t ptMin, Double_t delta_z, Double_t extRadius,
349     Double_t intRadius){
350    
351     if(!tracks) return 1.0;
352     if(tracks->GetEntries() <= 0) return 1.0;
353    
354     double Pt_jets_X = 0. ;
355     double Pt_jets_Y = 0. ;
356     double Pt_jets_X_tot = 0. ;
357     double Pt_jets_Y_tot = 0. ;
358    
359     for(int i=0;i<int(tracks->GetEntries());i++){
360     const Track *pTrack = tracks->At(i);
361    
362     if(pTrack && p->TrackerTrk() &&
363     pTrack == p->TrackerTrk()) continue;
364    
365     if(pTrack && p->GsfTrk() &&
366     pTrack == p->GsfTrk()) continue;
367    
368     if(pTrack->Pt() <= ptMin) continue;
369    
370     Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
371     if ( dr < extRadius && dr >= intRadius ) {
372     Pt_jets_X_tot += pTrack->Px();
373     Pt_jets_Y_tot += pTrack->Py();
374 ceballos 1.8 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
375 ceballos 1.7 if(pDz < delta_z){
376     Pt_jets_X += pTrack->Px();
377     Pt_jets_Y += pTrack->Py();
378     }
379     }
380     }
381    
382     if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
383     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);
384    
385     return 1.0;
386     }
387 fabstoec 1.10
388    
389     // method added by F.Stoeckli: computes the track isolation with NO constrint on the OV-track compatibility
390     Double_t IsolationTools::TrackIsolationNoPV(const mithep::Particle* p, const BaseVertex* bsp,
391     Double_t extRadius,
392     Double_t intRadius,
393     Double_t ptLow,
394     Double_t etaStrip,
395     Double_t maxD0,
396     mithep::TrackQuality::EQuality quality,
397 bendavid 1.11 const mithep::Collection<mithep::Track> *tracks,
398     UInt_t maxNExpectedHitsInner,
399     const mithep::DecayParticleCol *conversions) {
400 fabstoec 1.10
401     // loop over all tracks
402     Double_t tPt = 0.;
403     for(UInt_t i=0; i<tracks->GetEntries(); ++i) {
404     const Track* t = tracks->At(i);
405     if ( t->Pt() < ptLow ) continue;
406     if ( ! t->Quality().Quality(quality) ) continue;
407     // only check for beamspot if available, otherwise ignore cut
408     if ( bsp && fabs(t->D0Corrected( *bsp) ) > maxD0) continue;
409 bendavid 1.11 if (t->NExpectedHitsInner()>maxNExpectedHitsInner) continue;
410     if (conversions && PhotonTools::MatchedConversion(t,conversions,bsp)) continue;
411 fabstoec 1.10 Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
412     Double_t dEta = fabs(t->Eta()-p->Eta());
413     if(dR < extRadius && dR > intRadius && dEta > etaStrip) tPt += t->Pt();
414     }
415     return tPt;
416     }
417