ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.88
Committed: Sat Feb 23 14:51:45 2013 UTC (12 years, 2 months ago) by mingyang
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_029c, Mit_029b, Mit_029a
Changes since 1.87: +13 -1 lines
Log Message:
*** empty log message ***

File Contents

# Content
1 // $Id: MuonIDMod.cc,v 1.87 2013/01/07 22:14:56 mingyang 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 fPrintMVADebugInfo(kFALSE),
18 fMuonBranchName(Names::gkMuonBrn),
19 fCleanMuonsName(ModNames::gkCleanMuonsName),
20 fNonIsolatedMuonsName("random"),
21 fNonIsolatedElectronsName("random"),
22 fVertexName(ModNames::gkGoodVertexesName),
23 fBeamSpotName(Names::gkBeamSpotBrn),
24 fTrackName(Names::gkTrackBrn),
25 fPFCandidatesName(Names::gkPFCandidatesBrn),
26 fPFNoPileUpName("PFNoPileUp"),
27 fPFPileUpName("PFPileUp"),
28 fMuonIDType("WWMuIdV3"),
29 fMuonIsoType("PFIso"),
30 fMuonClassType("Global"),
31 fTrackIsolationCut(3.0),
32 fCaloIsolationCut(3.0),
33 fCombIsolationCut(0.15),
34 fCombRelativeIsolationCut(0.15),
35 fPFIsolationCut(-999.0),
36 fMuonPtMin(10),
37 fApplyD0Cut(kTRUE),
38 fApplyDZCut(kTRUE),
39 fD0Cut(0.020),
40 fDZCut(0.10),
41 fWhichVertex(-1),
42 fEtaCut(2.4),
43 fMuIDType(kIdUndef),
44 fMuIsoType(kIsoUndef),
45 fMuClassType(kClassUndef),
46 fMuons(0),
47 fVertices(0),
48 fBeamSpot(0),
49 fTracks(0),
50 fPFCandidates(0),
51 fPFNoPileUpCands(0),
52 fPFPileUpCands(0),
53 fIntRadius(0.0),
54 fNonIsolatedMuons(0),
55 fNonIsolatedElectrons(0),
56 fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
57 fPileupEnergyDensity(0),
58 fMuonTools(0),
59 fMuonIDMVA(0),
60 fTheRhoType(RhoUtilities::DEFAULT),
61 fPVName(Names::gkPVBeamSpotBrn)
62 {
63 // Constructor.
64 }
65
66 //--------------------------------------------------------------------------------------------------
67 void MuonIDMod::Process()
68 {
69 // Process entries of the tree.
70
71 if(fCleanMuonsName.CompareTo("HggLeptonTagMuons") == 0 ){
72 LoadEventObject(fPVName,fVertices);
73 }
74 else{
75 fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
76 }
77
78 if(fMuIsoType != kPFIsoNoL) {
79 LoadEventObject(fMuonBranchName, fMuons);
80 }
81 else {
82 fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
83 }
84 LoadEventObject(fBeamSpotName, fBeamSpot);
85 LoadEventObject(fTrackName, fTracks);
86 LoadEventObject(fPFCandidatesName, fPFCandidates);
87 if(fMuIsoType == kTrackCaloSliding ||
88 fMuIsoType == kCombinedRelativeConeAreaCorrected ||
89 fMuIsoType == kPFIsoEffectiveAreaCorrected ||
90 fMuIsoType == kMVAIso_BDTG_IDIso ||
91 fMuIsoType == kIsoRingsV0_BDTG_Iso ||
92 fMuIsoType == kIsoDeltaR
93 ) {
94 LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
95 }
96 if(fMuIsoType == kPFRadialIso || fMuIsoType == kIsoDeltaR || fMuIsoType == kPFIsoBetaPUCorrected){
97 // Name is hardcoded, can be changed if someone feels to do it *** did it--Heng
98 fPFNoPileUpCands = GetObjThisEvt<PFCandidateCol>(fPFNoPileUpName);
99 fPFPileUpCands = GetObjThisEvt<PFCandidateCol>(fPFPileUpName);
100 }
101
102 MuonOArr *CleanMuons = new MuonOArr;
103 CleanMuons->SetName(fCleanMuonsName);
104
105 for (UInt_t i=0; i<fMuons->GetEntries() && fVertices->GetEntries() > 0 ; ++i) {
106 const Muon *mu = fMuons->At(i);
107
108 Bool_t pass = kFALSE;
109 Double_t pt = 0; // make sure pt is taken from the correct track!
110 Double_t eta = 0; // make sure eta is taken from the correct track!
111 switch (fMuClassType) {
112 case kAll:
113 pass = kTRUE;
114 if (mu->HasTrk()) {
115 pt = mu->Pt();
116 eta = TMath::Abs(mu->Eta());
117 }
118 break;
119 case kGlobal:
120 pass = mu->HasGlobalTrk() && mu->IsTrackerMuon();
121 if (pass && mu->TrackerTrk()) {
122 pt = mu->TrackerTrk()->Pt();
123 eta = TMath::Abs(mu->TrackerTrk()->Eta());
124 }
125 else {
126 pt = mu->Pt();
127 eta = TMath::Abs(mu->Eta());
128 }
129 break;
130 case kGlobalTracker:
131 pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
132 (mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
133 (mu->IsTrackerMuon() &&
134 mu->Quality().Quality(MuonQuality::TMLastStationTight));
135 if (pass) {
136 pt = mu->TrackerTrk()->Pt();
137 eta = TMath::Abs(mu->TrackerTrk()->Eta());
138 }
139 else {
140 pt = mu->Pt();
141 eta = TMath::Abs(mu->Eta());
142 }
143 break;
144 case kSta:
145 pass = mu->HasStandaloneTrk();
146 if (pass) {
147 pt = mu->StandaloneTrk()->Pt();
148 eta = TMath::Abs(mu->StandaloneTrk()->Eta());
149 }
150 break;
151 case kTrackerMuon:
152 pass = mu->HasTrackerTrk() && mu->IsTrackerMuon() &&
153 mu->Quality().Quality(MuonQuality::TrackerMuonArbitrated);
154 if (pass) {
155 pt = mu->TrackerTrk()->Pt();
156 eta = TMath::Abs(mu->TrackerTrk()->Eta());
157 }
158 break;
159 case kCaloMuon:
160 pass = mu->HasTrackerTrk() && mu->IsCaloMuon();
161 if (pass) {
162 pt = mu->TrackerTrk()->Pt();
163 eta = TMath::Abs(mu->TrackerTrk()->Eta());
164 }
165 break;
166 case kTrackerBased:
167 pass = mu->HasTrackerTrk();
168 if (pass) {
169 pt = mu->TrackerTrk()->Pt();
170 eta = TMath::Abs(mu->TrackerTrk()->Eta());
171 }
172 break;
173 case kGlobalOnly:
174 pass = mu->HasGlobalTrk();
175 if (pass && mu->TrackerTrk()) {
176 pt = mu->TrackerTrk()->Pt();
177 eta = TMath::Abs(mu->TrackerTrk()->Eta());
178 }
179 else {
180 pt = mu->Pt();
181 eta = TMath::Abs(mu->Eta());
182 }
183 break;
184 default:
185 break;
186 }
187
188 if (!pass)
189 continue;
190
191 if (pt <= fMuonPtMin)
192 continue;
193
194 if (eta >= fEtaCut)
195 continue;
196
197 Double_t RChi2 = 0.0;
198 if (mu->HasGlobalTrk()) {
199 RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
200 }
201 else if(mu->BestTrk() != 0){
202 RChi2 = mu->BestTrk()->Chi2()/mu->BestTrk()->Ndof();
203 }
204 Bool_t idpass = kFALSE;
205 switch (fMuIDType) {
206 case kWMuId:
207 idpass = mu->BestTrk() != 0 &&
208 mu->BestTrk()->NHits() > 10 &&
209 RChi2 < 10.0 &&
210 (mu->NSegments() > 1 || mu->NMatches() > 1) &&
211 mu->BestTrk()->NPixelHits() > 0 &&
212 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
213 break;
214 case kZMuId:
215 idpass = mu->BestTrk() != 0 &&
216 mu->BestTrk()->NHits() > 10 &&
217 (mu->NSegments() > 1 || mu->NMatches() > 1) &&
218 mu->BestTrk()->NPixelHits() > 0 &&
219 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
220 break;
221 case kLoose:
222 idpass = mu->BestTrk() != 0 &&
223 mu->Quality().Quality(MuonQuality::TMOneStationLoose) &&
224 mu->Quality().Quality(MuonQuality::TM2DCompatibilityLoose) &&
225 mu->BestTrk()->NHits() > 10 &&
226 RChi2 < 10.0 &&
227 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
228 break;
229 case kTight:
230 idpass = mu->BestTrk() != 0 &&
231 mu->NTrkLayersHit() > 5 &&
232 mu->IsPFMuon() == kTRUE &&
233 mu->BestTrk()->NPixelHits() > 0 &&
234 RChi2 < 10.0;
235 break;
236 case kmuonPOG2012CutBasedIDTight:
237 idpass = mu->IsGlobalMuon() &&
238 mu->IsPFMuon() &&
239 mu->GlobalTrk()->RChi2() < 10 &&
240 mu->NValidHits() != 0 &&
241 mu->NMatches() > 1 &&
242 mu->BestTrk()->NPixelHits() != 0 &&
243 mu->NTrkLayersHit() > 5;
244 break;
245 // 2012 WW analysis for 42x (there is no PFMuon link)
246 case kWWMuIdV1:
247 idpass = mu->BestTrk() != 0 &&
248 mu->NTrkLayersHit() > 5 &&
249 mu->BestTrk()->NPixelHits() > 0 &&
250 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
251 mu->TrkKink() < 20.0;
252 break;
253 // 2010 WW analysis
254 case kWWMuIdV2:
255 idpass = mu->BestTrk() != 0 &&
256 mu->BestTrk()->NHits() > 10 &&
257 mu->BestTrk()->NPixelHits() > 0 &&
258 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
259 break;
260 // 2011 WW analysis
261 case kWWMuIdV3:
262 idpass = mu->BestTrk() != 0 &&
263 mu->BestTrk()->NHits() > 10 &&
264 mu->BestTrk()->NPixelHits() > 0 &&
265 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
266 mu->TrkKink() < 20.0;
267 break;
268 // 2012 WW analysis
269 case kWWMuIdV4:
270 idpass = mu->BestTrk() != 0 &&
271 mu->NTrkLayersHit() > 5 &&
272 mu->IsPFMuon() == kTRUE &&
273 mu->BestTrk()->NPixelHits() > 0 &&
274 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
275 mu->TrkKink() < 20.0;
276 break;
277 case kMVAID_BDTG_IDIso:
278 {
279 Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
280 mu->BestTrk()->NHits() > 10 &&
281 mu->BestTrk()->NPixelHits() > 0 &&
282 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
283 MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
284 MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
285 mu->TrkKink() < 20.0
286 );
287 idpass = passDenominatorM2;
288 //only evaluate MVA if muon passes M2 denominator to save time
289 if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
290 }
291 break;
292 case kNoId:
293 idpass = kTRUE;
294 break;
295 default:
296 break;
297 }
298
299 if (!idpass)
300 continue;
301
302 Double_t Rho = 0.;
303 if(fPileupEnergyDensity){
304 const PileupEnergyDensity *rho = fPileupEnergyDensity->At(0);
305
306 switch (fTheRhoType) {
307 case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
308 Rho = rho->RhoLowEta();
309 break;
310 case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
311 Rho = rho->Rho();
312 break;
313 case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
314 Rho = rho->RhoRandomLowEta();
315 break;
316 case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
317 Rho = rho->RhoRandom();
318 break;
319 case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
320 Rho = rho->RhoKt6PFJets();
321 break;
322 default:
323 Rho = rho->Rho();
324 }
325
326 if ((TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho()))) Rho = 0.;
327 }
328
329 Bool_t isocut = kFALSE;
330 switch (fMuIsoType) {
331 case kTrackCalo:
332 isocut = (mu->IsoR03SumPt() < fTrackIsolationCut) &&
333 (mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
334 break;
335 case kTrackCaloCombined:
336 isocut = (1.0 * mu->IsoR03SumPt() +
337 1.0 * mu->IsoR03EmEt() +
338 1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
339 break;
340 case kTrackCaloSliding:
341 {
342 Double_t totalIso = mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
343 // trick to change the signal region cut
344 double theIsoCut = fCombIsolationCut;
345 if(theIsoCut < 0.20){
346 if(mu->Pt() > 20.0) theIsoCut = 0.15;
347 else theIsoCut = 0.10;
348 }
349 if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
350 }
351 break;
352 case kTrackCaloSlidingNoCorrection:
353 {
354 Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
355 1.0 * mu->IsoR03EmEt() +
356 1.0 * mu->IsoR03HadEt();
357 // trick to change the signal region cut
358 double theIsoCut = fCombIsolationCut;
359 if(theIsoCut < 0.20){
360 if(mu->Pt() > 20.0) theIsoCut = 0.15;
361 else theIsoCut = 0.10;
362 }
363 if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
364 }
365 break;
366 case kCombinedRelativeConeAreaCorrected:
367 {
368 //const PileupEnergyDensity *rho = fPileupEnergyDensity->At(0); // Fabian: made Rho customable
369 Double_t totalIso = mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
370 double theIsoCut = fCombRelativeIsolationCut;
371 if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
372 }
373 break;
374 case kCombinedRelativeEffectiveAreaCorrected:
375 {
376 Double_t tmpRho = Rho; // Fabian: made the Rho type customable.
377 //if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))
378 //tmpRho = fPileupEnergyDensity->At(0)->Rho();
379
380 isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
381 - tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
382 - tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
383 ) < (mu->Pt()* 0.40);
384 }
385 break;
386 case kPFIso:
387 {
388 Double_t pfIsoCutValue = 9999;
389 if(fPFIsolationCut > 0){
390 pfIsoCutValue = fPFIsolationCut;
391 } else {
392 if (mu->AbsEta() < 1.479) {
393 if (mu->Pt() > 20) {
394 pfIsoCutValue = 0.13;
395 } else {
396 pfIsoCutValue = 0.06;
397 }
398 } else {
399 if (mu->Pt() > 20) {
400 pfIsoCutValue = 0.09;
401 } else {
402 pfIsoCutValue = 0.05;
403 }
404 }
405 }
406 Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
407 if (totalIso < (mu->Pt()*pfIsoCutValue) )
408 isocut = kTRUE;
409 }
410 break;
411 case kPFRadialIso:
412 {
413 Double_t pfIsoCutValue = 9999;
414 if(fPFIsolationCut > 0){
415 pfIsoCutValue = fPFIsolationCut;
416 } else {
417 if (mu->Pt() > 20) {
418 pfIsoCutValue = 0.10;
419 } else {
420 pfIsoCutValue = 0.05;
421 }
422 }
423 Double_t totalIso = IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
424 if (totalIso < (mu->Pt()*pfIsoCutValue) )
425 isocut = kTRUE;
426 }
427 break;
428 case kPFIsoBetaPUCorrected:
429 {
430 Double_t pfIsoCutValue = 9999;
431 if(fPFIsolationCut > 0){
432 pfIsoCutValue = fPFIsolationCut;
433 } else {
434 if (mu->Pt() > 20) {
435 pfIsoCutValue = 0.2;
436 } else {
437 pfIsoCutValue = 0.2;
438 }
439 }
440 Double_t totalIso = IsolationTools::BetaMwithPUCorrection(fPFNoPileUpCands, fPFPileUpCands, mu, 0.4);
441
442 if (totalIso < (mu->Pt()*pfIsoCutValue) )
443 isocut = kTRUE;
444 }
445 break;
446 case kPFIsoEffectiveAreaCorrected:
447 {
448 Double_t pfIsoCutValue = 9999;
449 if(fPFIsolationCut > 0){
450 pfIsoCutValue = fPFIsolationCut;
451 } else {
452 pfIsoCutValue = fPFIsolationCut; //leave it like this for now
453 }
454 Double_t EffectiveAreaCorrectedPFIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
455 - Rho * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
456 //- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta()); // Fabian: made Rho-type customable
457 isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
458 break;
459 }
460
461
462 case kPFIsoNoL:
463 {
464 fNonIsolatedMuons = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
465 fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
466
467 Double_t pfIsoCutValue = 9999;
468 if(fPFIsolationCut > 0){
469 pfIsoCutValue = fPFIsolationCut;
470 } else {
471 if (mu->AbsEta() < 1.479) {
472 if (mu->Pt() > 20) {
473 pfIsoCutValue = 0.13;
474 } else {
475 pfIsoCutValue = 0.06;
476 }
477 } else {
478 if (mu->Pt() > 20) {
479 pfIsoCutValue = 0.09;
480 } else {
481 pfIsoCutValue = 0.05;
482 }
483 }
484 }
485 Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
486 if (totalIso < (mu->Pt()*pfIsoCutValue) )
487 isocut = kTRUE;
488 }
489 break;
490 case kMVAIso_BDTG_IDIso:
491 {
492
493 Double_t totalIso = IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
494 isocut = (totalIso < (mu->Pt()*0.4));
495
496 }
497 break;
498 case kIsoRingsV0_BDTG_Iso:
499 {
500
501 isocut = PassMuonIsoRingsV0_BDTG_Iso(mu, fVertices->At(0), fPileupEnergyDensity);
502
503 }
504 break;
505 case kIsoDeltaR:
506 {
507
508 isocut = PassMuonIsoDeltaR(mu, fVertices->At(0), fPileupEnergyDensity);
509
510 }
511 break;
512 case kNoIso:
513 isocut = kTRUE;
514 break;
515 case kCustomIso:
516 default:
517 break;
518 }
519
520 if (isocut == kFALSE)
521 continue;
522
523 // apply d0 cut
524 if (fApplyD0Cut) {
525 Bool_t passD0cut = kTRUE;
526 if(fD0Cut < 0.05) { // trick to change the signal region cut
527 if (mu->Pt() > 20.0) fD0Cut = 0.02;
528 else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
529 }
530 if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
531 else passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
532 if (!passD0cut)
533 continue;
534 }
535
536 // apply dz cut
537 if (fApplyDZCut) {
538 Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
539 if (!passDZcut)
540 continue;
541 }
542
543 // add good muon
544 CleanMuons->Add(mu);
545 }
546
547 // sort according to pt
548 CleanMuons->Sort();
549
550 // add objects for other modules to use
551 AddObjThisEvt(CleanMuons);
552 }
553
554 //--------------------------------------------------------------------------------------------------
555 void MuonIDMod::SlaveBegin()
556 {
557 // Run startup code on the computer (slave) doing the actual analysis. Here,
558 // we just request the muon collection branch.
559
560 if(fCleanMuonsName.CompareTo("HggLeptonTagMuons") == 0 ){
561 ReqEventObject(fPVName,fVertices,true);
562 }
563
564 // In this case we cannot have a branch
565 if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
566 ReqEventObject(fMuonBranchName, fMuons, kTRUE);
567 }
568 ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
569 ReqEventObject(fTrackName, fTracks, kTRUE);
570 ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
571 if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
572 || fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
573 || fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
574 || fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
575 || fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
576 || fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0
577 || fMuonIsoType.CompareTo("IsoDeltaR") == 0
578 ) {
579 ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
580 }
581
582
583 if (fMuonIDType.CompareTo("WMuId") == 0)
584 fMuIDType = kWMuId;
585 else if (fMuonIDType.CompareTo("ZMuId") == 0)
586 fMuIDType = kZMuId;
587 else if (fMuonIDType.CompareTo("Tight") == 0)
588 fMuIDType = kTight;
589 else if (fMuonIDType.CompareTo("muonPOG2012CutBasedIDTight") == 0)
590 fMuIDType = kmuonPOG2012CutBasedIDTight;
591 else if (fMuonIDType.CompareTo("Loose") == 0)
592 fMuIDType = kLoose;
593 else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
594 fMuIDType = kWWMuIdV1;
595 else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
596 fMuIDType = kWWMuIdV2;
597 else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
598 fMuIDType = kWWMuIdV3;
599 else if (fMuonIDType.CompareTo("WWMuIdV4") == 0)
600 fMuIDType = kWWMuIdV4;
601 else if (fMuonIDType.CompareTo("NoId") == 0)
602 fMuIDType = kNoId;
603 else if (fMuonIDType.CompareTo("Custom") == 0) {
604 fMuIDType = kCustomId;
605 SendError(kWarning, "SlaveBegin",
606 "Custom muon identification is not yet implemented.");
607 } else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
608 fMuIDType = kMVAID_BDTG_IDIso;
609 } else {
610 SendError(kAbortAnalysis, "SlaveBegin",
611 "The specified muon identification %s is not defined.",
612 fMuonIDType.Data());
613 return;
614 }
615
616 if (fMuonIsoType.CompareTo("TrackCalo") == 0)
617 fMuIsoType = kTrackCalo;
618 else if (fMuonIsoType.CompareTo("TrackCaloCombined") == 0)
619 fMuIsoType = kTrackCaloCombined;
620 else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
621 fMuIsoType = kTrackCaloSliding;
622 else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
623 fMuIsoType = kTrackCaloSlidingNoCorrection;
624 else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
625 fMuIsoType = kCombinedRelativeConeAreaCorrected;
626 else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
627 fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
628 else if (fMuonIsoType.CompareTo("PFIso") == 0)
629 fMuIsoType = kPFIso;
630 else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
631 fMuIsoType = kPFRadialIso;
632 else if (fMuonIsoType.CompareTo("PFIsoBetaPUCorrected") == 0)
633 fMuIsoType = kPFIsoBetaPUCorrected;
634 else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
635 fMuIsoType = kPFIsoEffectiveAreaCorrected;
636 else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
637 fMuIsoType = kPFIsoNoL;
638 else if (fMuonIsoType.CompareTo("NoIso") == 0)
639 fMuIsoType = kNoIso;
640 else if (fMuonIsoType.CompareTo("Custom") == 0) {
641 fMuIsoType = kCustomIso;
642 SendError(kWarning, "SlaveBegin",
643 "Custom muon isolation is not yet implemented.");
644 } else if (fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0) {
645 fMuIsoType = kMVAIso_BDTG_IDIso;
646 } else if (fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0) {
647 fMuIsoType = kIsoRingsV0_BDTG_Iso;
648 } else if (fMuonIsoType.CompareTo("IsoDeltaR") == 0) {
649 fMuIsoType = kIsoDeltaR;
650 } else {
651 SendError(kAbortAnalysis, "SlaveBegin",
652 "The specified muon isolation %s is not defined.",
653 fMuonIsoType.Data());
654 return;
655 }
656
657 if (fMuonClassType.CompareTo("All") == 0)
658 fMuClassType = kAll;
659 else if (fMuonClassType.CompareTo("Global") == 0)
660 fMuClassType = kGlobal;
661 else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
662 fMuClassType = kGlobalTracker;
663 else if (fMuonClassType.CompareTo("Standalone") == 0)
664 fMuClassType = kSta;
665 else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
666 fMuClassType = kTrackerMuon;
667 else if (fMuonClassType.CompareTo("CaloMuon") == 0)
668 fMuClassType = kCaloMuon;
669 else if (fMuonClassType.CompareTo("TrackerBased") == 0)
670 fMuClassType = kTrackerBased;
671 else if (fMuonClassType.CompareTo("GlobalOnly") == 0)
672 fMuClassType = kGlobalOnly;
673 else {
674 SendError(kAbortAnalysis, "SlaveBegin",
675 "The specified muon class %s is not defined.",
676 fMuonClassType.Data());
677 return;
678 }
679
680
681 //If we use MVA ID, need to load MVA weights
682 if (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
683 fMuonTools = new MuonTools();
684 fMuonIDMVA = new MuonIDMVA();
685 fMuonIDMVA->Initialize("BDTG method",
686 string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin0_IDIsoCombined_BDTG.weights.xml"))),
687 string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin0_IDIsoCombined_BDTG.weights.xml"))),
688 string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin1_IDIsoCombined_BDTG.weights.xml"))),
689 string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin1_IDIsoCombined_BDTG.weights.xml"))),
690 string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin2_IDIsoCombined_BDTG.weights.xml"))),
691 string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin2_IDIsoCombined_BDTG.weights.xml"))),
692 MuonIDMVA::kIDIsoCombinedDetIso,
693 fTheRhoType);
694 }
695 else if(fMuIsoType == kIsoRingsV0_BDTG_Iso) {
696 std::vector<std::string> muonidiso_weightfiles;
697 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml"))));
698 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml"))));
699 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml"))));
700 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml"))));
701 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml"))));
702 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml"))));
703 fMuonTools = new MuonTools();
704 fMuonIDMVA = new MuonIDMVA();
705 fMuonIDMVA->Initialize("MuonIso_BDTG_IsoRings",
706 MuonIDMVA::kIsoRingsV0,
707 kTRUE,
708 muonidiso_weightfiles,
709 fTheRhoType);
710 }
711 else if(fMuIsoType == kIsoDeltaR) {
712 std::vector<std::string> muonidiso_weightfiles;
713 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LB_BDT.weights.xml"))));
714 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LE_BDT.weights.xml"))));
715 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HB_BDT.weights.xml"))));
716 muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HE_BDT.weights.xml"))));
717 fMuonTools = new MuonTools();
718 fMuonIDMVA = new MuonIDMVA();
719 fMuonIDMVA->Initialize("muonHZZ2012IsoDRMVA",
720 MuonIDMVA::kIsoDeltaR,
721 kTRUE,
722 muonidiso_weightfiles,
723 fTheRhoType);
724 }
725
726 }
727
728
729 //--------------------------------------------------------------------------------------------------
730 Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
731 const PileupEnergyDensityCol *PileupEnergyDensity) const
732 {
733
734 const Track *muTrk=0;
735 if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
736 else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
737
738 Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
739
740 Int_t subdet = 0;
741 if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
742 else subdet = 1;
743 Int_t ptBin = 0;
744 if (muTrk->Pt() > 14.5) ptBin = 1;
745 if (muTrk->Pt() > 20.0) ptBin = 2;
746
747 Int_t MVABin = -1;
748 if (subdet == 0 && ptBin == 0) MVABin = 0;
749 else if (subdet == 1 && ptBin == 0) MVABin = 1;
750 else if (subdet == 0 && ptBin == 1) MVABin = 2;
751 else if (subdet == 1 && ptBin == 1) MVABin = 3;
752 else if (subdet == 0 && ptBin == 2) MVABin = 4;
753 else if (subdet == 1 && ptBin == 2) MVABin = 5;
754
755 Double_t MVACut = -999;
756 if (MVABin == 0) MVACut = -0.5618;
757 else if (MVABin == 1) MVACut = -0.3002;
758 else if (MVABin == 2) MVACut = -0.4642;
759 else if (MVABin == 3) MVACut = -0.2478;
760 else if (MVABin == 4) MVACut = 0.1706;
761 else if (MVABin == 5) MVACut = 0.8146;
762
763 if (MVAValue > MVACut) return kTRUE;
764 return kFALSE;
765 }
766
767 //--------------------------------------------------------------------------------------------------
768 Bool_t MuonIDMod::PassMuonIsoRingsV0_BDTG_Iso(const Muon *mu, const Vertex *vertex,
769 const PileupEnergyDensityCol *PileupEnergyDensity) const
770 {
771
772 Bool_t isDebug = kFALSE;
773 const Track *muTrk=0;
774 if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
775 else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
776
777 ElectronOArr *tempElectrons = new ElectronOArr;
778 MuonOArr *tempMuons = new MuonOArr;
779 Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,
780 PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,isDebug);
781 delete tempElectrons;
782 delete tempMuons;
783
784 Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
785
786 Double_t MVACut = -1.0;
787 Double_t eta = mu->AbsEta();
788 if (mu->Pt() < 20 && eta < 1.479) MVACut = 0.86;
789 else if(mu->Pt() < 20 && eta >= 1.479) MVACut = 0.82;
790 else if(mu->Pt() >= 20 && eta < 1.479) MVACut = 0.82;
791 else if(mu->Pt() >= 20 && eta >= 1.479) MVACut = 0.86;
792
793 if(fPFIsolationCut > -1.0) MVACut = fPFIsolationCut;
794
795 if(isDebug == kTRUE){
796 printf("PassMuonIsoRingsV0_BDTG_IsoDebug: %d, pt, eta = %f, %f, rho = %f(%f) : RingsMVA = %f, bin: %d\n",
797 GetEventHeader()->EvtNum(),mu->Pt(), mu->Eta(),
798 fPileupEnergyDensity->At(0)->Rho(),fPileupEnergyDensity->At(0)->RhoKt6PFJets(),MVAValue,MVABin);
799 }
800
801 if (MVAValue > MVACut) return kTRUE;
802 return kFALSE;
803 }
804
805 //--------------------------------------------------------------------------------------------------
806 Bool_t MuonIDMod::PassMuonIsoDeltaR(const Muon *mu, const Vertex *vertex,
807 const PileupEnergyDensityCol *PileupEnergyDensity) const
808 {
809
810 const Track *muTrk=0;
811 if(mu->HasTrackerTrk()) { muTrk = mu->TrackerTrk(); }
812 else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
813
814 ElectronOArr *tempElectrons = new ElectronOArr;
815 MuonOArr *tempMuons = new MuonOArr;
816 Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFNoPileUpCands,
817 PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,kFALSE);
818 delete tempElectrons;
819 delete tempMuons;
820
821 Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
822
823 Double_t MVACut = -999;
824 if (MVABin == 0) MVACut = 0.000;
825 else if (MVABin == 1) MVACut = 0.000;
826 else if (MVABin == 2) MVACut = 0.000;
827 else if (MVABin == 3) MVACut = 0.000;
828
829 if (MVAValue > MVACut) return kTRUE;
830 return kFALSE;
831 }
832
833 //--------------------------------------------------------------------------------------------------
834 void MuonIDMod::Terminate()
835 {
836 // Run finishing code on the computer (slave) that did the analysis
837 delete fMuonIDMVA;
838
839 delete fMuonTools;
840 }