TrackState.h

Go to the documentation of this file.

#ifndef TRACKSTATE_H
#define TRACKSTATE_H

#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "lardata/RecoObjects/TrackingPlaneHelper.h"
#include "lardataobj/RecoBase/TrackingPlane.h"
#include "lardataobj/RecoBase/TrackingTypes.h"

namespace trkf {

  using SVector5 = recob::tracking::SVector5;
  using SVector6 = recob::tracking::SVector6;
  using SMatrix55 = recob::tracking::SMatrix55;
  using SMatrixSym55 = recob::tracking::SMatrixSym55;
  using SMatrixSym66 = recob::tracking::SMatrixSym66;
  using Plane = recob::tracking::Plane;
  using Point_t = recob::tracking::Point_t;
  using Vector_t = recob::tracking::Vector_t;

  namespace {
    const double elmass = 0.000510998; // Electron
    const double mumass = 0.105658367; // Muon
    const double pimass = 0.13957;     // Charged pion
    const double kmass = 0.493677;     // Charged kaon
    const double pmass = 0.938272;     // Proton
  }


  class HitState {
  public:
    HitState(double hitMeas,
             double hitMeasErr2,
             geo::WireID const& wireId,
             const geo::WireGeo& wgeom)
      : fHitMeas(hitMeas)
      , fHitMeasErr2(hitMeasErr2)
      , fWireId(wireId)
      , fPlane(recob::tracking::makePlane(wgeom))
    {}
    HitState(double hitMeas, double hitMeasErr2, geo::WireID&& wireId, const geo::WireGeo& wgeom)
      : fHitMeas(hitMeas)
      , fHitMeasErr2(hitMeasErr2)
      , fWireId(std::move(wireId))
      , fPlane(recob::tracking::makePlane(wgeom))
    {}
    double hitMeas() const { return fHitMeas; }
    double hitMeasErr2() const { return fHitMeasErr2; }
    const Plane& plane() const { return fPlane; }
    const geo::WireID& wireId() const { return fWireId; }
    std::ostream& dump(std::ostream& out = std::cout) const
    {
      out << "HitState with meas=" << hitMeas() << " err2=" << hitMeasErr2()
          << " plane=" << wireId().Plane << " wire=" << wireId().Wire
          << " on plane with pos=" << plane().position() << " and dir=" << plane().direction()
          << "\n";
      return out;
    }

  private:
    double fHitMeas;
    double fHitMeasErr2;
    const geo::WireID fWireId;
    Plane fPlane;
  };


  class TrackState {
  public:
    TrackState(const SVector5& trackStatePar,
               const SMatrixSym55& trackStateCov,
               const Plane& plane,
               bool trackAlongPlaneDir,
               int pid)
      : fTrackStatePar(trackStatePar), fTrackStateCov(trackStateCov), fPlane(plane), fPid(pid)
    {
      SVector6 par6d = fPlane.Local5DToGlobal6DParameters(fTrackStatePar, trackAlongPlaneDir);
      fPos = Point_t(par6d[0], par6d[1], par6d[2]);
      fMom = Point_t(par6d[3], par6d[4], par6d[5]);
    }
    //
    const SVector5& parameters() const { return fTrackStatePar; }
    const SMatrixSym55& covariance() const { return fTrackStateCov; }
    const Plane& plane() const { return fPlane; }
    const Point_t& position() const { return fPos; }
    const Vector_t& momentum() const { return fMom; }
    int pID() const { return fPid; }
    double mass() const
    {
      if (abs(fPid) == 11) { return elmass; }
      if (abs(fPid) == 13) { return mumass; }
      if (abs(fPid) == 211) { return pimass; }
      if (abs(fPid) == 321) { return kmass; }
      if (abs(fPid) == 2212) { return pmass; }
      return util::kBogusD;
    }
    SVector6 parameters6D() const
    {
      return SVector6(fPos.X(), fPos.Y(), fPos.Z(), fMom.X(), fMom.Y(), fMom.Z());
    }
    SMatrixSym66 covariance6D() const
    {
      return fPlane.Local5DToGlobal6DCovariance(fTrackStateCov, true, fMom);
    }
    //
    bool isTrackAlongPlaneDir() const { return fMom.Dot(fPlane.direction()) > 0; }
    //
    std::ostream& dump(std::ostream& out = std::cout) const
    {
      out << "TrackState with pID=" << pID() << " mass=" << mass() << "\npars=" << parameters()
          << " position=" << position() << " momentum=" << momentum() << "\ncov=\n"
          << covariance() << "\non plane with pos=" << plane().position()
          << " and dir=" << plane().direction() << " along=" << isTrackAlongPlaneDir() << "\n";
      return out;
    }
    //
    inline double residual(const HitState& hitstate) const
    {
      return hitstate.hitMeas() - fTrackStatePar(0);
    }

    inline double combinedError2(const HitState& hitstate) const
    {
      return hitstate.hitMeasErr2() + fTrackStateCov(0, 0);
    }

    inline double combinedError(const HitState& hitstate) const
    {
      return sqrt(combinedError2(hitstate));
    }

    inline double chi2(const HitState& hitstate) const
    {
      return residual(hitstate) * residual(hitstate) / combinedError2(hitstate);
    }

    void setCovariance(const SMatrixSym55& trackStateCov) { fTrackStateCov = trackStateCov; }

    void setParameters(const SVector5& trackStatePar)
    {
      fTrackStatePar = trackStatePar;
      SVector6 par6d = fPlane.Local5DToGlobal6DParameters(trackStatePar, isTrackAlongPlaneDir());
      fPos = Point_t(par6d[0], par6d[1], par6d[2]);
      fMom = Vector_t(par6d[3], par6d[4], par6d[5]);
    }
    //
  private:
    SVector5 fTrackStatePar;     
    SMatrixSym55 fTrackStateCov; 
    Plane fPlane;                
    int fPid;                    
    Point_t fPos;                
    Vector_t fMom;               
  };

}

#endif