Belle II Software prerelease-10-00-00a
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{
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 defined 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 defined 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 }
641
642}
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
Double32_t m_tanLambda
Memory for the slope of the track in the z coordinate over the two dimensional arc length (dz/ds)
Definition Helix.h:426
double getPerigeeX() const
Calculates the x coordinate of the perigee point.
Definition Helix.cc:54
double getSinPhi0() const
Getter for the cosine of the azimuth angle of travel direction at the perigee.
Definition Helix.h:384
double getMomentumZ(const double bZ) const
Calculates the z momentum of the particle at the perigee point.
Definition Helix.cc:84
Double32_t m_omega
Memory for the curvature of the signed curvature.
Definition Helix.h:420
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 getCosPhi0() const
Getter for the cosine of the azimuth angle of travel direction at the perigee.
Definition Helix.h:381
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
Double32_t m_z0
Memory for the z coordinate of the perigee.
Definition Helix.h:423
static double getAlpha(const double bZ)
Calculates the alpha value for a given magnetic field in Tesla.
Definition Helix.cc:109
double getOmega() const
Getter for omega, which is a signed curvature measure of the track.
Definition Helix.h:387
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
Double32_t m_d0
Memory for the signed distance to the perigee.
Definition Helix.h:414
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
double getD0() const
Getter for d0, which is the signed distance to the perigee in the r-phi plane.
Definition Helix.h:372
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
double getTanLambda() const
Getter for tan lambda, which is the z over two dimensional arc length slope of the track.
Definition Helix.h:393
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
Double32_t m_phi0
Memory for the azimuth angle between the transverse momentum and the x axis, which is in [-pi,...
Definition Helix.h:417
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 getZ0() const
Getter for z0, which is the z coordinate of the perigee.
Definition Helix.h:390
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
double getPhi0() const
Getter for phi0, which is the azimuth angle of the transverse momentum at the perigee.
Definition Helix.h:378
static double calcATanXDividedByX(const double &x)
Implementation of the function atan(x) / x which handles small x values smoothly.
Definition Helix.cc:477
std::ostream & operator<<(std::ostream &output, const IntervalOfValidity &iov)
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.