MuonCosmicsId.cc

Go to the documentation of this file.

#include "RecoMuon/MuonIdentification/interface/MuonCosmicsId.h"
#include "DataFormats/TrackReco/interface/Track.h"

bool directionAlongMomentum(const reco::Track& track) {
  // check is done in 2D
  return (track.innerPosition().x() - track.vx()) * track.px() + (track.innerPosition().y() - track.vy()) * track.py() >
         0;
}

// returns angle and dPt/Pt
std::pair<double, double> muonid::matchTracks(const reco::Track& ref, const reco::Track& probe) {
  std::pair<double, double> result(0, 0);
  // When both tracks are reconstructed as outside going, sign is -1
  // otherwise it's +1. There is also a crazy case of both are outside
  // going, then sign is -1 as well.
  int match_sign = directionAlongMomentum(ref) == directionAlongMomentum(probe) ? -1 : +1;
  double sprod = ref.px() * probe.px() + ref.py() * probe.py() + ref.pz() * probe.pz();
  double argCos = match_sign * (sprod / ref.p() / probe.p());
  if (argCos < -1.0)
    argCos = -1.0;
  if (argCos > 1.0)
    argCos = 1.0;
  result.first = acos(argCos);
  result.second = fabs(probe.pt() - ref.pt()) / sqrt(ref.pt() * probe.pt());  //SK: take a geom-mean pt
  return result;
}

reco::TrackRef muonid::findOppositeTrack(const edm::Handle<reco::TrackCollection>& tracks,
                                         const reco::Track& muonTrack,
                                         double angleMatch,
                                         double momentumMatch) {
  for (unsigned int i = 0; i < tracks->size(); ++i) {
    // When both tracks are reconstructed as outside going, sign is -1
    // otherwise it's +1. There is also a crazy case of both are outside
    // going, then sign is -1 as well.
    const std::pair<double, double>& match = matchTracks(muonTrack, tracks->at(i));
    if (match.first < angleMatch && match.second < momentumMatch)
      return reco::TrackRef(tracks, i);
  }
  return reco::TrackRef();
}