LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
trkf::PropYZLine Class Reference

#include "PropYZLine.h"

Inheritance diagram for trkf::PropYZLine:
trkf::Propagator

Public Types

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

Public Member Functions

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

Detailed Description

Definition at line 24 of file PropYZLine.h.

Member Enumeration Documentation

Propagation direction enum.

Enumerator
FORWARD 
BACKWARD 
UNKNOWN 

Definition at line 94 of file Propagator.h.

Constructor & Destructor Documentation

trkf::PropYZLine::PropYZLine ( 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 PropYZLine.cxx.

29  : Propagator(detProp,
30  tcut,
31  doDedx,
32  (tcut >= 0. ? std::make_shared<const InteractGeneral>(detProp, tcut) :
33  std::shared_ptr<const Interactor>{}))
34  {}
Propagator(detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx, const std::shared_ptr< const Interactor > &interactor)
Constructor.
Definition: Propagator.cxx:26

Member Function Documentation

Propagator* trkf::PropYZLine::clone ( ) const
inlineoverridevirtual

Clone method.

Implements trkf::Propagator.

Definition at line 28 of file PropYZLine.h.

References dir.

28 { return new PropYZLine(*this); }
PropYZLine(detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx)
Definition: PropYZLine.cxx:28
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 short_vec_prop(), trkf::PropXYZPlane::short_vec_prop(), and trkf::PropYZPlane::short_vec_prop().

451  {
452  // For infinite initial momentum, return with success status,
453  // still infinite momentum.
454 
455  if (pinv == 0.) return std::make_optional(0.);
456 
457  // Set the default return value to be uninitialized with value 0.
458 
459  std::optional<double> result{std::nullopt};
460 
461  // Calculate final energy.
462 
463  double p1 = 1. / std::abs(pinv);
464  double e1 = std::hypot(p1, mass);
465  double de = -0.001 * s * fDetProp.Eloss(p1, mass, fTcut);
466  double emid = e1 + 0.5 * de;
467  if (emid > mass) {
468  double pmid = std::sqrt(emid * emid - mass * mass);
469  double e2 = e1 - 0.001 * s * fDetProp.Eloss(pmid, mass, fTcut);
470  if (e2 > mass) {
471  double p2 = std::sqrt(e2 * e2 - mass * mass);
472  double pinv2 = 1. / p2;
473  if (pinv < 0.) pinv2 = -pinv2;
474 
475  // Calculation was successful, update result.
476 
477  result = std::make_optional(pinv2);
478 
479  // Also calculate derivative, if requested.
480 
481  if (deriv != 0) *deriv = pinv2 * pinv2 * pinv2 * e2 / (pinv * pinv * pinv * e1);
482  }
483  }
484 
485  // Done.
486 
487  return result;
488  }
constexpr auto abs(T v)
Returns the absolute value of the argument.
double fTcut
Maximum delta ray energy for dE/dx.
Definition: Propagator.h:165
detinfo::DetectorPropertiesData const & fDetProp
Definition: Propagator.h:164
double Eloss(double mom, double mass, double tcut) const
Restricted mean energy loss (dE/dx)
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().

349  {
350  // Propagate without error, get propagation matrix.
351 
352  TrackMatrix prop_temp;
353  if (prop_matrix == 0) prop_matrix = &prop_temp;
354  std::optional<double> result = lin_prop(tre, psurf, dir, doDedx, ref, prop_matrix, 0);
355 
356  // If propagation succeeded, update track error matrix.
357 
358  if (!!result) {
359  TrackMatrix temp = prod(tre.getError(), trans(*prop_matrix));
360  TrackMatrix temp2 = prod(*prop_matrix, temp);
361  TrackError newerr = ublas::symmetric_adaptor<TrackMatrix>(temp2);
362  tre.setError(newerr);
363  }
364 
365  // Done.
366 
367  return result;
368  }
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
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.
Definition: Propagator.cxx:249
TDirectory * dir
Definition: macro.C:5
bool trkf::Propagator::getDoDedx ( ) const
inlineinherited

Definition at line 108 of file Propagator.h.

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

108 { return fDoDedx; }
bool fDoDedx
Energy loss enable flag.
Definition: Propagator.h:166
const std::shared_ptr<const Interactor>& trkf::Propagator::getInteractor ( ) const
inlineinherited

Definition at line 109 of file Propagator.h.

References dir.

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

109 { return fInteractor; }
std::shared_ptr< const Interactor > fInteractor
Interactor (for calculating noise).
Definition: Propagator.h:167
double trkf::Propagator::getTcut ( ) const
inlineinherited

Definition at line 107 of file Propagator.h.

107 { return fTcut; }
double fTcut
Maximum delta ray energy for dE/dx.
Definition: Propagator.h:165
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().

256  {
257  // Default result.
258 
259  std::optional<double> result;
260 
261  if (ref == 0)
262  result = vec_prop(trk, psurf, dir, doDedx, prop_matrix, noise_matrix);
263  else {
264 
265  // A reference track has been provided.
266 
267  // It is an error (throw exception) if the reference track and
268  // the track to be propagted are not on the same surface.
269 
270  if (!trk.getSurface()->isEqual(*(ref->getSurface())))
271  throw cet::exception("Propagator")
272  << "Input track and reference track not on same surface.\n";
273 
274  // Remember the starting track and reference track.
275 
276  KTrack trk0(trk);
277  KTrack ref0(*ref);
278 
279  // Propagate the reference track. Make sure we calculate the
280  // propagation matrix.
281 
282  TrackMatrix prop_temp;
283  if (prop_matrix == 0) prop_matrix = &prop_temp;
284 
285  // Do the propgation. The returned result will be the result of
286  // this propagatrion.
287 
288  result = vec_prop(*ref, psurf, dir, doDedx, prop_matrix, noise_matrix);
289  if (!!result) {
290 
291  // Propagation of reference track succeeded. Update the track
292  // state vector and surface of the track to be propagated.
293 
294  TrackVector diff = trk.getSurface()->getDiff(trk.getVector(), ref0.getVector());
295  TrackVector newvec = ref->getVector() + prod(*prop_matrix, diff);
296 
297  // Store updated state vector and surface.
298 
299  trk.setVector(newvec);
300  trk.setSurface(psurf);
301  trk.setDirection(ref->getDirection());
302  if (!trk.getSurface()->isEqual(*(ref->getSurface())))
303  throw cet::exception("Propagator") << __func__ << ": surface mismatch";
304 
305  // Final validity check. In case of failure, restore the track
306  // and reference track to their starting values.
307 
308  if (!trk.isValid()) {
309  result = std::nullopt;
310  trk = trk0;
311  *ref = ref0;
312  }
313  }
314  else {
315 
316  // Propagation failed.
317  // Restore the reference track to its starting value, so that we ensure
318  // the reference track and the actual track remain on the same surface.
319 
320  trk = trk0;
321  *ref = ref0;
322  }
323  }
324 
325  // Done.
326 
327  return result;
328  }
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).
Definition: Propagator.cxx:52
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
KVector< 5 >::type TrackVector
Track state vector, dimension 5.
TDirectory * dir
Definition: macro.C:5
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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().

387  {
388  // Propagate without error, get propagation matrix and noise matrix.
389 
390  TrackMatrix prop_matrix;
391  TrackError noise_matrix;
392  std::optional<double> result =
393  lin_prop(tre, psurf, dir, doDedx, ref, &prop_matrix, &noise_matrix);
394 
395  // If propagation succeeded, update track error matrix.
396 
397  if (!!result) {
398  TrackMatrix temp = prod(tre.getError(), trans(prop_matrix));
399  TrackMatrix temp2 = prod(prop_matrix, temp);
400  TrackError newerr = ublas::symmetric_adaptor<TrackMatrix>(temp2);
401  newerr += noise_matrix;
402  tre.setError(newerr);
403  }
404 
405  // Done.
406 
407  return result;
408  }
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
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.
Definition: Propagator.cxx:249
TDirectory * dir
Definition: macro.C:5
std::optional< double > trkf::PropYZLine::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 265 of file PropYZLine.cxx.

References dir, trkf::KTrack::getDirection(), trkf::KTrack::getPosition(), trkf::KTrack::getSurface(), trkf::KTrack::getVector(), trkf::KTrack::isValid(), trkf::SurfYZLine::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().

268  {
269  // Set the default return value to be unitialized with value 0.
270 
271  std::optional<double> result{std::nullopt};
272 
273  // Remember starting track.
274 
275  KTrack trk0(trk);
276 
277  // Get initial track parameters and direction.
278  // Note the initial track can be on any type of surface.
279 
280  TrackVector vec = trk.getVector(); // Modifiable copy.
281  if (vec.size() != 5)
282  throw cet::exception("PropYZPlane")
283  << "Track state vector has wrong size" << vec.size() << "\n";
284  Surface::TrackDirection dir = trk.getDirection();
285 
286  // Get track position.
287 
288  double xyz[3];
289  trk.getPosition(xyz);
290  double x02 = xyz[0];
291  double y02 = xyz[1];
292  double z02 = xyz[2];
293 
294  // Generate the origin surface, which will be the destination surface.
295  // Return failure if orientation surface is the wrong type.
296 
297  const SurfYZLine* orient = dynamic_cast<const SurfYZLine*>(&*porient);
298  if (orient == 0) return result;
299  double phi2 = orient->phi();
300  std::shared_ptr<const Surface> porigin(new SurfYZLine(x02, y02, z02, phi2));
301 
302  // Test initial surface types.
303 
304  if (const SurfYZLine* from = dynamic_cast<const SurfYZLine*>(&*trk.getSurface())) {
305 
306  // Initial surface is SurfYZLine.
307  // Get surface paramters.
308 
309  double phi1 = from->phi();
310 
311  // Transform track to origin surface.
312 
313  bool ok = transformYZLine(phi1, phi2, vec, dir, prop_matrix);
314  result = std::make_optional(0.);
315  if (!ok) return std::nullopt;
316  }
317  else if (const SurfYZPlane* from = dynamic_cast<const SurfYZPlane*>(&*trk.getSurface())) {
318 
319  // Initial surface is SurfYZPlane.
320  // Get surface paramters.
321 
322  double phi1 = from->phi();
323 
324  // Transform track to origin surface.
325 
326  bool ok = transformYZPlane(phi1, phi2, vec, dir, prop_matrix);
327  result = std::make_optional(0.);
328  if (!ok) return std::nullopt;
329  }
330  else if (const SurfXYZPlane* from = dynamic_cast<const SurfXYZPlane*>(&*trk.getSurface())) {
331 
332  // Initial surface is SurfXYZPlane.
333  // Get surface paramters.
334 
335  double theta1 = from->theta();
336  double phi1 = from->phi();
337 
338  // Transform track to origin surface.
339 
340  bool ok = transformXYZPlane(theta1, phi1, phi2, vec, dir, prop_matrix);
341  result = std::make_optional(0.);
342  if (!ok) return std::nullopt;
343  }
344 
345  // Update track.
346 
347  trk.setSurface(porigin);
348  trk.setVector(vec);
349  trk.setDirection(dir);
350 
351  // Final validity check.
352 
353  if (!trk.isValid()) {
354  trk = trk0;
355  result = std::nullopt;
356  }
357 
358  // Done.
359 
360  return result;
361  }
TrackDirection
Track direction enum.
Definition: Surface.h:54
bool transformXYZPlane(double theta1, double phi1, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
Transform xyz plane -> yz line.
Definition: PropYZLine.cxx:728
bool transformYZLine(double phi1, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
Transform yz line -> yz line.
Definition: PropYZLine.cxx:365
bool transformYZPlane(double phi1, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
Transform yz plane -> yz line.
Definition: PropYZLine.cxx:554
KVector< 5 >::type TrackVector
Track state vector, dimension 5.
TDirectory * dir
Definition: macro.C:5
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::optional< double > trkf::PropYZLine::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 PropYZLine.cxx.

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

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

56  {
57  // Set the default return value to be unitialized with value 0.
58 
59  std::optional<double> result{std::nullopt};
60 
61  // Get destination surface and surface parameters.
62  // Return failure if wrong surface type.
63 
64  const SurfYZLine* to = dynamic_cast<const SurfYZLine*>(&*psurf);
65  if (to == 0) return result;
66  double x02 = to->x0();
67  double y02 = to->y0();
68  double z02 = to->z0();
69  double phi2 = to->phi();
70 
71  // Remember starting track.
72 
73  KTrack trk0(trk);
74 
75  // Get track position.
76 
77  double xyz[3];
78  trk.getPosition(xyz);
79  double x01 = xyz[0];
80  double y01 = xyz[1];
81  double z01 = xyz[2];
82 
83  // Propagate to origin surface.
84 
85  TrackMatrix local_prop_matrix;
86  TrackMatrix* plocal_prop_matrix = (prop_matrix == 0 ? 0 : &local_prop_matrix);
87  std::optional<double> result1 = origin_vec_prop(trk, psurf, plocal_prop_matrix);
88  if (!result1) return result1;
89 
90  // Get the intermediate track state vector and track parameters.
91 
92  const TrackVector& vec = trk.getVector();
93  if (vec.size() != 5)
94  throw cet::exception("PropYZLine")
95  << "Track state vector has wrong size" << vec.size() << "\n";
96  double r1 = vec(0);
97  double v1 = vec(1);
98  double phid1 = vec(2);
99  double eta1 = vec(3);
100  double pinv = vec(4);
101 
102  // Calculate transcendental functions.
103 
104  double sinphid1 = std::sin(phid1);
105  double cosphid1 = std::cos(phid1);
106  double sh1 = std::sinh(eta1);
107  double ch1 = std::cosh(eta1);
108  double sinphi2 = std::sin(phi2);
109  double cosphi2 = std::cos(phi2);
110 
111  // Calculate the initial position in the intermediate coordinate system.
112 
113  double u1 = -r1 * sinphid1;
114  double w1 = r1 * cosphid1;
115 
116  // Calculate initial position in the destination coordinate system.
117 
118  double u2 = x01 - x02 + u1;
119  double v2 = (y01 - y02) * cosphi2 + (z01 - z02) * sinphi2 + v1;
120  double w2 = -(y01 - y02) * sinphi2 + (z01 - z02) * cosphi2 + w1;
121 
122  // Calculate the impact parameter in the destination coordinate system.
123 
124  double r2 = w2 * cosphid1 - u2 * sinphid1;
125 
126  // Calculate the perpendicular propagation distance.
127 
128  double d2 = -(w2 * sinphid1 + u2 * cosphid1);
129 
130  // Calculate the final position in the destination coordinate system.
131 
132  //double u2p = -r2 * sinphid1;
133  double v2p = v2 + d2 * sh1;
134  //double w2p = r2 * cosphid1;
135 
136  // Calculate the signed propagation distance.
137 
138  double s = d2 * ch1;
139 
140  // Check if propagation was in the right direction.
141  // (Compare sign of s with requested direction).
142 
143  bool sok = (dir == Propagator::UNKNOWN || (dir == Propagator::FORWARD && s >= 0.) ||
144  (dir == Propagator::BACKWARD && s <= 0.));
145 
146  // If wrong direction, return failure without updating the track
147  // or propagation matrix.
148 
149  if (!sok) {
150  trk = trk0;
151  return result;
152  }
153 
154  // Find final momentum.
155 
156  double deriv = 1.;
157  auto pinv2 = std::make_optional(pinv);
158  if (getDoDedx() && doDedx && s != 0.) {
159  double* pderiv = (prop_matrix != 0 ? &deriv : 0);
160  pinv2 = dedx_prop(pinv, trk.Mass(), s, pderiv);
161  }
162 
163  // Return failure in case of range out.
164 
165  if (!pinv2) {
166  trk = trk0;
167  return result;
168  }
169 
170  // Update default result to success and store propagation distance.
171 
172  result = std::make_optional(s);
173 
174  // Update propagation matrix (if requested).
175 
176  if (prop_matrix != 0) {
177  TrackMatrix pm;
178  pm.resize(vec.size(), vec.size(), false);
179 
180  // Calculate partial derivatives.
181 
182  pm(0, 0) = 1.; // dr2/dr1
183  pm(1, 0) = 0.; // dv2/dr1
184  pm(2, 0) = 0.; // d(phi2)/dr1
185  pm(3, 0) = 0.; // d(eta2)/dr1
186  pm(4, 0) = 0.; // d(pinv2)/dr1
187 
188  pm(0, 1) = 0.; // dr2/dv1
189  pm(1, 1) = 1.; // dv2/dv1
190  pm(2, 1) = 0.; // d(phi2)/dv1
191  pm(3, 1) = 0.; // d(eta2)/dv1
192  pm(4, 1) = 0.; // d(pinv2)/dv1
193 
194  pm(0, 2) = d2; // dr2/d(phi1);
195  pm(1, 2) = -r2 * sh1; // dv2/d(phi1);
196  pm(2, 2) = 1.; // d(phi2)/d(phi1);
197  pm(3, 2) = 0.; // d(eta2)/d(phi1);
198  pm(4, 2) = 0.; // d(pinv2)/d(phi1);
199 
200  pm(0, 3) = 0.; // dr2/d(eta1);
201  pm(1, 3) = d2 * ch1; // dv2/d(eta1);
202  pm(2, 3) = 0.; // d(phi2)/d(eta1);
203  pm(3, 3) = 1.; // d(eta2)/d(eta1);
204  pm(4, 3) = 0.; // d(pinv2)/d(eta1);
205 
206  pm(0, 4) = 0.; // dr2/d(pinv1);
207  pm(1, 4) = 0.; // dv2/d(pinv1);
208  pm(2, 4) = 0.; // d(phi2)/d(pinv1);
209  pm(3, 4) = 0.; // d(eta2)/d(pinv1);
210  pm(4, 4) = deriv; // d(pinv2)/d(pinv1);
211 
212  // Compose the final propagation matrix from zero-distance propagation and
213  // parallel surface propagation.
214 
215  *prop_matrix = prod(pm, *plocal_prop_matrix);
216  }
217 
218  // Update noise matrix (if requested).
219 
220  if (noise_matrix != 0) {
221  noise_matrix->resize(vec.size(), vec.size(), false);
222  if (getInteractor().get() != 0) {
223  bool ok = getInteractor()->noise(trk, s, *noise_matrix);
224  if (!ok) {
225  trk = trk0;
226  return std::nullopt;
227  }
228  }
229  else
230  noise_matrix->clear();
231  }
232 
233  // Construct track vector at destination surface.
234 
235  TrackVector vec2(vec.size());
236  vec2(0) = r2;
237  vec2(1) = v2p;
238  vec2(2) = phid1;
239  vec2(3) = eta1;
240  vec2(4) = *pinv2;
241 
242  // Update track.
243 
244  trk.setSurface(psurf);
245  trk.setVector(vec2);
246 
247  // Done.
248 
249  return result;
250  }
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.
Definition: PropYZLine.cxx:265
const std::shared_ptr< const Interactor > & getInteractor() const
Definition: Propagator.h:109
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
KVector< 5 >::type TrackVector
Track state vector, dimension 5.
bool getDoDedx() const
Definition: Propagator.h:108
TDirectory * dir
Definition: macro.C:5
std::optional< double > dedx_prop(double pinv, double mass, double s, double *deriv=0) const
Method to calculate updated momentum due to dE/dx.
Definition: Propagator.cxx:447
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool trkf::PropYZLine::transformXYZPlane ( double  theta1,
double  phi1,
double  phi2,
TrackVector vec,
Surface::TrackDirection dir,
TrackMatrix prop_matrix 
) const
private

Transform xyz plane -> yz line.

Definition at line 728 of file PropYZLine.cxx.

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

Referenced by origin_vec_prop().

734  {
735  // Calculate surface transcendental functions.
736 
737  double sinth1 = std::sin(theta1);
738  double costh1 = std::cos(theta1);
739 
740  double sindphi = std::sin(phi2 - phi1);
741  double cosdphi = std::cos(phi2 - phi1);
742 
743  // Get the initial track parameters.
744 
745  double dudw1 = vec(2);
746  double dvdw1 = vec(3);
747 
748  // Make sure initial track has a valid direction.
749 
750  double dirf = 1.;
751  if (dir == Surface::BACKWARD)
752  dirf = -1.;
753  else if (dir != Surface::FORWARD)
754  return false;
755 
756  // Calculate elements of rotation matrix from initial coordinate
757  // system to destination coordinte system.
758 
759  double ruu = costh1;
760  double ruw = sinth1;
761 
762  double rvu = -sinth1 * sindphi;
763  double rvv = cosdphi;
764  double rvw = costh1 * sindphi;
765 
766  double rwu = -sinth1 * cosdphi;
767  double rwv = -sindphi;
768  double rww = costh1 * cosdphi;
769 
770  // Calculate direction in the starting coordinate system.
771 
772  double dw1 = dirf / std::sqrt(1. + dudw1 * dudw1 + dvdw1 * dvdw1);
773  double du1 = dudw1 * dw1;
774  double dv1 = dvdw1 * dw1;
775 
776  // Rotate direction vector into destination coordinate system.
777 
778  double du2 = ruu * du1 + ruw * dw1;
779  double dv2 = rvu * du1 + rvv * dv1 + rvw * dw1;
780  double dw2 = rwu * du1 + rwv * dv1 + rww * dw1;
781  double duw2 = std::hypot(du2, dw2);
782 
783  // Calculate final direction track parameters.
784 
785  double phid2 = atan2(dw2, du2);
786  double eta2 = std::asinh(dv2 / duw2);
787 
788  // Update propagation matrix (if requested).
789 
790  if (prop_matrix != 0) {
791  TrackMatrix& pm = *prop_matrix;
792  pm.resize(vec.size(), vec.size(), false);
793 
794  // Calculate partial derivatives.
795 
796  // Partials of initial positions and directions wrt initial t.p.'s.
797 
798  double ddu1ddudw1 = (1. + dvdw1 * dvdw1) * dw1 * dw1 * dw1;
799  double ddu1ddvdw1 = -dudw1 * dvdw1 * dw1 * dw1 * dw1;
800 
801  double ddv1ddudw1 = -dudw1 * dvdw1 * dw1 * dw1 * dw1;
802  double ddv1ddvdw1 = (1. + dudw1 * dudw1) * dw1 * dw1 * dw1;
803 
804  double ddw1ddudw1 = -dudw1 * dw1 * dw1 * dw1;
805  double ddw1ddvdw1 = -dvdw1 * dw1 * dw1 * dw1;
806 
807  // Rotate partials to destination coordinate system.
808 
809  double du2du1 = ruu;
810  double dv2du1 = rvu;
811  double dw2du1 = rwu;
812 
813  double dv2dv1 = rvv;
814  double dw2dv1 = rwv;
815 
816  double ddu2ddudw1 = ruu * ddu1ddudw1 + ruw * ddw1ddudw1;
817  double ddv2ddudw1 = rvu * ddu1ddudw1 + rvv * ddv1ddudw1 + rvw * ddw1ddudw1;
818  double ddw2ddudw1 = rwu * ddu1ddudw1 + rwv * ddv1ddudw1 + rww * ddw1ddudw1;
819 
820  double ddu2ddvdw1 = ruu * ddu1ddvdw1 + ruw * ddw1ddvdw1;
821  double ddv2ddvdw1 = rvu * ddu1ddvdw1 + rvv * ddv1ddvdw1 + rvw * ddw1ddvdw1;
822  double ddw2ddvdw1 = rwu * ddu1ddvdw1 + rwv * ddv1ddvdw1 + rww * ddw1ddvdw1;
823 
824  // Partials of final t.p. wrt final position and direction.
825 
826  double dr2du2 = -dw2 / duw2;
827  double dr2dw2 = du2 / duw2;
828 
829  double dphi2ddu2 = -dw2 / (duw2 * duw2);
830  double dphi2ddw2 = du2 / (duw2 * duw2);
831 
832  double deta2ddv2 = 1. / (duw2 * duw2);
833 
834  // Partials of final t.p. wrt initial t.p.
835 
836  double dr2du1 = dr2du2 * du2du1 + dr2dw2 * dw2du1;
837  double dr2dv1 = dr2dw2 * dw2dv1;
838 
839  double dphi2ddudw1 = dphi2ddu2 * ddu2ddudw1 + dphi2ddw2 * ddw2ddudw1;
840  double dphi2ddvdw1 = dphi2ddu2 * ddu2ddvdw1 + dphi2ddw2 * ddw2ddvdw1;
841 
842  double deta2ddudw1 = deta2ddv2 * ddv2ddudw1;
843  double deta2ddvdw1 = deta2ddv2 * ddv2ddvdw1;
844 
845  // We still need to calculate the corretion due to the dependence of the
846  // propagation distance on the initial track parameters. This correction is
847  // needed even though the actual propagation distance is zero.
848 
849  // This correction only effects the v track parameter, since the v parameter
850  // the only parameter that actually dependes on the propagation distance.
851 
852  // Partials of propagation distance wrt position and direction in the destination
853  // coordinate system.
854 
855  double dsdu2 = -du2 / (duw2 * duw2);
856  double dsdw2 = -dw2 / (duw2 * duw2);
857 
858  // Partials of propagation distance wrt initial t.p.
859 
860  double dsdu1 = dsdu2 * du2du1 + dsdw2 * dw2du1;
861  double dsdv1 = dsdw2 * dw2dv1;
862 
863  // Calculate correction to v parameter partials wrt initial t.p. due to path length.
864 
865  dv2du1 += dv2 * dsdu1;
866  dv2dv1 += dv2 * dsdv1;
867 
868  // Fill matrix.
869 
870  pm(0, 0) = dr2du1; // dr2/du1
871  pm(1, 0) = dv2du1; // dv2/du1
872  pm(2, 0) = 0.; // d(phi2)/du1
873  pm(3, 0) = 0.; // d(eta2)/du1
874  pm(4, 0) = 0.; // d(pinv2)/du1
875 
876  pm(0, 1) = dr2dv1; // dr2/dv1
877  pm(1, 1) = dv2dv1; // dv2/dv1
878  pm(2, 1) = 0.; // d(phi2)/dv1
879  pm(3, 1) = 0.; // d(eta2)/dv1
880  pm(4, 1) = 0.; // d(pinv2)/dv1
881 
882  pm(0, 2) = 0.; // dr2/d(dudw1);
883  pm(1, 2) = 0.; // dv2/d(dudw1);
884  pm(2, 2) = dphi2ddudw1; // d(dudw2)/d(dudw1);
885  pm(3, 2) = deta2ddudw1; // d(eta2)/d(dudw1);
886  pm(4, 2) = 0.; // d(pinv2)/d(dudw1);
887 
888  pm(0, 3) = 0.; // dr2/d(dvdw1);
889  pm(1, 3) = 0.; // dv2/d(dvdw1);
890  pm(2, 3) = dphi2ddvdw1; // d(phi2)/d(dvdw1);
891  pm(3, 3) = deta2ddvdw1; // d(eta2)/d(dvdw1);
892  pm(4, 3) = 0.; // d(pinv2)/d(dvdw1);
893 
894  pm(0, 4) = 0.; // dr2/d(pinv1);
895  pm(1, 4) = 0.; // dv2/d(pinv1);
896  pm(2, 4) = 0.; // d(phi2)/d(pinv1);
897  pm(3, 4) = 0.; // d(eta2)/d(pinv1);
898  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
899  }
900 
901  // Update track vector.
902 
903  vec(0) = 0.;
904  vec(1) = 0.;
905  vec(2) = phid2;
906  vec(3) = eta2;
907 
908  // Done (success).
909 
910  return true;
911  }
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
TDirectory * dir
Definition: macro.C:5
bool trkf::PropYZLine::transformYZLine ( double  phi1,
double  phi2,
TrackVector vec,
Surface::TrackDirection dir,
TrackMatrix prop_matrix 
) const
private

Transform yz line -> yz line.

The following methods transform the track parameters from initial surface to SurfYZLine 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 365 of file PropYZLine.cxx.

Referenced by origin_vec_prop().

370  {
371  // Calculate surface transcendental functions.
372 
373  double sindphi = std::sin(phi2 - phi1);
374  double cosdphi = std::cos(phi2 - phi1);
375 
376  // Get the initial track parameters.
377 
378  double r1 = vec(0);
379  double phid1 = vec(2);
380  double eta1 = vec(3);
381 
382  // Calculate elements of rotation matrix from initial coordinate
383  // system to destination coordinte system.
384 
385  double rvv = cosdphi;
386  double rvw = sindphi;
387 
388  double rwv = -sindphi;
389  double rww = cosdphi;
390 
391  // Calculate track transcendental functions.
392 
393  double sinphid1 = std::sin(phid1);
394  double cosphid1 = std::cos(phid1);
395  double sh1 = 1. / std::cosh(eta1); // sech(eta1)
396  double th1 = std::tanh(eta1);
397 
398  // Calculate initial position in Cartesian coordinates.
399 
400  double u1 = -r1 * sinphid1;
401  double w1 = r1 * cosphid1;
402 
403  // Calculate direction in destination coordinate system.
404 
405  double du2 = sh1 * cosphid1;
406  double dv2 = th1 * cosdphi + sh1 * sinphid1 * sindphi;
407  double dw2 = -th1 * sindphi + sh1 * sinphid1 * cosdphi;
408  double duw2 = std::hypot(du2, dw2);
409 
410  // Calculate final direction track parameters.
411 
412  double phid2 = atan2(dw2, du2);
413  double eta2 = std::asinh(dv2 / duw2);
414 
415  // Update propagation matrix (if requested).
416 
417  if (prop_matrix != 0) {
418  TrackMatrix& pm = *prop_matrix;
419  pm.resize(vec.size(), vec.size(), false);
420 
421  // Calculate partial derivatives.
422 
423  // Partials of initial positions and directions wrt initial t.p.'s.
424 
425  double du1dr1 = -sinphid1;
426  double du1dphi1 = -w1;
427 
428  double dw1dr1 = cosphid1;
429  double dw1dphi1 = u1;
430 
431  double ddu1dphi1 = -sinphid1 * sh1;
432  double ddu1deta1 = -cosphid1 * sh1 * th1;
433 
434  double ddv1deta1 = sh1 * sh1;
435 
436  double ddw1dphi1 = cosphid1 * sh1;
437  double ddw1deta1 = -sinphid1 * sh1 * th1;
438 
439  // Rotate partials to destination coordinate system.
440 
441  double du2dr1 = du1dr1;
442  double dv2dr1 = rvw * dw1dr1;
443  double dw2dr1 = rww * dw1dr1;
444 
445  double dv2dv1 = rvv;
446  double dw2dv1 = rwv;
447 
448  double du2dphi1 = du1dphi1;
449  double dv2dphi1 = rvw * dw1dphi1;
450  double dw2dphi1 = rww * dw1dphi1;
451 
452  double ddu2dphi1 = ddu1dphi1;
453  double ddv2dphi1 = rvw * ddw1dphi1;
454  double ddw2dphi1 = rww * ddw1dphi1;
455 
456  double ddu2deta1 = ddu1deta1;
457  double ddv2deta1 = rvv * ddv1deta1 + rvw * ddw1deta1;
458  double ddw2deta1 = rwv * ddv1deta1 + rww * ddw1deta1;
459 
460  // Partials of final t.p. wrt final position and direction.
461 
462  double dr2du2 = -dw2 / duw2;
463  double dr2dw2 = du2 / duw2;
464 
465  double dphi2ddu2 = -dw2 / (duw2 * duw2);
466  double dphi2ddw2 = du2 / (duw2 * duw2);
467 
468  double deta2ddv2 = 1. / (duw2 * duw2);
469 
470  // Partials of final t.p. wrt initial t.p.
471 
472  double dr2dr1 = dr2du2 * du2dr1 + dr2dw2 * dw2dr1;
473  double dr2dv1 = dr2dw2 * dw2dv1;
474  double dr2dphi1 = dr2du2 * du2dphi1 + dr2dw2 * dw2dphi1;
475 
476  double dphi2dphi1 = dphi2ddu2 * ddu2dphi1 + dphi2ddw2 * ddw2dphi1;
477  double dphi2deta1 = dphi2ddu2 * ddu2deta1 + dphi2ddw2 * ddw2deta1;
478 
479  double deta2dphi1 = deta2ddv2 * ddv2dphi1;
480  double deta2deta1 = deta2ddv2 * ddv2deta1;
481 
482  // We still need to calculate the corretion due to the dependence of the
483  // propagation distance on the initial track parameters. This correction is
484  // needed even though the actual propagation distance is zero.
485 
486  // This correction only effects the v track parameter, since the v parameter
487  // the only parameter that actually dependes on the propagation distance.
488 
489  // Partials of propagation distance wrt position and direction in the destination
490  // coordinate system.
491 
492  double dsdu2 = -du2 / (duw2 * duw2);
493  double dsdw2 = -dw2 / (duw2 * duw2);
494 
495  // Partials of propagation distance wrt initial t.p.
496 
497  double dsdr1 = dsdu2 * du2dr1 + dsdw2 * dw2dr1;
498  double dsdv1 = dsdw2 * dw2dv1;
499  double dsdphi1 = dsdu2 * du2dphi1 + dsdw2 * dw2dphi1;
500 
501  // Calculate correction to v parameter partials wrt initial t.p. due to path length.
502 
503  dv2dr1 += dv2 * dsdr1;
504  dv2dv1 += dv2 * dsdv1;
505  dv2dphi1 += dv2 * dsdphi1;
506 
507  // Fill derivative matrix.
508 
509  pm(0, 0) = dr2dr1; // dr2/dr1
510  pm(1, 0) = dv2dr1; // dv2/dr1
511  pm(2, 0) = 0.; // d(phi2)/dr1
512  pm(3, 0) = 0.; // d(eta2)/dr1
513  pm(4, 0) = 0.; // d(pinv2)/dr1
514 
515  pm(0, 1) = dr2dv1; // dr2/dv1
516  pm(1, 1) = dv2dv1; // dv2/dv1
517  pm(2, 1) = 0.; // d(phi2)/dv1
518  pm(3, 1) = 0.; // d(eta2)/dv1
519  pm(4, 1) = 0.; // d(pinv2)/dv1
520 
521  pm(0, 2) = dr2dphi1; // dr2/d(phi1);
522  pm(1, 2) = dv2dphi1; // dv2/d(phi1);
523  pm(2, 2) = dphi2dphi1; // d(phi2)/d(phi1);
524  pm(3, 2) = deta2dphi1; // d(eta2)/d(phi1);
525  pm(4, 2) = 0.; // d(pinv2)/d(phi1);
526 
527  pm(0, 3) = 0.; // dr2/d(eta1);
528  pm(1, 3) = 0.; // dv2/d(eta1);
529  pm(2, 3) = dphi2deta1; // d(phi2)/d(eta1);
530  pm(3, 3) = deta2deta1; // d(eta2)/d(eta1);
531  pm(4, 3) = 0.; // d(pinv2)/d(eta1);
532 
533  pm(0, 4) = 0.; // dr2/d(pinv1);
534  pm(1, 4) = 0.; // dv2/d(pinv1);
535  pm(2, 4) = 0.; // d(phi2)/d(pinv1);
536  pm(3, 4) = 0.; // d(eta2)/d(pinv1);
537  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
538  }
539 
540  // Update track vector.
541 
542  vec(0) = 0.;
543  vec(1) = 0.;
544  vec(2) = phid2;
545  vec(3) = eta2;
546 
547  // Done (success).
548 
549  return true;
550  }
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
bool trkf::PropYZLine::transformYZPlane ( double  phi1,
double  phi2,
TrackVector vec,
Surface::TrackDirection dir,
TrackMatrix prop_matrix 
) const
private

Transform yz plane -> yz line.

Definition at line 554 of file PropYZLine.cxx.

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

Referenced by origin_vec_prop().

559  {
560  // Calculate surface transcendental functions.
561 
562  double sindphi = std::sin(phi2 - phi1);
563  double cosdphi = std::cos(phi2 - phi1);
564 
565  // Get the initial track parameters.
566 
567  double dudw1 = vec(2);
568  double dvdw1 = vec(3);
569 
570  // Make sure initial track has a valid direction.
571 
572  double dirf = 1.;
573  if (dir == Surface::BACKWARD)
574  dirf = -1.;
575  else if (dir != Surface::FORWARD)
576  return false;
577 
578  // Calculate elements of rotation matrix from initial coordinate
579  // system to destination coordinte system.
580 
581  double rvv = cosdphi;
582  double rvw = sindphi;
583 
584  double rwv = -sindphi;
585  double rww = cosdphi;
586 
587  // Calculate direction in the starting coordinate system.
588 
589  double dw1 = dirf / std::sqrt(1. + dudw1 * dudw1 + dvdw1 * dvdw1);
590  double du1 = dudw1 * dw1;
591  double dv1 = dvdw1 * dw1;
592 
593  // Rotate direction vector into destination coordinate system.
594 
595  double du2 = du1;
596  double dv2 = rvv * dv1 + rvw * dw1;
597  double dw2 = rwv * dv1 + rww * dw1;
598  double duw2 = std::hypot(du2, dw2);
599 
600  // Calculate final direction track parameters.
601 
602  double phid2 = atan2(dw2, du2);
603  double eta2 = std::asinh(dv2 / duw2);
604 
605  // Update propagation matrix (if requested).
606 
607  if (prop_matrix != 0) {
608  TrackMatrix& pm = *prop_matrix;
609  pm.resize(vec.size(), vec.size(), false);
610 
611  // Calculate partial derivatives.
612 
613  // Partials of initial positions and directions wrt initial t.p.'s.
614 
615  double ddu1ddudw1 = (1. + dvdw1 * dvdw1) * dw1 * dw1 * dw1;
616  double ddu1ddvdw1 = -dudw1 * dvdw1 * dw1 * dw1 * dw1;
617 
618  double ddv1ddudw1 = -dudw1 * dvdw1 * dw1 * dw1 * dw1;
619  double ddv1ddvdw1 = (1. + dudw1 * dudw1) * dw1 * dw1 * dw1;
620 
621  double ddw1ddudw1 = -dudw1 * dw1 * dw1 * dw1;
622  double ddw1ddvdw1 = -dvdw1 * dw1 * dw1 * dw1;
623 
624  // Rotate partials to destination coordinate system.
625 
626  double dv2dv1 = rvv;
627  double dw2dv1 = rwv;
628 
629  double ddu2ddudw1 = ddu1ddudw1;
630  double ddv2ddudw1 = rvv * ddv1ddudw1 + rvw * ddw1ddudw1;
631  double ddw2ddudw1 = rwv * ddv1ddudw1 + rww * ddw1ddudw1;
632 
633  double ddu2ddvdw1 = ddu1ddvdw1;
634  double ddv2ddvdw1 = rvv * ddv1ddvdw1 + rvw * ddw1ddvdw1;
635  double ddw2ddvdw1 = rwv * ddv1ddvdw1 + rww * ddw1ddvdw1;
636 
637  // Partials of final t.p. wrt final position and direction.
638 
639  double dr2du2 = -dw2 / duw2;
640  double dr2dw2 = du2 / duw2;
641 
642  double dphi2ddu2 = -dw2 / (duw2 * duw2);
643  double dphi2ddw2 = du2 / (duw2 * duw2);
644 
645  double deta2ddv2 = 1. / (duw2 * duw2);
646 
647  // Partials of final t.p. wrt initial t.p.
648 
649  double dr2du1 = dr2du2;
650  double dr2dv1 = dr2dw2 * dw2dv1;
651 
652  double dphi2ddudw1 = dphi2ddu2 * ddu2ddudw1 + dphi2ddw2 * ddw2ddudw1;
653  double dphi2ddvdw1 = dphi2ddu2 * ddu2ddvdw1 + dphi2ddw2 * ddw2ddvdw1;
654 
655  double deta2ddudw1 = deta2ddv2 * ddv2ddudw1;
656  double deta2ddvdw1 = deta2ddv2 * ddv2ddvdw1;
657 
658  // We still need to calculate the corretion due to the dependence of the
659  // propagation distance on the initial track parameters. This correction is
660  // needed even though the actual propagation distance is zero.
661 
662  // This correction only effects the v track parameter, since the v parameter
663  // the only parameter that actually dependes on the propagation distance.
664 
665  // Partials of propagation distance wrt position and direction in the destination
666  // coordinate system.
667 
668  double dsdu2 = -du2 / (duw2 * duw2);
669  double dsdw2 = -dw2 / (duw2 * duw2);
670 
671  // Partials of propagation distance wrt initial t.p.
672 
673  double dsdu1 = dsdu2;
674  double dsdv1 = dsdw2 * dw2dv1;
675 
676  // Calculate correction to v parameter partials wrt initial t.p. due to path length.
677 
678  double dv2du1 = dv2 * dsdu1;
679  dv2dv1 += dv2 * dsdv1;
680 
681  // Fill matrix.
682 
683  pm(0, 0) = dr2du1; // dr2/du1
684  pm(1, 0) = dv2du1; // dv2/du1
685  pm(2, 0) = 0.; // d(phi2)/du1
686  pm(3, 0) = 0.; // d(eta2)/du1
687  pm(4, 0) = 0.; // d(pinv2)/du1
688 
689  pm(0, 1) = dr2dv1; // dr2/dv1
690  pm(1, 1) = dv2dv1; // dv2/dv1
691  pm(2, 1) = 0.; // d(phi2)/dv1
692  pm(3, 1) = 0.; // d(eta2)/dv1
693  pm(4, 1) = 0.; // d(pinv2)/dv1
694 
695  pm(0, 2) = 0.; // dr2/d(dudw1);
696  pm(1, 2) = 0.; // dv2/d(dudw1);
697  pm(2, 2) = dphi2ddudw1; // d(dudw2)/d(dudw1);
698  pm(3, 2) = deta2ddudw1; // d(eta2)/d(dudw1);
699  pm(4, 2) = 0.; // d(pinv2)/d(dudw1);
700 
701  pm(0, 3) = 0.; // dr2/d(dvdw1);
702  pm(1, 3) = 0.; // dv2/d(dvdw1);
703  pm(2, 3) = dphi2ddvdw1; // d(phi2)/d(dvdw1);
704  pm(3, 3) = deta2ddvdw1; // d(eta2)/d(dvdw1);
705  pm(4, 3) = 0.; // d(pinv2)/d(dvdw1);
706 
707  pm(0, 4) = 0.; // dr2/d(pinv1);
708  pm(1, 4) = 0.; // dv2/d(pinv1);
709  pm(2, 4) = 0.; // d(phi2)/d(pinv1);
710  pm(3, 4) = 0.; // d(eta2)/d(pinv1);
711  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
712  }
713 
714  // Update track vector.
715 
716  vec(0) = 0.;
717  vec(1) = 0.;
718  vec(2) = phid2;
719  vec(3) = eta2;
720 
721  // Done (success).
722 
723  return true;
724  }
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
TDirectory * dir
Definition: macro.C:5
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().

58  {
59  std::optional<double> result{std::nullopt};
60 
61  // Get the inverse momentum (assumed to be track parameter four).
62 
63  double pinv = trk.getVector()(4);
64 
65  // If dE/dx is not requested, or if inverse momentum is zero, then
66  // it is safe to propagate in one step. In this case, just pass
67  // the call to short_vec_prop with unlimited distance.
68 
69  bool dedx = getDoDedx() && doDedx;
70  if (!dedx || pinv == 0.)
71  result = short_vec_prop(trk, psurf, dir, dedx, prop_matrix, noise_matrix);
72 
73  else {
74 
75  // dE/dx is requested. In this case we limit the maximum
76  // propagation distance such that the kinetic energy of the
77  // particle should not change by more thatn 10%.
78 
79  // Initialize propagation matrix to unit matrix (if specified).
80 
81  int nvec = trk.getVector().size();
82  if (prop_matrix) *prop_matrix = ublas::identity_matrix<TrackVector::value_type>(nvec);
83 
84  // Initialize noise matrix to zero matrix (if specified).
85 
86  if (noise_matrix) {
87  noise_matrix->resize(nvec, nvec, false);
88  noise_matrix->clear();
89  }
90 
91  // Remember the starting track.
92 
93  KTrack trk0(trk);
94 
95  // Make pointer variables pointing to local versions of the
96  // propagation and noise matrices, or null if not specified.
97 
98  TrackMatrix local_prop_matrix;
99  TrackMatrix* plocal_prop_matrix = (prop_matrix == 0 ? 0 : &local_prop_matrix);
100  TrackError local_noise_matrix;
101  TrackError* plocal_noise_matrix = (noise_matrix == 0 ? 0 : &local_noise_matrix);
102 
103  // Cumulative propagation distance.
104 
105  double s = 0.;
106 
107  // Begin stepping loop.
108  // We put a maximum iteration count to prevent infinite loops caused by
109  // floating point pathologies. The iteration count is large enough to reach
110  // any point in the tpc using the minimum step size (for a reasonable tpc).
111 
112  bool done = false;
113  int nitmax = 10000; // Maximum number of iterations.
114  int nit = 0; // Iteration count.
115  while (!done) {
116 
117  // If the iteration count exceeds the maximum, return failure.
118 
119  ++nit;
120  if (nit > nitmax) {
121  trk = trk0;
122  result = std::nullopt;
123  return result;
124  }
125 
126  // Estimate maximum step distance according to the above
127  // stated principle.
128 
129  pinv = trk.getVector()(4);
130  double mass = trk.Mass();
131  double p = 1. / std::abs(pinv);
132  double e = std::hypot(p, mass);
133  double t = p * p / (e + mass);
134  double dedx = 0.001 * fDetProp.Eloss(p, mass, fTcut);
135  double smax = 0.1 * t / dedx;
136  if (smax <= 0.)
137  throw cet::exception("Propagator") << __func__ << ": maximum step " << smax << "\n";
138 
139  // Always allow a step of at least 0.3 cm (about one wire spacing).
140 
141  if (smax < 0.3) smax = 0.3;
142 
143  // First do a test propagation (without dE/dx and errors) to
144  // find the distance to the destination surface.
145 
146  KTrack trktest(trk);
147  std::optional<double> dist = short_vec_prop(trktest, psurf, dir, false, 0, 0);
148 
149  // If the test propagation failed, return failure.
150 
151  if (!dist) {
152  trk = trk0;
153  return dist;
154  }
155 
156  // Generate destionation surface for this step (either final
157  // destination, or some intermediate surface).
158 
159  std::shared_ptr<const Surface> pstep;
160  if (std::abs(*dist) <= smax) {
161  done = true;
162  pstep = psurf;
163  }
164  else {
165 
166  // Generate intermediate surface.
167  // First get point where track will intersect intermediate surface.
168 
169  double xyz0[3]; // Starting point.
170  trk.getPosition(xyz0);
171  double xyz1[3]; // Destination point.
172  trktest.getPosition(xyz1);
173  double frac = smax / std::abs(*dist);
174  double xyz[3]; // Intermediate point.
175  xyz[0] = xyz0[0] + frac * (xyz1[0] - xyz0[0]);
176  xyz[1] = xyz0[1] + frac * (xyz1[1] - xyz0[1]);
177  xyz[2] = xyz0[2] + frac * (xyz1[2] - xyz0[2]);
178 
179  // Choose orientation of intermediate surface perpendicular
180  // to track.
181 
182  double mom[3];
183  trk.getMomentum(mom);
184 
185  // Make intermediate surface object.
186 
187  pstep = std::shared_ptr<const Surface>(
188  new SurfXYZPlane(xyz[0], xyz[1], xyz[2], mom[0], mom[1], mom[2]));
189  }
190 
191  // Do the actual step propagation.
192 
193  dist = short_vec_prop(trk, pstep, dir, doDedx, plocal_prop_matrix, plocal_noise_matrix);
194 
195  // If the step propagation failed, return failure.
196 
197  if (!dist) {
198  trk = trk0;
199  return dist;
200  }
201 
202  // Update cumulative propagation distance.
203 
204  s += *dist;
205 
206  // Update cumulative propagation matrix (left-multiply).
207 
208  if (prop_matrix != 0) {
209  TrackMatrix temp = prod(*plocal_prop_matrix, *prop_matrix);
210  *prop_matrix = temp;
211  }
212 
213  // Update cumulative noise matrix.
214 
215  if (noise_matrix != 0) {
216  TrackMatrix temp = prod(*noise_matrix, trans(*plocal_prop_matrix));
217  TrackMatrix temp2 = prod(*plocal_prop_matrix, temp);
218  *noise_matrix = ublas::symmetric_adaptor<TrackMatrix>(temp2);
219  *noise_matrix += *plocal_noise_matrix;
220  }
221  }
222 
223  // Set the final result (distance + success).
224 
225  result = std::make_optional(s);
226  }
227 
228  // Done.
229 
230  return result;
231  }
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
constexpr auto abs(T v)
Returns the absolute value of the argument.
double fTcut
Maximum delta ray energy for dE/dx.
Definition: Propagator.h:165
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).
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
detinfo::DetectorPropertiesData const & fDetProp
Definition: Propagator.h:164
double Eloss(double mom, double mass, double tcut) const
Restricted mean energy loss (dE/dx)
bool getDoDedx() const
Definition: Propagator.h:108
TDirectory * dir
Definition: macro.C:5
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Float_t e
Definition: plot.C:35
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

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