8 #include <framework/dataobjects/Helix.h>
10 #include <framework/gearbox/Const.h>
12 #include <boost/math/special_functions/sign.hpp>
13 #include <boost/math/special_functions/sinc.hpp>
14 #include <boost/math/tools/precision.hpp>
21 using namespace HelixParameterIndex;
33 const TVector3& momentum,
34 const short int charge,
44 const double& tanLambda)
49 m_tanLambda(tanLambda)
115 return boost::math::sign(
getOmega());
134 const double dr =
getD0();
135 const double deltaCylindricalR = cylindricalR;
137 return absArcLength2D;
143 double arcLength2D = 0;
149 const double& nX,
const double& nY)
const
152 const double tX = nY;
153 const double tY = -nX;
159 const double d0 =
getD0();
163 const double deltaParallel = byX * cosPhi0 + byY * sinPhi0;
164 const double deltaOrthogonal = byY * cosPhi0 - byX * sinPhi0 + d0;
165 const double deltaCylindricalR = hypot(deltaOrthogonal, deltaParallel);
166 const double deltaCylindricalRSquared = deltaCylindricalR * deltaCylindricalR;
168 const double UOrthogonal = 1 + omega * deltaOrthogonal;
169 const double UParallel = omega * deltaParallel;
172 const double A = 2 * deltaOrthogonal + omega * deltaCylindricalRSquared;
174 const double tParallel = tX * cosPhi0 + tY * sinPhi0;
175 const double tOrthogonal = tY * cosPhi0 - tX * sinPhi0;
176 const double tCylindricalR = (tX * tX + tY * tY);
178 const double c = A / 2;
179 const double b = UParallel * tParallel + UOrthogonal * tOrthogonal;
180 const double a = omega / 2 * tCylindricalR;
182 const double discriminant = ((double)b) * b - 4 * a * c;
183 const double root = sqrt(discriminant);
184 const double bigSum = (b > 0) ? -b - root : -b + root;
186 const double bigOffset = bigSum / 2 / a;
187 const double smallOffset = 2 * c / bigSum;
189 const double distance1 = hypot(byX + bigOffset * tX, byY + bigOffset * tY);
190 const double distance2 = hypot(byX + smallOffset * tX, byY + smallOffset * tY);
192 if (distance1 < distance2) {
224 const double chi = -arcLength2D *
getOmega();
225 const double chiHalf = chi / 2.0;
227 using boost::math::sinc_pi;
229 const double x = arcLength2D * sinc_pi(chi);
230 const double y = arcLength2D * sinc_pi(chiHalf) * sin(chiHalf) -
getD0();
236 TVector3 position(x, y, z);
250 const double chi = - omega * arcLength2D;
252 const double tx = cos(chi + phi0);
253 const double ty = sin(chi + phi0);
254 const double tz = tanLambda;
256 TVector3 tangential(tx, ty, tz);
264 const double invNorm = 1 / norm;
265 unitTangential *= invNorm;
266 return unitTangential;
284 double arcLength2D = 0;
288 double chi = -arcLength2D *
getOmega();
289 double new_phi0 =
m_phi0 + chi;
302 TMatrixD jacobian(5, 5);
311 const double expandBelowChi)
const
315 jacobian.UnitMatrix();
316 assert(jacobian.GetNrows() == jacobian.GetNcols());
317 assert(jacobian.GetNrows() == 5 or jacobian.GetNrows() == 6);
323 const double d0 =
getD0();
328 const double deltaParallel = byX * cosPhi0 + byY * sinPhi0;
329 const double deltaOrthogonal = byY * cosPhi0 - byX * sinPhi0 + d0;
330 const double deltaCylindricalR = hypot(deltaOrthogonal, deltaParallel);
331 const double deltaCylindricalRSquared = deltaCylindricalR * deltaCylindricalR;
334 const double A = 2 * deltaOrthogonal + omega * deltaCylindricalRSquared;
335 const double USquared = 1 + omega * A;
336 const double U = sqrt(USquared);
337 const double UOrthogonal = 1 + omega * deltaOrthogonal;
338 const double UParallel = omega * deltaParallel;
340 const double u = 1 + omega * d0;
351 const double new_d0 = A / (1 + U);
353 jacobian(iD0, iOmega) = (deltaCylindricalR + new_d0) * (deltaCylindricalR - new_d0) / 2 / U;
354 jacobian(iD0, iPhi0) = - u * deltaParallel / U;
355 jacobian(iD0, iD0) = UOrthogonal / U;
360 const double dArcLength2D_dD0 = - UParallel / USquared;
361 const double dArcLength2D_dPhi0 = (omega * deltaCylindricalRSquared + deltaOrthogonal - d0 * UOrthogonal) / USquared;
363 jacobian(iPhi0, iD0) = - dArcLength2D_dD0 * omega;
364 jacobian(iPhi0, iPhi0) = u * UOrthogonal / USquared;
365 jacobian(iPhi0, iOmega) = -deltaParallel / USquared;
374 const double chi = -atan2(UParallel, UOrthogonal);
375 double arcLength2D = 0;
376 double dArcLength2D_dOmega = 0;
378 if (fabs(chi) < std::min(expandBelowChi, M_PI / 2.0)) {
381 double principleArcLength2D = deltaParallel / UOrthogonal;
382 const double dPrincipleArcLength2D_dOmega = - principleArcLength2D * deltaOrthogonal / UOrthogonal;
384 const double x = principleArcLength2D * omega;
388 arcLength2D = principleArcLength2D * f;
389 dArcLength2D_dOmega = (f + x * df_dx) * dPrincipleArcLength2D_dOmega + principleArcLength2D * df_dx * principleArcLength2D;
395 arcLength2D = - chi / omega;
396 dArcLength2D_dOmega = (-arcLength2D - jacobian(iPhi0, iOmega)) / omega;
403 jacobian(iZ0, iD0) = tanLambda * dArcLength2D_dD0;
404 jacobian(iZ0, iPhi0) = tanLambda * dArcLength2D_dPhi0;
405 jacobian(iZ0, iOmega) = tanLambda * dArcLength2D_dOmega;
406 jacobian(iZ0, iZ0) = 1.0;
407 jacobian(iZ0, iTanLambda) = arcLength2D;
409 if (jacobian.GetNrows() == 6) {
411 jacobian(iArcLength2D, iD0) = dArcLength2D_dD0;
412 jacobian(iArcLength2D, iPhi0) = dArcLength2D_dPhi0;
413 jacobian(iArcLength2D, iOmega) = dArcLength2D_dOmega;
414 jacobian(iArcLength2D, iZ0) = 0.0;
415 jacobian(iArcLength2D, iTanLambda) = 0;
416 jacobian(iArcLength2D, iArcLength2D) = 1.0;
422 return std::remainder(phi + M_PI, 2 * M_PI);
433 double x = secantLength *
getOmega() / 2.0;
442 BOOST_MATH_STD_USING;
444 auto const taylor_n_bound = boost::math::tools::forth_root_epsilon<double>();
446 if (abs(x) >= taylor_n_bound) {
459 auto const taylor_0_bound = boost::math::tools::epsilon<double>();
460 if (abs(x) >= taylor_0_bound) {
465 auto const taylor_2_bound = boost::math::tools::root_epsilon<double>();
466 if (abs(x) >= taylor_2_bound) {
468 result += x2 * x2 * (3.0 / 40.0);
481 BOOST_MATH_STD_USING;
483 auto const taylor_n_bound = boost::math::tools::forth_root_epsilon<double>();
485 if (abs(x) >= taylor_n_bound) {
492 auto const taylor_0_bound = boost::math::tools::epsilon<double>();
493 if (abs(x) >= taylor_0_bound) {
498 auto const taylor_2_bound = boost::math::tools::root_epsilon<double>();
499 if (abs(x) >= taylor_2_bound) {
501 result += x2 * x2 * (1.0 / 5.0);
513 BOOST_MATH_STD_USING;
515 auto const taylor_n_bound = boost::math::tools::forth_root_epsilon<double>();
517 const double x2 = x * x;
518 if (abs(x) >= taylor_n_bound) {
519 const double chi = atan(x);
520 return ((1 - chi / x) / x - chi) / (1 + x2);
526 auto const taylor_0_bound = boost::math::tools::epsilon<double>();
527 if (abs(x) >= taylor_0_bound) {
529 result -= 2.0 * x / 3.0;
531 auto const taylor_2_bound = boost::math::tools::root_epsilon<double>();
532 if (abs(x) >= taylor_2_bound) {
534 result += x2 * x * (4.0 / 5.0);
550 const double d0 =
getD0();
552 const double deltaParallel = x * cosPhi0 + y * sinPhi0;
553 const double deltaOrthogonal = y * cosPhi0 - x * sinPhi0 + d0;
554 const double deltaCylindricalR = hypot(deltaOrthogonal, deltaParallel);
555 const double deltaCylindricalRSquared = deltaCylindricalR * deltaCylindricalR;
557 const double A = 2 * deltaOrthogonal + omega * deltaCylindricalRSquared;
558 const double U = sqrt(1 + omega * A);
559 const double UOrthogonal = 1 + omega * deltaOrthogonal;
560 const double UParallel = omega * deltaParallel;
567 const double chi = -atan2(UParallel, UOrthogonal);
569 if (fabs(chi) < M_PI / 8) {
571 double principleArcLength2D = deltaParallel / UOrthogonal;
576 arcLength2D = -chi / omega;
583 double secantLength = sqrt((deltaCylindricalR + dr) * (deltaCylindricalR - dr) / (1 + dr * omega));
588 const TVector3& momentum,
589 const short int charge,
592 assert(abs(charge) <= 1);
600 const double x = position.X();
601 const double y = position.Y();
602 const double z = position.Z();
604 const double px = momentum.X();
605 const double py = momentum.Y();
606 const double pz = momentum.Z();
608 const double ptinv = 1 / hypot(px, py);
609 const double omega = charge * ptinv / alpha;
610 const double tanLambda = ptinv * pz;
611 const double phi0 = atan2(py, px);
634 <<
"d0=" << helix.
getD0() <<
", "
635 <<
"phi0=" << helix.
getPhi0() <<
", "
636 <<
"omega=" << helix.
getOmega() <<
", "
637 <<
"z0=" << helix.
getZ0() <<
", "
static const double speedOfLight
[cm/ns]
This class represents an ideal helix in perigee parameterization.
TMatrixD calcPassiveMoveByJacobian(const TVector3 &by, const double expandBelowChi=M_PI/8) const
Calculate the 5x5 jacobian matrix for the transport of the helix parameters, when moving the origin o...
Double32_t m_tanLambda
Memory for the slope of the track in the z coordinate over the two dimensional arc length (dz/ds)
double getPerigeeX() const
Calculates the x coordinate of the perigee point.
double getSinPhi0() const
Getter for the cosine of the azimuth angle of travel direction at the perigee.
double getMomentumZ(const double bZ) const
Calculates the z momentum of the particle at the perigee point.
Double32_t m_omega
Memory for the curvature of the signed curvature.
double getMomentumY(const double bZ) const
Calculates the y momentum of the particle at the perigee point.
void reverse()
Reverses the direction of travel of the helix in place.
double getCosPhi0() const
Getter for the cosine of the azimuth angle of travel direction at the perigee.
double getArcLength2DAtCylindricalR(const double &cylindricalR) const
Calculates the transverse travel distance at the point the helix first reaches the given cylindrical ...
static double calcASinXDividedByX(const double &x)
Implementation of the function asin(x) / x which handles small x values smoothly.
static double reversePhi(const double &phi)
Reverses an azimuthal angle to the opposite direction.
short getChargeSign() const
Return track charge sign (1, 0 or -1).
TVector3 getTangentialAtArcLength2D(const double &arcLength2D) const
Calculates the tangential vector to the helix curve at the given two dimensional arc length.
double getArcLength2DAtXY(const double &x, const double &y) const
Calculates the two dimensional arc length at which the circle in the xy projection is closest to the ...
TVector3 getMomentumAtArcLength2D(const double &arcLength2D, const double &bz) const
Calculates the momentum vector at the given two dimensional arc length.
Double32_t m_z0
Memory for the z coordinate of the perigee.
static double getAlpha(const double bZ)
Calculates the alpha value for a given magnetic field in Tesla.
TVector3 getPerigee() const
Getter for the perigee position.
double getOmega() const
Getter for omega, which is a signed curvature measure of the track.
Double32_t m_d0
Memory for the signed distance to the perigee.
double getPerigeeZ() const
Calculates the z coordinate of the perigee point.
double getD0() const
Getter for d0, which is the signed distance to the perigee in the r-phi plane.
void calcArcLength2DAndDrAtXY(const double &x, const double &y, double &arcLength2D, double &dr) const
Helper method to calculate the signed two dimensional arc length and the signed distance to the circl...
double getArcLength2DAtNormalPlane(const double &x, const double &y, const double &nX, const double &nY) const
Calculates the arc length to reach a plane parallel to the z axes.
double getMomentumX(const double bZ) const
Calculates the x momentum of the particle at the perigee point.
double getTanLambda() const
Getter for tan lambda, which is the z over two dimensional arc length slope of the track.
TVector3 getPositionAtArcLength2D(const double &arcLength2D) const
Calculates the position on the helix at the given two dimensional arc length.
double getPerigeeY() const
Calculates the y coordinate of the perigee point.
double calcSecantLengthToArcLength2DFactor(const double &secantLength2D) const
Helper function to calculate the factor between the dimensional secant length and the two dimensional...
Double32_t m_phi0
Memory for the azimuth angle between the transverse momentum and the x axis, which is in [-pi,...
TVector3 getUnitTangentialAtArcLength2D(const double &arcLength2D) const
Calculates the unit tangential vector to the helix curve at the given two dimensional arc length.
double calcArcLength2DAtDeltaCylindricalRAndDr(const double &deltaCylindricalR, const double &dr) const
Helper method to calculate the two dimensional arc length from the perigee to a point at cylindrical ...
double calcArcLength2DFromSecantLength(const double &secantLength2D) const
Helper function to calculate the two dimensional arc length from the length of a secant.
double getZ0() const
Getter for z0, which is the z coordinate of the perigee.
double passiveMoveBy(const TVector3 &by)
Moves origin of the coordinate system (passive transformation) by the given vector.
double getKappa(const double bZ) const
Getter for kappa, which is charge / transverse momentum or equivalently omega * alpha.
static double calcDerivativeOfATanXDividedByX(const double &x)
Implementation of the function d / dx (atan(x) / x) which handles small x values smoothly.
Helix()
Constructor initializing all perigee parameters to zero.
double getTransverseMomentum(const double bZ) const
Getter for the absolute value of the transverse momentum at the perigee.
void setCartesian(const TVector3 &position, const TVector3 &momentum, const short int charge, const double bZ)
Cartesian to Perigee conversion.
TVector3 getMomentum(const double bZ) const
Getter for vector of momentum at the perigee position.
double getPhi0() const
Getter for phi0, which is the azimuth angle of the transverse momentum at the perigee.
static double calcATanXDividedByX(const double &x)
Implementation of the function atan(x) / x which handles small x values smoothly.
TVector3 getDirection() const
Getter for unit vector of momentum at the perigee position.
std::ostream & operator<<(std::ostream &output, const IntervalOfValidity &iov)
Abstract base class for different kinds of events.