TrackBase.h

Go to the documentation of this file.

#ifndef TrackReco_TrackBase_h
#define TrackReco_TrackBase_h

#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "DataFormats/Math/interface/Vector.h"
#include "DataFormats/Math/interface/Error.h"
#include "DataFormats/Math/interface/Vector3D.h"
#include "DataFormats/Math/interface/Point3D.h"
#include "DataFormats/Math/interface/Error.h"
#include <bitset>

namespace reco {

  class TrackBase {
  public:
    enum { dimension = 5 };

    enum { covarianceSize = dimension * (dimension + 1) / 2 };

    typedef math::Vector<dimension>::type ParameterVector;

    typedef math::Error<dimension>::type CovarianceMatrix;

    typedef math::XYZVector Vector;

    typedef math::XYZPoint Point;

    enum { i_qoverp = 0, i_lambda, i_phi, i_dxy, i_dsz };

    typedef unsigned int index;

    enum TrackAlgorithm {
      undefAlgorithm = 0,
      ctf = 1,
      duplicateMerge = 2,
      cosmics = 3,
      initialStep = 4,
      lowPtTripletStep = 5,
      pixelPairStep = 6,
      detachedTripletStep = 7,
      mixedTripletStep = 8,
      pixelLessStep = 9,
      tobTecStep = 10,
      jetCoreRegionalStep = 11,
      conversionStep = 12,
      muonSeededStepInOut = 13,
      muonSeededStepOutIn = 14,
      outInEcalSeededConv = 15,
      inOutEcalSeededConv = 16,
      nuclInter = 17,
      standAloneMuon = 18,
      globalMuon = 19,
      cosmicStandAloneMuon = 20,
      cosmicGlobalMuon = 21,
      // Phase1
      highPtTripletStep = 22,
      lowPtQuadStep = 23,
      detachedQuadStep = 24,
      displacedGeneralStep = 25,
      displacedRegionalStep = 26,
      bTagGhostTracks = 27,
      beamhalo = 28,
      gsf = 29,
      // HLT algo name
      hltPixel = 30,
      // steps used by PF
      hltIter0 = 31,
      hltIter1 = 32,
      hltIter2 = 33,
      hltIter3 = 34,
      hltIter4 = 35,
      // steps used by all other objects @HLT
      hltIterX = 36,
      // steps used by HI muon regional iterative tracking
      hiRegitMuInitialStep = 37,
      hiRegitMuLowPtTripletStep = 38,
      hiRegitMuPixelPairStep = 39,
      hiRegitMuDetachedTripletStep = 40,
      hiRegitMuMixedTripletStep = 41,
      hiRegitMuPixelLessStep = 42,
      hiRegitMuTobTecStep = 43,
      hiRegitMuMuonSeededStepInOut = 44,
      hiRegitMuMuonSeededStepOutIn = 45,
      algoSize = 46
    };

    typedef std::bitset<algoSize> AlgoMask;

    static const std::string algoNames[];

    enum TrackQuality {
      undefQuality = -1,
      loose = 0,
      tight = 1,
      highPurity = 2,
      confirmed = 3,      // means found by more than one iteration
      goodIterative = 4,  // meaningless
      looseSetWithPV = 5,
      highPuritySetWithPV = 6,
      discarded = 7,  // because a better track found. kept in the collection for reference....
      qualitySize = 8
    };

    static const std::string qualityNames[];

    TrackBase();

    TrackBase(double chi2,
              double ndof,
              const Point &vertex,
              const Vector &momentum,
              int charge,
              const CovarianceMatrix &cov,
              TrackAlgorithm = undefAlgorithm,
              TrackQuality quality = undefQuality,
              signed char nloops = 0,
              uint8_t stopReason = 0,
              float t0 = 0.f,
              float beta = 0.f,
              float covt0t0 = -1.f,
              float covbetabeta = -1.f);

    virtual ~TrackBase();

    bool isTimeOk() const { return covt0t0_ > 0.f; }

    double chi2() const;

    double ndof() const;

    double normalizedChi2() const;

    int charge() const;

    double qoverp() const;

    double theta() const;

    double lambda() const;

    double dxy() const;

    double d0() const;

    double dsz() const;

    double dz() const;

    double p2() const;

    double p() const;

    double pt2() const;

    double pt() const;

    double px() const;

    double py() const;

    double pz() const;

    double phi() const;

    double eta() const;

    double vx() const;

    double vy() const;

    double vz() const;

    const Vector &momentum() const;

    const Point &referencePoint() const;

    double t0() const;

    double beta() const;

    const Point &vertex() const;
    //__attribute__((deprecated("This method is DEPRECATED, please use referencePoint() instead.")));

    double dxy(const Point &myBeamSpot) const;

    double dxy(const BeamSpot &theBeamSpot) const;

    double dsz(const Point &myBeamSpot) const;

    double dz(const Point &myBeamSpot) const;

    ParameterVector parameters() const;

    CovarianceMatrix covariance() const;

    double parameter(int i) const;

    double covariance(int i, int j) const;

    double covt0t0() const;

    double covBetaBeta() const;

    double error(int i) const;

    double qoverpError() const;

    double ptError2() const;

    double ptError() const;

    double thetaError() const;

    double lambdaError() const;

    double etaError() const;

    double phiError() const;

    double dxyError() const;

    double d0Error() const;

    double dszError() const;

    double dzError() const;

    double t0Error() const;

    double betaError() const;

    double dxyError(Point const &vtx, math::Error<3>::type const &vertexCov) const;

    double dxyError(const BeamSpot &theBeamSpot) const;

    CovarianceMatrix &fill(CovarianceMatrix &v) const;

    static index covIndex(index i, index j);

    const HitPattern &hitPattern() const;

    unsigned short numberOfValidHits() const;

    unsigned short numberOfLostHits() const;

    int missingInnerHits() const;

    int missingOuterHits() const;

    double validFraction() const;

    template <typename C>
    bool appendHits(const C &c, const TrackerTopology &ttopo);

    template <typename I>
    bool appendHits(const I &begin, const I &end, const TrackerTopology &ttopo);

    bool appendHitPattern(const TrackingRecHit &hit, const TrackerTopology &ttopo);
    bool appendHitPattern(const DetId &id, TrackingRecHit::Type hitType, const TrackerTopology &ttopo);

    bool appendTrackerHitPattern(uint16_t subdet, uint16_t layer, uint16_t stereo, TrackingRecHit::Type hitType);
    bool appendHitPattern(const uint16_t pattern, TrackingRecHit::Type hitType);

    bool appendMuonHitPattern(const DetId &id, TrackingRecHit::Type hitType);

    void resetHitPattern();

    void setAlgorithm(const TrackAlgorithm a);

    void setOriginalAlgorithm(const TrackAlgorithm a);

    void setAlgoMask(AlgoMask a) { algoMask_ = a; }

    AlgoMask algoMask() const { return algoMask_; }
    unsigned long long algoMaskUL() const { return algoMask().to_ullong(); }
    bool isAlgoInMask(TrackAlgorithm a) const { return algoMask()[a]; }

    TrackAlgorithm algo() const;
    TrackAlgorithm originalAlgo() const;

    std::string algoName() const;

    static std::string algoName(TrackAlgorithm);

    static TrackAlgorithm algoByName(const std::string &name);

    bool quality(const TrackQuality) const;

    void setQuality(const TrackQuality);

    static std::string qualityName(TrackQuality);

    static TrackQuality qualityByName(const std::string &name);

    int qualityMask() const;

    void setQualityMask(int qualMask);

    void setNLoops(signed char value);

    bool isLooper() const;

    signed char nLoops() const;

    void setStopReason(uint8_t value) { stopReason_ = value; }

    uint8_t stopReason() const { return stopReason_; }

  private:
    HitPattern hitPattern_;

    float covariance_[covarianceSize];

    float covt0t0_, covbetabeta_;

    float chi2_;

    Point vertex_;

    float t0_;

    Vector momentum_;

    float beta_;

    std::bitset<algoSize> algoMask_;

    float ndof_;

    char charge_;

    uint8_t algorithm_;

    uint8_t originalAlgorithm_;

    uint8_t quality_;

    // I use signed char because I don't expect more than 128 loops and I could use a negative value for a special purpose.
    signed char nLoops_;

    uint8_t stopReason_;
  };

  //  Access the hit pattern, indicating in which Tracker layers the track has hits.
  inline const HitPattern &TrackBase::hitPattern() const { return hitPattern_; }

  inline bool TrackBase::appendHitPattern(const DetId &id, TrackingRecHit::Type hitType, const TrackerTopology &ttopo) {
    return hitPattern_.appendHit(id, hitType, ttopo);
  }

  inline bool TrackBase::appendHitPattern(const TrackingRecHit &hit, const TrackerTopology &ttopo) {
    return hitPattern_.appendHit(hit, ttopo);
  }

  inline bool TrackBase::appendTrackerHitPattern(uint16_t subdet,
                                                 uint16_t layer,
                                                 uint16_t stereo,
                                                 TrackingRecHit::Type hitType) {
    return hitPattern_.appendTrackerHit(subdet, layer, stereo, hitType);
  }

  inline bool TrackBase::appendHitPattern(uint16_t pattern, TrackingRecHit::Type hitType) {
    return hitPattern_.appendHit(pattern, hitType);
  }

  inline bool TrackBase::appendMuonHitPattern(const DetId &id, TrackingRecHit::Type hitType) {
    return hitPattern_.appendMuonHit(id, hitType);
  }

  inline void TrackBase::resetHitPattern() { hitPattern_.clear(); }

  template <typename I>
  bool TrackBase::appendHits(const I &begin, const I &end, const TrackerTopology &ttopo) {
    return hitPattern_.appendHits(begin, end, ttopo);
  }

  template <typename C>
  bool TrackBase::appendHits(const C &c, const TrackerTopology &ttopo) {
    return hitPattern_.appendHits(c.begin(), c.end(), ttopo);
  }

  inline TrackBase::index TrackBase::covIndex(index i, index j) {
    int a = (i <= j ? i : j);
    int b = (i <= j ? j : i);
    return b * (b + 1) / 2 + a;
  }

  inline TrackBase::TrackAlgorithm TrackBase::algo() const { return (TrackAlgorithm)(algorithm_); }
  inline TrackBase::TrackAlgorithm TrackBase::originalAlgo() const { return (TrackAlgorithm)(originalAlgorithm_); }

  inline std::string TrackBase::algoName() const { return TrackBase::algoName(algo()); }

  inline bool TrackBase::quality(const TrackBase::TrackQuality q) const {
    switch (q) {
      case undefQuality:
        return quality_ == 0;
      case goodIterative:
        return (quality_ & (1 << TrackBase::highPurity)) >> TrackBase::highPurity;
      default:
        return (quality_ & (1 << q)) >> q;
    }
    return false;
  }

  inline void TrackBase::setQuality(const TrackBase::TrackQuality q) {
    if (q == undefQuality) {
      quality_ = 0;
    } else {
      quality_ |= (1 << q);
    }
  }

  inline std::string TrackBase::qualityName(TrackQuality q) {
    if (int(q) < int(qualitySize) && int(q) >= 0) {
      return qualityNames[int(q)];
    }
    return "undefQuality";
  }

  inline std::string TrackBase::algoName(TrackAlgorithm a) {
    if (int(a) < int(algoSize) && int(a) > 0) {
      return algoNames[int(a)];
    }
    return "undefAlgorithm";
  }

  // chi-squared of the fit
  inline double TrackBase::chi2() const { return chi2_; }

  // number of degrees of freedom of the fit
  inline double TrackBase::ndof() const { return ndof_; }

  // chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
  inline double TrackBase::normalizedChi2() const { return ndof_ != 0 ? chi2_ / ndof_ : chi2_ * 1e6; }

  // track electric charge
  inline int TrackBase::charge() const { return charge_; }

  // q / p
  inline double TrackBase::qoverp() const { return charge() / p(); }

  // polar angle
  inline double TrackBase::theta() const { return momentum_.theta(); }

  // Lambda angle
  inline double TrackBase::lambda() const { return M_PI_2 - momentum_.theta(); }

  // dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close to (0,0,0): see parametrization definition above for details). See also function dxy(myBeamSpot) below.
  inline double TrackBase::dxy() const { return (-vx() * py() + vy() * px()) / pt(); }

  // dxy parameter in perigee convention (d0 = -dxy)
  inline double TrackBase::d0() const { return -dxy(); }

  // dsz parameter (THIS IS NOT the SZ impact parameter to (0,0,0) if refPoint is far from (0,0,0): see parametrization definition above for details)
  inline double TrackBase::dsz() const {
    const auto thept = pt();
    const auto thepinv = 1 / p();
    const auto theptoverp = thept * thepinv;
    return vz() * theptoverp - (vx() * px() + vy() * py()) / thept * pz() * thepinv;
  }

  // dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to (0,0,0). See also function dz(myBeamSpot) below.
  inline double TrackBase::dz() const {
    const auto thept2inv = 1 / pt2();
    return vz() - (vx() * px() + vy() * py()) * pz() * thept2inv;
  }

  // momentum vector magnitude square
  inline double TrackBase::p2() const { return momentum_.Mag2(); }

  // momentum vector magnitude
  inline double TrackBase::p() const { return sqrt(p2()); }

  // track transverse momentum square
  inline double TrackBase::pt2() const { return momentum_.Perp2(); }

  // track transverse momentum
  inline double TrackBase::pt() const { return sqrt(pt2()); }

  // x coordinate of momentum vector
  inline double TrackBase::px() const { return momentum_.x(); }

  // y coordinate of momentum vector
  inline double TrackBase::py() const { return momentum_.y(); }

  // z coordinate of momentum vector
  inline double TrackBase::pz() const { return momentum_.z(); }

  // azimuthal angle of momentum vector
  inline double TrackBase::phi() const { return momentum_.Phi(); }

  // pseudorapidity of momentum vector
  inline double TrackBase::eta() const { return momentum_.Eta(); }

  // x coordinate of the reference point on track
  inline double TrackBase::vx() const { return vertex_.x(); }

  // y coordinate of the reference point on track
  inline double TrackBase::vy() const { return vertex_.y(); }

  // z coordinate of the reference point on track
  inline double TrackBase::vz() const { return vertex_.z(); }

  // track momentum vector
  inline const TrackBase::Vector &TrackBase::momentum() const { return momentum_; }

  // Reference point on the track
  inline const TrackBase::Point &TrackBase::referencePoint() const { return vertex_; }

  // Time at the reference point on the track
  inline double TrackBase::t0() const { return t0_; }

  // Velocity at the reference point on the track in natural units
  inline double TrackBase::beta() const { return beta_; }

  // reference point on the track. This method is DEPRECATED, please use referencePoint() instead
  inline const TrackBase::Point &TrackBase::vertex() const { return vertex_; }

  // dxy parameter with respect to a user-given beamSpot
  // (WARNING: this quantity can only be interpreted as a minimum transverse distance if beamSpot, if the beam spot is reasonably close to the refPoint, since linear approximations are involved).
  // This is a good approximation for Tracker tracks.
  inline double TrackBase::dxy(const Point &myBeamSpot) const {
    return (-(vx() - myBeamSpot.x()) * py() + (vy() - myBeamSpot.y()) * px()) / pt();
  }

  // dxy parameter with respect to the beamSpot taking into account the beamspot slopes
  // (WARNING: this quantity can only be interpreted as a minimum transverse distance if beamSpot, if the beam spot is reasonably close to the refPoint, since linear approximations are involved).
  // This is a good approximation for Tracker tracks.
  inline double TrackBase::dxy(const BeamSpot &theBeamSpot) const { return dxy(theBeamSpot.position(vz())); }

  // dsz parameter with respect to a user-given beamSpot
  // (WARNING: this quantity can only be interpreted as the distance in the S-Z plane to the beamSpot, if the beam spot is reasonably close to the refPoint, since linear approximations are involved).
  // This is a good approximation for Tracker tracks.
  inline double TrackBase::dsz(const Point &myBeamSpot) const {
    const auto thept = pt();
    const auto thepinv = 1 / p();
    const auto theptoverp = thept * thepinv;
    return (vz() - myBeamSpot.z()) * theptoverp -
           ((vx() - myBeamSpot.x()) * px() + (vy() - myBeamSpot.y()) * py()) / thept * pz() * thepinv;
  }

  // dz parameter with respect to a user-given beamSpot
  // (WARNING: this quantity can only be interpreted as the track z0, if the beamSpot is reasonably close to the refPoint, since linear approximations are involved).
  // This is a good approximation for Tracker tracks.
  inline double TrackBase::dz(const Point &myBeamSpot) const {
    const auto theptinv2 = 1 / pt2();
    return (vz() - myBeamSpot.z()) -
           ((vx() - myBeamSpot.x()) * px() + (vy() - myBeamSpot.y()) * py()) * pz() * theptinv2;
  }

  // Track parameters with one-to-one correspondence to the covariance matrix
  inline TrackBase::ParameterVector TrackBase::parameters() const {
    return TrackBase::ParameterVector(qoverp(), lambda(), phi(), dxy(), dsz());
  }

  // return track covariance matrix
  inline TrackBase::CovarianceMatrix TrackBase::covariance() const {
    CovarianceMatrix m;
    fill(m);
    return m;
  }

  // i-th parameter ( i = 0, ... 4 )
  inline double TrackBase::parameter(int i) const { return parameters()[i]; }

  // (i,j)-th element of covariance matrix (i, j = 0, ... 4)
  inline double TrackBase::covariance(int i, int j) const { return covariance_[covIndex(i, j)]; }

  // error on specified element
  inline double TrackBase::error(int i) const { return sqrt(covariance_[covIndex(i, i)]); }

  // error on signed transverse curvature
  inline double TrackBase::qoverpError() const { return error(i_qoverp); }

  // error on Pt (set to 1000**2 TeV**2 if charge==0 for safety)
  inline double TrackBase::ptError2() const {
    const auto thecharge = charge();

    if (thecharge != 0) {
      const auto thept2 = pt2();
      const auto thep2 = p2();
      const auto thepz = pz();
      const auto ptimespt = sqrt(thep2 * thept2);
      const auto oneovercharge = 1 / thecharge;

      return thept2 * thep2 * oneovercharge * oneovercharge * covariance(i_qoverp, i_qoverp) +
             2 * ptimespt * oneovercharge * thepz * covariance(i_qoverp, i_lambda) +
             thepz * thepz * covariance(i_lambda, i_lambda);
    }

    return 1.e12;
  }

  // error on Pt (set to 1000 TeV if charge==0 for safety)
  inline double TrackBase::ptError() const { return sqrt(ptError2()); }

  // error on theta
  inline double TrackBase::thetaError() const { return error(i_lambda); }

  // error on lambda
  inline double TrackBase::lambdaError() const { return error(i_lambda); }

  // error on eta
  inline double TrackBase::etaError() const { return error(i_lambda) * sqrt(p2() / pt2()); }

  // error on phi
  inline double TrackBase::phiError() const { return error(i_phi); }

  // error on dxy
  inline double TrackBase::dxyError() const { return error(i_dxy); }

  // error on d0
  inline double TrackBase::d0Error() const { return error(i_dxy); }

  // error on dsz
  inline double TrackBase::dszError() const { return error(i_dsz); }

  // error on dz
  inline double TrackBase::dzError() const { return error(i_dsz) * sqrt(p2() / pt2()); }

  // covariance of t0
  inline double TrackBase::covt0t0() const { return covt0t0_; }

  // covariance of beta
  inline double TrackBase::covBetaBeta() const { return covbetabeta_; }

  // error on t0
  inline double TrackBase::t0Error() const { return std::sqrt(covt0t0_); }

  // error on beta
  inline double TrackBase::betaError() const { return std::sqrt(covbetabeta_); }

  // error on dxy with respect to a given beamspot
  inline double TrackBase::dxyError(const BeamSpot &theBeamSpot) const {
    return dxyError(theBeamSpot.position(vz()), theBeamSpot.rotatedCovariance3D());
  }

  // number of valid hits found
  inline unsigned short TrackBase::numberOfValidHits() const { return hitPattern_.numberOfValidHits(); }

  // number of cases where track crossed a layer without getting a hit.
  inline unsigned short TrackBase::numberOfLostHits() const {
    return hitPattern_.numberOfLostHits(HitPattern::TRACK_HITS);
  }

  // number of hits expected from inner track extrapolation but missing
  inline int TrackBase::missingInnerHits() const {
    return hitPattern_.numberOfLostHits(HitPattern::MISSING_INNER_HITS);
  }

  // number of hits expected from outer track extrapolation but missing
  inline int TrackBase::missingOuterHits() const {
    return hitPattern_.numberOfLostHits(HitPattern::MISSING_OUTER_HITS);
  }

  // fraction of valid hits on the track
  inline double TrackBase::validFraction() const {
    int valid = hitPattern_.numberOfValidTrackerHits();
    int lost = hitPattern_.numberOfLostTrackerHits(HitPattern::TRACK_HITS);
    int lostIn = hitPattern_.numberOfLostTrackerHits(HitPattern::MISSING_INNER_HITS);
    int lostOut = hitPattern_.numberOfLostTrackerHits(HitPattern::MISSING_OUTER_HITS);

    const auto tot = valid + lost + lostIn + lostOut;

    if (tot == 0) {
      return -1;
    }

    return valid / (double)(tot);
  }

  //Track algorithm
  inline void TrackBase::setAlgorithm(const TrackBase::TrackAlgorithm a) {
    algorithm_ = a;
    algoMask_.reset();
    setOriginalAlgorithm(a);
  }

  inline void TrackBase::setOriginalAlgorithm(const TrackBase::TrackAlgorithm a) {
    originalAlgorithm_ = a;
    algoMask_.set(a);
  }

  inline int TrackBase::qualityMask() const { return quality_; }

  inline void TrackBase::setQualityMask(int qualMask) { quality_ = qualMask; }

  inline void TrackBase::setNLoops(signed char value) { nLoops_ = value; }

  inline bool TrackBase::isLooper() const { return (nLoops_ > 0); }

  inline signed char TrackBase::nLoops() const { return nLoops_; }

}  // namespace reco

#endif