ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/IsolationSelection.cc
Revision: 1.35
Committed: Mon Dec 17 17:12:27 2012 UTC (12 years, 5 months ago) by dkralph
Content type: text/plain
Branch: MAIN
CVS Tags: compiled
Changes since 1.34: +34 -304 lines
Log Message:
*** empty log message ***

File Contents

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