ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/IsolationTools.cc
Revision: 1.11
Committed: Thu Apr 14 22:05:42 2011 UTC (14 years ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_020d, TMit_020d, Mit_020c, Mit_021pre1, Mit_020b
Changes since 1.10: +7 -2 lines
Log Message:
add conversion veto option for photon isolation

File Contents

# User Rev Content
1 bendavid 1.11 // $Id: IsolationTools.cc,v 1.10 2011/04/06 18:03:48 fabstoec Exp $
2 loizides 1.1
3     #include "MitPhysics/Utils/interface/IsolationTools.h"
4 bendavid 1.11 #include "MitPhysics/Utils/interface/PhotonTools.h"
5 loizides 1.1 #include "MitCommon/MathTools/interface/MathUtils.h"
6    
7 loizides 1.4 ClassImp(mithep::IsolationTools)
8    
9 loizides 1.1 using namespace mithep;
10    
11     //--------------------------------------------------------------------------------------------------
12     Double_t IsolationTools::TrackIsolation(const Track *p, Double_t extRadius, Double_t intRadius,
13     Double_t ptLow, Double_t maxVtxZDist,
14 phedex 1.2 const Collection<Track> *tracks)
15 loizides 1.1 {
16     //Computes the Track Isolation: Summed Transverse Momentum of all tracks inside an
17     //annulus around the electron seed track.
18    
19     Double_t ptSum =0.;
20     for (UInt_t i=0; i<tracks->GetEntries();i++) {
21     Double_t tmpPt = tracks->At(i)->Pt();
22     Double_t deltaZ = fabs(p->Z0() - tracks->At(i)->Z0());
23    
24     //ignore the track if it is below the pt threshold
25     if (tmpPt < ptLow)
26     continue;
27     //ingore the track if it is too far away in Z
28     if (deltaZ > maxVtxZDist)
29     continue;
30    
31 phedex 1.2 Double_t dr = MathUtils::DeltaR(p->Phi(),p->Eta(),tracks->At(i)->Phi(), tracks->At(i)->Eta());
32 loizides 1.1 //add the track pt if it is inside the annulus
33     if ( dr < extRadius &&
34     dr >= intRadius ) {
35     ptSum += tmpPt;
36     }
37     }
38     return ptSum;
39     }
40    
41     //--------------------------------------------------------------------------------------------------
42     Double_t IsolationTools::EcalIsolation(const SuperCluster *sc, Double_t coneSize, Double_t etLow,
43 phedex 1.2 const Collection<BasicCluster> *basicClusters)
44 loizides 1.1 {
45     //Computes the Ecal Isolation: Summed Transverse Energy of all Basic Clusters inside a
46     //cone around the electron, excluding those that are inside the electron super cluster.
47    
48     Double_t ecalIsol=0.;
49     const BasicCluster *basicCluster= 0;
50     for (UInt_t i=0; i<basicClusters->GetEntries();i++) {
51     basicCluster = basicClusters->At(i);
52     Double_t basicClusterEnergy = basicCluster->Energy();
53     Double_t basicClusterEta = basicCluster->Eta();
54     Double_t basicClusterEt = basicClusterEnergy*sin(2*atan(exp(basicClusterEta)));
55    
56 bendavid 1.3 if (basicClusterEt > etLow) {
57 loizides 1.1 bool inSuperCluster = false;
58    
59     // loop over the basic clusters of the supercluster
60     // to make sure that the basic cluster is not inside
61     // the super cluster. We exclude those.
62     for (UInt_t j=0; j<sc->ClusterSize(); j++) {
63     const BasicCluster *tempBasicClusterInSuperCluster = sc->Cluster(j);
64     if (tempBasicClusterInSuperCluster == basicCluster) {
65     inSuperCluster = true;
66     }
67     }
68    
69     if (!inSuperCluster) {
70 phedex 1.2 Double_t dr = MathUtils::DeltaR(sc->Phi(), sc->Eta(),
71     basicCluster->Phi(),basicCluster->Eta());
72 loizides 1.1 if(dr < coneSize) {
73     ecalIsol += basicClusterEt;
74     }
75     }
76     }
77     }
78     return ecalIsol;
79     }
80    
81     //--------------------------------------------------------------------------------------------------
82     Double_t IsolationTools::CaloTowerHadIsolation(const ThreeVector *p, Double_t extRadius,
83     Double_t intRadius, Double_t etLow,
84 phedex 1.2 const Collection<CaloTower> *caloTowers)
85 loizides 1.1 {
86     //Computes the CaloTower Had Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
87     //inside an annulus around the electron super cluster position.
88    
89     Double_t sumEt = 0;
90     for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
91     Double_t caloTowerEt = caloTowers->At(i)->HadEt();
92     Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
93     p->Phi(), p->Eta());
94     if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
95     sumEt += caloTowerEt;
96     }
97     }
98     return sumEt;
99     }
100    
101     //--------------------------------------------------------------------------------------------------
102     Double_t IsolationTools::CaloTowerEmIsolation(const ThreeVector *p, Double_t extRadius,
103     Double_t intRadius, Double_t etLow,
104 phedex 1.2 const Collection<CaloTower> *caloTowers)
105 loizides 1.1 {
106     //Computes the CaloTower Em Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
107     //inside an annulus around the electron super cluster position.
108    
109     Double_t sumEt = 0;
110     for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
111     Double_t caloTowerEt = caloTowers->At(i)->EmEt();
112     Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
113     p->Phi(), p->Eta());
114     if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
115     sumEt += caloTowerEt;
116     }
117     }
118     return sumEt;
119     }
120 ceballos 1.5
121     //--------------------------------------------------------------------------------------------------
122     Double_t IsolationTools::PFMuonIsolation(const Muon *p, const Collection<PFCandidate> *PFCands,
123 ceballos 1.8 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
124     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 ceballos 1.6 if(pf->Pt() <= ptMin) continue;
152    
153 ceballos 1.5 if(pf->TrackerTrk() && p->TrackerTrk() &&
154     pf->TrackerTrk() == p->TrackerTrk()) continue;
155    
156     Double_t deltaZ = 0.0;
157     if(pf->BestTrk()) {
158 ceballos 1.8 deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
159 ceballos 1.5 }
160    
161     // ignore the pf candidate if it is too far away in Z
162 ceballos 1.6 if (deltaZ >= delta_z)
163 ceballos 1.5 continue;
164    
165     Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
166     // add the pf pt if it is inside the extRadius and outside the intRadius
167     if ( dr < extRadius &&
168     dr >= intRadius ) {
169 ceballos 1.8 Bool_t isLepton = kFALSE;
170     if(goodMuons && isoType == 3){
171     for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
172     const Muon *m = goodMuons->At(nl);
173     if(pf->TrackerTrk() && m->TrackerTrk() &&
174     pf->TrackerTrk() == m->TrackerTrk()) {
175     isLepton = kTRUE;
176     break;
177     }
178     }
179     }
180     if(goodElectrons && isLepton == kFALSE && isoType == 3){
181     for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
182     const Electron *e = goodElectrons->At(nl);
183     if(pf->TrackerTrk() && e->TrackerTrk() &&
184     pf->TrackerTrk() == e->TrackerTrk()) {
185     isLepton = kTRUE;
186     break;
187     }
188     if(pf->GsfTrk() && e->GsfTrk() &&
189     pf->GsfTrk() == e->GsfTrk()) {
190     isLepton = kTRUE;
191     break;
192     }
193     }
194     }
195     if(isLepton == kFALSE){
196     if(pf->BestTrk()) ptSum += pf->Pt();
197     else ptSum += pf->Pt()*beta;
198     }
199 ceballos 1.5 }
200     }
201     return ptSum;
202     }
203     //--------------------------------------------------------------------------------------------------
204     Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands,
205 ceballos 1.8 const Vertex *vertex, Double_t delta_z, Double_t ptMin,
206     Double_t extRadius, Double_t intRadius, int isoType,
207     Double_t beta, const MuonCol *goodMuons,
208     const ElectronCol *goodElectrons)
209 ceballos 1.5 {
210     //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
211     //annulus around the particle seed track.
212    
213     Double_t zLepton = 0.0;
214 ceballos 1.8 if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
215 ceballos 1.5
216     Double_t ptSum =0.;
217     for (UInt_t i=0; i<PFCands->GetEntries();i++) {
218     const PFCandidate *pf = PFCands->At(i);
219    
220     Bool_t isGoodType = kFALSE;
221     // all particles
222 ceballos 1.6 if (isoType == 0) isGoodType = kTRUE;
223 ceballos 1.5 // charged particles only
224 ceballos 1.6 else if(isoType == 1 && pf->BestTrk()) isGoodType = kTRUE;
225 ceballos 1.5 // charged particles and gammas only
226 ceballos 1.6 else if(isoType == 2 &&
227     (pf->BestTrk() || pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
228 ceballos 1.8 // all particles, rejecting good leptons
229     else if(isoType == 3) isGoodType = kTRUE;
230 ceballos 1.5
231     if(isGoodType == kFALSE) continue;
232    
233 ceballos 1.6 if(pf->Pt() <= ptMin) continue;
234    
235 ceballos 1.5 if(pf->TrackerTrk() && p->TrackerTrk() &&
236     pf->TrackerTrk() == p->TrackerTrk()) continue;
237    
238     if(pf->GsfTrk() && p->GsfTrk() &&
239     pf->GsfTrk() == p->GsfTrk()) continue;
240    
241     Double_t deltaZ = 0.0;
242     if(pf->BestTrk()) {
243 ceballos 1.8 deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
244 ceballos 1.5 }
245    
246     // ignore the pf candidate if it is too far away in Z
247 ceballos 1.6 if (deltaZ >= delta_z)
248 ceballos 1.5 continue;
249    
250     Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
251     // add the pf pt if it is inside the extRadius and outside the intRadius
252     if ( dr < extRadius &&
253     dr >= intRadius ) {
254 ceballos 1.8 Bool_t isLepton = kFALSE;
255     if(goodMuons && isoType == 3){
256     for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
257     const Muon *m = goodMuons->At(nl);
258     if(pf->TrackerTrk() && m->TrackerTrk() &&
259     pf->TrackerTrk() == m->TrackerTrk()) {
260     isLepton = kTRUE;
261     break;
262     }
263     }
264     }
265     if(goodElectrons && isLepton == kFALSE && isoType == 3){
266     for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
267     const Electron *e = goodElectrons->At(nl);
268     if(pf->TrackerTrk() && e->TrackerTrk() &&
269     pf->TrackerTrk() == e->TrackerTrk()) {
270     isLepton = kTRUE;
271     break;
272     }
273     if(pf->GsfTrk() && e->GsfTrk() &&
274     pf->GsfTrk() == e->GsfTrk()) {
275     isLepton = kTRUE;
276     break;
277     }
278     }
279     }
280     if(isLepton == kFALSE){
281     if(pf->BestTrk()) ptSum += pf->Pt();
282     else ptSum += pf->Pt()*beta;
283     }
284 ceballos 1.5 }
285     }
286     return ptSum;
287     }
288 ceballos 1.7 //--------------------------------------------------------------------------------------------------
289 ceballos 1.8 Double_t IsolationTools::BetaM(const TrackCol *tracks, const Muon *p, const Vertex *vertex,
290 ceballos 1.7 Double_t ptMin, Double_t delta_z, Double_t extRadius,
291     Double_t intRadius){
292    
293 ceballos 1.9 if(!tracks) return 1.0;
294 ceballos 1.7 if(tracks->GetEntries() <= 0) return 1.0;
295    
296     double Pt_jets_X = 0. ;
297     double Pt_jets_Y = 0. ;
298     double Pt_jets_X_tot = 0. ;
299     double Pt_jets_Y_tot = 0. ;
300    
301     for(int i=0;i<int(tracks->GetEntries());i++){
302     const Track *pTrack = tracks->At(i);
303    
304     if(pTrack && p->TrackerTrk() &&
305     pTrack == p->TrackerTrk()) continue;
306    
307     if(pTrack->Pt() <= ptMin) continue;
308    
309     Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
310     if ( dr < extRadius && dr >= intRadius ) {
311     Pt_jets_X_tot += pTrack->Px();
312     Pt_jets_Y_tot += pTrack->Py();
313 ceballos 1.8 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
314 ceballos 1.7 if(pDz < delta_z){
315     Pt_jets_X += pTrack->Px();
316     Pt_jets_Y += pTrack->Py();
317     }
318     }
319     }
320    
321     if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
322     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);
323    
324     return 1.0;
325     }
326    
327     //--------------------------------------------------------------------------------------------------
328 ceballos 1.8 Double_t IsolationTools::BetaE(const TrackCol *tracks, const Electron *p, const Vertex *vertex,
329 ceballos 1.7 Double_t ptMin, Double_t delta_z, Double_t extRadius,
330     Double_t intRadius){
331    
332     if(!tracks) return 1.0;
333     if(tracks->GetEntries() <= 0) return 1.0;
334    
335     double Pt_jets_X = 0. ;
336     double Pt_jets_Y = 0. ;
337     double Pt_jets_X_tot = 0. ;
338     double Pt_jets_Y_tot = 0. ;
339    
340     for(int i=0;i<int(tracks->GetEntries());i++){
341     const Track *pTrack = tracks->At(i);
342    
343     if(pTrack && p->TrackerTrk() &&
344     pTrack == p->TrackerTrk()) continue;
345    
346     if(pTrack && p->GsfTrk() &&
347     pTrack == p->GsfTrk()) continue;
348    
349     if(pTrack->Pt() <= ptMin) continue;
350    
351     Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
352     if ( dr < extRadius && dr >= intRadius ) {
353     Pt_jets_X_tot += pTrack->Px();
354     Pt_jets_Y_tot += pTrack->Py();
355 ceballos 1.8 double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
356 ceballos 1.7 if(pDz < delta_z){
357     Pt_jets_X += pTrack->Px();
358     Pt_jets_Y += pTrack->Py();
359     }
360     }
361     }
362    
363     if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
364     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);
365    
366     return 1.0;
367     }
368 fabstoec 1.10
369    
370     // method added by F.Stoeckli: computes the track isolation with NO constrint on the OV-track compatibility
371     Double_t IsolationTools::TrackIsolationNoPV(const mithep::Particle* p, const BaseVertex* bsp,
372     Double_t extRadius,
373     Double_t intRadius,
374     Double_t ptLow,
375     Double_t etaStrip,
376     Double_t maxD0,
377     mithep::TrackQuality::EQuality quality,
378 bendavid 1.11 const mithep::Collection<mithep::Track> *tracks,
379     UInt_t maxNExpectedHitsInner,
380     const mithep::DecayParticleCol *conversions) {
381 fabstoec 1.10
382     // loop over all tracks
383     Double_t tPt = 0.;
384     for(UInt_t i=0; i<tracks->GetEntries(); ++i) {
385     const Track* t = tracks->At(i);
386     if ( t->Pt() < ptLow ) continue;
387     if ( ! t->Quality().Quality(quality) ) continue;
388     // only check for beamspot if available, otherwise ignore cut
389     if ( bsp && fabs(t->D0Corrected( *bsp) ) > maxD0) continue;
390 bendavid 1.11 if (t->NExpectedHitsInner()>maxNExpectedHitsInner) continue;
391     if (conversions && PhotonTools::MatchedConversion(t,conversions,bsp)) continue;
392 fabstoec 1.10 Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
393     Double_t dEta = fabs(t->Eta()-p->Eta());
394     if(dR < extRadius && dR > intRadius && dEta > etaStrip) tPt += t->Pt();
395     }
396     return tPt;
397     }
398