MuonCocktails.cc

Go to the documentation of this file.

#include "TMath.h"
#include "DataFormats/MuonReco/interface/MuonCocktails.h"
#include "DataFormats/TrackReco/interface/Track.h"

//
// Return the TeV-optimized refit track (aka the cocktail or Tune P) or
// the tracker track if either the optimized pT or tracker pT is below the pT threshold
//
reco::Muon::MuonTrackTypePair  muon::tevOptimized(const reco::TrackRef& combinedTrack,
                          const reco::TrackRef& trackerTrack,
                          const reco::TrackRef& tpfmsTrack,
                          const reco::TrackRef& pickyTrack,
                          const reco::TrackRef& dytTrack,
                          const double ptThreshold,
                          const double tune1,
                          const double tune2,
                          double dptcut) {
  const unsigned int nAlgo=5; 

  // Array for convenience below.
  const reco::Muon::MuonTrackTypePair refit[nAlgo] = { 
    make_pair(trackerTrack, reco::Muon::InnerTrack), 
    make_pair(combinedTrack,reco::Muon::CombinedTrack),
    make_pair(tpfmsTrack,   reco::Muon::TPFMS),
    make_pair(pickyTrack,   reco::Muon::Picky),
    make_pair(dytTrack,     reco::Muon::DYT)
  }; 
  
  // Calculate the log(tail probabilities). If there's a problem,
  // signify this with prob == 0. The current problems recognized are:
  // the track being not available, whether the (re)fit failed or it's
  // just not in the event, or if the (re)fit ended up with no valid
  // hits.
  double prob[nAlgo] = {0.,0.,0.,0.,0.};
  bool valid[nAlgo] = {false,false,false,false,false};

  double dptmin = 1.;

  if (dptcut>0) {  
    for (unsigned int i = 0; i < nAlgo; ++i)
      if (refit[i].first.isNonnull())
        if (refit[i].first->ptError()/refit[i].first->pt()<dptmin) dptmin = refit[i].first->ptError()/refit[i].first->pt();
  
    if (dptmin>dptcut) dptcut = dptmin+0.15;
  }

  for (unsigned int i = 0; i < nAlgo; ++i) 
    if (refit[i].first.isNonnull()){ 
      valid[i] = true;
      if (refit[i].first->numberOfValidHits() && (refit[i].first->ptError()/refit[i].first->pt()<dptcut || dptcut<0)) 
    prob[i] = muon::trackProbability(refit[i].first); 
    }

  
  // Start with picky.
  int chosen = 3;
  
  // If there's a problem with picky, make the default one of the
  // other tracks. Try TPFMS first, then global, then tracker-only.
  if (prob[3] == 0.) { 

    // split so that passing dptcut<0 recreates EXACTLY the old tuneP behavior
    if (dptcut>0) {
      if      (prob[4] > 0.) chosen = 4;
      else if (prob[0] > 0.) chosen = 0;
      else if (prob[2] > 0.) chosen = 2;
      else if (prob[1] > 0.) chosen = 1;
    } else {
      if      (prob[2] > 0.) chosen = 2;
      else if (prob[1] > 0.) chosen = 1;
      else if (prob[0] > 0.) chosen = 0;
    }
  } 
  
  // Now the algorithm: switch from picky to dyt if the difference is lower than a tuned value. Then 
  //  switch from picky to tracker-only if the
  // difference, log(tail prob(picky)) - log(tail prob(tracker-only))
  // is greater than a tuned value. Then compare the
  // so-picked track to TPFMS in the same manner using another tuned
  // value.
  if (prob[4]>0. && prob[3]>0.) {
    if(refit[3].first->pt()>0 && refit[4].first->pt()>0 &&
       (refit[4].first->ptError()/refit[4].first->pt()-refit[3].first->ptError()/refit[3].first->pt())<=0) 
      chosen=4; // dyt
  }

  if (prob[0] > 0. && prob[chosen] > 0. && (prob[chosen] - prob[0]) > tune1)
    chosen = 0;
  if (prob[2] > 0. && (prob[chosen] - prob[2]) > tune2)
    chosen = 2;

  // Sanity checks 
  if (chosen == 4 && !valid[4] ) chosen = 3;
  if (chosen == 3 && !valid[3] ) chosen = 2;
  if (chosen == 2 && !valid[2] ) chosen = 1;
  if (chosen == 1 && !valid[1] ) chosen = 0; 

  // Done. If pT of the chosen track (or pT of the tracker track) is below the threshold value, return the tracker track.
  if (valid[chosen] && refit[chosen].first->pt() < ptThreshold && prob[0] > 0.) return make_pair(trackerTrack,reco::Muon::InnerTrack);    
  if (trackerTrack->pt() < ptThreshold && prob[0] > 0.) return make_pair(trackerTrack,reco::Muon::InnerTrack);  
  
  // Return the chosen track (which can be the global track in
  // very rare cases).
  return refit[chosen];
}

//
// calculate the tail probability (-ln(P)) of a fit
//
double muon::trackProbability(const reco::TrackRef track) {
  
  int nDOF = (int)track->ndof();
  if ( nDOF > 0 && track->chi2()> 0) { 
    return -log(TMath::Prob(track->chi2(), nDOF));
  } else { 
    return 0.0;
  }
  
}

reco::TrackRef muon::getTevRefitTrack(const reco::TrackRef& combinedTrack,
                             const reco::TrackToTrackMap& map) {
  reco::TrackToTrackMap::const_iterator it = map.find(combinedTrack);
  return it == map.end() ? reco::TrackRef() : it->val;
}


//
// Get the sigma-switch decision (tracker-only versus global).
//
reco::Muon::MuonTrackTypePair muon::sigmaSwitch(const reco::TrackRef& combinedTrack,
                        const reco::TrackRef& trackerTrack,
                        const double nSigma,
                        const double ptThreshold) {
  // If either the global or tracker-only fits have pT below threshold
  // (default 200 GeV), return the tracker-only fit.
  if (combinedTrack->pt() < ptThreshold || trackerTrack->pt() < ptThreshold)
    return make_pair(trackerTrack,reco::Muon::InnerTrack);
  
  // If both are above the pT threshold, compare the difference in
  // q/p: if less than two sigma of the tracker-only track, switch to
  // global. Otherwise, use tracker-only.
  const double delta = fabs(trackerTrack->qoverp() - combinedTrack->qoverp());
  const double threshold = nSigma * trackerTrack->qoverpError();
  return delta > threshold ? make_pair(trackerTrack,reco::Muon::InnerTrack) :  make_pair(combinedTrack,reco::Muon::CombinedTrack);
}

//
// Get the TMR decision (tracker-only versus TPFMS).
//
reco::Muon::MuonTrackTypePair muon::TMR(const reco::TrackRef& trackerTrack,
                    const reco::TrackRef& fmsTrack,
                    const double tune) {
  double probTK  = 0;
  double probFMS = 0;
  
  if (trackerTrack.isNonnull() && trackerTrack->numberOfValidHits())  
    probTK = muon::trackProbability(trackerTrack);
  if (fmsTrack.isNonnull() && fmsTrack->numberOfValidHits())
    probFMS = muon::trackProbability(fmsTrack);
  
  bool TKok  = probTK > 0;
  bool FMSok = probFMS > 0;
  
  if (TKok && FMSok) {
    if (probFMS - probTK > tune)
      return make_pair(trackerTrack,reco::Muon::InnerTrack);
    else
      return  make_pair(fmsTrack,reco::Muon::TPFMS);
  }
  else if (FMSok)
    return  make_pair(fmsTrack,reco::Muon::TPFMS);
  else if (TKok)
    return make_pair(trackerTrack,reco::Muon::InnerTrack);
  else
    return make_pair(reco::TrackRef(),reco::Muon::None); 
}