ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/IsolationSelection.cc
Revision: 1.6
Committed: Thu Apr 26 07:15:08 2012 UTC (13 years ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.5: +6 -6 lines
Log Message:
*** empty log message ***

File Contents

# Content
1 #include <math.h>
2
3 #include "IsolationSelection.h"
4 #include "IsolationSelectionDefs.h"
5
6 #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
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 //--------------------------------------------------------------------------------------------------
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 bool failiso = false;
185 if( isEB && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EB_HIGHPT ) {
186 failiso = true;
187 }
188 if( isEB && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EB_LOWPT ) {
189 failiso = true;
190 }
191 if( !(isEB) && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EE_HIGHPT ) {
192 failiso = true;
193 }
194 if( !(isEB) && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EE_LOWPT ) {
195 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 //--------------------------------------------------------------------------------------------------
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
213 bool failiso=false;
214
215 float reliso = computePFEleIso(ele,fVertex,fPFCandidates,0.4)/ele->Pt();
216
217 if( ele->IsEB() && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EB_HIGHPT ) {
218 failiso = true;
219 }
220 if( ele->IsEB() && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EB_LOWPT ) {
221 failiso = true;
222 }
223 if(ctrl.debug) cout << "before iso check ..." << endl;
224 if( !(ele->IsEB()) && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EE_HIGHPT ) {
225 if(ctrl.debug) cout << "\tit fails ..." << endl;
226 failiso = true;
227 }
228 if( !(ele->IsEB()) && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EE_LOWPT ) {
229 failiso = true;
230 }
231
232 SelectionStatus status;
233 if( !failiso ) {
234 status.orStatus(SelectionStatus::LOOSEISO);
235 status.orStatus(SelectionStatus::TIGHTISO);
236 }
237 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
238 return status;
239
240 }
241
242
243 bool noIso(ControlFlags &, vector<SimpleLepton> &, float rho) {
244
245 return true;
246 }
247
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 }