MuonReco/src/MuonCocktails.cc

0001 #include "TMath.h"
0002 #include "DataFormats/MuonReco/interface/MuonCocktails.h"
0003 #include "DataFormats/TrackReco/interface/Track.h"
0004
0005 //
0006 // Return the TeV-optimized refit track (aka the cocktail or Tune P) or
0007 // the tracker track if either the optimized pT or tracker pT is below the pT threshold
0008 //
0009 reco::Muon::MuonTrackTypePair muon::tevOptimized(const reco::TrackRef& combinedTrack,
0010                                                  const reco::TrackRef& trackerTrack,
0011                                                  const reco::TrackRef& tpfmsTrack,
0012                                                  const reco::TrackRef& pickyTrack,
0013                                                  const reco::TrackRef& dytTrack,
0014                                                  const double ptThreshold,
0015                                                  const double tune1,
0016                                                  const double tune2,
0017                                                  double dptcut) {
0018   const unsigned int nAlgo = 5;
0019
0020   // Array for convenience below.
0021   const reco::Muon::MuonTrackTypePair refit[nAlgo] = {make_pair(trackerTrack, reco::Muon::InnerTrack),
0022                                                       make_pair(combinedTrack, reco::Muon::CombinedTrack),
0023                                                       make_pair(tpfmsTrack, reco::Muon::TPFMS),
0024                                                       make_pair(pickyTrack, reco::Muon::Picky),
0025                                                       make_pair(dytTrack, reco::Muon::DYT)};
0026
0027   // Calculate the log(tail probabilities). If there's a problem,
0028   // signify this with prob == 0. The current problems recognized are:
0029   // the track being not available, whether the (re)fit failed or it's
0030   // just not in the event, or if the (re)fit ended up with no valid
0031   // hits.
0032   double prob[nAlgo] = {0., 0., 0., 0., 0.};
0033   bool valid[nAlgo] = {false, false, false, false, false};
0034
0035   double dptmin = 1.;
0036
0037   if (dptcut > 0) {
0038     for (unsigned int i = 0; i < nAlgo; ++i)
0039       if (refit[i].first.isNonnull())
0040         if (refit[i].first->ptError() / refit[i].first->pt() < dptmin)
0041           dptmin = refit[i].first->ptError() / refit[i].first->pt();
0042
0043     if (dptmin > dptcut)
0044       dptcut = dptmin + 0.15;
0045   }
0046
0047   for (unsigned int i = 0; i < nAlgo; ++i)
0048     if (refit[i].first.isNonnull()) {
0049       valid[i] = true;
0050       if (refit[i].first->numberOfValidHits() &&
0051           (refit[i].first->ptError() / refit[i].first->pt() < dptcut || dptcut < 0))
0052         prob[i] = muon::trackProbability(refit[i].first);
0053     }
0054
0055   // Start with picky.
0056   int chosen = 3;
0057
0058   // If there's a problem with picky, make the default one of the
0059   // other tracks. Try TPFMS first, then global, then tracker-only.
0060   if (prob[3] == 0.) {
0061     // split so that passing dptcut<0 recreates EXACTLY the old tuneP behavior
0062     if (dptcut > 0) {
0063       if (prob[4] > 0.)
0064         chosen = 4;
0065       else if (prob[0] > 0.)
0066         chosen = 0;
0067       else if (prob[2] > 0.)
0068         chosen = 2;
0069       else if (prob[1] > 0.)
0070         chosen = 1;
0071     } else {
0072       if (prob[2] > 0.)
0073         chosen = 2;
0074       else if (prob[1] > 0.)
0075         chosen = 1;
0076       else if (prob[0] > 0.)
0077         chosen = 0;
0078     }
0079   }
0080
0081   // Now the algorithm: switch from picky to dyt if the difference is lower than a tuned value. Then
0082   //  switch from picky to tracker-only if the
0083   // difference, log(tail prob(picky)) - log(tail prob(tracker-only))
0084   // is greater than a tuned value. Then compare the
0085   // so-picked track to TPFMS in the same manner using another tuned
0086   // value.
0087   if (prob[4] > 0. && prob[3] > 0.) {
0088     if (refit[3].first->pt() > 0 && refit[4].first->pt() > 0 &&
0089         (refit[4].first->ptError() / refit[4].first->pt() - refit[3].first->ptError() / refit[3].first->pt()) <= 0)
0090       chosen = 4;  // dyt
0091   }
0092
0093   if (prob[0] > 0. && prob[chosen] > 0. && (prob[chosen] - prob[0]) > tune1)
0094     chosen = 0;
0095   if (prob[2] > 0. && (prob[chosen] - prob[2]) > tune2)
0096     chosen = 2;
0097
0098   // Sanity checks
0099   if (chosen == 4 && !valid[4])
0100     chosen = 3;
0101   if (chosen == 3 && !valid[3])
0102     chosen = 2;
0103   if (chosen == 2 && !valid[2])
0104     chosen = 1;
0105   if (chosen == 1 && !valid[1])
0106     chosen = 0;
0107
0108   // Done. If pT of the chosen track (or pT of the tracker track) is below the threshold value, return the tracker track.
0109   if (valid[chosen] && refit[chosen].first->pt() < ptThreshold && prob[0] > 0.)
0110     return make_pair(trackerTrack, reco::Muon::InnerTrack);
0111   if (trackerTrack->pt() < ptThreshold && prob[0] > 0.)
0112     return make_pair(trackerTrack, reco::Muon::InnerTrack);
0113
0114   // Return the chosen track (which can be the global track in
0115   // very rare cases).
0116   return refit[chosen];
0117 }
0118
0119 //
0120 // calculate the tail probability (-ln(P)) of a fit
0121 //
0122 double muon::trackProbability(const reco::TrackRef track) {
0123   int nDOF = (int)track->ndof();
0124   if (nDOF > 0 && track->chi2() > 0) {
0125     return -log(TMath::Prob(track->chi2(), nDOF));
0126   } else {
0127     return 0.0;
0128   }
0129 }
0130
0131 reco::TrackRef muon::getTevRefitTrack(const reco::TrackRef& combinedTrack, const reco::TrackToTrackMap& map) {
0132   reco::TrackToTrackMap::const_iterator it = map.find(combinedTrack);
0133   return it == map.end() ? reco::TrackRef() : it->val;
0134 }
0135
0136 //
0137 // Get the sigma-switch decision (tracker-only versus global).
0138 //
0139 reco::Muon::MuonTrackTypePair muon::sigmaSwitch(const reco::TrackRef& combinedTrack,
0140                                                 const reco::TrackRef& trackerTrack,
0141                                                 const double nSigma,
0142                                                 const double ptThreshold) {
0143   // If either the global or tracker-only fits have pT below threshold
0144   // (default 200 GeV), return the tracker-only fit.
0145   if (combinedTrack->pt() < ptThreshold || trackerTrack->pt() < ptThreshold)
0146     return make_pair(trackerTrack, reco::Muon::InnerTrack);
0147
0148   // If both are above the pT threshold, compare the difference in
0149   // q/p: if less than two sigma of the tracker-only track, switch to
0150   // global. Otherwise, use tracker-only.
0151   const double delta = fabs(trackerTrack->qoverp() - combinedTrack->qoverp());
0152   const double threshold = nSigma * trackerTrack->qoverpError();
0153   return delta > threshold ? make_pair(trackerTrack, reco::Muon::InnerTrack)
0154                            : make_pair(combinedTrack, reco::Muon::CombinedTrack);
0155 }
0156
0157 //
0158 // Get the TMR decision (tracker-only versus TPFMS).
0159 //
0160 reco::Muon::MuonTrackTypePair muon::TMR(const reco::TrackRef& trackerTrack,
0161                                         const reco::TrackRef& fmsTrack,
0162                                         const double tune) {
0163   double probTK = 0;
0164   double probFMS = 0;
0165
0166   if (trackerTrack.isNonnull() && trackerTrack->numberOfValidHits())
0167     probTK = muon::trackProbability(trackerTrack);
0168   if (fmsTrack.isNonnull() && fmsTrack->numberOfValidHits())
0169     probFMS = muon::trackProbability(fmsTrack);
0170
0171   bool TKok = probTK > 0;
0172   bool FMSok = probFMS > 0;
0173
0174   if (TKok && FMSok) {
0175     if (probFMS - probTK > tune)
0176       return make_pair(trackerTrack, reco::Muon::InnerTrack);
0177     else
0178       return make_pair(fmsTrack, reco::Muon::TPFMS);
0179   } else if (FMSok)
0180     return make_pair(fmsTrack, reco::Muon::TPFMS);
0181   else if (TKok)
0182     return make_pair(trackerTrack, reco::Muon::InnerTrack);
0183   else
0184     return make_pair(reco::TrackRef(), reco::Muon::None);
0185 }