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.35 by ceballos, Wed Oct 20 02:44:16 2010 UTC vs.
Revision 1.73 by ceballos, Thu May 3 10:22:10 2012 UTC

# Line 1 | Line 1
1 <        // $Id$
1 > // $Id$
2  
3   #include "MitPhysics/Mods/interface/MuonIDMod.h"
4   #include "MitCommon/MathTools/interface/MathUtils.h"
5 < #include "MitAna/DataTree/interface/MuonCol.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  
# Line 13 | 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"),  
21 +  fNonIsolatedElectronsName("random"),  
22    fVertexName(ModNames::gkGoodVertexesName),
23 <  fMuonIDType("Minimal"),
24 <  fMuonIsoType("TrackCaloSliding"),  
23 >  fBeamSpotName(Names::gkBeamSpotBrn),
24 >  fTrackName(Names::gkTrackBrn),
25 >  fPFCandidatesName(Names::gkPFCandidatesBrn),
26 >  fMuonIDType("WWMuIdV3"),
27 >  fMuonIsoType("PFIso"),
28    fMuonClassType("Global"),  
29    fTrackIsolationCut(3.0),
30    fCaloIsolationCut(3.0),
31 <  fCombIsolationCut(5.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.10),
39 +  fWhichVertex(-1),
40    fEtaCut(2.4),
29  fReverseIsoCut(kFALSE),
30  fReverseD0Cut(kFALSE),
41    fMuIDType(kIdUndef),
42    fMuIsoType(kIsoUndef),
43    fMuClassType(kClassUndef),
44    fMuons(0),
45    fVertices(0),
46 <  fMuonTools(0)
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),
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 43 | Line 70 | void MuonIDMod::Process()
70   {
71    // Process entries of the tree.
72  
73 <  LoadEventObject(fMuonBranchName, fMuons);
73 >  if(fMuIsoType != kPFIsoNoL) {
74 >    LoadEventObject(fMuonBranchName, fMuons);
75 >  }
76 >  else {
77 >    fMuons = GetObjThisEvt<MuonOArr>(fMuonBranchName);
78 >  }
79 >  LoadEventObject(fBeamSpotName, fBeamSpot);
80 >  LoadEventObject(fTrackName, fTracks);
81 >  LoadEventObject(fPFCandidatesName, fPFCandidates);
82 >  if(fMuIsoType == kTrackCaloSliding ||
83 >     fMuIsoType == kCombinedRelativeConeAreaCorrected ||
84 >     fMuIsoType == kPFIsoEffectiveAreaCorrected ||
85 >     fMuIsoType == kMVAIso_BDTG_IDIso
86 >    ) {
87 >    LoadEventObject(fPileupEnergyDensityName, fPileupEnergyDensity);
88 >  }
89 >  if(fMuIsoType == kPFRadialIso){
90 >    // Name is hardcoded, can be changed if someone feels to do it
91 >    fPFNoPileUpCands = GetObjThisEvt<PFCandidateCol>("PFNoPileUp");    
92 >  }
93  
94    MuonOArr *CleanMuons = new MuonOArr;
95    CleanMuons->SetName(fCleanMuonsName);
96  
97 <  for (UInt_t i=0; i<fMuons->GetEntries(); ++i) {
97 >  fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
98 >
99 >  for (UInt_t i=0; i<fMuons->GetEntries() && fVertices->GetEntries() > 0 ; ++i) {
100      const Muon *mu = fMuons->At(i);
101  
102      Bool_t pass = kFALSE;
# Line 73 | Line 121 | void MuonIDMod::Process()
121            eta = TMath::Abs(mu->Eta());
122          }
123          break;
124 +      case kGlobalTracker:
125 +        pass = (mu->HasGlobalTrk() && mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof() < 10 &&
126 +               (mu->NSegments() > 1 || mu->NMatches() > 1) && mu->NValidHits() > 0) ||
127 +               (mu->IsTrackerMuon() &&
128 +                mu->Quality().Quality(MuonQuality::TMLastStationTight));
129 +        if (pass) {
130 +          pt  = mu->TrackerTrk()->Pt();
131 +          eta = TMath::Abs(mu->TrackerTrk()->Eta());
132 +        }
133 +        else {
134 +          pt  = mu->Pt();
135 +          eta = TMath::Abs(mu->Eta());
136 +        }
137 +        break;
138        case kSta:
139          pass = mu->HasStandaloneTrk();
140          if (pass) {
# Line 115 | Line 177 | void MuonIDMod::Process()
177      if (eta >= fEtaCut)
178        continue;
179  
180 +
181 +    //***********************************************************************************************
182 +    //Debug Info For Lepton MVA
183 +    //***********************************************************************************************
184 +    if( fPrintMVADebugInfo &&
185 +        (fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso)
186 +      ) {
187 +      cout << "Event: " << GetEventHeader()->RunNum() << " " << GetEventHeader()->LumiSec() << " "
188 +           << GetEventHeader()->EvtNum() << " : Rho = " << fPileupEnergyDensity->At(0)->Rho()
189 +           << " : Muon " << i << " "
190 +           << endl;
191 +      fMuonIDMVA->MVAValue(mu,fVertices->At(0),fMuonTools,fPFCandidates,fPileupEnergyDensity,kTRUE);
192 +    }
193 +    //***********************************************************************************************
194 +
195 +
196      Double_t RChi2 = 0.0;
197      if     (mu->HasGlobalTrk()) {
198        RChi2 = mu->GlobalTrk()->Chi2()/mu->GlobalTrk()->Ndof();
# Line 128 | Line 206 | void MuonIDMod::Process()
206          idpass = mu->BestTrk() != 0 &&
207                   mu->BestTrk()->NHits() > 10 &&
208                   RChi2 < 10.0 &&
209 <                 mu->NSegments() > 1 &&
209 >                (mu->NSegments() > 1 || mu->NMatches() > 1) &&
210                   mu->BestTrk()->NPixelHits() > 0 &&
211                   mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
212          break;
213        case kZMuId:
214          idpass = mu->BestTrk() != 0 &&
215                   mu->BestTrk()->NHits() > 10 &&
216 <                 mu->NSegments() > 1 &&
216 >                (mu->NSegments() > 1 || mu->NMatches() > 1) &&
217                   mu->BestTrk()->NPixelHits() > 0 &&
218                   mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
219          break;
# Line 155 | Line 233 | void MuonIDMod::Process()
233                   RChi2 < 10.0 &&
234                   mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
235          break;
236 <      case kMinimal:
236 >      case kWWMuIdV1:
237          idpass = mu->BestTrk() != 0 &&
238                   mu->BestTrk()->NHits() > 10 &&
239 +                 mu->BestTrk()->NPixelHits() > 0 &&
240 +                 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
241                   RChi2 < 10.0 &&
242 +                (mu->NSegments() > 1 || mu->NMatches() > 1) &&
243                   mu->Quality().Quality(MuonQuality::GlobalMuonPromptTight);
244          break;
245 +      case kWWMuIdV2:
246 +        idpass = mu->BestTrk() != 0 &&
247 +                 mu->BestTrk()->NHits() > 10 &&
248 +                 mu->BestTrk()->NPixelHits() > 0 &&
249 +                 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1;
250 +        break;
251 +      case kWWMuIdV3:
252 +        idpass = mu->BestTrk() != 0 &&
253 +                 mu->BestTrk()->NHits() > 10 &&
254 +                 mu->BestTrk()->NPixelHits() > 0 &&
255 +                 mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
256 +                 mu->TrkKink() < 20.0;
257 +        break;
258 +      case kMVAID_BDTG_IDIso:
259 +        {
260 +          Bool_t passDenominatorM2 = (mu->BestTrk() != 0 &&
261 +                                      mu->BestTrk()->NHits() > 10 &&
262 +                                      mu->BestTrk()->NPixelHits() > 0 &&
263 +                                      mu->BestTrk()->PtErr()/mu->BestTrk()->Pt() < 0.1 &&
264 +                                      MuonTools::PassD0Cut(mu, fVertices, 0.20, 0) &&
265 +                                      MuonTools::PassDZCut(mu, fVertices, 0.10, 0) &&
266 +                                      mu->TrkKink() < 20.0
267 +            );  
268 +          idpass =  passDenominatorM2;
269 +          //only evaluate MVA if muon passes M2 denominator to save time
270 +          if (idpass) idpass = PassMuonMVA_BDTG_IdIso(mu, fVertices->At(0), fPileupEnergyDensity);
271 +        }
272 +        break;
273        case kNoId:
274          idpass = kTRUE;
275          break;
# Line 171 | Line 280 | void MuonIDMod::Process()
280      if (!idpass)
281        continue;
282  
283 +    Double_t Rho = 0.;
284 +    if(fPileupEnergyDensity){
285 +      const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0);
286 +
287 +      switch (fTheRhoType) {
288 +      case RhoUtilities::MIT_RHO_VORONOI_LOW_ETA:
289 +        Rho = rho->RhoLowEta();
290 +        break;
291 +      case RhoUtilities::MIT_RHO_VORONOI_HIGH_ETA:
292 +        Rho = rho->Rho();
293 +        break;
294 +      case RhoUtilities::MIT_RHO_RANDOM_LOW_ETA:
295 +        Rho = rho->RhoRandomLowEta();
296 +        break;
297 +      case RhoUtilities::MIT_RHO_RANDOM_HIGH_ETA:
298 +        Rho = rho->RhoRandom();
299 +        break;
300 +      default:
301 +        Rho = rho->Rho();
302 +      }
303 +
304 +      if ((TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho()))) Rho = 0.;
305 +    }
306 +
307      Bool_t isocut = kFALSE;
308      switch (fMuIsoType) {
309        case kTrackCalo:
# Line 178 | Line 311 | void MuonIDMod::Process()
311            (mu->IsoR03EmEt() + mu->IsoR03HadEt() < fCaloIsolationCut);
312          break;
313        case kTrackCaloCombined:
314 <        isocut = (1.0 * mu->IsoR03SumPt() + 1.0 * mu->IsoR03EmEt() +
315 <                  1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
314 >        isocut = (1.0 * mu->IsoR03SumPt() +
315 >                  1.0 * mu->IsoR03EmEt()  +
316 >                  1.0 * mu->IsoR03HadEt() < fCombIsolationCut);
317          break;
318        case kTrackCaloSliding:
319          {
320 <          Double_t totalIso = 1.0 * mu->IsoR03SumPt() +
321 <                              1.0 * mu->IsoR03EmEt() +
322 <                              1.0 * mu->IsoR03HadEt();
323 <          if (totalIso < (mu->Pt()*0.15) )
320 >          Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
321 >          // trick to change the signal region cut
322 >          double theIsoCut = fCombIsolationCut;
323 >          if(theIsoCut < 0.20){
324 >            if(mu->Pt() >  20.0) theIsoCut = 0.15;
325 >            else                 theIsoCut = 0.10;
326 >          }
327 >          if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
328 >        }
329 >        break;
330 >      case kTrackCaloSlidingNoCorrection:
331 >        {
332 >          Double_t totalIso =  1.0 * mu->IsoR03SumPt() +
333 >                               1.0 * mu->IsoR03EmEt()  +
334 >                               1.0 * mu->IsoR03HadEt();
335 >          // trick to change the signal region cut
336 >          double theIsoCut = fCombIsolationCut;
337 >          if(theIsoCut < 0.20){
338 >            if(mu->Pt() >  20.0) theIsoCut = 0.15;
339 >            else                 theIsoCut = 0.10;
340 >          }
341 >          if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
342 >        }
343 >        break;
344 >      case kCombinedRelativeConeAreaCorrected:
345 >        {
346 >          //const PileupEnergyDensity *rho =  fPileupEnergyDensity->At(0); // Fabian: made Rho customable
347 >          Double_t totalIso =  mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt() - Rho * TMath::Pi() * 0.3 * 0.3 ;
348 >          double theIsoCut = fCombRelativeIsolationCut;
349 >          if (totalIso < (mu->Pt()*theIsoCut)) isocut = kTRUE;
350 >        }
351 >        break;          
352 >    case kCombinedRelativeEffectiveAreaCorrected:
353 >      {
354 >        Double_t tmpRho = Rho;   // Fabian: made the Rho type customable.
355 >        //if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || std::isinf(fPileupEnergyDensity->At(0)->Rho())))
356 >        //tmpRho = fPileupEnergyDensity->At(0)->Rho();
357 >        
358 >          isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
359 >                     -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
360 >                     -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
361 >            ) < (mu->Pt()* 0.40);
362 >        }
363 >        break;          
364 >      case kPFIso:
365 >        {
366 >          Double_t pfIsoCutValue = 9999;
367 >          if(fPFIsolationCut > 0){
368 >            pfIsoCutValue = fPFIsolationCut;
369 >          } else {
370 >            if (mu->AbsEta() < 1.479) {
371 >              if (mu->Pt() > 20) {
372 >                pfIsoCutValue = 0.13;
373 >              } else {
374 >                pfIsoCutValue = 0.06;
375 >              }
376 >            } else {
377 >              if (mu->Pt() > 20) {
378 >                pfIsoCutValue = 0.09;
379 >              } else {
380 >                pfIsoCutValue = 0.05;
381 >              }
382 >            }
383 >          }
384 >          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
385 >          if (totalIso < (mu->Pt()*pfIsoCutValue) )
386              isocut = kTRUE;
191
192          if     (fReverseIsoCut == kTRUE &&
193                  isocut == kFALSE && totalIso < 10)
194            isocut = kTRUE;
195          else if(fReverseIsoCut == kTRUE)
196            isocut = kFALSE;
387          }
388          break;
389 +      case kPFRadialIso:
390 +        {
391 +          Double_t pfIsoCutValue = 9999;
392 +          if(fPFIsolationCut > 0){
393 +            pfIsoCutValue = fPFIsolationCut;
394 +          } else {
395 +            if (mu->Pt() > 20) {
396 +              pfIsoCutValue = 0.10;
397 +            } else {
398 +              pfIsoCutValue = 0.05;
399 +            }
400 +          }
401 +          Double_t totalIso =  IsolationTools::PFRadialMuonIsolation(mu, fPFNoPileUpCands, 1.0, 0.3);
402 +          if (totalIso < (mu->Pt()*pfIsoCutValue) )
403 +            isocut = kTRUE;
404 +        }
405 +        break;
406 +      case kPFIsoEffectiveAreaCorrected:
407 +        {
408 +          Double_t pfIsoCutValue = 9999;
409 +          if(fPFIsolationCut > 0){
410 +            pfIsoCutValue = fPFIsolationCut;
411 +          } else {
412 +            pfIsoCutValue = fPFIsolationCut; //leave it like this for now
413 +          }
414 +          Double_t EffectiveAreaCorrectedPFIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius)
415 +            - Rho * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());
416 +          //- fPileupEnergyDensity->At(0)->Rho() * MuonTools::MuonEffectiveArea(MuonTools::kMuNeutralIso03, mu->Eta());  // Fabian: made Rho-type customable
417 +          isocut = EffectiveAreaCorrectedPFIso < (mu->Pt() * pfIsoCutValue);
418 +          break;
419 +        }
420 +      case kPFIsoNoL:
421 +        {
422 +          fNonIsolatedMuons     = GetObjThisEvt<MuonCol>(fNonIsolatedMuonsName);
423 +          fNonIsolatedElectrons = GetObjThisEvt<ElectronCol>(fNonIsolatedElectronsName);
424 +
425 +          Double_t pfIsoCutValue = 9999;
426 +          if(fPFIsolationCut > 0){
427 +            pfIsoCutValue = fPFIsolationCut;
428 +          } else {
429 +            if (mu->AbsEta() < 1.479) {
430 +              if (mu->Pt() > 20) {
431 +                pfIsoCutValue = 0.13;
432 +              } else {
433 +                pfIsoCutValue = 0.06;
434 +              }
435 +            } else {
436 +              if (mu->Pt() > 20) {
437 +                pfIsoCutValue = 0.09;
438 +              } else {
439 +                pfIsoCutValue = 0.05;
440 +              }
441 +            }
442 +          }
443 +          Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fNonIsolatedMuons, fNonIsolatedElectrons, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
444 +          if (totalIso < (mu->Pt()*pfIsoCutValue) )
445 +            isocut = kTRUE;
446 +        }
447 +        break;
448 +      case kMVAIso_BDTG_IDIso:
449 +      {
450 +
451 +        // **************************************************************************
452 +        // Don't use effective area correction denominator. Instead use the old one.
453 +        // **************************************************************************
454 +
455 +        //         Double_t tmpRho = 0;
456 +        //         if (!(TMath::IsNaN(fPileupEnergyDensity->At(0)->Rho()) || isinf(fPileupEnergyDensity->At(0)->Rho())))
457 +        //           tmpRho = fPileupEnergyDensity->At(0)->Rho();
458 +        
459 +        //         isocut = ( mu->IsoR03SumPt() + mu->IsoR03EmEt() + mu->IsoR03HadEt()
460 +        //                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuEMIso03, mu->Eta())
461 +        //                    -  tmpRho*MuonTools::MuonEffectiveArea(MuonTools::kMuHadIso03, mu->Eta())
462 +        //           ) < (mu->Pt()* 0.40);
463 +        
464 +        Double_t totalIso =  IsolationTools::PFMuonIsolation(mu, fPFCandidates, fVertices->At(0), 0.1, 1.0, 0.3, 0.0, fIntRadius);
465 +        isocut = (totalIso < (mu->Pt()*0.4));
466 +
467 +      }
468 +        break;
469        case kNoIso:
470          isocut = kTRUE;
471          break;
# Line 207 | Line 477 | void MuonIDMod::Process()
477      if (isocut == kFALSE)
478        continue;
479  
480 +    // apply d0 cut
481      if (fApplyD0Cut) {
482 <      fVertices = GetObjThisEvt<VertexOArr>(fVertexName);
483 <      Bool_t passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fReverseD0Cut);
482 >      Bool_t passD0cut = kTRUE;
483 >      if(fD0Cut < 0.05) { // trick to change the signal region cut
484 >        if      (mu->Pt() >  20.0) fD0Cut = 0.02;
485 >        else if (mu->Pt() <= 20.0) fD0Cut = 0.01;
486 >      }
487 >      if(fWhichVertex >= -1) passD0cut = MuonTools::PassD0Cut(mu, fVertices, fD0Cut, fWhichVertex);
488 >      else                   passD0cut = MuonTools::PassD0Cut(mu, fBeamSpot, fD0Cut);
489        if (!passD0cut)
490          continue;
491      }
492  
493 +    // apply dz cut
494 +    if (fApplyDZCut) {
495 +      Bool_t passDZcut = MuonTools::PassDZCut(mu, fVertices, fDZCut, fWhichVertex);
496 +      if (!passDZcut)
497 +        continue;
498 +    }
499 +
500      // add good muon
501      CleanMuons->Add(mu);
502    }
# Line 231 | Line 514 | void MuonIDMod::SlaveBegin()
514    // Run startup code on the computer (slave) doing the actual analysis. Here,
515    // we just request the muon collection branch.
516  
517 <  ReqEventObject(fMuonBranchName, fMuons, kTRUE);
517 >   // In this case we cannot have a branch
518 >  if (fMuonIsoType.CompareTo("PFIsoNoL") != 0) {
519 >    ReqEventObject(fMuonBranchName, fMuons, kTRUE);
520 >  }
521 >  ReqEventObject(fBeamSpotName, fBeamSpot, kTRUE);
522 >  ReqEventObject(fTrackName, fTracks, kTRUE);
523 >  ReqEventObject(fPFCandidatesName, fPFCandidates, kTRUE);
524 >  if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0
525 >      || fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0
526 >      || fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0
527 >       || fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0
528 >      || fMuonIsoType.CompareTo("MVA_BDTG_IDIso") == 0
529 >    ) {
530 >    ReqEventObject(fPileupEnergyDensityName, fPileupEnergyDensity, kTRUE);
531 >  }
532  
236  fMuonTools = new MuonTools;
533  
534    if (fMuonIDType.CompareTo("WMuId") == 0)
535      fMuIDType = kWMuId;
# Line 243 | Line 539 | void MuonIDMod::SlaveBegin()
539      fMuIDType = kTight;
540    else if (fMuonIDType.CompareTo("Loose") == 0)
541      fMuIDType = kLoose;
542 <  else if (fMuonIDType.CompareTo("Minimal") == 0)
543 <    fMuIDType = kMinimal;
542 >  else if (fMuonIDType.CompareTo("WWMuIdV1") == 0)
543 >    fMuIDType = kWWMuIdV1;
544 >  else if (fMuonIDType.CompareTo("WWMuIdV2") == 0)
545 >    fMuIDType = kWWMuIdV2;
546 >  else if (fMuonIDType.CompareTo("WWMuIdV3") == 0)
547 >    fMuIDType = kWWMuIdV3;
548    else if (fMuonIDType.CompareTo("NoId") == 0)
549      fMuIDType = kNoId;
550    else if (fMuonIDType.CompareTo("Custom") == 0) {
551      fMuIDType = kCustomId;
552      SendError(kWarning, "SlaveBegin",
553                "Custom muon identification is not yet implemented.");
554 +  } else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
555 +    fMuIDType = kMVAID_BDTG_IDIso;
556    } else {
557      SendError(kAbortAnalysis, "SlaveBegin",
558                "The specified muon identification %s is not defined.",
# Line 264 | Line 566 | void MuonIDMod::SlaveBegin()
566      fMuIsoType = kTrackCaloCombined;
567    else if (fMuonIsoType.CompareTo("TrackCaloSliding") == 0)
568      fMuIsoType = kTrackCaloSliding;
569 +  else if (fMuonIsoType.CompareTo("TrackCaloSlidingNoCorrection") == 0)
570 +    fMuIsoType = kTrackCaloSlidingNoCorrection;
571 +  else if (fMuonIsoType.CompareTo("CombinedRelativeConeAreaCorrected") == 0)
572 +    fMuIsoType = kCombinedRelativeConeAreaCorrected;
573 +  else if (fMuonIsoType.CompareTo("CombinedRelativeEffectiveAreaCorrected") == 0)
574 +    fMuIsoType = kCombinedRelativeEffectiveAreaCorrected;
575 +  else if (fMuonIsoType.CompareTo("PFIso") == 0)
576 +    fMuIsoType = kPFIso;
577 +  else if (fMuonIsoType.CompareTo("PFRadialIso") == 0)
578 +    fMuIsoType = kPFRadialIso;
579 +  else if (fMuonIsoType.CompareTo("PFIsoEffectiveAreaCorrected") == 0)
580 +    fMuIsoType = kPFIsoEffectiveAreaCorrected;
581 +  else if (fMuonIsoType.CompareTo("PFIsoNoL") == 0)
582 +    fMuIsoType = kPFIsoNoL;
583    else if (fMuonIsoType.CompareTo("NoIso") == 0)
584      fMuIsoType = kNoIso;
585    else if (fMuonIsoType.CompareTo("Custom") == 0) {
586      fMuIsoType = kCustomIso;
587      SendError(kWarning, "SlaveBegin",
588                "Custom muon isolation is not yet implemented.");
589 +  } else if (fMuonIDType.CompareTo("MVA_BDTG_IDIso") == 0) {
590 +    fMuIsoType = kMVAIso_BDTG_IDIso;
591    } else {
592      SendError(kAbortAnalysis, "SlaveBegin",
593                "The specified muon isolation %s is not defined.",
# Line 281 | Line 599 | void MuonIDMod::SlaveBegin()
599      fMuClassType = kAll;
600    else if (fMuonClassType.CompareTo("Global") == 0)
601      fMuClassType = kGlobal;
602 +  else if (fMuonClassType.CompareTo("GlobalTracker") == 0)
603 +    fMuClassType = kGlobalTracker;
604    else if (fMuonClassType.CompareTo("Standalone") == 0)
605      fMuClassType = kSta;
606    else if (fMuonClassType.CompareTo("TrackerMuon") == 0)
# Line 295 | Line 615 | void MuonIDMod::SlaveBegin()
615                fMuonClassType.Data());
616      return;
617    }
618 +
619 +
620 +  //If we use MVA ID, need to load MVA weights
621 +  if(fMuIsoType == kMVAIso_BDTG_IDIso || fMuIDType == kMVAID_BDTG_IDIso) {
622 +    fMuonTools = new MuonTools();
623 +    fMuonIDMVA = new MuonIDMVA();
624 +    fMuonIDMVA->Initialize("BDTG method",
625 +                           fMuonMVAWeights_Subdet0Pt10To14p5,
626 +                           fMuonMVAWeights_Subdet1Pt10To14p5,
627 +                           fMuonMVAWeights_Subdet0Pt14p5To20,
628 +                           fMuonMVAWeights_Subdet1Pt14p5To20,
629 +                           fMuonMVAWeights_Subdet0Pt20ToInf,
630 +                           fMuonMVAWeights_Subdet1Pt20ToInf,
631 +                           MuonIDMVA::kIDIsoCombinedDetIso);
632 +  }
633 +
634 + }
635 +
636 +
637 + //--------------------------------------------------------------------------------------------------
638 + Bool_t MuonIDMod::PassMuonMVA_BDTG_IdIso(const Muon *mu, const Vertex *vertex,
639 +                                         const PileupEnergyDensityCol *PileupEnergyDensity) const
640 + {
641 +
642 +  const Track *muTrk=0;
643 +  if(mu->HasTrackerTrk())         { muTrk = mu->TrackerTrk();    }
644 +  else if(mu->HasStandaloneTrk()) { muTrk = mu->StandaloneTrk(); }
645 +  
646 +  Double_t MVAValue = fMuonIDMVA->MVAValue(mu,vertex,fMuonTools,fPFCandidates,PileupEnergyDensity);
647 +
648 +  Int_t subdet = 0;
649 +  if (fabs(muTrk->Eta()) < 1.479) subdet = 0;
650 +  else subdet = 1;
651 +  Int_t ptBin = 0;
652 +  if (muTrk->Pt() > 14.5) ptBin = 1;
653 +  if (muTrk->Pt() > 20.0) ptBin = 2;
654 +
655 +  Int_t MVABin = -1;
656 +  if      (subdet == 0 && ptBin == 0) MVABin = 0;
657 +  else if (subdet == 1 && ptBin == 0) MVABin = 1;
658 +  else if (subdet == 0 && ptBin == 1) MVABin = 2;
659 +  else if (subdet == 1 && ptBin == 1) MVABin = 3;
660 +  else if (subdet == 0 && ptBin == 2) MVABin = 4;
661 +  else if (subdet == 1 && ptBin == 2) MVABin = 5;
662 +
663 +  Double_t MVACut = -999;
664 +  if      (MVABin == 0) MVACut = -0.5618;
665 +  else if (MVABin == 1) MVACut = -0.3002;
666 +  else if (MVABin == 2) MVACut = -0.4642;
667 +  else if (MVABin == 3) MVACut = -0.2478;
668 +  else if (MVABin == 4) MVACut =  0.1706;
669 +  else if (MVABin == 5) MVACut =  0.8146;
670 +
671 +  if (MVAValue > MVACut) return kTRUE;
672 +  return kFALSE;
673   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines