Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:23:58

0001 #include "PhysicsTools/PatAlgos/interface/SoftMuonMvaEstimator.h"
0002 
0003 #include "FWCore/ParameterSet/interface/FileInPath.h"
0004 #include "CommonTools/MVAUtils/interface/GBRForestTools.h"
0005 #include "DataFormats/Candidate/interface/Candidate.h"
0006 #include "DataFormats/MuonReco/interface/Muon.h"
0007 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0008 #include "DataFormats/PatCandidates/interface/Muon.h"
0009 
0010 #include "TMath.h"
0011 
0012 using namespace pat;
0013 
0014 SoftMuonMvaEstimator::SoftMuonMvaEstimator(const edm::FileInPath& weightsfile) {
0015   gbrForest_ = createGBRForest(weightsfile);
0016 }
0017 
0018 SoftMuonMvaEstimator::~SoftMuonMvaEstimator() {}
0019 
0020 namespace {
0021   enum inputIndexes {
0022     kSegmentCompatibility,
0023     kChi2LocalMomentum,
0024     kChi2LocalPosition,
0025     kGlbTrackProbability,
0026     kIValidFraction,
0027     kLayersWithMeasurement,
0028     kTrkKink,
0029     kLog2PlusGlbKink,
0030     kTimeAtIpInOutErr,
0031     kOuterChi2,
0032     kInnerChi2,
0033     kTrkRelChi2,
0034     kVMuonHitComb,
0035     kQProd,
0036     kPID,
0037     kPt,
0038     kEta,
0039     kMomID,
0040     kLast
0041   };
0042 }
0043 
0044 float SoftMuonMvaEstimator::computeMva(const pat::Muon& muon) const {
0045   float var[kLast]{};
0046 
0047   reco::TrackRef gTrack = muon.globalTrack();
0048   reco::TrackRef iTrack = muon.innerTrack();
0049   reco::TrackRef oTrack = muon.outerTrack();
0050 
0051   if (!(muon.innerTrack().isNonnull() and muon.outerTrack().isNonnull() and muon.globalTrack().isNonnull())) {
0052     return -1;
0053   }
0054 
0055   //VARIABLE EXTRACTION
0056   var[kPt] = muon.pt();
0057   var[kEta] = muon.eta();
0058   var[kMomID] = -1;
0059   var[kPID] = -1;
0060 
0061   var[kChi2LocalMomentum] = muon.combinedQuality().chi2LocalMomentum;
0062   var[kChi2LocalPosition] = muon.combinedQuality().chi2LocalPosition;
0063   var[kGlbTrackProbability] = muon.combinedQuality().glbTrackProbability;
0064   var[kTrkRelChi2] = muon.combinedQuality().trkRelChi2;
0065 
0066   var[kTrkKink] = muon.combinedQuality().trkKink;
0067   var[kLog2PlusGlbKink] = TMath::Log(2 + muon.combinedQuality().glbKink);
0068   var[kSegmentCompatibility] = muon.segmentCompatibility();
0069 
0070   var[kTimeAtIpInOutErr] = muon.time().timeAtIpInOutErr;
0071 
0072   //TRACK RELATED VARIABLES
0073 
0074   var[kIValidFraction] = iTrack->validFraction();
0075   var[kInnerChi2] = iTrack->normalizedChi2();
0076   var[kLayersWithMeasurement] = iTrack->hitPattern().trackerLayersWithMeasurement();
0077 
0078   var[kOuterChi2] = oTrack->normalizedChi2();
0079 
0080   var[kQProd] = iTrack->charge() * oTrack->charge();
0081 
0082   //vComb Calculation
0083 
0084   const reco::HitPattern& gMpattern = gTrack->hitPattern();
0085 
0086   std::vector<int> fvDThits{0, 0, 0, 0};
0087   std::vector<int> fvRPChits{0, 0, 0, 0};
0088   std::vector<int> fvCSChits{0, 0, 0, 0};
0089 
0090   var[kVMuonHitComb] = 0;
0091 
0092   for (int i = 0; i < gMpattern.numberOfAllHits(reco::HitPattern::TRACK_HITS); i++) {
0093     uint32_t hit = gMpattern.getHitPattern(reco::HitPattern::TRACK_HITS, i);
0094     if (!gMpattern.validHitFilter(hit))
0095       continue;
0096 
0097     int muStation0 = gMpattern.getMuonStation(hit) - 1;
0098     if (muStation0 >= 0 && muStation0 < 4) {
0099       if (gMpattern.muonDTHitFilter(hit))
0100         fvDThits[muStation0]++;
0101       if (gMpattern.muonRPCHitFilter(hit))
0102         fvRPChits[muStation0]++;
0103       if (gMpattern.muonCSCHitFilter(hit))
0104         fvCSChits[muStation0]++;
0105     }
0106   }
0107 
0108   for (unsigned int station = 0; station < 4; ++station) {
0109     var[kVMuonHitComb] += (fvDThits[station]) / 2.;
0110     var[kVMuonHitComb] += fvRPChits[station];
0111 
0112     if (fvCSChits[station] > 6) {
0113       var[kVMuonHitComb] += 6;
0114     } else {
0115       var[kVMuonHitComb] += fvCSChits[station];
0116     }
0117   }
0118 
0119   if (var[kChi2LocalMomentum] < 5000 and var[kChi2LocalPosition] < 2000 and var[kGlbTrackProbability] < 5000 and
0120       var[kTrkKink] < 900 and var[kLog2PlusGlbKink] < 50 and var[kTimeAtIpInOutErr] < 4 and var[kOuterChi2] < 1000 and
0121       var[kInnerChi2] < 10 and var[kTrkRelChi2] < 3) {
0122     return gbrForest_->GetAdaBoostClassifier(var);
0123   } else {
0124     return -1;
0125   }
0126 }