SoftMuonMvaEstimator.cc

Go to the documentation of this file.

#include "PhysicsTools/PatAlgos/interface/SoftMuonMvaEstimator.h"

#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "CommonTools/MVAUtils/interface/GBRForestTools.h"
#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonSelectors.h"
#include "DataFormats/PatCandidates/interface/Muon.h"

using namespace pat;

SoftMuonMvaEstimator::SoftMuonMvaEstimator(const edm::FileInPath& weightsfile)
{
  gbrForest_ = createGBRForest( weightsfile );
}

SoftMuonMvaEstimator::~SoftMuonMvaEstimator() { }

namespace {
  enum inputIndexes {
    kSegmentCompatibility,
    kChi2LocalMomentum,
    kChi2LocalPosition,
    kGlbTrackProbability,
    kIValidFraction,
    kLayersWithMeasurement,
    kTrkKink,
    kLog2PlusGlbKink,
    kTimeAtIpInOutErr,
    kOuterChi2,
    kInnerChi2,
    kTrkRelChi2,
    kVMuonHitComb,
    kQProd,
    kPID,
    kPt,
    kEta,
    kMomID,
    kLast
  };
}

float SoftMuonMvaEstimator::computeMva(const pat::Muon& muon) const
{
    float var[kLast]{};

    reco::TrackRef gTrack = muon.globalTrack();
    reco::TrackRef iTrack = muon.innerTrack();
    reco::TrackRef oTrack = muon.outerTrack();

    if(!(muon.innerTrack().isNonnull() and
         muon.outerTrack().isNonnull() and
         muon.globalTrack().isNonnull()))
    {
      return -1;
    }

    //VARIABLE EXTRACTION
    var[kPt] = muon.pt();
    var[kEta] = muon.eta();
    var[kMomID] = -1;
    var[kPID] = -1;

    var[kChi2LocalMomentum] = muon.combinedQuality().chi2LocalMomentum;
    var[kChi2LocalPosition] = muon.combinedQuality().chi2LocalPosition;
    var[kGlbTrackProbability] =  muon.combinedQuality().glbTrackProbability;
    var[kTrkRelChi2] =  muon.combinedQuality().trkRelChi2;

    var[kTrkKink] = muon.combinedQuality().trkKink;
    var[kLog2PlusGlbKink] =  TMath::Log(2+muon.combinedQuality().glbKink);
    var[kSegmentCompatibility] = muon.segmentCompatibility();

    var[kTimeAtIpInOutErr] = muon.time().timeAtIpInOutErr;

    //TRACK RELATED VARIABLES
    
    var[kIValidFraction] =  iTrack->validFraction();
    var[kInnerChi2] =  iTrack->normalizedChi2();
    var[kLayersWithMeasurement] = iTrack->hitPattern().trackerLayersWithMeasurement();

    var[kOuterChi2] = oTrack->normalizedChi2();

    var[kQProd] =  iTrack->charge()*oTrack->charge();

    //vComb Calculation

    const reco::HitPattern &gMpattern = gTrack->hitPattern();

    std::vector<int> fvDThits {0,0,0,0};
    std::vector<int> fvRPChits {0,0,0,0};
    std::vector<int> fvCSChits {0,0,0,0};

    var[kVMuonHitComb] = 0;

    for (int i=0;i<gMpattern.numberOfAllHits(reco::HitPattern::TRACK_HITS);i++){

      uint32_t hit = gMpattern.getHitPattern(reco::HitPattern::TRACK_HITS, i);
      if ( !gMpattern.validHitFilter(hit) ) continue;

      int muStation0 = gMpattern.getMuonStation(hit) - 1;
      if (muStation0 >=0 && muStation0 < 4){
        if (gMpattern.muonDTHitFilter(hit)) fvDThits[muStation0]++;
        if (gMpattern.muonRPCHitFilter(hit)) fvRPChits[muStation0]++;
        if (gMpattern.muonCSCHitFilter(hit)) fvCSChits[muStation0]++;
      }

    }
  

    for (unsigned int station = 0; station < 4; ++station) {

      var[kVMuonHitComb] += (fvDThits[station])/2.;
      var[kVMuonHitComb] += fvRPChits[station];

      if (fvCSChits[station] > 6){
        var[kVMuonHitComb] += 6;
      }else{
        var[kVMuonHitComb] += fvCSChits[station];
      }

    }

    if(var[kChi2LocalMomentum] < 5000 and var[kChi2LocalPosition] < 2000 and
       var[kGlbTrackProbability] < 5000 and var[kTrkKink] < 900 and
       var[kLog2PlusGlbKink] < 50 and var[kTimeAtIpInOutErr] < 4 and
       var[kOuterChi2] < 1000 and var[kInnerChi2] < 10 and var[kTrkRelChi2] < 3)
    {
      return gbrForest_->GetAdaBoostClassifier(var);
    } else {
      return -1;
        }
}