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