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.46 by ceballos, Tue Apr 5 06:37:07 2011 UTC vs.
Revision 1.75 by ceballos, Fri May 4 21:47:03 2012 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"),
27 <  fMuonIsoType("TrackCaloSliding"),
26 >  fMuonIDType("WWMuIdV3"),
27 >  fMuonIsoType("PFIso"),
28    fMuonClassType("Global"),  
29    fTrackIsolationCut(3.0),
30    fCaloIsolationCut(3.0),
31    fCombIsolationCut(0.15),
32 +  fCombRelativeIsolationCut(0.15),
33 +  fPFIsolationCut(-1.0),
34    fMuonPtMin(10),
35    fApplyD0Cut(kTRUE),
36    fApplyDZCut(kTRUE),
37    fD0Cut(0.020),
38 <  fDZCut(0.20),
38 >  fDZCut(0.10),
39    fWhichVertex(-1),
40    fEtaCut(2.4),
41    fMuIDType(kIdUndef),
# Line 43 | Line 46 | ClassImp(mithep::MuonIDMod)
46    fBeamSpot(0),
47    fTracks(0),
48    fPFCandidates(0),
49 +  fPFNoPileUpCands(0),
50 +  fIntRadius(0.0),
51    fNonIsolatedMuons(0),
52    fNonIsolatedElectrons(0),
53    fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn),
54 <  fPileupEnergyDensity(0)
54 >  fPileupEnergyDensity(0),
55 >  fMuonTools(0),
56 >  fMuonIDMVA(0),
57 >  fMuonMVAWeights_Subdet0Pt10To14p5(""),
58 >  fMuonMVAWeights_Subdet1Pt10To14p5(""),
59 >  fMuonMVAWeights_Subdet0Pt14p5To20(""),
60 >  fMuonMVAWeights_Subdet1Pt14p5To20(""),
61 >  fMuonMVAWeights_Subdet0Pt20ToInf(""),
62 >  fMuonMVAWeights_Subdet1Pt20ToInf(""),
63 >  fTheRhoType(RhoUtilities::DEFAULT)
64   {
65    // Constructor.
66   }
# Line 65 | Line 79 | void MuonIDMod::Process()
79    LoadEventObject(fBeamSpotName, fBeamSpot);
80    LoadEventObject(fTrackName, fTracks);
81    LoadEventObject(fPFCandidatesName, fPFCandidates);
82 <  if(fMuIsoType == kTrackCaloSliding) {
82 >  if(fMuIsoType == kTrackCaloSliding ||
83 >     fMuIsoType == kCombinedRelativeConeAreaCorrected ||
84 >     fMuIsoType == kPFIsoEffectiveAreaCorrected ||
85 >     fMuIsoType == kMVAIso_BDTG_IDIso ||
86 >     fMuIsoType == kIsoRingsV0_BDTG_Iso ||
87 >     fMuIsoType == kIsoDeltaR
88 >    ) {
89      LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
90    }
91 +  if(fMuIsoType == kPFRadialIso || fMuIsoType == kIsoDeltaR){
92 +    // Name is hardcoded, can be changed if someone feels to do it
93 +    fPFNoPileUpCands = GetObjThisEvt<PFCandidateCol>("PFNoPileUp");    
94 +  }
95 +
96    MuonOArr *CleanMuons = new MuonOArr;
97    CleanMuons->SetName(fCleanMuonsName);
98  
99    fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
100  
101 <  for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
101 >  for (UInt_t i=0; i<fMuons->GetEntries() && fVertices->GetEntries() > 0 ; ++i) {
102      const Muon *mu = fMuons->At(i);
103  
104      Bool_t pass = kFALSE;
# Line 98 | Line 123 | void MuonIDMod::Process()
123            eta = TMath::Abs(mu->Eta());
124          }
125          break;
126 +      case kGlobalTracker:
127 +        pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
128 +               (mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
129 +               (mu->IsTrackerMuon() &&
130 +                mu->Quality().Quality(MuonQuality::TMLastStationTight));
131 +        if (pass) {
132 +          pt  = mu->TrackerTrk()->Pt();
133 +          eta = TMath::Abs(mu->TrackerTrk()->Eta());
134 +        }
135 +        else {
136 +          pt  = mu->Pt();
137 +          eta = TMath::Abs(mu->Eta());
138 +        }
139 +        break;
140        case kSta:
141          pass = mu->HasStandaloneTrk();
142          if (pass) {
# Line 140 | Line 179 | void MuonIDMod::Process()
179      if (eta >= fEtaCut)
180        continue;
181  
182 +
183 +    //***********************************************************************************************
184 +    //Debug Info For Lepton MVA
185 +    //***********************************************************************************************
186 +    if( fPrintMVADebugInfo &&
187 +        (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso)
188 +      ) {
189 +      cout << "Event: " << GetEventHeader()->RunNum() << " " << GetEventHeader()->LumiSec() << " "
190 +           << GetEventHeader()->EvtNum() << " : Rho = " << fPileupEnergyDensity->At(0)->Rho()
191 +           << " : Muon " << i << " "
192 +           << endl;
193 +      fMuonIDMVA->MVAValue(mu,fVertices->At(0),fMuonTools,fPFCandidates,fPileupEnergyDensity,kTRUE);
194 +    }
195 +    //***********************************************************************************************
196 +
197 +
198      Double_t RChi2 = 0.0;
199      if     (mu->HasGlobalTrk()) {
200        RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
# Line 180 | Line 235 | void MuonIDMod::Process()
235                   RChi2 < 10.0 &&
236                   mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
237          break;
238 <      case kWWMuId:
238 >      case kWWMuIdV1:
239          idpass = mu->BestTrk() != 0 &&
240                   mu->BestTrk()->NHits() > 10 &&
241 +                 mu->BestTrk()->NPixelHits() > 0 &&
242 +                 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
243                   RChi2 < 10.0 &&
244                  (mu->NSegments() > 1 || mu->NMatches() > 1) &&
245 +                 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
246 +        break;
247 +      case kWWMuIdV2:
248 +        idpass = mu->BestTrk() != 0 &&
249 +                 mu->BestTrk()->NHits() > 10 &&
250                   mu->BestTrk()->NPixelHits() > 0 &&
189                 mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight) &&
251                   mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
252          break;
253 +      case kWWMuIdV3:
254 +        idpass = mu->BestTrk() != 0 &&
255 +                 mu->BestTrk()->NHits() > 10 &&
256 +                 mu->BestTrk()->NPixelHits() > 0 &&
257 +                 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
258 +                 mu->TrkKink() < 20.0;
259 +        break;
260 +      case kMVAID_BDTG_IDIso:
261 +        {
262 +          Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
263 +                                      mu->BestTrk()->NHits() > 10 &&
264 +                                      mu->BestTrk()->NPixelHits() > 0 &&
265 +                                      mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
266 +                                      MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
267 +                                      MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
268 +                                      mu->TrkKink() < 20.0
269 +            );  
270 +          idpass =  passDenominatorM2;
271 +          //only evaluate MVA if muon passes M2 denominator to save time
272 +          if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
273 +        }
274 +        break;
275        case kNoId:
276          idpass = kTRUE;
277          break;
# Line 199 | Line 282 | void MuonIDMod::Process()
282      if (!idpass)
283        continue;
284  
285 +    Double_t Rho = 0.;
286 +    if(fPileupEnergyDensity){
287 +      const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
288 +
289 +      switch (fTheRhoType) {
290 +      case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
291 +        Rho = rho->RhoLowEta();
292 +        break;
293 +      case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
294 +        Rho = rho->Rho();
295 +        break;
296 +      case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
297 +        Rho = rho->RhoRandomLowEta();
298 +        break;
299 +      case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
300 +        Rho = rho->RhoRandom();
301 +        break;
302 +      default:
303 +        Rho = rho->Rho();
304 +      }
305 +
306 +      if ((TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho()))) Rho = 0.;
307 +    }
308 +
309      Bool_t isocut = kFALSE;
310      switch (fMuIsoType) {
311        case kTrackCalo:
# Line 212 | Line 319 | void MuonIDMod::Process()
319          break;
320        case kTrackCaloSliding:
321          {
322 <          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);
322 >          Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
323            // trick to change the signal region cut
324            double theIsoCut = fCombIsolationCut;
325            if(theIsoCut < 0.20){
# Line 237 | Line 343 | void MuonIDMod::Process()
343            if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
344          }
345          break;
346 +      case kCombinedRelativeConeAreaCorrected:
347 +        {
348 +          //const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0); // Fabian: made Rho customable
349 +          Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
350 +          double theIsoCut = fCombRelativeIsolationCut;
351 +          if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
352 +        }
353 +        break;          
354 +    case kCombinedRelativeEffectiveAreaCorrected:
355 +      {
356 +        Double_t tmpRho = Rho;   // Fabian: made the Rho type customable.
357 +        //if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))
358 +        //tmpRho = fPileupEnergyDensity->At(0)->Rho();
359 +        
360 +          isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
361 +                     -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
362 +                     -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
363 +            ) < (mu->Pt()* 0.40);
364 +        }
365 +        break;          
366        case kPFIso:
367          {
368 <          Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
369 <          if(beta == 0) beta = 1.0;
370 <          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 0, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
371 <          if (totalIso < (mu->Pt()*fCombIsolationCut) )
368 >          Double_t pfIsoCutValue = 9999;
369 >          if(fPFIsolationCut > 0){
370 >            pfIsoCutValue = fPFIsolationCut;
371 >          } else {
372 >            if (mu->AbsEta() < 1.479) {
373 >              if (mu->Pt() > 20) {
374 >                pfIsoCutValue = 0.13;
375 >              } else {
376 >                pfIsoCutValue = 0.06;
377 >              }
378 >            } else {
379 >              if (mu->Pt() > 20) {
380 >                pfIsoCutValue = 0.09;
381 >              } else {
382 >                pfIsoCutValue = 0.05;
383 >              }
384 >            }
385 >          }
386 >          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
387 >          if (totalIso < (mu->Pt()*pfIsoCutValue) )
388 >            isocut = kTRUE;
389 >        }
390 >        break;
391 >      case kPFRadialIso:
392 >        {
393 >          Double_t pfIsoCutValue = 9999;
394 >          if(fPFIsolationCut > 0){
395 >            pfIsoCutValue = fPFIsolationCut;
396 >          } else {
397 >            if (mu->Pt() > 20) {
398 >              pfIsoCutValue = 0.10;
399 >            } else {
400 >              pfIsoCutValue = 0.05;
401 >            }
402 >          }
403 >          Double_t totalIso =  IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
404 >          if (totalIso < (mu->Pt()*pfIsoCutValue) )
405              isocut = kTRUE;
406          }
407          break;
408 +      case kPFIsoEffectiveAreaCorrected:
409 +        {
410 +          Double_t pfIsoCutValue = 9999;
411 +          if(fPFIsolationCut > 0){
412 +            pfIsoCutValue = fPFIsolationCut;
413 +          } else {
414 +            pfIsoCutValue = fPFIsolationCut; //leave it like this for now
415 +          }
416 +          Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
417 +            - Rho * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
418 +          //- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());  // Fabian: made Rho-type customable
419 +          isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
420 +          break;
421 +        }
422        case kPFIsoNoL:
423          {
424            fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
425            fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
426  
427 <          Double_t beta = IsolationTools::BetaM(fTracks, mu, fVertices->At(0), 0.0, 0.2, 0.3, 0.02);
428 <          if(beta == 0) beta = 1.0;
429 <          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.2, 0.5, 0.3, 0.02, 3, beta, fNonIsolatedMuons, fNonIsolatedElectrons);
430 <          if (totalIso < (mu->Pt()*fCombIsolationCut) )
427 >          Double_t pfIsoCutValue = 9999;
428 >          if(fPFIsolationCut > 0){
429 >            pfIsoCutValue = fPFIsolationCut;
430 >          } else {
431 >            if (mu->AbsEta() < 1.479) {
432 >              if (mu->Pt() > 20) {
433 >                pfIsoCutValue = 0.13;
434 >              } else {
435 >                pfIsoCutValue = 0.06;
436 >              }
437 >            } else {
438 >              if (mu->Pt() > 20) {
439 >                pfIsoCutValue = 0.09;
440 >              } else {
441 >                pfIsoCutValue = 0.05;
442 >              }
443 >            }
444 >          }
445 >          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
446 >          if (totalIso < (mu->Pt()*pfIsoCutValue) )
447              isocut = kTRUE;
448          }
449          break;
450 +      case kMVAIso_BDTG_IDIso:
451 +      {
452 +
453 +        // **************************************************************************
454 +        // Don't use effective area correction denominator. Instead use the old one.
455 +        // **************************************************************************
456 +
457 +        //         Double_t tmpRho = 0;
458 +        //         if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || isinf(fPileupEnergyDensity->At(0)->Rho())))
459 +        //           tmpRho = fPileupEnergyDensity->At(0)->Rho();
460 +        
461 +        //         isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
462 +        //                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
463 +        //                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
464 +        //           ) < (mu->Pt()* 0.40);
465 +        
466 +        Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
467 +        isocut = (totalIso < (mu->Pt()*0.4));
468 +
469 +      }
470 +        break;
471 +      case kIsoRingsV0_BDTG_Iso:
472 +      {
473 +        
474 +        isocut = PassMuonIsoRingsV0_BDTG_Iso(mu, fVertices->At(0), fPileupEnergyDensity);
475 +
476 +      }
477 +        break;
478 +      case kIsoDeltaR:
479 +      {
480 +        
481 +        isocut = PassMuonIsoDeltaR(mu, fVertices->At(0), fPileupEnergyDensity);
482 +
483 +      }
484 +        break;
485        case kNoIso:
486          isocut = kTRUE;
487          break;
# Line 313 | Line 537 | void MuonIDMod::SlaveBegin()
537    ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
538    ReqEventObject(fTrackName, fTracks, kTRUE);
539    ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
540 <  if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0) {
540 >  if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
541 >      || fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
542 >      || fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
543 >      || fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
544 >      || fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
545 >      || fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0
546 >      || fMuonIsoType.CompareTo("IsoDeltaR") == 0
547 >    ) {
548      ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
549    }
550  
551 +
552    if (fMuonIDType.CompareTo("WMuId") == 0)
553      fMuIDType = kWMuId;
554    else if (fMuonIDType.CompareTo("ZMuId") == 0)
# Line 325 | Line 557 | void MuonIDMod::SlaveBegin()
557      fMuIDType = kTight;
558    else if (fMuonIDType.CompareTo("Loose") == 0)
559      fMuIDType = kLoose;
560 <  else if (fMuonIDType.CompareTo("WWMuId") == 0)
561 <    fMuIDType = kWWMuId;
560 >  else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
561 >    fMuIDType = kWWMuIdV1;
562 >  else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
563 >    fMuIDType = kWWMuIdV2;
564 >  else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
565 >    fMuIDType = kWWMuIdV3;
566    else if (fMuonIDType.CompareTo("NoId") == 0)
567      fMuIDType = kNoId;
568    else if (fMuonIDType.CompareTo("Custom") == 0) {
569      fMuIDType = kCustomId;
570      SendError(kWarning, "SlaveBegin",
571                "Custom muon identification is not yet implemented.");
572 +  } else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
573 +    fMuIDType = kMVAID_BDTG_IDIso;
574    } else {
575      SendError(kAbortAnalysis, "SlaveBegin",
576                "The specified muon identification %s is not defined.",
# Line 348 | Line 586 | void MuonIDMod::SlaveBegin()
586      fMuIsoType = kTrackCaloSliding;
587    else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
588      fMuIsoType = kTrackCaloSlidingNoCorrection;
589 +  else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
590 +    fMuIsoType = kCombinedRelativeConeAreaCorrected;
591 +  else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
592 +    fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
593    else if (fMuonIsoType.CompareTo("PFIso") == 0)
594      fMuIsoType = kPFIso;
595 +  else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
596 +    fMuIsoType = kPFRadialIso;
597 +  else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
598 +    fMuIsoType = kPFIsoEffectiveAreaCorrected;
599    else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
600      fMuIsoType = kPFIsoNoL;
601    else if (fMuonIsoType.CompareTo("NoIso") == 0)
# Line 358 | Line 604 | void MuonIDMod::SlaveBegin()
604      fMuIsoType = kCustomIso;
605      SendError(kWarning, "SlaveBegin",
606                "Custom muon isolation is not yet implemented.");
607 +  } else if (fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0) {
608 +    fMuIsoType = kMVAIso_BDTG_IDIso;
609 +  } else if (fMuonIsoType.CompareTo("IsoRingsV0_BDTG_Iso") == 0) {
610 +    fMuIsoType = kIsoRingsV0_BDTG_Iso;
611 +  } else if (fMuonIsoType.CompareTo("IsoDeltaR") == 0) {
612 +    fMuIsoType = kIsoDeltaR;
613    } else {
614      SendError(kAbortAnalysis, "SlaveBegin",
615                "The specified muon isolation %s is not defined.",
# Line 369 | Line 621 | void MuonIDMod::SlaveBegin()
621      fMuClassType = kAll;
622    else if (fMuonClassType.CompareTo("Global") == 0)
623      fMuClassType = kGlobal;
624 +  else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
625 +    fMuClassType = kGlobalTracker;
626    else if (fMuonClassType.CompareTo("Standalone") == 0)
627      fMuClassType = kSta;
628    else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
# Line 383 | Line 637 | void MuonIDMod::SlaveBegin()
637                fMuonClassType.Data());
638      return;
639    }
640 +
641 +
642 +  //If we use MVA ID, need to load MVA weights
643 +  if     (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
644 +    fMuonTools = new MuonTools();
645 +    fMuonIDMVA = new MuonIDMVA();
646 +    fMuonIDMVA->Initialize("BDTG method",
647 +                           fMuonMVAWeights_Subdet0Pt10To14p5,
648 +                           fMuonMVAWeights_Subdet1Pt10To14p5,
649 +                           fMuonMVAWeights_Subdet0Pt14p5To20,
650 +                           fMuonMVAWeights_Subdet1Pt14p5To20,
651 +                           fMuonMVAWeights_Subdet0Pt20ToInf,
652 +                           fMuonMVAWeights_Subdet1Pt20ToInf,
653 +                           MuonIDMVA::kIDIsoCombinedDetIso);
654 +  }
655 +  else if(fMuIsoType == kIsoRingsV0_BDTG_Iso) {
656 +    std::vector<std::string> muonidiso_weightfiles;
657 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml"))));
658 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml"))));
659 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml"))));
660 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml"))));
661 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml"))));
662 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml"))));
663 +    fMuonTools = new MuonTools();
664 +    fMuonIDMVA = new MuonIDMVA();
665 +    fMuonIDMVA->Initialize("MuonIso_BDTG_IsoRings",
666 +                       MuonIDMVA::kIsoRingsV0,
667 +                       kTRUE,
668 +                       muonidiso_weightfiles);
669 +  }
670 +  else if(fMuIsoType == kIsoDeltaR) {
671 +    std::vector<std::string> muonidiso_weightfiles;
672 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LB_BDT.weights.xml"))));
673 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_LE_BDT.weights.xml"))));
674 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HB_BDT.weights.xml"))));
675 +    muonidiso_weightfiles.push_back(string((getenv("CMSSW_BASE")+string("/src/MitPhysics/data/MuonMVAWeights/MuonIsoMVA_santi-V1_HE_BDT.weights.xml"))));
676 +    fMuonTools = new MuonTools();
677 +    fMuonIDMVA = new MuonIDMVA();
678 +    fMuonIDMVA->Initialize("muonHZZ2012IsoDRMVA",
679 +                       MuonIDMVA::kIsoDeltaR,
680 +                       kTRUE,
681 +                       muonidiso_weightfiles);
682 +  }
683 +
684 + }
685 +
686 +
687 + //--------------------------------------------------------------------------------------------------
688 + Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
689 +                                         const PileupEnergyDensityCol *PileupEnergyDensity) const
690 + {
691 +
692 +  const Track *muTrk=0;
693 +  if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
694 +  else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
695 +  
696 +  Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
697 +
698 +  Int_t subdet = 0;
699 +  if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
700 +  else subdet = 1;
701 +  Int_t ptBin = 0;
702 +  if (muTrk->Pt() > 14.5) ptBin = 1;
703 +  if (muTrk->Pt() > 20.0) ptBin = 2;
704 +
705 +  Int_t MVABin = -1;
706 +  if      (subdet == 0 && ptBin == 0) MVABin = 0;
707 +  else if (subdet == 1 && ptBin == 0) MVABin = 1;
708 +  else if (subdet == 0 && ptBin == 1) MVABin = 2;
709 +  else if (subdet == 1 && ptBin == 1) MVABin = 3;
710 +  else if (subdet == 0 && ptBin == 2) MVABin = 4;
711 +  else if (subdet == 1 && ptBin == 2) MVABin = 5;
712 +
713 +  Double_t MVACut = -999;
714 +  if      (MVABin == 0) MVACut = -0.5618;
715 +  else if (MVABin == 1) MVACut = -0.3002;
716 +  else if (MVABin == 2) MVACut = -0.4642;
717 +  else if (MVABin == 3) MVACut = -0.2478;
718 +  else if (MVABin == 4) MVACut =  0.1706;
719 +  else if (MVABin == 5) MVACut =  0.8146;
720 +
721 +  if (MVAValue > MVACut) return kTRUE;
722 +  return kFALSE;
723 + }
724 +
725 + //--------------------------------------------------------------------------------------------------
726 + Bool_t MuonIDMod::PassMuonIsoRingsV0_BDTG_Iso(const Muon *mu, const Vertex *vertex,
727 +                                              const PileupEnergyDensityCol *PileupEnergyDensity) const
728 + {
729 +
730 +  const Track *muTrk=0;
731 +  if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
732 +  else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
733 +  
734 +  ElectronOArr *tempElectrons = new  ElectronOArr;
735 +  MuonOArr     *tempMuons     = new  MuonOArr;
736 +  Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,
737 +                      PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,kFALSE);
738 +  delete tempElectrons;
739 +  delete tempMuons;
740 +
741 +  Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
742 +
743 +  Double_t MVACut = -999;
744 +  if      (MVABin == 0) MVACut = -0.593;
745 +  else if (MVABin == 1) MVACut =  0.337;
746 +  else if (MVABin == 2) MVACut = -0.767;
747 +  else if (MVABin == 3) MVACut =  0.410;
748 +  else if (MVABin == 4) MVACut = -0.989;
749 +  else if (MVABin == 5) MVACut = -0.995;
750 +
751 +  if (MVAValue > MVACut) return kTRUE;
752 +  return kFALSE;
753 + }
754 +
755 + //--------------------------------------------------------------------------------------------------
756 + Bool_t MuonIDMod::PassMuonIsoDeltaR(const Muon *mu, const Vertex *vertex,
757 +                                    const PileupEnergyDensityCol *PileupEnergyDensity) const
758 + {
759 +
760 +  const Track *muTrk=0;
761 +  if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
762 +  else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
763 +  
764 +  ElectronOArr *tempElectrons = new  ElectronOArr;
765 +  MuonOArr     *tempMuons     = new  MuonOArr;
766 +  Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFNoPileUpCands,
767 +                      PileupEnergyDensity,MuonTools::kMuEAFall11MC,tempElectrons,tempMuons,kFALSE);
768 +  delete tempElectrons;
769 +  delete tempMuons;
770 +
771 +  Int_t MVABin = fMuonIDMVA->GetMVABin(muTrk->Eta(), muTrk->Pt(), mu->IsGlobalMuon(), mu->IsTrackerMuon());
772 +
773 +  Double_t MVACut = -999;
774 +  if      (MVABin == 0) MVACut =  0.000;
775 +  else if (MVABin == 1) MVACut =  0.000;
776 +  else if (MVABin == 2) MVACut =  0.000;
777 +  else if (MVABin == 3) MVACut =  0.000;
778 +
779 +  if (MVAValue > MVACut) return kTRUE;
780 +  return kFALSE;
781 + }
782 +
783 + //--------------------------------------------------------------------------------------------------
784 + void MuonIDMod::Terminate()
785 + {
786 +  // Run finishing code on the computer (slave) that did the analysis
787 +  delete fMuonIDMVA;
788 +  
789 +  delete fMuonTools;
790   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines