Public Member Functions | Private Member Functions | List of all members
trkf::PropYZPlane Class Reference

#include <PropYZPlane.h>

Inheritance diagram for trkf::PropYZPlane:
trkf::Propagator

Public Member Functions

 PropYZPlane (detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx)
 
Propagatorclone () const override
 Clone method. More...
 
std::optional< double > short_vec_prop (KTrack &trk, const std::shared_ptr< const Surface > &surf, Propagator::PropDirection dir, bool doDedx, TrackMatrix *prop_matrix=0, TrackError *noise_matrix=0) const override
 Propagate without error. More...
 
std::optional< double > origin_vec_prop (KTrack &trk, const std::shared_ptr< const Surface > &porient, TrackMatrix *prop_matrix=0) const override
 Propagate without error to surface whose origin parameters coincide with track position. More...
 
- Public Member Functions inherited from trkf::Propagator
 Propagator (detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx, const std::shared_ptr< const Interactor > &interactor)
 Constructor. More...
 
virtual ~Propagator ()
 Destructor. More...
 
double getTcut () const
 
bool getDoDedx () const
 
const std::shared_ptr< const Interactor > & getInteractor () const
 
std::optional< double > vec_prop (KTrack &trk, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, TrackMatrix *prop_matrix=0, TrackError *noise_matrix=0) const
 Propagate without error (long distance). More...
 
std::optional< double > lin_prop (KTrack &trk, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, KTrack *ref=0, TrackMatrix *prop_matrix=0, TrackError *noise_matrix=0) const
 Linearized propagate without error. More...
 
std::optional< double > err_prop (KETrack &tre, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, KTrack *ref=0, TrackMatrix *prop_matrix=0) const
 Propagate with error, but without noise. More...
 
std::optional< double > noise_prop (KETrack &tre, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, KTrack *ref=0) const
 Propagate with error and noise. More...
 
std::optional< double > dedx_prop (double pinv, double mass, double s, double *deriv=0) const
 Method to calculate updated momentum due to dE/dx. More...
 

Private Member Functions

bool transformYZLine (double phi1, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
 Transform yz line -> yz plane. More...
 
bool transformYZPlane (double phi1, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
 Transform yz plane -> yz plane. More...
 
bool transformXYZPlane (double theta1, double phi1, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
 Transform xyz plane -> yz plane. More...
 

Additional Inherited Members

- Public Types inherited from trkf::Propagator
enum  PropDirection { FORWARD, BACKWARD, UNKNOWN }
 Propagation direction enum. More...
 

Detailed Description

Definition at line 24 of file PropYZPlane.h.

Constructor & Destructor Documentation

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

Constructor.

Arguments.

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

Definition at line 28 of file PropYZPlane.cxx.

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

Member Function Documentation

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

Clone method.

Implements trkf::Propagator.

Definition at line 30 of file PropYZPlane.h.

31  {
32  return new PropYZPlane(*this);
33  }
PropYZPlane(detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx)
Definition: PropYZPlane.cxx:28
std::optional< double > trkf::PropYZPlane::origin_vec_prop ( KTrack trk,
const std::shared_ptr< const Surface > &  porient,
TrackMatrix prop_matrix = 0 
) const
overridevirtual

Propagate without error to surface whose origin parameters coincide with track position.

Propagate without error to dynamically generated origin surface. Optionally return propagation matrix.

Arguments:

trk - Track to propagate. porient - Orientation surface. prop_matrix - Pointer to optional propagation matrix.

Returned value: propagation distance + success flag.

Propagation distance is always zero after successful propagation.

Implements trkf::Propagator.

Definition at line 259 of file PropYZPlane.cxx.

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

Propagate without error.

Propagate without error. Optionally return propagation matrix and noise matrix.

Arguments:

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

Returned value: propagation distance + success flag.

Implements trkf::Propagator.

Definition at line 51 of file PropYZPlane.cxx.

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

Transform xyz plane -> yz plane.

Definition at line 648 of file PropYZPlane.cxx.

654  {
655  // Calculate transcendental functions.
656 
657  double sinth1 = std::sin(theta1);
658  double costh1 = std::cos(theta1);
659 
660  double sindphi = std::sin(phi2 - phi1);
661  double cosdphi = std::cos(phi2 - phi1);
662 
663  // Get the initial track state vector and track parameters.
664 
665  double dudw1 = vec(2);
666  double dvdw1 = vec(3);
667 
668  // Make sure initial track has a valid direction.
669 
670  if (dir == Surface::UNKNOWN) return false;
671 
672  // Calculate elements of rotation matrix from initial coordinate
673  // system to destination coordinte system.
674 
675  double ruu = costh1;
676  double ruw = sinth1;
677 
678  double rvu = -sinth1 * sindphi;
679  double rvv = cosdphi;
680  double rvw = costh1 * sindphi;
681 
682  double rwu = -sinth1 * cosdphi;
683  double rwv = -sindphi;
684  double rww = costh1 * cosdphi;
685 
686  // Calculate the derivative dw2/dw1;
687  // If dw2/dw1 == 0., that means the track is moving parallel
688  // to destination plane.
689  // In this case return propagation failure.
690 
691  double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
692  if (dw2dw1 == 0.) return false;
693 
694  // Calculate slope in destination plane coordinates.
695 
696  double dudw2 = (dudw1 * ruu + ruw) / dw2dw1;
697  double dvdw2 = (dudw1 * rvu + dvdw1 * rvv + rvw) / dw2dw1;
698 
699  // Calculate direction parameter at destination surface.
700  // Direction will flip if dw2dw1 < 0.;
701 
702  switch (dir) {
703  case Surface::FORWARD: dir = (dw2dw1 > 0.) ? Surface::FORWARD : Surface::BACKWARD; break;
704  case Surface::BACKWARD: dir = (dw2dw1 > 0.) ? Surface::BACKWARD : Surface::FORWARD; break;
705  default:
706  throw cet::exception("PropXYZPlane")
707  << __func__ << ": unexpected direction #" << ((int)dir) << "\n";
708  } // switch
709 
710  // Update propagation matrix (if requested).
711 
712  if (prop_matrix != 0) {
713  TrackMatrix& pm = *prop_matrix;
714  pm.resize(vec.size(), vec.size(), false);
715 
716  // Calculate partial derivatives.
717 
718  pm(0, 0) = ruu - dudw2 * rwu; // du2/du1
719  pm(1, 0) = rvu - dvdw2 * rwu; // dv2/du1
720  pm(2, 0) = 0.; // d(dudw2)/du1
721  pm(3, 0) = 0.; // d(dvdw2)/du1
722  pm(4, 0) = 0.; // d(pinv2)/du1
723 
724  pm(0, 1) = -dudw2 * rwv; // du2/dv1
725  pm(1, 1) = rvv - dvdw2 * rwv; // dv2/dv1
726  pm(2, 1) = 0.; // d(dudw2)/dv1
727  pm(3, 1) = 0.; // d(dvdw2)/dv1
728  pm(4, 1) = 0.; // d(pinv2)/dv1
729 
730  pm(0, 2) = 0.; // du2/d(dudw1);
731  pm(1, 2) = 0.; // dv2/d(dudw1);
732  pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1; // d(dudw2)/d(dudw1);
733  pm(3, 2) = (rvu - dvdw2 * rwu) / dw2dw1; // d(dvdw2)/d(dudw1);
734  pm(4, 2) = 0.; // d(pinv2)/d(dudw1);
735 
736  pm(0, 3) = 0.; // du2/d(dvdw1);
737  pm(1, 3) = 0.; // dv2/d(dvdw1);
738  pm(2, 3) = -dudw2 * rwv / dw2dw1; // d(dudw2)/d(dvdw1);
739  pm(3, 3) = (rvv - dvdw2 * rwv) / dw2dw1; // d(dvdw2)/d(dvdw1);
740  pm(4, 3) = 0.; // d(pinv2)/d(dvdw1);
741 
742  pm(0, 4) = 0.; // du2/d(pinv1);
743  pm(1, 4) = 0.; // dv2/d(pinv1);
744  pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
745  pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
746  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
747  }
748 
749  // Update track vector.
750 
751  vec(0) = 0.;
752  vec(1) = 0.;
753  vec(2) = dudw2;
754  vec(3) = dvdw2;
755 
756  // Done (success).
757 
758  return true;
759  }
string dir
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool trkf::PropYZPlane::transformYZLine ( double  phi1,
double  phi2,
TrackVector vec,
Surface::TrackDirection dir,
TrackMatrix prop_matrix 
) const
private

Transform yz line -> yz plane.

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

Definition at line 360 of file PropYZPlane.cxx.

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

Transform yz plane -> yz plane.

Definition at line 552 of file PropYZPlane.cxx.

557  {
558  // Calculate transcendental functions.
559 
560  double sindphi = std::sin(phi2 - phi1);
561  double cosdphi = std::cos(phi2 - phi1);
562 
563  // Get the initial track parameters.
564 
565  double dudw1 = vec(2);
566  double dvdw1 = vec(3);
567 
568  // Make sure initial track has a valid direction.
569 
570  if (dir == Surface::UNKNOWN) return false;
571 
572  // Calculate derivative dw2/dw1.
573  // If dw2/dw1 == 0., that means the track is moving parallel
574  // to destination plane.
575  // In this case return propagation failure.
576 
577  double dw2dw1 = cosdphi - dvdw1 * sindphi;
578  if (dw2dw1 == 0.) return false;
579 
580  // Calculate slope in destrination coordiante system.
581 
582  double dudw2 = dudw1 / dw2dw1;
583  double dvdw2 = (sindphi + dvdw1 * cosdphi) / dw2dw1;
584 
585  // Calculate direction parameter at destination surface.
586  // Direction will flip if dw2dw1 < 0.;
587 
588  switch (dir) {
589  case Surface::FORWARD: dir = (dw2dw1 > 0.) ? Surface::FORWARD : Surface::BACKWARD; break;
590  case Surface::BACKWARD: dir = (dw2dw1 > 0.) ? Surface::BACKWARD : Surface::FORWARD; break;
591  default: throw cet::exception("PropYZPlane") << "unexpected direction #" << ((int)dir) << "\n";
592  } // switch
593 
594  // Update propagation matrix (if requested).
595 
596  if (prop_matrix != 0) {
597  TrackMatrix& pm = *prop_matrix;
598  pm.resize(vec.size(), vec.size(), false);
599 
600  // Calculate partial derivatives.
601 
602  pm(0, 0) = 1.; // du2/du1
603  pm(1, 0) = 0.; // dv2/du1
604  pm(2, 0) = 0.; // d(dudw2)/du1
605  pm(3, 0) = 0.; // d(dvdw2)/du1
606  pm(4, 0) = 0.; // d(pinv2)/du1
607 
608  pm(0, 1) = dudw2 * sindphi; // du2/dv1
609  pm(1, 1) = cosdphi + dvdw2 * sindphi; // dv2/dv1
610  pm(2, 1) = 0.; // d(dudw2)/dv1
611  pm(3, 1) = 0.; // d(dvdw2)/dv1
612  pm(4, 1) = 0.; // d(pinv2)/dv1
613 
614  pm(0, 2) = 0.; // du2/d(dudw1);
615  pm(1, 2) = 0.; // dv2/d(dudw1);
616  pm(2, 2) = 1. / dw2dw1; // d(dudw2)/d(dudw1);
617  pm(3, 2) = 0.; // d(dvdw2)/d(dudw1);
618  pm(4, 2) = 0.; // d(pinv2)/d(dudw1);
619 
620  pm(0, 3) = 0.; // du2/d(dvdw1);
621  pm(1, 3) = 0.; // dv2/d(dvdw1);
622  pm(2, 3) = dudw1 * sindphi / (dw2dw1 * dw2dw1); // d(dudw2)/d(dvdw1);
623  pm(3, 3) = 1. / (dw2dw1 * dw2dw1); // d(dvdw2)/d(dvdw1);
624  pm(4, 3) = 0.; // d(pinv2)/d(dvdw1);
625 
626  pm(0, 4) = 0.; // du2/d(pinv1);
627  pm(1, 4) = 0.; // dv2/d(pinv1);
628  pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
629  pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
630  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
631  }
632 
633  // Update track vector.
634 
635  vec(0) = 0.;
636  vec(1) = 0.;
637  vec(2) = dudw2;
638  vec(3) = dvdw2;
639 
640  // Done (success).
641 
642  return true;
643  }
string dir
KMatrix< 5, 5 >::type TrackMatrix
General 5x5 matrix.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

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