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

# User Rev Content
1 khahn 1.1 #include <math.h>
2     #include <iostream>
3    
4     #include "MuonSelection.h"
5     #include "ParseArgs.h"
6     #include "SelectionStatus.h"
7    
8 khahn 1.4 #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 khahn 1.19 extern vector<bool> PFnoPUflag;
17 khahn 1.4 extern Float_t computePFMuonIso(const mithep::Muon *mu,
18 anlevin 1.21 const mithep::Vertex * vtx,
19 khahn 1.4 const mithep::Array<mithep::PFCandidate> * pfCandidates,
20     const Double_t dRMax);
21    
22 anlevin 1.32 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 dkralph 1.33 }
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 anlevin 1.27 {
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 dkralph 1.33
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 anlevin 1.27
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 dkralph 1.33 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 anlevin 1.27
143     return pass;
144    
145    
146     }
147    
148 anlevin 1.24 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 khahn 1.4 //--------------------------------------------------------------------------------------------------
162 khahn 1.8 SelectionStatus muonDummyVeto(ControlFlags &ctrl,
163     const mithep::Muon *muon,
164 anlevin 1.21 const mithep::Vertex * vtx)
165 khahn 1.8 //--------------------------------------------------------------------------------------------------
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 anlevin 1.21 const mithep::Vertex *vtx)
176 khahn 1.8 //--------------------------------------------------------------------------------------------------
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 khahn 1.4 SelectionStatus noPreSelection( ControlFlags &ctrl, const mithep::Muon * mu )
195    
196     //--------------------------------------------------------------------------------------------------
197     {
198 anlevin 1.3 SelectionStatus status;
199     status.setStatus(SelectionStatus::PRESELECTION);
200     if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
201     return status;
202     }
203 khahn 1.1
204 khahn 1.4 //--------------------------------------------------------------------------------------------------
205     SelectionStatus muonPreSelection( ControlFlags &ctrl,
206     const mithep::Muon * mu,
207 anlevin 1.21 const mithep::Vertex * vtx,
208 khahn 1.4 const mithep::Array<mithep::PFCandidate> * pfCandidates )
209     //--------------------------------------------------------------------------------------------------
210     {
211     bool ret = true;
212 khahn 1.20
213     ret &= ( fabs(mu->Ip3dPVSignificance()) < 100 );
214 khahn 1.14 ret &= ( mu->Pt() >= 5 );
215 khahn 1.7 ret &= ( fabs(mu->Eta()) <= 2.4 );
216 khahn 1.11 ret &= ( mu->IsTrackerMuon() || mu->IsGlobalMuon() );
217 khahn 1.20 // 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 khahn 1.1
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 khahn 1.4
228 khahn 1.1
229 khahn 1.13 //--------------------------------------------------------------------------------------------------
230     SelectionStatus muonReferencePreSelection( ControlFlags &ctrl,
231     const mithep::Muon * mu,
232 anlevin 1.21 const mithep::Vertex * vtx,
233 khahn 1.13 const mithep::Array<mithep::PFCandidate> * pfCandidates )
234     //--------------------------------------------------------------------------------------------------
235     {
236     bool ret = true;
237 khahn 1.20
238 khahn 1.13 if(ctrl.debug) cout << "inside muonpresel ... " << endl;
239 khahn 1.16 ret &= ( fabs(mu->Ip3dPVSignificance()) < 100 );
240 khahn 1.13 if(ctrl.debug) cout << "and pass IP (" << mu->Ip3dPVSignificance() << ") ? ... " << ret << endl;
241 khahn 1.14 ret &= ( mu->Pt() >= 5 );
242 khahn 1.13 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 anlevin 1.21 if(ctrl.debug) cout << "is Tracker or Global? ... " << ret << endl;
247 khahn 1.26 ret &= (mu->HasTrackerTrk() &&
248     fabs(mu->TrackerTrk()->D0Corrected(*vtx)) < 0.5 &&
249     fabs(mu->TrackerTrk()->DzCorrected(*vtx)) < 1.0);
250 anlevin 1.21 if(ctrl.debug) cout << "d0 and dz? ... " << ret << endl;
251    
252 khahn 1.13 // 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 dkralph 1.29 //--------------------------------------------------------------------------------------------------
265 khahn 1.25 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 khahn 1.13
295    
296 anlevin 1.3
297 khahn 1.4 //--------------------------------------------------------------------------------------------------
298     SelectionStatus muonIDMVASelection(ControlFlags &ctrl,
299     const mithep::Muon *mu,
300 anlevin 1.21 const mithep::Vertex * vtx )
301 khahn 1.4 //--------------------------------------------------------------------------------------------------
302     {
303     double global_rchi2 = (mu->IsGlobalMuon()) ? mu->GlobalTrk()->RChi2() : 0.;
304 anlevin 1.3
305 khahn 1.13 //
306     // WARNING !!! KH hacked out to sync
307     //
308     /*
309 khahn 1.4 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 khahn 1.11 (Double_t)(mu->TrackerTrk()->NHits()),
317 khahn 1.4 (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 khahn 1.13 */
328     double mvaval = 0.0;
329 anlevin 1.3
330 khahn 1.4 SelectionStatus status;
331     bool pass = false;
332    
333     if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
334 khahn 1.7 && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_IDMVA_CUT_BIN0) pass = true;
335 khahn 1.4 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
336 khahn 1.7 && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN1) pass = true;
337 khahn 1.4 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
338 khahn 1.7 && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_IDMVA_CUT_BIN2) pass = true;
339 khahn 1.4 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
340     && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN3) pass = true;
341 khahn 1.11 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 khahn 1.4
345    
346     if( pass ) {
347     status.orStatus(SelectionStatus::LOOSEID);
348     status.orStatus(SelectionStatus::TIGHTID);
349 anlevin 1.3 }
350    
351 khahn 1.4 return status;
352     }
353 anlevin 1.3
354 khahn 1.4 //--------------------------------------------------------------------------------------------------
355     void initMuonIDMVA()
356     //--------------------------------------------------------------------------------------------------
357     {
358     muIDMVA = new mithep::MuonIDMVA();
359     vector<string> weightFiles;
360 khahn 1.6 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 khahn 1.4 muIDMVA->Initialize( "MuonIDMVA",
367     mithep::MuonIDMVA::kIDV0,
368     kTRUE, weightFiles);
369     }
370 khahn 1.9
371     //--------------------------------------------------------------------------------------------------
372     SelectionStatus muonIDPFSelection(ControlFlags &ctrl,
373     const mithep::Muon *mu,
374 anlevin 1.21 const mithep::Vertex * vtx,
375 khahn 1.9 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 khahn 1.19 // tmp
384     if( !(PFnoPUflag[i]) ) continue; // my PF no PU hack
385    
386 khahn 1.9 const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*pfCandidates)[i]);
387 khahn 1.19 if( (pf->TrackerTrk() == mu->TrackerTrk()) && abs(pf->PFType()) == mithep::PFCandidate::eMuon ) {
388 khahn 1.17 pass = true;
389     break;
390 khahn 1.9 }
391     }
392 anlevin 1.22 //
393     // if(mu->IsPFMuon()) pass = true;
394     // else pass = false;
395 khahn 1.9
396 khahn 1.19 if( pass ) {
397     status.orStatus(SelectionStatus::LOOSEID);
398     status.orStatus(SelectionStatus::TIGHTID);
399     }
400 khahn 1.9
401     return status;
402     }
403 dkralph 1.28 //----------------------------------------------------------------------------------------
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 dkralph 1.29 if (pfiso > 0.2*mu->Pt())
437 dkralph 1.28 passIso = false;
438    
439     SelectionStatus status;
440     if(passID) status |= SelectionStatus::TIGHTIDANDPRE;
441     if(passIso) status |= SelectionStatus::TIGHTISOANDPRE;
442    
443     return status;
444    
445     }