ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Mods/src/MuonIDMod.cc
(Generate patch)

Comparing UserCode/MitPhysics/Mods/src/MuonIDMod.cc (file contents):
Revision 1.49 by ceballos, Thu May 12 23:15:48 2011 UTC vs.
Revision 1.87 by mingyang, Mon Jan 7 22:14:56 2013 UTC

# Line 14 | Line 14 | ClassImp(mithep::MuonIDMod)
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"),  
# Line 22 | Line 23 | ClassImp(mithep::MuonIDMod)
23    fBeamSpotName(Names::gkBeamSpotBrn),
24    fTrackName(Names::gkTrackBrn),
25    fPFCandidatesName(Names::gkPFCandidatesBrn),
26 <  fMuonIDType("WWMuId"),
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(-1.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.20),
40 >  fDZCut(0.10),
41    fWhichVertex(-1),
42    fEtaCut(2.4),
43    fMuIDType(kIdUndef),
# Line 43 | Line 48 | ClassImp(mithep::MuonIDMod)
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)
57 >  fPileupEnergyDensity(0),
58 >  fMuonTools(0),
59 >  fMuonIDMVA(0),
60 >  fTheRhoType(RhoUtilities::DEFAULT),
61 >  fPVName(Names::gkPVBeamSpotBrn)
62   {
63    // Constructor.
64   }
# Line 56 | Line 68 | 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    }
# Line 65 | Line 84 | void MuonIDMod::Process()
84    LoadEventObject(fBeamSpotName, fBeamSpot);
85    LoadEventObject(fTrackName, fTracks);
86    LoadEventObject(fPFCandidatesName, fPFCandidates);
87 <  if(fMuIsoType == kTrackCaloSliding) {
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 <  fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
75 <
76 <  for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
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;
# Line 98 | Line 127 | void MuonIDMod::Process()
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) {
# Line 127 | Line 170 | void MuonIDMod::Process()
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      }
# Line 173 | Line 227 | void MuonIDMod::Process()
227                   mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
228          break;
229        case kTight:
230 <        idpass = mu->BestTrk() !=  0 &&
231 <                 mu->Quality().Quality(MuonQuality::TMOneStationTight) &&
232 <                 mu->Quality().Quality(MuonQuality::TM2DCompatibilityTight) &&
233 <                 mu->BestTrk()->NHits() > 10 &&
234 <                 RChi2 < 10.0 &&
181 <                 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
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 kWWMuId:
236 >      // 2012 WW analysis for 42x (there is no PFMuon link)
237 >      case kWWMuIdV1:
238 >        idpass = mu->BestTrk() != 0 &&
239 >                 mu->NTrkLayersHit() > 5 &&
240 >                 mu->BestTrk()->NPixelHits() > 0 &&
241 >                 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
242 >                 mu->TrkKink() < 20.0;
243 >        break;
244 >      // 2010 WW analysis
245 >      case kWWMuIdV2:
246          idpass = mu->BestTrk() != 0 &&
247                   mu->BestTrk()->NHits() > 10 &&
186                 RChi2 < 10.0 &&
187                (mu->NSegments() > 1 || mu->NMatches() > 1) &&
248                   mu->BestTrk()->NPixelHits() > 0 &&
189                 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight) &&
249                   mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
250          break;
251 +      // 2011 WW analysis
252 +      case kWWMuIdV3:
253 +        idpass = mu->BestTrk() != 0 &&
254 +                 mu->BestTrk()->NHits() > 10 &&
255 +                 mu->BestTrk()->NPixelHits() > 0 &&
256 +                 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
257 +                 mu->TrkKink() < 20.0;
258 +        break;
259 +      // 2012 WW analysis
260 +      case kWWMuIdV4:
261 +        idpass = mu->BestTrk() != 0 &&
262 +                 mu->NTrkLayersHit() > 5 &&
263 +                 mu->IsPFMuon() == kTRUE &&
264 +                 mu->BestTrk()->NPixelHits() > 0 &&
265 +                 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
266 +                 mu->TrkKink() < 20.0;
267 +        break;
268 +      case kMVAID_BDTG_IDIso:
269 +        {
270 +          Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
271 +                                      mu->BestTrk()->NHits() > 10 &&
272 +                                      mu->BestTrk()->NPixelHits() > 0 &&
273 +                                      mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
274 +                                      MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
275 +                                      MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
276 +                                      mu->TrkKink() < 20.0
277 +            );  
278 +          idpass =  passDenominatorM2;
279 +          //only evaluate MVA if muon passes M2 denominator to save time
280 +          if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
281 +        }
282 +        break;
283        case kNoId:
284          idpass = kTRUE;
285          break;
# Line 199 | Line 290 | void MuonIDMod::Process()
290      if (!idpass)
291        continue;
292  
293 +    Double_t Rho = 0.;
294 +    if(fPileupEnergyDensity){
295 +      const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
296 +
297 +      switch (fTheRhoType) {
298 +      case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
299 +        Rho = rho->RhoLowEta();
300 +        break;
301 +      case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
302 +        Rho = rho->Rho();
303 +        break;
304 +      case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
305 +        Rho = rho->RhoRandomLowEta();
306 +        break;
307 +      case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
308 +        Rho = rho->RhoRandom();
309 +        break;
310 +      case RhoUtilities::CMS_RHO_RHOKT6PFJETS:
311 +        Rho = rho->RhoKt6PFJets();
312 +        break;
313 +      default:
314 +        Rho = rho->Rho();
315 +      }
316 +
317 +      if ((TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho()))) Rho = 0.;
318 +    }
319 +
320      Bool_t isocut = kFALSE;
321      switch (fMuIsoType) {
322        case kTrackCalo:
# Line 212 | Line 330 | void MuonIDMod::Process()
330          break;
331        case kTrackCaloSliding:
332          {
333 <          const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
216 <          Double_t totalIso =  mu->IsoR03SumPt() + TMath::Max(mu->IsoR03EmEt() + mu->IsoR03HadEt() - rho->Rho() * TMath::Pi() * 0.3 * 0.3, 0.0);
333 >          Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
334            // trick to change the signal region cut
335            double theIsoCut = fCombIsolationCut;
336            if(theIsoCut < 0.20){
# Line 223 | Line 340 | void MuonIDMod::Process()
340            if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
341          }
342          break;
343 <      case kTrackCaloSlidingNoCorrection:
343 >    case kTrackCaloSlidingNoCorrection:
344          {
345            Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
346                                 1.0 * mu->IsoR03EmEt()  +
# Line 237 | Line 354 | void MuonIDMod::Process()
354            if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
355          }
356          break;
357 <      case kPFIso:
358 <        {
359 <          Double_t pfIsoCutValue = 9999;
360 <          if(fCombIsolationCut > 0){
361 <            pfIsoCutValue = fCombIsolationCut;
362 <          } else {
363 <            if (mu->Pt() > 20) {
364 <              pfIsoCutValue = 0.17;
365 <            } else {
366 <              pfIsoCutValue = 0.13;
357 >    case kCombinedRelativeConeAreaCorrected:    
358 >      {          
359 >        //const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0); // Fabian: made Rho customable          
360 >        Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;          
361 >        double theIsoCut = fCombRelativeIsolationCut;    
362 >        if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;    
363 >      }          
364 >      break;    
365 >    case kCombinedRelativeEffectiveAreaCorrected:        
366 >      {          
367 >        Double_t tmpRho = Rho;   // Fabian: made the Rho type customable.        
368 >        //if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))    
369 >        //tmpRho = fPileupEnergyDensity->At(0)->Rho();  
370 >        
371 >        isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()      
372 >                   -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())      
373 >                   -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())    
374 >                   ) < (mu->Pt()* 0.40);        
375 >      }          
376 >      break;
377 >    case kPFIso:
378 >      {
379 >        Double_t pfIsoCutValue = 9999;
380 >        if(fPFIsolationCut > 0){
381 >          pfIsoCutValue = fPFIsolationCut;
382 >        } else {
383 >          if (mu->AbsEta() < 1.479) {
384 >            if (mu->Pt() > 20) {
385 >              pfIsoCutValue = 0.13;
386 >            } else {
387 >              pfIsoCutValue = 0.06;
388 >            }
389 >          } else {
390 >            if (mu->Pt() > 20) {
391 >              pfIsoCutValue = 0.09;
392 >            } else {
393 >              pfIsoCutValue = 0.05;
394 >            }
395 >          }
396 >        }
397 >        Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
398 >        if (totalIso < (mu->Pt()*pfIsoCutValue) )
399 >          isocut = kTRUE;
400 >      }
401 >      break;
402 >    case kPFRadialIso:
403 >      {
404 >        Double_t pfIsoCutValue = 9999;
405 >        if(fPFIsolationCut > 0){
406 >          pfIsoCutValue = fPFIsolationCut;
407 >        } else {
408 >          if (mu->Pt() > 20) {
409 >            pfIsoCutValue = 0.10;
410 >          } else {
411 >              pfIsoCutValue = 0.05;
412              }
413            }
414 <          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0),0.2, 1.0, 0.4, 0.0);
414 >          Double_t totalIso =  IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
415            if (totalIso < (mu->Pt()*pfIsoCutValue) )
416              isocut = kTRUE;
417 +      }
418 +      break;
419 +    case kPFIsoBetaPUCorrected:
420 +      {
421 +        Double_t pfIsoCutValue = 9999;
422 +        if(fPFIsolationCut > 0){
423 +          pfIsoCutValue = fPFIsolationCut;
424 +        } else {
425 +          if (mu->Pt() > 20) {
426 +            pfIsoCutValue = 0.2;
427 +          } else {
428 +            pfIsoCutValue = 0.2;
429 +          }
430          }
431 <        break;
432 <      case kPFIsoNoL:
431 >        Double_t totalIso =  IsolationTools::BetaMwithPUCorrection(fPFNoPileUpCands, fPFPileUpCands, mu, 0.4);
432 >        if (totalIso < (mu->Pt()*pfIsoCutValue) )
433 >          isocut = kTRUE;
434 >      }
435 >      break;
436 >    case kPFIsoEffectiveAreaCorrected:  
437 >      {          
438 >        Double_t pfIsoCutValue = 9999;  
439 >        if(fPFIsolationCut > 0){        
440 >          pfIsoCutValue = fPFIsolationCut;      
441 >        } else {        
442 >          pfIsoCutValue = fPFIsolationCut; //leave it like this for now          
443 >        }        
444 >        Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
445 >          - Rho * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
446 >        //- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());  // Fabian: made Rho-type customable      
447 >        isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);      
448 >        break;  
449 >      }
450 >      
451 >      
452 >    case kPFIsoNoL:
453          {
454            fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
455            fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
456  
457 <          Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
458 <          if(beta == 0) beta = 1.0;
459 <          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), fNonIsolatedMuons, fNonIsolatedElectrons, 0.2, 1.0, 0.4, 0.0, 3);
460 <          if (totalIso < (mu->Pt()*fCombIsolationCut) )
457 >          Double_t pfIsoCutValue = 9999;
458 >          if(fPFIsolationCut > 0){
459 >            pfIsoCutValue = fPFIsolationCut;
460 >          } else {
461 >            if (mu->AbsEta() < 1.479) {
462 >              if (mu->Pt() > 20) {
463 >                pfIsoCutValue = 0.13;
464 >              } else {
465 >                pfIsoCutValue = 0.06;
466 >              }
467 >            } else {
468 >              if (mu->Pt() > 20) {
469 >                pfIsoCutValue = 0.09;
470 >              } else {
471 >                pfIsoCutValue = 0.05;
472 >              }
473 >            }
474 >          }
475 >          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
476 >          if (totalIso < (mu->Pt()*pfIsoCutValue) )
477              isocut = kTRUE;
478          }
479          break;
480 +      case kMVAIso_BDTG_IDIso:
481 +      {
482 +
483 +        Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
484 +        isocut = (totalIso < (mu->Pt()*0.4));
485 +
486 +      }
487 +        break;
488 +      case kIsoRingsV0_BDTG_Iso:
489 +      {
490 +        
491 +        isocut = PassMuonIsoRingsV0_BDTG_Iso(mu, fVertices->At(0), fPileupEnergyDensity);
492 +
493 +      }
494 +        break;
495 +      case kIsoDeltaR:
496 +      {
497 +        
498 +        isocut = PassMuonIsoDeltaR(mu, fVertices->At(0), fPileupEnergyDensity);
499 +
500 +      }
501 +        break;
502        case kNoIso:
503          isocut = kTRUE;
504          break;
# Line 314 | Line 547 | void MuonIDMod::SlaveBegin()
547    // Run startup code on the computer (slave) doing the actual analysis. Here,
548    // we just request the muon collection branch.
549  
550 +  if(fCleanMuonsName.CompareTo("HggLeptonTagMuons") == 0 ){
551 +    ReqEventObject(fPVName,fVertices,true);
552 +  }
553 +
554     // In this case we cannot have a branch
555    if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
556      ReqEventObject(fMuonBranchName, fMuons, kTRUE);
# Line 321 | Line 558 | void MuonIDMod::SlaveBegin()
558    ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
559    ReqEventObject(fTrackName, fTracks, kTRUE);
560    ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
561 <  if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0) {
561 >  if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
562 >      || fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0        
563 >      || fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0  
564 >      || fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
565 >      || fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
566 >      || fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0
567 >      || fMuonIsoType.CompareTo("IsoDeltaR") == 0
568 >      ) {
569      ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
570    }
571  
572 +
573    if (fMuonIDType.CompareTo("WMuId") == 0)
574      fMuIDType = kWMuId;
575    else if (fMuonIDType.CompareTo("ZMuId") == 0)
# Line 333 | Line 578 | void MuonIDMod::SlaveBegin()
578      fMuIDType = kTight;
579    else if (fMuonIDType.CompareTo("Loose") == 0)
580      fMuIDType = kLoose;
581 <  else if (fMuonIDType.CompareTo("WWMuId") == 0)
582 <    fMuIDType = kWWMuId;
581 >  else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
582 >    fMuIDType = kWWMuIdV1;
583 >  else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
584 >    fMuIDType = kWWMuIdV2;
585 >  else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
586 >    fMuIDType = kWWMuIdV3;
587 >  else if (fMuonIDType.CompareTo("WWMuIdV4") == 0)
588 >    fMuIDType = kWWMuIdV4;
589    else if (fMuonIDType.CompareTo("NoId") == 0)
590      fMuIDType = kNoId;
591    else if (fMuonIDType.CompareTo("Custom") == 0) {
592      fMuIDType = kCustomId;
593      SendError(kWarning, "SlaveBegin",
594                "Custom muon identification is not yet implemented.");
595 +  } else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
596 +    fMuIDType = kMVAID_BDTG_IDIso;
597    } else {
598      SendError(kAbortAnalysis, "SlaveBegin",
599                "The specified muon identification %s is not defined.",
# Line 356 | Line 609 | void MuonIDMod::SlaveBegin()
609      fMuIsoType = kTrackCaloSliding;
610    else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
611      fMuIsoType = kTrackCaloSlidingNoCorrection;
612 +  else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)    
613 +    fMuIsoType = kCombinedRelativeConeAreaCorrected;    
614 +  else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)        
615 +    fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
616    else if (fMuonIsoType.CompareTo("PFIso") == 0)
617      fMuIsoType = kPFIso;
618 +  else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
619 +    fMuIsoType = kPFRadialIso;
620 +  else if (fMuonIsoType.CompareTo("PFIsoBetaPUCorrected") == 0)
621 +    fMuIsoType = kPFIsoBetaPUCorrected;
622 +  else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)  
623 +    fMuIsoType = kPFIsoEffectiveAreaCorrected;
624    else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
625      fMuIsoType = kPFIsoNoL;
626    else if (fMuonIsoType.CompareTo("NoIso") == 0)
# Line 366 | Line 629 | void MuonIDMod::SlaveBegin()
629      fMuIsoType = kCustomIso;
630      SendError(kWarning, "SlaveBegin",
631                "Custom muon isolation is not yet implemented.");
632 +  } else if (fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0) {
633 +    fMuIsoType = kMVAIso_BDTG_IDIso;
634 +  } else if (fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0) {
635 +    fMuIsoType = kIsoRingsV0_BDTG_Iso;
636 +  } else if (fMuonIsoType.CompareTo("IsoDeltaR") == 0) {
637 +    fMuIsoType = kIsoDeltaR;
638    } else {
639      SendError(kAbortAnalysis, "SlaveBegin",
640                "The specified muon isolation %s is not defined.",
# Line 377 | Line 646 | void MuonIDMod::SlaveBegin()
646      fMuClassType = kAll;
647    else if (fMuonClassType.CompareTo("Global") == 0)
648      fMuClassType = kGlobal;
649 +  else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
650 +    fMuClassType = kGlobalTracker;
651    else if (fMuonClassType.CompareTo("Standalone") == 0)
652      fMuClassType = kSta;
653    else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
# Line 385 | Line 656 | void MuonIDMod::SlaveBegin()
656      fMuClassType = kCaloMuon;
657    else if (fMuonClassType.CompareTo("TrackerBased") == 0)
658      fMuClassType = kTrackerBased;
659 +  else if (fMuonClassType.CompareTo("GlobalOnly") == 0)
660 +    fMuClassType = kGlobalOnly;
661    else {
662      SendError(kAbortAnalysis, "SlaveBegin",
663                "The specified muon class %s is not defined.",
664                fMuonClassType.Data());
665      return;
666    }
667 +
668 +
669 +  //If we use MVA ID, need to load MVA weights
670 +  if     (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
671 +    fMuonTools = new MuonTools();
672 +    fMuonIDMVA = new MuonIDMVA();
673 +    fMuonIDMVA->Initialize("BDTG method",
674 +                           string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin0_IDIsoCombined_BDTG.weights.xml"))),
675 +                           string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin0_IDIsoCombined_BDTG.weights.xml"))),
676 +                           string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin1_IDIsoCombined_BDTG.weights.xml"))),
677 +                           string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin1_IDIsoCombined_BDTG.weights.xml"))),
678 +                           string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/BarrelPtBin2_IDIsoCombined_BDTG.weights.xml"))),
679 +                           string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/EndcapPtBin2_IDIsoCombined_BDTG.weights.xml"))),
680 +                           MuonIDMVA::kIDIsoCombinedDetIso,
681 +                           fTheRhoType);
682 +  }
683 +  else if(fMuIsoType == kIsoRingsV0_BDTG_Iso) {
684 +    std::vector<std::string> muonidiso_weightfiles;
685 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml"))));
686 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml"))));
687 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml"))));
688 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml"))));
689 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml"))));
690 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml"))));
691 +    fMuonTools = new MuonTools();
692 +    fMuonIDMVA = new MuonIDMVA();
693 +    fMuonIDMVA->Initialize("MuonIso_BDTG_IsoRings",
694 +                       MuonIDMVA::kIsoRingsV0,
695 +                       kTRUE,
696 +                       muonidiso_weightfiles,
697 +                       fTheRhoType);
698 +  }
699 +  else if(fMuIsoType == kIsoDeltaR) {
700 +    std::vector<std::string> muonidiso_weightfiles;
701 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LB_BDT.weights.xml"))));
702 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LE_BDT.weights.xml"))));
703 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HB_BDT.weights.xml"))));
704 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HE_BDT.weights.xml"))));
705 +    fMuonTools = new MuonTools();
706 +    fMuonIDMVA = new MuonIDMVA();
707 +    fMuonIDMVA->Initialize("muonHZZ2012IsoDRMVA",
708 +                       MuonIDMVA::kIsoDeltaR,
709 +                       kTRUE,
710 +                       muonidiso_weightfiles,
711 +                       fTheRhoType);
712 +  }
713 +
714 + }
715 +
716 +
717 + //--------------------------------------------------------------------------------------------------
718 + Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
719 +                                         const PileupEnergyDensityCol *PileupEnergyDensity) const
720 + {
721 +
722 +  const Track *muTrk=0;
723 +  if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
724 +  else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
725 +  
726 +  Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
727 +
728 +  Int_t subdet = 0;
729 +  if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
730 +  else subdet = 1;
731 +  Int_t ptBin = 0;
732 +  if (muTrk->Pt() > 14.5) ptBin = 1;
733 +  if (muTrk->Pt() > 20.0) ptBin = 2;
734 +
735 +  Int_t MVABin = -1;
736 +  if      (subdet == 0 && ptBin == 0) MVABin = 0;
737 +  else if (subdet == 1 && ptBin == 0) MVABin = 1;
738 +  else if (subdet == 0 && ptBin == 1) MVABin = 2;
739 +  else if (subdet == 1 && ptBin == 1) MVABin = 3;
740 +  else if (subdet == 0 && ptBin == 2) MVABin = 4;
741 +  else if (subdet == 1 && ptBin == 2) MVABin = 5;
742 +
743 +  Double_t MVACut = -999;
744 +  if      (MVABin == 0) MVACut = -0.5618;
745 +  else if (MVABin == 1) MVACut = -0.3002;
746 +  else if (MVABin == 2) MVACut = -0.4642;
747 +  else if (MVABin == 3) MVACut = -0.2478;
748 +  else if (MVABin == 4) MVACut =  0.1706;
749 +  else if (MVABin == 5) MVACut =  0.8146;
750 +
751 +  if (MVAValue > MVACut) return kTRUE;
752 +  return kFALSE;
753 + }
754 +
755 + //--------------------------------------------------------------------------------------------------
756 + Bool_t MuonIDMod::PassMuonIsoRingsV0_BDTG_Iso(const Muon *mu, const Vertex *vertex,
757 +                                              const PileupEnergyDensityCol *PileupEnergyDensity) const
758 + {
759 +
760 +  Bool_t isDebug = kFALSE;
761 +  const Track *muTrk=0;
762 +  if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
763 +  else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
764 +  
765 +  ElectronOArr *tempElectrons = new  ElectronOArr;
766 +  MuonOArr     *tempMuons     = new  MuonOArr;
767 +  Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,
768 +                      PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,isDebug);
769 +  delete tempElectrons;
770 +  delete tempMuons;
771 +
772 +  Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
773 +
774 +  Double_t MVACut = -1.0;
775 +  Double_t eta = mu->AbsEta();
776 +  if     (mu->Pt() <  20 && eta <  1.479) MVACut = 0.86;
777 +  else if(mu->Pt() <  20 && eta >= 1.479) MVACut = 0.82;
778 +  else if(mu->Pt() >= 20 && eta <  1.479) MVACut = 0.82;
779 +  else if(mu->Pt() >= 20 && eta >= 1.479) MVACut = 0.86;
780 +
781 +  if(fPFIsolationCut > -1.0) MVACut = fPFIsolationCut;
782 +
783 +  if(isDebug == kTRUE){
784 +    printf("PassMuonIsoRingsV0_BDTG_IsoDebug: %d, pt, eta = %f, %f, rho = %f(%f) : RingsMVA = %f, bin: %d\n",
785 +           GetEventHeader()->EvtNum(),mu->Pt(), mu->Eta(),
786 +           fPileupEnergyDensity->At(0)->Rho(),fPileupEnergyDensity->At(0)->RhoKt6PFJets(),MVAValue,MVABin);
787 +  }
788 +
789 +  if (MVAValue > MVACut) return kTRUE;
790 +  return kFALSE;
791 + }
792 +
793 + //--------------------------------------------------------------------------------------------------
794 + Bool_t MuonIDMod::PassMuonIsoDeltaR(const Muon *mu, const Vertex *vertex,
795 +                                    const PileupEnergyDensityCol *PileupEnergyDensity) const
796 + {
797 +
798 +  const Track *muTrk=0;
799 +  if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
800 +  else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
801 +  
802 +  ElectronOArr *tempElectrons = new  ElectronOArr;
803 +  MuonOArr     *tempMuons     = new  MuonOArr;
804 +  Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFNoPileUpCands,
805 +                      PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,kFALSE);
806 +  delete tempElectrons;
807 +  delete tempMuons;
808 +
809 +  Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
810 +
811 +  Double_t MVACut = -999;
812 +  if      (MVABin == 0) MVACut =  0.000;
813 +  else if (MVABin == 1) MVACut =  0.000;
814 +  else if (MVABin == 2) MVACut =  0.000;
815 +  else if (MVABin == 3) MVACut =  0.000;
816 +
817 +  if (MVAValue > MVACut) return kTRUE;
818 +  return kFALSE;
819 + }
820 +
821 + //--------------------------------------------------------------------------------------------------
822 + void MuonIDMod::Terminate()
823 + {
824 +  // Run finishing code on the computer (slave) that did the analysis
825 +  delete fMuonIDMVA;
826 +  
827 +  delete fMuonTools;
828   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines