ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
Revision: 1.93
Committed: Fri Oct 18 14:10:13 2013 UTC (11 years, 6 months ago) by ceballos
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.92: +2 -2 lines
Error occurred while calculating annotation data.
Log Message:
minor fix

File Contents

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