ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/IsolationSelection.cc
Revision: 1.17
Committed: Fri May 11 20:26:43 2012 UTC (13 years ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.16: +3 -1 lines
Log Message:
fix presel check for allArb

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 // global hack to sync
20 double gChargedIso;
21 double gGammaIso;
22 double gNeutralIso;
23
24
25 //--------------------------------------------------------------------------------------------------
26 Float_t computePFMuonIso(const mithep::Muon *muon,
27 const mithep::Vertex & vtx,
28 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
29 const Double_t dRMax)
30 //--------------------------------------------------------------------------------------------------
31 {
32 const Double_t dRMin = 0;
33 const Double_t neuPtMin = 1.0;
34 const Double_t dzMax = 0.1;
35
36 Double_t zLepton = (muon->BestTrk()) ? muon->BestTrk()->DzCorrected(vtx) : 0.0;
37
38 Float_t iso=0;
39 for(UInt_t ipf=0; ipf<fPFCandidates->GetEntries(); ipf++) {
40 const PFCandidate *pfcand = fPFCandidates->At(ipf);
41
42 if(!pfcand->HasTrk() && (pfcand->Pt()<=neuPtMin)) continue; // pT cut on neutral particles
43
44 // exclude THE muon
45 if(pfcand->TrackerTrk() && muon->TrackerTrk() && (pfcand->TrackerTrk()==muon->TrackerTrk())) continue;
46
47 // dz cut
48 Double_t dz = (pfcand->BestTrk()) ? fabs(pfcand->BestTrk()->DzCorrected(vtx) - zLepton) : 0;
49 if(dz >= dzMax) continue;
50
51 // check iso cone
52 Double_t dr = MathUtils::DeltaR(muon->Mom(), pfcand->Mom());
53 if(dr<dRMax && dr>=dRMin)
54 iso += pfcand->Pt();
55 }
56
57 return iso;
58 }
59
60 //--------------------------------------------------------------------------------------------------
61 Float_t computePFEleIso(const mithep::Electron *electron,
62 const mithep::Vertex & fVertex,
63 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
64 const Double_t dRMax)
65 //--------------------------------------------------------------------------------------------------
66 {
67 const Double_t dRMin = 0;
68 const Double_t neuPtMin = 1.0;
69 const Double_t dzMax = 0.1;
70
71 Double_t zLepton = (electron->BestTrk()) ? electron->BestTrk()->DzCorrected(fVertex) : 0.0;
72
73 Float_t iso=0;
74 for(UInt_t ipf=0; ipf<fPFCandidates->GetEntries(); ipf++) {
75 const PFCandidate *pfcand = (PFCandidate*)(fPFCandidates->At(ipf));
76
77 if(!pfcand->HasTrk() && (pfcand->Pt()<=neuPtMin)) continue; // pT cut on neutral particles
78
79 // dz cut
80 Double_t dz = (pfcand->BestTrk()) ? fabs(pfcand->BestTrk()->DzCorrected(fVertex) - zLepton) : 0;
81 if(dz >= dzMax) continue;
82
83 // remove THE electron
84 if(pfcand->TrackerTrk() && electron->TrackerTrk() && (pfcand->TrackerTrk()==electron->TrackerTrk())) continue;
85 if(pfcand->GsfTrk() && electron->GsfTrk() && (pfcand->GsfTrk()==electron->GsfTrk())) continue;
86
87 // check iso cone
88 Double_t dr = MathUtils::DeltaR(electron->Mom(), pfcand->Mom());
89 if(dr<dRMax && dr>=dRMin) {
90 // eta-strip veto for photons
91 if((pfcand->PFType() == PFCandidate::eGamma) && fabs(electron->Eta() - pfcand->Eta()) < 0.025) continue;
92
93 // Inner cone (one tower = dR < 0.07) veto for non-photon neutrals
94 if(!pfcand->HasTrk() && (pfcand->PFType() == PFCandidate::eNeutralHadron) &&
95 (MathUtils::DeltaR(electron->Mom(), pfcand->Mom()) < 0.07)) continue;
96
97 iso += pfcand->Pt();
98 }
99 }
100
101 return iso;
102 };
103
104 //--------------------------------------------------------------------------------------------------
105 bool pairwiseIsoSelection( ControlFlags &ctrl,
106 vector<SimpleLepton> &lepvec,
107 float rho )
108 //--------------------------------------------------------------------------------------------------
109 {
110
111 bool passiso=true;
112
113 for( int i=0; i<lepvec.size(); i++ )
114 {
115
116 if( !(lepvec[i].is4l) ) continue;
117
118 float effArea_ecal_i, effArea_hcal_i;
119 if( lepvec[i].isEB ) {
120 if( lepvec[i].type == 11 ) {
121 effArea_ecal_i = 0.101;
122 effArea_hcal_i = 0.021;
123 } else {
124 effArea_ecal_i = 0.074;
125 effArea_hcal_i = 0.022;
126 }
127 } else {
128 if( lepvec[i].type == 11 ) {
129 effArea_ecal_i = 0.046;
130 effArea_hcal_i = 0.040;
131 } else {
132 effArea_ecal_i = 0.045;
133 effArea_hcal_i = 0.030;
134 }
135 }
136
137 float isoEcal_corr_i = lepvec[i].isoEcal - (effArea_ecal_i*rho);
138 float isoHcal_corr_i = lepvec[i].isoHcal - (effArea_hcal_i*rho);
139
140 for( int j=i+1; j<lepvec.size(); j++ )
141 {
142
143 if( !(lepvec[j].is4l) ) continue;
144
145 float effArea_ecal_j, effArea_hcal_j;
146 if( lepvec[j].isEB ) {
147 if( lepvec[j].type == 11 ) {
148 effArea_ecal_j = 0.101;
149 effArea_hcal_j = 0.021;
150 } else {
151 effArea_ecal_j = 0.074;
152 effArea_hcal_j = 0.022;
153 }
154 } else {
155 if( lepvec[j].type == 11 ) {
156 effArea_ecal_j = 0.046;
157 effArea_hcal_j = 0.040;
158 } else {
159 effArea_ecal_j = 0.045;
160 effArea_hcal_j = 0.030;
161 }
162 }
163
164 float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
165 float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
166 float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec.Pt();
167 float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec.Pt();
168 float comboIso = RIso_i + RIso_j;
169
170 if( comboIso > 0.35 ) {
171 if( ctrl.debug ) cout << "combo failing for indices: " << i << "," << j << endl;
172 passiso = false;
173 return passiso;
174 }
175 }
176 }
177
178 return passiso;
179 }
180
181 //--------------------------------------------------------------------------------------------------
182 SelectionStatus muonIsoSelection(ControlFlags &ctrl,
183 const mithep::Muon * mu,
184 const mithep::Vertex & vtx,
185 const mithep::Array<mithep::PFCandidate> * fPFCandidateCol )
186 //--------------------------------------------------------------------------------------------------
187 {
188 float reliso = computePFMuonIso(mu,vtx,fPFCandidateCol,0.3)/mu->Pt();
189 bool isEB = (fabs(mu->Eta()) < 1.479 ? 1 : 0 );
190 bool failiso = false;
191 if( isEB && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EB_HIGHPT ) {
192 failiso = true;
193 }
194 if( isEB && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EB_LOWPT ) {
195 failiso = true;
196 }
197 if( !(isEB) && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EE_HIGHPT ) {
198 failiso = true;
199 }
200 if( !(isEB) && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EE_LOWPT ) {
201 failiso = true;
202 }
203
204 SelectionStatus status;
205 if( !failiso ) status.setStatus(SelectionStatus::LOOSEISO);
206 if( !failiso ) status.setStatus(SelectionStatus::TIGHTISO);
207 return status;
208
209 };
210
211 //--------------------------------------------------------------------------------------------------
212 SelectionStatus electronIsoSelection(ControlFlags &ctrl,
213 const mithep::Electron * ele,
214 const mithep::Vertex &fVertex,
215 const mithep::Array<mithep::PFCandidate> * fPFCandidates)
216 //--------------------------------------------------------------------------------------------------
217 {
218
219 bool failiso=false;
220
221 float reliso = computePFEleIso(ele,fVertex,fPFCandidates,0.4)/ele->Pt();
222
223 if( ele->IsEB() && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EB_HIGHPT ) {
224 failiso = true;
225 }
226 if( ele->IsEB() && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EB_LOWPT ) {
227 failiso = true;
228 }
229 if(ctrl.debug) cout << "before iso check ..." << endl;
230 if( !(ele->IsEB()) && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EE_HIGHPT ) {
231 if(ctrl.debug) cout << "\tit fails ..." << endl;
232 failiso = true;
233 }
234 if( !(ele->IsEB()) && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EE_LOWPT ) {
235 failiso = true;
236 }
237
238 SelectionStatus status;
239 if( !failiso ) {
240 status.orStatus(SelectionStatus::LOOSEISO);
241 status.orStatus(SelectionStatus::TIGHTISO);
242 }
243 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
244 return status;
245
246 }
247
248
249 bool noIso(ControlFlags &, vector<SimpleLepton> &, float rho) {
250
251 return true;
252 }
253
254 //--------------------------------------------------------------------------------------------------
255 SelectionStatus muonIsoMVASelection(ControlFlags &ctrl,
256 const mithep::Muon * mu,
257 const mithep::Vertex & vtx,
258 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
259 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
260 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
261 vector<const mithep::Muon*> muonsToVeto,
262 vector<const mithep::Electron*> electronsToVeto)
263 //--------------------------------------------------------------------------------------------------
264 {
265
266 if( ctrl.debug ) {
267 cout << "muonIsoMVASelection :: muons to veto " << endl;
268 for( int i=0; i<muonsToVeto.size(); i++ ) {
269 const mithep::Muon * vmu = muonsToVeto[i];
270 cout << "\tpt: " << vmu->Pt()
271 << "\teta: " << vmu->Eta()
272 << "\tphi: " << vmu->Phi()
273 << endl;
274 }
275 cout << "muonIsoMVASelection :: electrson to veto " << endl;
276 for( int i=0; i<electronsToVeto.size(); i++ ) {
277 const mithep::Electron * vel = electronsToVeto[i];
278 cout << "\tpt: " << vel->Pt()
279 << "\teta: " << vel->Eta()
280 << "\tphi: " << vel->Phi()
281 << endl;
282 }
283 }
284 bool failiso=false;
285
286 //
287 // tmp iso rings
288 //
289 Double_t tmpChargedIso_DR0p0To0p1 = 0;
290 Double_t tmpChargedIso_DR0p1To0p2 = 0;
291 Double_t tmpChargedIso_DR0p2To0p3 = 0;
292 Double_t tmpChargedIso_DR0p3To0p4 = 0;
293 Double_t tmpChargedIso_DR0p4To0p5 = 0;
294 Double_t tmpChargedIso_DR0p5To0p7 = 0;
295
296 Double_t tmpGammaIso_DR0p0To0p1 = 0;
297 Double_t tmpGammaIso_DR0p1To0p2 = 0;
298 Double_t tmpGammaIso_DR0p2To0p3 = 0;
299 Double_t tmpGammaIso_DR0p3To0p4 = 0;
300 Double_t tmpGammaIso_DR0p4To0p5 = 0;
301 Double_t tmpGammaIso_DR0p5To0p7 = 0;
302
303 Double_t tmpNeutralHadronIso_DR0p0To0p1 = 0;
304 Double_t tmpNeutralHadronIso_DR0p1To0p2 = 0;
305 Double_t tmpNeutralHadronIso_DR0p2To0p3 = 0;
306 Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
307 Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
308 Double_t tmpNeutralHadronIso_DR0p5To0p7 = 0;
309
310
311
312 //
313 // final rings for the MVA
314 //
315 Double_t fChargedIso_DR0p0To0p1;
316 Double_t fChargedIso_DR0p1To0p2;
317 Double_t fChargedIso_DR0p2To0p3;
318 Double_t fChargedIso_DR0p3To0p4;
319 Double_t fChargedIso_DR0p4To0p5;
320 Double_t fChargedIso_DR0p5To0p7;
321
322 Double_t fGammaIso_DR0p0To0p1;
323 Double_t fGammaIso_DR0p1To0p2;
324 Double_t fGammaIso_DR0p2To0p3;
325 Double_t fGammaIso_DR0p3To0p4;
326 Double_t fGammaIso_DR0p4To0p5;
327 Double_t fGammaIso_DR0p5To0p7;
328
329 Double_t fNeutralHadronIso_DR0p0To0p1;
330 Double_t fNeutralHadronIso_DR0p1To0p2;
331 Double_t fNeutralHadronIso_DR0p2To0p3;
332 Double_t fNeutralHadronIso_DR0p3To0p4;
333 Double_t fNeutralHadronIso_DR0p4To0p5;
334 Double_t fNeutralHadronIso_DR0p5To0p7;
335
336
337 //
338 //Loop over PF Candidates
339 //
340 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
341 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
342
343 Double_t deta = (mu->Eta() - pf->Eta());
344 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
345 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
346 if (dr > 1.0) continue;
347
348 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
349
350 //
351 // Lepton Footprint Removal
352 //
353 Bool_t IsLeptonFootprint = kFALSE;
354 if (dr < 1.0) {
355
356 //
357 // Check for electrons
358 //
359 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
360 const mithep::Electron *tmpele = electronsToVeto[q];
361 // 4l electron
362 if( pf->HasTrackerTrk() ) {
363 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
364 IsLeptonFootprint = kTRUE;
365 }
366 if( pf->HasGsfTrk() ) {
367 if( pf->GsfTrk() == tmpele->GsfTrk() )
368 IsLeptonFootprint = kTRUE;
369 }
370 // PF charged
371 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
372 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
373 IsLeptonFootprint = kTRUE;
374 // PF gamma
375 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
376 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
377 IsLeptonFootprint = kTRUE;
378 } // loop over electrons
379
380 /* KH - commented for sync
381 //
382 // Check for muons
383 //
384 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
385 const mithep::Muon *tmpmu = muonsToVeto[q];
386 // 4l muon
387 if( pf->HasTrackerTrk() ) {
388 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
389 IsLeptonFootprint = kTRUE;
390 }
391 // PF charged
392 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
393 IsLeptonFootprint = kTRUE;
394 } // loop over muons
395 */
396
397 if (IsLeptonFootprint)
398 continue;
399
400 //
401 // Charged Iso Rings
402 //
403 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
404
405 if( dr < 0.01 ) continue; // only for muon iso mva?
406 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
407
408 if( pf->HasTrackerTrk() ) {
409 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
410 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
411 << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
412 << dr << endl;
413 }
414 if( pf->HasGsfTrk() ) {
415 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
416 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
417 << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
418 << dr << endl;
419 }
420
421 // Footprint Veto
422 if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
423 if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
424 if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
425 if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
426 if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
427 if (dr >= 0.5 && dr < 0.7) tmpChargedIso_DR0p5To0p7 += pf->Pt();
428 }
429
430 //
431 // Gamma Iso Rings
432 //
433 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
434 if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
435 if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
436 if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
437 if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
438 if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
439 if (dr >= 0.5 && dr < 0.7) tmpGammaIso_DR0p5To0p7 += pf->Pt();
440 }
441
442 //
443 // Other Neutral Iso Rings
444 //
445 else {
446 if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
447 if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
448 if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
449 if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
450 if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
451 if (dr >= 0.5 && dr < 0.7) tmpNeutralHadronIso_DR0p5To0p7 += pf->Pt();
452 }
453
454 }
455
456 }
457
458 fChargedIso_DR0p0To0p1 = min((tmpChargedIso_DR0p0To0p1)/mu->Pt(), 2.5);
459 fChargedIso_DR0p1To0p2 = min((tmpChargedIso_DR0p1To0p2)/mu->Pt(), 2.5);
460 fChargedIso_DR0p2To0p3 = min((tmpChargedIso_DR0p2To0p3)/mu->Pt(), 2.5);
461 fChargedIso_DR0p3To0p4 = min((tmpChargedIso_DR0p3To0p4)/mu->Pt(), 2.5);
462 fChargedIso_DR0p4To0p5 = min((tmpChargedIso_DR0p4To0p5)/mu->Pt(), 2.5);
463
464
465 double rho = 0;
466 if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
467 rho = fPUEnergyDensity->At(0)->Rho();
468
469
470 fGammaIso_DR0p0To0p1 = max(min((tmpGammaIso_DR0p0To0p1
471 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p0To0p1,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
472 ,2.5)
473 ,0.0);
474 fGammaIso_DR0p1To0p2 = max(min((tmpGammaIso_DR0p1To0p2
475 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p1To0p2,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
476 ,2.5)
477 ,0.0);
478 fGammaIso_DR0p2To0p3 = max(min((tmpGammaIso_DR0p2To0p3
479 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p2To0p3,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
480 ,2.5)
481 ,0.0);
482 fGammaIso_DR0p3To0p4 = max(min((tmpGammaIso_DR0p3To0p4
483 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p3To0p4,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
484 ,2.5)
485 ,0.0);
486 fGammaIso_DR0p4To0p5 = max(min((tmpGammaIso_DR0p4To0p5
487 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p4To0p5,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
488 ,2.5)
489 ,0.0);
490
491
492
493 fNeutralHadronIso_DR0p0To0p1 = max(min((tmpNeutralHadronIso_DR0p0To0p1
494 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p0To0p1,
495 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
496 , 2.5)
497 , 0.0);
498 fNeutralHadronIso_DR0p1To0p2 = max(min((tmpNeutralHadronIso_DR0p1To0p2
499 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p1To0p2,
500 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
501 , 2.5)
502 , 0.0);
503 fNeutralHadronIso_DR0p2To0p3 = max(min((tmpNeutralHadronIso_DR0p2To0p3
504 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p2To0p3,
505 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
506 , 2.5)
507 , 0.0);
508 fNeutralHadronIso_DR0p3To0p4 = max(min((tmpNeutralHadronIso_DR0p3To0p4
509 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p3To0p4,
510 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
511 , 2.5)
512 , 0.0);
513 fNeutralHadronIso_DR0p4To0p5 = max(min((tmpNeutralHadronIso_DR0p4To0p5
514 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p4To0p5,
515 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
516 , 2.5)
517 , 0.0);
518
519
520 double mvaval = muIsoMVA->MVAValue_IsoRings( mu->Pt(),
521 mu->Eta(),
522 fChargedIso_DR0p0To0p1,
523 fChargedIso_DR0p1To0p2,
524 fChargedIso_DR0p2To0p3,
525 fChargedIso_DR0p3To0p4,
526 fChargedIso_DR0p4To0p5,
527 fGammaIso_DR0p0To0p1,
528 fGammaIso_DR0p1To0p2,
529 fGammaIso_DR0p2To0p3,
530 fGammaIso_DR0p3To0p4,
531 fGammaIso_DR0p4To0p5,
532 fNeutralHadronIso_DR0p0To0p1,
533 fNeutralHadronIso_DR0p1To0p2,
534 fNeutralHadronIso_DR0p2To0p3,
535 fNeutralHadronIso_DR0p3To0p4,
536 fNeutralHadronIso_DR0p4To0p5,
537 ctrl.debug);
538
539 SelectionStatus status;
540 bool pass;
541
542 pass = false;
543 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
544 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN0) pass = true;
545 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
546 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN1) pass = true;
547 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
548 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN2) pass = true;
549 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
550 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN3) pass = true;
551 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN4) pass = true;
552 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN5) pass = true;
553 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
554
555 pass = false;
556 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
557 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN0) pass = true;
558 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
559 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN1) pass = true;
560 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
561 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN2) pass = true;
562 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
563 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN3) pass = true;
564 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN4) pass = true;
565 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN5) pass = true;
566 if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
567
568 // pass &= (fChargedIso_DR0p0To0p1 + fChargedIso_DR0p1To0p2 + fChargedIso_DR0p2To0p3 < 0.7);
569
570 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
571 return status;
572
573 }
574
575 //--------------------------------------------------------------------------------------------------
576 void initMuonIsoMVA() {
577 //--------------------------------------------------------------------------------------------------
578 muIsoMVA = new mithep::MuonIDMVA();
579 vector<string> weightFiles;
580 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml");
581 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml");
582 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml");
583 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml");
584 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml");
585 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml");
586 muIsoMVA->Initialize( "MuonIsoMVA",
587 mithep::MuonIDMVA::kIsoRingsV0,
588 kTRUE, weightFiles);
589 }
590
591
592 //--------------------------------------------------------------------------------------------------
593 double muonPFIso04(ControlFlags &ctrl,
594 const mithep::Muon * mu,
595 const mithep::Vertex & vtx,
596 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
597 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
598 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
599 vector<const mithep::Muon*> muonsToVeto,
600 vector<const mithep::Electron*> electronsToVeto)
601 //--------------------------------------------------------------------------------------------------
602 {
603
604 if( ctrl.debug ) {
605 cout << "muonIsoMVASelection :: muons to veto " << endl;
606 for( int i=0; i<muonsToVeto.size(); i++ ) {
607 const mithep::Muon * vmu = muonsToVeto[i];
608 cout << "\tpt: " << vmu->Pt()
609 << "\teta: " << vmu->Eta()
610 << "\tphi: " << vmu->Phi()
611 << endl;
612 }
613 cout << "muonIsoMVASelection :: electrson to veto " << endl;
614 for( int i=0; i<electronsToVeto.size(); i++ ) {
615 const mithep::Electron * vel = electronsToVeto[i];
616 cout << "\tpt: " << vel->Pt()
617 << "\teta: " << vel->Eta()
618 << "\tphi: " << vel->Phi()
619 << endl;
620 }
621 }
622
623 //
624 // final iso
625 //
626 Double_t fChargedIso = 0.0;
627 Double_t fGammaIso = 0.0;
628 Double_t fNeutralHadronIso = 0.0;
629
630 //
631 //Loop over PF Candidates
632 //
633 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
634 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
635
636 Double_t deta = (mu->Eta() - pf->Eta());
637 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
638 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
639 if (dr > 0.4) continue;
640
641 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
642
643 //
644 // Lepton Footprint Removal
645 //
646 Bool_t IsLeptonFootprint = kFALSE;
647 if (dr < 1.0) {
648
649 //
650 // Check for electrons
651 //
652 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
653 const mithep::Electron *tmpele = electronsToVeto[q];
654 // 4l electron
655 if( pf->HasTrackerTrk() ) {
656 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
657 IsLeptonFootprint = kTRUE;
658 }
659 if( pf->HasGsfTrk() ) {
660 if( pf->GsfTrk() == tmpele->GsfTrk() )
661 IsLeptonFootprint = kTRUE;
662 }
663 // PF charged
664 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
665 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
666 IsLeptonFootprint = kTRUE;
667 // PF gamma
668 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
669 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
670 IsLeptonFootprint = kTRUE;
671 } // loop over electrons
672
673 // KH, comment to sync
674 /*
675 //
676 // Check for muons
677 //
678 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
679 const mithep::Muon *tmpmu = muonsToVeto[q];
680 // 4l muon
681 if( pf->HasTrackerTrk() ) {
682 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
683 IsLeptonFootprint = kTRUE;
684 }
685 // PF charged
686 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
687 IsLeptonFootprint = kTRUE;
688 } // loop over muons
689 */
690
691 if (IsLeptonFootprint)
692 continue;
693
694 //
695 // Charged Iso
696 //
697 if (pf->Charge() != 0 ) {
698
699 //if( dr < 0.01 ) continue; // only for muon iso mva?
700 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
701
702 if( pf->HasTrackerTrk() ) {
703 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
704 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
705 << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
706 << dr << endl;
707 }
708 if( pf->HasGsfTrk() ) {
709 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
710 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
711 << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
712 << dr << endl;
713 }
714
715
716 fChargedIso += pf->Pt();
717 }
718
719 //
720 // Gamma Iso
721 //
722 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
723 // KH, add to sync
724 if( pf->Pt() > 0.5 )
725 fGammaIso += pf->Pt();
726 }
727
728 //
729 // Other Neutrals
730 //
731 else {
732 // KH, add to sync
733 if( pf->Pt() > 0.5 )
734 fNeutralHadronIso += pf->Pt();
735 }
736
737 }
738
739 }
740
741 double rho = 0;
742 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
743 // rho = fPUEnergyDensity->At(0)->Rho();
744 if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
745 rho = fPUEnergyDensity->At(0)->RhoLowEta();
746
747 // WARNING!!!!
748 // hardcode for sync ...
749 EffectiveAreaVersion = muT.kMuEAData2011;
750 // WARNING!!!!
751
752
753 double pfIso = fChargedIso + max(0.0,(fGammaIso + fNeutralHadronIso
754 -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
755 mu->Eta(),EffectiveAreaVersion)));
756
757 gChargedIso = fChargedIso;
758 gGammaIso = fGammaIso;
759 gNeutralIso = fNeutralHadronIso;
760 return pfIso;
761 }
762
763
764 //--------------------------------------------------------------------------------------------------
765 // hacked version
766 double muonPFIso04(ControlFlags &ctrl,
767 const mithep::Muon * mu,
768 const mithep::Vertex & vtx,
769 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
770 float rho,
771 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
772 vector<const mithep::Muon*> muonsToVeto,
773 vector<const mithep::Electron*> electronsToVeto)
774 //--------------------------------------------------------------------------------------------------
775 {
776
777 extern double gChargedIso;
778 extern double gGammaIso;
779 extern double gNeutralIso;
780
781 if( ctrl.debug ) {
782 cout << "muonIsoMVASelection :: muons to veto " << endl;
783 for( int i=0; i<muonsToVeto.size(); i++ ) {
784 const mithep::Muon * vmu = muonsToVeto[i];
785 cout << "\tpt: " << vmu->Pt()
786 << "\teta: " << vmu->Eta()
787 << "\tphi: " << vmu->Phi()
788 << endl;
789 }
790 cout << "muonIsoMVASelection :: electrson to veto " << endl;
791 for( int i=0; i<electronsToVeto.size(); i++ ) {
792 const mithep::Electron * vel = electronsToVeto[i];
793 cout << "\tpt: " << vel->Pt()
794 << "\teta: " << vel->Eta()
795 << "\tphi: " << vel->Phi()
796 << endl;
797 }
798 }
799
800 //
801 // final iso
802 //
803 Double_t fChargedIso = 0.0;
804 Double_t fGammaIso = 0.0;
805 Double_t fNeutralHadronIso = 0.0;
806
807 //
808 //Loop over PF Candidates
809 //
810 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
811 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
812
813 Double_t deta = (mu->Eta() - pf->Eta());
814 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
815 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
816 if (dr > 0.4) continue;
817
818 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
819
820 //
821 // Lepton Footprint Removal
822 //
823 Bool_t IsLeptonFootprint = kFALSE;
824 if (dr < 1.0) {
825
826 //
827 // Check for electrons
828 //
829 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
830 const mithep::Electron *tmpele = electronsToVeto[q];
831 // 4l electron
832 if( pf->HasTrackerTrk() ) {
833 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
834 IsLeptonFootprint = kTRUE;
835 }
836 if( pf->HasGsfTrk() ) {
837 if( pf->GsfTrk() == tmpele->GsfTrk() )
838 IsLeptonFootprint = kTRUE;
839 }
840 // PF charged
841 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
842 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
843 IsLeptonFootprint = kTRUE;
844 // PF gamma
845 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
846 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
847 IsLeptonFootprint = kTRUE;
848 } // loop over electrons
849
850 /* KH - comment for sync
851 //
852 // Check for muons
853 //
854 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
855 const mithep::Muon *tmpmu = muonsToVeto[q];
856 // 4l muon
857 if( pf->HasTrackerTrk() ) {
858 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
859 IsLeptonFootprint = kTRUE;
860 }
861 // PF charged
862 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
863 IsLeptonFootprint = kTRUE;
864 } // loop over muons
865 */
866
867 if (IsLeptonFootprint)
868 continue;
869
870 //
871 // Charged Iso
872 //
873 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
874
875 if( dr < 0.01 ) continue; // only for muon iso mva?
876 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
877
878 if( pf->HasTrackerTrk() ) {
879 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
880 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
881 << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
882 << dr << endl;
883 }
884 if( pf->HasGsfTrk() ) {
885 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
886 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
887 << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
888 << dr << endl;
889 }
890
891
892 fChargedIso += pf->Pt();
893 }
894
895 //
896 // Gamma Iso
897 //
898 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
899 fGammaIso += pf->Pt();
900 }
901
902 //
903 // Other Neutrals
904 //
905 else {
906 // KH, add to sync
907 if( pf->Pt() > 0.5 )
908 fNeutralHadronIso += pf->Pt();
909 }
910
911 }
912
913 }
914
915 // double rho = 0;
916 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
917 // rho = fPUEnergyDensity->At(0)->Rho();
918
919 // WARNING!!!!
920 // hardcode for sync ...
921 EffectiveAreaVersion = muT.kMuEAData2011;
922 // WARNING!!!!
923
924
925 double pfIso = fChargedIso + max(0.0,(fGammaIso + fNeutralHadronIso
926 -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
927 mu->Eta(),EffectiveAreaVersion)));
928 gChargedIso = fChargedIso;
929 gGammaIso = fGammaIso;
930 gNeutralIso = fNeutralHadronIso;
931
932 return pfIso;
933 }
934
935
936 //--------------------------------------------------------------------------------------------------
937 SelectionStatus muonReferenceIsoSelection(ControlFlags &ctrl,
938 const mithep::Muon * mu,
939 const mithep::Vertex & vtx,
940 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
941 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
942 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
943 vector<const mithep::Muon*> muonsToVeto,
944 vector<const mithep::Electron*> electronsToVeto)
945 //--------------------------------------------------------------------------------------------------
946 {
947
948 SelectionStatus status;
949
950 double pfIso = muonPFIso04( ctrl, mu, vtx, fPFCandidates, fPUEnergyDensity,
951 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
952 cout << "--------------> setting muon isoPF04 to" << pfIso << endl;
953 status.isoPF04 = pfIso;
954 status.chisoPF04 = gChargedIso;
955 status.gaisoPF04 = gGammaIso;
956 status.neisoPF04 = gNeutralIso;
957
958 bool pass = false;
959 if( (pfIso/mu->Pt()) < MUON_REFERENCE_PFISO_CUT ) pass = true;
960
961 if( pass ) {
962 status.orStatus(SelectionStatus::LOOSEISO);
963 status.orStatus(SelectionStatus::TIGHTISO);
964 }
965 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
966 return status;
967
968 }
969
970
971 //--------------------------------------------------------------------------------------------------
972 // hacked version
973 SelectionStatus muonReferenceIsoSelection(ControlFlags &ctrl,
974 const mithep::Muon * mu,
975 const mithep::Vertex & vtx,
976 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
977 float rho,
978 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
979 vector<const mithep::Muon*> muonsToVeto,
980 vector<const mithep::Electron*> electronsToVeto)
981 //--------------------------------------------------------------------------------------------------
982 {
983
984 SelectionStatus status;
985
986 double pfIso = muonPFIso04( ctrl, mu, vtx, fPFCandidates, rho,
987 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
988 bool pass = false;
989 if( (pfIso/mu->Pt()) < MUON_REFERENCE_PFISO_CUT ) pass = true;
990
991 if( pass ) {
992 status.orStatus(SelectionStatus::LOOSEISO);
993 status.orStatus(SelectionStatus::TIGHTISO);
994 }
995 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
996 return status;
997
998 }
999
1000
1001
1002 //--------------------------------------------------------------------------------------------------
1003 SelectionStatus electronIsoMVASelection(ControlFlags &ctrl,
1004 const mithep::Electron * ele,
1005 const mithep::Vertex & vtx,
1006 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1007 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1008 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1009 vector<const mithep::Muon*> muonsToVeto,
1010 vector<const mithep::Electron*> electronsToVeto)
1011 //--------------------------------------------------------------------------------------------------
1012 {
1013
1014 if( ctrl.debug ) {
1015 cout << "electronIsoMVASelection :: muons to veto " << endl;
1016 for( int i=0; i<muonsToVeto.size(); i++ ) {
1017 const mithep::Muon * vmu = muonsToVeto[i];
1018 cout << "\tpt: " << vmu->Pt()
1019 << "\teta: " << vmu->Eta()
1020 << "\tphi: " << vmu->Phi()
1021 << endl;
1022 }
1023 cout << "electronIsoMVASelection :: electrson to veto " << endl;
1024 for( int i=0; i<electronsToVeto.size(); i++ ) {
1025 const mithep::Electron * vel = electronsToVeto[i];
1026 cout << "\tpt: " << vel->Pt()
1027 << "\teta: " << vel->Eta()
1028 << "\tphi: " << vel->Phi()
1029 << "\ttrk: " << vel->TrackerTrk()
1030 << endl;
1031 }
1032 }
1033
1034 bool failiso=false;
1035
1036 //
1037 // tmp iso rings
1038 //
1039 Double_t tmpChargedIso_DR0p0To0p1 = 0;
1040 Double_t tmpChargedIso_DR0p1To0p2 = 0;
1041 Double_t tmpChargedIso_DR0p2To0p3 = 0;
1042 Double_t tmpChargedIso_DR0p3To0p4 = 0;
1043 Double_t tmpChargedIso_DR0p4To0p5 = 0;
1044 Double_t tmpChargedIso_DR0p5To0p7 = 0;
1045
1046 Double_t tmpGammaIso_DR0p0To0p1 = 0;
1047 Double_t tmpGammaIso_DR0p1To0p2 = 0;
1048 Double_t tmpGammaIso_DR0p2To0p3 = 0;
1049 Double_t tmpGammaIso_DR0p3To0p4 = 0;
1050 Double_t tmpGammaIso_DR0p4To0p5 = 0;
1051 Double_t tmpGammaIso_DR0p5To0p7 = 0;
1052
1053 Double_t tmpNeutralHadronIso_DR0p0To0p1 = 0;
1054 Double_t tmpNeutralHadronIso_DR0p1To0p2 = 0;
1055 Double_t tmpNeutralHadronIso_DR0p2To0p3 = 0;
1056 Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
1057 Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
1058 Double_t tmpNeutralHadronIso_DR0p5To0p7 = 0;
1059
1060
1061
1062 //
1063 // final rings for the MVA
1064 //
1065 Double_t fChargedIso_DR0p0To0p1;
1066 Double_t fChargedIso_DR0p1To0p2;
1067 Double_t fChargedIso_DR0p2To0p3;
1068 Double_t fChargedIso_DR0p3To0p4;
1069 Double_t fChargedIso_DR0p4To0p5;
1070
1071 Double_t fGammaIso_DR0p0To0p1;
1072 Double_t fGammaIso_DR0p1To0p2;
1073 Double_t fGammaIso_DR0p2To0p3;
1074 Double_t fGammaIso_DR0p3To0p4;
1075 Double_t fGammaIso_DR0p4To0p5;
1076
1077 Double_t fNeutralHadronIso_DR0p0To0p1;
1078 Double_t fNeutralHadronIso_DR0p1To0p2;
1079 Double_t fNeutralHadronIso_DR0p2To0p3;
1080 Double_t fNeutralHadronIso_DR0p3To0p4;
1081 Double_t fNeutralHadronIso_DR0p4To0p5;
1082
1083
1084 //
1085 //Loop over PF Candidates
1086 //
1087 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1088 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1089 Double_t deta = (ele->Eta() - pf->Eta());
1090 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
1091 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
1092 if (dr >= 0.5) continue;
1093 if(ctrl.debug) {
1094 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
1095 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(vtx);
1096 cout << endl;
1097 }
1098
1099
1100 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
1101 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
1102
1103
1104 //
1105 // Lepton Footprint Removal
1106 //
1107 Bool_t IsLeptonFootprint = kFALSE;
1108 if (dr < 1.0) {
1109
1110 //
1111 // Check for electrons
1112 //
1113 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1114 const mithep::Electron *tmpele = electronsToVeto[q];
1115 // 4l electron
1116 if( pf->HasTrackerTrk() ) {
1117 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
1118 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
1119 IsLeptonFootprint = kTRUE;
1120 }
1121 }
1122 if( pf->HasGsfTrk() ) {
1123 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
1124 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
1125 IsLeptonFootprint = kTRUE;
1126 }
1127 }
1128 // PF charged
1129 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
1130 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
1131 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
1132 IsLeptonFootprint = kTRUE;
1133 }
1134 // PF gamma
1135 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
1136 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
1137 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
1138 IsLeptonFootprint = kTRUE;
1139 }
1140 } // loop over electrons
1141
1142 /* KH - comment for sync
1143 //
1144 // Check for muons
1145 //
1146 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1147 const mithep::Muon *tmpmu = muonsToVeto[q];
1148 // 4l muon
1149 if( pf->HasTrackerTrk() ) {
1150 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
1151 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
1152 IsLeptonFootprint = kTRUE;
1153 }
1154 }
1155 // PF charged
1156 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
1157 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
1158 IsLeptonFootprint = kTRUE;
1159 }
1160 } // loop over muons
1161 */
1162
1163 if (IsLeptonFootprint)
1164 continue;
1165
1166 //
1167 // Charged Iso Rings
1168 //
1169 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
1170
1171 if( pf->HasTrackerTrk() )
1172 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1173 if( pf->HasGsfTrk() )
1174 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1175
1176 // Veto any PFmuon, or PFEle
1177 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1178
1179 // Footprint Veto
1180 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
1181
1182 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
1183 << "\ttype: " << pf->PFType()
1184 << "\ttrk: " << pf->TrackerTrk() << endl;
1185
1186 if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
1187 if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
1188 if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
1189 if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
1190 if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
1191 if (dr >= 0.5 && dr < 0.7) tmpChargedIso_DR0p5To0p7 += pf->Pt();
1192
1193 }
1194
1195 //
1196 // Gamma Iso Rings
1197 //
1198 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1199
1200 if (fabs(ele->SCluster()->Eta()) > 1.479) {
1201 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
1202 }
1203
1204 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
1205 << dr << endl;
1206
1207 if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
1208 if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
1209 if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
1210 if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
1211 if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
1212 if (dr >= 0.5 && dr < 0.7) tmpGammaIso_DR0p5To0p7 += pf->Pt();
1213
1214 }
1215
1216 //
1217 // Other Neutral Iso Rings
1218 //
1219 else {
1220 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
1221 << dr << endl;
1222 if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
1223 if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
1224 if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
1225 if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
1226 if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
1227 if (dr >= 0.5 && dr < 0.7) tmpNeutralHadronIso_DR0p5To0p7 += pf->Pt();
1228 }
1229
1230 }
1231
1232 }
1233
1234 fChargedIso_DR0p0To0p1 = min((tmpChargedIso_DR0p0To0p1)/ele->Pt(), 2.5);
1235 fChargedIso_DR0p1To0p2 = min((tmpChargedIso_DR0p1To0p2)/ele->Pt(), 2.5);
1236 fChargedIso_DR0p2To0p3 = min((tmpChargedIso_DR0p2To0p3)/ele->Pt(), 2.5);
1237 fChargedIso_DR0p3To0p4 = min((tmpChargedIso_DR0p3To0p4)/ele->Pt(), 2.5);
1238 fChargedIso_DR0p4To0p5 = min((tmpChargedIso_DR0p4To0p5)/ele->Pt(), 2.5);
1239
1240 double rho = 0;
1241 if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1242 rho = fPUEnergyDensity->At(0)->Rho();
1243
1244 if( ctrl.debug) {
1245 cout << "RHO: " << rho << endl;
1246 cout << "eta: " << ele->SCluster()->Eta() << endl;
1247 cout << "target: " << EffectiveAreaVersion << endl;
1248 cout << "effA 0-1: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
1249 ele->SCluster()->Eta(),
1250 EffectiveAreaVersion)
1251 << endl;
1252 cout << "effA 1-2: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
1253 ele->SCluster()->Eta(),
1254 EffectiveAreaVersion)
1255 << endl;
1256 cout << "effA 2-3: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
1257 ele->SCluster()->Eta(),
1258 EffectiveAreaVersion)
1259 << endl;
1260 cout << "effA 3-4: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
1261 ele->SCluster()->Eta(),
1262 EffectiveAreaVersion)
1263 << endl;
1264 }
1265
1266 fGammaIso_DR0p0To0p1 = max(min((tmpGammaIso_DR0p0To0p1
1267 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p0To0p1,
1268 ele->SCluster()->Eta(),
1269 EffectiveAreaVersion))/ele->Pt()
1270 ,2.5)
1271 ,0.0);
1272 fGammaIso_DR0p1To0p2 = max(min((tmpGammaIso_DR0p1To0p2
1273 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p1To0p2,
1274 ele->SCluster()->Eta(),
1275 EffectiveAreaVersion))/ele->Pt()
1276 ,2.5)
1277 ,0.0);
1278 fGammaIso_DR0p2To0p3 = max(min((tmpGammaIso_DR0p2To0p3
1279 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p2To0p3,
1280 ele->SCluster()->Eta()
1281 ,EffectiveAreaVersion))/ele->Pt()
1282 ,2.5)
1283 ,0.0);
1284 fGammaIso_DR0p3To0p4 = max(min((tmpGammaIso_DR0p3To0p4
1285 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p3To0p4,
1286 ele->SCluster()->Eta(),
1287 EffectiveAreaVersion))/ele->Pt()
1288 ,2.5)
1289 ,0.0);
1290 fGammaIso_DR0p4To0p5 = max(min((tmpGammaIso_DR0p4To0p5
1291 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p4To0p5,
1292 ele->SCluster()->Eta(),
1293 EffectiveAreaVersion))/ele->Pt()
1294 ,2.5)
1295 ,0.0);
1296
1297
1298 fNeutralHadronIso_DR0p0To0p1 = max(min((tmpNeutralHadronIso_DR0p0To0p1
1299 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
1300 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1301 , 2.5)
1302 , 0.0);
1303 fNeutralHadronIso_DR0p1To0p2 = max(min((tmpNeutralHadronIso_DR0p1To0p2
1304 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
1305 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1306 , 2.5)
1307 , 0.0);
1308 fNeutralHadronIso_DR0p2To0p3 = max(min((tmpNeutralHadronIso_DR0p2To0p3
1309 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
1310 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1311 , 2.5)
1312 , 0.0);
1313 fNeutralHadronIso_DR0p3To0p4 = max(min((tmpNeutralHadronIso_DR0p3To0p4
1314 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
1315 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
1316 , 2.5)
1317 , 0.0);
1318 fNeutralHadronIso_DR0p4To0p5 = max(min((tmpNeutralHadronIso_DR0p4To0p5
1319 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p4To0p5,
1320 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
1321 , 2.5)
1322 , 0.0);
1323
1324 double mvaval = eleIsoMVA->MVAValue_IsoRings( ele->Pt(),
1325 ele->SCluster()->Eta(),
1326 fChargedIso_DR0p0To0p1,
1327 fChargedIso_DR0p1To0p2,
1328 fChargedIso_DR0p2To0p3,
1329 fChargedIso_DR0p3To0p4,
1330 fChargedIso_DR0p4To0p5,
1331 fGammaIso_DR0p0To0p1,
1332 fGammaIso_DR0p1To0p2,
1333 fGammaIso_DR0p2To0p3,
1334 fGammaIso_DR0p3To0p4,
1335 fGammaIso_DR0p4To0p5,
1336 fNeutralHadronIso_DR0p0To0p1,
1337 fNeutralHadronIso_DR0p1To0p2,
1338 fNeutralHadronIso_DR0p2To0p3,
1339 fNeutralHadronIso_DR0p3To0p4,
1340 fNeutralHadronIso_DR0p4To0p5,
1341 ctrl.debug);
1342
1343 SelectionStatus status;
1344 bool pass = false;
1345
1346 Int_t subdet = 0;
1347 if (fabs(ele->SCluster()->Eta()) < 0.8) subdet = 0;
1348 else if (fabs(ele->SCluster()->Eta()) < 1.479) subdet = 1;
1349 else subdet = 2;
1350 Int_t ptBin = 0;
1351 if (ele->Pt() > 10.0) ptBin = 1;
1352
1353 Int_t MVABin = -1;
1354 if (subdet == 0 && ptBin == 0) MVABin = 0;
1355 if (subdet == 1 && ptBin == 0) MVABin = 1;
1356 if (subdet == 2 && ptBin == 0) MVABin = 2;
1357 if (subdet == 0 && ptBin == 1) MVABin = 3;
1358 if (subdet == 1 && ptBin == 1) MVABin = 4;
1359 if (subdet == 2 && ptBin == 1) MVABin = 5;
1360
1361 pass = false;
1362 if( MVABin == 0 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN0 ) pass = true;
1363 if( MVABin == 1 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN1 ) pass = true;
1364 if( MVABin == 2 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN2 ) pass = true;
1365 if( MVABin == 3 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN3 ) pass = true;
1366 if( MVABin == 4 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN4 ) pass = true;
1367 if( MVABin == 5 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN5 ) pass = true;
1368 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
1369
1370 pass = false;
1371 if( MVABin == 0 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN0 ) pass = true;
1372 if( MVABin == 1 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN1 ) pass = true;
1373 if( MVABin == 2 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN2 ) pass = true;
1374 if( MVABin == 3 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN3 ) pass = true;
1375 if( MVABin == 4 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN4 ) pass = true;
1376 if( MVABin == 5 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN5 ) pass = true;
1377 if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
1378
1379 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1380 return status;
1381
1382 }
1383
1384
1385 //--------------------------------------------------------------------------------------------------
1386 void initElectronIsoMVA() {
1387 //--------------------------------------------------------------------------------------------------
1388 eleIsoMVA = new mithep::ElectronIDMVA();
1389 vector<string> weightFiles;
1390 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt5To10.weights.xml");
1391 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt5To10.weights.xml");
1392 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt10ToInf.weights.xml");
1393 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt10ToInf.weights.xml");
1394 eleIsoMVA->Initialize( "ElectronIsoMVA",
1395 mithep::ElectronIDMVA::kIsoRingsV0,
1396 kTRUE, weightFiles);
1397 }
1398
1399
1400
1401 //--------------------------------------------------------------------------------------------------
1402 float electronPFIso04(ControlFlags &ctrl,
1403 const mithep::Electron * ele,
1404 const mithep::Vertex & vtx,
1405 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1406 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1407 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1408 vector<const mithep::Muon*> muonsToVeto,
1409 vector<const mithep::Electron*> electronsToVeto)
1410 //--------------------------------------------------------------------------------------------------
1411 {
1412
1413 if( ctrl.debug ) {
1414 cout << "electronIsoMVASelection :: muons to veto " << endl;
1415 for( int i=0; i<muonsToVeto.size(); i++ ) {
1416 const mithep::Muon * vmu = muonsToVeto[i];
1417 cout << "\tpt: " << vmu->Pt()
1418 << "\teta: " << vmu->Eta()
1419 << "\tphi: " << vmu->Phi()
1420 << endl;
1421 }
1422 cout << "electronIsoMVASelection :: electrson to veto " << endl;
1423 for( int i=0; i<electronsToVeto.size(); i++ ) {
1424 const mithep::Electron * vel = electronsToVeto[i];
1425 cout << "\tpt: " << vel->Pt()
1426 << "\teta: " << vel->Eta()
1427 << "\tphi: " << vel->Phi()
1428 << "\ttrk: " << vel->TrackerTrk()
1429 << endl;
1430 }
1431 }
1432
1433
1434 //
1435 // final iso
1436 //
1437 Double_t fChargedIso = 0.0;
1438 Double_t fGammaIso = 0.0;
1439 Double_t fNeutralHadronIso = 0.0;
1440
1441
1442 //
1443 //Loop over PF Candidates
1444 //
1445 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1446 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1447 Double_t deta = (ele->Eta() - pf->Eta());
1448 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
1449 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
1450 if (dr >= 0.4) continue;
1451 if(ctrl.debug) {
1452 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
1453 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(vtx);
1454 cout << endl;
1455 }
1456
1457
1458 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
1459 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
1460
1461
1462 //
1463 // Lepton Footprint Removal
1464 //
1465 Bool_t IsLeptonFootprint = kFALSE;
1466 if (dr < 1.0) {
1467
1468 //
1469 // Check for electrons
1470 //
1471 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1472 const mithep::Electron *tmpele = electronsToVeto[q];
1473 // 4l electron
1474 if( pf->HasTrackerTrk() ) {
1475 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
1476 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
1477 IsLeptonFootprint = kTRUE;
1478 }
1479 }
1480 if( pf->HasGsfTrk() ) {
1481 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
1482 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
1483 IsLeptonFootprint = kTRUE;
1484 }
1485 }
1486 // PF charged
1487 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
1488 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
1489 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
1490 IsLeptonFootprint = kTRUE;
1491 }
1492 // PF gamma
1493 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
1494 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
1495 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
1496 IsLeptonFootprint = kTRUE;
1497 }
1498 } // loop over electrons
1499
1500
1501 //
1502 // Check for muons
1503 //
1504 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1505 const mithep::Muon *tmpmu = muonsToVeto[q];
1506 // 4l muon
1507 if( pf->HasTrackerTrk() ) {
1508 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
1509 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
1510 IsLeptonFootprint = kTRUE;
1511 }
1512 }
1513 // PF charged
1514 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
1515 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
1516 IsLeptonFootprint = kTRUE;
1517 }
1518 } // loop over muons
1519
1520
1521 if (IsLeptonFootprint)
1522 continue;
1523
1524 //
1525 // Charged Iso
1526 //
1527 if (pf->Charge() != 0 ) {
1528
1529 if( pf->HasTrackerTrk() )
1530 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1531 if( pf->HasGsfTrk() )
1532 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1533
1534 // Veto any PFmuon, or PFEle
1535 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1536
1537 // Footprint Veto
1538 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
1539
1540 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
1541 << "\ttype: " << pf->PFType()
1542 << "\ttrk: " << pf->TrackerTrk() << endl;
1543
1544 fChargedIso += pf->Pt();
1545 }
1546
1547 //
1548 // Gamma Iso
1549 //
1550 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1551
1552 if (fabs(ele->SCluster()->Eta()) > 1.479) {
1553 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
1554 }
1555 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
1556 << dr << endl;
1557 if( pf->Pt() > 0.5 )
1558 fGammaIso += pf->Pt();
1559 }
1560
1561 //
1562 // Neutral Iso
1563 //
1564 else {
1565 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
1566 << dr << endl;
1567 if( pf->Pt() > 0.5 )
1568 fNeutralHadronIso += pf->Pt();
1569 }
1570
1571 }
1572
1573 }
1574
1575 double rho = 0;
1576 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1577 // rho = fPUEnergyDensity->At(0)->Rho();
1578 if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
1579 rho = fPUEnergyDensity->At(0)->RhoLowEta();
1580
1581 // WARNING!!!!
1582 // hardcode for sync ...
1583 EffectiveAreaVersion = eleT.kEleEAData2011;
1584 // WARNING!!!!
1585
1586
1587 double pfIso = fChargedIso +
1588 max(0.0,fGammaIso -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIso04,
1589 ele->Eta(),EffectiveAreaVersion)) +
1590 max(0.0,fNeutralHadronIso -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIso04,
1591 ele->Eta(),EffectiveAreaVersion)) ;
1592
1593 gChargedIso = fChargedIso;
1594 gGammaIso = fGammaIso;
1595 gNeutralIso = fNeutralHadronIso;
1596
1597 return pfIso;
1598 }
1599
1600 //--------------------------------------------------------------------------------------------------
1601 // hacked version
1602 float electronPFIso04(ControlFlags &ctrl,
1603 const mithep::Electron * ele,
1604 const mithep::Vertex & vtx,
1605 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1606 float rho,
1607 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1608 vector<const mithep::Muon*> muonsToVeto,
1609 vector<const mithep::Electron*> electronsToVeto)
1610 //--------------------------------------------------------------------------------------------------
1611 {
1612
1613 if( ctrl.debug ) {
1614 cout << "electronIsoMVASelection :: muons to veto " << endl;
1615 for( int i=0; i<muonsToVeto.size(); i++ ) {
1616 const mithep::Muon * vmu = muonsToVeto[i];
1617 cout << "\tpt: " << vmu->Pt()
1618 << "\teta: " << vmu->Eta()
1619 << "\tphi: " << vmu->Phi()
1620 << endl;
1621 }
1622 cout << "electronIsoMVASelection :: electrson to veto " << endl;
1623 for( int i=0; i<electronsToVeto.size(); i++ ) {
1624 const mithep::Electron * vel = electronsToVeto[i];
1625 cout << "\tpt: " << vel->Pt()
1626 << "\teta: " << vel->Eta()
1627 << "\tphi: " << vel->Phi()
1628 << "\ttrk: " << vel->TrackerTrk()
1629 << endl;
1630 }
1631 }
1632
1633
1634 //
1635 // final iso
1636 //
1637 Double_t fChargedIso = 0.0;
1638 Double_t fGammaIso = 0.0;
1639 Double_t fNeutralHadronIso = 0.0;
1640
1641
1642 //
1643 //Loop over PF Candidates
1644 //
1645 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1646 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1647 Double_t deta = (ele->Eta() - pf->Eta());
1648 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
1649 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
1650 if (dr >= 0.4) continue;
1651 if(ctrl.debug) {
1652 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
1653 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(vtx);
1654 cout << endl;
1655 }
1656
1657
1658 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
1659 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
1660
1661
1662 //
1663 // Lepton Footprint Removal
1664 //
1665 Bool_t IsLeptonFootprint = kFALSE;
1666 if (dr < 1.0) {
1667
1668 //
1669 // Check for electrons
1670 //
1671 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1672 const mithep::Electron *tmpele = electronsToVeto[q];
1673 // 4l electron
1674 if( pf->HasTrackerTrk() ) {
1675 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
1676 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
1677 IsLeptonFootprint = kTRUE;
1678 }
1679 }
1680 if( pf->HasGsfTrk() ) {
1681 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
1682 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
1683 IsLeptonFootprint = kTRUE;
1684 }
1685 }
1686 // PF charged
1687 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
1688 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
1689 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
1690 IsLeptonFootprint = kTRUE;
1691 }
1692 // PF gamma
1693 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
1694 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
1695 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
1696 IsLeptonFootprint = kTRUE;
1697 }
1698 } // loop over electrons
1699
1700 /* KH - comment for sync
1701 //
1702 // Check for muons
1703 //
1704 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1705 const mithep::Muon *tmpmu = muonsToVeto[q];
1706 // 4l muon
1707 if( pf->HasTrackerTrk() ) {
1708 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
1709 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
1710 IsLeptonFootprint = kTRUE;
1711 }
1712 }
1713 // PF charged
1714 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
1715 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
1716 IsLeptonFootprint = kTRUE;
1717 }
1718 } // loop over muons
1719 */
1720
1721 if (IsLeptonFootprint)
1722 continue;
1723
1724 //
1725 // Charged Iso
1726 //
1727 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
1728
1729 if( pf->HasTrackerTrk() )
1730 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1731 if( pf->HasGsfTrk() )
1732 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1733
1734 // Veto any PFmuon, or PFEle
1735 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1736
1737 // Footprint Veto
1738 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
1739
1740 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
1741 << "\ttype: " << pf->PFType()
1742 << "\ttrk: " << pf->TrackerTrk() << endl;
1743
1744 fChargedIso += pf->Pt();
1745 }
1746
1747 //
1748 // Gamma Iso
1749 //
1750 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1751
1752 if (fabs(ele->SCluster()->Eta()) > 1.479) {
1753 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
1754 }
1755 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
1756 << dr << endl;
1757 fGammaIso += pf->Pt();
1758 }
1759
1760 //
1761 // Neutral Iso
1762 //
1763 else {
1764 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
1765 << dr << endl;
1766 // KH, add to sync
1767 if( pf->Pt() > 0.5 )
1768 fNeutralHadronIso += pf->Pt();
1769 }
1770
1771 }
1772
1773 }
1774
1775 // double rho = 0;
1776 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1777 // rho = fPUEnergyDensity->At(0)->Rho();
1778
1779 // WARNING!!!!
1780 // hardcode for sync ...
1781 EffectiveAreaVersion = eleT.kEleEAData2011;
1782 // WARNING!!!!
1783
1784
1785 double pfIso = fChargedIso +
1786 max(0.0,fGammaIso -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIso04,
1787 ele->Eta(),EffectiveAreaVersion)) +
1788 max(0.0,fNeutralHadronIso -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIso04,
1789 ele->Eta(),EffectiveAreaVersion)) ;
1790 return pfIso;
1791 }
1792
1793
1794 //--------------------------------------------------------------------------------------------------
1795 SelectionStatus electronReferenceIsoSelection(ControlFlags &ctrl,
1796 const mithep::Electron * ele,
1797 const mithep::Vertex & vtx,
1798 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1799 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1800 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1801 vector<const mithep::Muon*> muonsToVeto,
1802 vector<const mithep::Electron*> electronsToVeto)
1803 //--------------------------------------------------------------------------------------------------
1804 {
1805
1806 SelectionStatus status;
1807
1808 double pfIso = electronPFIso04( ctrl, ele, vtx, fPFCandidates, fPUEnergyDensity,
1809 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
1810 cout << "--------------> setting electron isoPF04 to " << pfIso << endl;
1811 status.isoPF04 = pfIso;
1812 status.chisoPF04 = gChargedIso;
1813 status.gaisoPF04 = gGammaIso;
1814 status.neisoPF04 = gNeutralIso;
1815
1816 bool pass = false;
1817 if( (pfIso/ele->Pt()) < ELECTRON_REFERENCE_PFISO_CUT ) pass = true;
1818
1819 if( pass ) {
1820 status.orStatus(SelectionStatus::LOOSEISO);
1821 status.orStatus(SelectionStatus::TIGHTISO);
1822 }
1823 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1824 return status;
1825
1826 }
1827
1828
1829 //--------------------------------------------------------------------------------------------------
1830 // hacked version
1831 SelectionStatus electronReferenceIsoSelection(ControlFlags &ctrl,
1832 const mithep::Electron * ele,
1833 const mithep::Vertex & vtx,
1834 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1835 float rho,
1836 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1837 vector<const mithep::Muon*> muonsToVeto,
1838 vector<const mithep::Electron*> electronsToVeto)
1839 //--------------------------------------------------------------------------------------------------
1840 {
1841
1842 SelectionStatus status;
1843
1844 double pfIso = electronPFIso04( ctrl, ele, vtx, fPFCandidates, rho,
1845 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
1846 bool pass = false;
1847 if( (pfIso/ele->Pt()) < ELECTRON_REFERENCE_PFISO_CUT ) pass = true;
1848
1849 if( pass ) {
1850 status.orStatus(SelectionStatus::LOOSEISO);
1851 status.orStatus(SelectionStatus::TIGHTISO);
1852 }
1853 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1854 return status;
1855
1856 }