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 |
}
|