trkf::Propagator Class Reference

#include "Propagator.h"

Inheritance diagram for trkf::Propagator:

Public Types
enum	PropDirection { FORWARD, BACKWARD, UNKNOWN }
	Propagation direction enum. More...

Public Member Functions
	Propagator (detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx, const std::shared_ptr< const Interactor > &interactor)
	Constructor. More...
virtual	~Propagator ()
	Destructor. More...
double	getTcut () const
bool	getDoDedx () const
const std::shared_ptr< const Interactor > &	getInteractor () const
virtual Propagator *	clone () const =0
	Clone method. More...
virtual std::optional< double >	short_vec_prop (KTrack &trk, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, TrackMatrix *prop_matrix=0, TrackError *noise_matrix=0) const =0
	Propagate without error (short distance). More...
virtual std::optional< double >	origin_vec_prop (KTrack &trk, const std::shared_ptr< const Surface > &porient, TrackMatrix *prop_matrix=0) const =0
	Propagate without error to surface whose origin parameters coincide with track position. More...
std::optional< double >	vec_prop (KTrack &trk, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, TrackMatrix *prop_matrix=0, TrackError *noise_matrix=0) const
	Propagate without error (long distance). More...
std::optional< double >	lin_prop (KTrack &trk, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, KTrack *ref=0, TrackMatrix *prop_matrix=0, TrackError *noise_matrix=0) const
	Linearized propagate without error. More...
std::optional< double >	err_prop (KETrack &tre, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, KTrack *ref=0, TrackMatrix *prop_matrix=0) const
	Propagate with error, but without noise. More...
std::optional< double >	noise_prop (KETrack &tre, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, KTrack *ref=0) const
	Propagate with error and noise. More...
std::optional< double >	dedx_prop (double pinv, double mass, double s, double *deriv=0) const
	Method to calculate updated momentum due to dE/dx. More...

Private Attributes
detinfo::DetectorPropertiesData const &	fDetProp
double	fTcut
	Maximum delta ray energy for dE/dx. More...
bool	fDoDedx
	Energy loss enable flag. More...
std::shared_ptr< const Interactor >	fInteractor
	Interactor (for calculating noise). More...

Detailed Description

Definition at line 91 of file Propagator.h.

Member Enumeration Documentation

enum trkf::Propagator::PropDirection

Propagation direction enum.

Enumerator
FORWARD
BACKWARD
UNKNOWN

Definition at line 94 of file Propagator.h.

{ FORWARD, BACKWARD, UNKNOWN };

Constructor & Destructor Documentation

trkf::Propagator::Propagator	(	detinfo::DetectorPropertiesData const &	detProp,
		double	tcut,
		bool	doDedx,
		const std::shared_ptr< const Interactor > &	interactor
	)

Constructor.

Constructor.

Arguments:

tcut - Maximum delta ray energy. doDedx - dE/dx enable flag.

Definition at line 26 of file Propagator.cxx.

References fDoDedx, fInteractor, fTcut, and ~Propagator().

    : fDetProp{detProp}, fTcut(tcut), fDoDedx(doDedx), fInteractor(interactor)
  {}

trkf::Propagator::~Propagator ( )

virtualdefault

Destructor.

Referenced by Propagator().

Member Function Documentation

virtual Propagator* trkf::Propagator::clone ( ) const

pure virtual

Clone method.

Implemented in trkf::PropAny, trkf::PropXYZPlane, trkf::PropYZPlane, and trkf::PropYZLine.

std::optional< double > trkf::Propagator::dedx_prop	(	double	pinv,
		double	mass,
		double	s,
		double *	deriv = 0
	)		const

Method to calculate updated momentum due to dE/dx.

Method to calculate updated momentum due to dE/dx.

Arguments:

pinv - Initial inverse momentum (units c/GeV). mass - Particle mass (GeV/c^2). s - Path distance. deriv - Pointer to store derivative d(pinv2)/d(pinv1) if nonzero.

Returns: Final inverse momentum (pinv2) + success flag.

Failure is returned in case of range out.

Inverse momentum can be signed (q/p). Returned inverse momentum has the same sign as the input.

In this method, we are solving the differential equation in terms of energy.

dE/dx = -f(E)

where f(E) is the stopping power returned by method LArProperties::Eloss.

We expect that this method will be called exclusively for short distance propagation. The differential equation is solved using the midpoint method using a single step, which requires two evaluations of f(E).

dE = -s*f(E1) E2 = E1 - s*f(E1 + 0.5*dE)

The derivative is calculated assuming E2 = E1 + constant, giving

d(pinv2)/d(pinv1) = pinv2^3 E2 / (pinv1^3 E1).

Definition at line 447 of file Propagator.cxx.

References util::abs(), detinfo::DetectorPropertiesData::Eloss(), fDetProp, and fTcut.

Referenced by trkf::PropYZLine::short_vec_prop(), trkf::PropYZPlane::short_vec_prop(), and trkf::PropXYZPlane::short_vec_prop().

  {
    // For infinite initial momentum, return with success status,
    // still infinite momentum.

    if (pinv == 0.) return std::make_optional(0.);

    // Set the default return value to be uninitialized with value 0.

    std::optional<double> result{std::nullopt};

    // Calculate final energy.

    double p1 = 1. / std::abs(pinv);
    double e1 = std::hypot(p1, mass);
    double de = -0.001 * s * fDetProp.Eloss(p1, mass, fTcut);
    double emid = e1 + 0.5 * de;
    if (emid > mass) {
      double pmid = std::sqrt(emid * emid - mass * mass);
      double e2 = e1 - 0.001 * s * fDetProp.Eloss(pmid, mass, fTcut);
      if (e2 > mass) {
        double p2 = std::sqrt(e2 * e2 - mass * mass);
        double pinv2 = 1. / p2;
        if (pinv < 0.) pinv2 = -pinv2;

        // Calculation was successful, update result.

        result = std::make_optional(pinv2);

        // Also calculate derivative, if requested.

        if (deriv != 0) *deriv = pinv2 * pinv2 * pinv2 * e2 / (pinv * pinv * pinv * e1);
      }
    }

    // Done.

    return result;
  }

std::optional< double > trkf::Propagator::err_prop	(	KETrack &	tre,
		const std::shared_ptr< const Surface > &	psurf,
		PropDirection	dir,
		bool	doDedx,
		KTrack *	ref = 0,
		TrackMatrix *	prop_matrix = 0
	)		const

Propagate with error, but without noise.

Propagate with error, but without noise (i.e. reversibly).

Arguments:

tre - Track to propagate. psurf - Destination surface. dir - Propagation direction (FORWARD, BACKWARD, or UNKNOWN). doDedx - dE/dx enable/disable flag. ref - Reference track (for linearized propagation). Can be null. prop_matrix - Return propagation matrix if not null.

Returned value: propagation distance + success flag.

Definition at line 343 of file Propagator.cxx.

References trkf::KETrack::getError(), lin_prop(), and trkf::KETrack::setError().

Referenced by trkf::KGTrack::fillTrack(), trkf::KalmanFilterAlg::fitMomentumMS(), and trkf::KHit< N >::predict().

  {
    // Propagate without error, get propagation matrix.

    TrackMatrix prop_temp;
    if (prop_matrix == 0) prop_matrix = &prop_temp;
    std::optional<double> result = lin_prop(tre, psurf, dir, doDedx, ref, prop_matrix, 0);

    // If propagation succeeded, update track error matrix.

    if (!!result) {
      TrackMatrix temp = prod(tre.getError(), trans(*prop_matrix));
      TrackMatrix temp2 = prod(*prop_matrix, temp);
      TrackError newerr = ublas::symmetric_adaptor<TrackMatrix>(temp2);
      tre.setError(newerr);
    }

    // Done.

    return result;
  }

bool trkf::Propagator::getDoDedx ( ) const

inline

Definition at line 108 of file Propagator.h.

Referenced by trkf::PropYZLine::short_vec_prop(), trkf::PropYZPlane::short_vec_prop(), trkf::PropXYZPlane::short_vec_prop(), and vec_prop().

{ return fDoDedx; }

const std::shared_ptr<const Interactor>& trkf::Propagator::getInteractor ( ) const

inline

Definition at line 109 of file Propagator.h.

References dir.

Referenced by trkf::PropYZLine::short_vec_prop(), trkf::PropYZPlane::short_vec_prop(), and trkf::PropXYZPlane::short_vec_prop().

{ return fInteractor; }

double trkf::Propagator::getTcut ( ) const

inline

Definition at line 107 of file Propagator.h.

{ return fTcut; }

std::optional< double > trkf::Propagator::lin_prop	(	KTrack &	trk,
		const std::shared_ptr< const Surface > &	psurf,
		PropDirection	dir,
		bool	doDedx,
		KTrack *	ref = 0,
		TrackMatrix *	prop_matrix = 0,
		TrackError *	noise_matrix = 0
	)		const

Linearized propagate without error.

Linearized propagate without error.

Arguments:

trk - Track to propagate. psurf - Destination surface. dir - Propagation direction (FORWARD, BACKWARD, or UNKNOWN). doDedx - dE/dx enable/disable flag. ref - Reference track (for linearized propagation). Can be null. prop_matrix - Return propagation matrix if not null. noise_matrix - Return noise matrix if not null.

Returned value: Propagation distance & success flag.

If the reference track is null, this method simply calls vec_prop.

Definition at line 249 of file Propagator.cxx.

References trkf::KTrack::getDirection(), trkf::KTrack::getSurface(), trkf::KTrack::getVector(), trkf::KTrack::isValid(), trkf::KTrack::setDirection(), trkf::KTrack::setSurface(), trkf::KTrack::setVector(), and vec_prop().

Referenced by err_prop(), and noise_prop().

  {
    // Default result.

    std::optional<double> result;

    if (ref == 0)
      result = vec_prop(trk, psurf, dir, doDedx, prop_matrix, noise_matrix);
    else {

      // A reference track has been provided.

      // It is an error (throw exception) if the reference track and
      // the track to be propagted are not on the same surface.

      if (!trk.getSurface()->isEqual(*(ref->getSurface())))
        throw cet::exception("Propagator")
          << "Input track and reference track not on same surface.\n";

      // Remember the starting track and reference track.

      KTrack trk0(trk);
      KTrack ref0(*ref);

      // Propagate the reference track.  Make sure we calculate the
      // propagation matrix.

      TrackMatrix prop_temp;
      if (prop_matrix == 0) prop_matrix = &prop_temp;

      // Do the propgation.  The returned result will be the result of
      // this propagatrion.

      result = vec_prop(*ref, psurf, dir, doDedx, prop_matrix, noise_matrix);
      if (!!result) {

        // Propagation of reference track succeeded.  Update the track
        // state vector and surface of the track to be propagated.

        TrackVector diff = trk.getSurface()->getDiff(trk.getVector(), ref0.getVector());
        TrackVector newvec = ref->getVector() + prod(*prop_matrix, diff);

        // Store updated state vector and surface.

        trk.setVector(newvec);
        trk.setSurface(psurf);
        trk.setDirection(ref->getDirection());
        if (!trk.getSurface()->isEqual(*(ref->getSurface())))
          throw cet::exception("Propagator") << __func__ << ": surface mismatch";

        // Final validity check.  In case of failure, restore the track
        // and reference track to their starting values.

        if (!trk.isValid()) {
          result = std::nullopt;
          trk = trk0;
          *ref = ref0;
        }
      }
      else {

        // Propagation failed.
        // Restore the reference track to its starting value, so that we ensure
        // the reference track and the actual track remain on the same surface.

        trk = trk0;
        *ref = ref0;
      }
    }

    // Done.

    return result;
  }

std::optional< double > trkf::Propagator::noise_prop	(	KETrack &	tre,
		const std::shared_ptr< const Surface > &	psurf,
		PropDirection	dir,
		bool	doDedx,
		KTrack *	ref = 0
	)		const

Propagate with error and noise.

Propagate with error and noise.

Arguments:

tre - Track to propagate. psurf - Destination surface. dir - Propagation direction (FORWARD, BACKWARD, or UNKNOWN). doDedx - dE/dx enable/disable flag. ref - Reference track (for linearized propagation). Can be null.

Returned value: propagation distance + success flag.

Definition at line 382 of file Propagator.cxx.

References trkf::KETrack::getError(), lin_prop(), and trkf::KETrack::setError().

Referenced by trkf::KalmanFilterAlg::buildTrack(), trkf::KalmanFilterAlg::extendTrack(), trkf::KalmanFilterAlg::fitMomentumMS(), trkf::KalmanFilterAlg::smoothTrack(), and trkf::KalmanFilterAlg::updateMomentum().

  {
    // Propagate without error, get propagation matrix and noise matrix.

    TrackMatrix prop_matrix;
    TrackError noise_matrix;
    std::optional<double> result =
      lin_prop(tre, psurf, dir, doDedx, ref, &prop_matrix, &noise_matrix);

    // If propagation succeeded, update track error matrix.

    if (!!result) {
      TrackMatrix temp = prod(tre.getError(), trans(prop_matrix));
      TrackMatrix temp2 = prod(prop_matrix, temp);
      TrackError newerr = ublas::symmetric_adaptor<TrackMatrix>(temp2);
      newerr += noise_matrix;
      tre.setError(newerr);
    }

    // Done.

    return result;
  }

virtual std::optional<double> trkf::Propagator::origin_vec_prop ( KTrack &  trk, const std::shared_ptr< const Surface > &  porient, TrackMatrix *  prop_matrix = 0  ) const

pure virtual

Propagate without error to surface whose origin parameters coincide with track position.

Implemented in trkf::PropAny, trkf::PropXYZPlane, trkf::PropYZPlane, and trkf::PropYZLine.

virtual std::optional<double> trkf::Propagator::short_vec_prop ( KTrack &  trk, const std::shared_ptr< const Surface > &  psurf, PropDirection  dir, bool  doDedx, TrackMatrix *  prop_matrix = 0, TrackError *  noise_matrix = 0  ) const

pure virtual

Propagate without error (short distance).

Implemented in trkf::PropAny, trkf::PropXYZPlane, trkf::PropYZPlane, and trkf::PropYZLine.

Referenced by vec_prop().

std::optional< double > trkf::Propagator::vec_prop	(	KTrack &	trk,
		const std::shared_ptr< const Surface > &	psurf,
		PropDirection	dir,
		bool	doDedx,
		TrackMatrix *	prop_matrix = 0,
		TrackError *	noise_matrix = 0
	)		const

Propagate without error (long distance).

Propagate without error (long distance).

Arguments:

trk - Track to propagate. psurf - Destination surface. dir - Propagation direction (FORWARD, BACKWARD, or UNKNOWN). doDedx - dE/dx enable/disable flag. prop_matrix - Return propagation matrix if not null. noise_matrix - Return noise matrix if not null.

Returned value: Propagation distance & success flag.

This method calls virtual method short_vec_prop in steps of some maximum size.

Definition at line 52 of file Propagator.cxx.

References util::abs(), larg4::dist(), e, detinfo::DetectorPropertiesData::Eloss(), fDetProp, fTcut, getDoDedx(), trkf::KTrack::getMomentum(), trkf::KTrack::getPosition(), trkf::KTrack::getVector(), trkf::KTrack::Mass(), and short_vec_prop().

Referenced by trkf::KalmanFilterAlg::fitMomentumMS(), lin_prop(), trkf::KHitContainer::sort(), and trkf::KalmanFilterAlg::updateMomentum().

  {
    std::optional<double> result{std::nullopt};

    // Get the inverse momentum (assumed to be track parameter four).

    double pinv = trk.getVector()(4);

    // If dE/dx is not requested, or if inverse momentum is zero, then
    // it is safe to propagate in one step.  In this case, just pass
    // the call to short_vec_prop with unlimited distance.

    bool dedx = getDoDedx() && doDedx;
    if (!dedx || pinv == 0.)
      result = short_vec_prop(trk, psurf, dir, dedx, prop_matrix, noise_matrix);

    else {

      // dE/dx is requested.  In this case we limit the maximum
      // propagation distance such that the kinetic energy of the
      // particle should not change by more thatn 10%.

      // Initialize propagation matrix to unit matrix (if specified).

      int nvec = trk.getVector().size();
      if (prop_matrix) *prop_matrix = ublas::identity_matrix<TrackVector::value_type>(nvec);

      // Initialize noise matrix to zero matrix (if specified).

      if (noise_matrix) {
        noise_matrix->resize(nvec, nvec, false);
        noise_matrix->clear();
      }

      // Remember the starting track.

      KTrack trk0(trk);

      // Make pointer variables pointing to local versions of the
      // propagation and noise matrices, or null if not specified.

      TrackMatrix local_prop_matrix;
      TrackMatrix* plocal_prop_matrix = (prop_matrix == 0 ? 0 : &local_prop_matrix);
      TrackError local_noise_matrix;
      TrackError* plocal_noise_matrix = (noise_matrix == 0 ? 0 : &local_noise_matrix);

      // Cumulative propagation distance.

      double s = 0.;

      // Begin stepping loop.
      // We put a maximum iteration count to prevent infinite loops caused by
      // floating point pathologies.  The iteration count is large enough to reach
      // any point in the tpc using the minimum step size (for a reasonable tpc).

      bool done = false;
      int nitmax = 10000; // Maximum number of iterations.
      int nit = 0;        // Iteration count.
      while (!done) {

        // If the iteration count exceeds the maximum, return failure.

        ++nit;
        if (nit > nitmax) {
          trk = trk0;
          result = std::nullopt;
          return result;
        }

        // Estimate maximum step distance according to the above
        // stated principle.

        pinv = trk.getVector()(4);
        double mass = trk.Mass();
        double p = 1. / std::abs(pinv);
        double e = std::hypot(p, mass);
        double t = p * p / (e + mass);
        double dedx = 0.001 * fDetProp.Eloss(p, mass, fTcut);
        double smax = 0.1 * t / dedx;
        if (smax <= 0.)
          throw cet::exception("Propagator") << __func__ << ": maximum step " << smax << "\n";

        // Always allow a step of at least 0.3 cm (about one wire spacing).

        if (smax < 0.3) smax = 0.3;

        // First do a test propagation (without dE/dx and errors) to
        // find the distance to the destination surface.

        KTrack trktest(trk);
        std::optional<double> dist = short_vec_prop(trktest, psurf, dir, false, 0, 0);

        // If the test propagation failed, return failure.

        if (!dist) {
          trk = trk0;
          return dist;
        }

        // Generate destionation surface for this step (either final
        // destination, or some intermediate surface).

        std::shared_ptr<const Surface> pstep;
        if (std::abs(*dist) <= smax) {
          done = true;
          pstep = psurf;
        }
        else {

          // Generate intermediate surface.
          // First get point where track will intersect intermediate surface.

          double xyz0[3]; // Starting point.
          trk.getPosition(xyz0);
          double xyz1[3]; // Destination point.
          trktest.getPosition(xyz1);
          double frac = smax / std::abs(*dist);
          double xyz[3]; // Intermediate point.
          xyz[0] = xyz0[0] + frac * (xyz1[0] - xyz0[0]);
          xyz[1] = xyz0[1] + frac * (xyz1[1] - xyz0[1]);
          xyz[2] = xyz0[2] + frac * (xyz1[2] - xyz0[2]);

          // Choose orientation of intermediate surface perpendicular
          // to track.

          double mom[3];
          trk.getMomentum(mom);

          // Make intermediate surface object.

          pstep = std::shared_ptr<const Surface>(
            new SurfXYZPlane(xyz[0], xyz[1], xyz[2], mom[0], mom[1], mom[2]));
        }

        // Do the actual step propagation.

        dist = short_vec_prop(trk, pstep, dir, doDedx, plocal_prop_matrix, plocal_noise_matrix);

        // If the step propagation failed, return failure.

        if (!dist) {
          trk = trk0;
          return dist;
        }

        // Update cumulative propagation distance.

        s += *dist;

        // Update cumulative propagation matrix (left-multiply).

        if (prop_matrix != 0) {
          TrackMatrix temp = prod(*plocal_prop_matrix, *prop_matrix);
          *prop_matrix = temp;
        }

        // Update cumulative noise matrix.

        if (noise_matrix != 0) {
          TrackMatrix temp = prod(*noise_matrix, trans(*plocal_prop_matrix));
          TrackMatrix temp2 = prod(*plocal_prop_matrix, temp);
          *noise_matrix = ublas::symmetric_adaptor<TrackMatrix>(temp2);
          *noise_matrix += *plocal_noise_matrix;
        }
      }

      // Set the final result (distance + success).

      result = std::make_optional(s);
    }

    // Done.

    return result;
  }

Member Data Documentation

detinfo::DetectorPropertiesData const& trkf::Propagator::fDetProp

private

Definition at line 164 of file Propagator.h.

Referenced by dedx_prop(), and vec_prop().

bool trkf::Propagator::fDoDedx

private

Energy loss enable flag.

Definition at line 166 of file Propagator.h.

Referenced by Propagator().

std::shared_ptr<const Interactor> trkf::Propagator::fInteractor

private

Interactor (for calculating noise).

Definition at line 167 of file Propagator.h.

Referenced by Propagator().

double trkf::Propagator::fTcut

private

Maximum delta ray energy for dE/dx.

Definition at line 165 of file Propagator.h.

Referenced by dedx_prop(), Propagator(), and vec_prop().

---

The documentation for this class was generated from the following files:

• lardata/v09_16_03/source/lardata/RecoObjects/Propagator.h
• lardata/v09_16_03/source/lardata/RecoObjects/Propagator.cxx