ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/IsolationSelection.cc
Revision: 1.4
Committed: Thu Apr 26 06:56:18 2012 UTC (13 years ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.3: +741 -20 lines
Log Message:
first pass port to bambu

File Contents

# User Rev Content
1 khahn 1.1 #include <math.h>
2    
3     #include "IsolationSelection.h"
4     #include "IsolationSelectionDefs.h"
5    
6 khahn 1.4 #include "MathUtils.h"
7     #include "MuonTools.h"
8     #include "MuonIDMVA.h"
9     #include "ElectronTools.h"
10     #include "ElectronIDMVA.h"
11    
12     using namespace mithep;
13    
14     mithep::MuonIDMVA * muIsoMVA;
15     mithep::MuonTools muT;
16     mithep::ElectronIDMVA * eleIsoMVA;
17     mithep::ElectronTools eleT;
18    
19     //--------------------------------------------------------------------------------------------------
20     Float_t computePFMuonIso(const mithep::Muon *muon,
21     const mithep::Vertex & vtx,
22     const mithep::Array<mithep::PFCandidate> * fPFCandidates,
23     const Double_t dRMax)
24     //--------------------------------------------------------------------------------------------------
25     {
26     const Double_t dRMin = 0;
27     const Double_t neuPtMin = 1.0;
28     const Double_t dzMax = 0.1;
29    
30     Double_t zLepton = (muon->BestTrk()) ? muon->BestTrk()->DzCorrected(vtx) : 0.0;
31    
32     Float_t iso=0;
33     for(UInt_t ipf=0; ipf<fPFCandidates->GetEntries(); ipf++) {
34     const PFCandidate *pfcand = fPFCandidates->At(ipf);
35    
36     if(!pfcand->HasTrk() && (pfcand->Pt()<=neuPtMin)) continue; // pT cut on neutral particles
37    
38     // exclude THE muon
39     if(pfcand->TrackerTrk() && muon->TrackerTrk() && (pfcand->TrackerTrk()==muon->TrackerTrk())) continue;
40    
41     // dz cut
42     Double_t dz = (pfcand->BestTrk()) ? fabs(pfcand->BestTrk()->DzCorrected(vtx) - zLepton) : 0;
43     if(dz >= dzMax) continue;
44    
45     // check iso cone
46     Double_t dr = MathUtils::DeltaR(muon->Mom(), pfcand->Mom());
47     if(dr<dRMax && dr>=dRMin)
48     iso += pfcand->Pt();
49     }
50    
51     return iso;
52     }
53    
54     //--------------------------------------------------------------------------------------------------
55     Float_t computePFEleIso(const mithep::Electron *electron,
56     const mithep::Vertex & fVertex,
57     const mithep::Array<mithep::PFCandidate> * fPFCandidates,
58     const Double_t dRMax)
59     //--------------------------------------------------------------------------------------------------
60     {
61     const Double_t dRMin = 0;
62     const Double_t neuPtMin = 1.0;
63     const Double_t dzMax = 0.1;
64    
65     Double_t zLepton = (electron->BestTrk()) ? electron->BestTrk()->DzCorrected(fVertex) : 0.0;
66    
67     Float_t iso=0;
68     for(UInt_t ipf=0; ipf<fPFCandidates->GetEntries(); ipf++) {
69     const PFCandidate *pfcand = (PFCandidate*)(fPFCandidates->At(ipf));
70    
71     if(!pfcand->HasTrk() && (pfcand->Pt()<=neuPtMin)) continue; // pT cut on neutral particles
72    
73     // dz cut
74     Double_t dz = (pfcand->BestTrk()) ? fabs(pfcand->BestTrk()->DzCorrected(fVertex) - zLepton) : 0;
75     if(dz >= dzMax) continue;
76    
77     // remove THE electron
78     if(pfcand->TrackerTrk() && electron->TrackerTrk() && (pfcand->TrackerTrk()==electron->TrackerTrk())) continue;
79     if(pfcand->GsfTrk() && electron->GsfTrk() && (pfcand->GsfTrk()==electron->GsfTrk())) continue;
80    
81     // check iso cone
82     Double_t dr = MathUtils::DeltaR(electron->Mom(), pfcand->Mom());
83     if(dr<dRMax && dr>=dRMin) {
84     // eta-strip veto for photons
85     if((pfcand->PFType() == PFCandidate::eGamma) && fabs(electron->Eta() - pfcand->Eta()) < 0.025) continue;
86    
87     // Inner cone (one tower = dR < 0.07) veto for non-photon neutrals
88     if(!pfcand->HasTrk() && (pfcand->PFType() == PFCandidate::eNeutralHadron) &&
89     (MathUtils::DeltaR(electron->Mom(), pfcand->Mom()) < 0.07)) continue;
90    
91     iso += pfcand->Pt();
92     }
93     }
94    
95     return iso;
96     };
97    
98     //--------------------------------------------------------------------------------------------------
99     bool pairwiseIsoSelection( ControlFlags &ctrl,
100     vector<SimpleLepton> &lepvec,
101     float rho )
102     //--------------------------------------------------------------------------------------------------
103     {
104 khahn 1.2
105     bool passiso=true;
106    
107     for( int i=0; i<lepvec.size(); i++ )
108     {
109    
110     if( !(lepvec[i].is4l) ) continue;
111    
112     float effArea_ecal_i, effArea_hcal_i;
113     if( lepvec[i].isEB ) {
114     if( lepvec[i].type == 11 ) {
115     effArea_ecal_i = 0.101;
116     effArea_hcal_i = 0.021;
117     } else {
118     effArea_ecal_i = 0.074;
119     effArea_hcal_i = 0.022;
120     }
121     } else {
122     if( lepvec[i].type == 11 ) {
123     effArea_ecal_i = 0.046;
124     effArea_hcal_i = 0.040;
125     } else {
126     effArea_ecal_i = 0.045;
127     effArea_hcal_i = 0.030;
128     }
129     }
130    
131     float isoEcal_corr_i = lepvec[i].isoEcal - (effArea_ecal_i*rho);
132     float isoHcal_corr_i = lepvec[i].isoHcal - (effArea_hcal_i*rho);
133    
134     for( int j=i+1; j<lepvec.size(); j++ )
135     {
136    
137     if( !(lepvec[j].is4l) ) continue;
138    
139     float effArea_ecal_j, effArea_hcal_j;
140     if( lepvec[j].isEB ) {
141     if( lepvec[j].type == 11 ) {
142     effArea_ecal_j = 0.101;
143     effArea_hcal_j = 0.021;
144     } else {
145     effArea_ecal_j = 0.074;
146     effArea_hcal_j = 0.022;
147     }
148     } else {
149     if( lepvec[j].type == 11 ) {
150     effArea_ecal_j = 0.046;
151     effArea_hcal_j = 0.040;
152     } else {
153     effArea_ecal_j = 0.045;
154     effArea_hcal_j = 0.030;
155     }
156     }
157    
158     float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
159     float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
160     float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec->Pt();
161     float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec->Pt();
162     float comboIso = RIso_i + RIso_j;
163    
164     if( comboIso > 0.35 ) {
165     if( ctrl.debug ) cout << "combo failing for indices: " << i << "," << j << endl;
166     passiso = false;
167     return passiso;
168     }
169     }
170     }
171    
172     return passiso;
173     }
174    
175 khahn 1.4 //--------------------------------------------------------------------------------------------------
176     SelectionStatus muonIsoSelection(ControlFlags &ctrl,
177     const mithep::Muon * mu,
178     const mithep::Vertex & vtx,
179     const mithep::Array<mithep::PFCandidate> * fPFCandidateCol )
180     //--------------------------------------------------------------------------------------------------
181     {
182     float reliso = computePFMuonIso(mu,vtx,fPFCandidateCol,0.3)/mu->Pt();
183     bool isEB = (fabs(mu->Eta()) < 1.479 ? 1 : 0 );
184 khahn 1.1 bool failiso = false;
185 khahn 1.4 if( isEB && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EB_HIGHPT ) {
186 khahn 1.1 failiso = true;
187     }
188 khahn 1.4 if( isEB && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EB_LOWPT ) {
189 khahn 1.1 failiso = true;
190     }
191 khahn 1.4 if( !(isEB) && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EE_HIGHPT ) {
192 khahn 1.1 failiso = true;
193     }
194 khahn 1.4 if( !(isEB) && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EE_LOWPT ) {
195 khahn 1.1 failiso = true;
196     }
197    
198     SelectionStatus status;
199     if( !failiso ) status.setStatus(SelectionStatus::LOOSEISO);
200     if( !failiso ) status.setStatus(SelectionStatus::TIGHTISO);
201     return status;
202    
203     };
204    
205 khahn 1.4 //--------------------------------------------------------------------------------------------------
206     SelectionStatus electronIsoSelection(ControlFlags &ctrl,
207     const mithep::Electron * ele,
208     const mithep::Vertex &fVertex,
209     const mithep::Array<mithep::PFCandidate> * fPFCandidates)
210     //--------------------------------------------------------------------------------------------------
211     {
212 khahn 1.1
213 khahn 1.4 bool failiso=false;
214 khahn 1.1
215 khahn 1.4 float reliso = computePFEleIso(ele,fVertex,fPFCandidates,0.4)/ele->Pt();
216 khahn 1.1
217 khahn 1.4 if( ele->IsEB() && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EB_HIGHPT ) {
218 khahn 1.1 failiso = true;
219     }
220 khahn 1.4 if( ele->IsEB() && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EB_LOWPT ) {
221 khahn 1.1 failiso = true;
222     }
223     if(ctrl.debug) cout << "before iso check ..." << endl;
224 khahn 1.4 if( !(ele->IsEB()) && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EE_HIGHPT ) {
225 khahn 1.1 if(ctrl.debug) cout << "\tit fails ..." << endl;
226     failiso = true;
227     }
228 khahn 1.4 if( !(ele->IsEB()) && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EE_LOWPT ) {
229 khahn 1.1 failiso = true;
230     }
231    
232     SelectionStatus status;
233     if( !failiso ) {
234 khahn 1.4 status.orStatus(SelectionStatus::LOOSEISO);
235     status.orStatus(SelectionStatus::TIGHTISO);
236 khahn 1.1 }
237     if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
238     return status;
239    
240     }
241 anlevin 1.3
242 khahn 1.4
243 anlevin 1.3 bool noIso(ControlFlags &, vector<SimpleLepton> &, float rho) {
244    
245     return true;
246     }
247 khahn 1.4
248     //--------------------------------------------------------------------------------------------------
249     SelectionStatus muonIsoMVASelection(ControlFlags &ctrl,
250     const mithep::Muon * mu,
251     const mithep::Vertex & vtx,
252     const mithep::Array<mithep::PFCandidate> * fPFCandidates,
253     const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
254     mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
255     vector<const mithep::Muon*> muonsToVeto,
256     vector<const mithep::Electron*> electronsToVeto)
257     //--------------------------------------------------------------------------------------------------
258     {
259    
260     bool failiso=false;
261    
262     //
263     // tmp iso rings
264     //
265     Double_t tmpChargedIso_DR0p0To0p05 = 0;
266     Double_t tmpChargedIso_DR0p05To0p1 = 0;
267     Double_t tmpChargedIso_DR0p1To0p15 = 0;
268     Double_t tmpChargedIso_DR0p15To0p2 = 0;
269     Double_t tmpChargedIso_DR0p2To0p25 = 0;
270     Double_t tmpChargedIso_DR0p25To0p3 = 0;
271     Double_t tmpChargedIso_DR0p3To0p4 = 0;
272     Double_t tmpChargedIso_DR0p4To0p5 = 0;
273    
274     Double_t tmpGammaIso_DR0p0To0p05 = 0;
275     Double_t tmpGammaIso_DR0p05To0p1 = 0;
276     Double_t tmpGammaIso_DR0p1To0p15 = 0;
277     Double_t tmpGammaIso_DR0p15To0p2 = 0;
278     Double_t tmpGammaIso_DR0p2To0p25 = 0;
279     Double_t tmpGammaIso_DR0p25To0p3 = 0;
280     Double_t tmpGammaIso_DR0p3To0p4 = 0;
281     Double_t tmpGammaIso_DR0p4To0p5 = 0;
282    
283     Double_t tmpNeutralHadronIso_DR0p0To0p05 = 0;
284     Double_t tmpNeutralHadronIso_DR0p05To0p1 = 0;
285     Double_t tmpNeutralHadronIso_DR0p1To0p15 = 0;
286     Double_t tmpNeutralHadronIso_DR0p15To0p2 = 0;
287     Double_t tmpNeutralHadronIso_DR0p2To0p25 = 0;
288     Double_t tmpNeutralHadronIso_DR0p25To0p3 = 0;
289     Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
290     Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
291    
292     Double_t tmp2ChargedIso_DR0p5To1p0 = 0;
293    
294     //
295     // final rings for the MVA
296     //
297     Double_t fChargedIso_DR0p0To0p1;
298     Double_t fChargedIso_DR0p1To0p2;
299     Double_t fChargedIso_DR0p2To0p3;
300     Double_t fChargedIso_DR0p3To0p4;
301     Double_t fChargedIso_DR0p4To0p5;
302    
303     Double_t fGammaIso_DR0p0To0p1;
304     Double_t fGammaIso_DR0p1To0p2;
305     Double_t fGammaIso_DR0p2To0p3;
306     Double_t fGammaIso_DR0p3To0p4;
307     Double_t fGammaIso_DR0p4To0p5;
308    
309     Double_t fNeutralHadronIso_DR0p0To0p1;
310     Double_t fNeutralHadronIso_DR0p1To0p2;
311     Double_t fNeutralHadronIso_DR0p2To0p3;
312     Double_t fNeutralHadronIso_DR0p3To0p4;
313     Double_t fNeutralHadronIso_DR0p4To0p5;
314    
315    
316     //
317     //Loop over PF Candidates
318     //
319     for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
320     const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
321    
322     Double_t deta = (mu->Eta() - pf->Eta());
323     Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
324     Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
325     if (dr > 1.0) continue;
326    
327     if (pf->PFType() == PFCandidate::eMuon && pf->TrackerTrk() == mu->TrackerTrk() ) continue;
328    
329     //
330     // Lepton Footprint Removal
331     //
332     Bool_t IsLeptonFootprint = kFALSE;
333     if (dr < 1.0) {
334    
335     //
336     // Check for electrons
337     //
338     for (Int_t q=0; q < electronsToVeto.size(); ++q) {
339     const mithep::Electron *tmpele = electronsToVeto[q];
340     // PF electron
341     if( pf->PFType() == PFCandidate::eElectron && pf->TrackerTrk() == tmpele->TrackerTrk() )
342     IsLeptonFootprint = kTRUE;
343     // PF charged
344     if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
345     && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
346     IsLeptonFootprint = kTRUE;
347     // PF gamma
348     if (pf->PFType() == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
349     && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
350     IsLeptonFootprint = kTRUE;
351     } // loop over electrons
352    
353     //
354     // Check for muons
355     //
356     for (Int_t q=0; q < muonsToVeto.size(); ++q) {
357     const mithep::Muon *tmpmu = muonsToVeto[q];
358     // PF muons
359     if (pf->PFType() == PFCandidate::eMuon && pf->TrackerTrk() == tmpmu->TrackerTrk() )
360     IsLeptonFootprint = kTRUE;
361     // PF charged
362     if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
363     IsLeptonFootprint = kTRUE;
364     } // loop over muons
365    
366    
367     if (IsLeptonFootprint)
368     continue;
369    
370     //
371     // Charged Iso Rings
372     //
373     if (pf->Charge() != 0 && pf->HasTrackerTrk() ) {
374    
375     if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
376    
377     // Veto any PFmuon, or PFEle
378     if (pf->PFType() == PFCandidate::eElectron || pf->PFType() == PFCandidate::eMuon) continue;
379    
380     // Footprint Veto
381     if (dr < 0.01) continue;
382    
383     if (dr < 0.05) tmpChargedIso_DR0p0To0p05 += pf->Pt();
384     if (dr >= 0.05 && dr < 0.10) tmpChargedIso_DR0p05To0p1 += pf->Pt();
385     if (dr >= 0.10 && dr < 0.15) tmpChargedIso_DR0p1To0p15 += pf->Pt();
386     if (dr >= 0.15 && dr < 0.20) tmpChargedIso_DR0p15To0p2 += pf->Pt();
387     if (dr >= 0.20 && dr < 0.25) tmpChargedIso_DR0p2To0p25 += pf->Pt();
388     if (dr >= 0.25 && dr < 0.3) tmpChargedIso_DR0p25To0p3 += pf->Pt();
389     if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
390     if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
391     }
392    
393     //
394     // Gamma Iso Rings
395     //
396     else if (pf->PFType() == PFCandidate::eGamma) {
397    
398     if (dr < 0.05) tmpGammaIso_DR0p0To0p05 += pf->Pt();
399     if (dr >= 0.05 && dr < 0.10) tmpGammaIso_DR0p05To0p1 += pf->Pt();
400     if (dr >= 0.10 && dr < 0.15) tmpGammaIso_DR0p1To0p15 += pf->Pt();
401     if (dr >= 0.15 && dr < 0.20) tmpGammaIso_DR0p15To0p2 += pf->Pt();
402     if (dr >= 0.20 && dr < 0.25) tmpGammaIso_DR0p2To0p25 += pf->Pt();
403     if (dr >= 0.25 && dr < 0.3) tmpGammaIso_DR0p25To0p3 += pf->Pt();
404     if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
405     if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
406     }
407    
408     //
409     // Other Neutral Iso Rings
410     //
411     else {
412     if (dr < 0.05) tmpNeutralHadronIso_DR0p0To0p05 += pf->Pt();
413     if (dr >= 0.05 && dr < 0.10) tmpNeutralHadronIso_DR0p05To0p1 += pf->Pt();
414     if (dr >= 0.10 && dr < 0.15) tmpNeutralHadronIso_DR0p1To0p15 += pf->Pt();
415     if (dr >= 0.15 && dr < 0.20) tmpNeutralHadronIso_DR0p15To0p2 += pf->Pt();
416     if (dr >= 0.20 && dr < 0.25) tmpNeutralHadronIso_DR0p2To0p25 += pf->Pt();
417     if (dr >= 0.25 && dr < 0.3) tmpNeutralHadronIso_DR0p25To0p3 += pf->Pt();
418     if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
419     if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
420     }
421    
422     }
423    
424     }
425    
426     fChargedIso_DR0p0To0p1 = min((tmpChargedIso_DR0p0To0p05 + tmpChargedIso_DR0p05To0p1 )/mu->Pt(), 2.5);
427     fChargedIso_DR0p1To0p2 = min((tmpChargedIso_DR0p1To0p15 + tmpChargedIso_DR0p15To0p2)/mu->Pt(), 2.5);
428     fChargedIso_DR0p2To0p3 = min((tmpChargedIso_DR0p2To0p25 + tmpChargedIso_DR0p25To0p3)/mu->Pt(), 2.5);
429     fChargedIso_DR0p3To0p4 = min((tmpChargedIso_DR0p3To0p4)/mu->Pt(), 2.5);
430     fChargedIso_DR0p4To0p5 = min((tmpChargedIso_DR0p4To0p5)/mu->Pt(), 2.5);
431    
432     double rho = 0;
433     if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
434     rho = fPUEnergyDensity->At(0)->Rho();
435    
436    
437     fGammaIso_DR0p0To0p1 = max(min((tmpGammaIso_DR0p0To0p05 + tmpGammaIso_DR0p05To0p1
438     -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p0To0p1,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
439     ,2.5)
440     ,0.0);
441     fGammaIso_DR0p1To0p2 = max(min((tmpGammaIso_DR0p1To0p15 + tmpGammaIso_DR0p15To0p2
442     -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p1To0p2,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
443     ,2.5)
444     ,0.0);
445     fGammaIso_DR0p2To0p3 = max(min((tmpGammaIso_DR0p2To0p25 + tmpGammaIso_DR0p25To0p3
446     -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p2To0p3,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
447     ,2.5)
448     ,0.0);
449     fGammaIso_DR0p3To0p4 = max(min((tmpGammaIso_DR0p3To0p4
450     -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p3To0p4,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
451     ,2.5)
452     ,0.0);
453     fGammaIso_DR0p4To0p5 = max(min((tmpGammaIso_DR0p4To0p5
454     -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p4To0p5,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
455     ,2.5)
456     ,0.0);
457    
458    
459     fNeutralHadronIso_DR0p0To0p1 = max(min((tmpNeutralHadronIso_DR0p0To0p05 + tmpNeutralHadronIso_DR0p05To0p1
460     -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p0To0p1,
461     mu->Eta(),EffectiveAreaVersion))/mu->Pt()
462     , 2.5)
463     , 0.0);
464     fNeutralHadronIso_DR0p1To0p2 = max(min((tmpNeutralHadronIso_DR0p1To0p15 + tmpNeutralHadronIso_DR0p15To0p2
465     -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p1To0p2,
466     mu->Eta(),EffectiveAreaVersion))/mu->Pt()
467     , 2.5)
468     , 0.0);
469     fNeutralHadronIso_DR0p2To0p3 = max(min((tmpNeutralHadronIso_DR0p2To0p25 + tmpNeutralHadronIso_DR0p25To0p3
470     -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p2To0p3,
471     mu->Eta(),EffectiveAreaVersion))/mu->Pt()
472     , 2.5)
473     , 0.0);
474     fNeutralHadronIso_DR0p3To0p4 = max(min((tmpNeutralHadronIso_DR0p3To0p4
475     -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p3To0p4,
476     mu->Eta(), EffectiveAreaVersion))/mu->Pt()
477     , 2.5)
478     , 0.0);
479     fNeutralHadronIso_DR0p4To0p5 = max(min((tmpNeutralHadronIso_DR0p4To0p5
480     -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p4To0p5,
481     mu->Eta(), EffectiveAreaVersion))/mu->Pt()
482     , 2.5)
483     , 0.0);
484    
485     double mvaval = muIsoMVA->MVAValue_IsoRings( mu->Pt(),
486     mu->Eta(),
487     fChargedIso_DR0p0To0p1,
488     fChargedIso_DR0p1To0p2,
489     fChargedIso_DR0p2To0p3,
490     fChargedIso_DR0p3To0p4,
491     fChargedIso_DR0p4To0p5,
492     fGammaIso_DR0p0To0p1,
493     fGammaIso_DR0p1To0p2,
494     fGammaIso_DR0p2To0p3,
495     fGammaIso_DR0p3To0p4,
496     fGammaIso_DR0p4To0p5,
497     fNeutralHadronIso_DR0p0To0p1,
498     fNeutralHadronIso_DR0p1To0p2,
499     fNeutralHadronIso_DR0p2To0p3,
500     fNeutralHadronIso_DR0p3To0p4,
501     fNeutralHadronIso_DR0p4To0p5,
502     ctrl.debug);
503    
504     SelectionStatus status;
505     bool pass = false;
506    
507     if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
508     && fabs(mu->Eta()) < 1.5 && mu->Pt() < 10 && mvaval >= MUON_ISOMVA_CUT_BIN0) pass = true;
509     else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
510     && fabs(mu->Eta()) < 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_CUT_BIN1) pass = true;
511     else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
512     && fabs(mu->Eta()) > 1.5 && mu->Pt() < 10 && mvaval >= MUON_ISOMVA_CUT_BIN2) pass = true;
513     else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
514     && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_CUT_BIN3) pass = true;
515     else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon()
516     && (mu->Quality().QualityMask().Mask() & mithep::MuonQuality::AllArbitrated) && mvaval >= MUON_ISOMVA_CUT_BIN4)
517     pass = true;
518    
519    
520     if( pass ) {
521     status.orStatus(SelectionStatus::LOOSEISO);
522     status.orStatus(SelectionStatus::TIGHTISO);
523     }
524     if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
525     return status;
526    
527     }
528    
529     //--------------------------------------------------------------------------------------------------
530     void initMuonIsoMVA() {
531     //--------------------------------------------------------------------------------------------------
532     muIsoMVA = new mithep::MuonIDMVA();
533     vector<string> weightFiles;
534     weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml");
535     weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml");
536     weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml");
537     weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml");
538     weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml");
539     weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml");
540     muIsoMVA->Initialize( "MuonIsoMVA",
541     mithep::MuonIDMVA::kIsoRingsV0,
542     kTRUE, weightFiles);
543     }
544    
545    
546    
547     //--------------------------------------------------------------------------------------------------
548     SelectionStatus electronIsoMVASelection(ControlFlags &ctrl,
549     const mithep::Electron * ele,
550     const mithep::Vertex & vtx,
551     const mithep::Array<mithep::PFCandidate> * fPFCandidates,
552     const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
553     mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
554     vector<const mithep::Muon*> muonsToVeto,
555     vector<const mithep::Electron*> electronsToVeto)
556     //--------------------------------------------------------------------------------------------------
557     {
558    
559     bool failiso=false;
560    
561     //
562     // tmp iso rings
563     //
564     Double_t tmpChargedIso_DR0p0To0p05 = 0;
565     Double_t tmpChargedIso_DR0p05To0p1 = 0;
566     Double_t tmpChargedIso_DR0p1To0p15 = 0;
567     Double_t tmpChargedIso_DR0p15To0p2 = 0;
568     Double_t tmpChargedIso_DR0p2To0p25 = 0;
569     Double_t tmpChargedIso_DR0p25To0p3 = 0;
570     Double_t tmpChargedIso_DR0p3To0p4 = 0;
571     Double_t tmpChargedIso_DR0p4To0p5 = 0;
572    
573     Double_t tmpGammaIso_DR0p0To0p05 = 0;
574     Double_t tmpGammaIso_DR0p05To0p1 = 0;
575     Double_t tmpGammaIso_DR0p1To0p15 = 0;
576     Double_t tmpGammaIso_DR0p15To0p2 = 0;
577     Double_t tmpGammaIso_DR0p2To0p25 = 0;
578     Double_t tmpGammaIso_DR0p25To0p3 = 0;
579     Double_t tmpGammaIso_DR0p3To0p4 = 0;
580     Double_t tmpGammaIso_DR0p4To0p5 = 0;
581    
582     Double_t tmpNeutralHadronIso_DR0p0To0p05 = 0;
583     Double_t tmpNeutralHadronIso_DR0p05To0p1 = 0;
584     Double_t tmpNeutralHadronIso_DR0p1To0p15 = 0;
585     Double_t tmpNeutralHadronIso_DR0p15To0p2 = 0;
586     Double_t tmpNeutralHadronIso_DR0p2To0p25 = 0;
587     Double_t tmpNeutralHadronIso_DR0p25To0p3 = 0;
588     Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
589     Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
590    
591     Double_t tmp2ChargedIso_DR0p5To1p0 = 0;
592    
593     //
594     // final rings for the MVA
595     //
596     Double_t fChargedIso_DR0p0To0p1;
597     Double_t fChargedIso_DR0p1To0p2;
598     Double_t fChargedIso_DR0p2To0p3;
599     Double_t fChargedIso_DR0p3To0p4;
600     Double_t fChargedIso_DR0p4To0p5;
601    
602     Double_t fGammaIso_DR0p0To0p1;
603     Double_t fGammaIso_DR0p1To0p2;
604     Double_t fGammaIso_DR0p2To0p3;
605     Double_t fGammaIso_DR0p3To0p4;
606     Double_t fGammaIso_DR0p4To0p5;
607    
608     Double_t fNeutralHadronIso_DR0p0To0p1;
609     Double_t fNeutralHadronIso_DR0p1To0p2;
610     Double_t fNeutralHadronIso_DR0p2To0p3;
611     Double_t fNeutralHadronIso_DR0p3To0p4;
612     Double_t fNeutralHadronIso_DR0p4To0p5;
613    
614    
615     //
616     //Loop over PF Candidates
617     //
618     for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
619     const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
620    
621     Double_t deta = (ele->Eta() - pf->Eta());
622     Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
623     Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
624     if (dr > 1.0) continue;
625    
626     if (pf->PFType() == PFCandidate::eElectron && pf->TrackerTrk() == ele->TrackerTrk() ) continue;
627    
628     //
629     // Lepton Footprint Removal
630     //
631     Bool_t IsLeptonFootprint = kFALSE;
632     if (dr < 1.0) {
633    
634     //
635     // Check for electrons
636     //
637     for (Int_t q=0; q < electronsToVeto.size(); ++q) {
638     const mithep::Electron *tmpele = electronsToVeto[q];
639     // PF electron
640     if( pf->PFType() == PFCandidate::eElectron && pf->TrackerTrk() == tmpele->TrackerTrk() )
641     IsLeptonFootprint = kTRUE;
642     // PF charged
643     if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
644     && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
645     IsLeptonFootprint = kTRUE;
646     // PF gamma
647     if (pf->PFType() == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
648     && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
649     IsLeptonFootprint = kTRUE;
650     } // loop over electrons
651    
652     //
653     // Check for muons
654     //
655     for (Int_t q=0; q < muonsToVeto.size(); ++q) {
656     const mithep::Muon *tmpmu = muonsToVeto[q];
657     // PF muons
658     if (pf->PFType() == PFCandidate::eMuon && pf->TrackerTrk() == tmpmu->TrackerTrk() )
659     IsLeptonFootprint = kTRUE;
660     // PF charged
661     if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
662     IsLeptonFootprint = kTRUE;
663     } // loop over muons
664    
665    
666     if (IsLeptonFootprint)
667     continue;
668    
669     //
670     // Charged Iso Rings
671     //
672     if (pf->Charge() != 0 && pf->HasTrackerTrk() ) {
673    
674     if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
675    
676     // Veto any PFmuon, or PFEle
677     if (pf->PFType() == PFCandidate::eElectron || pf->PFType() == PFCandidate::eMuon) continue;
678    
679     // Footprint Veto
680     if (dr < 0.01) continue;
681    
682     if (dr < 0.05) tmpChargedIso_DR0p0To0p05 += pf->Pt();
683     if (dr >= 0.05 && dr < 0.10) tmpChargedIso_DR0p05To0p1 += pf->Pt();
684     if (dr >= 0.10 && dr < 0.15) tmpChargedIso_DR0p1To0p15 += pf->Pt();
685     if (dr >= 0.15 && dr < 0.20) tmpChargedIso_DR0p15To0p2 += pf->Pt();
686     if (dr >= 0.20 && dr < 0.25) tmpChargedIso_DR0p2To0p25 += pf->Pt();
687     if (dr >= 0.25 && dr < 0.3) tmpChargedIso_DR0p25To0p3 += pf->Pt();
688     if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
689     if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
690     }
691    
692     //
693     // Gamma Iso Rings
694     //
695     else if (pf->PFType() == PFCandidate::eGamma) {
696    
697     if (dr < 0.05) tmpGammaIso_DR0p0To0p05 += pf->Pt();
698     if (dr >= 0.05 && dr < 0.10) tmpGammaIso_DR0p05To0p1 += pf->Pt();
699     if (dr >= 0.10 && dr < 0.15) tmpGammaIso_DR0p1To0p15 += pf->Pt();
700     if (dr >= 0.15 && dr < 0.20) tmpGammaIso_DR0p15To0p2 += pf->Pt();
701     if (dr >= 0.20 && dr < 0.25) tmpGammaIso_DR0p2To0p25 += pf->Pt();
702     if (dr >= 0.25 && dr < 0.3) tmpGammaIso_DR0p25To0p3 += pf->Pt();
703     if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
704     if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
705     }
706    
707     //
708     // Other Neutral Iso Rings
709     //
710     else {
711     if (dr < 0.05) tmpNeutralHadronIso_DR0p0To0p05 += pf->Pt();
712     if (dr >= 0.05 && dr < 0.10) tmpNeutralHadronIso_DR0p05To0p1 += pf->Pt();
713     if (dr >= 0.10 && dr < 0.15) tmpNeutralHadronIso_DR0p1To0p15 += pf->Pt();
714     if (dr >= 0.15 && dr < 0.20) tmpNeutralHadronIso_DR0p15To0p2 += pf->Pt();
715     if (dr >= 0.20 && dr < 0.25) tmpNeutralHadronIso_DR0p2To0p25 += pf->Pt();
716     if (dr >= 0.25 && dr < 0.3) tmpNeutralHadronIso_DR0p25To0p3 += pf->Pt();
717     if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
718     if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
719     }
720    
721     }
722    
723     }
724    
725     fChargedIso_DR0p0To0p1 = min((tmpChargedIso_DR0p0To0p05 + tmpChargedIso_DR0p05To0p1 )/ele->Pt(), 2.5);
726     fChargedIso_DR0p1To0p2 = min((tmpChargedIso_DR0p1To0p15 + tmpChargedIso_DR0p15To0p2)/ele->Pt(), 2.5);
727     fChargedIso_DR0p2To0p3 = min((tmpChargedIso_DR0p2To0p25 + tmpChargedIso_DR0p25To0p3)/ele->Pt(), 2.5);
728     fChargedIso_DR0p3To0p4 = min((tmpChargedIso_DR0p3To0p4)/ele->Pt(), 2.5);
729     fChargedIso_DR0p4To0p5 = min((tmpChargedIso_DR0p4To0p5)/ele->Pt(), 2.5);
730    
731     double rho = 0;
732     if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
733     rho = fPUEnergyDensity->At(0)->Rho();
734    
735    
736     fGammaIso_DR0p0To0p1 = max(min((tmpGammaIso_DR0p0To0p05 + tmpGammaIso_DR0p05To0p1
737     -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p0To0p1,
738     ele->Eta(),
739     EffectiveAreaVersion))/ele->Pt()
740     ,2.5)
741     ,0.0);
742     fGammaIso_DR0p1To0p2 = max(min((tmpGammaIso_DR0p1To0p15 + tmpGammaIso_DR0p15To0p2
743     -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p1To0p2,
744     ele->Eta(),
745     EffectiveAreaVersion))/ele->Pt()
746     ,2.5)
747     ,0.0);
748     fGammaIso_DR0p2To0p3 = max(min((tmpGammaIso_DR0p2To0p25 + tmpGammaIso_DR0p25To0p3
749     -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p2To0p3,
750     ele->Eta()
751     ,EffectiveAreaVersion))/ele->Pt()
752     ,2.5)
753     ,0.0);
754     fGammaIso_DR0p3To0p4 = max(min((tmpGammaIso_DR0p3To0p4
755     -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p3To0p4,
756     ele->Eta(),
757     EffectiveAreaVersion))/ele->Pt()
758     ,2.5)
759     ,0.0);
760     fGammaIso_DR0p4To0p5 = max(min((tmpGammaIso_DR0p4To0p5
761     -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p4To0p5,
762     ele->Eta(),
763     EffectiveAreaVersion))/ele->Pt()
764     ,2.5)
765     ,0.0);
766    
767    
768     fNeutralHadronIso_DR0p0To0p1 = max(min((tmpNeutralHadronIso_DR0p0To0p05 + tmpNeutralHadronIso_DR0p05To0p1
769     -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
770     ele->Eta(),EffectiveAreaVersion))/ele->Pt()
771     , 2.5)
772     , 0.0);
773     fNeutralHadronIso_DR0p1To0p2 = max(min((tmpNeutralHadronIso_DR0p1To0p15 + tmpNeutralHadronIso_DR0p15To0p2
774     -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
775     ele->Eta(),EffectiveAreaVersion))/ele->Pt()
776     , 2.5)
777     , 0.0);
778     fNeutralHadronIso_DR0p2To0p3 = max(min((tmpNeutralHadronIso_DR0p2To0p25 + tmpNeutralHadronIso_DR0p25To0p3
779     -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
780     ele->Eta(),EffectiveAreaVersion))/ele->Pt()
781     , 2.5)
782     , 0.0);
783     fNeutralHadronIso_DR0p3To0p4 = max(min((tmpNeutralHadronIso_DR0p3To0p4
784     -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
785     ele->Eta(), EffectiveAreaVersion))/ele->Pt()
786     , 2.5)
787     , 0.0);
788     fNeutralHadronIso_DR0p4To0p5 = max(min((tmpNeutralHadronIso_DR0p4To0p5
789     -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p4To0p5,
790     ele->Eta(), EffectiveAreaVersion))/ele->Pt()
791     , 2.5)
792     , 0.0);
793    
794     double mvaval = eleIsoMVA->MVAValue_IsoRings( ele->Pt(),
795     ele->Eta(),
796     fChargedIso_DR0p0To0p1,
797     fChargedIso_DR0p1To0p2,
798     fChargedIso_DR0p2To0p3,
799     fChargedIso_DR0p3To0p4,
800     fChargedIso_DR0p4To0p5,
801     fGammaIso_DR0p0To0p1,
802     fGammaIso_DR0p1To0p2,
803     fGammaIso_DR0p2To0p3,
804     fGammaIso_DR0p3To0p4,
805     fGammaIso_DR0p4To0p5,
806     fNeutralHadronIso_DR0p0To0p1,
807     fNeutralHadronIso_DR0p1To0p2,
808     fNeutralHadronIso_DR0p2To0p3,
809     fNeutralHadronIso_DR0p3To0p4,
810     fNeutralHadronIso_DR0p4To0p5,
811     ctrl.debug);
812    
813     SelectionStatus status;
814     bool pass = false;
815    
816     Int_t subdet = 0;
817     if (fabs(ele->SCluster()->Eta()) < 0.8) subdet = 0;
818     else if (fabs(ele->SCluster()->Eta()) < 1.479) subdet = 1;
819     else subdet = 2;
820     Int_t ptBin = 0;
821     if (ele->Pt() > 10.0) ptBin = 1;
822    
823     Int_t MVABin = -1;
824     if (subdet == 0 && ptBin == 0) MVABin = 0;
825     if (subdet == 1 && ptBin == 0) MVABin = 1;
826     if (subdet == 2 && ptBin == 0) MVABin = 2;
827     if (subdet == 0 && ptBin == 1) MVABin = 3;
828     if (subdet == 1 && ptBin == 1) MVABin = 4;
829     if (subdet == 2 && ptBin == 1) MVABin = 5;
830    
831     if( MVABin == 0 && mvaval > ELECTRON_ISOMVA_CUT_BIN0 ) pass = true;
832     if( MVABin == 1 && mvaval > ELECTRON_ISOMVA_CUT_BIN1 ) pass = true;
833     if( MVABin == 2 && mvaval > ELECTRON_ISOMVA_CUT_BIN2 ) pass = true;
834     if( MVABin == 3 && mvaval > ELECTRON_ISOMVA_CUT_BIN3 ) pass = true;
835     if( MVABin == 4 && mvaval > ELECTRON_ISOMVA_CUT_BIN4 ) pass = true;
836     if( MVABin == 5 && mvaval > ELECTRON_ISOMVA_CUT_BIN5 ) pass = true;
837    
838     if( pass ) {
839     status.orStatus(SelectionStatus::LOOSEISO);
840     status.orStatus(SelectionStatus::TIGHTISO);
841     }
842     if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
843     return status;
844    
845     }
846    
847    
848     //--------------------------------------------------------------------------------------------------
849     void initElectronIsoMVA() {
850     //--------------------------------------------------------------------------------------------------
851     eleIsoMVA = new mithep::ElectronIDMVA();
852     vector<string> weightFiles;
853     weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt5To10.weights.xml");
854     weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt5To10.weights.xml");
855     weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt10ToInf.weights.xml");
856     weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt10ToInf.weights.xml");
857     eleIsoMVA->Initialize( "ElectronIsoMVA",
858     mithep::ElectronIDMVA::kIsoRingsV0,
859     kTRUE, weightFiles);
860     }