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