TrackStatePropagator.h

Go to the documentation of this file.

#ifndef TRACKSTATEPROPAGATOR_H
#define TRACKSTATEPROPAGATOR_H

#include <memory>

#include "lardata/RecoObjects/TrackState.h"
#include "fhiclcpp/ParameterSet.h"
#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Table.h"

namespace detinfo {
  class DetectorProperties;
  class LArProperties;
}

namespace trkf {


  class TrackStatePropagator
  {
  public:

    using Plane      = recob::tracking::Plane;
    using Point_t    = recob::tracking::Point_t;
    using Vector_t   = recob::tracking::Vector_t;
    using SVector6   = recob::tracking::SVector6;

    struct Config {
      using Name = fhicl::Name;
      using Comment = fhicl::Comment;
      fhicl::Atom<double> minStep {
        Name("minStep"),
        Comment("Minimum propagation step length guaranteed."),
        1.0
       };
      fhicl::Atom<double> maxElossFrac {
        Name("maxElossFrac"),
        Comment("Maximum propagation step length based on fraction of energy loss."),
        0.1
       };
      fhicl::Atom<int> maxNit {
        Name("maxNit"),
        Comment("Maximum number of iterations when applying material effects."),
        10
       };
      fhicl::Atom<double> tcut {
        Name("tcut"),
        Comment("Maximum delta ray energy for dE/dx."),
        10.
       };
      fhicl::Atom<double> wrongDirDistTolerance {
        Name("wrongDirDistTolerance"),
        Comment("Allowed propagation distance in the wrong direction."),
        0.01
       };
      fhicl::Atom<bool> propPinvErr {
        Name("propPinvErr"),
        Comment("Propagate error on 1/p or not (in order to avoid infs, it should be set to false when 1/p not updated)."),
        false
       };
    };
    using Parameters = fhicl::Table<Config>;

    enum PropDirection {FORWARD=0, BACKWARD=1, UNKNOWN=2};

    TrackStatePropagator(double minStep, double maxElossFrac, int maxNit, double tcut, double wrongDirDistTolerance, bool propPinvErr);

    explicit TrackStatePropagator(Parameters const & p) : TrackStatePropagator(p().minStep(),p().maxElossFrac(),p().maxNit(),p().tcut(),p().wrongDirDistTolerance(),p().propPinvErr()) {}

    virtual ~TrackStatePropagator();

    TrackState propagateToPlane(bool& success, const TrackState& origin, const Plane& target, bool dodedx, bool domcs, PropDirection dir = FORWARD) const;

    inline TrackState rotateToPlane(bool& success, const TrackState& origin, const Plane& target) const { double dw2dw1 = 0.; return rotateToPlane(success, origin, target, dw2dw1);}

    inline Point_t propagatedPosByDistance(const Point_t& origpos, const Vector_t& origdir, double distance) const { return origpos+distance*origdir; }

    double distanceToPlane(bool& success, const Point_t& origpos, const Vector_t& origdir, const Plane& target) const;
    inline double distanceToPlane(bool& success, const TrackState& origin, const Plane& target) const {
      return distanceToPlane(success, origin.position(), origin.momentum().Unit(), target);
    }

    double perpDistanceToPlane(bool& success, const Point_t& origpos, const Plane& target) const;
    inline double perpDistanceToPlane(bool& success, const TrackState& origin, const Plane& target) const {
      return perpDistanceToPlane(success, origin.position(), target);
    }

    std::pair<double, double> distancePairToPlane(bool& success, const Point_t& origpos, const Vector_t& origdir, const Plane& target) const;
    inline std::pair<double, double> distancePairToPlane(bool& success, const TrackState& origin, const Plane& target) const {
      return distancePairToPlane(success, origin.position(), origin.momentum().Unit(), target);
    }

    void apply_dedx(double& pinv, double dedx, double e1, double mass, double s, double& deriv) const;

    bool apply_mcs(double dudw, double dvdw, double pinv, double mass, double s, double range, double p, double e2, bool flipSign, SMatrixSym55& noise_matrix) const;

    double getTcut() const {return fTcut;}

  private:

    TrackState rotateToPlane(bool& success, const TrackState& origin, const Plane& target, double& dw2dw1) const;

    double fMinStep;               
    double fMaxElossFrac;          
    int    fMaxNit;                
    double fTcut;                  
    double fWrongDirDistTolerance; 
    bool   fPropPinvErr;           
    const detinfo::DetectorProperties* detprop;
    const detinfo::LArProperties* larprop;
  };
}

#endif