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,
        reservedForUpgrades1 = 25, reservedForUpgrades2 = 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 p() 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 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;

    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_;

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

    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
{
    return vz() * pt() / p() - (vx() * px() + vy() * py()) / pt() * pz() / p();
}

// 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
{
    return vz() - (vx() * px() + vy() * py()) / pt() * (pz() / pt());
}

// momentum vector magnitude
inline double TrackBase::p() const
{
    return momentum_.R();
}

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

// 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
{
    return (vz() - myBeamSpot.z()) * pt() / p() - ((vx() - myBeamSpot.x()) * px() + (vy() - myBeamSpot.y()) * py()) / pt() * pz() / p();
}

// 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
{
    return (vz() - myBeamSpot.z()) - ((vx() - myBeamSpot.x()) * px() + (vy() - myBeamSpot.y()) * py()) / pt() * pz() / pt();
}

// 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 TeV if charge==0 for safety)
inline double TrackBase::ptError() const
{
    return (charge() != 0) ?  sqrt(
               pt() * pt() * p() * p() / charge() / charge() * covariance(i_qoverp, i_qoverp)
               + 2 * pt() * p() / charge() * pz() * covariance(i_qoverp, i_lambda)
               + pz() * pz() * covariance(i_lambda, i_lambda)) : 1.e6;
}

// 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) * p() / pt();
}

// 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) * p() / pt();
}

// 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);
}

// 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);

    if ((valid + lost + lostIn + lostOut) == 0) {
        return -1;
    }

    return valid / (double)(valid + lost + lostIn + lostOut);
}

//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