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
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
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
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 }