ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.57
Committed: Tue Oct 11 17:00:24 2011 UTC (13 years, 6 months ago) by ceballos
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025c, Mit_025b, Mit_025a
Changes since 1.56: +11 -3 lines
Log Message:
adding new id

File Contents

# Content
1 // $Id: MuonIDMod.cc,v 1.56 2011/09/16 14:09:18 ceballos Exp $
2
3 #include "MitPhysics/Mods/interface/MuonIDMod.h"
4 #include "MitCommon/MathTools/interface/MathUtils.h"
5 #include "MitAna/DataTree/interface/MuonFwd.h"
6 #include "MitAna/DataTree/interface/ElectronFwd.h"
7 #include "MitAna/DataTree/interface/VertexCol.h"
8 #include "MitPhysics/Init/interface/ModNames.h"
9
10 using namespace mithep;
11
12 ClassImp(mithep::MuonIDMod)
13
14 //--------------------------------------------------------------------------------------------------
15 MuonIDMod::MuonIDMod(const char *name, const char *title) :
16 BaseMod(name,title),
17 fMuonBranchName(Names::gkMuonBrn),
18 fCleanMuonsName(ModNames::gkCleanMuonsName),
19 fNonIsolatedMuonsName("random"),
20 fNonIsolatedElectronsName("random"),
21 fVertexName(ModNames::gkGoodVertexesName),
22 fBeamSpotName(Names::gkBeamSpotBrn),
23 fTrackName(Names::gkTrackBrn),
24 fPFCandidatesName(Names::gkPFCandidatesBrn),
25 fMuonIDType("WWMuIdV3"),
26 fMuonIsoType("PFIso"),
27 fMuonClassType("Global"),
28 fTrackIsolationCut(3.0),
29 fCaloIsolationCut(3.0),
30 fCombIsolationCut(0.15),
31 fCombRelativeIsolationCut(0.15),
32 fPFIsolationCut(-1.0),
33 fMuonPtMin(10),
34 fApplyD0Cut(kTRUE),
35 fApplyDZCut(kTRUE),
36 fD0Cut(0.020),
37 fDZCut(0.10),
38 fWhichVertex(-1),
39 fEtaCut(2.4),
40 fMuIDType(kIdUndef),
41 fMuIsoType(kIsoUndef),
42 fMuClassType(kClassUndef),
43 fMuons(0),
44 fVertices(0),
45 fBeamSpot(0),
46 fTracks(0),
47 fPFCandidates(0),
48 fIntRadius(0.0),
49 fNonIsolatedMuons(0),
50 fNonIsolatedElectrons(0),
51 fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
52 fPileupEnergyDensity(0)
53 {
54 // Constructor.
55 }
56
57 //--------------------------------------------------------------------------------------------------
58 void MuonIDMod::Process()
59 {
60 // Process entries of the tree.
61
62 if(fMuIsoType != kPFIsoNoL) {
63 LoadEventObject(fMuonBranchName, fMuons);
64 }
65 else {
66 fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
67 }
68 LoadEventObject(fBeamSpotName, fBeamSpot);
69 LoadEventObject(fTrackName, fTracks);
70 LoadEventObject(fPFCandidatesName, fPFCandidates);
71 if(fMuIsoType == kTrackCaloSliding || fMuIsoType == kCombinedRelativeConeAreaCorrected) {
72 LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
73 }
74 MuonOArr *CleanMuons = new MuonOArr;
75 CleanMuons->SetName(fCleanMuonsName);
76
77 fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
78
79 for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
80 const Muon *mu = fMuons->At(i);
81
82 Bool_t pass = kFALSE;
83 Double_t pt = 0; // make sure pt is taken from the correct track!
84 Double_t eta = 0; // make sure eta is taken from the correct track!
85 switch (fMuClassType) {
86 case kAll:
87 pass = kTRUE;
88 if (mu->HasTrk()) {
89 pt = mu->Pt();
90 eta = TMath::Abs(mu->Eta());
91 }
92 break;
93 case kGlobal:
94 pass = mu->HasGlobalTrk() && mu->IsTrackerMuon();
95 if (pass && mu->TrackerTrk()) {
96 pt = mu->TrackerTrk()->Pt();
97 eta = TMath::Abs(mu->TrackerTrk()->Eta());
98 }
99 else {
100 pt = mu->Pt();
101 eta = TMath::Abs(mu->Eta());
102 }
103 break;
104 case kGlobalTracker:
105 pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
106 (mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
107 (mu->IsTrackerMuon() &&
108 mu->Quality().Quality(MuonQuality::TMLastStationTight));
109 if (pass) {
110 pt = mu->TrackerTrk()->Pt();
111 eta = TMath::Abs(mu->TrackerTrk()->Eta());
112 }
113 else {
114 pt = mu->Pt();
115 eta = TMath::Abs(mu->Eta());
116 }
117 break;
118 case kSta:
119 pass = mu->HasStandaloneTrk();
120 if (pass) {
121 pt = mu->StandaloneTrk()->Pt();
122 eta = TMath::Abs(mu->StandaloneTrk()->Eta());
123 }
124 break;
125 case kTrackerMuon:
126 pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
127 mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
128 if (pass) {
129 pt = mu->TrackerTrk()->Pt();
130 eta = TMath::Abs(mu->TrackerTrk()->Eta());
131 }
132 break;
133 case kCaloMuon:
134 pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
135 if (pass) {
136 pt = mu->TrackerTrk()->Pt();
137 eta = TMath::Abs(mu->TrackerTrk()->Eta());
138 }
139 break;
140 case kTrackerBased:
141 pass = mu->HasTrackerTrk();
142 if (pass) {
143 pt = mu->TrackerTrk()->Pt();
144 eta = TMath::Abs(mu->TrackerTrk()->Eta());
145 }
146 break;
147 default:
148 break;
149 }
150
151 if (!pass)
152 continue;
153
154 if (pt <= fMuonPtMin)
155 continue;
156
157 if (eta >= fEtaCut)
158 continue;
159
160 Double_t RChi2 = 0.0;
161 if (mu->HasGlobalTrk()) {
162 RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
163 }
164 else if(mu->BestTrk() != 0){
165 RChi2 = mu->BestTrk()->Chi2()/mu->BestTrk()->Ndof();
166 }
167 Bool_t idpass = kFALSE;
168 switch (fMuIDType) {
169 case kWMuId:
170 idpass = mu->BestTrk() != 0 &&
171 mu->BestTrk()->NHits() > 10 &&
172 RChi2 < 10.0 &&
173 (mu->NSegments() > 1 || mu->NMatches() > 1) &&
174 mu->BestTrk()->NPixelHits() > 0 &&
175 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
176 break;
177 case kZMuId:
178 idpass = mu->BestTrk() != 0 &&
179 mu->BestTrk()->NHits() > 10 &&
180 (mu->NSegments() > 1 || mu->NMatches() > 1) &&
181 mu->BestTrk()->NPixelHits() > 0 &&
182 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
183 break;
184 case kLoose:
185 idpass = mu->BestTrk() != 0 &&
186 mu->Quality().Quality(MuonQuality::TMOneStationLoose) &&
187 mu->Quality().Quality(MuonQuality::TM2DCompatibilityLoose) &&
188 mu->BestTrk()->NHits() > 10 &&
189 RChi2 < 10.0 &&
190 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
191 break;
192 case kTight:
193 idpass = mu->BestTrk() != 0 &&
194 mu->Quality().Quality(MuonQuality::TMOneStationTight) &&
195 mu->Quality().Quality(MuonQuality::TM2DCompatibilityTight) &&
196 mu->BestTrk()->NHits() > 10 &&
197 RChi2 < 10.0 &&
198 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
199 break;
200 case kWWMuIdV1:
201 idpass = mu->BestTrk() != 0 &&
202 mu->BestTrk()->NHits() > 10 &&
203 mu->BestTrk()->NPixelHits() > 0 &&
204 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
205 RChi2 < 10.0 &&
206 (mu->NSegments() > 1 || mu->NMatches() > 1) &&
207 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
208 break;
209 case kWWMuIdV2:
210 idpass = mu->BestTrk() != 0 &&
211 mu->BestTrk()->NHits() > 10 &&
212 mu->BestTrk()->NPixelHits() > 0 &&
213 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
214 break;
215 case kWWMuIdV3:
216 idpass = mu->BestTrk() != 0 &&
217 mu->BestTrk()->NHits() > 10 &&
218 mu->BestTrk()->NPixelHits() > 0 &&
219 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
220 mu->TrkKink() < 20.0;
221 break;
222 case kNoId:
223 idpass = kTRUE;
224 break;
225 default:
226 break;
227 }
228
229 if (!idpass)
230 continue;
231
232 Bool_t isocut = kFALSE;
233 switch (fMuIsoType) {
234 case kTrackCalo:
235 isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
236 (mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
237 break;
238 case kTrackCaloCombined:
239 isocut = (1.0 * mu->IsoR03SumPt() +
240 1.0 * mu->IsoR03EmEt() +
241 1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
242 break;
243 case kTrackCaloSliding:
244 {
245 const PileupEnergyDensity *rho = fPileupEnergyDensity->At(0);
246 Double_t totalIso = mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
247 // trick to change the signal region cut
248 double theIsoCut = fCombIsolationCut;
249 if(theIsoCut < 0.20){
250 if(mu->Pt() > 20.0) theIsoCut = 0.15;
251 else theIsoCut = 0.10;
252 }
253 if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
254 }
255 break;
256 case kTrackCaloSlidingNoCorrection:
257 {
258 Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
259 1.0 * mu->IsoR03EmEt() +
260 1.0 * mu->IsoR03HadEt();
261 // trick to change the signal region cut
262 double theIsoCut = fCombIsolationCut;
263 if(theIsoCut < 0.20){
264 if(mu->Pt() > 20.0) theIsoCut = 0.15;
265 else theIsoCut = 0.10;
266 }
267 if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
268 }
269 break;
270 case kCombinedRelativeConeAreaCorrected:
271 {
272 const PileupEnergyDensity *rho = fPileupEnergyDensity->At(0);
273 Double_t totalIso = mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3 ;
274 double theIsoCut = fCombRelativeIsolationCut;
275 if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
276 }
277 break;
278 case kPFIso:
279 {
280 Double_t pfIsoCutValue = 9999;
281 if(fPFIsolationCut > 0){
282 pfIsoCutValue = fPFIsolationCut;
283 } else {
284 if (mu->AbsEta() < 1.479) {
285 if (mu->Pt() > 20) {
286 pfIsoCutValue = 0.13;
287 } else {
288 pfIsoCutValue = 0.06;
289 }
290 } else {
291 if (mu->Pt() > 20) {
292 pfIsoCutValue = 0.09;
293 } else {
294 pfIsoCutValue = 0.05;
295 }
296 }
297 }
298 Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
299 if (totalIso < (mu->Pt()*pfIsoCutValue) )
300 isocut = kTRUE;
301 }
302 break;
303 case kPFIsoNoL:
304 {
305 fNonIsolatedMuons = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
306 fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
307
308 Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
309 if(beta == 0) beta = 1.0;
310 Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), fNonIsolatedMuons, fNonIsolatedElectrons, 0.2, 1.0, 0.4, 0.0, 3);
311 if (totalIso < (mu->Pt()*fPFIsolationCut) )
312 isocut = kTRUE;
313 }
314 break;
315 case kNoIso:
316 isocut = kTRUE;
317 break;
318 case kCustomIso:
319 default:
320 break;
321 }
322
323 if (isocut == kFALSE)
324 continue;
325
326 // apply d0 cut
327 if (fApplyD0Cut) {
328 Bool_t passD0cut = kTRUE;
329 if(fD0Cut < 0.05) { // trick to change the signal region cut
330 if (mu->Pt() > 20.0) fD0Cut = 0.02;
331 else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
332 }
333 if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
334 else passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
335 if (!passD0cut)
336 continue;
337 }
338
339 // apply dz cut
340 if (fApplyDZCut) {
341 Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
342 if (!passDZcut)
343 continue;
344 }
345
346 // add good muon
347 CleanMuons->Add(mu);
348 }
349
350 // sort according to pt
351 CleanMuons->Sort();
352
353 // add objects for other modules to use
354 AddObjThisEvt(CleanMuons);
355 }
356
357 //--------------------------------------------------------------------------------------------------
358 void MuonIDMod::SlaveBegin()
359 {
360 // Run startup code on the computer (slave) doing the actual analysis. Here,
361 // we just request the muon collection branch.
362
363 // In this case we cannot have a branch
364 if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
365 ReqEventObject(fMuonBranchName, fMuons, kTRUE);
366 }
367 ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
368 ReqEventObject(fTrackName, fTracks, kTRUE);
369 ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
370 if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
371 || fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0) {
372 ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
373 }
374
375 if (fMuonIDType.CompareTo("WMuId") == 0)
376 fMuIDType = kWMuId;
377 else if (fMuonIDType.CompareTo("ZMuId") == 0)
378 fMuIDType = kZMuId;
379 else if (fMuonIDType.CompareTo("Tight") == 0)
380 fMuIDType = kTight;
381 else if (fMuonIDType.CompareTo("Loose") == 0)
382 fMuIDType = kLoose;
383 else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
384 fMuIDType = kWWMuIdV1;
385 else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
386 fMuIDType = kWWMuIdV2;
387 else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
388 fMuIDType = kWWMuIdV3;
389 else if (fMuonIDType.CompareTo("NoId") == 0)
390 fMuIDType = kNoId;
391 else if (fMuonIDType.CompareTo("Custom") == 0) {
392 fMuIDType = kCustomId;
393 SendError(kWarning, "SlaveBegin",
394 "Custom muon identification is not yet implemented.");
395 } else {
396 SendError(kAbortAnalysis, "SlaveBegin",
397 "The specified muon identification %s is not defined.",
398 fMuonIDType.Data());
399 return;
400 }
401
402 if (fMuonIsoType.CompareTo("TrackCalo") == 0)
403 fMuIsoType = kTrackCalo;
404 else if (fMuonIsoType.CompareTo("TrackCaloCombined") == 0)
405 fMuIsoType = kTrackCaloCombined;
406 else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
407 fMuIsoType = kTrackCaloSliding;
408 else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
409 fMuIsoType = kTrackCaloSlidingNoCorrection;
410 else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
411 fMuIsoType = kCombinedRelativeConeAreaCorrected;
412 else if (fMuonIsoType.CompareTo("PFIso") == 0)
413 fMuIsoType = kPFIso;
414 else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
415 fMuIsoType = kPFIsoNoL;
416 else if (fMuonIsoType.CompareTo("NoIso") == 0)
417 fMuIsoType = kNoIso;
418 else if (fMuonIsoType.CompareTo("Custom") == 0) {
419 fMuIsoType = kCustomIso;
420 SendError(kWarning, "SlaveBegin",
421 "Custom muon isolation is not yet implemented.");
422 } else {
423 SendError(kAbortAnalysis, "SlaveBegin",
424 "The specified muon isolation %s is not defined.",
425 fMuonIsoType.Data());
426 return;
427 }
428
429 if (fMuonClassType.CompareTo("All") == 0)
430 fMuClassType = kAll;
431 else if (fMuonClassType.CompareTo("Global") == 0)
432 fMuClassType = kGlobal;
433 else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
434 fMuClassType = kGlobalTracker;
435 else if (fMuonClassType.CompareTo("Standalone") == 0)
436 fMuClassType = kSta;
437 else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
438 fMuClassType = kTrackerMuon;
439 else if (fMuonClassType.CompareTo("CaloMuon") == 0)
440 fMuClassType = kCaloMuon;
441 else if (fMuonClassType.CompareTo("TrackerBased") == 0)
442 fMuClassType = kTrackerBased;
443 else {
444 SendError(kAbortAnalysis, "SlaveBegin",
445 "The specified muon class %s is not defined.",
446 fMuonClassType.Data());
447 return;
448 }
449 }