12 #include "cetlib_except/exception.h" 34 (tcut >= 0. ? std::make_shared<InteractPlane const>(detProp, tcut) :
35 std::shared_ptr<Interactor const>{})}
53 const std::shared_ptr<const Surface>& psurf,
61 std::optional<double> result{std::nullopt};
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();
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;
98 <<
"Track state vector has wrong size" << vec.size() <<
"\n";
101 double dudw1 = vec(2);
102 double dvdw1 = vec(3);
103 double pinv = vec(4);
115 double sinth2 = std::sin(theta2);
116 double costh2 = std::cos(theta2);
117 double sinphi2 = std::sin(phi2);
118 double cosphi2 = std::cos(phi2);
124 double ruy = sinth2 * sinphi2;
125 double ruz = -sinth2 * cosphi2;
127 double rvy = cosphi2;
128 double rvz = sinphi2;
131 double rwy = -costh2 * sinphi2;
132 double rwz = costh2 * cosphi2;
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;
143 double u2p = u2 - w2 * dudw1;
144 double v2p = v2 - w2 * dvdw1;
148 double s = -w2 * std::sqrt(1. + dudw1 * dudw1 + dvdw1 * dvdw1);
167 auto pinv2 = std::make_optional(pinv);
169 double* pderiv = (prop_matrix != 0 ? &deriv : 0);
182 result = std::make_optional(s);
186 if (prop_matrix != 0) {
188 pm.resize(vec.size(), vec.size(),
false);
225 *prop_matrix = prod(pm, *plocal_prop_matrix);
230 if (noise_matrix != 0) {
231 noise_matrix->resize(vec.size(), vec.size(),
false);
240 noise_matrix->clear();
276 const std::shared_ptr<const Surface>& porient,
281 std::optional<double> result{std::nullopt};
293 <<
"Track state vector has wrong size" << vec.size() <<
"\n";
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));
320 double phi1 = from->phi();
325 result = std::make_optional(0.);
326 if (!ok)
return std::nullopt;
333 double phi1 = from->phi();
338 result = std::make_optional(0.);
339 if (!ok)
return std::nullopt;
346 double theta1 = from->theta();
347 double phi1 = from->phi();
352 result = std::make_optional(0.);
353 if (!ok)
return std::nullopt;
366 result = std::nullopt;
385 double sinth2 = std::sin(theta2);
386 double costh2 = std::cos(theta2);
388 double sindphi = std::sin(phi2 - phi1);
389 double cosdphi = std::cos(phi2 - phi1);
394 double phid1 = vec(2);
395 double eta1 = vec(3);
401 double ruv = sinth2 * sindphi;
402 double ruw = -sinth2 * cosdphi;
404 double rvv = cosdphi;
405 double rvw = sindphi;
408 double rwv = -costh2 * sindphi;
409 double rww = costh2 * cosdphi;
413 double sinphid1 = std::sin(phid1);
414 double cosphid1 = std::cos(phid1);
415 double sh1 = 1. / std::cosh(eta1);
416 double th1 = std::tanh(eta1);
420 double u1 = -r1 * sinphid1;
421 double w1 = r1 * cosphid1;
425 double du1 = sh1 * cosphid1;
427 double dw1 = sh1 * sinphid1;
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;
441 dir = Surface::TrackDirection::FORWARD;
443 dir = Surface::TrackDirection::BACKWARD;
449 double dudw2 = du2 / dw2;
450 double dvdw2 = dv2 / dw2;
454 if (prop_matrix != 0) {
456 pm.resize(vec.size(), vec.size(),
false);
462 double du1dr1 = -sinphid1;
463 double du1dphi1 = -w1;
465 double dw1dr1 = cosphid1;
466 double dw1dphi1 = u1;
468 double ddu1dphi1 = -sinphid1 * sh1;
469 double ddu1deta1 = -cosphid1 * sh1 * th1;
471 double ddv1deta1 = sh1 * sh1;
473 double ddw1dphi1 = cosphid1 * sh1;
474 double ddw1deta1 = -sinphid1 * sh1 * th1;
478 double du2dr1 = ruu * du1dr1 + ruw * dw1dr1;
479 double dv2dr1 = rvw * dw1dr1;
480 double dw2dr1 = rwu * du1dr1 + rww * dw1dr1;
486 double du2dphi1 = ruu * du1dphi1 + ruw * dw1dphi1;
487 double dv2dphi1 = rvw * dw1dphi1;
488 double dw2dphi1 = rwu * du1dphi1 + rww * dw1dphi1;
490 double ddu2dphi1 = ruu * ddu1dphi1 + ruw * ddw1dphi1;
491 double ddv2dphi1 = rvw * ddw1dphi1;
492 double ddw2dphi1 = rwu * ddu1dphi1 + rww * ddw1dphi1;
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;
500 double ddudw2ddu2 = 1. / dw2;
501 double ddudw2ddw2 = -dudw2 / dw2;
503 double ddvdw2ddv2 = 1. / dw2;
504 double ddvdw2ddw2 = -dvdw2 / dw2;
508 double ddudw2dphi1 = ddudw2ddu2 * ddu2dphi1 + ddudw2ddw2 * ddw2dphi1;
509 double ddudw2deta1 = ddudw2ddu2 * ddu2deta1 + ddudw2ddw2 * ddw2deta1;
511 double ddvdw2dphi1 = ddvdw2ddv2 * ddv2dphi1 + ddvdw2ddw2 * ddw2dphi1;
512 double ddvdw2deta1 = ddvdw2ddv2 * ddv2deta1 + ddvdw2ddw2 * ddw2deta1;
522 double dstdr1 = -dw2dr1;
523 double dstdv1 = -dw2dv1;
524 double dstdphi1 = -dw2dphi1;
528 du2dr1 += dstdr1 * dudw2;
529 du2dv1 += dstdv1 * dudw2;
530 du2dphi1 += dstdphi1 * dudw2;
532 dv2dr1 += dstdr1 * dvdw2;
533 dv2dv1 += dstdv1 * dvdw2;
534 dv2dphi1 += dstdphi1 * dvdw2;
552 pm(2, 2) = ddudw2dphi1;
553 pm(3, 2) = ddvdw2dphi1;
558 pm(2, 3) = ddudw2deta1;
559 pm(3, 3) = ddvdw2deta1;
592 double sinth2 = std::sin(theta2);
593 double costh2 = std::cos(theta2);
595 double sindphi = std::sin(phi2 - phi1);
596 double cosdphi = std::cos(phi2 - phi1);
600 double dudw1 = vec(2);
601 double dvdw1 = vec(3);
611 double ruv = sinth2 * sindphi;
612 double ruw = -sinth2 * cosdphi;
614 double rvv = cosdphi;
615 double rvw = sindphi;
618 double rwv = -costh2 * sindphi;
619 double rww = costh2 * cosdphi;
626 double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
627 if (dw2dw1 == 0.)
return false;
631 double dudw2 = (dudw1 * ruu + dvdw1 * ruv + ruw) / dw2dw1;
632 double dvdw2 = (dvdw1 * rvv + rvw) / dw2dw1;
642 << __func__ <<
": unexpected direction #" << ((int)dir) <<
"\n";
647 if (prop_matrix != 0) {
649 pm.resize(vec.size(), vec.size(),
false);
653 pm(0, 0) = ruu - dudw2 * rwu;
654 pm(1, 0) = -dvdw2 * rwu;
659 pm(0, 1) = ruv - dudw2 * rwv;
660 pm(1, 1) = rvv - dvdw2 * rwv;
667 pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1;
668 pm(3, 2) = -dvdw2 * rwu / dw2dw1;
673 pm(2, 3) = (ruv - dudw2 * rwv) / dw2dw1;
674 pm(3, 3) = (rvv - dvdw2 * rwv) / dw2dw1;
708 double sinth1 = std::sin(theta1);
709 double costh1 = std::cos(theta1);
710 double sinth2 = std::sin(theta2);
711 double costh2 = std::cos(theta2);
713 double sindphi = std::sin(phi2 - phi1);
714 double cosdphi = std::cos(phi2 - phi1);
718 double dudw1 = vec(2);
719 double dvdw1 = vec(3);
728 double ruu = costh1 * costh2 + sinth1 * sinth2 * cosdphi;
729 double ruv = sinth2 * sindphi;
730 double ruw = sinth1 * costh2 - costh1 * sinth2 * cosdphi;
732 double rvu = -sinth1 * sindphi;
733 double rvv = cosdphi;
734 double rvw = costh1 * sindphi;
736 double rwu = costh1 * sinth2 - sinth1 * costh2 * cosdphi;
737 double rwv = -costh2 * sindphi;
738 double rww = sinth1 * sinth2 + costh1 * costh2 * cosdphi;
745 double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
746 if (dw2dw1 == 0.)
return false;
750 double dudw2 = (dudw1 * ruu + dvdw1 * ruv + ruw) / dw2dw1;
751 double dvdw2 = (dudw1 * rvu + dvdw1 * rvv + rvw) / dw2dw1;
761 << __func__ <<
": unexpected direction #" << ((int)dir) <<
"\n";
766 if (prop_matrix != 0) {
768 pm.resize(vec.size(), vec.size(),
false);
772 pm(0, 0) = ruu - dudw2 * rwu;
773 pm(1, 0) = rvu - dvdw2 * rwu;
778 pm(0, 1) = ruv - dudw2 * rwv;
779 pm(1, 1) = rvv - dvdw2 * rwv;
786 pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1;
787 pm(3, 2) = (rvu - dvdw2 * rwu) / dw2dw1;
792 pm(2, 3) = (ruv - dudw2 * rwv) / dw2dw1;
793 pm(3, 3) = (rvv - dvdw2 * rwv) / dw2dw1;
TrackDirection
Track direction enum.
double Mass() const
Based on pdg code.
const std::shared_ptr< const Surface > & getSurface() const
Surface.
double y0() const
Y origin.
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.
void setDirection(Surface::TrackDirection dir)
Set direction.
Planar surface parallel to x-axis.
const std::shared_ptr< const Interactor > & getInteractor() const
PropXYZPlane(detinfo::DetectorPropertiesData const &detProp, double tcut, bool doDedx)
Constructor.
void setSurface(const std::shared_ptr< const Surface > &psurf)
Set surface.
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.
Propagate to SurfXYZPlane surface.
double phi() const
Rot. angle about x-axis (wire angle).
double x0() const
X origin.
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'-axis (projected Lorentz angle).
const TrackVector & getVector() const
Track state vector.
std::optional< double > dedx_prop(double pinv, double mass, double s, double *deriv=0) const
Method to calculate updated momentum due to dE/dx.
Interactor for planar surfaces.
Surface::TrackDirection getDirection() const
Track direction.
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.
Line surface perpendicular to x-axis.
cet::coded_exception< error, detail::translate > exception
bool isValid() const
Test if track is valid.
PropDirection
Propagation direction enum.