PatAlgos/src/SoftMuonMvaRun3Estimator.cc

0001 #include "PhysicsTools/PatAlgos/interface/SoftMuonMvaRun3Estimator.h"
0002 #include "DataFormats/PatCandidates/interface/Muon.h"
0003 #include "PhysicsTools/XGBoost/interface/XGBooster.h"
0004
0005 typedef std::pair<const reco::MuonChamberMatch*, const reco::MuonSegmentMatch*> MatchPair;
0006
0007 const MatchPair& getBetterMatch(const MatchPair& match1, const MatchPair& match2) {
0008   // Prefer DT over CSC simply because it's closer to IP
0009   // and will have less multiple scattering (at least for
0010   // RB1 vs ME1/3 case). RB1 & ME1/2 overlap is tiny
0011   if (match2.first->detector() == MuonSubdetId::DT and match1.first->detector() != MuonSubdetId::DT)
0012     return match2;
0013
0014   // For the rest compare local x match. We expect that
0015   // segments belong to the muon, so the difference in
0016   // local x is a reflection on how well we can measure it
0017   if (abs(match1.first->x - match1.second->x) > abs(match2.first->x - match2.second->x))
0018     return match2;
0019
0020   return match1;
0021 }
0022
0023 float dX(const MatchPair& match) {
0024   if (match.first and match.second->hasPhi())
0025     return (match.first->x - match.second->x);
0026   else
0027     return 9999.;
0028 }
0029
0030 float pullX(const MatchPair& match) {
0031   if (match.first and match.second->hasPhi())
0032     return dX(match) / sqrt(pow(match.first->xErr, 2) + pow(match.second->xErr, 2));
0033   else
0034     return 9999.;
0035 }
0036
0037 float pullDxDz(const MatchPair& match) {
0038   if (match.first and match.second->hasPhi())
0039     return (match.first->dXdZ - match.second->dXdZ) /
0040            sqrt(pow(match.first->dXdZErr, 2) + pow(match.second->dXdZErr, 2));
0041   else
0042     return 9999.;
0043 }
0044
0045 float dY(const MatchPair& match) {
0046   if (match.first and match.second->hasZed())
0047     return (match.first->y - match.second->y);
0048   else
0049     return 9999.;
0050 }
0051
0052 float pullY(const MatchPair& match) {
0053   if (match.first and match.second->hasZed())
0054     return dY(match) / sqrt(pow(match.first->yErr, 2) + pow(match.second->yErr, 2));
0055   else
0056     return 9999.;
0057 }
0058
0059 float pullDyDz(const MatchPair& match) {
0060   if (match.first and match.second->hasZed())
0061     return (match.first->dYdZ - match.second->dYdZ) /
0062            sqrt(pow(match.first->dYdZErr, 2) + pow(match.second->dYdZErr, 2));
0063   else
0064     return 9999.;
0065 }
0066
0067 void fillMatchInfoForStation(std::string prefix, pat::XGBooster& booster, const MatchPair& match) {
0068   booster.set(prefix + "_dX", dX(match));
0069   booster.set(prefix + "_pullX", pullX(match));
0070   booster.set(prefix + "_pullDxDz", pullDxDz(match));
0071   booster.set(prefix + "_dY", dY(match));
0072   booster.set(prefix + "_pullY", pullY(match));
0073   booster.set(prefix + "_pullDyDz", pullDyDz(match));
0074 }
0075
0076 void fillMatchInfo(pat::XGBooster& booster, const pat::Muon& muon) {
0077   // Initiate containter for results
0078   const int n_stations = 2;
0079   std::vector<MatchPair> matches;
0080   for (unsigned int i = 0; i < n_stations; ++i)
0081     matches.push_back(std::pair(nullptr, nullptr));
0082
0083   // Find best matches
0084   for (auto& chamberMatch : muon.matches()) {
0085     unsigned int station = chamberMatch.station() - 1;
0086     if (station >= n_stations)
0087       continue;
0088
0089     // Find best segment match.
0090     // We could consider all segments, but we will restrict to segments
0091     // that match to this candidate better than to other muon candidates
0092     for (auto& segmentMatch : chamberMatch.segmentMatches) {
0093       if (not segmentMatch.isMask(reco::MuonSegmentMatch::BestInStationByDR) ||
0094           not segmentMatch.isMask(reco::MuonSegmentMatch::BelongsToTrackByDR))
0095         continue;
0096
0097       // Multiple segment matches are possible in different
0098       // chambers that are either overlapping or belong to
0099       // different detectors. We need to select one.
0100       auto match_pair = MatchPair(&chamberMatch, &segmentMatch);
0101
0102       if (matches[station].first)
0103         matches[station] = getBetterMatch(matches[station], match_pair);
0104       else
0105         matches[station] = match_pair;
0106     }
0107   }
0108
0109   // Fill matching information
0110   fillMatchInfoForStation("match1", booster, matches[0]);
0111   fillMatchInfoForStation("match2", booster, matches[1]);
0112 }
0113
0114 float pat::computeSoftMvaRun3(pat::XGBooster& booster, const pat::Muon& muon) {
0115   if (!muon.isTrackerMuon() && !muon.isGlobalMuon())
0116     return 0;
0117
0118   fillMatchInfo(booster, muon);
0119
0120   booster.set("pt", muon.pt());
0121   booster.set("eta", muon.eta());
0122   booster.set("trkValidFrac", muon.innerTrack()->validFraction());
0123   booster.set("glbTrackProbability", muon.combinedQuality().glbTrackProbability);
0124   booster.set("nLostHitsInner",
0125               muon.innerTrack()->hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_INNER_HITS));
0126   booster.set("nLostHitsOuter",
0127               muon.innerTrack()->hitPattern().numberOfLostTrackerHits(reco::HitPattern::MISSING_OUTER_HITS));
0128   booster.set("trkKink", muon.combinedQuality().trkKink);
0129   booster.set("chi2LocalPosition", muon.combinedQuality().chi2LocalPosition);
0130   booster.set("nPixels", muon.innerTrack()->hitPattern().numberOfValidPixelHits());
0131   booster.set("nValidHits", muon.innerTrack()->hitPattern().numberOfValidTrackerHits());
0132   booster.set("nLostHitsOn", muon.innerTrack()->hitPattern().numberOfLostTrackerHits(reco::HitPattern::TRACK_HITS));
0133   booster.set("glbNormChi2", muon.isGlobalMuon() ? muon.globalTrack()->normalizedChi2() : 9999.);
0134   booster.set("trkLayers", muon.innerTrack()->hitPattern().trackerLayersWithMeasurement());
0135   booster.set("highPurity", muon.innerTrack()->quality(reco::Track::highPurity));
0136
0137   return booster.predict();
0138 }