TrackStatePropagator.h

Go to the documentation of this file.

#ifndef TRACKSTATEPROPAGATOR_H
#define TRACKSTATEPROPAGATOR_H

#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Table.h"
#include "lardata/RecoObjects/TrackState.h"
#include "lardataobj/RecoBase/TrackingPlane.h"
#include "lardataobj/RecoBase/TrackingTypes.h"

#include <utility>

namespace detinfo {
  class DetectorPropertiesData;
  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())
    {}

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

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

    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;
    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;
    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;
    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,
                    detinfo::DetectorPropertiesData const& detProp,
                    double dedx,
                    double e1,
                    double mass,
                    double s,
                    double& deriv) const;

    bool apply_mcs(detinfo::DetectorPropertiesData const& detProp,
                   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::LArProperties* larprop;
  };
}

#endif