ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/IsolationSelection.cc
Revision: 1.12
Committed: Sun May 6 12:44:45 2012 UTC (13 years ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.11: +6 -6 lines
Log Message:
initialze isos for PF04

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 if( ctrl.debug ) {
261 cout << "muonIsoMVASelection :: muons to veto " << endl;
262 for( int i=0; i<muonsToVeto.size(); i++ ) {
263 const mithep::Muon * vmu = muonsToVeto[i];
264 cout << "\tpt: " << vmu->Pt()
265 << "\teta: " << vmu->Eta()
266 << "\tphi: " << vmu->Phi()
267 << endl;
268 }
269 cout << "muonIsoMVASelection :: electrson to veto " << endl;
270 for( int i=0; i<electronsToVeto.size(); i++ ) {
271 const mithep::Electron * vel = electronsToVeto[i];
272 cout << "\tpt: " << vel->Pt()
273 << "\teta: " << vel->Eta()
274 << "\tphi: " << vel->Phi()
275 << endl;
276 }
277 }
278 bool failiso=false;
279
280 //
281 // tmp iso rings
282 //
283 Double_t tmpChargedIso_DR0p0To0p1 = 0;
284 Double_t tmpChargedIso_DR0p1To0p2 = 0;
285 Double_t tmpChargedIso_DR0p2To0p3 = 0;
286 Double_t tmpChargedIso_DR0p3To0p4 = 0;
287 Double_t tmpChargedIso_DR0p4To0p5 = 0;
288 Double_t tmpChargedIso_DR0p5To0p7 = 0;
289
290 Double_t tmpGammaIso_DR0p0To0p1 = 0;
291 Double_t tmpGammaIso_DR0p1To0p2 = 0;
292 Double_t tmpGammaIso_DR0p2To0p3 = 0;
293 Double_t tmpGammaIso_DR0p3To0p4 = 0;
294 Double_t tmpGammaIso_DR0p4To0p5 = 0;
295 Double_t tmpGammaIso_DR0p5To0p7 = 0;
296
297 Double_t tmpNeutralHadronIso_DR0p0To0p1 = 0;
298 Double_t tmpNeutralHadronIso_DR0p1To0p2 = 0;
299 Double_t tmpNeutralHadronIso_DR0p2To0p3 = 0;
300 Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
301 Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
302 Double_t tmpNeutralHadronIso_DR0p5To0p7 = 0;
303
304
305
306 //
307 // final rings for the MVA
308 //
309 Double_t fChargedIso_DR0p0To0p1;
310 Double_t fChargedIso_DR0p1To0p2;
311 Double_t fChargedIso_DR0p2To0p3;
312 Double_t fChargedIso_DR0p3To0p4;
313 Double_t fChargedIso_DR0p4To0p5;
314 Double_t fChargedIso_DR0p5To0p7;
315
316 Double_t fGammaIso_DR0p0To0p1;
317 Double_t fGammaIso_DR0p1To0p2;
318 Double_t fGammaIso_DR0p2To0p3;
319 Double_t fGammaIso_DR0p3To0p4;
320 Double_t fGammaIso_DR0p4To0p5;
321 Double_t fGammaIso_DR0p5To0p7;
322
323 Double_t fNeutralHadronIso_DR0p0To0p1;
324 Double_t fNeutralHadronIso_DR0p1To0p2;
325 Double_t fNeutralHadronIso_DR0p2To0p3;
326 Double_t fNeutralHadronIso_DR0p3To0p4;
327 Double_t fNeutralHadronIso_DR0p4To0p5;
328 Double_t fNeutralHadronIso_DR0p5To0p7;
329
330
331 //
332 //Loop over PF Candidates
333 //
334 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
335 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
336
337 Double_t deta = (mu->Eta() - pf->Eta());
338 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
339 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
340 if (dr > 1.0) continue;
341
342 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
343
344 //
345 // Lepton Footprint Removal
346 //
347 Bool_t IsLeptonFootprint = kFALSE;
348 if (dr < 1.0) {
349
350 //
351 // Check for electrons
352 //
353 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
354 const mithep::Electron *tmpele = electronsToVeto[q];
355 // 4l electron
356 if( pf->HasTrackerTrk() ) {
357 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
358 IsLeptonFootprint = kTRUE;
359 }
360 if( pf->HasGsfTrk() ) {
361 if( pf->GsfTrk() == tmpele->GsfTrk() )
362 IsLeptonFootprint = kTRUE;
363 }
364 // PF charged
365 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
366 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
367 IsLeptonFootprint = kTRUE;
368 // PF gamma
369 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
370 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
371 IsLeptonFootprint = kTRUE;
372 } // loop over electrons
373
374 //
375 // Check for muons
376 //
377 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
378 const mithep::Muon *tmpmu = muonsToVeto[q];
379 // 4l muon
380 if( pf->HasTrackerTrk() ) {
381 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
382 IsLeptonFootprint = kTRUE;
383 }
384 // PF charged
385 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
386 IsLeptonFootprint = kTRUE;
387 } // loop over muons
388
389
390 if (IsLeptonFootprint)
391 continue;
392
393 //
394 // Charged Iso Rings
395 //
396 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
397
398 if( dr < 0.01 ) continue; // only for muon iso mva?
399 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
400
401 if( pf->HasTrackerTrk() ) {
402 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
403 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
404 << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
405 << dr << endl;
406 }
407 if( pf->HasGsfTrk() ) {
408 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
409 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
410 << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
411 << dr << endl;
412 }
413
414 // Footprint Veto
415 if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
416 if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
417 if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
418 if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
419 if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
420 if (dr >= 0.5 && dr < 0.7) tmpChargedIso_DR0p5To0p7 += pf->Pt();
421 }
422
423 //
424 // Gamma Iso Rings
425 //
426 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
427 if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
428 if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
429 if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
430 if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
431 if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
432 if (dr >= 0.5 && dr < 0.7) tmpGammaIso_DR0p5To0p7 += pf->Pt();
433 }
434
435 //
436 // Other Neutral Iso Rings
437 //
438 else {
439 if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
440 if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
441 if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
442 if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
443 if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
444 if (dr >= 0.5 && dr < 0.7) tmpNeutralHadronIso_DR0p5To0p7 += pf->Pt();
445 }
446
447 }
448
449 }
450
451 fChargedIso_DR0p0To0p1 = min((tmpChargedIso_DR0p0To0p1)/mu->Pt(), 2.5);
452 fChargedIso_DR0p1To0p2 = min((tmpChargedIso_DR0p1To0p2)/mu->Pt(), 2.5);
453 fChargedIso_DR0p2To0p3 = min((tmpChargedIso_DR0p2To0p3)/mu->Pt(), 2.5);
454 fChargedIso_DR0p3To0p4 = min((tmpChargedIso_DR0p3To0p4)/mu->Pt(), 2.5);
455 fChargedIso_DR0p4To0p5 = min((tmpChargedIso_DR0p4To0p5)/mu->Pt(), 2.5);
456
457
458 double rho = 0;
459 if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
460 rho = fPUEnergyDensity->At(0)->Rho();
461
462
463 fGammaIso_DR0p0To0p1 = max(min((tmpGammaIso_DR0p0To0p1
464 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p0To0p1,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
465 ,2.5)
466 ,0.0);
467 fGammaIso_DR0p1To0p2 = max(min((tmpGammaIso_DR0p1To0p2
468 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p1To0p2,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
469 ,2.5)
470 ,0.0);
471 fGammaIso_DR0p2To0p3 = max(min((tmpGammaIso_DR0p2To0p3
472 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p2To0p3,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
473 ,2.5)
474 ,0.0);
475 fGammaIso_DR0p3To0p4 = max(min((tmpGammaIso_DR0p3To0p4
476 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p3To0p4,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
477 ,2.5)
478 ,0.0);
479 fGammaIso_DR0p4To0p5 = max(min((tmpGammaIso_DR0p4To0p5
480 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p4To0p5,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
481 ,2.5)
482 ,0.0);
483
484
485
486 fNeutralHadronIso_DR0p0To0p1 = max(min((tmpNeutralHadronIso_DR0p0To0p1
487 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p0To0p1,
488 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
489 , 2.5)
490 , 0.0);
491 fNeutralHadronIso_DR0p1To0p2 = max(min((tmpNeutralHadronIso_DR0p1To0p2
492 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p1To0p2,
493 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
494 , 2.5)
495 , 0.0);
496 fNeutralHadronIso_DR0p2To0p3 = max(min((tmpNeutralHadronIso_DR0p2To0p3
497 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p2To0p3,
498 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
499 , 2.5)
500 , 0.0);
501 fNeutralHadronIso_DR0p3To0p4 = max(min((tmpNeutralHadronIso_DR0p3To0p4
502 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p3To0p4,
503 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
504 , 2.5)
505 , 0.0);
506 fNeutralHadronIso_DR0p4To0p5 = max(min((tmpNeutralHadronIso_DR0p4To0p5
507 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p4To0p5,
508 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
509 , 2.5)
510 , 0.0);
511
512
513 double mvaval = muIsoMVA->MVAValue_IsoRings( mu->Pt(),
514 mu->Eta(),
515 fChargedIso_DR0p0To0p1,
516 fChargedIso_DR0p1To0p2,
517 fChargedIso_DR0p2To0p3,
518 fChargedIso_DR0p3To0p4,
519 fChargedIso_DR0p4To0p5,
520 fGammaIso_DR0p0To0p1,
521 fGammaIso_DR0p1To0p2,
522 fGammaIso_DR0p2To0p3,
523 fGammaIso_DR0p3To0p4,
524 fGammaIso_DR0p4To0p5,
525 fNeutralHadronIso_DR0p0To0p1,
526 fNeutralHadronIso_DR0p1To0p2,
527 fNeutralHadronIso_DR0p2To0p3,
528 fNeutralHadronIso_DR0p3To0p4,
529 fNeutralHadronIso_DR0p4To0p5,
530 ctrl.debug);
531
532 SelectionStatus status;
533 bool pass;
534
535 pass = false;
536 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
537 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN0) pass = true;
538 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
539 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN1) pass = true;
540 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
541 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN2) pass = true;
542 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
543 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN3) pass = true;
544 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN4) pass = true;
545 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN5) pass = true;
546 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
547
548 pass = false;
549 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
550 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN0) pass = true;
551 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
552 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN1) pass = true;
553 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
554 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN2) pass = true;
555 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
556 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN3) pass = true;
557 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN4) pass = true;
558 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN5) pass = true;
559 if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
560
561 // pass &= (fChargedIso_DR0p0To0p1 + fChargedIso_DR0p1To0p2 + fChargedIso_DR0p2To0p3 < 0.7);
562
563 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
564 return status;
565
566 }
567
568 //--------------------------------------------------------------------------------------------------
569 void initMuonIsoMVA() {
570 //--------------------------------------------------------------------------------------------------
571 muIsoMVA = new mithep::MuonIDMVA();
572 vector<string> weightFiles;
573 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml");
574 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml");
575 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml");
576 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml");
577 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml");
578 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml");
579 muIsoMVA->Initialize( "MuonIsoMVA",
580 mithep::MuonIDMVA::kIsoRingsV0,
581 kTRUE, weightFiles);
582 }
583
584
585 //--------------------------------------------------------------------------------------------------
586 double muonPFIso04(ControlFlags &ctrl,
587 const mithep::Muon * mu,
588 const mithep::Vertex & vtx,
589 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
590 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
591 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
592 vector<const mithep::Muon*> muonsToVeto,
593 vector<const mithep::Electron*> electronsToVeto)
594 //--------------------------------------------------------------------------------------------------
595 {
596
597 if( ctrl.debug ) {
598 cout << "muonIsoMVASelection :: muons to veto " << endl;
599 for( int i=0; i<muonsToVeto.size(); i++ ) {
600 const mithep::Muon * vmu = muonsToVeto[i];
601 cout << "\tpt: " << vmu->Pt()
602 << "\teta: " << vmu->Eta()
603 << "\tphi: " << vmu->Phi()
604 << endl;
605 }
606 cout << "muonIsoMVASelection :: electrson to veto " << endl;
607 for( int i=0; i<electronsToVeto.size(); i++ ) {
608 const mithep::Electron * vel = electronsToVeto[i];
609 cout << "\tpt: " << vel->Pt()
610 << "\teta: " << vel->Eta()
611 << "\tphi: " << vel->Phi()
612 << endl;
613 }
614 }
615
616 //
617 // final iso
618 //
619 Double_t fChargedIso = 0.0;
620 Double_t fGammaIso = 0.0;
621 Double_t fNeutralHadronIso = 0.0;
622
623 //
624 //Loop over PF Candidates
625 //
626 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
627 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
628
629 Double_t deta = (mu->Eta() - pf->Eta());
630 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
631 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
632 if (dr > 0.4) continue;
633
634 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
635
636 //
637 // Lepton Footprint Removal
638 //
639 Bool_t IsLeptonFootprint = kFALSE;
640 if (dr < 1.0) {
641
642 //
643 // Check for electrons
644 //
645 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
646 const mithep::Electron *tmpele = electronsToVeto[q];
647 // 4l electron
648 if( pf->HasTrackerTrk() ) {
649 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
650 IsLeptonFootprint = kTRUE;
651 }
652 if( pf->HasGsfTrk() ) {
653 if( pf->GsfTrk() == tmpele->GsfTrk() )
654 IsLeptonFootprint = kTRUE;
655 }
656 // PF charged
657 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
658 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
659 IsLeptonFootprint = kTRUE;
660 // PF gamma
661 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
662 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
663 IsLeptonFootprint = kTRUE;
664 } // loop over electrons
665
666 //
667 // Check for muons
668 //
669 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
670 const mithep::Muon *tmpmu = muonsToVeto[q];
671 // 4l muon
672 if( pf->HasTrackerTrk() ) {
673 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
674 IsLeptonFootprint = kTRUE;
675 }
676 // PF charged
677 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
678 IsLeptonFootprint = kTRUE;
679 } // loop over muons
680
681
682 if (IsLeptonFootprint)
683 continue;
684
685 //
686 // Charged Iso Rings
687 //
688 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
689
690 if( dr < 0.01 ) continue; // only for muon iso mva?
691 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
692
693 if( pf->HasTrackerTrk() ) {
694 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
695 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
696 << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
697 << dr << endl;
698 }
699 if( pf->HasGsfTrk() ) {
700 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
701 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
702 << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
703 << dr << endl;
704 }
705
706
707 fChargedIso += pf->Pt();
708 }
709
710 //
711 // Gamma Iso
712 //
713 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
714 fGammaIso += pf->Pt();
715 }
716
717 //
718 // Other Neutral Iso Rings
719 //
720 else {
721 fNeutralHadronIso += pf->Pt();
722 }
723
724 }
725
726 }
727
728 double rho = 0;
729 if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
730 rho = fPUEnergyDensity->At(0)->Rho();
731
732 // WARNING!!!!
733 // hardcode for sync ...
734 EffectiveAreaVersion = muT.kMuEAData2011;
735 // WARNING!!!!
736
737
738 double pfIso = fChargedIso + max(0.0,(fGammaIso + fNeutralHadronIso
739 -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
740 mu->Eta(),EffectiveAreaVersion)));
741
742 return pfIso;
743 }
744
745 //--------------------------------------------------------------------------------------------------
746 SelectionStatus muonReferenceIsoSelection(ControlFlags &ctrl,
747 const mithep::Muon * mu,
748 const mithep::Vertex & vtx,
749 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
750 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
751 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
752 vector<const mithep::Muon*> muonsToVeto,
753 vector<const mithep::Electron*> electronsToVeto)
754 //--------------------------------------------------------------------------------------------------
755 {
756
757 SelectionStatus status;
758
759 double pfIso = muonPFIso04( ctrl, mu, vtx, fPFCandidates, fPUEnergyDensity,
760 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
761 bool pass = false;
762 if( (pfIso/mu->Pt()) < MUON_REFERENCE_PFISO_CUT ) pass = true;
763
764 if( pass ) {
765 status.orStatus(SelectionStatus::LOOSEISO);
766 status.orStatus(SelectionStatus::TIGHTISO);
767 }
768 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
769 return status;
770
771 }
772
773
774
775 //--------------------------------------------------------------------------------------------------
776 SelectionStatus electronIsoMVASelection(ControlFlags &ctrl,
777 const mithep::Electron * ele,
778 const mithep::Vertex & vtx,
779 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
780 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
781 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
782 vector<const mithep::Muon*> muonsToVeto,
783 vector<const mithep::Electron*> electronsToVeto)
784 //--------------------------------------------------------------------------------------------------
785 {
786
787 if( ctrl.debug ) {
788 cout << "electronIsoMVASelection :: muons to veto " << endl;
789 for( int i=0; i<muonsToVeto.size(); i++ ) {
790 const mithep::Muon * vmu = muonsToVeto[i];
791 cout << "\tpt: " << vmu->Pt()
792 << "\teta: " << vmu->Eta()
793 << "\tphi: " << vmu->Phi()
794 << endl;
795 }
796 cout << "electronIsoMVASelection :: electrson to veto " << endl;
797 for( int i=0; i<electronsToVeto.size(); i++ ) {
798 const mithep::Electron * vel = electronsToVeto[i];
799 cout << "\tpt: " << vel->Pt()
800 << "\teta: " << vel->Eta()
801 << "\tphi: " << vel->Phi()
802 << "\ttrk: " << vel->TrackerTrk()
803 << endl;
804 }
805 }
806
807 bool failiso=false;
808
809 //
810 // tmp iso rings
811 //
812 Double_t tmpChargedIso_DR0p0To0p1 = 0;
813 Double_t tmpChargedIso_DR0p1To0p2 = 0;
814 Double_t tmpChargedIso_DR0p2To0p3 = 0;
815 Double_t tmpChargedIso_DR0p3To0p4 = 0;
816 Double_t tmpChargedIso_DR0p4To0p5 = 0;
817 Double_t tmpChargedIso_DR0p5To0p7 = 0;
818
819 Double_t tmpGammaIso_DR0p0To0p1 = 0;
820 Double_t tmpGammaIso_DR0p1To0p2 = 0;
821 Double_t tmpGammaIso_DR0p2To0p3 = 0;
822 Double_t tmpGammaIso_DR0p3To0p4 = 0;
823 Double_t tmpGammaIso_DR0p4To0p5 = 0;
824 Double_t tmpGammaIso_DR0p5To0p7 = 0;
825
826 Double_t tmpNeutralHadronIso_DR0p0To0p1 = 0;
827 Double_t tmpNeutralHadronIso_DR0p1To0p2 = 0;
828 Double_t tmpNeutralHadronIso_DR0p2To0p3 = 0;
829 Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
830 Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
831 Double_t tmpNeutralHadronIso_DR0p5To0p7 = 0;
832
833
834
835 //
836 // final rings for the MVA
837 //
838 Double_t fChargedIso_DR0p0To0p1;
839 Double_t fChargedIso_DR0p1To0p2;
840 Double_t fChargedIso_DR0p2To0p3;
841 Double_t fChargedIso_DR0p3To0p4;
842 Double_t fChargedIso_DR0p4To0p5;
843
844 Double_t fGammaIso_DR0p0To0p1;
845 Double_t fGammaIso_DR0p1To0p2;
846 Double_t fGammaIso_DR0p2To0p3;
847 Double_t fGammaIso_DR0p3To0p4;
848 Double_t fGammaIso_DR0p4To0p5;
849
850 Double_t fNeutralHadronIso_DR0p0To0p1;
851 Double_t fNeutralHadronIso_DR0p1To0p2;
852 Double_t fNeutralHadronIso_DR0p2To0p3;
853 Double_t fNeutralHadronIso_DR0p3To0p4;
854 Double_t fNeutralHadronIso_DR0p4To0p5;
855
856
857 //
858 //Loop over PF Candidates
859 //
860 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
861 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
862 Double_t deta = (ele->Eta() - pf->Eta());
863 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
864 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
865 if (dr >= 0.5) continue;
866 if(ctrl.debug) {
867 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
868 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(vtx);
869 cout << endl;
870 }
871
872
873 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
874 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
875
876
877 //
878 // Lepton Footprint Removal
879 //
880 Bool_t IsLeptonFootprint = kFALSE;
881 if (dr < 1.0) {
882
883 //
884 // Check for electrons
885 //
886 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
887 const mithep::Electron *tmpele = electronsToVeto[q];
888 // 4l electron
889 if( pf->HasTrackerTrk() ) {
890 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
891 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
892 IsLeptonFootprint = kTRUE;
893 }
894 }
895 if( pf->HasGsfTrk() ) {
896 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
897 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
898 IsLeptonFootprint = kTRUE;
899 }
900 }
901 // PF charged
902 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
903 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
904 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
905 IsLeptonFootprint = kTRUE;
906 }
907 // PF gamma
908 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
909 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
910 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
911 IsLeptonFootprint = kTRUE;
912 }
913 } // loop over electrons
914
915 //
916 // Check for muons
917 //
918 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
919 const mithep::Muon *tmpmu = muonsToVeto[q];
920 // 4l muon
921 if( pf->HasTrackerTrk() ) {
922 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
923 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
924 IsLeptonFootprint = kTRUE;
925 }
926 }
927 // PF charged
928 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
929 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
930 IsLeptonFootprint = kTRUE;
931 }
932 } // loop over muons
933
934
935 if (IsLeptonFootprint)
936 continue;
937
938 //
939 // Charged Iso Rings
940 //
941 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
942
943 if( pf->HasTrackerTrk() )
944 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
945 if( pf->HasGsfTrk() )
946 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
947
948 // Veto any PFmuon, or PFEle
949 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
950
951 // Footprint Veto
952 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
953
954 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
955 << "\ttype: " << pf->PFType()
956 << "\ttrk: " << pf->TrackerTrk() << endl;
957
958 if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
959 if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
960 if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
961 if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
962 if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
963 if (dr >= 0.5 && dr < 0.7) tmpChargedIso_DR0p5To0p7 += pf->Pt();
964
965 }
966
967 //
968 // Gamma Iso Rings
969 //
970 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
971
972 if (fabs(ele->SCluster()->Eta()) > 1.479) {
973 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
974 }
975
976 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
977 << dr << endl;
978
979 if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
980 if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
981 if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
982 if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
983 if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
984 if (dr >= 0.5 && dr < 0.7) tmpGammaIso_DR0p5To0p7 += pf->Pt();
985
986 }
987
988 //
989 // Other Neutral Iso Rings
990 //
991 else {
992 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
993 << dr << endl;
994 if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
995 if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
996 if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
997 if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
998 if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
999 if (dr >= 0.5 && dr < 0.7) tmpNeutralHadronIso_DR0p5To0p7 += pf->Pt();
1000 }
1001
1002 }
1003
1004 }
1005
1006 fChargedIso_DR0p0To0p1 = min((tmpChargedIso_DR0p0To0p1)/ele->Pt(), 2.5);
1007 fChargedIso_DR0p1To0p2 = min((tmpChargedIso_DR0p1To0p2)/ele->Pt(), 2.5);
1008 fChargedIso_DR0p2To0p3 = min((tmpChargedIso_DR0p2To0p3)/ele->Pt(), 2.5);
1009 fChargedIso_DR0p3To0p4 = min((tmpChargedIso_DR0p3To0p4)/ele->Pt(), 2.5);
1010 fChargedIso_DR0p4To0p5 = min((tmpChargedIso_DR0p4To0p5)/ele->Pt(), 2.5);
1011
1012 double rho = 0;
1013 if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1014 rho = fPUEnergyDensity->At(0)->Rho();
1015
1016 if( ctrl.debug) {
1017 cout << "RHO: " << rho << endl;
1018 cout << "eta: " << ele->SCluster()->Eta() << endl;
1019 cout << "target: " << EffectiveAreaVersion << endl;
1020 cout << "effA 0-1: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
1021 ele->SCluster()->Eta(),
1022 EffectiveAreaVersion)
1023 << endl;
1024 cout << "effA 1-2: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
1025 ele->SCluster()->Eta(),
1026 EffectiveAreaVersion)
1027 << endl;
1028 cout << "effA 2-3: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
1029 ele->SCluster()->Eta(),
1030 EffectiveAreaVersion)
1031 << endl;
1032 cout << "effA 3-4: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
1033 ele->SCluster()->Eta(),
1034 EffectiveAreaVersion)
1035 << endl;
1036 }
1037
1038 fGammaIso_DR0p0To0p1 = max(min((tmpGammaIso_DR0p0To0p1
1039 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p0To0p1,
1040 ele->SCluster()->Eta(),
1041 EffectiveAreaVersion))/ele->Pt()
1042 ,2.5)
1043 ,0.0);
1044 fGammaIso_DR0p1To0p2 = max(min((tmpGammaIso_DR0p1To0p2
1045 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p1To0p2,
1046 ele->SCluster()->Eta(),
1047 EffectiveAreaVersion))/ele->Pt()
1048 ,2.5)
1049 ,0.0);
1050 fGammaIso_DR0p2To0p3 = max(min((tmpGammaIso_DR0p2To0p3
1051 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p2To0p3,
1052 ele->SCluster()->Eta()
1053 ,EffectiveAreaVersion))/ele->Pt()
1054 ,2.5)
1055 ,0.0);
1056 fGammaIso_DR0p3To0p4 = max(min((tmpGammaIso_DR0p3To0p4
1057 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p3To0p4,
1058 ele->SCluster()->Eta(),
1059 EffectiveAreaVersion))/ele->Pt()
1060 ,2.5)
1061 ,0.0);
1062 fGammaIso_DR0p4To0p5 = max(min((tmpGammaIso_DR0p4To0p5
1063 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p4To0p5,
1064 ele->SCluster()->Eta(),
1065 EffectiveAreaVersion))/ele->Pt()
1066 ,2.5)
1067 ,0.0);
1068
1069
1070 fNeutralHadronIso_DR0p0To0p1 = max(min((tmpNeutralHadronIso_DR0p0To0p1
1071 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
1072 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1073 , 2.5)
1074 , 0.0);
1075 fNeutralHadronIso_DR0p1To0p2 = max(min((tmpNeutralHadronIso_DR0p1To0p2
1076 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
1077 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1078 , 2.5)
1079 , 0.0);
1080 fNeutralHadronIso_DR0p2To0p3 = max(min((tmpNeutralHadronIso_DR0p2To0p3
1081 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
1082 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1083 , 2.5)
1084 , 0.0);
1085 fNeutralHadronIso_DR0p3To0p4 = max(min((tmpNeutralHadronIso_DR0p3To0p4
1086 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
1087 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
1088 , 2.5)
1089 , 0.0);
1090 fNeutralHadronIso_DR0p4To0p5 = max(min((tmpNeutralHadronIso_DR0p4To0p5
1091 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p4To0p5,
1092 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
1093 , 2.5)
1094 , 0.0);
1095
1096 double mvaval = eleIsoMVA->MVAValue_IsoRings( ele->Pt(),
1097 ele->SCluster()->Eta(),
1098 fChargedIso_DR0p0To0p1,
1099 fChargedIso_DR0p1To0p2,
1100 fChargedIso_DR0p2To0p3,
1101 fChargedIso_DR0p3To0p4,
1102 fChargedIso_DR0p4To0p5,
1103 fGammaIso_DR0p0To0p1,
1104 fGammaIso_DR0p1To0p2,
1105 fGammaIso_DR0p2To0p3,
1106 fGammaIso_DR0p3To0p4,
1107 fGammaIso_DR0p4To0p5,
1108 fNeutralHadronIso_DR0p0To0p1,
1109 fNeutralHadronIso_DR0p1To0p2,
1110 fNeutralHadronIso_DR0p2To0p3,
1111 fNeutralHadronIso_DR0p3To0p4,
1112 fNeutralHadronIso_DR0p4To0p5,
1113 ctrl.debug);
1114
1115 SelectionStatus status;
1116 bool pass = false;
1117
1118 Int_t subdet = 0;
1119 if (fabs(ele->SCluster()->Eta()) < 0.8) subdet = 0;
1120 else if (fabs(ele->SCluster()->Eta()) < 1.479) subdet = 1;
1121 else subdet = 2;
1122 Int_t ptBin = 0;
1123 if (ele->Pt() > 10.0) ptBin = 1;
1124
1125 Int_t MVABin = -1;
1126 if (subdet == 0 && ptBin == 0) MVABin = 0;
1127 if (subdet == 1 && ptBin == 0) MVABin = 1;
1128 if (subdet == 2 && ptBin == 0) MVABin = 2;
1129 if (subdet == 0 && ptBin == 1) MVABin = 3;
1130 if (subdet == 1 && ptBin == 1) MVABin = 4;
1131 if (subdet == 2 && ptBin == 1) MVABin = 5;
1132
1133 pass = false;
1134 if( MVABin == 0 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN0 ) pass = true;
1135 if( MVABin == 1 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN1 ) pass = true;
1136 if( MVABin == 2 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN2 ) pass = true;
1137 if( MVABin == 3 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN3 ) pass = true;
1138 if( MVABin == 4 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN4 ) pass = true;
1139 if( MVABin == 5 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN5 ) pass = true;
1140 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
1141
1142 pass = false;
1143 if( MVABin == 0 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN0 ) pass = true;
1144 if( MVABin == 1 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN1 ) pass = true;
1145 if( MVABin == 2 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN2 ) pass = true;
1146 if( MVABin == 3 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN3 ) pass = true;
1147 if( MVABin == 4 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN4 ) pass = true;
1148 if( MVABin == 5 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN5 ) pass = true;
1149 if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
1150
1151 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1152 return status;
1153
1154 }
1155
1156
1157 //--------------------------------------------------------------------------------------------------
1158 void initElectronIsoMVA() {
1159 //--------------------------------------------------------------------------------------------------
1160 eleIsoMVA = new mithep::ElectronIDMVA();
1161 vector<string> weightFiles;
1162 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt5To10.weights.xml");
1163 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt5To10.weights.xml");
1164 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt10ToInf.weights.xml");
1165 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt10ToInf.weights.xml");
1166 eleIsoMVA->Initialize( "ElectronIsoMVA",
1167 mithep::ElectronIDMVA::kIsoRingsV0,
1168 kTRUE, weightFiles);
1169 }
1170
1171
1172
1173 //--------------------------------------------------------------------------------------------------
1174 float electronPFIso04(ControlFlags &ctrl,
1175 const mithep::Electron * ele,
1176 const mithep::Vertex & vtx,
1177 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1178 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1179 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1180 vector<const mithep::Muon*> muonsToVeto,
1181 vector<const mithep::Electron*> electronsToVeto)
1182 //--------------------------------------------------------------------------------------------------
1183 {
1184
1185 if( ctrl.debug ) {
1186 cout << "electronIsoMVASelection :: muons to veto " << endl;
1187 for( int i=0; i<muonsToVeto.size(); i++ ) {
1188 const mithep::Muon * vmu = muonsToVeto[i];
1189 cout << "\tpt: " << vmu->Pt()
1190 << "\teta: " << vmu->Eta()
1191 << "\tphi: " << vmu->Phi()
1192 << endl;
1193 }
1194 cout << "electronIsoMVASelection :: electrson to veto " << endl;
1195 for( int i=0; i<electronsToVeto.size(); i++ ) {
1196 const mithep::Electron * vel = electronsToVeto[i];
1197 cout << "\tpt: " << vel->Pt()
1198 << "\teta: " << vel->Eta()
1199 << "\tphi: " << vel->Phi()
1200 << "\ttrk: " << vel->TrackerTrk()
1201 << endl;
1202 }
1203 }
1204
1205
1206 //
1207 // final iso
1208 //
1209 Double_t fChargedIso = 0.0;
1210 Double_t fGammaIso = 0.0;
1211 Double_t fNeutralHadronIso = 0.0;
1212
1213
1214 //
1215 //Loop over PF Candidates
1216 //
1217 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1218 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1219 Double_t deta = (ele->Eta() - pf->Eta());
1220 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
1221 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
1222 if (dr >= 0.4) continue;
1223 if(ctrl.debug) {
1224 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
1225 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(vtx);
1226 cout << endl;
1227 }
1228
1229
1230 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
1231 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
1232
1233
1234 //
1235 // Lepton Footprint Removal
1236 //
1237 Bool_t IsLeptonFootprint = kFALSE;
1238 if (dr < 1.0) {
1239
1240 //
1241 // Check for electrons
1242 //
1243 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1244 const mithep::Electron *tmpele = electronsToVeto[q];
1245 // 4l electron
1246 if( pf->HasTrackerTrk() ) {
1247 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
1248 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
1249 IsLeptonFootprint = kTRUE;
1250 }
1251 }
1252 if( pf->HasGsfTrk() ) {
1253 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
1254 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
1255 IsLeptonFootprint = kTRUE;
1256 }
1257 }
1258 // PF charged
1259 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
1260 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
1261 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
1262 IsLeptonFootprint = kTRUE;
1263 }
1264 // PF gamma
1265 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
1266 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
1267 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
1268 IsLeptonFootprint = kTRUE;
1269 }
1270 } // loop over electrons
1271
1272 //
1273 // Check for muons
1274 //
1275 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1276 const mithep::Muon *tmpmu = muonsToVeto[q];
1277 // 4l muon
1278 if( pf->HasTrackerTrk() ) {
1279 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
1280 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
1281 IsLeptonFootprint = kTRUE;
1282 }
1283 }
1284 // PF charged
1285 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
1286 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
1287 IsLeptonFootprint = kTRUE;
1288 }
1289 } // loop over muons
1290
1291
1292 if (IsLeptonFootprint)
1293 continue;
1294
1295 //
1296 // Charged Iso Rings
1297 //
1298 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
1299
1300 if( pf->HasTrackerTrk() )
1301 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1302 if( pf->HasGsfTrk() )
1303 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1304
1305 // Veto any PFmuon, or PFEle
1306 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1307
1308 // Footprint Veto
1309 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
1310
1311 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
1312 << "\ttype: " << pf->PFType()
1313 << "\ttrk: " << pf->TrackerTrk() << endl;
1314
1315 fChargedIso += pf->Pt();
1316 }
1317
1318 //
1319 // Gamma Iso
1320 //
1321 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1322
1323 if (fabs(ele->SCluster()->Eta()) > 1.479) {
1324 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
1325 }
1326 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
1327 << dr << endl;
1328 fGammaIso += pf->Pt();
1329 }
1330
1331 //
1332 // Neutral Iso
1333 //
1334 else {
1335 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
1336 << dr << endl;
1337 fNeutralHadronIso += pf->Pt();
1338 }
1339
1340 }
1341
1342 }
1343
1344 double rho = 0;
1345 if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1346 rho = fPUEnergyDensity->At(0)->Rho();
1347
1348 // WARNING!!!!
1349 // hardcode for sync ...
1350 EffectiveAreaVersion = eleT.kEleEAData2011;
1351 // WARNING!!!!
1352
1353
1354 double pfIso = fChargedIso +
1355 max(0.0,fGammaIso -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIso04,
1356 ele->Eta(),EffectiveAreaVersion)) +
1357 max(0.0,fNeutralHadronIso -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIso04,
1358 ele->Eta(),EffectiveAreaVersion)) ;
1359 return pfIso;
1360 }
1361
1362 //--------------------------------------------------------------------------------------------------
1363 SelectionStatus electronReferenceIsoSelection(ControlFlags &ctrl,
1364 const mithep::Electron * ele,
1365 const mithep::Vertex & vtx,
1366 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1367 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1368 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1369 vector<const mithep::Muon*> muonsToVeto,
1370 vector<const mithep::Electron*> electronsToVeto)
1371 //--------------------------------------------------------------------------------------------------
1372 {
1373
1374 SelectionStatus status;
1375
1376 double pfIso = electronPFIso04( ctrl, ele, vtx, fPFCandidates, fPUEnergyDensity,
1377 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
1378 bool pass = false;
1379 if( (pfIso/ele->Pt()) < ELECTRON_REFERENCE_PFISO_CUT ) pass = true;
1380
1381 if( pass ) {
1382 status.orStatus(SelectionStatus::LOOSEISO);
1383 status.orStatus(SelectionStatus::TIGHTISO);
1384 }
1385 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1386 return status;
1387
1388 }