ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/MuonSelection.cc
Revision: 1.33
Committed: Mon Jul 9 08:53:53 2012 UTC (12 years, 10 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Changes since 1.32: +30 -52 lines
Log Message:
Change names on lepton tag ID functions; change rootcint flags

File Contents

# Content
1 #include <math.h>
2 #include <iostream>
3
4 #include "MuonSelection.h"
5 #include "ParseArgs.h"
6 #include "SelectionStatus.h"
7
8 #include "Track.h"
9 #include "MuonTools.h"
10 #include "MuonQuality.h"
11 #include "MuonIDMVA.h"
12
13 mithep::MuonIDMVA * muIDMVA;
14 mithep::MuonTools muTools;
15
16 extern vector<bool> PFnoPUflag;
17 extern Float_t computePFMuonIso(const mithep::Muon *mu,
18 const mithep::Vertex * vtx,
19 const mithep::Array<mithep::PFCandidate> * pfCandidates,
20 const Double_t dRMax);
21
22 bool muonPOG2012CutBasedIDTight(const mithep::Muon *mu, const mithep::Vertex * vtx, const mithep::Array<mithep::PFCandidate> * pfCandidates, const mithep::Array<mithep::PileupEnergyDensity> * puEnergyDensity)
23 {
24
25 Float_t charged_iso = 0;
26
27 Float_t sumChargedHadronPt = 0;
28 Float_t sumNeutralHadronPt = 0;
29 Float_t sumPhotonPt = 0;
30 Float_t sumPUPt = 0;
31
32 for(unsigned k=0; k<pfCandidates->GetEntries(); ++k) {
33 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfCandidates)[k]);
34 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
35 if(dr > 0.4) continue;
36 if(pf->TrackerTrk() && mu->TrackerTrk() && pf->TrackerTrk() == mu->TrackerTrk()) continue;
37
38 if(PFnoPUflag[k]){ //pf no pileup
39 if(abs(pf->PFType()) == mithep::PFCandidate::eHadron)
40 sumChargedHadronPt += pf->Pt();
41 else if (pf->Pt() > 0.5 && abs(pf->PFType()) == mithep::PFCandidate::eNeutralHadron)
42 sumNeutralHadronPt += pf->Pt();
43 else if (pf->Pt() > 0.5 && abs(pf->PFType()) == mithep::PFCandidate::eGamma)
44 sumPhotonPt += pf->Pt();
45 }
46 else {//pf pile up
47 if(pf->Charge() != 0)
48 sumPUPt += pf->Pt();
49 }
50
51 }
52 double iso_sum = sumChargedHadronPt+ max(0.,sumNeutralHadronPt+sumPhotonPt-0.5*sumPUPt);
53
54 bool pass = kTRUE;
55 if(!mu->IsGlobalMuon()) pass = kFALSE;
56
57 //there is not necessarily a global track or a tracker track if is global muon is false
58 if(!pass) return pass;
59
60 if(!mu->IsPFMuon()) pass = kFALSE;
61 if(mu->GlobalTrk()->RChi2() > 10) pass= kFALSE;
62 if(mu->NValidHits() == 0 ) pass= kFALSE;
63 if(mu->NMatches() <= 1) pass = kFALSE;
64 if(fabs(mu->TrackerTrk()->D0Corrected(*vtx)) > 0.2) pass = kFALSE;
65 if(fabs(mu->TrackerTrk()->DzCorrected(*vtx)) > 0.5) pass = kFALSE;
66 if(mu->BestTrk()->NPixelHits() == 0) pass = kFALSE;
67 if(mu->NTrkLayersHit() <= 5) pass = kFALSE;
68 if(iso_sum/mu->Pt() > 0.12 ) pass = kFALSE;
69
70 return pass;
71
72
73 }
74 //----------------------------------------------------------------------------------------
75 bool muon2012CutBasedIDTightVersionWithOldIsoThatWorksOn2011(ControlFlags &ctrl, const mithep::Muon *mu, const mithep::Vertex * vtx, const mithep::Array<mithep::PFCandidate> * pfCandidates,
76 const mithep::Array<mithep::PileupEnergyDensity> * puEnergyDensity, mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion
77 )
78 {
79
80 mithep::MuonTools muT;
81
82 Float_t charged_iso = 0;
83 Float_t neutral_iso = 0;
84 Float_t gamma_iso = 0;
85
86 for(int k=0; k<pfCandidates->GetEntries(); ++k) {
87
88 if( !(PFnoPUflag[k]) ) continue;
89 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfCandidates)[k]);
90
91 Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
92 if (dr > 0.4) continue;
93
94 if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {
95
96 if (abs(pf->PFType()) == mithep::PFCandidate::eElectron || abs(pf->PFType()) == mithep::PFCandidate::eMuon) continue;
97
98 charged_iso += pf->Pt();
99 }
100
101 else if (abs(pf->PFType()) == mithep::PFCandidate::eGamma) {
102
103 gamma_iso += pf->Pt();
104 }
105
106 else {
107 neutral_iso += pf->Pt();
108 }
109 }
110
111 double iso_sum=999;
112 double rho;
113 if(ctrl.era == 2012) {
114 EffectiveAreaVersion = mithep::MuonTools::kMuEAData2012;
115 rho = puEnergyDensity->At(0)->RhoKt6PFJetsCentralNeutral();
116
117 iso_sum = charged_iso + fmax(0.0,(gamma_iso + neutral_iso
118 -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
119 mu->Eta(),EffectiveAreaVersion)));
120 } else if(ctrl.era == 2011) {
121 setMuonRhoAndEaVersion(rho,EffectiveAreaVersion,puEnergyDensity);
122 iso_sum = charged_iso + fmax(0.0,(gamma_iso + neutral_iso
123 -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
124 mu->Eta(),EffectiveAreaVersion)));
125 } else assert(0);
126
127 bool pass = kTRUE;
128 if(!mu->IsGlobalMuon()) pass = kFALSE;
129
130 //there is not necessarily a global track or a tracker track if is global muon is false
131 if(!pass) return pass;
132
133 if(!mu->IsPFMuon() && ctrl.era==2012) pass = kFALSE;
134 if(mu->GlobalTrk()->RChi2() > 10) pass = kFALSE;
135 if(mu->NValidHits() == 0 ) pass = kFALSE;
136 if(mu->NMatches() <= 1) pass = kFALSE;
137 if(fabs(mu->TrackerTrk()->D0Corrected(*vtx)) > 0.2) pass = kFALSE;
138 if(fabs(mu->TrackerTrk()->DzCorrected(*vtx)) > 0.5) pass = kFALSE;
139 if(mu->BestTrk()->NPixelHits() == 0) pass = kFALSE;
140 if(mu->NTrkLayersHit() <= 5) pass = kFALSE;
141 if(iso_sum/mu->Pt() > 0.12 ) pass = kFALSE;
142
143 return pass;
144
145
146 }
147
148 bool trackToPFMuonSelection( const mithep::Track *track, const mithep::MuonCol * muoncol )
149 {
150 bool pass = false;
151 for(UInt_t i=0; i<muoncol->GetEntries(); ++i) {
152 const mithep::Muon *mu = muoncol->At(i);
153 if(mu->TrackerTrk() == track && (mu->IsGlobalMuon() || mu->IsTrackerMuon()) && mu->IsPFMuon())
154 pass = true;
155 }
156
157 return pass;
158
159 }
160
161 //--------------------------------------------------------------------------------------------------
162 SelectionStatus muonDummyVeto(ControlFlags &ctrl,
163 const mithep::Muon *muon,
164 const mithep::Vertex * vtx)
165 //--------------------------------------------------------------------------------------------------
166 {
167 SelectionStatus status;
168 status.setStatus(SelectionStatus::PRESELECTION);
169 return status;
170 }
171
172 //--------------------------------------------------------------------------------------------------
173 SelectionStatus muonCutBasedVeto(ControlFlags &ctrl,
174 const mithep::Muon *muon,
175 const mithep::Vertex *vtx)
176 //--------------------------------------------------------------------------------------------------
177 {
178 //
179 // Loose cut-based ID for isolation veto
180 //
181 bool ret = true;
182
183 if(!(muon->IsGlobalMuon() || muon->IsTrackerMuon())) ret = false;
184 if( muon->NValidHits() < 11 ) ret = false;
185 if( fabs(muon->Ip3dPVSignificance()) >= 4 ) ret = false;
186
187 SelectionStatus status;
188 if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
189 return status;
190 }
191
192
193 //--------------------------------------------------------------------------------------------------
194 SelectionStatus noPreSelection( ControlFlags &ctrl, const mithep::Muon * mu )
195
196 //--------------------------------------------------------------------------------------------------
197 {
198 SelectionStatus status;
199 status.setStatus(SelectionStatus::PRESELECTION);
200 if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
201 return status;
202 }
203
204 //--------------------------------------------------------------------------------------------------
205 SelectionStatus muonPreSelection( ControlFlags &ctrl,
206 const mithep::Muon * mu,
207 const mithep::Vertex * vtx,
208 const mithep::Array<mithep::PFCandidate> * pfCandidates )
209 //--------------------------------------------------------------------------------------------------
210 {
211 bool ret = true;
212
213 ret &= ( fabs(mu->Ip3dPVSignificance()) < 100 );
214 ret &= ( mu->Pt() >= 5 );
215 ret &= ( fabs(mu->Eta()) <= 2.4 );
216 ret &= ( mu->IsTrackerMuon() || mu->IsGlobalMuon() );
217 // if( mu->IsTrackerMuon() && mu->HasTrackerTrk() && !mu->IsGlobalMuon() )
218 // ret &= ( mu->Quality().QualityMask().Mask() & mithep::MuonQuality::AllArbitrated);
219 ret &= ( mu->IsoR03SumPt()/mu->Pt() < 0.7 );
220
221 SelectionStatus status;
222 if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
223 if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
224 return status;
225 }
226
227
228
229 //--------------------------------------------------------------------------------------------------
230 SelectionStatus muonReferencePreSelection( ControlFlags &ctrl,
231 const mithep::Muon * mu,
232 const mithep::Vertex * vtx,
233 const mithep::Array<mithep::PFCandidate> * pfCandidates )
234 //--------------------------------------------------------------------------------------------------
235 {
236 bool ret = true;
237
238 if(ctrl.debug) cout << "inside muonpresel ... " << endl;
239 ret &= ( fabs(mu->Ip3dPVSignificance()) < 100 );
240 if(ctrl.debug) cout << "and pass IP (" << mu->Ip3dPVSignificance() << ") ? ... " << ret << endl;
241 ret &= ( mu->Pt() >= 5 );
242 if(ctrl.debug) cout << "and >5 GeV (" << mu->Pt() << ") ? ... " << ret << endl;
243 ret &= ( fabs(mu->Eta()) <= 2.4 );
244 if(ctrl.debug) cout << "and < 2.4 eta (" << mu->Eta() << ")? ... " << ret << endl;
245 ret &= ( mu->IsTrackerMuon() || mu->IsGlobalMuon() );
246 if(ctrl.debug) cout << "is Tracker or Global? ... " << ret << endl;
247 ret &= (mu->HasTrackerTrk() &&
248 fabs(mu->TrackerTrk()->D0Corrected(*vtx)) < 0.5 &&
249 fabs(mu->TrackerTrk()->DzCorrected(*vtx)) < 1.0);
250 if(ctrl.debug) cout << "d0 and dz? ... " << ret << endl;
251
252 // if( mu->IsTrackerMuon() && mu->HasTrackerTrk() && !mu->IsGlobalMuon() )
253 // ret &= ( mu->Quality().QualityMask().Mask() & mithep::MuonQuality::AllArbitrated);
254 // if(ctrl.debug) cout << "if isTrackerMuon, arbitrated ? ... " << ret << endl;
255 // ret &= ( mu->BestTrk()->DzCorrected(vtx) < 0.1);
256 // if( ctrl.debug ) cout << "elepresel : ret after dZcorr : " << ret << endl;
257
258
259 SelectionStatus status;
260 if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
261 if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
262 return status;
263 }
264 //--------------------------------------------------------------------------------------------------
265 SelectionStatus muonPreSelectionNoD0DzIP( ControlFlags &ctrl,
266 const mithep::Muon * mu,
267 const mithep::Vertex * vtx,
268 const mithep::Array<mithep::PFCandidate> * pfCandidates )
269 //--------------------------------------------------------------------------------------------------
270 {
271 bool ret = true;
272
273 if(ctrl.debug) cout << "inside muonpresel ... " << endl;
274 ret &= ( mu->Pt() >= 5 );
275 if(ctrl.debug) cout << "and >5 GeV (" << mu->Pt() << ") ? ... " << ret << endl;
276 ret &= ( fabs(mu->Eta()) <= 2.4 );
277 if(ctrl.debug) cout << "and < 2.4 eta (" << mu->Eta() << ")? ... " << ret << endl;
278 ret &= ( mu->IsTrackerMuon() || mu->IsGlobalMuon() );
279 if(ctrl.debug) cout << "is Tracker or Global? ... " << ret << endl;
280
281 // if( mu->IsTrackerMuon() && mu->HasTrackerTrk() && !mu->IsGlobalMuon() )
282 // ret &= ( mu->Quality().QualityMask().Mask() & mithep::MuonQuality::AllArbitrated);
283 // if(ctrl.debug) cout << "if isTrackerMuon, arbitrated ? ... " << ret << endl;
284 // ret &= ( mu->BestTrk()->DzCorrected(vtx) < 0.1);
285 // if( ctrl.debug ) cout << "elepresel : ret after dZcorr : " << ret << endl;
286
287
288 SelectionStatus status;
289 if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
290 if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
291 return status;
292 }
293
294
295
296
297 //--------------------------------------------------------------------------------------------------
298 SelectionStatus muonIDMVASelection(ControlFlags &ctrl,
299 const mithep::Muon *mu,
300 const mithep::Vertex * vtx )
301 //--------------------------------------------------------------------------------------------------
302 {
303 double global_rchi2 = (mu->IsGlobalMuon()) ? mu->GlobalTrk()->RChi2() : 0.;
304
305 //
306 // WARNING !!! KH hacked out to sync
307 //
308 /*
309 double mvaval = muIDMVA->MVAValue_ID(mu->Pt(),
310 mu->Eta(),
311 mu->IsGlobalMuon(),
312 mu->IsTrackerMuon(),
313 mu->TrackerTrk()->RChi2(),
314 global_rchi2,
315 (Double_t)(mu->NValidHits()),
316 (Double_t)(mu->TrackerTrk()->NHits()),
317 (Double_t)(mu->TrackerTrk()->NPixelHits()),
318 (Double_t)(mu->NMatches()),
319 mu->TrkKink(),
320 muTools.GetSegmentCompatability(mu),
321 muTools.GetCaloCompatability(mu,true,true),
322 mu->HadEnergy(),
323 mu->EmEnergy(),
324 mu->HadS9Energy(),
325 mu->EmS9Energy(),
326 (Bool_t)(ctrl.debug) );
327 */
328 double mvaval = 0.0;
329
330 SelectionStatus status;
331 bool pass = false;
332
333 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
334 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_IDMVA_CUT_BIN0) pass = true;
335 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
336 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN1) pass = true;
337 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
338 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_IDMVA_CUT_BIN2) pass = true;
339 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
340 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN3) pass = true;
341 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_IDMVA_CUT_BIN4 ) pass = true;
342 else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_IDMVA_CUT_BIN5 ) pass = true;
343
344
345
346 if( pass ) {
347 status.orStatus(SelectionStatus::LOOSEID);
348 status.orStatus(SelectionStatus::TIGHTID);
349 }
350
351 return status;
352 }
353
354 //--------------------------------------------------------------------------------------------------
355 void initMuonIDMVA()
356 //--------------------------------------------------------------------------------------------------
357 {
358 muIDMVA = new mithep::MuonIDMVA();
359 vector<string> weightFiles;
360 weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_barrel_lowpt_V2.weights.xml");
361 weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_barrel_highpt_V2.weights.xml");
362 weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_endcap_lowpt_V2.weights.xml");
363 weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_endcap_highpt_V2.weights.xml");
364 weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_tracker_V2.weights.xml");
365 weightFiles.push_back("./data/MuonIDMVAWeights/MuonIDMVA_BDTG_global_V2.weights.xml");
366 muIDMVA->Initialize( "MuonIDMVA",
367 mithep::MuonIDMVA::kIDV0,
368 kTRUE, weightFiles);
369 }
370
371 //--------------------------------------------------------------------------------------------------
372 SelectionStatus muonIDPFSelection(ControlFlags &ctrl,
373 const mithep::Muon *mu,
374 const mithep::Vertex * vtx,
375 const mithep::Array<mithep::PFCandidate> * pfCandidates )
376 //--------------------------------------------------------------------------------------------------
377 {
378 bool pass = false;
379 SelectionStatus status; // init'ed to FAIL
380
381 // check that it matches to a PF muon
382 for( int i=0; i<pfCandidates->GetEntries(); i++ ) {
383 // tmp
384 if( !(PFnoPUflag[i]) ) continue; // my PF no PU hack
385
386 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfCandidates)[i]);
387 if( (pf->TrackerTrk() == mu->TrackerTrk()) && abs(pf->PFType()) == mithep::PFCandidate::eMuon ) {
388 pass = true;
389 break;
390 }
391 }
392 //
393 // if(mu->IsPFMuon()) pass = true;
394 // else pass = false;
395
396 if( pass ) {
397 status.orStatus(SelectionStatus::LOOSEID);
398 status.orStatus(SelectionStatus::TIGHTID);
399 }
400
401 return status;
402 }
403 //----------------------------------------------------------------------------------------
404 SelectionStatus PassWwMuonSel(const mithep::Muon *mu, const mithep::Vertex *vtx,
405 const mithep::Array<mithep::PFCandidate> *pfArr
406 )
407 {
408 bool passID=true;
409 bool passIso=true;
410
411 if (mu->Pt() < 5)
412 passID = false;
413 if (fabs(mu->Eta()) > 2.4)
414 passID = false;
415 if (mu->BestTrk()->PtErr()/mu->Pt() > 0.1)
416 passID = false;
417 if (fabs(mu->BestTrk()->DzCorrected(*vtx)) > 0.1)
418 passID = false;
419 Bool_t isGlobal = (mu->IsGlobalMuon()) && (mu->BestTrk()->RChi2() < 10) &&
420 (mu->NMatches() > 1) && (mu->NValidHits() > 0);
421 Bool_t isTracker = (mu->IsTrackerMuon()) &&
422 (mu->Quality().Quality(mithep::MuonQuality::TMLastStationTight));
423 if (!isGlobal && !isTracker)
424 passID = false;
425 int ntrkhits = (mu->HasTrackerTrk()) ? mu->TrackerTrk()->NHits() : 0;
426 if (ntrkhits < 10)
427 passID = false;
428 if (mu->BestTrk()->NPixelHits() < 1)
429 passID = false;
430 if (fabs(mu->BestTrk()->D0Corrected(*vtx)) > 0.02)
431 passID = false;
432
433 // note: this isn't really ww muon isolation
434 float pfiso = computePFMuonIso(mu,vtx,pfArr,0.3);
435
436 if (pfiso > 0.2*mu->Pt())
437 passIso = false;
438
439 SelectionStatus status;
440 if(passID) status |= SelectionStatus::TIGHTIDANDPRE;
441 if(passIso) status |= SelectionStatus::TIGHTISOANDPRE;
442
443 return status;
444
445 }