12 #include "cetlib_except/exception.h" 34 (tcut >= 0. ? std::make_shared<InteractPlane const>(detProp, tcut) :
35 std::shared_ptr<Interactor const>{})}
54 const std::shared_ptr<const Surface>& psurf,
62 std::optional<double>
result{std::nullopt};
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();
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;
99 <<
"Track state vector has wrong size" << vec.size() <<
"\n";
102 double dudw1 = vec(2);
103 double dvdw1 = vec(3);
104 double pinv = vec(4);
116 double sinth2 = std::sin(theta2);
117 double costh2 = std::cos(theta2);
118 double sinphi2 = std::sin(phi2);
119 double cosphi2 = std::cos(phi2);
125 double ruy = sinth2 * sinphi2;
126 double ruz = -sinth2 * cosphi2;
128 double rvy = cosphi2;
129 double rvz = sinphi2;
132 double rwy = -costh2 * sinphi2;
133 double rwz = costh2 * cosphi2;
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;
144 double u2p = u2 - w2 * dudw1;
145 double v2p = v2 - w2 * dvdw1;
149 double s = -w2 * std::sqrt(1. + dudw1 * dudw1 + dvdw1 * dvdw1);
168 auto pinv2 = std::make_optional(pinv);
170 double* pderiv = (prop_matrix != 0 ? &deriv : 0);
183 result = std::make_optional(s);
187 if (prop_matrix != 0) {
189 pm.resize(vec.size(), vec.size(),
false);
226 *prop_matrix = prod(pm, *plocal_prop_matrix);
231 if (noise_matrix != 0) {
232 noise_matrix->resize(vec.size(), vec.size(),
false);
241 noise_matrix->clear();
276 std::optional<double>
278 const std::shared_ptr<const Surface>& porient,
283 std::optional<double>
result{std::nullopt};
295 <<
"Track state vector has wrong size" << vec.size() <<
"\n";
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));
322 double phi1 = from->phi();
327 result = std::make_optional(0.);
328 if (!ok)
return std::nullopt;
335 double phi1 = from->phi();
340 result = std::make_optional(0.);
341 if (!ok)
return std::nullopt;
348 double theta1 = from->theta();
349 double phi1 = from->phi();
354 result = std::make_optional(0.);
355 if (!ok)
return std::nullopt;
388 double sinth2 = std::sin(theta2);
389 double costh2 = std::cos(theta2);
391 double sindphi = std::sin(phi2 - phi1);
392 double cosdphi = std::cos(phi2 - phi1);
397 double phid1 = vec(2);
398 double eta1 = vec(3);
404 double ruv = sinth2 * sindphi;
405 double ruw = -sinth2 * cosdphi;
407 double rvv = cosdphi;
408 double rvw = sindphi;
411 double rwv = -costh2 * sindphi;
412 double rww = costh2 * cosdphi;
416 double sinphid1 = std::sin(phid1);
417 double cosphid1 = std::cos(phid1);
418 double sh1 = 1. / std::cosh(eta1);
419 double th1 = std::tanh(eta1);
423 double u1 = -r1 * sinphid1;
424 double w1 = r1 * cosphid1;
428 double du1 = sh1 * cosphid1;
430 double dw1 = sh1 * sinphid1;
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;
444 dir = Surface::TrackDirection::FORWARD;
446 dir = Surface::TrackDirection::BACKWARD;
452 double dudw2 = du2 / dw2;
453 double dvdw2 = dv2 / dw2;
457 if (prop_matrix != 0) {
459 pm.resize(vec.size(), vec.size(),
false);
465 double du1dr1 = -sinphid1;
466 double du1dphi1 = -w1;
468 double dw1dr1 = cosphid1;
469 double dw1dphi1 = u1;
471 double ddu1dphi1 = -sinphid1 * sh1;
472 double ddu1deta1 = -cosphid1 * sh1 * th1;
474 double ddv1deta1 = sh1 * sh1;
476 double ddw1dphi1 = cosphid1 * sh1;
477 double ddw1deta1 = -sinphid1 * sh1 * th1;
481 double du2dr1 = ruu * du1dr1 + ruw * dw1dr1;
482 double dv2dr1 = rvw * dw1dr1;
483 double dw2dr1 = rwu * du1dr1 + rww * dw1dr1;
489 double du2dphi1 = ruu * du1dphi1 + ruw * dw1dphi1;
490 double dv2dphi1 = rvw * dw1dphi1;
491 double dw2dphi1 = rwu * du1dphi1 + rww * dw1dphi1;
493 double ddu2dphi1 = ruu * ddu1dphi1 + ruw * ddw1dphi1;
494 double ddv2dphi1 = rvw * ddw1dphi1;
495 double ddw2dphi1 = rwu * ddu1dphi1 + rww * ddw1dphi1;
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;
503 double ddudw2ddu2 = 1. / dw2;
504 double ddudw2ddw2 = -dudw2 / dw2;
506 double ddvdw2ddv2 = 1. / dw2;
507 double ddvdw2ddw2 = -dvdw2 / dw2;
511 double ddudw2dphi1 = ddudw2ddu2 * ddu2dphi1 + ddudw2ddw2 * ddw2dphi1;
512 double ddudw2deta1 = ddudw2ddu2 * ddu2deta1 + ddudw2ddw2 * ddw2deta1;
514 double ddvdw2dphi1 = ddvdw2ddv2 * ddv2dphi1 + ddvdw2ddw2 * ddw2dphi1;
515 double ddvdw2deta1 = ddvdw2ddv2 * ddv2deta1 + ddvdw2ddw2 * ddw2deta1;
525 double dstdr1 = -dw2dr1;
526 double dstdv1 = -dw2dv1;
527 double dstdphi1 = -dw2dphi1;
531 du2dr1 += dstdr1 * dudw2;
532 du2dv1 += dstdv1 * dudw2;
533 du2dphi1 += dstdphi1 * dudw2;
535 dv2dr1 += dstdr1 * dvdw2;
536 dv2dv1 += dstdv1 * dvdw2;
537 dv2dphi1 += dstdphi1 * dvdw2;
555 pm(2, 2) = ddudw2dphi1;
556 pm(3, 2) = ddvdw2dphi1;
561 pm(2, 3) = ddudw2deta1;
562 pm(3, 3) = ddvdw2deta1;
596 double sinth2 = std::sin(theta2);
597 double costh2 = std::cos(theta2);
599 double sindphi = std::sin(phi2 - phi1);
600 double cosdphi = std::cos(phi2 - phi1);
604 double dudw1 = vec(2);
605 double dvdw1 = vec(3);
615 double ruv = sinth2 * sindphi;
616 double ruw = -sinth2 * cosdphi;
618 double rvv = cosdphi;
619 double rvw = sindphi;
622 double rwv = -costh2 * sindphi;
623 double rww = costh2 * cosdphi;
630 double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
631 if (dw2dw1 == 0.)
return false;
635 double dudw2 = (dudw1 * ruu + dvdw1 * ruv + ruw) / dw2dw1;
636 double dvdw2 = (dvdw1 * rvv + rvw) / dw2dw1;
646 << __func__ <<
": unexpected direction #" << ((
int)dir) <<
"\n";
651 if (prop_matrix != 0) {
653 pm.resize(vec.size(), vec.size(),
false);
657 pm(0, 0) = ruu - dudw2 * rwu;
658 pm(1, 0) = -dvdw2 * rwu;
663 pm(0, 1) = ruv - dudw2 * rwv;
664 pm(1, 1) = rvv - dvdw2 * rwv;
671 pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1;
672 pm(3, 2) = -dvdw2 * rwu / dw2dw1;
677 pm(2, 3) = (ruv - dudw2 * rwv) / dw2dw1;
678 pm(3, 3) = (rvv - dvdw2 * rwv) / dw2dw1;
713 double sinth1 = std::sin(theta1);
714 double costh1 = std::cos(theta1);
715 double sinth2 = std::sin(theta2);
716 double costh2 = std::cos(theta2);
718 double sindphi = std::sin(phi2 - phi1);
719 double cosdphi = std::cos(phi2 - phi1);
723 double dudw1 = vec(2);
724 double dvdw1 = vec(3);
733 double ruu = costh1 * costh2 + sinth1 * sinth2 * cosdphi;
734 double ruv = sinth2 * sindphi;
735 double ruw = sinth1 * costh2 - costh1 * sinth2 * cosdphi;
737 double rvu = -sinth1 * sindphi;
738 double rvv = cosdphi;
739 double rvw = costh1 * sindphi;
741 double rwu = costh1 * sinth2 - sinth1 * costh2 * cosdphi;
742 double rwv = -costh2 * sindphi;
743 double rww = sinth1 * sinth2 + costh1 * costh2 * cosdphi;
750 double dw2dw1 = dudw1 * rwu + dvdw1 * rwv + rww;
751 if (dw2dw1 == 0.)
return false;
755 double dudw2 = (dudw1 * ruu + dvdw1 * ruv + ruw) / dw2dw1;
756 double dvdw2 = (dudw1 * rvu + dvdw1 * rvv + rvw) / dw2dw1;
766 << __func__ <<
": unexpected direction #" << ((
int)dir) <<
"\n";
771 if (prop_matrix != 0) {
773 pm.resize(vec.size(), vec.size(),
false);
777 pm(0, 0) = ruu - dudw2 * rwu;
778 pm(1, 0) = rvu - dvdw2 * rwu;
783 pm(0, 1) = ruv - dudw2 * rwv;
784 pm(1, 1) = rvv - dvdw2 * rwv;
791 pm(2, 2) = (ruu - dudw2 * rwu) / dw2dw1;
792 pm(3, 2) = (rvu - dvdw2 * rwu) / dw2dw1;
797 pm(2, 3) = (ruv - dudw2 * rwv) / dw2dw1;
798 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.
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.
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.
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.