ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/LeptonSelection/src/MuonSelection.cc
Revision: 1.5
Committed: Thu Apr 26 07:02:28 2012 UTC (13 years ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.4: +6 -6 lines
Log Message:
first pass port to bambu

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 Float_t computePFMuonIso(const mithep::Muon *mu,
17 const mithep::Vertex & vtx,
18 const mithep::Array<mithep::PFCandidate> * pfCandidates,
19 const Double_t dRMax);
20
21 //--------------------------------------------------------------------------------------------------
22 SelectionStatus noPreSelection( ControlFlags &ctrl, const mithep::Muon * mu )
23
24 //--------------------------------------------------------------------------------------------------
25 {
26 SelectionStatus status;
27 status.setStatus(SelectionStatus::PRESELECTION);
28 if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
29 return status;
30 }
31
32 //--------------------------------------------------------------------------------------------------
33 SelectionStatus muonPreSelection( ControlFlags &ctrl,
34 const mithep::Muon * mu,
35 const mithep::Vertex & vtx,
36 const mithep::Array<mithep::PFCandidate> * pfCandidates )
37 //--------------------------------------------------------------------------------------------------
38 {
39 bool ret = true;
40 if(ctrl.debug) cout << "inside muonpresel ... " << endl;
41 // bool ret = isMuFO(mu,vtx,pfCandidates);
42 if(ctrl.debug) cout << "iS fo? ... " << ret << endl;
43 ret &= ( fabs(mu->Ip3dPVSignificance()) < 4 );
44 if(ctrl.debug) cout << "and pass IP (" << mu->Ip3dPVSignificance() << ") ? ... " << ret << endl;
45 ret &= ( mu->Pt() > 5 );
46 if(ctrl.debug) cout << "and >5 GeV ? ... " << ret << endl;
47 ret &= ( fabs(mu->Eta()) < 2.4 );
48 if(ctrl.debug) cout << "and < 2.4 eta ? ... " << ret << endl;
49 ret &= (mu->IsTrackerMuon() && mu->HasTrackerTrk());
50 if(ctrl.debug) cout << "and isTrackerMuon ? ... " << ret << endl;
51 ret &= (mu->IsoR03SumPt()/mu->Pt() < 0.7 );
52 if(ctrl.debug) cout << "and loose trk iso ? ... " << ret << endl;
53
54 SelectionStatus status;
55 if( ret ) status.setStatus(SelectionStatus::PRESELECTION);
56 if(ctrl.debug) cout << "muon presel returning status : " << status.getStatus() << endl;
57 return status;
58 }
59
60
61 //--------------------------------------------------------------------------------------------------
62 bool isMuFO( const mithep::Muon * mu,
63 const mithep::Vertex & vtx,
64 const mithep::Array<mithep::PFCandidate> * pfCandidates )
65 //--------------------------------------------------------------------------------------------------
66 {
67 bool isgood=true;
68 float reliso = computePFMuonIso(mu,vtx,pfCandidates,0.3)/mu->Pt();
69
70 if( mu->Pt() < 5 ) isgood=false;
71 if ( fabs(mu->BestTrk()->DzCorrected(vtx)) > 0.2 ) isgood=false; // needed to remove cosmics in HF sample
72
73 //
74
75 // HF seems to not want tkquality, SS does
76 if( mu->Pt()>20 ) {
77 if ( !((mu->IsGlobalMuon()) && mu->NSegments()>0 ) ) isgood=false; //&& mu->nValidHits>0
78 } else {
79 if ( !((mu->IsGlobalMuon()) && mu->NSegments()>0 ) //&& //&& mu->nValidHits>0
80 && !( (mu->IsTrackerMuon()) && (mu->Quality().QualityMask().Mask() & muTools.kTMOneStationLoose)
81 )
82 ) isgood=false;
83 }
84
85
86 //
87 // for HF MC to agree w/ data, cannot put any Tk quality bits
88 // but to get to 5% closure in the SS test, we need at least Tk oneStationLoose
89
90 /*
91 && ( (mu->qualityBits & kTMOneStationLoose) ||
92 (mu->qualityBits & kTMLastStationOptimizedBarrelLowPtLoose) ||
93 (mu->qualityBits & kTMLastStationOptimizedLowPtLoose)
94 )
95 */
96
97
98 // comment for more stats for mu fake shape
99 if( mu->Pt() < 20 ) {
100 if( reliso > 3 ) isgood=false;
101 } else {
102 if( reliso > 5 ) isgood=false;
103 }
104
105 return isgood;
106 };
107
108 //--------------------------------------------------------------------------------------------------
109 SelectionStatus passSoftMuonSelection( ControlFlags &ctrl,
110 const mithep::Muon * mu ,
111 const mithep::Vertex & vtx)
112 //--------------------------------------------------------------------------------------------------
113 {
114
115 int level=0;
116 unsigned failmask=0x0;
117
118 // Use tracker track when available
119 const mithep::Track *muTrk=0;
120 if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
121 else if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
122 else if(mu->HasGlobalTrk()) { muTrk = mu->GlobalTrk(); }
123 assert(muTrk);
124
125 if(mu->TrackerTrk()->NHits() < 11 )
126 failmask |= (1<<level);
127 level++;
128
129 if(fabs(muTrk->D0Corrected(vtx)) > 0.2)
130 failmask |= (1<<level);
131 level++;
132
133 if(fabs(muTrk->DzCorrected(vtx)) > 0.1)
134 failmask |= (1<<level);
135 level++;
136
137 if(!(mu->IsTrackerMuon()))
138 failmask |= (1<<level);
139 level++;
140
141 // if(!(mu->Quality().QualityMask().Mask() & muTools.kTMLastStationAngTight))
142 if(!(mu->Quality().QualityMask().Mask() & muTools.kTMLastStationTight))
143 failmask |= (1<<level);
144 level++;
145
146 Double_t iso = (mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt())/mu->Pt();
147 if(mu->Pt()>20 && iso<0.1)
148 failmask |= (1<<level);
149
150 bool passtight = !(failmask);
151 bool passloose = !(failmask);
152 return SelectionStatus( passtight | passloose );
153
154 }
155
156 //--------------------------------------------------------------------------------------------------
157 //unsigned passMuonSelectionZZ( const mithep::TMuon * mu ) {
158 SelectionStatus passMuonSelectionZZ( ControlFlags &ctrl, const mithep::Muon * mu )
159 //--------------------------------------------------------------------------------------------------
160 {
161 int level=0;
162 unsigned failmask=0x0;
163
164 if(mu->Pt() < 5) {
165 failmask |= (1<<level);
166 }
167
168 level++;
169 if(fabs(mu->Eta()) > 2.4) {
170 failmask |= (1<<level);
171 }
172
173 level++;
174 if(!(mu->IsGlobalMuon())) {
175 failmask |= (1<<level);
176 }
177
178 level++;
179 if(mu->TrackerTrk()->NHits() < 10) {
180 failmask |= (1<<level);
181 }
182
183
184 level++;
185 if( (mu->IsoR03SumPt()/mu->Pt()) > 0.7 ) {
186 failmask |= (1<<level);
187 }
188
189 // if(muon->nPixHits < 1) return false;
190 // if(muon->muNchi2 > 10) return false;
191 // if(muon->nMatch < 2) return false;
192 // if(muon->nValidHits < 1) return false;
193 // if(muon->pterr/muon->pt > 0.1) return false;
194 // if(fabs(muon->dz) > 0.1) return false;
195
196 SelectionStatus status;
197 if( !failmask ) status.orStatus(SelectionStatus::LOOSEID);
198 if( !failmask ) status.orStatus(SelectionStatus::TIGHTID);
199 return status;
200
201
202
203 };
204
205
206
207
208 //
209 // Kevin's WW selection
210 //
211 //--------------------------------------------------------------------------------------------------
212 SelectionStatus passMuonSelection( ControlFlags &ctrl,
213 const mithep::Muon * mu,
214 const mithep::Vertex & vtx )
215 //--------------------------------------------------------------------------------------------------
216 {
217 SelectionStatus status;
218 int level=0;
219 unsigned failmask=0x0;
220
221 // Use tracker track when available
222 const mithep::Track *muTrk=0;
223 if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
224 else if(mu->HasGlobalTrk()) { muTrk = mu->GlobalTrk(); }
225 else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
226 assert(muTrk);
227
228
229
230 // 0x1
231 if(muTrk->Pt() < 5) {
232 failmask |= (1<<level);
233 }
234
235 // 0x2
236 level++;
237 if(fabs(muTrk->Eta()) > 2.4) {
238 failmask |= (1<<level);
239 }
240
241 // 0x4
242 level++;
243 if(muTrk->PtErr()/muTrk->Pt() > 0.1) {
244 failmask |= (1<<level);
245 }
246
247 // 0x8
248 level++;
249 if( fabs(muTrk->DzCorrected(vtx)) > 0.1 ) {
250 failmask |= (1<<level);
251 }
252
253 double muNchi2;
254 if(mu->HasStandaloneTrk()) { muNchi2 = mu->StandaloneTrk()->RChi2(); }
255 else if(mu->HasGlobalTrk()) { muNchi2 = mu->GlobalTrk()->RChi2(); }
256 else if(mu->HasTrackerTrk()) { muNchi2 = mu->TrackerTrk()->RChi2(); }
257
258 unsigned qualityBits = mu->Quality().QualityMask().Mask();
259
260 Bool_t isGlobal = (mu->IsGlobalMuon()) && (muNchi2 < 10) && (mu->NMatches() > 1) && (mu->NValidHits() > 0);
261 Bool_t isTracker = (mu->IsTrackerMuon() ) && (qualityBits & muTools.kTMLastStationTight);
262
263 // 0x10
264 level++;
265 if((!isGlobal && !isTracker) || !(mu->HasTrackerTrk())) {
266 failmask |= (1<<level);
267 status.setStatus(SelectionStatus::FAIL);
268 return status;
269 }
270
271 // 0x20
272 level++;
273 assert(mu->HasTrackerTrk());
274 assert(mu->TrackerTrk());
275 if( mu->TrackerTrk()->NHits() < 10 ) {
276 failmask |= (1<<level);
277 }
278
279 level++;
280 if(muTrk->NPixelHits() < 1) {
281 failmask |= (1<<level);
282 }
283
284
285 level++;
286 if(fabs(muTrk->D0Corrected(vtx))>0.02) {
287 failmask |= (1<<level);
288 }
289
290 /*
291 if(mu->pt>20) {
292 if(fabs(mu->d0)>0.02) {
293 failmask |= (1<<level);
294 }
295 } else {
296 if(fabs(mu->d0)>0.01) {
297 failmask |= (1<<level);
298 }
299 }
300 */
301
302 if( !failmask ) status.orStatus(SelectionStatus::LOOSEID);
303 if( !failmask ) status.orStatus(SelectionStatus::TIGHTID);
304 return status;
305
306
307 };
308
309
310 //--------------------------------------------------------------------------------------------------
311 SelectionStatus muonIDMVASelection(ControlFlags &ctrl,
312 const mithep::Muon *mu,
313 const mithep::Vertex & vtx )
314 //--------------------------------------------------------------------------------------------------
315 {
316 assert(mu->IsTrackerMuon()); // must be
317 assert(mu->HasTrackerTrk());
318 assert(mu->TrackerTrk());
319
320 /*
321 cerr << "TrackerTrk: " << mu->TrackerTrk() << endl;
322 flush(cerr);
323 cerr << "Rchi2: " << mu->TrackerTrk()->RChi2() << endl;
324 flush(cerr);
325
326 cerr << "calling MVa ... "<< endl;
327 flush(cerr);
328 */
329
330 double global_rchi2 = (mu->IsGlobalMuon()) ? mu->GlobalTrk()->RChi2() : 0.;
331
332 double mvaval = muIDMVA->MVAValue_ID(mu->Pt(),
333 mu->Eta(),
334 mu->IsGlobalMuon(),
335 mu->IsTrackerMuon(),
336 mu->TrackerTrk()->RChi2(),
337 global_rchi2,
338 (Double_t)(mu->NValidHits()),
339 (Double_t)(mu->TrackerTrk()->NHits()),
340 (Double_t)(mu->TrackerTrk()->NPixelHits()),
341 (Double_t)(mu->NMatches()),
342 mu->TrkKink(),
343 muTools.GetSegmentCompatability(mu),
344 muTools.GetCaloCompatability(mu,true,true),
345 mu->HadEnergy(),
346 mu->EmEnergy(),
347 mu->HadS9Energy(),
348 mu->EmS9Energy(),
349 (Bool_t)(ctrl.debug) );
350
351
352 SelectionStatus status;
353 bool pass = false;
354
355 if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
356 && fabs(mu->Eta()) < 1.5 && mu->Pt() < 10 && mvaval >= MUON_IDMVA_CUT_BIN0) pass = true;
357 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
358 && fabs(mu->Eta()) < 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN1) pass = true;
359 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
360 && fabs(mu->Eta()) > 1.5 && mu->Pt() < 10 && mvaval >= MUON_IDMVA_CUT_BIN2) pass = true;
361 else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
362 && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_IDMVA_CUT_BIN3) pass = true;
363 else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon()
364 && (mu->Quality().QualityMask().Mask() & mithep::MuonQuality::AllArbitrated) && mvaval >= MUON_IDMVA_CUT_BIN4)
365 pass = true;
366
367
368 if( pass ) {
369 status.orStatus(SelectionStatus::LOOSEID);
370 status.orStatus(SelectionStatus::TIGHTID);
371 }
372
373 return status;
374 }
375
376 //--------------------------------------------------------------------------------------------------
377 void initMuonIDMVA()
378 //--------------------------------------------------------------------------------------------------
379 {
380 muIDMVA = new mithep::MuonIDMVA();
381 vector<string> weightFiles;
382 weightFiles.push_back("../MitPhysics/data/MuonIDMVAWeights/MuonIDMVA_BDTG_barrel_lowpt_V2.weights.xml");
383 weightFiles.push_back("../MitPhysics/data/MuonIDMVAWeights/MuonIDMVA_BDTG_barrel_highpt_V2.weights.xml");
384 weightFiles.push_back("../MitPhysics/data/MuonIDMVAWeights/MuonIDMVA_BDTG_endcap_lowpt_V2.weights.xml");
385 weightFiles.push_back("../MitPhysics/data/MuonIDMVAWeights/MuonIDMVA_BDTG_endcap_highpt_V2.weights.xml");
386 weightFiles.push_back("../MitPhysics/data/MuonIDMVAWeights/MuonIDMVA_BDTG_tracker_V2.weights.xml");
387 weightFiles.push_back("../MitPhysics/data/MuonIDMVAWeights/MuonIDMVA_BDTG_global_V2.weights.xml");
388 muIDMVA->Initialize( "MuonIDMVA",
389 mithep::MuonIDMVA::kIDV0,
390 kTRUE, weightFiles);
391 }