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

#include "PropXYZPlane.h"

Inheritance diagram for trkf::PropXYZPlane:
trkf::Propagator

Public Types

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

Public Member Functions

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

Detailed Description

Definition at line 24 of file PropXYZPlane.h.

Member Enumeration Documentation

Propagation direction enum.

Enumerator
FORWARD 
BACKWARD 
UNKNOWN 

Definition at line 94 of file Propagator.h.

Constructor & Destructor Documentation

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

Constructor.

Constructor.

Arguments.

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

Definition at line 28 of file PropXYZPlane.cxx.

31  : Propagator{detProp,
32  tcut,
33  doDedx,
34  (tcut >= 0. ? std::make_shared<InteractPlane const>(detProp, tcut) :
35  std::shared_ptr<Interactor const>{})}
36  {}
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::PropXYZPlane::clone ( ) const
inlineoverridevirtual

Clone method.

Implements trkf::Propagator.

Definition at line 29 of file PropXYZPlane.h.

References dir.

29 { return new PropXYZPlane(*this); }
PropXYZPlane(detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx)
Constructor.
std::optional< double > trkf::Propagator::dedx_prop ( double  pinv,
double  mass,
double  s,
double *  deriv = 0 
) const
inherited

Method to calculate updated momentum due to dE/dx.

Method to calculate updated momentum due to dE/dx.

Arguments:

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

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

Failure is returned in case of range out.

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

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

dE/dx = -f(E)

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

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

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

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

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

Definition at line 447 of file Propagator.cxx.

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

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

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

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

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

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

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

Transform xyz plane -> xyz plane.

Definition at line 698 of file PropXYZPlane.cxx.

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

Referenced by origin_vec_prop().

705  {
706  // Calculate transcendental functions.
707 
708  double sinth1 = std::sin(theta1);
709  double costh1 = std::cos(theta1);
710  double sinth2 = std::sin(theta2);
711  double costh2 = std::cos(theta2);
712 
713  double sindphi = std::sin(phi2 - phi1);
714  double cosdphi = std::cos(phi2 - phi1);
715 
716  // Get the initial track state vector and track parameters.
717 
718  double dudw1 = vec(2);
719  double dvdw1 = vec(3);
720 
721  // Make sure initial track has a valid direction.
722 
723  if (dir == Surface::UNKNOWN) return false;
724 
725  // Calculate elements of rotation matrix from initial coordinate
726  // system to destination coordinte system.
727 
728  double ruu = costh1 * costh2 + sinth1 * sinth2 * cosdphi;
729  double ruv = sinth2 * sindphi;
730  double ruw = sinth1 * costh2 - costh1 * sinth2 * cosdphi;
731 
732  double rvu = -sinth1 * sindphi;
733  double rvv = cosdphi;
734  double rvw = costh1 * sindphi;
735 
736  double rwu = costh1 * sinth2 - sinth1 * costh2 * cosdphi;
737  double rwv = -costh2 * sindphi;
738  double rww = sinth1 * sinth2 + costh1 * costh2 * cosdphi;
739 
740  // Calculate the derivative dw2/dw1;
741  // If dw2/dw1 == 0., that means the track is moving parallel
742  // to destination plane.
743  // In this case return propagation failure.
744 
745  double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
746  if (dw2dw1 == 0.) return false;
747 
748  // Calculate slope in destination plane coordinates.
749 
750  double dudw2 = (dudw1 * ruu + dvdw1 * ruv + ruw) / dw2dw1;
751  double dvdw2 = (dudw1 * rvu + dvdw1 * rvv + rvw) / dw2dw1;
752 
753  // Calculate direction parameter at destination surface.
754  // Direction will flip if dw2dw1 < 0.;
755 
756  switch (dir) {
757  case Surface::FORWARD: dir = (dw2dw1 > 0.) ? Surface::FORWARD : Surface::BACKWARD; break;
758  case Surface::BACKWARD: dir = (dw2dw1 > 0.) ? Surface::BACKWARD : Surface::FORWARD; break;
759  default:
760  throw cet::exception("PropXYZPlane")
761  << __func__ << ": unexpected direction #" << ((int)dir) << "\n";
762  } // switch
763 
764  // Update propagation matrix (if requested).
765 
766  if (prop_matrix != 0) {
767  TrackMatrix& pm = *prop_matrix;
768  pm.resize(vec.size(), vec.size(), false);
769 
770  // Calculate partial derivatives.
771 
772  pm(0, 0) = ruu - dudw2 * rwu; // du2/du1
773  pm(1, 0) = rvu - dvdw2 * rwu; // dv2/du1
774  pm(2, 0) = 0.; // d(dudw2)/du1
775  pm(3, 0) = 0.; // d(dvdw2)/du1
776  pm(4, 0) = 0.; // d(pinv2)/du1
777 
778  pm(0, 1) = ruv - dudw2 * rwv; // du2/dv1
779  pm(1, 1) = rvv - dvdw2 * rwv; // dv2/dv1
780  pm(2, 1) = 0.; // d(dudw2)/dv1
781  pm(3, 1) = 0.; // d(dvdw2)/dv1
782  pm(4, 1) = 0.; // d(pinv2)/dv1
783 
784  pm(0, 2) = 0.; // du2/d(dudw1);
785  pm(1, 2) = 0.; // dv2/d(dudw1);
786  pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1; // d(dudw2)/d(dudw1);
787  pm(3, 2) = (rvu - dvdw2 * rwu) / dw2dw1; // d(dvdw2)/d(dudw1);
788  pm(4, 2) = 0.; // d(pinv2)/d(dudw1);
789 
790  pm(0, 3) = 0.; // du2/d(dvdw1);
791  pm(1, 3) = 0.; // dv2/d(dvdw1);
792  pm(2, 3) = (ruv - dudw2 * rwv) / dw2dw1; // d(dudw2)/d(dvdw1);
793  pm(3, 3) = (rvv - dvdw2 * rwv) / dw2dw1; // d(dvdw2)/d(dvdw1);
794  pm(4, 3) = 0.; // d(pinv2)/d(dvdw1);
795 
796  pm(0, 4) = 0.; // du2/d(pinv1);
797  pm(1, 4) = 0.; // dv2/d(pinv1);
798  pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
799  pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
800  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
801  }
802 
803  // Update track vector.
804 
805  vec(0) = 0.;
806  vec(1) = 0.;
807  vec(2) = dudw2;
808  vec(3) = dvdw2;
809 
810  // Done (success).
811 
812  return true;
813  }
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
TDirectory * dir
Definition: macro.C:5
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool trkf::PropXYZPlane::transformYZLine ( double  phi1,
double  theta2,
double  phi2,
TrackVector vec,
Surface::TrackDirection dir,
TrackMatrix prop_matrix 
) const
private

Transform yz line -> xyz plane.

The following methods transform the track parameters from initial surface to SurfXYZPlane 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 376 of file PropXYZPlane.cxx.

Referenced by origin_vec_prop().

382  {
383  // Calculate surface transcendental functions.
384 
385  double sinth2 = std::sin(theta2);
386  double costh2 = std::cos(theta2);
387 
388  double sindphi = std::sin(phi2 - phi1);
389  double cosdphi = std::cos(phi2 - phi1);
390 
391  // Get the initial track parameters.
392 
393  double r1 = vec(0);
394  double phid1 = vec(2);
395  double eta1 = vec(3);
396 
397  // Calculate elements of rotation matrix from initial coordinate
398  // system to destination coordinte system.
399 
400  double ruu = costh2;
401  double ruv = sinth2 * sindphi;
402  double ruw = -sinth2 * cosdphi;
403 
404  double rvv = cosdphi;
405  double rvw = sindphi;
406 
407  double rwu = sinth2;
408  double rwv = -costh2 * sindphi;
409  double rww = costh2 * cosdphi;
410 
411  // Calculate track transcendental functions.
412 
413  double sinphid1 = std::sin(phid1);
414  double cosphid1 = std::cos(phid1);
415  double sh1 = 1. / std::cosh(eta1); // sech(eta1)
416  double th1 = std::tanh(eta1);
417 
418  // Calculate initial position in Cartesian coordinates.
419 
420  double u1 = -r1 * sinphid1;
421  double w1 = r1 * cosphid1;
422 
423  // Calculate direction in source coordinate system.
424 
425  double du1 = sh1 * cosphid1;
426  double dv1 = th1;
427  double dw1 = sh1 * sinphid1;
428  //double duw2 = std::hypot(du2, dw2);
429 
430  // Rotate direction to destination coordinate system.
431 
432  double du2 = ruu * du1 + ruv * dv1 + ruw * dw1;
433  double dv2 = rvv * dv1 + rvw * dw1;
434  double dw2 = rwu * du1 + rwv * dv1 + rww * dw1;
435 
436  // Calculate the track direction relative to the destination surface.
437  // The track direction comes from the sign of dw2 (=dw/ds).
438  // If dw2 is zero, the destionation surface is unreachable, return failure.
439 
440  if (dw2 > 0.)
441  dir = Surface::TrackDirection::FORWARD;
442  else if (dw2 < 0.)
443  dir = Surface::TrackDirection::BACKWARD;
444  else
445  return false;
446 
447  // Calculate final track slope track parameters.
448 
449  double dudw2 = du2 / dw2;
450  double dvdw2 = dv2 / dw2;
451 
452  // Update propagation matrix (if requested).
453 
454  if (prop_matrix != 0) {
455  TrackMatrix& pm = *prop_matrix;
456  pm.resize(vec.size(), vec.size(), false);
457 
458  // Calculate partial derivatives.
459 
460  // Partials of initial positions and directions wrt initial t.p.'s.
461 
462  double du1dr1 = -sinphid1;
463  double du1dphi1 = -w1;
464 
465  double dw1dr1 = cosphid1;
466  double dw1dphi1 = u1;
467 
468  double ddu1dphi1 = -sinphid1 * sh1;
469  double ddu1deta1 = -cosphid1 * sh1 * th1;
470 
471  double ddv1deta1 = sh1 * sh1;
472 
473  double ddw1dphi1 = cosphid1 * sh1;
474  double ddw1deta1 = -sinphid1 * sh1 * th1;
475 
476  // Rotate partials to destination coordinate system.
477 
478  double du2dr1 = ruu * du1dr1 + ruw * dw1dr1;
479  double dv2dr1 = rvw * dw1dr1;
480  double dw2dr1 = rwu * du1dr1 + rww * dw1dr1;
481 
482  double du2dv1 = ruv;
483  double dv2dv1 = rvv;
484  double dw2dv1 = rwv;
485 
486  double du2dphi1 = ruu * du1dphi1 + ruw * dw1dphi1;
487  double dv2dphi1 = rvw * dw1dphi1;
488  double dw2dphi1 = rwu * du1dphi1 + rww * dw1dphi1;
489 
490  double ddu2dphi1 = ruu * ddu1dphi1 + ruw * ddw1dphi1;
491  double ddv2dphi1 = rvw * ddw1dphi1;
492  double ddw2dphi1 = rwu * ddu1dphi1 + rww * ddw1dphi1;
493 
494  double ddu2deta1 = ruu * ddu1deta1 + ruv * ddv1deta1 + ruw * ddw1deta1;
495  double ddv2deta1 = rvv * ddv1deta1 + rvw * ddw1deta1;
496  double ddw2deta1 = rwu * ddu1deta1 + rwv * ddv1deta1 + rww * ddw1deta1;
497 
498  // Partials of final slope t.p. wrt final position and direction.
499 
500  double ddudw2ddu2 = 1. / dw2;
501  double ddudw2ddw2 = -dudw2 / dw2;
502 
503  double ddvdw2ddv2 = 1. / dw2;
504  double ddvdw2ddw2 = -dvdw2 / dw2;
505 
506  // Partials of final slope t.p. wrt initial t.p.
507 
508  double ddudw2dphi1 = ddudw2ddu2 * ddu2dphi1 + ddudw2ddw2 * ddw2dphi1;
509  double ddudw2deta1 = ddudw2ddu2 * ddu2deta1 + ddudw2ddw2 * ddw2deta1;
510 
511  double ddvdw2dphi1 = ddvdw2ddv2 * ddv2dphi1 + ddvdw2ddw2 * ddw2dphi1;
512  double ddvdw2deta1 = ddvdw2ddv2 * ddv2deta1 + ddvdw2ddw2 * ddw2deta1;
513 
514  // We still need to calculate the corretion due to the dependence of the
515  // propagation distance on the initial track parameters. This correction is
516  // needed even though the actual propagation distance is zero.
517 
518  // This correction effects the u and v track parameters.
519 
520  // Partials of perpendicular propagation distance wrt initial t.p.
521 
522  double dstdr1 = -dw2dr1;
523  double dstdv1 = -dw2dv1;
524  double dstdphi1 = -dw2dphi1;
525 
526  // Calculate correction to u and v parameter partials wrt initial t.p. due to path length.
527 
528  du2dr1 += dstdr1 * dudw2;
529  du2dv1 += dstdv1 * dudw2;
530  du2dphi1 += dstdphi1 * dudw2;
531 
532  dv2dr1 += dstdr1 * dvdw2;
533  dv2dv1 += dstdv1 * dvdw2;
534  dv2dphi1 += dstdphi1 * dvdw2;
535 
536  // Fill derivative matrix.
537 
538  pm(0, 0) = du2dr1; // du2/dr1
539  pm(1, 0) = dv2dr1; // dv2/dr1
540  pm(2, 0) = 0.; // d(dudw2)/dr1
541  pm(3, 0) = 0.; // d(dvdw2)/dr1
542  pm(4, 0) = 0.; // d(pinv2)/dr1
543 
544  pm(0, 1) = du2dv1; // du2/dv1
545  pm(1, 1) = dv2dv1; // dv2/dv1
546  pm(2, 1) = 0.; // d(dudw2)/dv1
547  pm(3, 1) = 0.; // d(dvdw2)/dv1
548  pm(4, 1) = 0.; // d(pinv2)/dv1
549 
550  pm(0, 2) = du2dphi1; // du2/d(phi1);
551  pm(1, 2) = dv2dphi1; // dv2/d(phi1);
552  pm(2, 2) = ddudw2dphi1; // d(dudw2)/d(phi1);
553  pm(3, 2) = ddvdw2dphi1; // d(dvdw2)/d(phi1);
554  pm(4, 2) = 0.; // d(pinv2)/d(phi1);
555 
556  pm(0, 3) = 0.; // du2/d(eta1);
557  pm(1, 3) = 0.; // dv2/d(eta1);
558  pm(2, 3) = ddudw2deta1; // d(dudw2)/d(eta1);
559  pm(3, 3) = ddvdw2deta1; // d(dvdw2)/d(eta1);
560  pm(4, 3) = 0.; // d(pinv2)/d(eta1);
561 
562  pm(0, 4) = 0.; // du2/d(pinv1);
563  pm(1, 4) = 0.; // dv2/d(pinv1);
564  pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
565  pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
566  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
567  }
568 
569  // Update track vector.
570 
571  vec(0) = 0.;
572  vec(1) = 0.;
573  vec(2) = dudw2;
574  vec(3) = dvdw2;
575 
576  // Done (success).
577 
578  return true;
579  }
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
TDirectory * dir
Definition: macro.C:5
bool trkf::PropXYZPlane::transformYZPlane ( double  phi1,
double  theta2,
double  phi2,
TrackVector vec,
Surface::TrackDirection dir,
TrackMatrix prop_matrix 
) const
private

Transform yz plane -> xyz plane.

Definition at line 583 of file PropXYZPlane.cxx.

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

Referenced by origin_vec_prop().

589  {
590  // Calculate transcendental functions.
591 
592  double sinth2 = std::sin(theta2);
593  double costh2 = std::cos(theta2);
594 
595  double sindphi = std::sin(phi2 - phi1);
596  double cosdphi = std::cos(phi2 - phi1);
597 
598  // Get the initial track state vector and track parameters.
599 
600  double dudw1 = vec(2);
601  double dvdw1 = vec(3);
602 
603  // Make sure initial track has a valid direction.
604 
605  if (dir == Surface::UNKNOWN) return false;
606 
607  // Calculate elements of rotation matrix from initial coordinate
608  // system to destination coordinte system.
609 
610  double ruu = costh2;
611  double ruv = sinth2 * sindphi;
612  double ruw = -sinth2 * cosdphi;
613 
614  double rvv = cosdphi;
615  double rvw = sindphi;
616 
617  double rwu = sinth2;
618  double rwv = -costh2 * sindphi;
619  double rww = costh2 * cosdphi;
620 
621  // Calculate the derivative dw2/dw1;
622  // If dw2/dw1 == 0., that means the track is moving parallel
623  // to destination plane.
624  // In this case return propagation failure.
625 
626  double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
627  if (dw2dw1 == 0.) return false;
628 
629  // Calculate slope in destination plane coordinates.
630 
631  double dudw2 = (dudw1 * ruu + dvdw1 * ruv + ruw) / dw2dw1;
632  double dvdw2 = (dvdw1 * rvv + rvw) / dw2dw1;
633 
634  // Calculate direction parameter at destination surface.
635  // Direction will flip if dw2dw1 < 0.;
636 
637  switch (dir) {
638  case Surface::FORWARD: dir = (dw2dw1 > 0.) ? Surface::FORWARD : Surface::BACKWARD; break;
639  case Surface::BACKWARD: dir = (dw2dw1 > 0.) ? Surface::BACKWARD : Surface::FORWARD; break;
640  default:
641  throw cet::exception("PropXYZPlane")
642  << __func__ << ": unexpected direction #" << ((int)dir) << "\n";
643  } // switch
644 
645  // Update propagation matrix (if requested).
646 
647  if (prop_matrix != 0) {
648  TrackMatrix& pm = *prop_matrix;
649  pm.resize(vec.size(), vec.size(), false);
650 
651  // Calculate partial derivatives.
652 
653  pm(0, 0) = ruu - dudw2 * rwu; // du2/du1
654  pm(1, 0) = -dvdw2 * rwu; // dv2/du1
655  pm(2, 0) = 0.; // d(dudw2)/du1
656  pm(3, 0) = 0.; // d(dvdw2)/du1
657  pm(4, 0) = 0.; // d(pinv2)/du1
658 
659  pm(0, 1) = ruv - dudw2 * rwv; // du2/dv1
660  pm(1, 1) = rvv - dvdw2 * rwv; // dv2/dv1
661  pm(2, 1) = 0.; // d(dudw2)/dv1
662  pm(3, 1) = 0.; // d(dvdw2)/dv1
663  pm(4, 1) = 0.; // d(pinv2)/dv1
664 
665  pm(0, 2) = 0.; // du2/d(dudw1);
666  pm(1, 2) = 0.; // dv2/d(dudw1);
667  pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1; // d(dudw2)/d(dudw1);
668  pm(3, 2) = -dvdw2 * rwu / dw2dw1; // d(dvdw2)/d(dudw1);
669  pm(4, 2) = 0.; // d(pinv2)/d(dudw1);
670 
671  pm(0, 3) = 0.; // du2/d(dvdw1);
672  pm(1, 3) = 0.; // dv2/d(dvdw1);
673  pm(2, 3) = (ruv - dudw2 * rwv) / dw2dw1; // d(dudw2)/d(dvdw1);
674  pm(3, 3) = (rvv - dvdw2 * rwv) / dw2dw1; // d(dvdw2)/d(dvdw1);
675  pm(4, 3) = 0.; // d(pinv2)/d(dvdw1);
676 
677  pm(0, 4) = 0.; // du2/d(pinv1);
678  pm(1, 4) = 0.; // dv2/d(pinv1);
679  pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
680  pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
681  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
682  }
683 
684  // Update track vector.
685 
686  vec(0) = 0.;
687  vec(1) = 0.;
688  vec(2) = dudw2;
689  vec(3) = dvdw2;
690 
691  // Done (success).
692 
693  return true;
694  }
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
TDirectory * dir
Definition: macro.C:5
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
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: