LArSoft  v07_13_02
Liquid Argon Software toolkit - http://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 (double tcut, bool doDedx)
 Constructor. More...
 
virtual ~PropXYZPlane ()
 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 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 20 of file PropXYZPlane.h.

Member Enumeration Documentation

Propagation direction enum.

Enumerator
FORWARD 
BACKWARD 
UNKNOWN 

Definition at line 92 of file Propagator.h.

Constructor & Destructor Documentation

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

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  {}
Propagator(double tcut, bool doDedx, const std::shared_ptr< const Interactor > &interactor)
Constructor.
Definition: Propagator.cxx:25
trkf::PropXYZPlane::~PropXYZPlane ( )
virtual

Destructor.

Definition at line 35 of file PropXYZPlane.cxx.

36  {}

Member Function Documentation

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

Clone method.

Implements trkf::Propagator.

Definition at line 33 of file PropXYZPlane.h.

References dir, origin_vec_prop(), PropXYZPlane(), short_vec_prop(), transformXYZPlane(), transformYZLine(), and transformYZPlane().

33 {return new PropXYZPlane(*this);}
PropXYZPlane(double tcut, bool doDedx)
Constructor.
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::PropYZPlane::short_vec_prop(), trkf::PropYZLine::short_vec_prop(), and 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

Definition at line 103 of file Propagator.h.

References trkf::Propagator::fDoDedx.

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

103 {return fDoDedx;}
bool fDoDedx
Energy loss enable flag.
Definition: Propagator.h:164
const std::shared_ptr<const Interactor>& trkf::Propagator::getInteractor ( ) 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::PropXYZPlane::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 281 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 clone(), trkf::PropAny::origin_vec_prop(), and short_vec_prop().

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

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 53 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(), s, 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 clone(), and trkf::PropAny::short_vec_prop().

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

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

Referenced by clone(), and origin_vec_prop().

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

Referenced by clone(), and origin_vec_prop().

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

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

Referenced by clone(), and origin_vec_prop().

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

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