ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.56
Committed: Fri Sep 16 14:09:18 2011 UTC (13 years, 7 months ago) by ceballos
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_025, Mit_025pre2
Changes since 1.55: +3 -2 lines
Log Message: 
fixing iso

File Contents

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