ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/IsolationSelection.cc
Revision: 1.18
Committed: Sat May 12 03:00:14 2012 UTC (13 years ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.17: +47 -30 lines
Log Message:
*** empty log message ***

File Contents

# Content
1 #include <math.h>
2
3 #include "IsolationSelection.h"
4 #include "IsolationSelectionDefs.h"
5
6 #include "MathUtils.h"
7 #include "MuonTools.h"
8 #include "MuonIDMVA.h"
9 #include "ElectronTools.h"
10 #include "ElectronIDMVA.h"
11
12 using namespace mithep;
13
14 mithep::MuonIDMVA * muIsoMVA;
15 mithep::MuonTools muT;
16 mithep::ElectronIDMVA * eleIsoMVA;
17 mithep::ElectronTools eleT;
18
19 // 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 // if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
469 // rho = fPUEnergyDensity->At(0)->RhoLowEta();
470
471 // WARNING!!!!
472 // hardcode for sync ...
473 EffectiveAreaVersion = muT.kMuEAData2011;
474 // WARNING!!!!
475
476
477 fGammaIso_DR0p0To0p1 = max(min((tmpGammaIso_DR0p0To0p1
478 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p0To0p1,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
479 ,2.5)
480 ,0.0);
481 fGammaIso_DR0p1To0p2 = max(min((tmpGammaIso_DR0p1To0p2
482 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p1To0p2,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
483 ,2.5)
484 ,0.0);
485 fGammaIso_DR0p2To0p3 = max(min((tmpGammaIso_DR0p2To0p3
486 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p2To0p3,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
487 ,2.5)
488 ,0.0);
489 fGammaIso_DR0p3To0p4 = max(min((tmpGammaIso_DR0p3To0p4
490 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p3To0p4,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
491 ,2.5)
492 ,0.0);
493 fGammaIso_DR0p4To0p5 = max(min((tmpGammaIso_DR0p4To0p5
494 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p4To0p5,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
495 ,2.5)
496 ,0.0);
497
498
499
500 fNeutralHadronIso_DR0p0To0p1 = max(min((tmpNeutralHadronIso_DR0p0To0p1
501 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p0To0p1,
502 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
503 , 2.5)
504 , 0.0);
505 fNeutralHadronIso_DR0p1To0p2 = max(min((tmpNeutralHadronIso_DR0p1To0p2
506 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p1To0p2,
507 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
508 , 2.5)
509 , 0.0);
510 fNeutralHadronIso_DR0p2To0p3 = max(min((tmpNeutralHadronIso_DR0p2To0p3
511 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p2To0p3,
512 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
513 , 2.5)
514 , 0.0);
515 fNeutralHadronIso_DR0p3To0p4 = max(min((tmpNeutralHadronIso_DR0p3To0p4
516 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p3To0p4,
517 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
518 , 2.5)
519 , 0.0);
520 fNeutralHadronIso_DR0p4To0p5 = max(min((tmpNeutralHadronIso_DR0p4To0p5
521 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p4To0p5,
522 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
523 , 2.5)
524 , 0.0);
525
526
527 double mvaval = muIsoMVA->MVAValue_IsoRings( mu->Pt(),
528 mu->Eta(),
529 mu->IsGlobalMuon(),
530 mu->IsTrackerMuon(),
531 fChargedIso_DR0p0To0p1,
532 fChargedIso_DR0p1To0p2,
533 fChargedIso_DR0p2To0p3,
534 fChargedIso_DR0p3To0p4,
535 fChargedIso_DR0p4To0p5,
536 fGammaIso_DR0p0To0p1,
537 fGammaIso_DR0p1To0p2,
538 fGammaIso_DR0p2To0p3,
539 fGammaIso_DR0p3To0p4,
540 fGammaIso_DR0p4To0p5,
541 fNeutralHadronIso_DR0p0To0p1,
542 fNeutralHadronIso_DR0p1To0p2,
543 fNeutralHadronIso_DR0p2To0p3,
544 fNeutralHadronIso_DR0p3To0p4,
545 fNeutralHadronIso_DR0p4To0p5,
546 ctrl.debug);
547
548 SelectionStatus status;
549 bool pass;
550
551 pass = false;
552 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
553 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN0) pass = true;
554 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
555 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN1) pass = true;
556 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
557 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN2) pass = true;
558 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
559 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN3) pass = true;
560 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN4) pass = true;
561 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN5) pass = true;
562 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
563
564 pass = false;
565 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
566 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN0) pass = true;
567 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
568 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN1) pass = true;
569 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
570 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN2) pass = true;
571 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
572 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN3) pass = true;
573 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN4) pass = true;
574 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN5) pass = true;
575 if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
576
577 // pass &= (fChargedIso_DR0p0To0p1 + fChargedIso_DR0p1To0p2 + fChargedIso_DR0p2To0p3 < 0.7);
578
579 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
580 return status;
581
582 }
583
584 //--------------------------------------------------------------------------------------------------
585 void initMuonIsoMVA() {
586 //--------------------------------------------------------------------------------------------------
587 muIsoMVA = new mithep::MuonIDMVA();
588 vector<string> weightFiles;
589 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml");
590 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml");
591 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml");
592 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml");
593 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml");
594 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml");
595 muIsoMVA->Initialize( "MuonIsoMVA",
596 mithep::MuonIDMVA::kIsoRingsV0,
597 kTRUE, weightFiles);
598 }
599
600
601 //--------------------------------------------------------------------------------------------------
602 double muonPFIso04(ControlFlags &ctrl,
603 const mithep::Muon * mu,
604 const mithep::Vertex & vtx,
605 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
606 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
607 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
608 vector<const mithep::Muon*> muonsToVeto,
609 vector<const mithep::Electron*> electronsToVeto)
610 //--------------------------------------------------------------------------------------------------
611 {
612
613 if( ctrl.debug ) {
614 cout << "muonIsoMVASelection :: muons to veto " << endl;
615 for( int i=0; i<muonsToVeto.size(); i++ ) {
616 const mithep::Muon * vmu = muonsToVeto[i];
617 cout << "\tpt: " << vmu->Pt()
618 << "\teta: " << vmu->Eta()
619 << "\tphi: " << vmu->Phi()
620 << endl;
621 }
622 cout << "muonIsoMVASelection :: electrson to veto " << endl;
623 for( int i=0; i<electronsToVeto.size(); i++ ) {
624 const mithep::Electron * vel = electronsToVeto[i];
625 cout << "\tpt: " << vel->Pt()
626 << "\teta: " << vel->Eta()
627 << "\tphi: " << vel->Phi()
628 << endl;
629 }
630 }
631
632 //
633 // final iso
634 //
635 Double_t fChargedIso = 0.0;
636 Double_t fGammaIso = 0.0;
637 Double_t fNeutralHadronIso = 0.0;
638
639 //
640 //Loop over PF Candidates
641 //
642 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
643 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
644
645 Double_t deta = (mu->Eta() - pf->Eta());
646 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
647 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
648 if (dr > 0.4) continue;
649
650 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
651
652 //
653 // Lepton Footprint Removal
654 //
655 Bool_t IsLeptonFootprint = kFALSE;
656 if (dr < 1.0) {
657
658 //
659 // Check for electrons
660 //
661 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
662 const mithep::Electron *tmpele = electronsToVeto[q];
663 // 4l electron
664 if( pf->HasTrackerTrk() ) {
665 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
666 IsLeptonFootprint = kTRUE;
667 }
668 if( pf->HasGsfTrk() ) {
669 if( pf->GsfTrk() == tmpele->GsfTrk() )
670 IsLeptonFootprint = kTRUE;
671 }
672 // PF charged
673 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
674 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
675 IsLeptonFootprint = kTRUE;
676 // PF gamma
677 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
678 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
679 IsLeptonFootprint = kTRUE;
680 } // loop over electrons
681
682 // KH, comment to sync
683 /*
684 //
685 // Check for muons
686 //
687 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
688 const mithep::Muon *tmpmu = muonsToVeto[q];
689 // 4l muon
690 if( pf->HasTrackerTrk() ) {
691 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
692 IsLeptonFootprint = kTRUE;
693 }
694 // PF charged
695 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
696 IsLeptonFootprint = kTRUE;
697 } // loop over muons
698 */
699
700 if (IsLeptonFootprint)
701 continue;
702
703 //
704 // Charged Iso
705 //
706 if (pf->Charge() != 0 ) {
707
708 //if( dr < 0.01 ) continue; // only for muon iso mva?
709 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
710
711 if( pf->HasTrackerTrk() ) {
712 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
713 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
714 << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
715 << dr << endl;
716 }
717 if( pf->HasGsfTrk() ) {
718 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
719 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
720 << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
721 << dr << endl;
722 }
723
724
725 fChargedIso += pf->Pt();
726 }
727
728 //
729 // Gamma Iso
730 //
731 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
732 // KH, add to sync
733 if( pf->Pt() > 0.5 )
734 fGammaIso += pf->Pt();
735 }
736
737 //
738 // Other Neutrals
739 //
740 else {
741 // KH, add to sync
742 if( pf->Pt() > 0.5 )
743 fNeutralHadronIso += pf->Pt();
744 }
745
746 }
747
748 }
749
750 double rho = 0;
751 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
752 // rho = fPUEnergyDensity->At(0)->Rho();
753 if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
754 rho = fPUEnergyDensity->At(0)->RhoLowEta();
755
756 // WARNING!!!!
757 // hardcode for sync ...
758 EffectiveAreaVersion = muT.kMuEAData2011;
759 // WARNING!!!!
760
761
762 double pfIso = fChargedIso + max(0.0,(fGammaIso + fNeutralHadronIso
763 -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
764 mu->Eta(),EffectiveAreaVersion)));
765
766 gChargedIso = fChargedIso;
767 gGammaIso = fGammaIso;
768 gNeutralIso = fNeutralHadronIso;
769 return pfIso;
770 }
771
772
773 //--------------------------------------------------------------------------------------------------
774 // hacked version
775 double muonPFIso04(ControlFlags &ctrl,
776 const mithep::Muon * mu,
777 const mithep::Vertex & vtx,
778 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
779 float rho,
780 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
781 vector<const mithep::Muon*> muonsToVeto,
782 vector<const mithep::Electron*> electronsToVeto)
783 //--------------------------------------------------------------------------------------------------
784 {
785
786 extern double gChargedIso;
787 extern double gGammaIso;
788 extern double gNeutralIso;
789
790 if( ctrl.debug ) {
791 cout << "muonIsoMVASelection :: muons to veto " << endl;
792 for( int i=0; i<muonsToVeto.size(); i++ ) {
793 const mithep::Muon * vmu = muonsToVeto[i];
794 cout << "\tpt: " << vmu->Pt()
795 << "\teta: " << vmu->Eta()
796 << "\tphi: " << vmu->Phi()
797 << endl;
798 }
799 cout << "muonIsoMVASelection :: electrson to veto " << endl;
800 for( int i=0; i<electronsToVeto.size(); i++ ) {
801 const mithep::Electron * vel = electronsToVeto[i];
802 cout << "\tpt: " << vel->Pt()
803 << "\teta: " << vel->Eta()
804 << "\tphi: " << vel->Phi()
805 << endl;
806 }
807 }
808
809 //
810 // final iso
811 //
812 Double_t fChargedIso = 0.0;
813 Double_t fGammaIso = 0.0;
814 Double_t fNeutralHadronIso = 0.0;
815
816 //
817 //Loop over PF Candidates
818 //
819 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
820 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
821
822 Double_t deta = (mu->Eta() - pf->Eta());
823 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
824 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
825 if (dr > 0.4) continue;
826
827 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
828
829 //
830 // Lepton Footprint Removal
831 //
832 Bool_t IsLeptonFootprint = kFALSE;
833 if (dr < 1.0) {
834
835 //
836 // Check for electrons
837 //
838 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
839 const mithep::Electron *tmpele = electronsToVeto[q];
840 // 4l electron
841 if( pf->HasTrackerTrk() ) {
842 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
843 IsLeptonFootprint = kTRUE;
844 }
845 if( pf->HasGsfTrk() ) {
846 if( pf->GsfTrk() == tmpele->GsfTrk() )
847 IsLeptonFootprint = kTRUE;
848 }
849 // PF charged
850 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
851 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
852 IsLeptonFootprint = kTRUE;
853 // PF gamma
854 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
855 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
856 IsLeptonFootprint = kTRUE;
857 } // loop over electrons
858
859 /* KH - comment for sync
860 //
861 // Check for muons
862 //
863 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
864 const mithep::Muon *tmpmu = muonsToVeto[q];
865 // 4l muon
866 if( pf->HasTrackerTrk() ) {
867 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
868 IsLeptonFootprint = kTRUE;
869 }
870 // PF charged
871 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
872 IsLeptonFootprint = kTRUE;
873 } // loop over muons
874 */
875
876 if (IsLeptonFootprint)
877 continue;
878
879 //
880 // Charged Iso
881 //
882 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
883
884 if( dr < 0.01 ) continue; // only for muon iso mva?
885 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
886
887 if( pf->HasTrackerTrk() ) {
888 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
889 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
890 << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
891 << dr << endl;
892 }
893 if( pf->HasGsfTrk() ) {
894 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
895 if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
896 << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
897 << dr << endl;
898 }
899
900
901 fChargedIso += pf->Pt();
902 }
903
904 //
905 // Gamma Iso
906 //
907 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
908 if( pf->Pt() > 0.5 )
909 fGammaIso += pf->Pt();
910 }
911
912 //
913 // Other Neutrals
914 //
915 else {
916 // KH, add to sync
917 if( pf->Pt() > 0.5 )
918 fNeutralHadronIso += pf->Pt();
919 }
920
921 }
922
923 }
924
925 // double rho = 0;
926 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
927 // rho = fPUEnergyDensity->At(0)->Rho();
928
929 // WARNING!!!!
930 // hardcode for sync ...
931 EffectiveAreaVersion = muT.kMuEAData2011;
932 // WARNING!!!!
933
934
935 double pfIso = fChargedIso + max(0.0,(fGammaIso + fNeutralHadronIso
936 -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
937 mu->Eta(),EffectiveAreaVersion)));
938 gChargedIso = fChargedIso;
939 gGammaIso = fGammaIso;
940 gNeutralIso = fNeutralHadronIso;
941
942 return pfIso;
943 }
944
945
946 //--------------------------------------------------------------------------------------------------
947 SelectionStatus muonReferenceIsoSelection(ControlFlags &ctrl,
948 const mithep::Muon * mu,
949 const mithep::Vertex & vtx,
950 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
951 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
952 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
953 vector<const mithep::Muon*> muonsToVeto,
954 vector<const mithep::Electron*> electronsToVeto)
955 //--------------------------------------------------------------------------------------------------
956 {
957
958 SelectionStatus status;
959
960 double pfIso = muonPFIso04( ctrl, mu, vtx, fPFCandidates, fPUEnergyDensity,
961 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
962 cout << "--------------> setting muon isoPF04 to" << pfIso << endl;
963 status.isoPF04 = pfIso;
964 status.chisoPF04 = gChargedIso;
965 status.gaisoPF04 = gGammaIso;
966 status.neisoPF04 = gNeutralIso;
967
968 bool pass = false;
969 if( (pfIso/mu->Pt()) < MUON_REFERENCE_PFISO_CUT ) pass = true;
970
971 if( pass ) {
972 status.orStatus(SelectionStatus::LOOSEISO);
973 status.orStatus(SelectionStatus::TIGHTISO);
974 }
975 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
976 return status;
977
978 }
979
980
981 //--------------------------------------------------------------------------------------------------
982 // hacked version
983 SelectionStatus muonReferenceIsoSelection(ControlFlags &ctrl,
984 const mithep::Muon * mu,
985 const mithep::Vertex & vtx,
986 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
987 float rho,
988 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
989 vector<const mithep::Muon*> muonsToVeto,
990 vector<const mithep::Electron*> electronsToVeto)
991 //--------------------------------------------------------------------------------------------------
992 {
993
994 SelectionStatus status;
995
996 double pfIso = muonPFIso04( ctrl, mu, vtx, fPFCandidates, rho,
997 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
998 bool pass = false;
999 if( (pfIso/mu->Pt()) < MUON_REFERENCE_PFISO_CUT ) pass = true;
1000
1001 if( pass ) {
1002 status.orStatus(SelectionStatus::LOOSEISO);
1003 status.orStatus(SelectionStatus::TIGHTISO);
1004 }
1005 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1006 return status;
1007
1008 }
1009
1010
1011
1012 //--------------------------------------------------------------------------------------------------
1013 SelectionStatus electronIsoMVASelection(ControlFlags &ctrl,
1014 const mithep::Electron * ele,
1015 const mithep::Vertex & vtx,
1016 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1017 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1018 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1019 vector<const mithep::Muon*> muonsToVeto,
1020 vector<const mithep::Electron*> electronsToVeto)
1021 //--------------------------------------------------------------------------------------------------
1022 {
1023
1024 if( ctrl.debug ) {
1025 cout << "electronIsoMVASelection :: muons to veto " << endl;
1026 for( int i=0; i<muonsToVeto.size(); i++ ) {
1027 const mithep::Muon * vmu = muonsToVeto[i];
1028 cout << "\tpt: " << vmu->Pt()
1029 << "\teta: " << vmu->Eta()
1030 << "\tphi: " << vmu->Phi()
1031 << endl;
1032 }
1033 cout << "electronIsoMVASelection :: electrson to veto " << endl;
1034 for( int i=0; i<electronsToVeto.size(); i++ ) {
1035 const mithep::Electron * vel = electronsToVeto[i];
1036 cout << "\tpt: " << vel->Pt()
1037 << "\teta: " << vel->Eta()
1038 << "\tphi: " << vel->Phi()
1039 << "\ttrk: " << vel->TrackerTrk()
1040 << endl;
1041 }
1042 }
1043
1044 bool failiso=false;
1045
1046 //
1047 // tmp iso rings
1048 //
1049 Double_t tmpChargedIso_DR0p0To0p1 = 0;
1050 Double_t tmpChargedIso_DR0p1To0p2 = 0;
1051 Double_t tmpChargedIso_DR0p2To0p3 = 0;
1052 Double_t tmpChargedIso_DR0p3To0p4 = 0;
1053 Double_t tmpChargedIso_DR0p4To0p5 = 0;
1054 Double_t tmpChargedIso_DR0p5To0p7 = 0;
1055
1056 Double_t tmpGammaIso_DR0p0To0p1 = 0;
1057 Double_t tmpGammaIso_DR0p1To0p2 = 0;
1058 Double_t tmpGammaIso_DR0p2To0p3 = 0;
1059 Double_t tmpGammaIso_DR0p3To0p4 = 0;
1060 Double_t tmpGammaIso_DR0p4To0p5 = 0;
1061 Double_t tmpGammaIso_DR0p5To0p7 = 0;
1062
1063 Double_t tmpNeutralHadronIso_DR0p0To0p1 = 0;
1064 Double_t tmpNeutralHadronIso_DR0p1To0p2 = 0;
1065 Double_t tmpNeutralHadronIso_DR0p2To0p3 = 0;
1066 Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
1067 Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
1068 Double_t tmpNeutralHadronIso_DR0p5To0p7 = 0;
1069
1070
1071
1072 //
1073 // final rings for the MVA
1074 //
1075 Double_t fChargedIso_DR0p0To0p1;
1076 Double_t fChargedIso_DR0p1To0p2;
1077 Double_t fChargedIso_DR0p2To0p3;
1078 Double_t fChargedIso_DR0p3To0p4;
1079 Double_t fChargedIso_DR0p4To0p5;
1080
1081 Double_t fGammaIso_DR0p0To0p1;
1082 Double_t fGammaIso_DR0p1To0p2;
1083 Double_t fGammaIso_DR0p2To0p3;
1084 Double_t fGammaIso_DR0p3To0p4;
1085 Double_t fGammaIso_DR0p4To0p5;
1086
1087 Double_t fNeutralHadronIso_DR0p0To0p1;
1088 Double_t fNeutralHadronIso_DR0p1To0p2;
1089 Double_t fNeutralHadronIso_DR0p2To0p3;
1090 Double_t fNeutralHadronIso_DR0p3To0p4;
1091 Double_t fNeutralHadronIso_DR0p4To0p5;
1092
1093
1094 //
1095 //Loop over PF Candidates
1096 //
1097 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1098 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1099 Double_t deta = (ele->Eta() - pf->Eta());
1100 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
1101 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
1102 if (dr >= 0.5) continue;
1103 if(ctrl.debug) {
1104 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
1105 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(vtx);
1106 cout << endl;
1107 }
1108
1109
1110 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
1111 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
1112
1113
1114 //
1115 // Lepton Footprint Removal
1116 //
1117 Bool_t IsLeptonFootprint = kFALSE;
1118 if (dr < 1.0) {
1119
1120 //
1121 // Check for electrons
1122 //
1123 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1124 const mithep::Electron *tmpele = electronsToVeto[q];
1125 // 4l electron
1126 if( pf->HasTrackerTrk() ) {
1127 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
1128 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
1129 IsLeptonFootprint = kTRUE;
1130 }
1131 }
1132 if( pf->HasGsfTrk() ) {
1133 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
1134 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
1135 IsLeptonFootprint = kTRUE;
1136 }
1137 }
1138 // PF charged
1139 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
1140 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
1141 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
1142 IsLeptonFootprint = kTRUE;
1143 }
1144 // PF gamma
1145 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
1146 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
1147 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
1148 IsLeptonFootprint = kTRUE;
1149 }
1150 } // loop over electrons
1151
1152 /* KH - comment for sync
1153 //
1154 // Check for muons
1155 //
1156 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1157 const mithep::Muon *tmpmu = muonsToVeto[q];
1158 // 4l muon
1159 if( pf->HasTrackerTrk() ) {
1160 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
1161 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
1162 IsLeptonFootprint = kTRUE;
1163 }
1164 }
1165 // PF charged
1166 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
1167 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
1168 IsLeptonFootprint = kTRUE;
1169 }
1170 } // loop over muons
1171 */
1172
1173 if (IsLeptonFootprint)
1174 continue;
1175
1176 //
1177 // Charged Iso Rings
1178 //
1179 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
1180
1181 if( pf->HasTrackerTrk() )
1182 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1183 if( pf->HasGsfTrk() )
1184 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1185
1186 // Veto any PFmuon, or PFEle
1187 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1188
1189 // Footprint Veto
1190 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
1191
1192 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
1193 << "\ttype: " << pf->PFType()
1194 << "\ttrk: " << pf->TrackerTrk() << endl;
1195
1196 if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
1197 if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
1198 if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
1199 if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
1200 if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
1201 if (dr >= 0.5 && dr < 0.7) tmpChargedIso_DR0p5To0p7 += pf->Pt();
1202
1203 }
1204
1205 //
1206 // Gamma Iso Rings
1207 //
1208 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1209
1210 if (fabs(ele->SCluster()->Eta()) > 1.479) {
1211 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
1212 }
1213
1214 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
1215 << dr << endl;
1216
1217 if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
1218 if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
1219 if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
1220 if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
1221 if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
1222 if (dr >= 0.5 && dr < 0.7) tmpGammaIso_DR0p5To0p7 += pf->Pt();
1223
1224 }
1225
1226 //
1227 // Other Neutral Iso Rings
1228 //
1229 else {
1230 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
1231 << dr << endl;
1232 if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
1233 if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
1234 if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
1235 if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
1236 if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
1237 if (dr >= 0.5 && dr < 0.7) tmpNeutralHadronIso_DR0p5To0p7 += pf->Pt();
1238 }
1239
1240 }
1241
1242 }
1243
1244 fChargedIso_DR0p0To0p1 = min((tmpChargedIso_DR0p0To0p1)/ele->Pt(), 2.5);
1245 fChargedIso_DR0p1To0p2 = min((tmpChargedIso_DR0p1To0p2)/ele->Pt(), 2.5);
1246 fChargedIso_DR0p2To0p3 = min((tmpChargedIso_DR0p2To0p3)/ele->Pt(), 2.5);
1247 fChargedIso_DR0p3To0p4 = min((tmpChargedIso_DR0p3To0p4)/ele->Pt(), 2.5);
1248 fChargedIso_DR0p4To0p5 = min((tmpChargedIso_DR0p4To0p5)/ele->Pt(), 2.5);
1249
1250 double rho = 0;
1251 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1252 // rho = fPUEnergyDensity->At(0)->Rho();
1253 if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
1254 rho = fPUEnergyDensity->At(0)->RhoLowEta();
1255
1256 // WARNING!!!!
1257 // hardcode for sync ...
1258 EffectiveAreaVersion = eleT.kEleEAData2011;
1259 // WARNING!!!!
1260
1261 if( ctrl.debug) {
1262 cout << "RHO: " << rho << endl;
1263 cout << "eta: " << ele->SCluster()->Eta() << endl;
1264 cout << "target: " << EffectiveAreaVersion << endl;
1265 cout << "effA 0-1: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
1266 ele->SCluster()->Eta(),
1267 EffectiveAreaVersion)
1268 << endl;
1269 cout << "effA 1-2: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
1270 ele->SCluster()->Eta(),
1271 EffectiveAreaVersion)
1272 << endl;
1273 cout << "effA 2-3: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
1274 ele->SCluster()->Eta(),
1275 EffectiveAreaVersion)
1276 << endl;
1277 cout << "effA 3-4: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
1278 ele->SCluster()->Eta(),
1279 EffectiveAreaVersion)
1280 << endl;
1281 }
1282
1283 fGammaIso_DR0p0To0p1 = max(min((tmpGammaIso_DR0p0To0p1
1284 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p0To0p1,
1285 ele->SCluster()->Eta(),
1286 EffectiveAreaVersion))/ele->Pt()
1287 ,2.5)
1288 ,0.0);
1289 fGammaIso_DR0p1To0p2 = max(min((tmpGammaIso_DR0p1To0p2
1290 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p1To0p2,
1291 ele->SCluster()->Eta(),
1292 EffectiveAreaVersion))/ele->Pt()
1293 ,2.5)
1294 ,0.0);
1295 fGammaIso_DR0p2To0p3 = max(min((tmpGammaIso_DR0p2To0p3
1296 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p2To0p3,
1297 ele->SCluster()->Eta()
1298 ,EffectiveAreaVersion))/ele->Pt()
1299 ,2.5)
1300 ,0.0);
1301 fGammaIso_DR0p3To0p4 = max(min((tmpGammaIso_DR0p3To0p4
1302 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p3To0p4,
1303 ele->SCluster()->Eta(),
1304 EffectiveAreaVersion))/ele->Pt()
1305 ,2.5)
1306 ,0.0);
1307 fGammaIso_DR0p4To0p5 = max(min((tmpGammaIso_DR0p4To0p5
1308 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p4To0p5,
1309 ele->SCluster()->Eta(),
1310 EffectiveAreaVersion))/ele->Pt()
1311 ,2.5)
1312 ,0.0);
1313
1314
1315 fNeutralHadronIso_DR0p0To0p1 = max(min((tmpNeutralHadronIso_DR0p0To0p1
1316 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
1317 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1318 , 2.5)
1319 , 0.0);
1320 fNeutralHadronIso_DR0p1To0p2 = max(min((tmpNeutralHadronIso_DR0p1To0p2
1321 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
1322 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1323 , 2.5)
1324 , 0.0);
1325 fNeutralHadronIso_DR0p2To0p3 = max(min((tmpNeutralHadronIso_DR0p2To0p3
1326 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
1327 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1328 , 2.5)
1329 , 0.0);
1330 fNeutralHadronIso_DR0p3To0p4 = max(min((tmpNeutralHadronIso_DR0p3To0p4
1331 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
1332 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
1333 , 2.5)
1334 , 0.0);
1335 fNeutralHadronIso_DR0p4To0p5 = max(min((tmpNeutralHadronIso_DR0p4To0p5
1336 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p4To0p5,
1337 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
1338 , 2.5)
1339 , 0.0);
1340
1341 double mvaval = eleIsoMVA->MVAValue_IsoRings( ele->Pt(),
1342 ele->SCluster()->Eta(),
1343 fChargedIso_DR0p0To0p1,
1344 fChargedIso_DR0p1To0p2,
1345 fChargedIso_DR0p2To0p3,
1346 fChargedIso_DR0p3To0p4,
1347 fChargedIso_DR0p4To0p5,
1348 fGammaIso_DR0p0To0p1,
1349 fGammaIso_DR0p1To0p2,
1350 fGammaIso_DR0p2To0p3,
1351 fGammaIso_DR0p3To0p4,
1352 fGammaIso_DR0p4To0p5,
1353 fNeutralHadronIso_DR0p0To0p1,
1354 fNeutralHadronIso_DR0p1To0p2,
1355 fNeutralHadronIso_DR0p2To0p3,
1356 fNeutralHadronIso_DR0p3To0p4,
1357 fNeutralHadronIso_DR0p4To0p5,
1358 ctrl.debug);
1359
1360 SelectionStatus status;
1361 bool pass = false;
1362
1363 Int_t subdet = 0;
1364 if (fabs(ele->SCluster()->Eta()) < 0.8) subdet = 0;
1365 else if (fabs(ele->SCluster()->Eta()) < 1.479) subdet = 1;
1366 else subdet = 2;
1367 Int_t ptBin = 0;
1368 if (ele->Pt() > 10.0) ptBin = 1;
1369
1370 Int_t MVABin = -1;
1371 if (subdet == 0 && ptBin == 0) MVABin = 0;
1372 if (subdet == 1 && ptBin == 0) MVABin = 1;
1373 if (subdet == 2 && ptBin == 0) MVABin = 2;
1374 if (subdet == 0 && ptBin == 1) MVABin = 3;
1375 if (subdet == 1 && ptBin == 1) MVABin = 4;
1376 if (subdet == 2 && ptBin == 1) MVABin = 5;
1377
1378 pass = false;
1379 if( MVABin == 0 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN0 ) pass = true;
1380 if( MVABin == 1 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN1 ) pass = true;
1381 if( MVABin == 2 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN2 ) pass = true;
1382 if( MVABin == 3 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN3 ) pass = true;
1383 if( MVABin == 4 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN4 ) pass = true;
1384 if( MVABin == 5 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN5 ) pass = true;
1385 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
1386
1387 // pass = false;
1388 // if( MVABin == 0 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN0 ) pass = true;
1389 // if( MVABin == 1 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN1 ) pass = true;
1390 // if( MVABin == 2 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN2 ) pass = true;
1391 // if( MVABin == 3 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN3 ) pass = true;
1392 // if( MVABin == 4 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN4 ) pass = true;
1393 // if( MVABin == 5 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN5 ) pass = true;
1394 // if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
1395
1396 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1397 return status;
1398
1399 }
1400
1401
1402 //--------------------------------------------------------------------------------------------------
1403 void initElectronIsoMVA() {
1404 //--------------------------------------------------------------------------------------------------
1405 eleIsoMVA = new mithep::ElectronIDMVA();
1406 vector<string> weightFiles;
1407 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt5To10.weights.xml");
1408 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt5To10.weights.xml");
1409 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt10ToInf.weights.xml");
1410 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt10ToInf.weights.xml");
1411 eleIsoMVA->Initialize( "ElectronIsoMVA",
1412 mithep::ElectronIDMVA::kIsoRingsV0,
1413 kTRUE, weightFiles);
1414 }
1415
1416
1417
1418 //--------------------------------------------------------------------------------------------------
1419 float electronPFIso04(ControlFlags &ctrl,
1420 const mithep::Electron * ele,
1421 const mithep::Vertex & vtx,
1422 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1423 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1424 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1425 vector<const mithep::Muon*> muonsToVeto,
1426 vector<const mithep::Electron*> electronsToVeto)
1427 //--------------------------------------------------------------------------------------------------
1428 {
1429
1430 if( ctrl.debug ) {
1431 cout << "electronIsoMVASelection :: muons to veto " << endl;
1432 for( int i=0; i<muonsToVeto.size(); i++ ) {
1433 const mithep::Muon * vmu = muonsToVeto[i];
1434 cout << "\tpt: " << vmu->Pt()
1435 << "\teta: " << vmu->Eta()
1436 << "\tphi: " << vmu->Phi()
1437 << endl;
1438 }
1439 cout << "electronIsoMVASelection :: electrson to veto " << endl;
1440 for( int i=0; i<electronsToVeto.size(); i++ ) {
1441 const mithep::Electron * vel = electronsToVeto[i];
1442 cout << "\tpt: " << vel->Pt()
1443 << "\teta: " << vel->Eta()
1444 << "\tphi: " << vel->Phi()
1445 << "\ttrk: " << vel->TrackerTrk()
1446 << endl;
1447 }
1448 }
1449
1450
1451 //
1452 // final iso
1453 //
1454 Double_t fChargedIso = 0.0;
1455 Double_t fGammaIso = 0.0;
1456 Double_t fNeutralHadronIso = 0.0;
1457
1458
1459 //
1460 //Loop over PF Candidates
1461 //
1462 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1463 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1464 Double_t deta = (ele->Eta() - pf->Eta());
1465 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
1466 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
1467 if (dr >= 0.4) continue;
1468 if(ctrl.debug) {
1469 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
1470 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(vtx);
1471 cout << endl;
1472 }
1473
1474
1475 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
1476 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
1477
1478
1479 //
1480 // Lepton Footprint Removal
1481 //
1482 Bool_t IsLeptonFootprint = kFALSE;
1483 if (dr < 1.0) {
1484
1485 //
1486 // Check for electrons
1487 //
1488 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1489 const mithep::Electron *tmpele = electronsToVeto[q];
1490 // 4l electron
1491 if( pf->HasTrackerTrk() ) {
1492 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
1493 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
1494 IsLeptonFootprint = kTRUE;
1495 }
1496 }
1497 if( pf->HasGsfTrk() ) {
1498 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
1499 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
1500 IsLeptonFootprint = kTRUE;
1501 }
1502 }
1503 // PF charged
1504 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
1505 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
1506 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
1507 IsLeptonFootprint = kTRUE;
1508 }
1509 // PF gamma
1510 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
1511 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
1512 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
1513 IsLeptonFootprint = kTRUE;
1514 }
1515 } // loop over electrons
1516
1517
1518 //
1519 // Check for muons
1520 //
1521 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1522 const mithep::Muon *tmpmu = muonsToVeto[q];
1523 // 4l muon
1524 if( pf->HasTrackerTrk() ) {
1525 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
1526 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
1527 IsLeptonFootprint = kTRUE;
1528 }
1529 }
1530 // PF charged
1531 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
1532 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
1533 IsLeptonFootprint = kTRUE;
1534 }
1535 } // loop over muons
1536
1537
1538 if (IsLeptonFootprint)
1539 continue;
1540
1541 //
1542 // Charged Iso
1543 //
1544 if (pf->Charge() != 0 ) {
1545
1546 if( pf->HasTrackerTrk() )
1547 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1548 if( pf->HasGsfTrk() )
1549 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1550
1551 // Veto any PFmuon, or PFEle
1552 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1553
1554 // Footprint Veto
1555 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
1556
1557 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
1558 << "\ttype: " << pf->PFType()
1559 << "\ttrk: " << pf->TrackerTrk() << endl;
1560
1561 fChargedIso += pf->Pt();
1562 }
1563
1564 //
1565 // Gamma Iso
1566 //
1567 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1568
1569 if (fabs(ele->SCluster()->Eta()) > 1.479) {
1570 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
1571 }
1572 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
1573 << dr << endl;
1574 if( pf->Pt() > 0.5 )
1575 fGammaIso += pf->Pt();
1576 }
1577
1578 //
1579 // Neutral Iso
1580 //
1581 else {
1582 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
1583 << dr << endl;
1584 if( pf->Pt() > 0.5 )
1585 fNeutralHadronIso += pf->Pt();
1586 }
1587
1588 }
1589
1590 }
1591
1592 double rho = 0;
1593 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1594 // rho = fPUEnergyDensity->At(0)->Rho();
1595 if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
1596 rho = fPUEnergyDensity->At(0)->RhoLowEta();
1597
1598 // WARNING!!!!
1599 // hardcode for sync ...
1600 EffectiveAreaVersion = eleT.kEleEAData2011;
1601 // WARNING!!!!
1602
1603
1604 double pfIso = fChargedIso +
1605 max(0.0,fGammaIso -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIso04,
1606 ele->Eta(),EffectiveAreaVersion)) +
1607 max(0.0,fNeutralHadronIso -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIso04,
1608 ele->Eta(),EffectiveAreaVersion)) ;
1609
1610 gChargedIso = fChargedIso;
1611 gGammaIso = fGammaIso;
1612 gNeutralIso = fNeutralHadronIso;
1613
1614 return pfIso;
1615 }
1616
1617 //--------------------------------------------------------------------------------------------------
1618 // hacked version
1619 float electronPFIso04(ControlFlags &ctrl,
1620 const mithep::Electron * ele,
1621 const mithep::Vertex & vtx,
1622 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1623 float rho,
1624 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1625 vector<const mithep::Muon*> muonsToVeto,
1626 vector<const mithep::Electron*> electronsToVeto)
1627 //--------------------------------------------------------------------------------------------------
1628 {
1629
1630 if( ctrl.debug ) {
1631 cout << "electronIsoMVASelection :: muons to veto " << endl;
1632 for( int i=0; i<muonsToVeto.size(); i++ ) {
1633 const mithep::Muon * vmu = muonsToVeto[i];
1634 cout << "\tpt: " << vmu->Pt()
1635 << "\teta: " << vmu->Eta()
1636 << "\tphi: " << vmu->Phi()
1637 << endl;
1638 }
1639 cout << "electronIsoMVASelection :: electrson to veto " << endl;
1640 for( int i=0; i<electronsToVeto.size(); i++ ) {
1641 const mithep::Electron * vel = electronsToVeto[i];
1642 cout << "\tpt: " << vel->Pt()
1643 << "\teta: " << vel->Eta()
1644 << "\tphi: " << vel->Phi()
1645 << "\ttrk: " << vel->TrackerTrk()
1646 << endl;
1647 }
1648 }
1649
1650
1651 //
1652 // final iso
1653 //
1654 Double_t fChargedIso = 0.0;
1655 Double_t fGammaIso = 0.0;
1656 Double_t fNeutralHadronIso = 0.0;
1657
1658
1659 //
1660 //Loop over PF Candidates
1661 //
1662 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1663 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1664 Double_t deta = (ele->Eta() - pf->Eta());
1665 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
1666 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
1667 if (dr >= 0.4) continue;
1668 if(ctrl.debug) {
1669 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
1670 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(vtx);
1671 cout << endl;
1672 }
1673
1674
1675 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
1676 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
1677
1678
1679 //
1680 // Lepton Footprint Removal
1681 //
1682 Bool_t IsLeptonFootprint = kFALSE;
1683 if (dr < 1.0) {
1684
1685 //
1686 // Check for electrons
1687 //
1688 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1689 const mithep::Electron *tmpele = electronsToVeto[q];
1690 // 4l electron
1691 if( pf->HasTrackerTrk() ) {
1692 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
1693 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
1694 IsLeptonFootprint = kTRUE;
1695 }
1696 }
1697 if( pf->HasGsfTrk() ) {
1698 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
1699 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
1700 IsLeptonFootprint = kTRUE;
1701 }
1702 }
1703 // PF charged
1704 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
1705 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
1706 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
1707 IsLeptonFootprint = kTRUE;
1708 }
1709 // PF gamma
1710 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
1711 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
1712 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
1713 IsLeptonFootprint = kTRUE;
1714 }
1715 } // loop over electrons
1716
1717 /* KH - comment for sync
1718 //
1719 // Check for muons
1720 //
1721 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1722 const mithep::Muon *tmpmu = muonsToVeto[q];
1723 // 4l muon
1724 if( pf->HasTrackerTrk() ) {
1725 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
1726 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
1727 IsLeptonFootprint = kTRUE;
1728 }
1729 }
1730 // PF charged
1731 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
1732 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
1733 IsLeptonFootprint = kTRUE;
1734 }
1735 } // loop over muons
1736 */
1737
1738 if (IsLeptonFootprint)
1739 continue;
1740
1741 //
1742 // Charged Iso
1743 //
1744 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
1745
1746 if( pf->HasTrackerTrk() )
1747 if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1748 if( pf->HasGsfTrk() )
1749 if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1750
1751 // Veto any PFmuon, or PFEle
1752 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1753
1754 // Footprint Veto
1755 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
1756
1757 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
1758 << "\ttype: " << pf->PFType()
1759 << "\ttrk: " << pf->TrackerTrk() << endl;
1760
1761 fChargedIso += pf->Pt();
1762 }
1763
1764 //
1765 // Gamma Iso
1766 //
1767 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1768
1769 if (fabs(ele->SCluster()->Eta()) > 1.479) {
1770 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
1771 }
1772 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
1773 << dr << endl;
1774 fGammaIso += pf->Pt();
1775 }
1776
1777 //
1778 // Neutral Iso
1779 //
1780 else {
1781 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
1782 << dr << endl;
1783 // KH, add to sync
1784 if( pf->Pt() > 0.5 )
1785 fNeutralHadronIso += pf->Pt();
1786 }
1787
1788 }
1789
1790 }
1791
1792 // double rho = 0;
1793 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1794 // rho = fPUEnergyDensity->At(0)->Rho();
1795
1796 // WARNING!!!!
1797 // hardcode for sync ...
1798 EffectiveAreaVersion = eleT.kEleEAData2011;
1799 // WARNING!!!!
1800
1801
1802 double pfIso = fChargedIso +
1803 max(0.0,fGammaIso -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIso04,
1804 ele->Eta(),EffectiveAreaVersion)) +
1805 max(0.0,fNeutralHadronIso -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIso04,
1806 ele->Eta(),EffectiveAreaVersion)) ;
1807 return pfIso;
1808 }
1809
1810
1811 //--------------------------------------------------------------------------------------------------
1812 SelectionStatus electronReferenceIsoSelection(ControlFlags &ctrl,
1813 const mithep::Electron * ele,
1814 const mithep::Vertex & vtx,
1815 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1816 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1817 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1818 vector<const mithep::Muon*> muonsToVeto,
1819 vector<const mithep::Electron*> electronsToVeto)
1820 //--------------------------------------------------------------------------------------------------
1821 {
1822
1823 SelectionStatus status;
1824
1825 double pfIso = electronPFIso04( ctrl, ele, vtx, fPFCandidates, fPUEnergyDensity,
1826 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
1827 cout << "--------------> setting electron isoPF04 to " << pfIso << endl;
1828 status.isoPF04 = pfIso;
1829 status.chisoPF04 = gChargedIso;
1830 status.gaisoPF04 = gGammaIso;
1831 status.neisoPF04 = gNeutralIso;
1832
1833 bool pass = false;
1834 if( (pfIso/ele->Pt()) < ELECTRON_REFERENCE_PFISO_CUT ) pass = true;
1835
1836 if( pass ) {
1837 status.orStatus(SelectionStatus::LOOSEISO);
1838 status.orStatus(SelectionStatus::TIGHTISO);
1839 }
1840 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1841 return status;
1842
1843 }
1844
1845
1846 //--------------------------------------------------------------------------------------------------
1847 // hacked version
1848 SelectionStatus electronReferenceIsoSelection(ControlFlags &ctrl,
1849 const mithep::Electron * ele,
1850 const mithep::Vertex & vtx,
1851 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1852 float rho,
1853 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1854 vector<const mithep::Muon*> muonsToVeto,
1855 vector<const mithep::Electron*> electronsToVeto)
1856 //--------------------------------------------------------------------------------------------------
1857 {
1858
1859 SelectionStatus status;
1860
1861 double pfIso = electronPFIso04( ctrl, ele, vtx, fPFCandidates, rho,
1862 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
1863 bool pass = false;
1864 if( (pfIso/ele->Pt()) < ELECTRON_REFERENCE_PFISO_CUT ) pass = true;
1865
1866 if( pass ) {
1867 status.orStatus(SelectionStatus::LOOSEISO);
1868 status.orStatus(SelectionStatus::TIGHTISO);
1869 }
1870 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1871 return status;
1872
1873 }