trkf::PropYZPlane Class Reference

#include "PropYZPlane.h"

Inheritance diagram for trkf::PropYZPlane:

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

Public Member Functions
	PropYZPlane (detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx)
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...
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 Member Functions
bool	transformYZLine (double phi1, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
	Transform yz line -> yz plane. More...
bool	transformYZPlane (double phi1, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
	Transform yz plane -> yz plane. More...
bool	transformXYZPlane (double theta1, double phi1, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
	Transform xyz plane -> yz plane. More...

Detailed Description

Definition at line 24 of file PropYZPlane.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::PropYZPlane::PropYZPlane	(	detinfo::DetectorPropertiesData const &	detProp,
		double	tcut,
		bool	doDedx
	)

Constructor.

Arguments.

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

Definition at line 28 of file PropYZPlane.cxx.

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

Member Function Documentation

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

inlineoverridevirtual

Clone method.

Implements trkf::Propagator.

Definition at line 29 of file PropYZPlane.h.

References dir.

{ return new PropYZPlane(*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(), 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

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(), 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(), short_vec_prop(), and trkf::PropXYZPlane::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::PropYZPlane::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 257 of file PropYZPlane.cxx.

References dir, trkf::KTrack::getDirection(), trkf::KTrack::getPosition(), trkf::KTrack::getSurface(), trkf::KTrack::getVector(), trkf::KTrack::isValid(), trkf::SurfYZPlane::phi(), trkf::KTrack::setDirection(), trkf::KTrack::setSurface(), trkf::KTrack::setVector(), transformXYZPlane(), transformYZLine(), and transformYZPlane().

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

  {
    // Set the default return value to be unitialized with value 0.

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

    // Remember starting track.

    KTrack trk0(trk);

    // Get initial track parameters and direction.
    // Note the initial track can be on any type of surface.

    TrackVector vec = trk.getVector(); // Modifiable copy.
    if (vec.size() != 5)
      throw cet::exception("PropYZPlane")
        << "Track state vector has wrong size" << vec.size() << "\n";
    Surface::TrackDirection dir = trk.getDirection();

    // Get track position.

    double xyz[3];
    trk.getPosition(xyz);
    double x02 = xyz[0];
    double y02 = xyz[1];
    double z02 = xyz[2];

    // Generate the origin surface, which will be the destination surface.
    // Return failure if orientation surface is the wrong type.

    const SurfYZPlane* orient = dynamic_cast<const SurfYZPlane*>(&*porient);
    if (orient == 0) return result;
    double phi2 = orient->phi();
    std::shared_ptr<const Surface> porigin(new SurfYZPlane(x02, y02, z02, phi2));

    // Test initial surface types.

    if (const SurfYZLine* from = dynamic_cast<const SurfYZLine*>(&*trk.getSurface())) {

      // Initial surface is SurfYZLine.
      // Get surface paramters.

      double phi1 = from->phi();

      // Transform track to origin surface.

      bool ok = transformYZLine(phi1, phi2, vec, dir, prop_matrix);
      result = std::make_optional(0.);
      if (!ok) return std::nullopt;
    }
    else if (const SurfYZPlane* from = dynamic_cast<const SurfYZPlane*>(&*trk.getSurface())) {

      // Initial surface is SurfYZPlane.
      // Get surface paramters.

      double phi1 = from->phi();

      // Transform track to origin surface.

      bool ok = transformYZPlane(phi1, phi2, vec, dir, prop_matrix);
      result = std::make_optional(0.);
      if (!ok) return std::nullopt;
    }
    else if (const SurfXYZPlane* from = dynamic_cast<const SurfXYZPlane*>(&*trk.getSurface())) {

      // Initial surface is SurfXYZPlane.
      // Get surface paramters.

      double theta1 = from->theta();
      double phi1 = from->phi();

      // Transform track to origin surface.

      bool ok = transformXYZPlane(theta1, phi1, phi2, vec, dir, prop_matrix);
      result = std::make_optional(0.);
      if (!ok) return std::nullopt;
    }

    // Update track.

    trk.setSurface(porigin);
    trk.setVector(vec);
    trk.setDirection(dir);

    // Final validity check.

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

    // Done.

    return result;
  }

std::optional< double > trkf::PropYZPlane::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.

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 50 of file PropYZPlane.cxx.

References trkf::Surface::BACKWARD, trkf::Propagator::BACKWARD, trkf::Propagator::dedx_prop(), trkf::Propagator::FORWARD, trkf::KTrack::getDirection(), trkf::Propagator::getDoDedx(), trkf::Propagator::getInteractor(), trkf::KTrack::getPosition(), trkf::KTrack::getVector(), trkf::KTrack::Mass(), origin_vec_prop(), trkf::SurfYZPlane::phi(), trkf::KTrack::setSurface(), trkf::KTrack::setVector(), trkf::Surface::UNKNOWN, trkf::Propagator::UNKNOWN, trkf::SurfYZPlane::x0(), trkf::SurfYZPlane::y0(), and trkf::SurfYZPlane::z0().

Referenced by trkf::PropAny::short_vec_prop().

  {
    // Set the default return value to be unitialized with value 0.

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

    // Get destination surface and surface parameters.
    // Return failure if wrong surface type.

    const SurfYZPlane* to = dynamic_cast<const SurfYZPlane*>(&*psurf);
    if (to == 0) return result;
    double x02 = to->x0();
    double y02 = to->y0();
    double z02 = to->z0();
    double phi2 = to->phi();

    // Remember starting track.

    KTrack trk0(trk);

    // Get track position.

    double xyz[3];
    trk.getPosition(xyz);
    double x01 = xyz[0];
    double y01 = xyz[1];
    double z01 = xyz[2];

    // Propagate to origin surface.

    TrackMatrix local_prop_matrix;
    TrackMatrix* plocal_prop_matrix = (prop_matrix == 0 ? 0 : &local_prop_matrix);
    std::optional<double> result1 = origin_vec_prop(trk, psurf, plocal_prop_matrix);
    if (!result1) return result1;

    // Get the intermediate track state vector and track parameters.

    const TrackVector& vec = trk.getVector();
    if (vec.size() != 5)
      throw cet::exception("PropYZPlane")
        << "Track state vector has wrong size" << vec.size() << "\n";
    double u1 = vec(0);
    double v1 = vec(1);
    double dudw1 = vec(2);
    double dvdw1 = vec(3);
    double pinv = vec(4);
    Surface::TrackDirection dir1 = trk.getDirection();

    // Make sure intermediate track has a valid direction.

    if (dir1 == Surface::UNKNOWN) {
      trk = trk0;
      return result;
    }

    // Calculate transcendental functions.

    double sinphi2 = std::sin(phi2);
    double cosphi2 = std::cos(phi2);

    // Calculate initial position in the destination coordinate
    // system.

    double u2 = x01 - x02 + u1;
    double v2 = (y01 - y02) * cosphi2 + (z01 - z02) * sinphi2 + v1;
    double w2 = -(y01 - y02) * sinphi2 + (z01 - z02) * cosphi2;

    // Calculate position at destination surface (propagate distance -w2).

    double u2p = u2 - w2 * dudw1;
    double v2p = v2 - w2 * dvdw1;

    // Calculate the signed propagation distance.

    double s = -w2 * std::sqrt(1. + dudw1 * dudw1 + dvdw1 * dvdw1);
    if (dir1 == Surface::BACKWARD) s = -s;

    // Check if propagation was in the right direction.
    // (Compare sign of s with requested direction).

    bool sok = (dir == Propagator::UNKNOWN || (dir == Propagator::FORWARD && s >= 0.) ||
                (dir == Propagator::BACKWARD && s <= 0.));

    // If wrong direction, return failure without updating the track
    // or propagation matrix.

    if (!sok) {
      trk = trk0;
      return result;
    }

    // Find final momentum.

    double deriv = 1.;
    auto pinv2 = std::make_optional(pinv);
    if (getDoDedx() && doDedx && s != 0.) {
      double* pderiv = (prop_matrix != 0 ? &deriv : 0);
      pinv2 = dedx_prop(pinv, trk.Mass(), s, pderiv);
    }

    // Return failure in case of range out.

    if (!pinv2) {
      trk = trk0;
      return result;
    }

    // Update default result to success and store propagation distance.

    result = std::make_optional(s);

    // Update propagation matrix (if requested).

    if (prop_matrix != 0) {
      TrackMatrix pm;
      pm.resize(vec.size(), vec.size(), false);

      // Calculate partial derivatives.

      pm(0, 0) = 1.; // du2/du1
      pm(1, 0) = 0.; // dv2/du1
      pm(2, 0) = 0.; // d(dudw2)/du1
      pm(3, 0) = 0.; // d(dvdw2)/du1
      pm(4, 0) = 0.; // d(pinv2)/du1

      pm(0, 1) = 0.; // du2/dv1
      pm(1, 1) = 1.; // dv2/dv1
      pm(2, 1) = 0.; // d(dudw2)/dv1
      pm(3, 1) = 0.; // d(dvdw2)/dv1
      pm(4, 1) = 0.; // d(pinv2)/dv1

      pm(0, 2) = -w2; // du2/d(dudw1);
      pm(1, 2) = 0.;  // dv2/d(dudw1);
      pm(2, 2) = 1.;  // d(dudw2)/d(dudw1);
      pm(3, 2) = 0.;  // d(dvdw2)/d(dudw1);
      pm(4, 2) = 0.;  // d(pinv2)/d(dudw1);

      pm(0, 3) = 0.;  // du2/d(dvdw1);
      pm(1, 3) = -w2; // dv2/d(dvdw1);
      pm(2, 3) = 0.;  // d(dudw2)/d(dvdw1);
      pm(3, 3) = 1.;  // d(dvdw2)/d(dvdw1);
      pm(4, 3) = 0.;  // d(pinv2)/d(dvdw1);

      pm(0, 4) = 0.;    // du2/d(pinv1);
      pm(1, 4) = 0.;    // dv2/d(pinv1);
      pm(2, 4) = 0.;    // d(dudw2)/d(pinv1);
      pm(3, 4) = 0.;    // d(dvdw2)/d(pinv1);
      pm(4, 4) = deriv; // d(pinv2)/d(pinv1);

      // Compose the final propagation matrix from zero-distance propagation and
      // parallel surface propagation.

      *prop_matrix = prod(pm, *plocal_prop_matrix);
    }

    // Update noise matrix (if requested).

    if (noise_matrix != 0) {
      noise_matrix->resize(vec.size(), vec.size(), false);
      if (getInteractor().get() != 0) {
        bool ok = getInteractor()->noise(trk, s, *noise_matrix);
        if (!ok) {
          trk = trk0;
          return std::nullopt;
        }
      }
      else
        noise_matrix->clear();
    }

    // Construct track vector at destination surface.

    TrackVector vec2(vec.size());
    vec2(0) = u2p;
    vec2(1) = v2p;
    vec2(2) = dudw1;
    vec2(3) = dvdw1;
    vec2(4) = *pinv2;

    // Update track.

    trk.setSurface(psurf);
    trk.setVector(vec2);

    // Done.

    return result;
  }

bool trkf::PropYZPlane::transformXYZPlane ( double  theta1, double  phi1, double  phi2, TrackVector &  vec, Surface::TrackDirection &  dir, TrackMatrix *  prop_matrix  ) const

private

Transform xyz plane -> yz plane.

Definition at line 643 of file PropYZPlane.cxx.

References trkf::Surface::BACKWARD, trkf::Surface::FORWARD, and trkf::Surface::UNKNOWN.

Referenced by origin_vec_prop().

  {
    // Calculate transcendental functions.

    double sinth1 = std::sin(theta1);
    double costh1 = std::cos(theta1);

    double sindphi = std::sin(phi2 - phi1);
    double cosdphi = std::cos(phi2 - phi1);

    // Get the initial track state vector and track parameters.

    double dudw1 = vec(2);
    double dvdw1 = vec(3);

    // Make sure initial track has a valid direction.

    if (dir == Surface::UNKNOWN) return false;

    // Calculate elements of rotation matrix from initial coordinate
    // system to destination coordinte system.

    double ruu = costh1;
    double ruw = sinth1;

    double rvu = -sinth1 * sindphi;
    double rvv = cosdphi;
    double rvw = costh1 * sindphi;

    double rwu = -sinth1 * cosdphi;
    double rwv = -sindphi;
    double rww = costh1 * cosdphi;

    // Calculate the derivative dw2/dw1;
    // If dw2/dw1 == 0., that means the track is moving parallel
    // to destination plane.
    // In this case return propagation failure.

    double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
    if (dw2dw1 == 0.) return false;

    // Calculate slope in destination plane coordinates.

    double dudw2 = (dudw1 * ruu + ruw) / dw2dw1;
    double dvdw2 = (dudw1 * rvu + dvdw1 * rvv + rvw) / dw2dw1;

    // Calculate direction parameter at destination surface.
    // Direction will flip if dw2dw1 < 0.;

    switch (dir) {
    case Surface::FORWARD: dir = (dw2dw1 > 0.) ? Surface::FORWARD : Surface::BACKWARD; break;
    case Surface::BACKWARD: dir = (dw2dw1 > 0.) ? Surface::BACKWARD : Surface::FORWARD; break;
    default:
      throw cet::exception("PropXYZPlane")
        << __func__ << ": unexpected direction #" << ((int)dir) << "\n";
    } // switch

    // Update propagation matrix (if requested).

    if (prop_matrix != 0) {
      TrackMatrix& pm = *prop_matrix;
      pm.resize(vec.size(), vec.size(), false);

      // Calculate partial derivatives.

      pm(0, 0) = ruu - dudw2 * rwu; // du2/du1
      pm(1, 0) = rvu - dvdw2 * rwu; // dv2/du1
      pm(2, 0) = 0.;                // d(dudw2)/du1
      pm(3, 0) = 0.;                // d(dvdw2)/du1
      pm(4, 0) = 0.;                // d(pinv2)/du1

      pm(0, 1) = -dudw2 * rwv;      // du2/dv1
      pm(1, 1) = rvv - dvdw2 * rwv; // dv2/dv1
      pm(2, 1) = 0.;                // d(dudw2)/dv1
      pm(3, 1) = 0.;                // d(dvdw2)/dv1
      pm(4, 1) = 0.;                // d(pinv2)/dv1

      pm(0, 2) = 0.;                           // du2/d(dudw1);
      pm(1, 2) = 0.;                           // dv2/d(dudw1);
      pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1; // d(dudw2)/d(dudw1);
      pm(3, 2) = (rvu - dvdw2 * rwu) / dw2dw1; // d(dvdw2)/d(dudw1);
      pm(4, 2) = 0.;                           // d(pinv2)/d(dudw1);

      pm(0, 3) = 0.;                           // du2/d(dvdw1);
      pm(1, 3) = 0.;                           // dv2/d(dvdw1);
      pm(2, 3) = -dudw2 * rwv / dw2dw1;        // d(dudw2)/d(dvdw1);
      pm(3, 3) = (rvv - dvdw2 * rwv) / dw2dw1; // d(dvdw2)/d(dvdw1);
      pm(4, 3) = 0.;                           // d(pinv2)/d(dvdw1);

      pm(0, 4) = 0.; // du2/d(pinv1);
      pm(1, 4) = 0.; // dv2/d(pinv1);
      pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
      pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
      pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
    }

    // Update track vector.

    vec(0) = 0.;
    vec(1) = 0.;
    vec(2) = dudw2;
    vec(3) = dvdw2;

    // Done (success).

    return true;
  }

bool trkf::PropYZPlane::transformYZLine ( double  phi1, double  phi2, TrackVector &  vec, Surface::TrackDirection &  dir, TrackMatrix *  prop_matrix  ) const

private

Transform yz line -> yz plane.

The following methods transform the track parameters from initial surface to SurfYZPlane origin surface, and generate a propagation matrix. The first group of function parameters are the orientation surface parameters of the initial surface. The second group of function parameters are the orientation parameters of the of the destination surface. The origin parameters of the destination surface are assumed to match the position of the track.

Definition at line 357 of file PropYZPlane.cxx.

Referenced by origin_vec_prop().

  {
    // Calculate surface transcendental functions.

    double sindphi = std::sin(phi2 - phi1);
    double cosdphi = std::cos(phi2 - phi1);

    // Get the initial track parameters.

    double r1 = vec(0);
    double phid1 = vec(2);
    double eta1 = vec(3);

    // Calculate elements of rotation matrix from initial coordinate
    // system to destination coordinte system.

    double rvv = cosdphi;
    double rvw = sindphi;

    double rwv = -sindphi;
    double rww = cosdphi;

    // Calculate track transcendental functions.

    double sinphid1 = std::sin(phid1);
    double cosphid1 = std::cos(phid1);
    double sh1 = 1. / std::cosh(eta1); // sech(eta1)
    double th1 = std::tanh(eta1);

    // Calculate initial position in Cartesian coordinates.

    double u1 = -r1 * sinphid1;
    double w1 = r1 * cosphid1;

    // Calculate direction in destination coordinate system.

    double du2 = sh1 * cosphid1;
    double dv2 = th1 * cosdphi + sh1 * sinphid1 * sindphi;
    double dw2 = -th1 * sindphi + sh1 * sinphid1 * cosdphi;
    //double duw2 = std::hypot(du2, dw2);

    // Calculate the track direction relative to the destination surface.
    // The track direction comes from the sign of dw2 (=dw/ds).
    // If dw2 is zero, the destionation surface is unreachable, return failure.

    if (dw2 > 0.)
      dir = Surface::TrackDirection::FORWARD;
    else if (dw2 < 0.)
      dir = Surface::TrackDirection::BACKWARD;
    else
      return false;

    // Calculate final track slope track parameters.

    double dudw2 = du2 / dw2;
    double dvdw2 = dv2 / dw2;

    // Update propagation matrix (if requested).

    if (prop_matrix != 0) {
      TrackMatrix& pm = *prop_matrix;
      pm.resize(vec.size(), vec.size(), false);

      // Calculate partial derivatives.

      // Partials of initial positions and directions wrt initial t.p.'s.

      double du1dr1 = -sinphid1;
      double du1dphi1 = -w1;

      double dw1dr1 = cosphid1;
      double dw1dphi1 = u1;

      double ddu1dphi1 = -sinphid1 * sh1;
      double ddu1deta1 = -cosphid1 * sh1 * th1;

      double ddv1deta1 = sh1 * sh1;

      double ddw1dphi1 = cosphid1 * sh1;
      double ddw1deta1 = -sinphid1 * sh1 * th1;

      // Rotate partials to destination coordinate system.

      double du2dr1 = du1dr1;
      double dv2dr1 = rvw * dw1dr1;
      double dw2dr1 = rww * dw1dr1;

      double dv2dv1 = rvv;
      double dw2dv1 = rwv;

      double du2dphi1 = du1dphi1;
      double dv2dphi1 = rvw * dw1dphi1;
      double dw2dphi1 = rww * dw1dphi1;

      double ddu2dphi1 = ddu1dphi1;
      double ddv2dphi1 = rvw * ddw1dphi1;
      double ddw2dphi1 = rww * ddw1dphi1;

      double ddu2deta1 = ddu1deta1;
      double ddv2deta1 = rvv * ddv1deta1 + rvw * ddw1deta1;
      double ddw2deta1 = rwv * ddv1deta1 + rww * ddw1deta1;

      // Partials of final slope t.p. wrt final position and direction.

      double ddudw2ddu2 = 1. / dw2;
      double ddudw2ddw2 = -dudw2 / dw2;

      double ddvdw2ddv2 = 1. / dw2;
      double ddvdw2ddw2 = -dvdw2 / dw2;

      // Partials of final slope t.p. wrt initial t.p.

      double ddudw2dphi1 = ddudw2ddu2 * ddu2dphi1 + ddudw2ddw2 * ddw2dphi1;
      double ddudw2deta1 = ddudw2ddu2 * ddu2deta1 + ddudw2ddw2 * ddw2deta1;

      double ddvdw2dphi1 = ddvdw2ddv2 * ddv2dphi1 + ddvdw2ddw2 * ddw2dphi1;
      double ddvdw2deta1 = ddvdw2ddv2 * ddv2deta1 + ddvdw2ddw2 * ddw2deta1;

      // We still need to calculate the corretion due to the dependence of the
      // propagation distance on the initial track parameters.  This correction is
      // needed even though the actual propagation distance is zero.

      // This correction effects the u and v track parameters.

      // Partials of perpendicular propagation distance wrt initial t.p.

      double dstdr1 = -dw2dr1;
      double dstdv1 = -dw2dv1;
      double dstdphi1 = -dw2dphi1;

      // Calculate correction to u and v parameter partials wrt initial t.p. due to path length.

      du2dr1 += dstdr1 * dudw2;
      double du2dv1 = dstdv1 * dudw2;
      du2dphi1 += dstdphi1 * dudw2;

      dv2dr1 += dstdr1 * dvdw2;
      dv2dv1 += dstdv1 * dvdw2;
      dv2dphi1 += dstdphi1 * dvdw2;

      // Fill derivative matrix.

      pm(0, 0) = du2dr1; // du2/dr1
      pm(1, 0) = dv2dr1; // dv2/dr1
      pm(2, 0) = 0.;     // d(dudw2)/dr1
      pm(3, 0) = 0.;     // d(dvdw2)/dr1
      pm(4, 0) = 0.;     // d(pinv2)/dr1

      pm(0, 1) = du2dv1; // du2/dv1
      pm(1, 1) = dv2dv1; // dv2/dv1
      pm(2, 1) = 0.;     // d(dudw2)/dv1
      pm(3, 1) = 0.;     // d(dvdw2)/dv1
      pm(4, 1) = 0.;     // d(pinv2)/dv1

      pm(0, 2) = du2dphi1;    // du2/d(phi1);
      pm(1, 2) = dv2dphi1;    // dv2/d(phi1);
      pm(2, 2) = ddudw2dphi1; // d(dudw2)/d(phi1);
      pm(3, 2) = ddvdw2dphi1; // d(dvdw2)/d(phi1);
      pm(4, 2) = 0.;          // d(pinv2)/d(phi1);

      pm(0, 3) = 0.;          // du2/d(eta1);
      pm(1, 3) = 0.;          // dv2/d(eta1);
      pm(2, 3) = ddudw2deta1; // d(dudw2)/d(eta1);
      pm(3, 3) = ddvdw2deta1; // d(dvdw2)/d(eta1);
      pm(4, 3) = 0.;          // d(pinv2)/d(eta1);

      pm(0, 4) = 0.; // du2/d(pinv1);
      pm(1, 4) = 0.; // dv2/d(pinv1);
      pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
      pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
      pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
    }

    // Update track vector.

    vec(0) = 0.;
    vec(1) = 0.;
    vec(2) = dudw2;
    vec(3) = dvdw2;

    // Done (success).

    return true;
  }

bool trkf::PropYZPlane::transformYZPlane ( double  phi1, double  phi2, TrackVector &  vec, Surface::TrackDirection &  dir, TrackMatrix *  prop_matrix  ) const

private

Transform yz plane -> yz plane.

Definition at line 548 of file PropYZPlane.cxx.

References trkf::Surface::BACKWARD, trkf::Surface::FORWARD, and trkf::Surface::UNKNOWN.

Referenced by origin_vec_prop().

  {
    // Calculate transcendental functions.

    double sindphi = std::sin(phi2 - phi1);
    double cosdphi = std::cos(phi2 - phi1);

    // Get the initial track parameters.

    double dudw1 = vec(2);
    double dvdw1 = vec(3);

    // Make sure initial track has a valid direction.

    if (dir == Surface::UNKNOWN) return false;

    // Calculate derivative dw2/dw1.
    // If dw2/dw1 == 0., that means the track is moving parallel
    // to destination plane.
    // In this case return propagation failure.

    double dw2dw1 = cosdphi - dvdw1 * sindphi;
    if (dw2dw1 == 0.) return false;

    // Calculate slope in destrination coordiante system.

    double dudw2 = dudw1 / dw2dw1;
    double dvdw2 = (sindphi + dvdw1 * cosdphi) / dw2dw1;

    // Calculate direction parameter at destination surface.
    // Direction will flip if dw2dw1 < 0.;

    switch (dir) {
    case Surface::FORWARD: dir = (dw2dw1 > 0.) ? Surface::FORWARD : Surface::BACKWARD; break;
    case Surface::BACKWARD: dir = (dw2dw1 > 0.) ? Surface::BACKWARD : Surface::FORWARD; break;
    default: throw cet::exception("PropYZPlane") << "unexpected direction #" << ((int)dir) << "\n";
    } // switch

    // Update propagation matrix (if requested).

    if (prop_matrix != 0) {
      TrackMatrix& pm = *prop_matrix;
      pm.resize(vec.size(), vec.size(), false);

      // Calculate partial derivatives.

      pm(0, 0) = 1.; // du2/du1
      pm(1, 0) = 0.; // dv2/du1
      pm(2, 0) = 0.; // d(dudw2)/du1
      pm(3, 0) = 0.; // d(dvdw2)/du1
      pm(4, 0) = 0.; // d(pinv2)/du1

      pm(0, 1) = dudw2 * sindphi;           // du2/dv1
      pm(1, 1) = cosdphi + dvdw2 * sindphi; // dv2/dv1
      pm(2, 1) = 0.;                        // d(dudw2)/dv1
      pm(3, 1) = 0.;                        // d(dvdw2)/dv1
      pm(4, 1) = 0.;                        // d(pinv2)/dv1

      pm(0, 2) = 0.;          // du2/d(dudw1);
      pm(1, 2) = 0.;          // dv2/d(dudw1);
      pm(2, 2) = 1. / dw2dw1; // d(dudw2)/d(dudw1);
      pm(3, 2) = 0.;          // d(dvdw2)/d(dudw1);
      pm(4, 2) = 0.;          // d(pinv2)/d(dudw1);

      pm(0, 3) = 0.;                                  // du2/d(dvdw1);
      pm(1, 3) = 0.;                                  // dv2/d(dvdw1);
      pm(2, 3) = dudw1 * sindphi / (dw2dw1 * dw2dw1); // d(dudw2)/d(dvdw1);
      pm(3, 3) = 1. / (dw2dw1 * dw2dw1);              // d(dvdw2)/d(dvdw1);
      pm(4, 3) = 0.;                                  // d(pinv2)/d(dvdw1);

      pm(0, 4) = 0.; // du2/d(pinv1);
      pm(1, 4) = 0.; // dv2/d(pinv1);
      pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
      pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
      pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
    }

    // Update track vector.

    vec(0) = 0.;
    vec(1) = 0.;
    vec(2) = dudw2;
    vec(3) = dvdw2;

    // Done (success).

    return true;
  }

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;
  }

---

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

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