ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/IsolationSelection.cc
Revision: 1.24
Committed: Wed May 23 22:06:39 2012 UTC (12 years, 11 months ago) by anlevin
Content type: text/plain
Branch: MAIN
Changes since 1.23: +10 -7 lines
Log Message: 
made changes for sync on the summer 12 file

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 extern vector<bool> PFnoPUflag;
25
26 //--------------------------------------------------------------------------------------------------
27 Float_t computePFMuonIso(const mithep::Muon *muon,
28 const mithep::Vertex * vtx,
29 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
30 const Double_t dRMax)
31 //--------------------------------------------------------------------------------------------------
32 {
33 const Double_t dRMin = 0;
34 const Double_t neuPtMin = 1.0;
35 const Double_t dzMax = 0.1;
36
37 Double_t zLepton = (muon->BestTrk()) ? muon->BestTrk()->DzCorrected(*vtx) : 0.0;
38
39 Float_t iso=0;
40 for(UInt_t ipf=0; ipf<fPFCandidates->GetEntries(); ipf++) {
41 const PFCandidate *pfcand = fPFCandidates->At(ipf);
42
43 if(!pfcand->HasTrk() && (pfcand->Pt()<=neuPtMin)) continue; // pT cut on neutral particles
44
45 // exclude THE muon
46 if(pfcand->TrackerTrk() && muon->TrackerTrk() && (pfcand->TrackerTrk()==muon->TrackerTrk())) continue;
47
48 // dz cut
49 Double_t dz = (pfcand->BestTrk()) ? fabs(pfcand->BestTrk()->DzCorrected(*vtx) - zLepton) : 0;
50 if(dz >= dzMax) continue;
51
52 // check iso cone
53 Double_t dr = MathUtils::DeltaR(muon->Mom(), pfcand->Mom());
54 if(dr<dRMax && dr>=dRMin)
55 iso += pfcand->Pt();
56 }
57
58 return iso;
59 }
60
61 //--------------------------------------------------------------------------------------------------
62 Float_t computePFEleIso(const mithep::Electron *electron,
63 const mithep::Vertex * fVertex,
64 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
65 const Double_t dRMax)
66 //--------------------------------------------------------------------------------------------------
67 {
68 const Double_t dRMin = 0;
69 const Double_t neuPtMin = 1.0;
70 const Double_t dzMax = 0.1;
71
72 Double_t zLepton = (electron->BestTrk()) ? electron->BestTrk()->DzCorrected(*fVertex) : 0.0;
73
74 Float_t iso=0;
75 for(UInt_t ipf=0; ipf<fPFCandidates->GetEntries(); ipf++) {
76 const PFCandidate *pfcand = (PFCandidate*)(fPFCandidates->At(ipf));
77
78 if(!pfcand->HasTrk() && (pfcand->Pt()<=neuPtMin)) continue; // pT cut on neutral particles
79
80 // dz cut
81 Double_t dz = (pfcand->BestTrk()) ? fabs(pfcand->BestTrk()->DzCorrected(*fVertex) - zLepton) : 0;
82 if(dz >= dzMax) continue;
83
84 // remove THE electron
85 if(pfcand->TrackerTrk() && electron->TrackerTrk() && (pfcand->TrackerTrk()==electron->TrackerTrk())) continue;
86 if(pfcand->GsfTrk() && electron->GsfTrk() && (pfcand->GsfTrk()==electron->GsfTrk())) continue;
87
88 // check iso cone
89 Double_t dr = MathUtils::DeltaR(electron->Mom(), pfcand->Mom());
90 if(dr<dRMax && dr>=dRMin) {
91 // eta-strip veto for photons
92 if((pfcand->PFType() == PFCandidate::eGamma) && fabs(electron->Eta() - pfcand->Eta()) < 0.025) continue;
93
94 // Inner cone (one tower = dR < 0.07) veto for non-photon neutrals
95 if(!pfcand->HasTrk() && (pfcand->PFType() == PFCandidate::eNeutralHadron) &&
96 (MathUtils::DeltaR(electron->Mom(), pfcand->Mom()) < 0.07)) continue;
97
98 iso += pfcand->Pt();
99 }
100 }
101
102 return iso;
103 };
104
105 //--------------------------------------------------------------------------------------------------
106 bool pairwiseIsoSelection( ControlFlags &ctrl,
107 vector<SimpleLepton> &lepvec,
108 float rho )
109 //--------------------------------------------------------------------------------------------------
110 {
111
112 bool passiso=true;
113
114 for( int i=0; i<lepvec.size(); i++ )
115 {
116
117 if( !(lepvec[i].is4l) ) continue;
118
119 float effArea_ecal_i, effArea_hcal_i;
120 if( lepvec[i].isEB ) {
121 if( lepvec[i].type == 11 ) {
122 effArea_ecal_i = 0.101;
123 effArea_hcal_i = 0.021;
124 } else {
125 effArea_ecal_i = 0.074;
126 effArea_hcal_i = 0.022;
127 }
128 } else {
129 if( lepvec[i].type == 11 ) {
130 effArea_ecal_i = 0.046;
131 effArea_hcal_i = 0.040;
132 } else {
133 effArea_ecal_i = 0.045;
134 effArea_hcal_i = 0.030;
135 }
136 }
137
138 float isoEcal_corr_i = lepvec[i].isoEcal - (effArea_ecal_i*rho);
139 float isoHcal_corr_i = lepvec[i].isoHcal - (effArea_hcal_i*rho);
140
141 for( int j=i+1; j<lepvec.size(); j++ )
142 {
143
144 if( !(lepvec[j].is4l) ) continue;
145
146 float effArea_ecal_j, effArea_hcal_j;
147 if( lepvec[j].isEB ) {
148 if( lepvec[j].type == 11 ) {
149 effArea_ecal_j = 0.101;
150 effArea_hcal_j = 0.021;
151 } else {
152 effArea_ecal_j = 0.074;
153 effArea_hcal_j = 0.022;
154 }
155 } else {
156 if( lepvec[j].type == 11 ) {
157 effArea_ecal_j = 0.046;
158 effArea_hcal_j = 0.040;
159 } else {
160 effArea_ecal_j = 0.045;
161 effArea_hcal_j = 0.030;
162 }
163 }
164
165 float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
166 float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
167 float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec.Pt();
168 float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec.Pt();
169 float comboIso = RIso_i + RIso_j;
170
171 if( comboIso > 0.35 ) {
172 if( ctrl.debug ) cout << "combo failing for indices: " << i << "," << j << endl;
173 passiso = false;
174 return passiso;
175 }
176 }
177 }
178
179 return passiso;
180 }
181
182 //--------------------------------------------------------------------------------------------------
183 SelectionStatus muonIsoSelection(ControlFlags &ctrl,
184 const mithep::Muon * mu,
185 const mithep::Vertex * vtx,
186 const mithep::Array<mithep::PFCandidate> * fPFCandidateCol )
187 //--------------------------------------------------------------------------------------------------
188 {
189 float reliso = computePFMuonIso(mu,vtx,fPFCandidateCol,0.3)/mu->Pt();
190 bool isEB = (fabs(mu->Eta()) < 1.479 ? 1 : 0 );
191 bool failiso = false;
192 if( isEB && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EB_HIGHPT ) {
193 failiso = true;
194 }
195 if( isEB && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EB_LOWPT ) {
196 failiso = true;
197 }
198 if( !(isEB) && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EE_HIGHPT ) {
199 failiso = true;
200 }
201 if( !(isEB) && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EE_LOWPT ) {
202 failiso = true;
203 }
204
205 SelectionStatus status;
206 if( !failiso ) status.setStatus(SelectionStatus::LOOSEISO);
207 if( !failiso ) status.setStatus(SelectionStatus::TIGHTISO);
208 return status;
209
210 };
211
212 //--------------------------------------------------------------------------------------------------
213 SelectionStatus electronIsoSelection(ControlFlags &ctrl,
214 const mithep::Electron * ele,
215 const mithep::Vertex *fVertex,
216 const mithep::Array<mithep::PFCandidate> * fPFCandidates)
217 //--------------------------------------------------------------------------------------------------
218 {
219
220 bool failiso=false;
221
222 float reliso = computePFEleIso(ele,fVertex,fPFCandidates,0.4)/ele->Pt();
223
224 if( ele->IsEB() && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EB_HIGHPT ) {
225 failiso = true;
226 }
227 if( ele->IsEB() && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EB_LOWPT ) {
228 failiso = true;
229 }
230 if( !(ele->IsEB()) && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EE_HIGHPT ) {
231 failiso = true;
232 }
233 if( !(ele->IsEB()) && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EE_LOWPT ) {
234 failiso = true;
235 }
236
237 SelectionStatus status;
238 if( !failiso ) {
239 status.orStatus(SelectionStatus::LOOSEISO);
240 status.orStatus(SelectionStatus::TIGHTISO);
241 }
242 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
243 return status;
244
245 }
246
247
248 bool noIso(ControlFlags &, vector<SimpleLepton> &, float rho) {
249
250 return true;
251 }
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
342 if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
343
344 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
345
346 Double_t deta = (mu->Eta() - pf->Eta());
347 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
348 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
349 if (dr > 1.0) continue;
350
351 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
352
353 //
354 // Lepton Footprint Removal
355 //
356 Bool_t IsLeptonFootprint = kFALSE;
357 if (dr < 1.0) {
358
359 //
360 // Check for electrons
361 //
362 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
363 const mithep::Electron *tmpele = electronsToVeto[q];
364 // 4l electron
365 if( pf->HasTrackerTrk() ) {
366 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
367 IsLeptonFootprint = kTRUE;
368 }
369 if( pf->HasGsfTrk() ) {
370 if( pf->GsfTrk() == tmpele->GsfTrk() )
371 IsLeptonFootprint = kTRUE;
372 }
373 // PF charged
374 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) >= 1.479
375 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
376 IsLeptonFootprint = kTRUE;
377 // PF gamma
378 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) >= 1.479
379 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
380 IsLeptonFootprint = kTRUE;
381 } // loop over electrons
382
383 /* KH - commented for sync
384 //
385 // Check for muons
386 //
387 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
388 const mithep::Muon *tmpmu = muonsToVeto[q];
389 // 4l muon
390 if( pf->HasTrackerTrk() ) {
391 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
392 IsLeptonFootprint = kTRUE;
393 }
394 // PF charged
395 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
396 IsLeptonFootprint = kTRUE;
397 } // loop over muons
398 */
399
400 if (IsLeptonFootprint)
401 continue;
402
403 //
404 // Charged Iso Rings
405 //
406 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
407
408 if( dr < 0.01 ) continue; // only for muon iso mva?
409 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
410
411 // if( pf->HasTrackerTrk() ) {
412 // if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
413 // if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
414 // << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
415 // << dr << endl;
416 // }
417 // if( pf->HasGsfTrk() ) {
418 // if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
419 // if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
420 // << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
421 // << dr << endl;
422 // }
423
424 // Footprint Veto
425 if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
426 if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
427 if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
428 if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
429 if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
430 if (dr >= 0.5 && dr < 0.7) tmpChargedIso_DR0p5To0p7 += pf->Pt();
431 }
432
433 //
434 // Gamma Iso Rings
435 //
436 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
437 if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
438 if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
439 if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
440 if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
441 if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
442 if (dr >= 0.5 && dr < 0.7) tmpGammaIso_DR0p5To0p7 += pf->Pt();
443 }
444
445 //
446 // Other Neutral Iso Rings
447 //
448 else {
449 if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
450 if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
451 if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
452 if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
453 if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
454 if (dr >= 0.5 && dr < 0.7) tmpNeutralHadronIso_DR0p5To0p7 += pf->Pt();
455 }
456
457 }
458
459 }
460
461 fChargedIso_DR0p0To0p1 = fmin((tmpChargedIso_DR0p0To0p1)/mu->Pt(), 2.5);
462 fChargedIso_DR0p1To0p2 = fmin((tmpChargedIso_DR0p1To0p2)/mu->Pt(), 2.5);
463 fChargedIso_DR0p2To0p3 = fmin((tmpChargedIso_DR0p2To0p3)/mu->Pt(), 2.5);
464 fChargedIso_DR0p3To0p4 = fmin((tmpChargedIso_DR0p3To0p4)/mu->Pt(), 2.5);
465 fChargedIso_DR0p4To0p5 = fmin((tmpChargedIso_DR0p4To0p5)/mu->Pt(), 2.5);
466
467
468 double rho = 0;
469 if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
470 rho = fPUEnergyDensity->At(0)->Rho();
471 // if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
472 // rho = fPUEnergyDensity->At(0)->RhoLowEta();
473
474 // WARNING!!!!
475 // hardcode for sync ...
476 EffectiveAreaVersion = muT.kMuEAData2011;
477 // WARNING!!!!
478
479
480 fGammaIso_DR0p0To0p1 = fmax(fmin((tmpGammaIso_DR0p0To0p1
481 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p0To0p1,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
482 ,2.5)
483 ,0.0);
484 fGammaIso_DR0p1To0p2 = fmax(fmin((tmpGammaIso_DR0p1To0p2
485 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p1To0p2,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
486 ,2.5)
487 ,0.0);
488 fGammaIso_DR0p2To0p3 = fmax(fmin((tmpGammaIso_DR0p2To0p3
489 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p2To0p3,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
490 ,2.5)
491 ,0.0);
492 fGammaIso_DR0p3To0p4 = fmax(fmin((tmpGammaIso_DR0p3To0p4
493 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p3To0p4,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
494 ,2.5)
495 ,0.0);
496 fGammaIso_DR0p4To0p5 = fmax(fmin((tmpGammaIso_DR0p4To0p5
497 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p4To0p5,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
498 ,2.5)
499 ,0.0);
500
501
502
503 fNeutralHadronIso_DR0p0To0p1 = fmax(fmin((tmpNeutralHadronIso_DR0p0To0p1
504 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p0To0p1,
505 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
506 , 2.5)
507 , 0.0);
508 fNeutralHadronIso_DR0p1To0p2 = fmax(fmin((tmpNeutralHadronIso_DR0p1To0p2
509 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p1To0p2,
510 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
511 , 2.5)
512 , 0.0);
513 fNeutralHadronIso_DR0p2To0p3 = fmax(fmin((tmpNeutralHadronIso_DR0p2To0p3
514 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p2To0p3,
515 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
516 , 2.5)
517 , 0.0);
518 fNeutralHadronIso_DR0p3To0p4 = fmax(fmin((tmpNeutralHadronIso_DR0p3To0p4
519 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p3To0p4,
520 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
521 , 2.5)
522 , 0.0);
523 fNeutralHadronIso_DR0p4To0p5 = fmax(fmin((tmpNeutralHadronIso_DR0p4To0p5
524 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p4To0p5,
525 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
526 , 2.5)
527 , 0.0);
528
529
530 double mvaval = muIsoMVA->MVAValue_IsoRings( mu->Pt(),
531 mu->Eta(),
532 mu->IsGlobalMuon(),
533 mu->IsTrackerMuon(),
534 fChargedIso_DR0p0To0p1,
535 fChargedIso_DR0p1To0p2,
536 fChargedIso_DR0p2To0p3,
537 fChargedIso_DR0p3To0p4,
538 fChargedIso_DR0p4To0p5,
539 fGammaIso_DR0p0To0p1,
540 fGammaIso_DR0p1To0p2,
541 fGammaIso_DR0p2To0p3,
542 fGammaIso_DR0p3To0p4,
543 fGammaIso_DR0p4To0p5,
544 fNeutralHadronIso_DR0p0To0p1,
545 fNeutralHadronIso_DR0p1To0p2,
546 fNeutralHadronIso_DR0p2To0p3,
547 fNeutralHadronIso_DR0p3To0p4,
548 fNeutralHadronIso_DR0p4To0p5,
549 ctrl.debug);
550
551 SelectionStatus status;
552 bool pass;
553
554 pass = false;
555 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
556 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN0) pass = true;
557 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
558 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN1) pass = true;
559 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
560 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN2) pass = true;
561 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
562 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN3) pass = true;
563 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN4) pass = true;
564 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN5) pass = true;
565 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
566
567 /*
568 pass = false;
569 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
570 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN0) pass = true;
571 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
572 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN1) pass = true;
573 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
574 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN2) pass = true;
575 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
576 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN3) pass = true;
577 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN4) pass = true;
578 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN5) pass = true;
579 if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
580 */
581
582 // pass &= (fChargedIso_DR0p0To0p1 + fChargedIso_DR0p1To0p2 + fChargedIso_DR0p2To0p3 < 0.7);
583
584 status.isoMVA = mvaval;
585
586 if(ctrl.debug) {
587 cout << "returning status : " << hex << status.getStatus() << dec << endl;
588 cout << "MVAVAL : " << status.isoMVA << endl;
589 }
590 return status;
591
592 }
593
594
595 //--------------------------------------------------------------------------------------------------
596 SelectionStatus muonIsoMVASelection(ControlFlags &ctrl,
597 const mithep::Muon * mu,
598 const mithep::Vertex * vtx,
599 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
600 float rho,
601 //const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
602 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
603 vector<const mithep::Muon*> muonsToVeto,
604 vector<const mithep::Electron*> electronsToVeto)
605 //--------------------------------------------------------------------------------------------------
606 // hacked version
607 {
608
609 if( ctrl.debug ) {
610 cout << "muonIsoMVASelection :: muons to veto " << endl;
611 for( int i=0; i<muonsToVeto.size(); i++ ) {
612 const mithep::Muon * vmu = muonsToVeto[i];
613 cout << "\tpt: " << vmu->Pt()
614 << "\teta: " << vmu->Eta()
615 << "\tphi: " << vmu->Phi()
616 << endl;
617 }
618 cout << "muonIsoMVASelection :: electrson to veto " << endl;
619 for( int i=0; i<electronsToVeto.size(); i++ ) {
620 const mithep::Electron * vel = electronsToVeto[i];
621 cout << "\tpt: " << vel->Pt()
622 << "\teta: " << vel->Eta()
623 << "\tphi: " << vel->Phi()
624 << endl;
625 }
626 }
627 bool failiso=false;
628
629 //
630 // tmp iso rings
631 //
632 Double_t tmpChargedIso_DR0p0To0p1 = 0;
633 Double_t tmpChargedIso_DR0p1To0p2 = 0;
634 Double_t tmpChargedIso_DR0p2To0p3 = 0;
635 Double_t tmpChargedIso_DR0p3To0p4 = 0;
636 Double_t tmpChargedIso_DR0p4To0p5 = 0;
637 Double_t tmpChargedIso_DR0p5To0p7 = 0;
638
639 Double_t tmpGammaIso_DR0p0To0p1 = 0;
640 Double_t tmpGammaIso_DR0p1To0p2 = 0;
641 Double_t tmpGammaIso_DR0p2To0p3 = 0;
642 Double_t tmpGammaIso_DR0p3To0p4 = 0;
643 Double_t tmpGammaIso_DR0p4To0p5 = 0;
644 Double_t tmpGammaIso_DR0p5To0p7 = 0;
645
646 Double_t tmpNeutralHadronIso_DR0p0To0p1 = 0;
647 Double_t tmpNeutralHadronIso_DR0p1To0p2 = 0;
648 Double_t tmpNeutralHadronIso_DR0p2To0p3 = 0;
649 Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
650 Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
651 Double_t tmpNeutralHadronIso_DR0p5To0p7 = 0;
652
653
654
655 //
656 // final rings for the MVA
657 //
658 Double_t fChargedIso_DR0p0To0p1;
659 Double_t fChargedIso_DR0p1To0p2;
660 Double_t fChargedIso_DR0p2To0p3;
661 Double_t fChargedIso_DR0p3To0p4;
662 Double_t fChargedIso_DR0p4To0p5;
663 Double_t fChargedIso_DR0p5To0p7;
664
665 Double_t fGammaIso_DR0p0To0p1;
666 Double_t fGammaIso_DR0p1To0p2;
667 Double_t fGammaIso_DR0p2To0p3;
668 Double_t fGammaIso_DR0p3To0p4;
669 Double_t fGammaIso_DR0p4To0p5;
670 Double_t fGammaIso_DR0p5To0p7;
671
672 Double_t fNeutralHadronIso_DR0p0To0p1;
673 Double_t fNeutralHadronIso_DR0p1To0p2;
674 Double_t fNeutralHadronIso_DR0p2To0p3;
675 Double_t fNeutralHadronIso_DR0p3To0p4;
676 Double_t fNeutralHadronIso_DR0p4To0p5;
677 Double_t fNeutralHadronIso_DR0p5To0p7;
678
679
680 //
681 //Loop over PF Candidates
682 //
683 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
684
685 if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
686
687 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
688
689 Double_t deta = (mu->Eta() - pf->Eta());
690 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
691 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
692 if (dr > 1.0) continue;
693
694 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
695
696 //
697 // Lepton Footprint Removal
698 //
699 Bool_t IsLeptonFootprint = kFALSE;
700 if (dr < 1.0) {
701
702 //
703 // Check for electrons
704 //
705 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
706 const mithep::Electron *tmpele = electronsToVeto[q];
707 // 4l electron
708 if( pf->HasTrackerTrk() ) {
709 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
710 IsLeptonFootprint = kTRUE;
711 }
712 if( pf->HasGsfTrk() ) {
713 if( pf->GsfTrk() == tmpele->GsfTrk() )
714 IsLeptonFootprint = kTRUE;
715 }
716 // PF charged
717 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) >= 1.479
718 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
719 IsLeptonFootprint = kTRUE;
720 // PF gamma
721 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) >= 1.479
722 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
723 IsLeptonFootprint = kTRUE;
724 } // loop over electrons
725
726 /* KH - commented for sync
727 //
728 // Check for muons
729 //
730 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
731 const mithep::Muon *tmpmu = muonsToVeto[q];
732 // 4l muon
733 if( pf->HasTrackerTrk() ) {
734 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
735 IsLeptonFootprint = kTRUE;
736 }
737 // PF charged
738 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
739 IsLeptonFootprint = kTRUE;
740 } // loop over muons
741 */
742
743 if (IsLeptonFootprint)
744 continue;
745
746 //
747 // Charged Iso Rings
748 //
749 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
750
751 if( dr < 0.01 ) continue; // only for muon iso mva?
752 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
753
754 // if( pf->HasTrackerTrk() ) {
755 // if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
756 // if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
757 // << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
758 // << dr << endl;
759 // }
760 // if( pf->HasGsfTrk() ) {
761 // if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
762 // if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
763 // << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
764 // << dr << endl;
765 // }
766
767 // Footprint Veto
768 if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
769 if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
770 if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
771 if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
772 if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
773 if (dr >= 0.5 && dr < 0.7) tmpChargedIso_DR0p5To0p7 += pf->Pt();
774 }
775
776 //
777 // Gamma Iso Rings
778 //
779 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
780 if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
781 if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
782 if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
783 if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
784 if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
785 if (dr >= 0.5 && dr < 0.7) tmpGammaIso_DR0p5To0p7 += pf->Pt();
786 }
787
788 //
789 // Other Neutral Iso Rings
790 //
791 else {
792 if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
793 if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
794 if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
795 if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
796 if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
797 if (dr >= 0.5 && dr < 0.7) tmpNeutralHadronIso_DR0p5To0p7 += pf->Pt();
798 }
799
800 }
801
802 }
803
804 fChargedIso_DR0p0To0p1 = fmin((tmpChargedIso_DR0p0To0p1)/mu->Pt(), 2.5);
805 fChargedIso_DR0p1To0p2 = fmin((tmpChargedIso_DR0p1To0p2)/mu->Pt(), 2.5);
806 fChargedIso_DR0p2To0p3 = fmin((tmpChargedIso_DR0p2To0p3)/mu->Pt(), 2.5);
807 fChargedIso_DR0p3To0p4 = fmin((tmpChargedIso_DR0p3To0p4)/mu->Pt(), 2.5);
808 fChargedIso_DR0p4To0p5 = fmin((tmpChargedIso_DR0p4To0p5)/mu->Pt(), 2.5);
809
810
811 // double rho = 0;
812 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
813 // rho = fPUEnergyDensity->At(0)->Rho();
814 // if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
815 // rho = fPUEnergyDensity->At(0)->RhoLowEta();
816
817 // WARNING!!!!
818 // hardcode for sync ...
819 EffectiveAreaVersion = muT.kMuEAData2011;
820 // WARNING!!!!
821
822
823 fGammaIso_DR0p0To0p1 = fmax(fmin((tmpGammaIso_DR0p0To0p1
824 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p0To0p1,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
825 ,2.5)
826 ,0.0);
827 fGammaIso_DR0p1To0p2 = fmax(fmin((tmpGammaIso_DR0p1To0p2
828 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p1To0p2,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
829 ,2.5)
830 ,0.0);
831 fGammaIso_DR0p2To0p3 = fmax(fmin((tmpGammaIso_DR0p2To0p3
832 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p2To0p3,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
833 ,2.5)
834 ,0.0);
835 fGammaIso_DR0p3To0p4 = fmax(fmin((tmpGammaIso_DR0p3To0p4
836 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p3To0p4,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
837 ,2.5)
838 ,0.0);
839 fGammaIso_DR0p4To0p5 = fmax(fmin((tmpGammaIso_DR0p4To0p5
840 -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p4To0p5,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
841 ,2.5)
842 ,0.0);
843
844
845
846 fNeutralHadronIso_DR0p0To0p1 = fmax(fmin((tmpNeutralHadronIso_DR0p0To0p1
847 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p0To0p1,
848 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
849 , 2.5)
850 , 0.0);
851 fNeutralHadronIso_DR0p1To0p2 = fmax(fmin((tmpNeutralHadronIso_DR0p1To0p2
852 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p1To0p2,
853 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
854 , 2.5)
855 , 0.0);
856 fNeutralHadronIso_DR0p2To0p3 = fmax(fmin((tmpNeutralHadronIso_DR0p2To0p3
857 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p2To0p3,
858 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
859 , 2.5)
860 , 0.0);
861 fNeutralHadronIso_DR0p3To0p4 = fmax(fmin((tmpNeutralHadronIso_DR0p3To0p4
862 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p3To0p4,
863 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
864 , 2.5)
865 , 0.0);
866 fNeutralHadronIso_DR0p4To0p5 = fmax(fmin((tmpNeutralHadronIso_DR0p4To0p5
867 -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p4To0p5,
868 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
869 , 2.5)
870 , 0.0);
871
872
873 double mvaval = muIsoMVA->MVAValue_IsoRings( mu->Pt(),
874 mu->Eta(),
875 mu->IsGlobalMuon(),
876 mu->IsTrackerMuon(),
877 fChargedIso_DR0p0To0p1,
878 fChargedIso_DR0p1To0p2,
879 fChargedIso_DR0p2To0p3,
880 fChargedIso_DR0p3To0p4,
881 fChargedIso_DR0p4To0p5,
882 fGammaIso_DR0p0To0p1,
883 fGammaIso_DR0p1To0p2,
884 fGammaIso_DR0p2To0p3,
885 fGammaIso_DR0p3To0p4,
886 fGammaIso_DR0p4To0p5,
887 fNeutralHadronIso_DR0p0To0p1,
888 fNeutralHadronIso_DR0p1To0p2,
889 fNeutralHadronIso_DR0p2To0p3,
890 fNeutralHadronIso_DR0p3To0p4,
891 fNeutralHadronIso_DR0p4To0p5,
892 ctrl.debug);
893
894 SelectionStatus status;
895 bool pass;
896
897 pass = false;
898 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
899 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN0) pass = true;
900 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
901 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN1) pass = true;
902 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
903 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN2) pass = true;
904 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
905 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN3) pass = true;
906 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN4) pass = true;
907 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN5) pass = true;
908 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
909
910 /*
911 pass = false;
912 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
913 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN0) pass = true;
914 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
915 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN1) pass = true;
916 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
917 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN2) pass = true;
918 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
919 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN3) pass = true;
920 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN4) pass = true;
921 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN5) pass = true;
922 if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
923 */
924
925 // pass &= (fChargedIso_DR0p0To0p1 + fChargedIso_DR0p1To0p2 + fChargedIso_DR0p2To0p3 < 0.7);
926
927 status.isoMVA = mvaval;
928
929 if(ctrl.debug) {
930 cout << "returning status : " << hex << status.getStatus() << dec << endl;
931 cout << "MVAVAL : " << status.isoMVA << endl;
932 }
933 return status;
934
935 }
936
937
938 //--------------------------------------------------------------------------------------------------
939 void initMuonIsoMVA() {
940 //--------------------------------------------------------------------------------------------------
941 muIsoMVA = new mithep::MuonIDMVA();
942 vector<string> weightFiles;
943 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml");
944 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml");
945 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml");
946 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml");
947 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml");
948 weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml");
949 muIsoMVA->Initialize( "MuonIsoMVA",
950 mithep::MuonIDMVA::kIsoRingsV0,
951 kTRUE, weightFiles);
952 }
953
954
955
956
957 //--------------------------------------------------------------------------------------------------
958 double muonPFIso04(ControlFlags &ctrl,
959 const mithep::Muon * mu,
960 const mithep::Vertex * vtx,
961 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
962 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
963 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
964 vector<const mithep::Muon*> muonsToVeto,
965 vector<const mithep::Electron*> electronsToVeto)
966 //--------------------------------------------------------------------------------------------------
967 {
968
969 extern double gChargedIso;
970 extern double gGammaIso;
971 extern double gNeutralIso;
972
973 if( ctrl.debug ) {
974 cout << "muonIsoMVASelection :: muons to veto " << endl;
975 for( int i=0; i<muonsToVeto.size(); i++ ) {
976 const mithep::Muon * vmu = muonsToVeto[i];
977 cout << "\tpt: " << vmu->Pt()
978 << "\teta: " << vmu->Eta()
979 << "\tphi: " << vmu->Phi()
980 << endl;
981 }
982 cout << "muonIsoMVASelection :: electrson to veto " << endl;
983 for( int i=0; i<electronsToVeto.size(); i++ ) {
984 const mithep::Electron * vel = electronsToVeto[i];
985 cout << "\tpt: " << vel->Pt()
986 << "\teta: " << vel->Eta()
987 << "\tphi: " << vel->Phi()
988 << endl;
989 }
990 }
991
992 //
993 // final iso
994 //
995 Double_t fChargedIso = 0.0;
996 Double_t fGammaIso = 0.0;
997 Double_t fNeutralHadronIso = 0.0;
998
999 //
1000 //Loop over PF Candidates
1001 //
1002 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1003
1004 if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
1005 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1006
1007 Double_t deta = (mu->Eta() - pf->Eta());
1008 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
1009 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
1010 if (dr > 0.4) continue;
1011
1012 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
1013
1014 //
1015 // Lepton Footprint Removal
1016 //
1017 Bool_t IsLeptonFootprint = kFALSE;
1018 if (dr < 1.0) {
1019
1020 //
1021 // Check for electrons
1022 //
1023 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1024 const mithep::Electron *tmpele = electronsToVeto[q];
1025 // 4l electron
1026 if( pf->HasTrackerTrk() ) {
1027 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
1028 IsLeptonFootprint = kTRUE;
1029 }
1030 if( pf->HasGsfTrk() ) {
1031 if( pf->GsfTrk() == tmpele->GsfTrk() )
1032 IsLeptonFootprint = kTRUE;
1033 }
1034 // PF charged
1035 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
1036 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
1037 if( ctrl.debug) cout << "\tcharged trk, dR ("
1038 << mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta())
1039 << " matches 4L ele ..." << endl;
1040 IsLeptonFootprint = kTRUE;
1041 }
1042 // PF gamma
1043 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
1044 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
1045 IsLeptonFootprint = kTRUE;
1046 } // loop over electrons
1047
1048 /* KH - comment for sync
1049 //
1050 // Check for muons
1051 //
1052 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1053 const mithep::Muon *tmpmu = muonsToVeto[q];
1054 // 4l muon
1055 if( pf->HasTrackerTrk() ) {
1056 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
1057 IsLeptonFootprint = kTRUE;
1058 }
1059 // PF charged
1060 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
1061 IsLeptonFootprint = kTRUE;
1062 } // loop over muons
1063 */
1064
1065 if (IsLeptonFootprint)
1066 continue;
1067
1068 //
1069 // Charged Iso
1070 //
1071 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
1072
1073 //if( dr < 0.01 ) continue; // only for muon iso mva?
1074 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1075
1076
1077 // if( pf->HasTrackerTrk() ) {
1078 // if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1079 // if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
1080 // << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
1081 // << dr << endl;
1082 // }
1083 // if( pf->HasGsfTrk() ) {
1084 // if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1085 // if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
1086 // << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
1087 // << dr << endl;
1088 // }
1089
1090
1091 fChargedIso += pf->Pt();
1092 }
1093
1094 //
1095 // Gamma Iso
1096 //
1097 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1098 // KH, add to sync
1099 if( pf->Pt() > 0.5 )
1100 fGammaIso += pf->Pt();
1101 }
1102
1103 //
1104 // Other Neutrals
1105 //
1106 else {
1107 // KH, add to sync
1108 if( pf->Pt() > 0.5 )
1109 fNeutralHadronIso += pf->Pt();
1110 }
1111
1112 }
1113
1114 }
1115
1116 double rho=0;
1117 if (!(isnan(fPUEnergyDensity->At(0)->RhoKt6PFJetsCentralNeutral()) || isinf(fPUEnergyDensity->At(0)->RhoKt6PFJetsCentralNeutral())))
1118 rho = fPUEnergyDensity->At(0)->RhoKt6PFJetsCentralNeutral();
1119
1120 // WARNING!!!!
1121 // hardcode for sync ...
1122 EffectiveAreaVersion = muT.kMuEAData2012;
1123 // WARNING!!!!
1124
1125
1126
1127 double pfIso = fChargedIso + fmax(0.0,(fGammaIso + fNeutralHadronIso
1128 -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
1129 mu->Eta(),EffectiveAreaVersion)));
1130 gChargedIso = fChargedIso;
1131 gGammaIso = fGammaIso;
1132 gNeutralIso = fNeutralHadronIso;
1133
1134 return pfIso;
1135 }
1136
1137
1138
1139
1140 //--------------------------------------------------------------------------------------------------
1141 // hacked version
1142 double muonPFIso04(ControlFlags &ctrl,
1143 const mithep::Muon * mu,
1144 const mithep::Vertex * vtx,
1145 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1146 float rho,
1147 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
1148 vector<const mithep::Muon*> muonsToVeto,
1149 vector<const mithep::Electron*> electronsToVeto)
1150 //--------------------------------------------------------------------------------------------------
1151 {
1152
1153 extern double gChargedIso;
1154 extern double gGammaIso;
1155 extern double gNeutralIso;
1156
1157 if( ctrl.debug ) {
1158 cout << "muonIsoMVASelection :: muons to veto " << endl;
1159 for( int i=0; i<muonsToVeto.size(); i++ ) {
1160 const mithep::Muon * vmu = muonsToVeto[i];
1161 cout << "\tpt: " << vmu->Pt()
1162 << "\teta: " << vmu->Eta()
1163 << "\tphi: " << vmu->Phi()
1164 << endl;
1165 }
1166 cout << "muonIsoMVASelection :: electrson to veto " << endl;
1167 for( int i=0; i<electronsToVeto.size(); i++ ) {
1168 const mithep::Electron * vel = electronsToVeto[i];
1169 cout << "\tpt: " << vel->Pt()
1170 << "\teta: " << vel->Eta()
1171 << "\tphi: " << vel->Phi()
1172 << endl;
1173 }
1174 }
1175
1176 //
1177 // final iso
1178 //
1179 Double_t fChargedIso = 0.0;
1180 Double_t fGammaIso = 0.0;
1181 Double_t fNeutralHadronIso = 0.0;
1182
1183 //
1184 //Loop over PF Candidates
1185 //
1186 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1187
1188 if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
1189 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1190
1191 Double_t deta = (mu->Eta() - pf->Eta());
1192 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
1193 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
1194 if (dr > 0.4) continue;
1195
1196 if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;
1197
1198 //
1199 // Lepton Footprint Removal
1200 //
1201 Bool_t IsLeptonFootprint = kFALSE;
1202 if (dr < 1.0) {
1203
1204 //
1205 // Check for electrons
1206 //
1207 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1208 const mithep::Electron *tmpele = electronsToVeto[q];
1209 // 4l electron
1210 if( pf->HasTrackerTrk() ) {
1211 if( pf->TrackerTrk() == tmpele->TrackerTrk() )
1212 IsLeptonFootprint = kTRUE;
1213 }
1214 if( pf->HasGsfTrk() ) {
1215 if( pf->GsfTrk() == tmpele->GsfTrk() )
1216 IsLeptonFootprint = kTRUE;
1217 }
1218 // PF charged
1219 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
1220 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015)
1221 IsLeptonFootprint = kTRUE;
1222 // PF gamma
1223 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
1224 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08)
1225 IsLeptonFootprint = kTRUE;
1226 } // loop over electrons
1227
1228 /* KH - comment for sync
1229 //
1230 // Check for muons
1231 //
1232 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1233 const mithep::Muon *tmpmu = muonsToVeto[q];
1234 // 4l muon
1235 if( pf->HasTrackerTrk() ) {
1236 if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
1237 IsLeptonFootprint = kTRUE;
1238 }
1239 // PF charged
1240 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01)
1241 IsLeptonFootprint = kTRUE;
1242 } // loop over muons
1243 */
1244
1245 if (IsLeptonFootprint)
1246 continue;
1247
1248 //
1249 // Charged Iso
1250 //
1251 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
1252
1253 //if( dr < 0.01 ) continue; // only for muon iso mva?
1254 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1255
1256
1257 // if( pf->HasTrackerTrk() ) {
1258 // if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1259 // if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
1260 // << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
1261 // << dr << endl;
1262 // }
1263 // if( pf->HasGsfTrk() ) {
1264 // if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1265 // if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " "
1266 // << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
1267 // << dr << endl;
1268 // }
1269
1270
1271 fChargedIso += pf->Pt();
1272 }
1273
1274 //
1275 // Gamma Iso
1276 //
1277 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1278 // KH, add to sync
1279 if( pf->Pt() > 0.5 )
1280 fGammaIso += pf->Pt();
1281 }
1282
1283 //
1284 // Other Neutrals
1285 //
1286 else {
1287 // KH, add to sync
1288 if( pf->Pt() > 0.5 )
1289 fNeutralHadronIso += pf->Pt();
1290 }
1291
1292 }
1293
1294 }
1295
1296 // double rho = 0;
1297 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1298 // rho = fPUEnergyDensity->At(0)->Rho();
1299
1300 // WARNING!!!!
1301 // hardcode for sync ...
1302 EffectiveAreaVersion = muT.kMuEAData2011;
1303 // WARNING!!!!
1304
1305
1306 double pfIso = fChargedIso + fmax(0.0,(fGammaIso + fNeutralHadronIso
1307 -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
1308 mu->Eta(),EffectiveAreaVersion)));
1309 gChargedIso = fChargedIso;
1310 gGammaIso = fGammaIso;
1311 gNeutralIso = fNeutralHadronIso;
1312
1313 return pfIso;
1314 }
1315
1316
1317 //--------------------------------------------------------------------------------------------------
1318 SelectionStatus muonReferenceIsoSelection(ControlFlags &ctrl,
1319 const mithep::Muon * mu,
1320 const mithep::Vertex * vtx,
1321 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1322 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1323 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
1324 vector<const mithep::Muon*> muonsToVeto,
1325 vector<const mithep::Electron*> electronsToVeto)
1326 //--------------------------------------------------------------------------------------------------
1327 {
1328
1329 SelectionStatus status;
1330
1331 double pfIso = muonPFIso04( ctrl, mu, vtx, fPFCandidates, fPUEnergyDensity,
1332 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
1333 // cout << "--------------> setting muon isoPF04 to" << pfIso << endl;
1334 status.isoPF04 = pfIso;
1335 status.chisoPF04 = gChargedIso;
1336 status.gaisoPF04 = gGammaIso;
1337 status.neisoPF04 = gNeutralIso;
1338
1339 bool pass = false;
1340 if( (pfIso/mu->Pt()) < MUON_REFERENCE_PFISO_CUT ) pass = true;
1341
1342 if( pass ) {
1343 status.orStatus(SelectionStatus::LOOSEISO);
1344 status.orStatus(SelectionStatus::TIGHTISO);
1345 }
1346 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1347 return status;
1348
1349 }
1350
1351
1352 //--------------------------------------------------------------------------------------------------
1353 // hacked version
1354 SelectionStatus muonReferenceIsoSelection(ControlFlags &ctrl,
1355 const mithep::Muon * mu,
1356 const mithep::Vertex * vtx,
1357 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1358 float rho,
1359 mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
1360 vector<const mithep::Muon*> muonsToVeto,
1361 vector<const mithep::Electron*> electronsToVeto)
1362 //--------------------------------------------------------------------------------------------------
1363 {
1364
1365 SelectionStatus status;
1366
1367 double pfIso = muonPFIso04( ctrl, mu, vtx, fPFCandidates, rho,
1368 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
1369
1370 status.isoPF04 = pfIso;
1371 status.chisoPF04 = gChargedIso;
1372 status.gaisoPF04 = gGammaIso;
1373 status.neisoPF04 = gNeutralIso;
1374
1375 bool pass = false;
1376 if( (pfIso/mu->Pt()) < MUON_REFERENCE_PFISO_CUT ) pass = true;
1377
1378 if( pass ) {
1379 status.orStatus(SelectionStatus::LOOSEISO);
1380 status.orStatus(SelectionStatus::TIGHTISO);
1381 }
1382 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1383 return status;
1384
1385 }
1386
1387
1388
1389
1390 //--------------------------------------------------------------------------------------------------
1391 SelectionStatus electronIsoMVASelection(ControlFlags &ctrl,
1392 const mithep::Electron * ele,
1393 const mithep::Vertex * vtx,
1394 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1395 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1396 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1397 vector<const mithep::Muon*> muonsToVeto,
1398 vector<const mithep::Electron*> electronsToVeto)
1399 //--------------------------------------------------------------------------------------------------
1400 {
1401
1402 if( ctrl.debug ) {
1403 cout << "electronIsoMVASelection :: muons to veto " << endl;
1404 for( int i=0; i<muonsToVeto.size(); i++ ) {
1405 const mithep::Muon * vmu = muonsToVeto[i];
1406 cout << "\tpt: " << vmu->Pt()
1407 << "\teta: " << vmu->Eta()
1408 << "\tphi: " << vmu->Phi()
1409 << endl;
1410 }
1411 cout << "electronIsoMVASelection :: electrson to veto " << endl;
1412 for( int i=0; i<electronsToVeto.size(); i++ ) {
1413 const mithep::Electron * vel = electronsToVeto[i];
1414 cout << "\tpt: " << vel->Pt()
1415 << "\teta: " << vel->Eta()
1416 << "\tphi: " << vel->Phi()
1417 << "\ttrk: " << vel->TrackerTrk()
1418 << endl;
1419 }
1420 }
1421
1422 bool failiso=false;
1423
1424 //
1425 // tmp iso rings
1426 //
1427 Double_t tmpChargedIso_DR0p0To0p1 = 0;
1428 Double_t tmpChargedIso_DR0p1To0p2 = 0;
1429 Double_t tmpChargedIso_DR0p2To0p3 = 0;
1430 Double_t tmpChargedIso_DR0p3To0p4 = 0;
1431 Double_t tmpChargedIso_DR0p4To0p5 = 0;
1432
1433 Double_t tmpGammaIso_DR0p0To0p1 = 0;
1434 Double_t tmpGammaIso_DR0p1To0p2 = 0;
1435 Double_t tmpGammaIso_DR0p2To0p3 = 0;
1436 Double_t tmpGammaIso_DR0p3To0p4 = 0;
1437 Double_t tmpGammaIso_DR0p4To0p5 = 0;
1438
1439
1440 Double_t tmpNeutralHadronIso_DR0p0To0p1 = 0;
1441 Double_t tmpNeutralHadronIso_DR0p1To0p2 = 0;
1442 Double_t tmpNeutralHadronIso_DR0p2To0p3 = 0;
1443 Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
1444 Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
1445
1446
1447
1448 //
1449 // final rings for the MVA
1450 //
1451 Double_t fChargedIso_DR0p0To0p1;
1452 Double_t fChargedIso_DR0p1To0p2;
1453 Double_t fChargedIso_DR0p2To0p3;
1454 Double_t fChargedIso_DR0p3To0p4;
1455 Double_t fChargedIso_DR0p4To0p5;
1456
1457 Double_t fGammaIso_DR0p0To0p1;
1458 Double_t fGammaIso_DR0p1To0p2;
1459 Double_t fGammaIso_DR0p2To0p3;
1460 Double_t fGammaIso_DR0p3To0p4;
1461 Double_t fGammaIso_DR0p4To0p5;
1462
1463 Double_t fNeutralHadronIso_DR0p0To0p1;
1464 Double_t fNeutralHadronIso_DR0p1To0p2;
1465 Double_t fNeutralHadronIso_DR0p2To0p3;
1466 Double_t fNeutralHadronIso_DR0p3To0p4;
1467 Double_t fNeutralHadronIso_DR0p4To0p5;
1468
1469
1470 //
1471 //Loop over PF Candidates
1472 //
1473 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1474
1475 if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
1476
1477 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1478 Double_t deta = (ele->Eta() - pf->Eta());
1479 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
1480 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
1481 if (dr > 1.0) continue;
1482
1483 if(ctrl.debug) {
1484 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
1485 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(*vtx);
1486 cout << endl;
1487 }
1488
1489
1490 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
1491 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
1492
1493
1494 //
1495 // Lepton Footprint Removal
1496 //
1497 Bool_t IsLeptonFootprint = kFALSE;
1498 if (dr < 1.0) {
1499
1500
1501 //
1502 // Check for electrons
1503 //
1504
1505 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1506 const mithep::Electron *tmpele = electronsToVeto[q];
1507 double tmpdr = mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta());
1508
1509 // 4l electron
1510 if( pf->HasTrackerTrk() ) {
1511 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
1512 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
1513 IsLeptonFootprint = kTRUE;
1514 }
1515 }
1516 if( pf->HasGsfTrk() ) {
1517 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
1518 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
1519 IsLeptonFootprint = kTRUE;
1520 }
1521 }
1522 // PF charged
1523 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) >= 1.479 && tmpdr < 0.015) {
1524 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
1525 IsLeptonFootprint = kTRUE;
1526 }
1527 // PF gamma
1528 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) >= 1.479
1529 && tmpdr < 0.08) {
1530 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
1531 IsLeptonFootprint = kTRUE;
1532 }
1533 } // loop over electrons
1534
1535
1536 /* KH - comment for sync
1537 //
1538 // Check for muons
1539 //
1540 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1541 const mithep::Muon *tmpmu = muonsToVeto[q];
1542 // 4l muon
1543 if( pf->HasTrackerTrk() ) {
1544 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
1545 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
1546 IsLeptonFootprint = kTRUE;
1547 }
1548 }
1549 // PF charged
1550 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
1551 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
1552 IsLeptonFootprint = kTRUE;
1553 }
1554 } // loop over muons
1555 */
1556
1557 if (IsLeptonFootprint)
1558 continue;
1559
1560 //
1561 // Charged Iso Rings
1562 //
1563 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
1564
1565 // if( pf->HasGsfTrk() ) {
1566 // if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1567 // } else if( pf->HasTrackerTrk() ){
1568 // if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1569 // }
1570
1571 // Veto any PFmuon, or PFEle
1572 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1573
1574 // Footprint Veto
1575 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
1576
1577 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
1578 << "\ttype: " << pf->PFType()
1579 << "\ttrk: " << pf->TrackerTrk() << endl;
1580
1581 if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
1582 if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
1583 if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
1584 if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
1585 if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
1586
1587 }
1588
1589 //
1590 // Gamma Iso Rings
1591 //
1592 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
1593
1594 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.08) continue;
1595
1596 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
1597 << dr << endl;
1598
1599 if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
1600 if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
1601 if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
1602 if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
1603 if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
1604 }
1605
1606 //
1607 // Other Neutral Iso Rings
1608 //
1609 else {
1610 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
1611 << dr << endl;
1612 if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
1613 if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
1614 if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
1615 if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
1616 if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
1617 }
1618
1619 }
1620
1621 }
1622
1623 fChargedIso_DR0p0To0p1 = fmin((tmpChargedIso_DR0p0To0p1)/ele->Pt(), 2.5);
1624 fChargedIso_DR0p1To0p2 = fmin((tmpChargedIso_DR0p1To0p2)/ele->Pt(), 2.5);
1625 fChargedIso_DR0p2To0p3 = fmin((tmpChargedIso_DR0p2To0p3)/ele->Pt(), 2.5);
1626 fChargedIso_DR0p3To0p4 = fmin((tmpChargedIso_DR0p3To0p4)/ele->Pt(), 2.5);
1627 fChargedIso_DR0p4To0p5 = fmin((tmpChargedIso_DR0p4To0p5)/ele->Pt(), 2.5);
1628
1629 if(ctrl.debug) {
1630 cout << "fChargedIso_DR0p0To0p1 : " << fChargedIso_DR0p0To0p1 << endl;
1631 cout << "fChargedIso_DR0p1To0p2 : " << fChargedIso_DR0p1To0p2 << endl;
1632 cout << "fChargedIso_DR0p2To0p3 : " << fChargedIso_DR0p2To0p3 << endl;
1633 cout << "fChargedIso_DR0p3To0p4 : " << fChargedIso_DR0p3To0p4 << endl;
1634 cout << "fChargedIso_DR0p4To0p5 : " << fChargedIso_DR0p4To0p5 << endl;
1635 }
1636
1637
1638 double rho = 0;
1639 if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
1640 rho = fPUEnergyDensity->At(0)->Rho();
1641 // if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
1642 // rho = fPUEnergyDensity->At(0)->RhoLowEta();
1643
1644 // WARNING!!!!
1645 // hardcode for sync ...
1646 EffectiveAreaVersion = eleT.kEleEAData2011;
1647 // WARNING!!!!
1648
1649 if( ctrl.debug) {
1650 cout << "RHO: " << rho << endl;
1651 cout << "eta: " << ele->SCluster()->Eta() << endl;
1652 cout << "target: " << EffectiveAreaVersion << endl;
1653 cout << "effA 0-1: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
1654 ele->SCluster()->Eta(),
1655 EffectiveAreaVersion)
1656 << endl;
1657 cout << "effA 1-2: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
1658 ele->SCluster()->Eta(),
1659 EffectiveAreaVersion)
1660 << endl;
1661 cout << "effA 2-3: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
1662 ele->SCluster()->Eta(),
1663 EffectiveAreaVersion)
1664 << endl;
1665 cout << "effA 3-4: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
1666 ele->SCluster()->Eta(),
1667 EffectiveAreaVersion)
1668 << endl;
1669 }
1670
1671 fGammaIso_DR0p0To0p1 = fmax(fmin((tmpGammaIso_DR0p0To0p1
1672 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p0To0p1,
1673 ele->SCluster()->Eta(),
1674 EffectiveAreaVersion))/ele->Pt()
1675 ,2.5)
1676 ,0.0);
1677 fGammaIso_DR0p1To0p2 = fmax(fmin((tmpGammaIso_DR0p1To0p2
1678 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p1To0p2,
1679 ele->SCluster()->Eta(),
1680 EffectiveAreaVersion))/ele->Pt()
1681 ,2.5)
1682 ,0.0);
1683 fGammaIso_DR0p2To0p3 = fmax(fmin((tmpGammaIso_DR0p2To0p3
1684 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p2To0p3,
1685 ele->SCluster()->Eta()
1686 ,EffectiveAreaVersion))/ele->Pt()
1687 ,2.5)
1688 ,0.0);
1689 fGammaIso_DR0p3To0p4 = fmax(fmin((tmpGammaIso_DR0p3To0p4
1690 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p3To0p4,
1691 ele->SCluster()->Eta(),
1692 EffectiveAreaVersion))/ele->Pt()
1693 ,2.5)
1694 ,0.0);
1695 fGammaIso_DR0p4To0p5 = fmax(fmin((tmpGammaIso_DR0p4To0p5
1696 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p4To0p5,
1697 ele->SCluster()->Eta(),
1698 EffectiveAreaVersion))/ele->Pt()
1699 ,2.5)
1700 ,0.0);
1701
1702
1703 if( ctrl.debug) {
1704 cout << "fGammaIso_DR0p0To0p1: " << fGammaIso_DR0p0To0p1 << endl;
1705 cout << "fGammaIso_DR0p1To0p2: " << fGammaIso_DR0p1To0p2 << endl;
1706 cout << "fGammaIso_DR0p2To0p3: " << fGammaIso_DR0p2To0p3 << endl;
1707 cout << "fGammaIso_DR0p3To0p4: " << fGammaIso_DR0p3To0p4 << endl;
1708 cout << "fGammaIso_DR0p4To0p5: " << fGammaIso_DR0p4To0p5 << endl;
1709 }
1710
1711 fNeutralHadronIso_DR0p0To0p1 = fmax(fmin((tmpNeutralHadronIso_DR0p0To0p1
1712 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
1713 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1714 , 2.5)
1715 , 0.0);
1716 fNeutralHadronIso_DR0p1To0p2 = fmax(fmin((tmpNeutralHadronIso_DR0p1To0p2
1717 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
1718 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1719 , 2.5)
1720 , 0.0);
1721 fNeutralHadronIso_DR0p2To0p3 = fmax(fmin((tmpNeutralHadronIso_DR0p2To0p3
1722 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
1723 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
1724 , 2.5)
1725 , 0.0);
1726 fNeutralHadronIso_DR0p3To0p4 = fmax(fmin((tmpNeutralHadronIso_DR0p3To0p4
1727 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
1728 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
1729 , 2.5)
1730 , 0.0);
1731 fNeutralHadronIso_DR0p4To0p5 = fmax(fmin((tmpNeutralHadronIso_DR0p4To0p5
1732 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p4To0p5,
1733 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
1734 , 2.5)
1735 , 0.0);
1736
1737 if( ctrl.debug) {
1738 cout << "fNeutralHadronIso_DR0p0To0p1: " << fNeutralHadronIso_DR0p0To0p1 << endl;
1739 cout << "fNeutralHadronIso_DR0p1To0p2: " << fNeutralHadronIso_DR0p1To0p2 << endl;
1740 cout << "fNeutralHadronIso_DR0p2To0p3: " << fNeutralHadronIso_DR0p2To0p3 << endl;
1741 cout << "fNeutralHadronIso_DR0p3To0p4: " << fNeutralHadronIso_DR0p3To0p4 << endl;
1742 cout << "fNeutralHadronIso_DR0p4To0p5: " << fNeutralHadronIso_DR0p4To0p5 << endl;
1743 }
1744
1745 double mvaval = eleIsoMVA->MVAValue_IsoRings( ele->Pt(),
1746 ele->SCluster()->Eta(),
1747 fChargedIso_DR0p0To0p1,
1748 fChargedIso_DR0p1To0p2,
1749 fChargedIso_DR0p2To0p3,
1750 fChargedIso_DR0p3To0p4,
1751 fChargedIso_DR0p4To0p5,
1752 fGammaIso_DR0p0To0p1,
1753 fGammaIso_DR0p1To0p2,
1754 fGammaIso_DR0p2To0p3,
1755 fGammaIso_DR0p3To0p4,
1756 fGammaIso_DR0p4To0p5,
1757 fNeutralHadronIso_DR0p0To0p1,
1758 fNeutralHadronIso_DR0p1To0p2,
1759 fNeutralHadronIso_DR0p2To0p3,
1760 fNeutralHadronIso_DR0p3To0p4,
1761 fNeutralHadronIso_DR0p4To0p5,
1762 ctrl.debug);
1763
1764 SelectionStatus status;
1765 status.isoMVA = mvaval;
1766 bool pass = false;
1767
1768 Int_t subdet = 0;
1769 if (fabs(ele->SCluster()->Eta()) < 0.8) subdet = 0;
1770 else if (fabs(ele->SCluster()->Eta()) < 1.479) subdet = 1;
1771 else subdet = 2;
1772
1773 Int_t ptBin = 0;
1774 if (ele->Pt() >= 10.0) ptBin = 1;
1775
1776 Int_t MVABin = -1;
1777 if (subdet == 0 && ptBin == 0) MVABin = 0;
1778 if (subdet == 1 && ptBin == 0) MVABin = 1;
1779 if (subdet == 2 && ptBin == 0) MVABin = 2;
1780 if (subdet == 0 && ptBin == 1) MVABin = 3;
1781 if (subdet == 1 && ptBin == 1) MVABin = 4;
1782 if (subdet == 2 && ptBin == 1) MVABin = 5;
1783
1784 pass = false;
1785 if( MVABin == 0 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN0 ) pass = true;
1786 if( MVABin == 1 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN1 ) pass = true;
1787 if( MVABin == 2 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN2 ) pass = true;
1788 if( MVABin == 3 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN3 ) pass = true;
1789 if( MVABin == 4 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN4 ) pass = true;
1790 if( MVABin == 5 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN5 ) pass = true;
1791 // pass &= (fChargedIso_DR0p0To0p1 + fChargedIso_DR0p1To0p2 + fChargedIso_DR0p2To0p3 < 0.7);
1792 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
1793
1794 // pass = false;
1795 // if( MVABin == 0 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN0 ) pass = true;
1796 // if( MVABin == 1 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN1 ) pass = true;
1797 // if( MVABin == 2 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN2 ) pass = true;
1798 // if( MVABin == 3 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN3 ) pass = true;
1799 // if( MVABin == 4 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN4 ) pass = true;
1800 // if( MVABin == 5 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN5 ) pass = true;
1801 // if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
1802
1803 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
1804 return status;
1805
1806 }
1807
1808
1809 //--------------------------------------------------------------------------------------------------
1810 SelectionStatus electronIsoMVASelection(ControlFlags &ctrl,
1811 const mithep::Electron * ele,
1812 const mithep::Vertex * vtx,
1813 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
1814 float rho,
1815 //const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
1816 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
1817 vector<const mithep::Muon*> muonsToVeto,
1818 vector<const mithep::Electron*> electronsToVeto)
1819 //--------------------------------------------------------------------------------------------------
1820 // hacked version
1821 {
1822 if( ctrl.debug ) {
1823 cout << "================> hacked ele Iso MVA <======================" << endl;
1824 }
1825
1826 if( ctrl.debug ) {
1827 cout << "electronIsoMVASelection :: muons to veto " << endl;
1828 for( int i=0; i<muonsToVeto.size(); i++ ) {
1829 const mithep::Muon * vmu = muonsToVeto[i];
1830 cout << "\tpt: " << vmu->Pt()
1831 << "\teta: " << vmu->Eta()
1832 << "\tphi: " << vmu->Phi()
1833 << endl;
1834 }
1835 cout << "electronIsoMVASelection :: electrson to veto " << endl;
1836 for( int i=0; i<electronsToVeto.size(); i++ ) {
1837 const mithep::Electron * vel = electronsToVeto[i];
1838 cout << "\tpt: " << vel->Pt()
1839 << "\teta: " << vel->Eta()
1840 << "\tphi: " << vel->Phi()
1841 << "\ttrk: " << vel->TrackerTrk()
1842 << endl;
1843 }
1844 }
1845
1846 bool failiso=false;
1847
1848 //
1849 // tmp iso rings
1850 //
1851 Double_t tmpChargedIso_DR0p0To0p1 = 0;
1852 Double_t tmpChargedIso_DR0p1To0p2 = 0;
1853 Double_t tmpChargedIso_DR0p2To0p3 = 0;
1854 Double_t tmpChargedIso_DR0p3To0p4 = 0;
1855 Double_t tmpChargedIso_DR0p4To0p5 = 0;
1856
1857 Double_t tmpGammaIso_DR0p0To0p1 = 0;
1858 Double_t tmpGammaIso_DR0p1To0p2 = 0;
1859 Double_t tmpGammaIso_DR0p2To0p3 = 0;
1860 Double_t tmpGammaIso_DR0p3To0p4 = 0;
1861 Double_t tmpGammaIso_DR0p4To0p5 = 0;
1862
1863
1864 Double_t tmpNeutralHadronIso_DR0p0To0p1 = 0;
1865 Double_t tmpNeutralHadronIso_DR0p1To0p2 = 0;
1866 Double_t tmpNeutralHadronIso_DR0p2To0p3 = 0;
1867 Double_t tmpNeutralHadronIso_DR0p3To0p4 = 0;
1868 Double_t tmpNeutralHadronIso_DR0p4To0p5 = 0;
1869
1870
1871
1872 //
1873 // final rings for the MVA
1874 //
1875 Double_t fChargedIso_DR0p0To0p1;
1876 Double_t fChargedIso_DR0p1To0p2;
1877 Double_t fChargedIso_DR0p2To0p3;
1878 Double_t fChargedIso_DR0p3To0p4;
1879 Double_t fChargedIso_DR0p4To0p5;
1880
1881 Double_t fGammaIso_DR0p0To0p1;
1882 Double_t fGammaIso_DR0p1To0p2;
1883 Double_t fGammaIso_DR0p2To0p3;
1884 Double_t fGammaIso_DR0p3To0p4;
1885 Double_t fGammaIso_DR0p4To0p5;
1886
1887 Double_t fNeutralHadronIso_DR0p0To0p1;
1888 Double_t fNeutralHadronIso_DR0p1To0p2;
1889 Double_t fNeutralHadronIso_DR0p2To0p3;
1890 Double_t fNeutralHadronIso_DR0p3To0p4;
1891 Double_t fNeutralHadronIso_DR0p4To0p5;
1892
1893
1894 //
1895 //Loop over PF Candidates
1896 //
1897 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
1898
1899 if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
1900
1901 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
1902 Double_t deta = (ele->Eta() - pf->Eta());
1903 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
1904 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
1905 if (dr > 1.0) continue;
1906
1907 if(ctrl.debug) {
1908 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
1909 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(*vtx);
1910 cout << endl;
1911 }
1912
1913
1914 if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
1915 (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
1916
1917
1918 //
1919 // Lepton Footprint Removal
1920 //
1921 Bool_t IsLeptonFootprint = kFALSE;
1922 if (dr < 1.0) {
1923
1924
1925 //
1926 // Check for electrons
1927 //
1928
1929 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
1930 const mithep::Electron *tmpele = electronsToVeto[q];
1931 double tmpdr = mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta());
1932
1933 // 4l electron
1934 if( pf->HasTrackerTrk() ) {
1935 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
1936 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
1937 IsLeptonFootprint = kTRUE;
1938 }
1939 }
1940 if( pf->HasGsfTrk() ) {
1941 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
1942 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
1943 IsLeptonFootprint = kTRUE;
1944 }
1945 }
1946 // PF charged
1947 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) >= 1.479 && tmpdr < 0.015) {
1948 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
1949 IsLeptonFootprint = kTRUE;
1950 }
1951 // PF gamma
1952 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) >= 1.479
1953 && tmpdr < 0.08) {
1954 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
1955 IsLeptonFootprint = kTRUE;
1956 }
1957 } // loop over electrons
1958
1959
1960 /* KH - comment for sync
1961 //
1962 // Check for muons
1963 //
1964 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
1965 const mithep::Muon *tmpmu = muonsToVeto[q];
1966 // 4l muon
1967 if( pf->HasTrackerTrk() ) {
1968 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
1969 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
1970 IsLeptonFootprint = kTRUE;
1971 }
1972 }
1973 // PF charged
1974 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
1975 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
1976 IsLeptonFootprint = kTRUE;
1977 }
1978 } // loop over muons
1979 */
1980
1981 if (IsLeptonFootprint)
1982 continue;
1983
1984 //
1985 // Charged Iso Rings
1986 //
1987 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
1988
1989 // if( pf->HasGsfTrk() ) {
1990 // if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
1991 // } else if( pf->HasTrackerTrk() ){
1992 // if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
1993 // }
1994
1995 // Veto any PFmuon, or PFEle
1996 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;
1997
1998 // Footprint Veto
1999 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
2000
2001 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
2002 << "\ttype: " << pf->PFType()
2003 << "\ttrk: " << pf->TrackerTrk() << endl;
2004
2005 if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
2006 if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
2007 if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
2008 if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
2009 if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
2010
2011 }
2012
2013 //
2014 // Gamma Iso Rings
2015 //
2016 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
2017
2018 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.08) continue;
2019
2020 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
2021 << dr << endl;
2022
2023 if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
2024 if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
2025 if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
2026 if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
2027 if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
2028 }
2029
2030 //
2031 // Other Neutral Iso Rings
2032 //
2033 else {
2034 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
2035 << dr << endl;
2036 if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
2037 if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
2038 if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
2039 if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
2040 if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
2041 }
2042
2043 }
2044
2045 }
2046
2047 fChargedIso_DR0p0To0p1 = fmin((tmpChargedIso_DR0p0To0p1)/ele->Pt(), 2.5);
2048 fChargedIso_DR0p1To0p2 = fmin((tmpChargedIso_DR0p1To0p2)/ele->Pt(), 2.5);
2049 fChargedIso_DR0p2To0p3 = fmin((tmpChargedIso_DR0p2To0p3)/ele->Pt(), 2.5);
2050 fChargedIso_DR0p3To0p4 = fmin((tmpChargedIso_DR0p3To0p4)/ele->Pt(), 2.5);
2051 fChargedIso_DR0p4To0p5 = fmin((tmpChargedIso_DR0p4To0p5)/ele->Pt(), 2.5);
2052
2053 if(ctrl.debug) {
2054 cout << "fChargedIso_DR0p0To0p1 : " << fChargedIso_DR0p0To0p1 << endl;
2055 cout << "fChargedIso_DR0p1To0p2 : " << fChargedIso_DR0p1To0p2 << endl;
2056 cout << "fChargedIso_DR0p2To0p3 : " << fChargedIso_DR0p2To0p3 << endl;
2057 cout << "fChargedIso_DR0p3To0p4 : " << fChargedIso_DR0p3To0p4 << endl;
2058 cout << "fChargedIso_DR0p4To0p5 : " << fChargedIso_DR0p4To0p5 << endl;
2059 }
2060
2061
2062 // rho=0;
2063 // double rho = 0;
2064 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
2065 // rho = fPUEnergyDensity->At(0)->Rho();
2066 // if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta())))
2067 // rho = fPUEnergyDensity->At(0)->RhoLowEta();
2068
2069 // WARNING!!!!
2070 // hardcode for sync ...
2071 EffectiveAreaVersion = eleT.kEleEAData2011;
2072 // WARNING!!!!
2073
2074 if( ctrl.debug) {
2075 cout << "RHO: " << rho << endl;
2076 cout << "eta: " << ele->SCluster()->Eta() << endl;
2077 cout << "target: " << EffectiveAreaVersion << endl;
2078 cout << "effA 0-1: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
2079 ele->SCluster()->Eta(),
2080 EffectiveAreaVersion)
2081 << endl;
2082 cout << "effA 1-2: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
2083 ele->SCluster()->Eta(),
2084 EffectiveAreaVersion)
2085 << endl;
2086 cout << "effA 2-3: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
2087 ele->SCluster()->Eta(),
2088 EffectiveAreaVersion)
2089 << endl;
2090 cout << "effA 3-4: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
2091 ele->SCluster()->Eta(),
2092 EffectiveAreaVersion)
2093 << endl;
2094 }
2095
2096 fGammaIso_DR0p0To0p1 = fmax(fmin((tmpGammaIso_DR0p0To0p1
2097 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p0To0p1,
2098 ele->SCluster()->Eta(),
2099 EffectiveAreaVersion))/ele->Pt()
2100 ,2.5)
2101 ,0.0);
2102 fGammaIso_DR0p1To0p2 = fmax(fmin((tmpGammaIso_DR0p1To0p2
2103 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p1To0p2,
2104 ele->SCluster()->Eta(),
2105 EffectiveAreaVersion))/ele->Pt()
2106 ,2.5)
2107 ,0.0);
2108 fGammaIso_DR0p2To0p3 = fmax(fmin((tmpGammaIso_DR0p2To0p3
2109 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p2To0p3,
2110 ele->SCluster()->Eta()
2111 ,EffectiveAreaVersion))/ele->Pt()
2112 ,2.5)
2113 ,0.0);
2114 fGammaIso_DR0p3To0p4 = fmax(fmin((tmpGammaIso_DR0p3To0p4
2115 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p3To0p4,
2116 ele->SCluster()->Eta(),
2117 EffectiveAreaVersion))/ele->Pt()
2118 ,2.5)
2119 ,0.0);
2120 fGammaIso_DR0p4To0p5 = fmax(fmin((tmpGammaIso_DR0p4To0p5
2121 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p4To0p5,
2122 ele->SCluster()->Eta(),
2123 EffectiveAreaVersion))/ele->Pt()
2124 ,2.5)
2125 ,0.0);
2126
2127
2128 if( ctrl.debug) {
2129 cout << "fGammaIso_DR0p0To0p1: " << fGammaIso_DR0p0To0p1 << endl;
2130 cout << "fGammaIso_DR0p1To0p2: " << fGammaIso_DR0p1To0p2 << endl;
2131 cout << "fGammaIso_DR0p2To0p3: " << fGammaIso_DR0p2To0p3 << endl;
2132 cout << "fGammaIso_DR0p3To0p4: " << fGammaIso_DR0p3To0p4 << endl;
2133 cout << "fGammaIso_DR0p4To0p5: " << fGammaIso_DR0p4To0p5 << endl;
2134 }
2135
2136 fNeutralHadronIso_DR0p0To0p1 = fmax(fmin((tmpNeutralHadronIso_DR0p0To0p1
2137 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
2138 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
2139 , 2.5)
2140 , 0.0);
2141 fNeutralHadronIso_DR0p1To0p2 = fmax(fmin((tmpNeutralHadronIso_DR0p1To0p2
2142 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
2143 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
2144 , 2.5)
2145 , 0.0);
2146 fNeutralHadronIso_DR0p2To0p3 = fmax(fmin((tmpNeutralHadronIso_DR0p2To0p3
2147 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
2148 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
2149 , 2.5)
2150 , 0.0);
2151 fNeutralHadronIso_DR0p3To0p4 = fmax(fmin((tmpNeutralHadronIso_DR0p3To0p4
2152 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
2153 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
2154 , 2.5)
2155 , 0.0);
2156 fNeutralHadronIso_DR0p4To0p5 = fmax(fmin((tmpNeutralHadronIso_DR0p4To0p5
2157 -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p4To0p5,
2158 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
2159 , 2.5)
2160 , 0.0);
2161
2162 if( ctrl.debug) {
2163 cout << "fNeutralHadronIso_DR0p0To0p1: " << fNeutralHadronIso_DR0p0To0p1 << endl;
2164 cout << "fNeutralHadronIso_DR0p1To0p2: " << fNeutralHadronIso_DR0p1To0p2 << endl;
2165 cout << "fNeutralHadronIso_DR0p2To0p3: " << fNeutralHadronIso_DR0p2To0p3 << endl;
2166 cout << "fNeutralHadronIso_DR0p3To0p4: " << fNeutralHadronIso_DR0p3To0p4 << endl;
2167 cout << "fNeutralHadronIso_DR0p4To0p5: " << fNeutralHadronIso_DR0p4To0p5 << endl;
2168 }
2169
2170 double mvaval = eleIsoMVA->MVAValue_IsoRings( ele->Pt(),
2171 ele->SCluster()->Eta(),
2172 fChargedIso_DR0p0To0p1,
2173 fChargedIso_DR0p1To0p2,
2174 fChargedIso_DR0p2To0p3,
2175 fChargedIso_DR0p3To0p4,
2176 fChargedIso_DR0p4To0p5,
2177 fGammaIso_DR0p0To0p1,
2178 fGammaIso_DR0p1To0p2,
2179 fGammaIso_DR0p2To0p3,
2180 fGammaIso_DR0p3To0p4,
2181 fGammaIso_DR0p4To0p5,
2182 fNeutralHadronIso_DR0p0To0p1,
2183 fNeutralHadronIso_DR0p1To0p2,
2184 fNeutralHadronIso_DR0p2To0p3,
2185 fNeutralHadronIso_DR0p3To0p4,
2186 fNeutralHadronIso_DR0p4To0p5,
2187 ctrl.debug);
2188
2189 SelectionStatus status;
2190 status.isoMVA = mvaval;
2191 bool pass = false;
2192
2193 Int_t subdet = 0;
2194 if (fabs(ele->SCluster()->Eta()) < 0.8) subdet = 0;
2195 else if (fabs(ele->SCluster()->Eta()) < 1.479) subdet = 1;
2196 else subdet = 2;
2197
2198 Int_t ptBin = 0;
2199 if (ele->Pt() >= 10.0) ptBin = 1;
2200
2201 Int_t MVABin = -1;
2202 if (subdet == 0 && ptBin == 0) MVABin = 0;
2203 if (subdet == 1 && ptBin == 0) MVABin = 1;
2204 if (subdet == 2 && ptBin == 0) MVABin = 2;
2205 if (subdet == 0 && ptBin == 1) MVABin = 3;
2206 if (subdet == 1 && ptBin == 1) MVABin = 4;
2207 if (subdet == 2 && ptBin == 1) MVABin = 5;
2208
2209 pass = false;
2210 if( MVABin == 0 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN0 ) pass = true;
2211 if( MVABin == 1 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN1 ) pass = true;
2212 if( MVABin == 2 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN2 ) pass = true;
2213 if( MVABin == 3 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN3 ) pass = true;
2214 if( MVABin == 4 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN4 ) pass = true;
2215 if( MVABin == 5 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN5 ) pass = true;
2216 if( pass ) status.orStatus(SelectionStatus::LOOSEISO);
2217
2218 // pass = false;
2219 // if( MVABin == 0 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN0 ) pass = true;
2220 // if( MVABin == 1 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN1 ) pass = true;
2221 // if( MVABin == 2 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN2 ) pass = true;
2222 // if( MVABin == 3 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN3 ) pass = true;
2223 // if( MVABin == 4 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN4 ) pass = true;
2224 // if( MVABin == 5 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN5 ) pass = true;
2225 // if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
2226
2227 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
2228 return status;
2229
2230 }
2231
2232
2233 //--------------------------------------------------------------------------------------------------
2234 void initElectronIsoMVA() {
2235 //--------------------------------------------------------------------------------------------------
2236 eleIsoMVA = new mithep::ElectronIDMVA();
2237 vector<string> weightFiles;
2238 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt5To10.weights.xml");
2239 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt5To10.weights.xml");
2240 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt10ToInf.weights.xml");
2241 weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt10ToInf.weights.xml");
2242 eleIsoMVA->Initialize( "ElectronIsoMVA",
2243 mithep::ElectronIDMVA::kIsoRingsV0,
2244 kTRUE, weightFiles);
2245 }
2246
2247
2248
2249
2250 //--------------------------------------------------------------------------------------------------
2251 float electronPFIso04(ControlFlags &ctrl,
2252 const mithep::Electron * ele,
2253 const mithep::Vertex * vtx,
2254 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
2255 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
2256 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
2257 vector<const mithep::Muon*> muonsToVeto,
2258 vector<const mithep::Electron*> electronsToVeto)
2259 //--------------------------------------------------------------------------------------------------
2260 {
2261 /*
2262 if( ctrl.debug ) {
2263 cout << "electronIsoMVASelection :: muons to veto " << endl;
2264 for( int i=0; i<muonsToVeto.size(); i++ ) {
2265 const mithep::Muon * vmu = muonsToVeto[i];
2266 cout << "\tpt: " << vmu->Pt()
2267 << "\teta: " << vmu->Eta()
2268 << "\tphi: " << vmu->Phi()
2269 << endl;
2270 }
2271 cout << "electronIsoMVASelection :: electrons to veto " << endl;
2272 for( int i=0; i<electronsToVeto.size(); i++ ) {
2273 const mithep::Electron * vel = electronsToVeto[i];
2274 cout << "\tpt: " << vel->Pt()
2275 << "\teta: " << vel->Eta()
2276 << "\tphi: " << vel->Phi()
2277 << "\ttrk: " << vel->TrackerTrk()
2278 << endl;
2279 }
2280 }
2281 */
2282
2283 //
2284 // final iso
2285 //
2286 Double_t fChargedIso = 0.0;
2287 Double_t fGammaIso = 0.0;
2288 Double_t fNeutralHadronIso = 0.0;
2289
2290
2291 //
2292 //Loop over PF Candidates
2293 //
2294 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
2295
2296
2297 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
2298 Double_t deta = (ele->Eta() - pf->Eta());
2299 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
2300 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
2301
2302 if (dr > 0.4) continue;
2303 if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
2304
2305 if(ctrl.debug) {
2306 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt() << "\tdR: " << dr;
2307 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(*vtx)
2308 << "\ttrk: " << pf->HasTrackerTrk()
2309 << "\tgsf: " << pf->HasGsfTrk();
2310
2311 cout << endl;
2312 }
2313
2314
2315 //
2316 // sync : I don't think theyre doing this ...
2317 //
2318 // if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
2319 // (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) {
2320 // if( ctrl.debug ) cout << "\tskipping, matches to the electron ..." << endl;
2321 // continue;
2322 // }
2323
2324
2325 //
2326 // Lepton Footprint Removal
2327 //
2328 Bool_t IsLeptonFootprint = kFALSE;
2329 if (dr < 1.0) {
2330
2331 //
2332 // Check for electrons
2333 //
2334 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
2335 const mithep::Electron *tmpele = electronsToVeto[q];
2336 /*
2337 // 4l electron
2338 if( pf->HasTrackerTrk() ) {
2339 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
2340 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
2341 IsLeptonFootprint = kTRUE;
2342 }
2343 }
2344 if( pf->HasGsfTrk() ) {
2345 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
2346 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
2347 IsLeptonFootprint = kTRUE;
2348 }
2349 }
2350 */
2351 // PF charged
2352 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
2353 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
2354 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
2355 IsLeptonFootprint = kTRUE;
2356 }
2357 // PF gamma
2358 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
2359 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
2360 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
2361 IsLeptonFootprint = kTRUE;
2362 }
2363 } // loop over electrons
2364
2365 /* KH - comment for sync
2366 //
2367 // Check for muons
2368 //
2369 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
2370 const mithep::Muon *tmpmu = muonsToVeto[q];
2371 // 4l muon
2372 if( pf->HasTrackerTrk() ) {
2373 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
2374 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
2375 IsLeptonFootprint = kTRUE;
2376 }
2377 }
2378 // PF charged
2379 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
2380 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
2381 IsLeptonFootprint = kTRUE;
2382 }
2383 } // loop over muons
2384 */
2385
2386 if (IsLeptonFootprint)
2387 continue;
2388
2389 //
2390 // Charged Iso
2391 //
2392 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
2393
2394 // if( pf->HasTrackerTrk() )
2395 // if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
2396 // if( pf->HasGsfTrk() )
2397 // if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
2398
2399 // Veto any PFmuon, or PFEle
2400 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) {
2401 cout << "\t skipping, pf is and ele or mu .." <<endl;
2402 continue;
2403 }
2404
2405 // Footprint Veto
2406 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
2407
2408 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
2409 << "\ttype: " << pf->PFType()
2410 << "\ttrk: " << pf->TrackerTrk() << endl;
2411
2412 fChargedIso += pf->Pt();
2413 }
2414
2415 //
2416 // Gamma Iso
2417 //
2418 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
2419
2420 if (fabs(ele->SCluster()->Eta()) > 1.479) {
2421 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
2422 }
2423 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
2424 << dr << endl;
2425 // KH, add to sync
2426 // if( pf->Pt() > 0.5 )
2427 fGammaIso += pf->Pt();
2428 }
2429
2430 //
2431 // Neutral Iso
2432 //
2433 else {
2434 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
2435 << dr << endl;
2436 // KH, add to sync
2437 // if( pf->Pt() > 0.5 )
2438 fNeutralHadronIso += pf->Pt();
2439 }
2440
2441 }
2442
2443 }
2444
2445
2446 double rho=0;
2447 if (!(isnan(fPUEnergyDensity->At(0)->RhoKt6PFJets()) || isinf(fPUEnergyDensity->At(0)->RhoKt6PFJets())))
2448 rho = fPUEnergyDensity->At(0)->RhoKt6PFJets();
2449
2450
2451
2452 // WARNING!!!!
2453 // hardcode for sync ...
2454 EffectiveAreaVersion = eleT.kEleEAData2012;
2455 // WARNING!!!!
2456
2457 eleT.ElectronEffectiveArea(eleT.kEleGammaAndNeutralHadronIso04,ele->Eta(),EffectiveAreaVersion) << endl;
2458
2459
2460 double pfIso = fChargedIso + fmax(0.0,(fGammaIso + fNeutralHadronIso
2461 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaAndNeutralHadronIso04,
2462 ele->Eta(),EffectiveAreaVersion)));
2463
2464
2465 gChargedIso = fChargedIso;
2466 gGammaIso = fGammaIso;
2467 gNeutralIso = fNeutralHadronIso;
2468 return pfIso;
2469 }
2470
2471
2472
2473 //--------------------------------------------------------------------------------------------------
2474 // hacked version
2475 float electronPFIso04(ControlFlags &ctrl,
2476 const mithep::Electron * ele,
2477 const mithep::Vertex * vtx,
2478 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
2479 float rho,
2480 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
2481 vector<const mithep::Muon*> muonsToVeto,
2482 vector<const mithep::Electron*> electronsToVeto)
2483 //--------------------------------------------------------------------------------------------------
2484 {
2485
2486 if( ctrl.debug ) {
2487 cout << "electronIsoMVASelection :: muons to veto " << endl;
2488 for( int i=0; i<muonsToVeto.size(); i++ ) {
2489 const mithep::Muon * vmu = muonsToVeto[i];
2490 cout << "\tpt: " << vmu->Pt()
2491 << "\teta: " << vmu->Eta()
2492 << "\tphi: " << vmu->Phi()
2493 << endl;
2494 }
2495 cout << "electronIsoMVASelection :: electrons to veto " << endl;
2496 for( int i=0; i<electronsToVeto.size(); i++ ) {
2497 const mithep::Electron * vel = electronsToVeto[i];
2498 cout << "\tpt: " << vel->Pt()
2499 << "\teta: " << vel->Eta()
2500 << "\tphi: " << vel->Phi()
2501 << "\ttrk: " << vel->TrackerTrk()
2502 << endl;
2503 }
2504 }
2505
2506
2507 //
2508 // final iso
2509 //
2510 Double_t fChargedIso = 0.0;
2511 Double_t fGammaIso = 0.0;
2512 Double_t fNeutralHadronIso = 0.0;
2513
2514
2515 //
2516 //Loop over PF Candidates
2517 //
2518 for(int k=0; k<fPFCandidates->GetEntries(); ++k) {
2519
2520
2521 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
2522 Double_t deta = (ele->Eta() - pf->Eta());
2523 Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
2524 Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
2525
2526 if (dr > 0.4) continue;
2527 if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
2528
2529 if(ctrl.debug) {
2530 cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt() << "\tdR: " << dr;
2531 if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(*vtx)
2532 << "\ttrk: " << pf->HasTrackerTrk()
2533 << "\tgsf: " << pf->HasGsfTrk();
2534
2535 cout << endl;
2536 }
2537
2538
2539 //
2540 // sync : I don't think theyre doing this ...
2541 //
2542 // if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
2543 // (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) {
2544 // if( ctrl.debug ) cout << "\tskipping, matches to the electron ..." << endl;
2545 // continue;
2546 // }
2547
2548
2549 //
2550 // Lepton Footprint Removal
2551 //
2552 Bool_t IsLeptonFootprint = kFALSE;
2553 if (dr < 1.0) {
2554
2555 //
2556 // Check for electrons
2557 //
2558 for (Int_t q=0; q < electronsToVeto.size(); ++q) {
2559 const mithep::Electron *tmpele = electronsToVeto[q];
2560 /*
2561 // 4l electron
2562 if( pf->HasTrackerTrk() ) {
2563 if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
2564 if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
2565 IsLeptonFootprint = kTRUE;
2566 }
2567 }
2568 if( pf->HasGsfTrk() ) {
2569 if( pf->GsfTrk() == tmpele->GsfTrk() ) {
2570 if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
2571 IsLeptonFootprint = kTRUE;
2572 }
2573 }
2574 */
2575 // PF charged
2576 if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479
2577 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
2578 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
2579 IsLeptonFootprint = kTRUE;
2580 }
2581 // PF gamma
2582 if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479
2583 && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
2584 if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
2585 IsLeptonFootprint = kTRUE;
2586 }
2587 } // loop over electrons
2588
2589 /* KH - comment for sync
2590 //
2591 // Check for muons
2592 //
2593 for (Int_t q=0; q < muonsToVeto.size(); ++q) {
2594 const mithep::Muon *tmpmu = muonsToVeto[q];
2595 // 4l muon
2596 if( pf->HasTrackerTrk() ) {
2597 if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
2598 if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
2599 IsLeptonFootprint = kTRUE;
2600 }
2601 }
2602 // PF charged
2603 if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
2604 if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
2605 IsLeptonFootprint = kTRUE;
2606 }
2607 } // loop over muons
2608 */
2609
2610 if (IsLeptonFootprint)
2611 continue;
2612
2613 //
2614 // Charged Iso
2615 //
2616 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
2617
2618 // if( pf->HasTrackerTrk() )
2619 // if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
2620 // if( pf->HasGsfTrk() )
2621 // if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
2622
2623 // Veto any PFmuon, or PFEle
2624 if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) {
2625 cout << "\t skipping, pf is and ele or mu .." <<endl;
2626 continue;
2627 }
2628
2629 // Footprint Veto
2630 if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;
2631
2632 if( ctrl.debug) cout << "charged:: pt: " << pf->Pt()
2633 << "\ttype: " << pf->PFType()
2634 << "\ttrk: " << pf->TrackerTrk() << endl;
2635
2636 fChargedIso += pf->Pt();
2637 }
2638
2639 //
2640 // Gamma Iso
2641 //
2642 else if (abs(pf->PFType()) == PFCandidate::eGamma) {
2643
2644 if (fabs(ele->SCluster()->Eta()) > 1.479) {
2645 if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
2646 }
2647 if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " "
2648 << dr << endl;
2649 // KH, add to sync
2650 // if( pf->Pt() > 0.5 )
2651 fGammaIso += pf->Pt();
2652 }
2653
2654 //
2655 // Neutral Iso
2656 //
2657 else {
2658 if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " "
2659 << dr << endl;
2660 // KH, add to sync
2661 // if( pf->Pt() > 0.5 )
2662 fNeutralHadronIso += pf->Pt();
2663 }
2664
2665 }
2666
2667 }
2668
2669 // double rho = 0;
2670 // if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho())))
2671 // rho = fPUEnergyDensity->At(0)->Rho();
2672
2673 // WARNING!!!!
2674 // hardcode for sync ...
2675 EffectiveAreaVersion = eleT.kEleEAData2011;
2676 // WARNING!!!!
2677
2678
2679 double pfIso = fChargedIso + fmax(0.0,(fGammaIso + fNeutralHadronIso
2680 -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaAndNeutralHadronIso04,
2681 ele->Eta(),EffectiveAreaVersion)));
2682
2683
2684 gChargedIso = fChargedIso;
2685 gGammaIso = fGammaIso;
2686 gNeutralIso = fNeutralHadronIso;
2687 return pfIso;
2688 }
2689
2690
2691 //--------------------------------------------------------------------------------------------------
2692 SelectionStatus electronReferenceIsoSelection(ControlFlags &ctrl,
2693 const mithep::Electron * ele,
2694 const mithep::Vertex * vtx,
2695 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
2696 const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity,
2697 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
2698 vector<const mithep::Muon*> muonsToVeto,
2699 vector<const mithep::Electron*> electronsToVeto)
2700 //--------------------------------------------------------------------------------------------------
2701 {
2702
2703 SelectionStatus status;
2704
2705 double pfIso = electronPFIso04( ctrl, ele, vtx, fPFCandidates, fPUEnergyDensity,
2706 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
2707 // cout << "--------------> setting electron isoPF04 to " << pfIso << endl;
2708 status.isoPF04 = pfIso;
2709 status.chisoPF04 = gChargedIso;
2710 status.gaisoPF04 = gGammaIso;
2711 status.neisoPF04 = gNeutralIso;
2712
2713 bool pass = false;
2714 if( (pfIso/ele->Pt()) < ELECTRON_REFERENCE_PFISO_CUT ) pass = true;
2715
2716 if( pass ) {
2717 status.orStatus(SelectionStatus::LOOSEISO);
2718 status.orStatus(SelectionStatus::TIGHTISO);
2719 }
2720 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
2721 return status;
2722
2723 }
2724
2725
2726 //--------------------------------------------------------------------------------------------------
2727 // hacked version
2728 SelectionStatus electronReferenceIsoSelection(ControlFlags &ctrl,
2729 const mithep::Electron * ele,
2730 const mithep::Vertex * vtx,
2731 const mithep::Array<mithep::PFCandidate> * fPFCandidates,
2732 float rho,
2733 mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
2734 vector<const mithep::Muon*> muonsToVeto,
2735 vector<const mithep::Electron*> electronsToVeto)
2736 //--------------------------------------------------------------------------------------------------
2737 {
2738
2739 SelectionStatus status;
2740
2741 double pfIso = electronPFIso04( ctrl, ele, vtx, fPFCandidates, rho,
2742 EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
2743 status.isoPF04 = pfIso;
2744 status.chisoPF04 = gChargedIso;
2745 status.gaisoPF04 = gGammaIso;
2746 status.neisoPF04 = gNeutralIso;
2747 bool pass = false;
2748 if( (pfIso/ele->Pt()) < ELECTRON_REFERENCE_PFISO_CUT ) pass = true;
2749
2750 if( pass ) {
2751 status.orStatus(SelectionStatus::LOOSEISO);
2752 status.orStatus(SelectionStatus::TIGHTISO);
2753 }
2754 if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
2755 return status;
2756
2757 }