PropXYZPlane.cxx
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 ///
3 /// \file PropXYZPlane.cxx
4 ///
5 /// \brief Propagate to SurfXYZPlane surface.
6 ///
7 /// \author H. Greenlee
8 ///
9 ////////////////////////////////////////////////////////////////////////
10 
12 #include "cetlib_except/exception.h"
17 #include <cmath>
18 
19 namespace trkf {
20 
21  /// Constructor.
22  ///
23  /// Arguments.
24  ///
25  /// tcut - Delta ray energy cutoff for calculating dE/dx.
26  /// doDedx - dE/dx enable flag.
27  ///
29  double tcut,
30  bool doDedx)
31  : Propagator{detProp,
32  tcut,
33  doDedx,
34  (tcut >= 0. ? std::make_shared<InteractPlane const>(detProp, tcut) :
35  std::shared_ptr<Interactor const>{})}
36  {}
37 
38  /// Propagate without error.
39  /// Optionally return propagation matrix and noise matrix.
40  ///
41  /// Arguments:
42  ///
43  /// trk - Track to propagate.
44  /// psurf - Destination surface.
45  /// dir - Propagation direction (FORWARD, BACKWARD, or UNKNOWN).
46  /// doDedx - dE/dx enable/disable flag.
47  /// prop_matrix - Pointer to optional propagation matrix.
48  /// noise_matrix - Pointer to optional noise matrix.
49  ///
50  /// Returned value: propagation distance + success flag.
51  ///
52  std::optional<double>
54  const std::shared_ptr<const Surface>& psurf,
56  bool doDedx,
57  TrackMatrix* prop_matrix,
58  TrackError* noise_matrix) const
59  {
60  // Set the default return value to be unitialized with value 0.
61 
62  std::optional<double> result{std::nullopt};
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) return result;
69  double x02 = to->x0();
70  double y02 = to->y0();
71  double z02 = to->z0();
72  double theta2 = to->theta();
73  double phi2 = to->phi();
74 
75  // Remember starting track.
76 
77  KTrack trk0(trk);
78 
79  // Get track position.
80 
81  double xyz[3];
82  trk.getPosition(xyz);
83  double x01 = xyz[0];
84  double y01 = xyz[1];
85  double z01 = xyz[2];
86 
87  // Propagate to origin surface.
88 
89  TrackMatrix local_prop_matrix;
90  TrackMatrix* plocal_prop_matrix = (prop_matrix == 0 ? 0 : &local_prop_matrix);
91  std::optional<double> result1 = origin_vec_prop(trk, psurf, plocal_prop_matrix);
92  if (!result1) return result1;
93 
94  // Get the intermediate track state vector and track parameters.
95 
96  const TrackVector& vec = trk.getVector();
97  if (vec.size() != 5)
98  throw cet::exception("PropXYZPlane")
99  << "Track state vector has wrong size" << vec.size() << "\n";
100  double u1 = vec(0);
101  double v1 = vec(1);
102  double dudw1 = vec(2);
103  double dvdw1 = vec(3);
104  double pinv = vec(4);
106 
107  // Make sure intermediate track has a valid direction.
108 
109  if (dir1 == Surface::UNKNOWN) {
110  trk = trk0;
111  return result;
112  }
113 
114  // Calculate transcendental functions.
115 
116  double sinth2 = std::sin(theta2);
117  double costh2 = std::cos(theta2);
118  double sinphi2 = std::sin(phi2);
119  double cosphi2 = std::cos(phi2);
120 
121  // Calculate elements of rotation matrix from global coordinate
122  // system to destination coordinate system.
123 
124  double rux = costh2;
125  double ruy = sinth2 * sinphi2;
126  double ruz = -sinth2 * cosphi2;
127 
128  double rvy = cosphi2;
129  double rvz = sinphi2;
130 
131  double rwx = sinth2;
132  double rwy = -costh2 * sinphi2;
133  double rwz = costh2 * cosphi2;
134 
135  // Calculate the initial position in the destination coordinate
136  // system.
137 
138  double u2 = (x01 - x02) * rux + (y01 - y02) * ruy + (z01 - z02) * ruz + u1;
139  double v2 = (y01 - y02) * rvy + (z01 - z02) * rvz + v1;
140  double w2 = (x01 - x02) * rwx + (y01 - y02) * rwy + (z01 - z02) * rwz;
141 
142  // Calculate position at destination surface (propagate distance -w2).
143 
144  double u2p = u2 - w2 * dudw1;
145  double v2p = v2 - w2 * dvdw1;
146 
147  // Calculate the signed propagation distance.
148 
149  double s = -w2 * std::sqrt(1. + dudw1 * dudw1 + dvdw1 * dvdw1);
150  if (dir1 == Surface::BACKWARD) s = -s;
151 
152  // Check if propagation was in the right direction.
153  // (Compare sign of s with requested direction).
154 
155  bool sok = (dir == Propagator::UNKNOWN || (dir == Propagator::FORWARD && s >= 0.) ||
156  (dir == Propagator::BACKWARD && s <= 0.));
157 
158  // If wrong direction, return failure without updating the track
159  // or propagation matrix.
160 
161  if (!sok) {
162  trk = trk0;
163  return result;
164  }
165 
166  // Find final momentum.
167  double deriv = 1.;
168  auto pinv2 = std::make_optional(pinv);
169  if (getDoDedx() && doDedx && s != 0.) {
170  double* pderiv = (prop_matrix != 0 ? &deriv : 0);
171  pinv2 = dedx_prop(pinv, trk.Mass(), s, pderiv);
172  }
173 
174  // Return failure in case of range out.
175 
176  if (!pinv2) {
177  trk = trk0;
178  return result;
179  }
180 
181  // Update default result to success and store propagation distance.
182 
183  result = std::make_optional(s);
184 
185  // Update propagation matrix (if requested).
186 
187  if (prop_matrix != 0) {
188  TrackMatrix pm;
189  pm.resize(vec.size(), vec.size(), false);
190 
191  // Calculate partial derivatives.
192 
193  pm(0, 0) = 1.; // du2/du1
194  pm(1, 0) = 0.; // dv2/du1
195  pm(2, 0) = 0.; // d(dudw2)/du1
196  pm(3, 0) = 0.; // d(dvdw2)/du1
197  pm(4, 0) = 0.; // d(pinv2)/du1
198 
199  pm(0, 1) = 0.; // du2/dv1
200  pm(1, 1) = 1.; // dv2/dv1
201  pm(2, 1) = 0.; // d(dudw2)/dv1
202  pm(3, 1) = 0.; // d(dvdw2)/dv1
203  pm(4, 1) = 0.; // d(pinv2)/dv1
204 
205  pm(0, 2) = -w2; // du2/d(dudw1);
206  pm(1, 2) = 0.; // dv2/d(dudw1);
207  pm(2, 2) = 1.; // d(dudw2)/d(dudw1);
208  pm(3, 2) = 0.; // d(dvdw2)/d(dudw1);
209  pm(4, 2) = 0.; // d(pinv2)/d(dudw1);
210 
211  pm(0, 3) = 0.; // du2/d(dvdw1);
212  pm(1, 3) = -w2; // dv2/d(dvdw1);
213  pm(2, 3) = 0.; // d(dudw2)/d(dvdw1);
214  pm(3, 3) = 1.; // d(dvdw2)/d(dvdw1);
215  pm(4, 3) = 0.; // d(pinv2)/d(dvdw1);
216 
217  pm(0, 4) = 0.; // du2/d(pinv1);
218  pm(1, 4) = 0.; // dv2/d(pinv1);
219  pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
220  pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
221  pm(4, 4) = deriv; // d(pinv2)/d(pinv1);
222 
223  // Compose the final propagation matrix from zero-distance propagation and
224  // parallel surface propagation.
225 
226  *prop_matrix = prod(pm, *plocal_prop_matrix);
227  }
228 
229  // Update noise matrix (if requested).
230 
231  if (noise_matrix != 0) {
232  noise_matrix->resize(vec.size(), vec.size(), false);
233  if (getInteractor().get() != 0) {
234  bool ok = getInteractor()->noise(trk, s, *noise_matrix);
235  if (!ok) {
236  trk = trk0;
237  return std::nullopt;
238  }
239  }
240  else
241  noise_matrix->clear();
242  }
243 
244  // Construct track vector at destination surface.
245 
246  TrackVector vec2(vec.size());
247  vec2(0) = u2p;
248  vec2(1) = v2p;
249  vec2(2) = dudw1;
250  vec2(3) = dvdw1;
251  vec2(4) = *pinv2;
252 
253  // Update track.
254 
255  trk.setSurface(psurf);
256  trk.setVector(vec2);
257 
258  // Done.
259 
260  return result;
261  }
262 
263  /// Propagate without error to dynamically generated origin surface.
264  /// Optionally return propagation matrix.
265  ///
266  /// Arguments:
267  ///
268  /// trk - Track to propagate.
269  /// porient - Orientation surface.
270  /// prop_matrix - Pointer to optional propagation matrix.
271  ///
272  /// Returned value: propagation distance + success flag.
273  ///
274  /// Propagation distance is always zero after successful propagation.
275  ///
276  std::optional<double>
278  const std::shared_ptr<const Surface>& porient,
279  TrackMatrix* prop_matrix) const
280  {
281  // Set the default return value to be unitialized with value 0.
282 
283  std::optional<double> result{std::nullopt};
284 
285  // Remember starting track.
286 
287  KTrack trk0(trk);
288 
289  // Get initial track parameters and direction.
290  // Note the initial track can be on any type of surface.
291 
292  TrackVector vec = trk.getVector(); // Modifiable copy.
293  if (vec.size() != 5)
294  throw cet::exception("PropYZPlane")
295  << "Track state vector has wrong size" << vec.size() << "\n";
297 
298  // Get track position.
299 
300  double xyz[3];
301  trk.getPosition(xyz);
302  double x02 = xyz[0];
303  double y02 = xyz[1];
304  double z02 = xyz[2];
305 
306  // Generate the origin surface, which will be the destination surface.
307  // Return failure if orientation surface is the wrong type.
308 
309  const SurfXYZPlane* orient = dynamic_cast<const SurfXYZPlane*>(&*porient);
310  if (orient == 0) return result;
311  double theta2 = orient->theta();
312  double phi2 = orient->phi();
313  std::shared_ptr<const Surface> porigin(new SurfXYZPlane(x02, y02, z02, phi2, theta2));
314 
315  // Test initial surface types.
316 
317  if (const SurfYZLine* from = dynamic_cast<const SurfYZLine*>(&*trk.getSurface())) {
318 
319  // Initial surface is SurfYZLine.
320  // Get surface paramters.
321 
322  double phi1 = from->phi();
323 
324  // Transform track to origin surface.
325 
326  bool ok = transformYZLine(phi1, theta2, phi2, vec, dir, prop_matrix);
327  result = std::make_optional(0.);
328  if (!ok) return std::nullopt;
329  }
330  else if (const SurfYZPlane* from = dynamic_cast<const SurfYZPlane*>(&*trk.getSurface())) {
331 
332  // Initial surface is SurfYZPlane.
333  // Get surface paramters.
334 
335  double phi1 = from->phi();
336 
337  // Transform track to origin surface.
338 
339  bool ok = transformYZPlane(phi1, theta2, phi2, vec, dir, prop_matrix);
340  result = std::make_optional(0.);
341  if (!ok) return std::nullopt;
342  }
343  else if (const SurfXYZPlane* from = dynamic_cast<const SurfXYZPlane*>(&*trk.getSurface())) {
344 
345  // Initial surface is SurfXYZPlane.
346  // Get surface paramters.
347 
348  double theta1 = from->theta();
349  double phi1 = from->phi();
350 
351  // Transform track to origin surface.
352 
353  bool ok = transformXYZPlane(theta1, phi1, theta2, phi2, vec, dir, prop_matrix);
354  result = std::make_optional(0.);
355  if (!ok) return std::nullopt;
356  }
357 
358  // Update track.
359 
360  trk.setSurface(porigin);
361  trk.setVector(vec);
362  trk.setDirection(dir);
363 
364  // Final validity check.
365 
366  if (!trk.isValid()) {
367  trk = trk0;
368  result = std::nullopt;
369  }
370 
371  // Done.
372 
373  return result;
374  }
375 
376  // Transform track parameters from SurfXYZLine to SurfXYZPlane.
377 
378  bool
380  double theta2,
381  double phi2,
382  TrackVector& vec,
384  TrackMatrix* prop_matrix) const
385  {
386  // Calculate surface transcendental functions.
387 
388  double sinth2 = std::sin(theta2);
389  double costh2 = std::cos(theta2);
390 
391  double sindphi = std::sin(phi2 - phi1);
392  double cosdphi = std::cos(phi2 - phi1);
393 
394  // Get the initial track parameters.
395 
396  double r1 = vec(0);
397  double phid1 = vec(2);
398  double eta1 = vec(3);
399 
400  // Calculate elements of rotation matrix from initial coordinate
401  // system to destination coordinte system.
402 
403  double ruu = costh2;
404  double ruv = sinth2 * sindphi;
405  double ruw = -sinth2 * cosdphi;
406 
407  double rvv = cosdphi;
408  double rvw = sindphi;
409 
410  double rwu = sinth2;
411  double rwv = -costh2 * sindphi;
412  double rww = costh2 * cosdphi;
413 
414  // Calculate track transcendental functions.
415 
416  double sinphid1 = std::sin(phid1);
417  double cosphid1 = std::cos(phid1);
418  double sh1 = 1. / std::cosh(eta1); // sech(eta1)
419  double th1 = std::tanh(eta1);
420 
421  // Calculate initial position in Cartesian coordinates.
422 
423  double u1 = -r1 * sinphid1;
424  double w1 = r1 * cosphid1;
425 
426  // Calculate direction in source coordinate system.
427 
428  double du1 = sh1 * cosphid1;
429  double dv1 = th1;
430  double dw1 = sh1 * sinphid1;
431  //double duw2 = std::hypot(du2, dw2);
432 
433  // Rotate direction to destination coordinate system.
434 
435  double du2 = ruu * du1 + ruv * dv1 + ruw * dw1;
436  double dv2 = rvv * dv1 + rvw * dw1;
437  double dw2 = rwu * du1 + rwv * dv1 + rww * dw1;
438 
439  // Calculate the track direction relative to the destination surface.
440  // The track direction comes from the sign of dw2 (=dw/ds).
441  // If dw2 is zero, the destionation surface is unreachable, return failure.
442 
443  if (dw2 > 0.)
444  dir = Surface::TrackDirection::FORWARD;
445  else if (dw2 < 0.)
446  dir = Surface::TrackDirection::BACKWARD;
447  else
448  return false;
449 
450  // Calculate final track slope track parameters.
451 
452  double dudw2 = du2 / dw2;
453  double dvdw2 = dv2 / dw2;
454 
455  // Update propagation matrix (if requested).
456 
457  if (prop_matrix != 0) {
458  TrackMatrix& pm = *prop_matrix;
459  pm.resize(vec.size(), vec.size(), false);
460 
461  // Calculate partial derivatives.
462 
463  // Partials of initial positions and directions wrt initial t.p.'s.
464 
465  double du1dr1 = -sinphid1;
466  double du1dphi1 = -w1;
467 
468  double dw1dr1 = cosphid1;
469  double dw1dphi1 = u1;
470 
471  double ddu1dphi1 = -sinphid1 * sh1;
472  double ddu1deta1 = -cosphid1 * sh1 * th1;
473 
474  double ddv1deta1 = sh1 * sh1;
475 
476  double ddw1dphi1 = cosphid1 * sh1;
477  double ddw1deta1 = -sinphid1 * sh1 * th1;
478 
479  // Rotate partials to destination coordinate system.
480 
481  double du2dr1 = ruu * du1dr1 + ruw * dw1dr1;
482  double dv2dr1 = rvw * dw1dr1;
483  double dw2dr1 = rwu * du1dr1 + rww * dw1dr1;
484 
485  double du2dv1 = ruv;
486  double dv2dv1 = rvv;
487  double dw2dv1 = rwv;
488 
489  double du2dphi1 = ruu * du1dphi1 + ruw * dw1dphi1;
490  double dv2dphi1 = rvw * dw1dphi1;
491  double dw2dphi1 = rwu * du1dphi1 + rww * dw1dphi1;
492 
493  double ddu2dphi1 = ruu * ddu1dphi1 + ruw * ddw1dphi1;
494  double ddv2dphi1 = rvw * ddw1dphi1;
495  double ddw2dphi1 = rwu * ddu1dphi1 + rww * ddw1dphi1;
496 
497  double ddu2deta1 = ruu * ddu1deta1 + ruv * ddv1deta1 + ruw * ddw1deta1;
498  double ddv2deta1 = rvv * ddv1deta1 + rvw * ddw1deta1;
499  double ddw2deta1 = rwu * ddu1deta1 + rwv * ddv1deta1 + rww * ddw1deta1;
500 
501  // Partials of final slope t.p. wrt final position and direction.
502 
503  double ddudw2ddu2 = 1. / dw2;
504  double ddudw2ddw2 = -dudw2 / dw2;
505 
506  double ddvdw2ddv2 = 1. / dw2;
507  double ddvdw2ddw2 = -dvdw2 / dw2;
508 
509  // Partials of final slope t.p. wrt initial t.p.
510 
511  double ddudw2dphi1 = ddudw2ddu2 * ddu2dphi1 + ddudw2ddw2 * ddw2dphi1;
512  double ddudw2deta1 = ddudw2ddu2 * ddu2deta1 + ddudw2ddw2 * ddw2deta1;
513 
514  double ddvdw2dphi1 = ddvdw2ddv2 * ddv2dphi1 + ddvdw2ddw2 * ddw2dphi1;
515  double ddvdw2deta1 = ddvdw2ddv2 * ddv2deta1 + ddvdw2ddw2 * ddw2deta1;
516 
517  // We still need to calculate the corretion due to the dependence of the
518  // propagation distance on the initial track parameters. This correction is
519  // needed even though the actual propagation distance is zero.
520 
521  // This correction effects the u and v track parameters.
522 
523  // Partials of perpendicular propagation distance wrt initial t.p.
524 
525  double dstdr1 = -dw2dr1;
526  double dstdv1 = -dw2dv1;
527  double dstdphi1 = -dw2dphi1;
528 
529  // Calculate correction to u and v parameter partials wrt initial t.p. due to path length.
530 
531  du2dr1 += dstdr1 * dudw2;
532  du2dv1 += dstdv1 * dudw2;
533  du2dphi1 += dstdphi1 * dudw2;
534 
535  dv2dr1 += dstdr1 * dvdw2;
536  dv2dv1 += dstdv1 * dvdw2;
537  dv2dphi1 += dstdphi1 * dvdw2;
538 
539  // Fill derivative matrix.
540 
541  pm(0, 0) = du2dr1; // du2/dr1
542  pm(1, 0) = dv2dr1; // dv2/dr1
543  pm(2, 0) = 0.; // d(dudw2)/dr1
544  pm(3, 0) = 0.; // d(dvdw2)/dr1
545  pm(4, 0) = 0.; // d(pinv2)/dr1
546 
547  pm(0, 1) = du2dv1; // du2/dv1
548  pm(1, 1) = dv2dv1; // dv2/dv1
549  pm(2, 1) = 0.; // d(dudw2)/dv1
550  pm(3, 1) = 0.; // d(dvdw2)/dv1
551  pm(4, 1) = 0.; // d(pinv2)/dv1
552 
553  pm(0, 2) = du2dphi1; // du2/d(phi1);
554  pm(1, 2) = dv2dphi1; // dv2/d(phi1);
555  pm(2, 2) = ddudw2dphi1; // d(dudw2)/d(phi1);
556  pm(3, 2) = ddvdw2dphi1; // d(dvdw2)/d(phi1);
557  pm(4, 2) = 0.; // d(pinv2)/d(phi1);
558 
559  pm(0, 3) = 0.; // du2/d(eta1);
560  pm(1, 3) = 0.; // dv2/d(eta1);
561  pm(2, 3) = ddudw2deta1; // d(dudw2)/d(eta1);
562  pm(3, 3) = ddvdw2deta1; // d(dvdw2)/d(eta1);
563  pm(4, 3) = 0.; // d(pinv2)/d(eta1);
564 
565  pm(0, 4) = 0.; // du2/d(pinv1);
566  pm(1, 4) = 0.; // dv2/d(pinv1);
567  pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
568  pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
569  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
570  }
571 
572  // Update track vector.
573 
574  vec(0) = 0.;
575  vec(1) = 0.;
576  vec(2) = dudw2;
577  vec(3) = dvdw2;
578 
579  // Done (success).
580 
581  return true;
582  }
583 
584  // Transform track parameters from SurfYZPlane to SurfXYZPlane.
585 
586  bool
588  double theta2,
589  double phi2,
590  TrackVector& vec,
592  TrackMatrix* prop_matrix) const
593  {
594  // Calculate transcendental functions.
595 
596  double sinth2 = std::sin(theta2);
597  double costh2 = std::cos(theta2);
598 
599  double sindphi = std::sin(phi2 - phi1);
600  double cosdphi = std::cos(phi2 - phi1);
601 
602  // Get the initial track state vector and track parameters.
603 
604  double dudw1 = vec(2);
605  double dvdw1 = vec(3);
606 
607  // Make sure initial track has a valid direction.
608 
609  if (dir == Surface::UNKNOWN) return false;
610 
611  // Calculate elements of rotation matrix from initial coordinate
612  // system to destination coordinte system.
613 
614  double ruu = costh2;
615  double ruv = sinth2 * sindphi;
616  double ruw = -sinth2 * cosdphi;
617 
618  double rvv = cosdphi;
619  double rvw = sindphi;
620 
621  double rwu = sinth2;
622  double rwv = -costh2 * sindphi;
623  double rww = costh2 * cosdphi;
624 
625  // Calculate the derivative dw2/dw1;
626  // If dw2/dw1 == 0., that means the track is moving parallel
627  // to destination plane.
628  // In this case return propagation failure.
629 
630  double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
631  if (dw2dw1 == 0.) return false;
632 
633  // Calculate slope in destination plane coordinates.
634 
635  double dudw2 = (dudw1 * ruu + dvdw1 * ruv + ruw) / dw2dw1;
636  double dvdw2 = (dvdw1 * rvv + rvw) / dw2dw1;
637 
638  // Calculate direction parameter at destination surface.
639  // Direction will flip if dw2dw1 < 0.;
640 
641  switch (dir) {
642  case Surface::FORWARD: dir = (dw2dw1 > 0.) ? Surface::FORWARD : Surface::BACKWARD; break;
643  case Surface::BACKWARD: dir = (dw2dw1 > 0.) ? Surface::BACKWARD : Surface::FORWARD; break;
644  default:
645  throw cet::exception("PropXYZPlane")
646  << __func__ << ": unexpected direction #" << ((int)dir) << "\n";
647  } // switch
648 
649  // Update propagation matrix (if requested).
650 
651  if (prop_matrix != 0) {
652  TrackMatrix& pm = *prop_matrix;
653  pm.resize(vec.size(), vec.size(), false);
654 
655  // Calculate partial derivatives.
656 
657  pm(0, 0) = ruu - dudw2 * rwu; // du2/du1
658  pm(1, 0) = -dvdw2 * rwu; // dv2/du1
659  pm(2, 0) = 0.; // d(dudw2)/du1
660  pm(3, 0) = 0.; // d(dvdw2)/du1
661  pm(4, 0) = 0.; // d(pinv2)/du1
662 
663  pm(0, 1) = ruv - dudw2 * rwv; // du2/dv1
664  pm(1, 1) = rvv - dvdw2 * rwv; // dv2/dv1
665  pm(2, 1) = 0.; // d(dudw2)/dv1
666  pm(3, 1) = 0.; // d(dvdw2)/dv1
667  pm(4, 1) = 0.; // d(pinv2)/dv1
668 
669  pm(0, 2) = 0.; // du2/d(dudw1);
670  pm(1, 2) = 0.; // dv2/d(dudw1);
671  pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1; // d(dudw2)/d(dudw1);
672  pm(3, 2) = -dvdw2 * rwu / dw2dw1; // d(dvdw2)/d(dudw1);
673  pm(4, 2) = 0.; // d(pinv2)/d(dudw1);
674 
675  pm(0, 3) = 0.; // du2/d(dvdw1);
676  pm(1, 3) = 0.; // dv2/d(dvdw1);
677  pm(2, 3) = (ruv - dudw2 * rwv) / dw2dw1; // d(dudw2)/d(dvdw1);
678  pm(3, 3) = (rvv - dvdw2 * rwv) / dw2dw1; // d(dvdw2)/d(dvdw1);
679  pm(4, 3) = 0.; // d(pinv2)/d(dvdw1);
680 
681  pm(0, 4) = 0.; // du2/d(pinv1);
682  pm(1, 4) = 0.; // dv2/d(pinv1);
683  pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
684  pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
685  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
686  }
687 
688  // Update track vector.
689 
690  vec(0) = 0.;
691  vec(1) = 0.;
692  vec(2) = dudw2;
693  vec(3) = dvdw2;
694 
695  // Done (success).
696 
697  return true;
698  }
699 
700  // Transform track parameters from SurfXYZPlane to SurfXYZPlane.
701 
702  bool
704  double phi1,
705  double theta2,
706  double phi2,
707  TrackVector& vec,
709  TrackMatrix* prop_matrix) const
710  {
711  // Calculate transcendental functions.
712 
713  double sinth1 = std::sin(theta1);
714  double costh1 = std::cos(theta1);
715  double sinth2 = std::sin(theta2);
716  double costh2 = std::cos(theta2);
717 
718  double sindphi = std::sin(phi2 - phi1);
719  double cosdphi = std::cos(phi2 - phi1);
720 
721  // Get the initial track state vector and track parameters.
722 
723  double dudw1 = vec(2);
724  double dvdw1 = vec(3);
725 
726  // Make sure initial track has a valid direction.
727 
728  if (dir == Surface::UNKNOWN) return false;
729 
730  // Calculate elements of rotation matrix from initial coordinate
731  // system to destination coordinte system.
732 
733  double ruu = costh1 * costh2 + sinth1 * sinth2 * cosdphi;
734  double ruv = sinth2 * sindphi;
735  double ruw = sinth1 * costh2 - costh1 * sinth2 * cosdphi;
736 
737  double rvu = -sinth1 * sindphi;
738  double rvv = cosdphi;
739  double rvw = costh1 * sindphi;
740 
741  double rwu = costh1 * sinth2 - sinth1 * costh2 * cosdphi;
742  double rwv = -costh2 * sindphi;
743  double rww = sinth1 * sinth2 + costh1 * costh2 * cosdphi;
744 
745  // Calculate the derivative dw2/dw1;
746  // If dw2/dw1 == 0., that means the track is moving parallel
747  // to destination plane.
748  // In this case return propagation failure.
749 
750  double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
751  if (dw2dw1 == 0.) return false;
752 
753  // Calculate slope in destination plane coordinates.
754 
755  double dudw2 = (dudw1 * ruu + dvdw1 * ruv + ruw) / dw2dw1;
756  double dvdw2 = (dudw1 * rvu + dvdw1 * rvv + rvw) / dw2dw1;
757 
758  // Calculate direction parameter at destination surface.
759  // Direction will flip if dw2dw1 < 0.;
760 
761  switch (dir) {
762  case Surface::FORWARD: dir = (dw2dw1 > 0.) ? Surface::FORWARD : Surface::BACKWARD; break;
763  case Surface::BACKWARD: dir = (dw2dw1 > 0.) ? Surface::BACKWARD : Surface::FORWARD; break;
764  default:
765  throw cet::exception("PropXYZPlane")
766  << __func__ << ": unexpected direction #" << ((int)dir) << "\n";
767  } // switch
768 
769  // Update propagation matrix (if requested).
770 
771  if (prop_matrix != 0) {
772  TrackMatrix& pm = *prop_matrix;
773  pm.resize(vec.size(), vec.size(), false);
774 
775  // Calculate partial derivatives.
776 
777  pm(0, 0) = ruu - dudw2 * rwu; // du2/du1
778  pm(1, 0) = rvu - dvdw2 * rwu; // dv2/du1
779  pm(2, 0) = 0.; // d(dudw2)/du1
780  pm(3, 0) = 0.; // d(dvdw2)/du1
781  pm(4, 0) = 0.; // d(pinv2)/du1
782 
783  pm(0, 1) = ruv - dudw2 * rwv; // du2/dv1
784  pm(1, 1) = rvv - dvdw2 * rwv; // dv2/dv1
785  pm(2, 1) = 0.; // d(dudw2)/dv1
786  pm(3, 1) = 0.; // d(dvdw2)/dv1
787  pm(4, 1) = 0.; // d(pinv2)/dv1
788 
789  pm(0, 2) = 0.; // du2/d(dudw1);
790  pm(1, 2) = 0.; // dv2/d(dudw1);
791  pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1; // d(dudw2)/d(dudw1);
792  pm(3, 2) = (rvu - dvdw2 * rwu) / dw2dw1; // d(dvdw2)/d(dudw1);
793  pm(4, 2) = 0.; // d(pinv2)/d(dudw1);
794 
795  pm(0, 3) = 0.; // du2/d(dvdw1);
796  pm(1, 3) = 0.; // dv2/d(dvdw1);
797  pm(2, 3) = (ruv - dudw2 * rwv) / dw2dw1; // d(dudw2)/d(dvdw1);
798  pm(3, 3) = (rvv - dvdw2 * rwv) / dw2dw1; // d(dvdw2)/d(dvdw1);
799  pm(4, 3) = 0.; // d(pinv2)/d(dvdw1);
800 
801  pm(0, 4) = 0.; // du2/d(pinv1);
802  pm(1, 4) = 0.; // dv2/d(pinv1);
803  pm(2, 4) = 0.; // d(dudw2)/d(pinv1);
804  pm(3, 4) = 0.; // d(dvdw2)/d(pinv1);
805  pm(4, 4) = 1.; // d(pinv2)/d(pinv1);
806  }
807 
808  // Update track vector.
809 
810  vec(0) = 0.;
811  vec(1) = 0.;
812  vec(2) = dudw2;
813  vec(3) = dvdw2;
814 
815  // Done (success).
816 
817  return true;
818  }
819 } // end namespace trkf
TrackDirection
Track direction enum.
Definition: Surface.h:56
double Mass() const
Based on pdg code.
Definition: KTrack.cxx:129
const std::shared_ptr< const Surface > & getSurface() const
Surface.
Definition: KTrack.h:55
static QCString result
General planar surface.
double y0() const
Y origin.
Definition: SurfXYZPlane.h:102
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.
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
void setVector(const TrackVector &vec)
Set state vector.
Definition: KTrack.h:67
void setDirection(Surface::TrackDirection dir)
Set direction.
Definition: KTrack.h:68
Planar surface parallel to x-axis.
const std::shared_ptr< const Interactor > & getInteractor() const
Definition: Propagator.h:118
string dir
PropXYZPlane(detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx)
Constructor.
void setSurface(const std::shared_ptr< const Surface > &psurf)
Set surface.
Definition: KTrack.h:66
bool transformXYZPlane(double theta1, double phi1, double theta2, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
Transform xyz plane -> xyz plane.
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.
void getPosition(double xyz[3]) const
Get position of track.
Definition: KTrack.cxx:171
Propagate to SurfXYZPlane surface.
double phi() const
Rot. angle about x-axis (wire angle).
Definition: SurfXYZPlane.h:104
double x0() const
X origin.
Definition: SurfXYZPlane.h:101
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.
double theta() const
Rot. angle about y&#39;-axis (projected Lorentz angle).
Definition: SurfXYZPlane.h:105
bool getDoDedx() const
Definition: Propagator.h:113
const TrackVector & getVector() const
Track state vector.
Definition: KTrack.h:56
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
Interactor for planar surfaces.
Surface::TrackDirection getDirection() const
Track direction.
Definition: KTrack.cxx:73
bool transformYZPlane(double phi1, double theta2, double phi2, TrackVector &vec, Surface::TrackDirection &dir, TrackMatrix *prop_matrix) const
Transform yz plane -> xyz plane.
double z0() const
Z origin.
Definition: SurfXYZPlane.h:103
Line surface perpendicular to x-axis.
static QCString * s
Definition: config.cpp:1042
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool isValid() const
Test if track is valid.
Definition: KTrack.cxx:91
PropDirection
Propagation direction enum.
Definition: Propagator.h:94