ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/ElectronIDMod.cc
Revision: 1.48
Committed: Sun Nov 8 13:52:29 2009 UTC (15 years, 5 months ago) by sixie
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_012d, Mit_012c, Mit_012b
Changes since 1.47: +3 -1 lines
Log Message:
fix bug with conversion filter. was failing incorrectly if we did not use conversion filter.

File Contents

# User Rev Content
1 sixie 1.48 // $Id: ElectronIDMod.cc,v 1.47 2009/11/02 13:52:55 ceballos Exp $
2 loizides 1.1
3     #include "MitPhysics/Mods/interface/ElectronIDMod.h"
4 loizides 1.26 #include "MitAna/DataTree/interface/StableData.h"
5     #include "MitAna/DataTree/interface/ElectronCol.h"
6     #include "MitAna/DataTree/interface/VertexCol.h"
7     #include "MitAna/DataTree/interface/DecayParticleCol.h"
8 loizides 1.5 #include "MitPhysics/Init/interface/ModNames.h"
9 loizides 1.1
10     using namespace mithep;
11    
12     ClassImp(mithep::ElectronIDMod)
13    
14     //--------------------------------------------------------------------------------------------------
15 loizides 1.5 ElectronIDMod::ElectronIDMod(const char *name, const char *title) :
16 loizides 1.1 BaseMod(name,title),
17 loizides 1.5 fElectronBranchName(Names::gkElectronBrn),
18 ceballos 1.12 fConversionBranchName(Names::gkMvfConversionBrn),
19 loizides 1.5 fGoodElectronsName(ModNames::gkGoodElectronsName),
20 ceballos 1.12 fVertexName(string("PrimaryVertexesBeamSpot").c_str()),
21 ceballos 1.31 fElectronIDType("CustomTight"),
22 ceballos 1.8 fElectronIsoType("TrackJuraSliding"),
23 loizides 1.1 fElectronPtMin(10),
24     fIDLikelihoodCut(0.9),
25     fTrackIsolationCut(5.0),
26     fCaloIsolationCut(5.0),
27 loizides 1.5 fEcalJuraIsoCut(5.0),
28     fHcalIsolationCut(5.0),
29 loizides 1.14 fApplyConvFilter(kTRUE),
30 bendavid 1.38 fWrongHitsRequirement(kTRUE),
31 ceballos 1.18 fApplyD0Cut(kTRUE),
32 ceballos 1.45 fChargeFilter(kTRUE),
33 loizides 1.14 fD0Cut(0.025),
34 loizides 1.25 fReverseIsoCut(kFALSE),
35     fReverseD0Cut(kFALSE),
36 loizides 1.5 fElIdType(kIdUndef),
37 ceballos 1.12 fElIsoType(kIsoUndef),
38 loizides 1.14 fElectrons(0),
39     fConversions(0),
40 loizides 1.25 fVertices(0)
41 loizides 1.1 {
42     // Constructor.
43     }
44    
45     //--------------------------------------------------------------------------------------------------
46 loizides 1.30 Bool_t ElectronIDMod::PassCustomID(const Electron *ele) const
47     {
48 loizides 1.32 // Based on RecoEgamma/ElectronIdentification/src/CutBasedElectronID.cc.
49 loizides 1.30
50     Double_t eOverP = ele->ESuperClusterOverP();
51     Double_t fBrem = ele->FBrem();
52    
53     if ((eOverP < fCuts[5][0]) && (fBrem < fCuts[5][1]))
54     return kFALSE;
55    
56     if (eOverP < fCuts[5][2]*(1-fBrem))
57     return kFALSE;
58    
59     Int_t cat = 2;
60     if ((ele->IsEB() && fBrem<0.06) || (ele->IsEE() && fBrem<0.1))
61     cat=1;
62     else if (eOverP < 1.2 && eOverP > 0.8)
63     cat=0;
64    
65     Double_t eSeedOverPin = ele->ESeedClusterOverPIn();
66     Double_t hOverE = ele->HadronicOverEm();
67     Double_t sigmaee = ele->CoviEtaiEta();
68     Double_t deltaPhiIn = TMath::Abs(ele->DeltaPhiSuperClusterTrackAtVtx());
69     Double_t deltaEtaIn = TMath::Abs(ele->DeltaEtaSuperClusterTrackAtVtx());
70    
71     Int_t eb = 1;
72     if (ele->IsEB())
73     eb = 0;
74    
75     if (hOverE>fCuts[0][cat+4*eb])
76     return kFALSE;
77    
78     if (sigmaee>fCuts[1][cat+4*eb])
79     return kFALSE;
80    
81 ceballos 1.40 if (deltaPhiIn>fCuts[2][cat+4*eb])
82     return kFALSE;
83 loizides 1.30
84     if(deltaEtaIn>fCuts[3][cat+4*eb])
85     return kFALSE;
86    
87     if(eSeedOverPin<fCuts[4][cat+4*eb])
88     return kFALSE;
89    
90     return kTRUE;
91     }
92    
93     //--------------------------------------------------------------------------------------------------
94 sixie 1.42 Bool_t ElectronIDMod::PassIDCut(const Electron *ele, EElIdType idType) const
95     {
96    
97     Bool_t idcut = kFALSE;
98     switch (idType) {
99     case kTight:
100     idcut = ele->PassTightID();
101     break;
102     case kLoose:
103     idcut = ele->PassLooseID();
104     break;
105     case kLikelihood:
106     idcut = (ele->IDLikelihood() > fIDLikelihoodCut);
107     break;
108     case kNoId:
109     idcut = kTRUE;
110     break;
111     case kCustomIdLoose:
112     idcut = ElectronIDMod::PassCustomID(ele);
113     break;
114     case kCustomIdTight:
115     idcut = ElectronIDMod::PassCustomID(ele);
116     break;
117     case kZeeId:
118     if (ele->IsEB()) {
119     idcut = (ele->CoviEtaiEta() < 0.01 && ele->DeltaEtaSuperClusterTrackAtVtx() < 0.0071);
120     } else {
121     idcut = (ele->CoviEtaiEta() < 0.028 && ele->DeltaEtaSuperClusterTrackAtVtx() < 0.0066);
122     }
123     break;
124     default:
125     break;
126     }
127    
128     return idcut;
129     }
130    
131    
132     //--------------------------------------------------------------------------------------------------
133     Bool_t ElectronIDMod::PassIsolationCut(const Electron *ele, EElIsoType isoType) const
134     {
135    
136     Bool_t isocut = kFALSE;
137     switch (isoType) {
138     case kTrackCalo:
139     isocut = (ele->TrackIsolationDr03() < fTrackIsolationCut) &&
140     (ele->CaloIsolation() < fCaloIsolationCut);
141     break;
142     case kTrackJura:
143     isocut = (ele->TrackIsolationDr03() < fTrackIsolationCut) &&
144     (ele->EcalRecHitIsoDr04() < fEcalJuraIsoCut) &&
145     (ele->HcalIsolation() < fHcalIsolationCut);
146     break;
147     case kTrackJuraSliding:
148     {
149     Double_t totalIso = ele->TrackIsolationDr03() + ele->EcalRecHitIsoDr04() - 1.5;
150     if (totalIso < (ele->Pt()-10.0)*4.5/20.0 ||
151     totalIso <= 0)
152     isocut = kTRUE;
153    
154     if (fReverseIsoCut == kTRUE &&
155     isocut == kFALSE && totalIso < 10)
156     isocut = kTRUE;
157     else if(fReverseIsoCut == kTRUE)
158     isocut = kFALSE;
159     }
160     break;
161     case kNoIso:
162     isocut = kTRUE;
163     break;
164     case kZeeIso:
165     if (ele->IsEB()) {
166     isocut = (ele->TrackIsolationDr04() < 7.2 && ele->EcalRecHitIsoDr04() < 5.7 && ele->HcalTowerSumEtDr04() < 8.1);
167     } else {
168     isocut = (ele->TrackIsolationDr04() < 5.1 && ele->EcalRecHitIsoDr04() < 5.0 && ele->HcalTowerSumEtDr04() < 3.4);
169     }
170     break;
171     case kCustomIso:
172     default:
173     break;
174     }
175    
176     return isocut;
177     }
178    
179    
180     //--------------------------------------------------------------------------------------------------
181     Bool_t ElectronIDMod::PassConversionFilter(const Electron *ele, const DecayParticleCol *conversions) const
182     {
183     Bool_t isGoodConversion = kFALSE;
184    
185     for (UInt_t ifc=0; ifc<conversions->GetEntries(); ifc++) {
186    
187     Bool_t ConversionMatchFound = kFALSE;
188     for (UInt_t d=0; d<conversions->At(ifc)->NDaughters(); d++) {
189     const Track *trk = dynamic_cast<const ChargedParticle*>
190     (conversions->At(ifc)->Daughter(d))->Trk();
191     if (ele->GsfTrk() == trk) {
192     ConversionMatchFound = kTRUE;
193     break;
194     }
195     }
196    
197     // if match between the e-track and one of the conversion legs
198     if (ConversionMatchFound == kTRUE){
199     isGoodConversion = (conversions->At(ifc)->Prob() > 1e-6) &&
200     (conversions->At(ifc)->Lxy() > 0) &&
201     (conversions->At(ifc)->Lz() > 0) &&
202     (conversions->At(ifc)->Position().Rho() > 2.0);
203    
204     if (isGoodConversion == kTRUE) {
205     for (UInt_t d=0; d<conversions->At(ifc)->NDaughters(); d++) {
206     const Track *trk = dynamic_cast<const ChargedParticle*>
207     (conversions->At(ifc)->Daughter(d))->Trk();
208    
209     if (trk) {
210     // These requirements are not used for the GSF track
211     if (!(trk->NHits() >= 3 && trk->Prob() > 1e-6) && trk!=ele->GsfTrk())
212     isGoodConversion = kFALSE;
213    
214     const StableData *sd = dynamic_cast<const StableData*>
215     (conversions->At(ifc)->DaughterDat(d));
216     if (fWrongHitsRequirement && sd->NWrongHits() != 0)
217     isGoodConversion = kFALSE;
218    
219     } else {
220     isGoodConversion = kFALSE;
221     }
222     }
223     }
224     }
225    
226     if (isGoodConversion == kTRUE) break;
227    
228     } // loop over all conversions
229    
230 ceballos 1.46 return !isGoodConversion;
231 sixie 1.42 }
232    
233     //--------------------------------------------------------------------------------------------------
234     Bool_t ElectronIDMod::PassD0Cut(const Electron *ele, const VertexCol *vertices) const
235     {
236     Bool_t d0cut = kFALSE;
237     // d0 cut
238     Double_t d0_real = 99999;
239     for(UInt_t i0 = 0; i0 < vertices->GetEntries(); i0++) {
240     Double_t pD0 = ele->GsfTrk()->D0Corrected(*vertices->At(i0));
241     if(TMath::Abs(pD0) < TMath::Abs(d0_real)) d0_real = TMath::Abs(pD0);
242     }
243     if(d0_real < fD0Cut) d0cut = kTRUE;
244    
245     if (fReverseD0Cut == kTRUE &&
246     d0cut == kFALSE && d0_real < 0.05)
247     d0cut = kTRUE;
248     else if(fReverseD0Cut == kTRUE)
249     d0cut = kFALSE;
250    
251     return d0cut;
252     }
253    
254     //--------------------------------------------------------------------------------------------------
255     Bool_t ElectronIDMod::PassChargeFilter(const Electron *ele) const
256     {
257     Bool_t passChargeFilter = kTRUE;
258     if(ele->TrackerTrk() &&
259     ele->TrackerTrk()->Charge() != ele->Charge()) passChargeFilter = kFALSE;
260    
261     return passChargeFilter;
262     }
263    
264     //--------------------------------------------------------------------------------------------------
265 loizides 1.1 void ElectronIDMod::Process()
266     {
267     // Process entries of the tree.
268    
269 loizides 1.23 LoadEventObject(fElectronBranchName, fElectrons);
270 loizides 1.1
271 loizides 1.6 ElectronOArr *GoodElectrons = new ElectronOArr;
272     GoodElectrons->SetName(fGoodElectronsName);
273 loizides 1.1
274 ceballos 1.18 for (UInt_t i=0; i<fElectrons->GetEntries(); ++i) {
275 loizides 1.5 const Electron *e = fElectrons->At(i);
276 loizides 1.1
277 loizides 1.5 if (e->Pt() <= fElectronPtMin)
278     continue;
279 loizides 1.1
280 sixie 1.42 //apply id cut
281     Bool_t idcut = PassIDCut(e, fElIdType);
282 loizides 1.5 if (!idcut)
283     continue;
284    
285 sixie 1.42 //apply Isolation Cut
286     Bool_t isocut = PassIsolationCut(e, fElIsoType);
287     if (!isocut)
288 loizides 1.5 continue;
289    
290 loizides 1.14 // apply conversion filter
291 ceballos 1.47 Bool_t passConvVeto = kFALSE;
292 loizides 1.14 if (fApplyConvFilter) {
293 sixie 1.42 LoadEventObject(fConversionBranchName, fConversions);
294 ceballos 1.47 passConvVeto = PassConversionFilter(e, fConversions);
295 sixie 1.48 } else {
296     passConvVeto = kTRUE;
297 ceballos 1.12 }
298 ceballos 1.47 if (passConvVeto == kFALSE) continue;
299 sixie 1.42
300     // apply d0 cut
301 ceballos 1.15 if (fApplyD0Cut) {
302 sixie 1.42 LoadEventObject(fVertexName, fVertices);
303 ceballos 1.46 Bool_t passD0cut = PassD0Cut(e, fVertices);
304 sixie 1.42 if (!passD0cut)
305 ceballos 1.24 continue;
306 ceballos 1.12 }
307    
308 sixie 1.42 //apply charge filter
309     if(fChargeFilter == kTRUE) {
310     Bool_t passChargeFilter = PassChargeFilter(e);
311     if (!passChargeFilter) continue;
312 ceballos 1.45 }
313    
314 loizides 1.5 // add good electron
315 ceballos 1.12 GoodElectrons->Add(e);
316 loizides 1.5 }
317 loizides 1.1
318 loizides 1.9 // sort according to pt
319     GoodElectrons->Sort();
320    
321 loizides 1.5 // add to event for other modules to use
322 loizides 1.6 AddObjThisEvt(GoodElectrons);
323 loizides 1.1 }
324    
325     //--------------------------------------------------------------------------------------------------
326     void ElectronIDMod::SlaveBegin()
327     {
328     // Run startup code on the computer (slave) doing the actual analysis. Here,
329 loizides 1.5 // we just request the electron collection branch.
330 loizides 1.1
331 loizides 1.23 ReqEventObject(fElectronBranchName, fElectrons, kTRUE);
332 loizides 1.17
333 ceballos 1.16 if (fApplyConvFilter)
334 loizides 1.23 ReqEventObject(fConversionBranchName, fConversions, kTRUE);
335 loizides 1.17
336 ceballos 1.15 if (fApplyD0Cut)
337 loizides 1.23 ReqEventObject(fVertexName, fVertices, kTRUE);
338 loizides 1.1
339 sixie 1.42 Setup();
340     }
341    
342     //--------------------------------------------------------------------------------------------------
343     void ElectronIDMod::Setup()
344     {
345     // Set all options properly before execution.
346    
347 loizides 1.5 if (fElectronIDType.CompareTo("Tight") == 0)
348     fElIdType = kTight;
349     else if (fElectronIDType.CompareTo("Loose") == 0)
350     fElIdType = kLoose;
351     else if (fElectronIDType.CompareTo("Likelihood") == 0)
352     fElIdType = kLikelihood;
353 loizides 1.10 else if (fElectronIDType.CompareTo("NoId") == 0)
354     fElIdType = kNoId;
355 sixie 1.42 else if (fElectronIDType.CompareTo("ZeeId") == 0)
356     fElIdType = kZeeId;
357 loizides 1.30 else if (fElectronIDType.CompareTo("CustomLoose") == 0) {
358 peveraer 1.29 fElIdType = kCustomIdLoose;
359 loizides 1.30 } else if (fElectronIDType.CompareTo("CustomTight") == 0) {
360 peveraer 1.29 fElIdType = kCustomIdTight;
361 loizides 1.30 }
362 peveraer 1.29 else {
363 loizides 1.5 SendError(kAbortAnalysis, "SlaveBegin",
364     "The specified electron identification %s is not defined.",
365     fElectronIDType.Data());
366     return;
367     }
368    
369 loizides 1.30 SetCustomIDCuts(fElIdType);
370    
371 loizides 1.5 if (fElectronIsoType.CompareTo("TrackCalo") == 0 )
372     fElIsoType = kTrackCalo;
373     else if (fElectronIsoType.CompareTo("TrackJura") == 0)
374     fElIsoType = kTrackJura;
375     else if(fElectronIsoType.CompareTo("TrackJuraSliding") == 0)
376     fElIsoType = kTrackJuraSliding;
377     else if (fElectronIsoType.CompareTo("NoIso") == 0 )
378     fElIsoType = kNoIso;
379 sixie 1.42 else if (fElectronIsoType.CompareTo("ZeeIso") == 0 )
380     fElIsoType = kZeeIso;
381 loizides 1.5 else if (fElectronIsoType.CompareTo("Custom") == 0 ) {
382     fElIsoType = kCustomIso;
383     SendError(kWarning, "SlaveBegin",
384     "Custom electron isolation is not yet implemented.");
385     } else {
386     SendError(kAbortAnalysis, "SlaveBegin",
387     "The specified electron isolation %s is not defined.",
388     fElectronIsoType.Data());
389     return;
390     }
391 sixie 1.42
392 loizides 1.1 }
393 loizides 1.30
394     //--------------------------------------------------------------------------------------------------
395     void ElectronIDMod::SetCustomIDCuts(EElIdType idt)
396     {
397 loizides 1.32 // Set cut values based on RecoEgamma/ElectronIdentification/python/electronIdCutBasedExt_cfi.py.
398     // The following changes are in sigmaetaeta for endcups and deltaetain.
399 loizides 1.30
400     Double_t tightcuts[6][8]={
401 ceballos 1.40 {0.086, 0.1, 0.052, 0.0, 0.050, 0.059, 0.061, 0.0}, //hovere
402 sixie 1.44 {0.011, 0.011, 0.011, 0.0, 0.033, 0.029, 0.030, 0.0}, //sigmaetaeta
403     {0.038, 0.024, 0.045, 0.0, 0.034, 0.017, 0.026, 0.0}, //deltaphiin
404     {0.0081, 0.0029, 0.0051, 0.0, 0.0070, 0.0062, 0.0088, 0.0}, //deltaetain
405 ceballos 1.40 {0.0, 0.9, 0.0, 0.0, 0.0, 0.78, 0.0, 0.0}, //eoverp
406     {0.8,0.2,0.9,0,0,0,0,0}}; //extra cuts fbrem and E_Over_P
407 loizides 1.30
408     Double_t loosecuts[6][8]={
409 loizides 1.39 {0.076, 0.033, 0.07, 0.0, 0.083,0.148, 0.033, 0.0}, //hovere
410     {0.0101, 0.0095, 0.0097, 0.0, 0.03, 0.03, 0.03, 0.0}, //sigmaetaeta
411     {0.053, 0.0189, 0.059, 0.099, 0.0278,0.0157, 0.042, 0.080}, //deltaphiin
412     {0.0078, 0.00259, 0.0062, 0.0, 0.0078,0.0061, 0.0061, 0.0}, //deltaetain
413     {0.3, 0.92, 0.211, 0.0, 0.42, 0.88, 0.68, 0.0}, //eoverp
414     {0.8,0.2,0,0,0,0,0,0}}; //extra cuts fbrem and E_Over_P
415 loizides 1.30
416 sixie 1.42
417 loizides 1.30 switch (idt) {
418     case kCustomIdTight:
419     memcpy(fCuts,tightcuts,sizeof(fCuts));
420     break;
421     case kCustomIdLoose:
422     memcpy(fCuts,loosecuts,sizeof(fCuts));
423     break;
424     default:
425     memset(fCuts,0,sizeof(fCuts));
426     break;
427     }
428     }