LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
trkf::PropAny Class Reference

#include "PropAny.h"

Inheritance diagram for trkf::PropAny:
trkf::Propagator

Public Types

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

Public Member Functions

 PropAny (double tcut, bool doDedx)
 Constructor. More...
 
virtual ~PropAny ()
 Destructor. More...
 
Propagatorclone () const
 Clone method. More...
 
boost::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
 Propagate without error. More...
 
virtual boost::optional< double > origin_vec_prop (KTrack &trk, const std::shared_ptr< const Surface > &porient, TrackMatrix *prop_matrix=0) const
 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
 
boost::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...
 
boost::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...
 
boost::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...
 
boost::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...
 
boost::optional< double > dedx_prop (double pinv, double mass, double s, double *deriv=0) const
 Method to calculate updated momentum due to dE/dx. More...
 

Private Attributes

PropYZLine fPropYZLine
 Underlying propagators. More...
 
PropYZPlane fPropYZPlane
 
PropXYZPlane fPropXYZPlane
 

Detailed Description

Definition at line 25 of file PropAny.h.

Member Enumeration Documentation

Propagation direction enum.

Enumerator
FORWARD 
BACKWARD 
UNKNOWN 

Definition at line 92 of file Propagator.h.

Constructor & Destructor Documentation

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

Referenced by clone().

28  :
29  Propagator(tcut, doDedx, (tcut >= 0 ?
30  std::shared_ptr<const Interactor>(new InteractPlane(tcut)) :
31  std::shared_ptr<const Interactor>())),
32  fPropYZLine(tcut, doDedx),
33  fPropYZPlane(tcut, doDedx),
34  fPropXYZPlane(tcut, doDedx)
35  {}
PropYZPlane fPropYZPlane
Definition: PropAny.h:60
PropYZLine fPropYZLine
Underlying propagators.
Definition: PropAny.h:59
Propagator(double tcut, bool doDedx, const std::shared_ptr< const Interactor > &interactor)
Constructor.
Definition: Propagator.cxx:25
PropXYZPlane fPropXYZPlane
Definition: PropAny.h:61
trkf::PropAny::~PropAny ( )
virtual

Destructor.

Definition at line 38 of file PropAny.cxx.

39  {}

Member Function Documentation

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

Clone method.

Implements trkf::Propagator.

Definition at line 38 of file PropAny.h.

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

38 {return new PropAny(*this);}
PropAny(double tcut, bool doDedx)
Constructor.
Definition: PropAny.cxx:28
boost::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 459 of file Propagator.cxx.

References trkf::Propagator::fTcut.

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

461  {
462  // For infinite initial momentum, return with success status,
463  // still infinite momentum.
464 
465  if(pinv == 0.)
466  return boost::optional<double>(true, 0.);
467 
468  // Set the default return value to be uninitialized with value 0.
469 
470  boost::optional<double> result(false, 0.);
471 
472  // Get LAr service.
473 
474  auto const * detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
475 
476  // Calculate final energy.
477 
478  double p1 = 1./std::abs(pinv);
479  double e1 = std::hypot(p1, mass);
480  double de = -0.001 * s * detprop->Eloss(p1, mass, fTcut);
481  double emid = e1 + 0.5 * de;
482  if(emid > mass) {
483  double pmid = std::sqrt(emid*emid - mass*mass);
484  double e2 = e1 - 0.001 * s * detprop->Eloss(pmid, mass, fTcut);
485  if(e2 > mass) {
486  double p2 = std::sqrt(e2*e2 - mass*mass);
487  double pinv2 = 1./p2;
488  if(pinv < 0.)
489  pinv2 = -pinv2;
490 
491  // Calculation was successful, update result.
492 
493  result = boost::optional<double>(true, pinv2);
494 
495  // Also calculate derivative, if requested.
496 
497  if(deriv != 0)
498  *deriv = pinv2*pinv2*pinv2 * e2 / (pinv*pinv*pinv * e1);
499  }
500  }
501 
502  // Done.
503 
504  return result;
505  }
Float_t s
Definition: plot.C:23
double fTcut
Maximum delta ray energy for dE/dx.
Definition: Propagator.h:163
boost::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 354 of file Propagator.cxx.

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

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

360  {
361  // Propagate without error, get propagation matrix.
362 
363  TrackMatrix prop_temp;
364  if(prop_matrix == 0)
365  prop_matrix = &prop_temp;
366  boost::optional<double> result = lin_prop(tre, psurf, dir, doDedx, ref, prop_matrix, 0);
367 
368  // If propagation succeeded, update track error matrix.
369 
370  if(!!result) {
371  TrackMatrix temp = prod(tre.getError(), trans(*prop_matrix));
372  TrackMatrix temp2 = prod(*prop_matrix, temp);
373  TrackError newerr = ublas::symmetric_adaptor<TrackMatrix>(temp2);
374  tre.setError(newerr);
375  }
376 
377  // Done.
378 
379  return result;
380  }
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
TDirectory * dir
Definition: macro.C:5
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
boost::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:258
bool trkf::Propagator::getDoDedx ( ) const
inlineinherited
double trkf::Propagator::getTcut ( ) const
inlineinherited

Definition at line 102 of file Propagator.h.

References trkf::Propagator::fTcut.

102 {return fTcut;}
double fTcut
Maximum delta ray energy for dE/dx.
Definition: Propagator.h:163
boost::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 258 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(), trkf::Propagator::getInteractor(), and trkf::Propagator::noise_prop().

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

399  {
400  // Propagate without error, get propagation matrix and noise matrix.
401 
402  TrackMatrix prop_matrix;
403  TrackError noise_matrix;
404  boost::optional<double> result = lin_prop(tre, psurf, dir, doDedx, ref,
405  &prop_matrix, &noise_matrix);
406 
407  // If propagation succeeded, update track error matrix.
408 
409  if(!!result) {
410  TrackMatrix temp = prod(tre.getError(), trans(prop_matrix));
411  TrackMatrix temp2 = prod(prop_matrix, temp);
412  TrackError newerr = ublas::symmetric_adaptor<TrackMatrix>(temp2);
413  newerr += noise_matrix;
414  tre.setError(newerr);
415  }
416 
417  // Done.
418 
419  return result;
420  }
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
TDirectory * dir
Definition: macro.C:5
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
boost::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:258
boost::optional< double > trkf::PropAny::origin_vec_prop ( KTrack trk,
const std::shared_ptr< const Surface > &  porient,
TrackMatrix prop_matrix = 0 
) const
virtual

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 99 of file PropAny.cxx.

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

Referenced by clone().

102  {
103  // Default result.
104 
105  boost::optional<double> result(false, 0.);
106 
107  // Test the type of the destination surface.
108 
109  if(dynamic_cast<const SurfYZLine*>(&*porient))
110  result = fPropYZLine.origin_vec_prop(trk, porient, prop_matrix);
111  else if(dynamic_cast<const SurfYZPlane*>(&*porient))
112  result = fPropYZPlane.origin_vec_prop(trk, porient, prop_matrix);
113  else if(dynamic_cast<const SurfXYZPlane*>(&*porient))
114  result = fPropXYZPlane.origin_vec_prop(trk, porient, prop_matrix);
115  else
116  throw cet::exception("PropAny") << "Destination surface has unknown type.\n";
117 
118  // Done.
119 
120  return result;
121  }
virtual boost::optional< double > origin_vec_prop(KTrack &trk, const std::shared_ptr< const Surface > &porient, TrackMatrix *prop_matrix=0) const
Propagate without error to surface whose origin parameters coincide with track position.
Definition: PropYZLine.cxx:272
PropYZPlane fPropYZPlane
Definition: PropAny.h:60
PropYZLine fPropYZLine
Underlying propagators.
Definition: PropAny.h:59
PropXYZPlane fPropXYZPlane
Definition: PropAny.h:61
virtual boost::optional< double > origin_vec_prop(KTrack &trk, const std::shared_ptr< const Surface > &porient, TrackMatrix *prop_matrix=0) const
Propagate without error to surface whose origin parameters coincide with track position.
virtual boost::optional< double > origin_vec_prop(KTrack &trk, const std::shared_ptr< const Surface > &porient, TrackMatrix *prop_matrix=0) const
Propagate without error to surface whose origin parameters coincide with track position.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
boost::optional< double > trkf::PropAny::short_vec_prop ( KTrack trk,
const std::shared_ptr< const Surface > &  psurf,
Propagator::PropDirection  dir,
bool  doDedx,
TrackMatrix prop_matrix = 0,
TrackError noise_matrix = 0 
) const
virtual

Propagate without error.

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

Arguments:

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

Returned value: propagation distance + success flag.

Implements trkf::Propagator.

Definition at line 58 of file PropAny.cxx.

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

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

64  {
65  // Default result.
66 
67  boost::optional<double> result(false, 0.);
68 
69  // Test the type of the destination surface.
70 
71  if(dynamic_cast<const SurfYZLine*>(&*psurf))
72  result = fPropYZLine.short_vec_prop(trk, psurf, dir, doDedx, prop_matrix, noise_matrix);
73  else if(dynamic_cast<const SurfYZPlane*>(&*psurf))
74  result = fPropYZPlane.short_vec_prop(trk, psurf, dir, doDedx, prop_matrix, noise_matrix);
75  else if(dynamic_cast<const SurfXYZPlane*>(&*psurf))
76  result = fPropXYZPlane.short_vec_prop(trk, psurf, dir, doDedx, prop_matrix, noise_matrix);
77  else
78  throw cet::exception("PropAny") << "Destination surface has unknown type.\n";
79 
80  // Done.
81 
82  return result;
83  }
boost::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
Propagate without error.
PropYZPlane fPropYZPlane
Definition: PropAny.h:60
boost::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
Propagate without error.
Definition: PropYZLine.cxx:53
PropYZLine fPropYZLine
Underlying propagators.
Definition: PropAny.h:59
PropXYZPlane fPropXYZPlane
Definition: PropAny.h:61
TDirectory * dir
Definition: macro.C:5
boost::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
Propagate without error.
Definition: PropYZPlane.cxx:53
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
boost::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 e, trkf::Propagator::fTcut, trkf::Propagator::getDoDedx(), trkf::KTrack::getMomentum(), trkf::KTrack::getPosition(), trkf::KTrack::getVector(), trkf::KTrack::Mass(), s, and trkf::Propagator::short_vec_prop().

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

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

Member Data Documentation

PropXYZPlane trkf::PropAny::fPropXYZPlane
private

Definition at line 61 of file PropAny.h.

Referenced by origin_vec_prop(), and short_vec_prop().

PropYZLine trkf::PropAny::fPropYZLine
private

Underlying propagators.

Definition at line 59 of file PropAny.h.

Referenced by origin_vec_prop(), and short_vec_prop().

PropYZPlane trkf::PropAny::fPropYZPlane
private

Definition at line 60 of file PropAny.h.

Referenced by origin_vec_prop(), and short_vec_prop().


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