trkf::PropAny Class Reference

#include "PropAny.h"

Inheritance diagram for trkf::PropAny:

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

Public Member Functions
	PropAny (detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx)
	Constructor. More...
Propagator *	clone () const override
	Clone method. More...
std::optional< double >	short_vec_prop (KTrack &trk, const std::shared_ptr< const Surface > &surf, Propagator::PropDirection dir, bool doDedx, TrackMatrix *prop_matrix=0, TrackError *noise_matrix=0) const override
	Propagate without error. More...
virtual std::optional< double >	origin_vec_prop (KTrack &trk, const std::shared_ptr< const Surface > &porient, TrackMatrix *prop_matrix=0) const override
	Propagate without error to surface whose origin parameters coincide with track position. More...
double	getTcut () const
bool	getDoDedx () const
const std::shared_ptr< const Interactor > &	getInteractor () const
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
PropYZLine	fPropYZLine
	Underlying propagators. More...
PropYZPlane	fPropYZPlane
PropXYZPlane	fPropXYZPlane

Detailed Description

Definition at line 25 of file PropAny.h.

Member Enumeration Documentation

enum trkf::Propagator::PropDirection

inherited

Propagation direction enum.

Enumerator
FORWARD
BACKWARD
UNKNOWN

Definition at line 94 of file Propagator.h.

{ FORWARD, BACKWARD, UNKNOWN };

Constructor & Destructor Documentation

trkf::PropAny::PropAny	(	detinfo::DetectorPropertiesData const &	detProp,
		double	tcut,
		bool	doDedx
	)

Constructor.

Constructor.

Arguments.

tcut - Delta ray energy cutoff for calculating dE/dx. doDedx - dE/dx enable flag.

Definition at line 27 of file PropAny.cxx.

References fPropXYZPlane, fPropYZLine, and fPropYZPlane.

Referenced by clone().

    : Propagator(detProp,
                 tcut,
                 doDedx,
                 (tcut >= 0 ? std::make_shared<InteractPlane const>(detProp, tcut) :
                              std::shared_ptr<Interactor const>{}))
    , fPropYZLine(detProp, tcut, doDedx)
    , fPropYZPlane(detProp, tcut, doDedx)
    , fPropXYZPlane(detProp, tcut, doDedx)
  {}

Member Function Documentation

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

inlineoverridevirtual

Clone method.

Implements trkf::Propagator.

Definition at line 30 of file PropAny.h.

References dir, origin_vec_prop(), PropAny(), and short_vec_prop().

{ return new PropAny(*this); }

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

inherited

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(), trkf::Propagator::fDetProp, and trkf::Propagator::fTcut.

Referenced by trkf::PropYZLine::short_vec_prop(), trkf::PropXYZPlane::short_vec_prop(), and trkf::PropYZPlane::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

inherited

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(), trkf::Propagator::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

inlineinherited

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 trkf::Propagator::vec_prop().

{ return fDoDedx; }

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

inlineinherited

Definition at line 109 of file Propagator.h.

References dir.

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

{ return fInteractor; }

double trkf::Propagator::getTcut ( ) const

inlineinherited

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

inherited

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 trkf::Propagator::vec_prop().

Referenced by trkf::Propagator::err_prop(), and trkf::Propagator::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

inherited

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(), trkf::Propagator::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;
  }

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

overridevirtual

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

Propagate without error to dynamically generated origin surface. Optionally return propagation matrix.

Arguments:

trk - Track to propagate. porient - Orientation surface. prop_matrix - Pointer to optional propagation matrix.

Returned value: propagation distance + success flag.

Propagation distance is always zero after successful propagation.

Implements trkf::Propagator.

Definition at line 94 of file PropAny.cxx.

References fPropXYZPlane, fPropYZLine, fPropYZPlane, trkf::PropYZLine::origin_vec_prop(), trkf::PropXYZPlane::origin_vec_prop(), and trkf::PropYZPlane::origin_vec_prop().

Referenced by clone().

  {
    // Default result.

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

    // Test the type of the destination surface.

    if (dynamic_cast<const SurfYZLine*>(&*porient))
      result = fPropYZLine.origin_vec_prop(trk, porient, prop_matrix);
    else if (dynamic_cast<const SurfYZPlane*>(&*porient))
      result = fPropYZPlane.origin_vec_prop(trk, porient, prop_matrix);
    else if (dynamic_cast<const SurfXYZPlane*>(&*porient))
      result = fPropXYZPlane.origin_vec_prop(trk, porient, prop_matrix);
    else
      throw cet::exception("PropAny") << "Destination surface has unknown type.\n";

    // Done.

    return result;
  }

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

overridevirtual

Propagate without error.

Propagate without error. Optionally return propagation matrix and noise matrix. This method tests the type of the destination surface, and calls the corresponding typed propagator.

Arguments:

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

Returned value: propagation distance + success flag.

Implements trkf::Propagator.

Definition at line 54 of file PropAny.cxx.

References fPropXYZPlane, fPropYZLine, fPropYZPlane, trkf::PropYZLine::short_vec_prop(), trkf::PropXYZPlane::short_vec_prop(), and trkf::PropYZPlane::short_vec_prop().

Referenced by clone(), and trkf::InteractGeneral::noise().

  {
    // Default result.

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

    // Test the type of the destination surface.

    if (dynamic_cast<const SurfYZLine*>(&*psurf))
      result = fPropYZLine.short_vec_prop(trk, psurf, dir, doDedx, prop_matrix, noise_matrix);
    else if (dynamic_cast<const SurfYZPlane*>(&*psurf))
      result = fPropYZPlane.short_vec_prop(trk, psurf, dir, doDedx, prop_matrix, noise_matrix);
    else if (dynamic_cast<const SurfXYZPlane*>(&*psurf))
      result = fPropXYZPlane.short_vec_prop(trk, psurf, dir, doDedx, prop_matrix, noise_matrix);
    else
      throw cet::exception("PropAny") << "Destination surface has unknown type.\n";

    // Done.

    return result;
  }

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

inherited

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(), trkf::Propagator::fDetProp, trkf::Propagator::fTcut, trkf::Propagator::getDoDedx(), trkf::KTrack::getMomentum(), trkf::KTrack::getPosition(), trkf::KTrack::getVector(), trkf::KTrack::Mass(), and trkf::Propagator::short_vec_prop().

Referenced by trkf::KalmanFilterAlg::fitMomentumMS(), trkf::Propagator::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

PropXYZPlane trkf::PropAny::fPropXYZPlane

private

Definition at line 50 of file PropAny.h.

Referenced by origin_vec_prop(), PropAny(), and short_vec_prop().

PropYZLine trkf::PropAny::fPropYZLine

private

Underlying propagators.

Definition at line 48 of file PropAny.h.

Referenced by origin_vec_prop(), PropAny(), and short_vec_prop().

PropYZPlane trkf::PropAny::fPropYZPlane

private

Definition at line 49 of file PropAny.h.

Referenced by origin_vec_prop(), PropAny(), and short_vec_prop().

---

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

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