Belle II Software development
Helix.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include <framework/dataobjects/Helix.h>
9
10#include <framework/gearbox/Const.h>
11
12#include <boost/math/special_functions/sign.hpp>
13#include <boost/math/special_functions/sinc.hpp>
14#include <boost/math/tools/precision.hpp>
15
16#include <Math/VectorUtil.h>
17#include <TMatrixD.h>
18
19#include <cassert>
20
21using namespace Belle2;
22using namespace HelixParameterIndex;
23
25 m_d0(0.0),
26 m_phi0(0.0),
27 m_omega(0.0),
28 m_z0(0.0),
29 m_tanLambda(0.0)
30{
31}
32
33Helix::Helix(const ROOT::Math::XYZVector& position,
34 const ROOT::Math::XYZVector& momentum,
35 const short int charge,
36 const double bZ)
37{
38 setCartesian(position, momentum, charge, bZ);
39}
40
41Helix::Helix(const double& d0,
42 const double& phi0,
43 const double& omega,
44 const double& z0,
45 const double& tanLambda)
46 : m_d0(d0),
47 m_phi0(phi0),
48 m_omega(omega),
49 m_z0(z0),
50 m_tanLambda(tanLambda)
51{
52}
53
54double Helix::getPerigeeX() const
55{
56 return getD0() * getSinPhi0();
57}
58
59double Helix::getPerigeeY() const
60{
61 return -getD0() * getCosPhi0();
62}
63
64double Helix::getPerigeeZ() const
65{
66 return getZ0();
67}
68
69ROOT::Math::XYZVector Helix::getPerigee() const
70{
71 return ROOT::Math::XYZVector(getPerigeeX(), getPerigeeY(), getPerigeeZ());
72}
73
74double Helix::getMomentumX(const double bZ) const
75{
76 return getCosPhi0() * getTransverseMomentum(bZ);
77}
78
79double Helix::getMomentumY(const double bZ) const
80{
81 return getSinPhi0() * getTransverseMomentum(bZ);
82}
83
84double Helix::getMomentumZ(const double bZ) const
85{
86 return getTanLambda() * getTransverseMomentum(bZ);
87}
88
89ROOT::Math::XYZVector Helix:: getMomentum(const double bZ) const
90{
91 return ROOT::Math::XYZVector(getMomentumX(bZ), getMomentumY(bZ), getMomentumZ(bZ));
92}
93
94ROOT::Math::XYZVector Helix::getDirection() const
95{
96 return ROOT::Math::XYZVector(getCosPhi0(), getSinPhi0(), getTanLambda());
97}
98
99double Helix::getTransverseMomentum(const double bZ) const
100{
101 return 1 / std::fabs(getAlpha(bZ) * getOmega());
102}
103
104double Helix::getKappa(const double bZ) const
105{
106 return getOmega() * getAlpha(bZ);
107}
108
109double Helix::getAlpha(const double bZ)
110{
111 return 1.0 / (bZ * Const::speedOfLight) * 1E4;
112}
113
115{
116 return boost::math::sign(getOmega());
117}
118
119void Helix::reverse()
120{
121 // All except z0 have to be taken to their opposites
122 m_d0 = -m_d0; //d0
123 m_phi0 = reversePhi(m_phi0);
124 m_omega = -m_omega;
125 m_tanLambda = -m_tanLambda;
126}
127
128double Helix::getArcLength2DAtCylindricalR(const double& cylindricalR) const
129{
130 // Slight trick here
131 // Since the sought point is on the helix we treat it as the perigee
132 // and the origin as the point to extrapolate to.
133 // We know the distance of the origin to the circle, which is just d0
134 // The direct distance from the origin to the imaginary perigee is just the given cylindricalR.
135 const double dr = getD0();
136 const double deltaCylindricalR = cylindricalR;
137 const double absArcLength2D = calcArcLength2DAtDeltaCylindricalRAndDr(deltaCylindricalR, dr);
138 return absArcLength2D;
139}
140
141double Helix::getArcLength2DAtXY(const double& x, const double& y) const
142{
143 double dr = 0;
144 double arcLength2D = 0;
145 calcArcLength2DAndDrAtXY(x, y, arcLength2D, dr);
146 return arcLength2D;
147}
148
149double Helix::getArcLength2DAtNormalPlane(const double& byX, const double& byY,
150 const double& nX, const double& nY) const
151{
152 // Construct the tangential vector to the plan in xy space
153 const double tX = nY;
154 const double tY = -nX;
155
156 // Fetch the helix parameters
157 const double omega = getOmega();
158 const double cosPhi0 = getCosPhi0();
159 const double sinPhi0 = getSinPhi0();
160 const double d0 = getD0();
161
162 // Prepare a delta vector, which is the vector from the perigee point to the support point of the plane
163 // Split it in component parallel and a component orthogonal to tangent at the perigee.
164 const double deltaParallel = byX * cosPhi0 + byY * sinPhi0;
165 const double deltaOrthogonal = byY * cosPhi0 - byX * sinPhi0 + d0;
166 const double deltaCylindricalR = hypot(deltaOrthogonal, deltaParallel);
167 const double deltaCylindricalRSquared = deltaCylindricalR * deltaCylindricalR;
168
169 const double UOrthogonal = 1 + omega * deltaOrthogonal; // called nu in the Karimaki paper.
170 const double UParallel = omega * deltaParallel;
171
172 // Some commonly used terms - compare Karimaki 1990
173 const double A = 2 * deltaOrthogonal + omega * deltaCylindricalRSquared;
174
175 const double tParallel = tX * cosPhi0 + tY * sinPhi0;
176 const double tOrthogonal = tY * cosPhi0 - tX * sinPhi0;
177 const double tCylindricalR = (tX * tX + tY * tY);
178
179 const double c = A / 2;
180 const double b = UParallel * tParallel + UOrthogonal * tOrthogonal;
181 const double a = omega / 2 * tCylindricalR;
182
183 const double discriminant = ((double)b) * b - 4 * a * c;
184 const double root = sqrt(discriminant);
185 const double bigSum = (b > 0) ? -b - root : -b + root;
186
187 const double bigOffset = bigSum / 2 / a;
188 const double smallOffset = 2 * c / bigSum;
189
190 const double distance1 = hypot(byX + bigOffset * tX, byY + bigOffset * tY);
191 const double distance2 = hypot(byX + smallOffset * tX, byY + smallOffset * tY);
192
193 if (distance1 < distance2) {
194 return getArcLength2DAtXY(byX + bigOffset * tX, byY + bigOffset * tY);
195 } else {
196 return getArcLength2DAtXY(byX + smallOffset * tX, byY + smallOffset * tY);
197 }
198}
199
200
201ROOT::Math::XYZVector Helix::getPositionAtArcLength2D(const double& arcLength2D) const
202{
203 /*
204 / \ / \ / \
205 | x | | cos phi0 -sin phi0 | | - sin(chi) / omega |
206 | | = | | * | |
207 | y | | sin phi0 cos phi0 | | -(1 - cos(chi)) / omega - d0 |
208 \ / \ / \ /
209
210 and
211
212 z = tanLambda * arclength + z0;
213
214 where chi = -arcLength2D * omega
215
216 Here arcLength2D means the arc length of the circle in the xy projection
217 traversed in the forward direction.
218 */
219
220 // First calculating the position assuming the circle center lies on the y axes (phi0 = 0)
221 // Rotate to the right phi position afterwards
222 // Using the sinus cardinalis yields expressions that are smooth in the limit of omega -> 0
223
224 // Do calculations in double
225 const double chi = -arcLength2D * getOmega();
226 const double chiHalf = chi / 2.0;
227
228 using boost::math::sinc_pi;
229
230 const double x = arcLength2D * sinc_pi(chi);
231 const double y = arcLength2D * sinc_pi(chiHalf) * sin(chiHalf) - getD0();
232
233 // const double z = s * tan lambda + z0
234 const double z = fma((double)arcLength2D, getTanLambda(), getZ0());
235
236 // Unrotated position
237 ROOT::Math::XYZVector position(x, y, z);
238
239 // Rotate to the right phi0 position
240 position = ROOT::Math::VectorUtil::RotateZ(position, getPhi0());
241
242 return position;
243}
244
245ROOT::Math::XYZVector Helix::getTangentialAtArcLength2D(const double& arcLength2D) const
246{
247 const double omega = getOmega();
248 const double phi0 = getPhi0();
249 const double tanLambda = getTanLambda();
250
251 const double chi = - omega * arcLength2D;
252
253 const double tx = cos(chi + phi0);
254 const double ty = sin(chi + phi0);
255 const double tz = tanLambda;
256
257 ROOT::Math::XYZVector tangential(tx, ty, tz);
258 return tangential;
259}
260
261ROOT::Math::XYZVector Helix::getUnitTangentialAtArcLength2D(const double& arcLength2D) const
262{
263 ROOT::Math::XYZVector unitTangential = getTangentialAtArcLength2D(arcLength2D);
264 const double norm = hypot(1, getTanLambda());
265 const double invNorm = 1 / norm;
266 unitTangential *= invNorm;
267 return unitTangential;
268}
269
270ROOT::Math::XYZVector Helix::getMomentumAtArcLength2D(const double& arcLength2D, const double& bz) const
271{
272 ROOT::Math::XYZVector momentum = getTangentialAtArcLength2D(arcLength2D);
273 const double pr = getTransverseMomentum(bz);
274 momentum *= pr;
275
276 return momentum;
277}
278
279double Helix::passiveMoveBy(const double& byX,
280 const double& byY,
281 const double& byZ)
282{
283 // First calculate the distance of the new origin to the helix in the xy projection
284 double new_d0 = 0;
285 double arcLength2D = 0;
286 calcArcLength2DAndDrAtXY(byX, byY, arcLength2D, new_d0);
287
288 // Third the new phi0 and z0 can be calculated from the arc length
289 double chi = -arcLength2D * getOmega();
290 double new_phi0 = m_phi0 + chi;
291 double new_z0 = getZ0() - byZ + getTanLambda() * arcLength2D;
292
294 m_d0 = new_d0;
295 m_phi0 = new_phi0;
296 m_z0 = new_z0;
297
298 return arcLength2D;
299}
300
301TMatrixD Helix::calcPassiveMoveByJacobian(const ROOT::Math::XYZVector& by, const double expandBelowChi) const
302{
303 TMatrixD jacobian(5, 5);
304 calcPassiveMoveByJacobian(by.X(), by.Y(), jacobian, expandBelowChi);
305 return jacobian;
306}
307
308
309void Helix::calcPassiveMoveByJacobian(const double& byX,
310 const double& byY,
311 TMatrixD& jacobian,
312 const double expandBelowChi) const
313{
314 // 0. Preparations
315 // Initialise the return value to a unit matrix
316 jacobian.UnitMatrix();
317 assert(jacobian.GetNrows() == jacobian.GetNcols());
318 assert(jacobian.GetNrows() == 5 or jacobian.GetNrows() == 6);
319
320 // Fetch the helix parameters
321 const double omega = getOmega();
322 const double cosPhi0 = getCosPhi0();
323 const double sinPhi0 = getSinPhi0();
324 const double d0 = getD0();
325 const double tanLambda = getTanLambda();
326
327 // Prepare a delta vector, which is the vector from the perigee point to the new origin
328 // Split it in component parallel and a component orthogonal to tangent at the perigee.
329 const double deltaParallel = byX * cosPhi0 + byY * sinPhi0;
330 const double deltaOrthogonal = byY * cosPhi0 - byX * sinPhi0 + d0;
331 const double deltaCylindricalR = hypot(deltaOrthogonal, deltaParallel);
332 const double deltaCylindricalRSquared = deltaCylindricalR * deltaCylindricalR;
333
334 // Some commonly used terms - compare Karimaki 1990
335 const double A = 2 * deltaOrthogonal + omega * deltaCylindricalRSquared;
336 const double USquared = 1 + omega * A;
337 const double U = sqrt(USquared);
338 const double UOrthogonal = 1 + omega * deltaOrthogonal; // called nu in the Karimaki paper.
339 const double UParallel = omega * deltaParallel;
340 // Note U is a vector pointing from the middle of the projected circle scaled by a factor omega.
341 const double u = 1 + omega * d0; // just called u in the Karimaki paper.
342
343 // ---------------------------------------------------------------------------------------
344 // 1. Set the parts related to the xy coordinates
345 // Fills the upper left 3x3 corner of the jacobain
346
347 // a. Calculate the row related to omega
348 // The row related to omega is not different from the unit jacobian initialised above.
349 // jacobian(iOmega, iOmega) = 1.0; // From UnitMatrix above.
350
351 // b. Calculate the row related to d0
352 const double new_d0 = A / (1 + U);
353
354 jacobian(iD0, iOmega) = (deltaCylindricalR + new_d0) * (deltaCylindricalR - new_d0) / 2 / U;
355 jacobian(iD0, iPhi0) = - u * deltaParallel / U;
356 jacobian(iD0, iD0) = UOrthogonal / U;
357
358 // c. Calculate the row related to phi0
359 // Also calculate the derivatives of the arc length
360 // which are need for the row related to z.
361 const double dArcLength2D_dD0 = - UParallel / USquared;
362 const double dArcLength2D_dPhi0 = (omega * deltaCylindricalRSquared + deltaOrthogonal - d0 * UOrthogonal) / USquared;
363
364 jacobian(iPhi0, iD0) = - dArcLength2D_dD0 * omega;
365 jacobian(iPhi0, iPhi0) = u * UOrthogonal / USquared;
366 jacobian(iPhi0, iOmega) = -deltaParallel / USquared;
367
368 // For jacobian(iPhi0, iOmega) we have to dig deeper
369 // since the normal equations have a divergence for omega -> 0.
370 // This hinders not only the straight line case but extrapolations between points which are close together,
371 // like it happens when the current perigee is very close the new perigee.
372 // To have a smooth transition in this limit we have to carefully
373 // factor out the divergent quotients and approximate them with their taylor series.
374
375 const double chi = -atan2(UParallel, UOrthogonal);
376 double arcLength2D = 0;
377 double dArcLength2D_dOmega = 0;
378
379 if (fabs(chi) < std::min(expandBelowChi, M_PI / 2.0)) {
380 // Never expand for the far side of the circle by limiting the expandBelow to maximally half a pi.
381 // Close side of the circle
382 double principleArcLength2D = deltaParallel / UOrthogonal;
383 const double dPrincipleArcLength2D_dOmega = - principleArcLength2D * deltaOrthogonal / UOrthogonal;
384
385 const double x = principleArcLength2D * omega;
386 const double f = calcATanXDividedByX(x);
387 const double df_dx = calcDerivativeOfATanXDividedByX(x);
388
389 arcLength2D = principleArcLength2D * f;
390 dArcLength2D_dOmega = (f + x * df_dx) * dPrincipleArcLength2D_dOmega + principleArcLength2D * df_dx * principleArcLength2D;
391 } else {
392 // Far side of the circle
393 // If the far side of the circle is a well definied concept, omega is high enough that
394 // we can divide by it.
395 // Otherwise nothing can rescue us since the far side of the circle is so far away that no reasonable extrapolation can be made.
396 arcLength2D = - chi / omega;
397 dArcLength2D_dOmega = (-arcLength2D - jacobian(iPhi0, iOmega)) / omega;
398 }
399
400 // ---------------------------------------------------------------------------------------
401 // 2. Set the parts related to the z coordinate
402 // Since tanLambda stays the same there are no entries for tanLambda in the jacobian matrix
403 // For the new z0' = z0 + arcLength2D * tanLambda
404 jacobian(iZ0, iD0) = tanLambda * dArcLength2D_dD0;
405 jacobian(iZ0, iPhi0) = tanLambda * dArcLength2D_dPhi0;
406 jacobian(iZ0, iOmega) = tanLambda * dArcLength2D_dOmega;
407 jacobian(iZ0, iZ0) = 1.0; // From UnitMatrix above.
408 jacobian(iZ0, iTanLambda) = arcLength2D;
409
410 if (jacobian.GetNrows() == 6) {
411 // Also write the derivates of arcLength2D to the jacobian if the matrix size has space for them.
412 jacobian(iArcLength2D, iD0) = dArcLength2D_dD0;
413 jacobian(iArcLength2D, iPhi0) = dArcLength2D_dPhi0;
414 jacobian(iArcLength2D, iOmega) = dArcLength2D_dOmega;
415 jacobian(iArcLength2D, iZ0) = 0.0;
416 jacobian(iArcLength2D, iTanLambda) = 0;
417 jacobian(iArcLength2D, iArcLength2D) = 1.0;
418 }
419}
420
421double Helix::reversePhi(const double& phi)
422{
423 return std::remainder(phi + M_PI, 2 * M_PI);
424}
425
426double Helix::calcArcLength2DFromSecantLength(const double& secantLength) const
427{
428 return secantLength * calcSecantLengthToArcLength2DFactor(secantLength);
429}
430
431
432double Helix::calcSecantLengthToArcLength2DFactor(const double& secantLength) const
433{
434 double x = secantLength * getOmega() / 2.0;
435 return calcASinXDividedByX(x);
436}
437
438double Helix::calcASinXDividedByX(const double& x)
439{
440 // Approximation of asin(x) / x
441 // Inspired by BOOST's sinc
442
443 BOOST_MATH_STD_USING;
444
445 auto const taylor_n_bound = boost::math::tools::forth_root_epsilon<double>();
446
447 if (abs(x) >= taylor_n_bound) {
448 if (fabs(x) == 1) {
449 return M_PI / 2.0;
450
451 } else {
452 return asin(x) / x;
453
454 }
455
456 } else {
457 // approximation by taylor series in x at 0 up to order 0
458 double result = 1.0;
459
460 auto const taylor_0_bound = boost::math::tools::epsilon<double>();
461 if (abs(x) >= taylor_0_bound) {
462 double x2 = x * x;
463 // approximation by taylor series in x at 0 up to order 2
464 result += x2 / 6.0;
465
466 auto const taylor_2_bound = boost::math::tools::root_epsilon<double>();
467 if (abs(x) >= taylor_2_bound) {
468 // approximation by taylor series in x at 0 up to order 4
469 result += x2 * x2 * (3.0 / 40.0);
470 }
471 }
472 return result;
473 }
474
475}
476
477double Helix::calcATanXDividedByX(const double& x)
478{
479 // Approximation of atan(x) / x
480 // Inspired by BOOST's sinc
481
482 BOOST_MATH_STD_USING;
483
484 auto const taylor_n_bound = boost::math::tools::forth_root_epsilon<double>();
485
486 if (abs(x) >= taylor_n_bound) {
487 return atan(x) / x;
488
489 } else {
490 // approximation by taylor series in x at 0 up to order 0
491 double result = 1.0;
492
493 auto const taylor_0_bound = boost::math::tools::epsilon<double>();
494 if (abs(x) >= taylor_0_bound) {
495 double x2 = x * x;
496 // approximation by taylor series in x at 0 up to order 2
497 result -= x2 / 3.0;
498
499 auto const taylor_2_bound = boost::math::tools::root_epsilon<double>();
500 if (abs(x) >= taylor_2_bound) {
501 // approximation by taylor series in x at 0 up to order 4
502 result += x2 * x2 * (1.0 / 5.0);
503 }
504 }
505 return result;
506 }
507}
508
510{
511 // Approximation of atan(x) / x
512 // Inspired by BOOST's sinc
513
514 BOOST_MATH_STD_USING;
515
516 auto const taylor_n_bound = boost::math::tools::forth_root_epsilon<double>();
517
518 const double x2 = x * x;
519 if (abs(x) >= taylor_n_bound) {
520 const double chi = atan(x);
521 return ((1 - chi / x) / x - chi) / (1 + x2);
522
523 } else {
524 // approximation by taylor series in x at 0 up to order 0
525 double result = 1.0;
526
527 auto const taylor_0_bound = boost::math::tools::epsilon<double>();
528 if (abs(x) >= taylor_0_bound) {
529 // approximation by taylor series in x at 0 up to order 2
530 result -= 2.0 * x / 3.0;
531
532 auto const taylor_2_bound = boost::math::tools::root_epsilon<double>();
533 if (abs(x) >= taylor_2_bound) {
534 // approximation by taylor series in x at 0 up to order 4
535 result += x2 * x * (4.0 / 5.0);
536 }
537 }
538 return result;
539 }
540}
541
542
543
544
545void Helix::calcArcLength2DAndDrAtXY(const double& x, const double& y, double& arcLength2D, double& dr) const
546{
547 // Prepare common variables
548 const double omega = getOmega();
549 const double cosPhi0 = getCosPhi0();
550 const double sinPhi0 = getSinPhi0();
551 const double d0 = getD0();
552
553 const double deltaParallel = x * cosPhi0 + y * sinPhi0;
554 const double deltaOrthogonal = y * cosPhi0 - x * sinPhi0 + d0;
555 const double deltaCylindricalR = hypot(deltaOrthogonal, deltaParallel);
556 const double deltaCylindricalRSquared = deltaCylindricalR * deltaCylindricalR;
557
558 const double A = 2 * deltaOrthogonal + omega * deltaCylindricalRSquared;
559 const double U = sqrt(1 + omega * A);
560 const double UOrthogonal = 1 + omega * deltaOrthogonal; // called nu in the Karimaki paper.
561 const double UParallel = omega * deltaParallel;
562 // Note U is a vector pointing from the middle of the projected circle scaled by a factor omega.
563
564 // Calculate dr
565 dr = A / (1 + U);
566
567 // Calculate the absolute value of the arc length
568 const double chi = -atan2(UParallel, UOrthogonal);
569
570 if (fabs(chi) < M_PI / 8) { // Rough guess where the critical zone for approximations begins
571 // Close side of the circle
572 double principleArcLength2D = deltaParallel / UOrthogonal;
573 arcLength2D = principleArcLength2D * calcATanXDividedByX(principleArcLength2D * omega);
574 } else {
575 // Far side of the circle
576 // If the far side of the circle is a well definied concept meaning that we have big enough omega.
577 arcLength2D = -chi / omega;
578 }
579}
580
581double Helix::calcArcLength2DAtDeltaCylindricalRAndDr(const double& deltaCylindricalR, const double& dr) const
582{
583 const double omega = getOmega();
584 double secantLength = sqrt((deltaCylindricalR + dr) * (deltaCylindricalR - dr) / (1 + dr * omega));
585 return calcArcLength2DFromSecantLength(secantLength);
586}
587
588void Helix::setCartesian(const ROOT::Math::XYZVector& position,
589 const ROOT::Math::XYZVector& momentum,
590 const short int charge,
591 const double bZ)
592{
593 assert(abs(charge) <= 1); // Not prepared for doubly-charged particles.
594 const double alpha = getAlpha(bZ);
595
596 // We allow for the case that position, momentum are not given
597 // exactly in the perigee. Therefore we have to transform the momentum
598 // with the position as the reference point and then move the coordinate system
599 // to the origin.
600
601 const double x = position.X();
602 const double y = position.Y();
603 const double z = position.Z();
604
605 const double px = momentum.X();
606 const double py = momentum.Y();
607 const double pz = momentum.Z();
608
609 const double ptinv = 1 / hypot(px, py);
610 const double omega = charge * ptinv / alpha;
611 const double tanLambda = ptinv * pz;
612 const double phi0 = atan2(py, px);
613 const double z0 = z;
614 const double d0 = 0;
615
616 m_omega = omega;
617 m_phi0 = phi0;
618 m_d0 = d0;
619 m_tanLambda = tanLambda;
620 m_z0 = z0;
621
622 passiveMoveBy(-x, -y, 0);
623}
624
625namespace Belle2 {
631 std::ostream& operator<<(std::ostream& output, const Helix& helix)
632 {
633 return output
634 << "Helix("
635 << "d0=" << helix.getD0() << ", "
636 << "phi0=" << helix.getPhi0() << ", "
637 << "omega=" << helix.getOmega() << ", "
638 << "z0=" << helix.getZ0() << ", "
639 << "tanLambda=" << helix.getTanLambda() << ")";
640 }
642}
Helix parameter class.
Definition: Helix.h:48
double phi0(void) const
Return helix parameter phi0.
Definition: Helix.h:388
const HepVector & a(void) const
Returns helix parameters.
Definition: Helix.h:428
double dr(void) const
Return helix parameter dr.
Definition: Helix.h:380
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition: Helix.cc:259
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Helix.cc:197
double cosPhi0(void) const
Return cos phi0.
Definition: Helix.h:500
double sinPhi0(void) const
Return sin phi0.
Definition: Helix.h:492
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
double getPerigeeX() const
Calculates the x coordinate of the perigee point.
Definition: Helix.cc:54
double getMomentumZ(const double bZ) const
Calculates the z momentum of the particle at the perigee point.
Definition: Helix.cc:84
double getMomentumY(const double bZ) const
Calculates the y momentum of the particle at the perigee point.
Definition: Helix.cc:79
void reverse()
Reverses the direction of travel of the helix in place.
double getArcLength2DAtCylindricalR(const double &cylindricalR) const
Calculates the transverse travel distance at the point the helix first reaches the given cylindrical ...
Definition: Helix.cc:128
static double calcASinXDividedByX(const double &x)
Implementation of the function asin(x) / x which handles small x values smoothly.
Definition: Helix.cc:438
static double reversePhi(const double &phi)
Reverses an azimuthal angle to the opposite direction.
Definition: Helix.cc:421
short getChargeSign() const
Return track charge sign (1, 0 or -1).
Definition: Helix.cc:114
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 ...
Definition: Helix.cc:141
static double getAlpha(const double bZ)
Calculates the alpha value for a given magnetic field in Tesla.
Definition: Helix.cc:109
ROOT::Math::XYZVector getDirection() const
Getter for unit vector of momentum at the perigee position.
Definition: Helix.cc:94
ROOT::Math::XYZVector getTangentialAtArcLength2D(const double &arcLength2D) const
Calculates the tangential vector to the helix curve at the given two dimensional arc length.
Definition: Helix.cc:245
ROOT::Math::XYZVector getPositionAtArcLength2D(const double &arcLength2D) const
Calculates the position on the helix at the given two dimensional arc length.
Definition: Helix.cc:201
ROOT::Math::XYZVector getUnitTangentialAtArcLength2D(const double &arcLength2D) const
Calculates the unit tangential vector to the helix curve at the given two dimensional arc length.
Definition: Helix.cc:261
double getPerigeeZ() const
Calculates the z coordinate of the perigee point.
Definition: Helix.cc:64
void setCartesian(const ROOT::Math::XYZVector &position, const ROOT::Math::XYZVector &momentum, const short int charge, const double bZ)
Cartesian to Perigee conversion.
Definition: Helix.cc:588
ROOT::Math::XYZVector getMomentum(const double bZ) const
Getter for vector of momentum at the perigee position.
Definition: Helix.cc:89
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...
Definition: Helix.cc:545
double passiveMoveBy(const ROOT::Math::XYZVector &by)
Moves origin of the coordinate system (passive transformation) by the given vector.
Definition: Helix.h:245
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.
Definition: Helix.cc:149
double getMomentumX(const double bZ) const
Calculates the x momentum of the particle at the perigee point.
Definition: Helix.cc:74
TMatrixD calcPassiveMoveByJacobian(const ROOT::Math::XYZVector &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...
Definition: Helix.cc:301
double getPerigeeY() const
Calculates the y coordinate of the perigee point.
Definition: Helix.cc:59
double calcSecantLengthToArcLength2DFactor(const double &secantLength2D) const
Helper function to calculate the factor between the dimensional secant length and the two dimensional...
Definition: Helix.cc:432
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 ...
Definition: Helix.cc:581
double calcArcLength2DFromSecantLength(const double &secantLength2D) const
Helper function to calculate the two dimensional arc length from the length of a secant.
Definition: Helix.cc:426
double getKappa(const double bZ) const
Getter for kappa, which is charge / transverse momentum or equivalently omega * alpha.
Definition: Helix.cc:104
static double calcDerivativeOfATanXDividedByX(const double &x)
Implementation of the function d / dx (atan(x) / x) which handles small x values smoothly.
Definition: Helix.cc:509
Helix()
Constructor initializing all perigee parameters to zero.
ROOT::Math::XYZVector getMomentumAtArcLength2D(const double &arcLength2D, const double &bz) const
Calculates the momentum vector at the given two dimensional arc length.
Definition: Helix.cc:270
ROOT::Math::XYZVector getPerigee() const
Getter for the perigee position.
Definition: Helix.cc:69
double getTransverseMomentum(const double bZ) const
Getter for the absolute value of the transverse momentum at the perigee.
Definition: Helix.cc:99
static double calcATanXDividedByX(const double &x)
Implementation of the function atan(x) / x which handles small x values smoothly.
Definition: Helix.cc:477
double atan(double a)
atan for double
Definition: beamHelpers.h:34
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Abstract base class for different kinds of events.