ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.55
Committed: Fri Jul 22 15:45:38 2011 UTC (13 years, 9 months ago) by sixie
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_024b, Mit_025pre1, Mit_024a, Mit_024
Changes since 1.54: +4 -3 lines
Log Message: 
fix bug

File Contents

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