1 #include <framework/dataobjects/Helix.h>
6 #include <framework/utilities/TestHelpers.h>
8 #include <gtest/gtest.h>
12 using namespace HelixParameterIndex;
15 std::ostream&
operator<<(::std::ostream& output,
const TVector3& tVector3)
19 << tVector3.X() <<
", "
20 << tVector3.Y() <<
", "
21 << tVector3.Z() <<
")";
35 bool d0Near = fabs(expected.getD0() - actual.getD0()) < tolerance;
36 bool phi0Near =
angleNear(expected.getPhi0(), actual.getPhi0(), tolerance);
37 bool omegaNear = fabs(expected.getOmega() - actual.getOmega()) < tolerance;
38 bool z0Near = fabs(expected.getZ0() - actual.getZ0()) < tolerance;
39 bool tanLambdaNear = fabs(expected.getTanLambda() - actual.getTanLambda()) < tolerance;
41 return d0Near and phi0Near and omegaNear and z0Near and tanLambdaNear;
46 Belle2::Helix helixFromCenter(
const TVector3& center,
const float& radius,
const float& tanLambda)
48 double omega = 1 / radius;
49 double phi0 = center.Phi() + copysign(M_PI / 2.0, radius);
50 double d0 = copysign(center.Perp(), radius) - radius;
51 double z0 = center.Z();
57 vector<float> linspace(
const float& start,
const float& end,
const int n)
59 std::vector<float> result(n);
63 for (
int i = 1; i < n - 1; ++i) {
64 float start_weight = (float)(n - 1 - i) / (n - 1);
65 float end_weight = 1 - start_weight;
66 result[i] = start * start_weight + end * end_weight;
73 class HelixTest :
public ::testing::Test {
77 double absError = 1
e-6;
78 double nominalBz = 1.5;
80 std::vector<float> omegas { -1, 0, 1};
82 std::vector<float> phi0s = linspace(-M_PI, M_PI, 11);
83 std::vector<float> d0s { -0.5, -0.2, 0, 0.2, 0.5};
85 std::vector<float> chis = linspace(-5 * M_PI / 6, 5 * M_PI / 6, 11);
87 std::vector<TVector3> bys = {
88 TVector3(-2.0, 0.5, 0.0),
89 TVector3(-1.0, 0.5, 0.0),
90 TVector3(0.0, 0.5, 0.0),
91 TVector3(1.0, 0.5, 0.0),
92 TVector3(2.0, 0.5, 0.0),
94 TVector3(-2.0, 0.0, 0.0),
95 TVector3(-1.0, 0.0, 0.0),
96 TVector3(0.0, 0.0, 0.0),
97 TVector3(1.0, 0.0, 0.0),
98 TVector3(2.0, 0.0, 0.0),
100 TVector3(-2.0, -1.0, 0.0),
101 TVector3(-1.0, -1.0, 0.0),
102 TVector3(0.0, -1.0, 0.0),
103 TVector3(1.0, -1.0, 0.0),
104 TVector3(2.0, -1.0, 0.0),
110 TEST_F(HelixTest, Getters)
113 unsigned int nCases = 1;
114 double bField = nominalBz;
116 for (
unsigned int i = 0; i < nCases; ++i) {
118 short int charge = generator.Uniform(-1, 1) > 0 ? 1 : -1;
121 TVector2 d(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
122 TVector2 pt(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
123 d.Set(d.X(), -d.X() * pt.Px() / pt.Py());
126 TVector3 position(d.X(), d.Y(), generator.Uniform(-1, 1));
127 TVector3 momentum(pt.Px(), pt.Py(), generator.Uniform(-1, 1));
130 Helix helix(position, momentum, charge, bField);
133 EXPECT_ALL_NEAR(position, helix.getPerigee(), absError);
134 EXPECT_ALL_NEAR(momentum, helix.getMomentum(bField), absError);
137 EXPECT_NEAR(momentum.Perp(), helix.getTransverseMomentum(bField), absError);
140 EXPECT_NEAR(charge / momentum.Perp(), helix.getKappa(bField), absError);
144 EXPECT_NEAR(position.Perp(), fabs(helix.getD0()), absError);
147 EXPECT_SAME_SIGN(position.Cross(momentum).Z(), helix.getD0());
150 EXPECT_NEAR(momentum.Phi(), helix.getPhi0(), absError);
153 EXPECT_NEAR(position.Z(), helix.getZ0(), absError);
156 EXPECT_NEAR(1 / tan(momentum.Theta()), helix.getTanLambda(), absError);
159 EXPECT_NEAR(1 / tan(momentum.Theta()), helix.getCotTheta(), absError);
162 EXPECT_EQ(charge, helix.getChargeSign());
167 TEST_F(HelixTest, SignOfD0)
171 const TVector3 position(1, 0, 0);
172 const TVector3 momentum(0, 1, 0);
173 const TVector3 oppositeMomentum(0, -1, 0);
175 const float bField = nominalBz;
177 Helix helix(position, momentum, charge, bField);
178 EXPECT_NEAR(1, helix.getD0(), absError);
181 Helix helix2(position, momentum, -charge, bField);
182 EXPECT_NEAR(1, helix2.getD0(), absError);
185 Helix oppositeMomentumHelix(position, oppositeMomentum, charge, bField);
186 EXPECT_NEAR(-1, oppositeMomentumHelix.getD0(), absError);
188 Helix oppositeMomentumHelix2(position, oppositeMomentum, -charge, bField);
189 EXPECT_NEAR(-1, oppositeMomentumHelix2.getD0(), absError);
201 TVector3 center(0.0, -2.0, 0.0);
206 Helix helix = helixFromCenter(center, radius, tanLambda);
208 TVector3 actualCenter = helix.getPerigee() * ((1 / helix.getOmega() + helix.getD0()) / helix.getD0());
209 EXPECT_NEAR(center.X(), actualCenter.X(), absError);
210 EXPECT_NEAR(center.Y(), actualCenter.Y(), absError);
213 TEST_F(HelixTest, Explicit)
221 TVector3 center(0.0, -2.0, 0.0);
226 Helix helix = helixFromCenter(center, radius, tanLambda);
227 EXPECT_NEAR(-1, helix.getD0(), absError);
228 EXPECT_ANGLE_NEAR(-M_PI, helix.getPhi0(), absError);
229 EXPECT_NEAR(-1, helix.getOmega(), absError);
234 float arcLength2D = 0;
235 TVector3 position = helix.getPositionAtArcLength2D(arcLength2D);
236 TVector3 tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
238 EXPECT_ALL_NEAR(TVector3(0.0, -1.0, 0.0), position, absError);
239 EXPECT_ANGLE_NEAR(-M_PI, tangential.Phi(), absError);
243 float arcLength2D = M_PI / 2;
244 TVector3 position = helix.getPositionAtArcLength2D(arcLength2D);
245 TVector3 tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
246 EXPECT_ALL_NEAR(TVector3(-1.0, -2.0, 0.0), position, absError);
247 EXPECT_ANGLE_NEAR(-M_PI / 2, tangential.Phi(), absError);
251 float arcLength2D = M_PI;
252 TVector3 position = helix.getPositionAtArcLength2D(arcLength2D);
253 TVector3 tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
254 EXPECT_ALL_NEAR(TVector3(0.0, -3.0, 0.0), position, absError);
255 EXPECT_ANGLE_NEAR(0, tangential.Phi(), absError);
259 float arcLength2D = 3 * M_PI / 2 ;
260 TVector3 position = helix.getPositionAtArcLength2D(arcLength2D);
261 TVector3 tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
262 EXPECT_ALL_NEAR(TVector3(1.0, -2.0, 0.0), position, absError);
263 EXPECT_ANGLE_NEAR(M_PI / 2, tangential.Phi(), absError);
267 TEST_F(HelixTest, Tangential)
272 for (
const float d0 : d0s) {
273 for (
const float phi0 : phi0s) {
274 for (
const float omega : omegas) {
276 Helix helix(d0, phi0, omega, z0, tanLambda);
277 TEST_CONTEXT(
"Failed for " << helix);
279 TVector3 tangentialAtPerigee = helix.getUnitTangentialAtArcLength2D(0.0);
281 EXPECT_ANGLE_NEAR(phi0, tangentialAtPerigee.Phi(), absError);
282 EXPECT_FLOAT_EQ(1.0, tangentialAtPerigee.Mag());
283 EXPECT_FLOAT_EQ(tanLambda, 1 / tan(tangentialAtPerigee.Theta()));
285 for (
const float chi : chis) {
289 float arcLength2D = chi;
292 TVector3 tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
293 EXPECT_ALL_NEAR(tangentialAtPerigee, tangential, absError);
296 float arcLength2D = -chi / omega;
297 TVector3 tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
299 float actualChi = tangential.DeltaPhi(tangentialAtPerigee);
300 EXPECT_ANGLE_NEAR(chi, actualChi, absError);
301 EXPECT_FLOAT_EQ(tangentialAtPerigee.Theta(), tangential.Theta());
302 EXPECT_FLOAT_EQ(1, tangential.Mag());
313 TEST_F(HelixTest, MomentumExtrapolation)
316 float tanLambda = -2;
318 for (
const float d0 : d0s) {
319 for (
const float phi0 : phi0s) {
320 for (
const float omega : omegas) {
323 Helix helix(d0, phi0, omega, z0, tanLambda);
324 TVector3 momentumAtPerigee = helix.getMomentum(nominalBz);
325 for (
const float chi : chis) {
327 float arcLength2D = -chi / omega;
328 TVector3 extrapolatedMomentum = helix.getMomentumAtArcLength2D(arcLength2D, nominalBz);
330 float actualChi = extrapolatedMomentum.DeltaPhi(momentumAtPerigee);
331 EXPECT_ANGLE_NEAR(chi, actualChi, absError);
332 EXPECT_FLOAT_EQ(momentumAtPerigee.Theta(), extrapolatedMomentum.Theta());
333 EXPECT_FLOAT_EQ(momentumAtPerigee.Mag(), extrapolatedMomentum.Mag());
343 TEST_F(HelixTest, Extrapolation)
350 for (
const float d0 : d0s) {
351 for (
const float phi0 : phi0s) {
352 for (
const float omega : omegas) {
354 Helix helix(d0, phi0, omega, z0, tanLambda);
355 TVector3 perigee = helix.getPerigee();
357 TVector3 tangentialAtPerigee = helix.getUnitTangentialAtArcLength2D(0.0);
358 TEST_CONTEXT(
"Failed for " << helix);
362 for (
const float chi : chis) {
363 TEST_CONTEXT(
"Failed for chi = " << chi);
367 float expectedArcLength2D = omega != 0 ? -chi / omega : chi;
368 TVector3 pointOnHelix = helix.getPositionAtArcLength2D(expectedArcLength2D);
370 float cylindricalR = pointOnHelix.Perp();
371 float arcLength2D = helix.getArcLength2DAtCylindricalR(cylindricalR);
374 EXPECT_NEAR(fabs(expectedArcLength2D), arcLength2D, absError);
377 TVector3 secantVector = pointOnHelix - perigee;
379 if (expectedArcLength2D == 0) {
380 EXPECT_NEAR(0, secantVector.Mag(), absError);
382 TVector2 secantVectorXY = secantVector.XYvector();
384 TVector2 tangentialXY = tangentialAtPerigee.XYvector();
385 float coalignment = secantVectorXY * tangentialXY ;
387 bool extrapolationIsForward = coalignment > 0;
388 bool expectedIsForward = expectedArcLength2D > 0;
389 EXPECT_EQ(expectedIsForward, extrapolationIsForward);
398 TEST_F(HelixTest, ExtrapolationToNormalPlane)
401 TVector3 center(0.0, 2.0, 0.0);
405 Helix helix = helixFromCenter(center, radius, tanLambda);
412 double arcLength2D = helix.getArcLength2DAtNormalPlane(x, y, nx, ny);
413 EXPECT_NEAR(-M_PI / 2, arcLength2D, absError);
417 TVector3 center(0.0, 2.0, 0.0);
421 Helix helix = helixFromCenter(center, radius, tanLambda);
428 double arcLength2D = helix.getArcLength2DAtNormalPlane(x, y, nx, ny);
429 EXPECT_NEAR(M_PI / 2, arcLength2D, absError);
433 TVector3 center(0.0, 2.0, 0.0);
437 Helix helix = helixFromCenter(center, radius, tanLambda);
444 double arcLength2D = helix.getArcLength2DAtNormalPlane(x, y, nx, ny);
445 EXPECT_NEAR(M_PI / 2, arcLength2D, absError);
449 TVector3 center(0.0, 2.0, 0.0);
453 Helix helix = helixFromCenter(center, radius, tanLambda);
460 double arcLength2D = helix.getArcLength2DAtNormalPlane(x, y, nx, ny);
461 EXPECT_NEAR(-M_PI / 2, arcLength2D, absError);
468 TEST_F(HelixTest, PerigeeExtrapolateRoundTrip)
471 float tanLambda = -2;
473 for (
const float d0 : d0s) {
474 for (
const float phi0 : phi0s) {
475 for (
const float omega : omegas) {
479 Helix expectedHelix(d0, phi0, omega, z0, tanLambda);
481 for (
const float chi : chis) {
482 float arcLength2D = -chi / omega;
483 TVector3 position = expectedHelix.getPositionAtArcLength2D(arcLength2D);
484 TVector3 momentum = expectedHelix.getMomentumAtArcLength2D(arcLength2D, nominalBz);
485 int chargeSign = expectedHelix.getChargeSign();
487 EXPECT_NEAR(tanLambda, 1 / tan(momentum.Theta()), absError);
488 EXPECT_ANGLE_NEAR(phi0 + chi, momentum.Phi(), absError);
489 EXPECT_NEAR(z0 + tanLambda * arcLength2D, position.Z(), absError);
492 Helix helix(position, momentum, chargeSign, nominalBz);
494 EXPECT_NEAR(expectedHelix.getOmega(), helix.getOmega(), absError);
495 EXPECT_ANGLE_NEAR(expectedHelix.getPhi0(), helix.getPhi0(), absError);
496 EXPECT_NEAR(expectedHelix.getD0(), helix.getD0(), absError);
497 EXPECT_NEAR(expectedHelix.getTanLambda(), helix.getTanLambda(), absError);
498 EXPECT_NEAR(expectedHelix.getZ0(), helix.getZ0(), absError);
592 TEST_F(HelixTest, PassiveMoveExplicit)
594 TVector3 center(0.0, 1.0, 0.0);
598 Helix helix = helixFromCenter(center, radius, tanLambda);
600 ASSERT_NEAR(0, helix.getD0(), absError);
601 ASSERT_ANGLE_NEAR(0, helix.getPhi0(), absError);
602 ASSERT_NEAR(-1, helix.getOmega(), absError);
605 Helix expectedHelix(helix);
609 TVector3 by(1.0, 1.0, 0.0);
611 float arcLength2D = helix.passiveMoveBy(by);
615 ASSERT_NEAR(M_PI / 2, arcLength2D, absError);
617 ASSERT_NEAR(0, helix.getD0(), absError);
618 ASSERT_ANGLE_NEAR(M_PI / 2, helix.getPhi0(), absError);
619 ASSERT_NEAR(-1, helix.getOmega(), absError);
621 ASSERT_NEAR(3 * M_PI / 2, helix.getZ0(), absError);
622 ASSERT_NEAR(3, helix.getTanLambda(), absError);
625 float arcLength2DBackward = helix.passiveMoveBy(-by);
627 ASSERT_NEAR(arcLength2D, -arcLength2DBackward, absError);
628 ASSERT_ALL_NEAR(expectedHelix, helix, absError);
632 TEST_F(HelixTest, PassiveMove)
637 for (
const float phi0 : phi0s) {
638 for (
const float omega : omegas) {
639 for (
const float d0 : d0s) {
640 for (
const float chi : chis) {
641 for (
const float newD0 : d0s) {
642 Helix helix(d0, phi0, omega, z0, tanLambda);
643 TEST_CONTEXT(
"Failed for " << helix);
647 float expectedArcLength2D = omega == 0 ? chi : -chi / omega;
648 TVector3 positionOnHelix = helix.getPositionAtArcLength2D(expectedArcLength2D);
649 TVector3 tangentialToHelix = helix.getUnitTangentialAtArcLength2D(expectedArcLength2D);
651 TVector3 perpendicularToHelix = tangentialToHelix;
652 perpendicularToHelix.RotateZ(M_PI / 2.0);
654 perpendicularToHelix *= 1 / perpendicularToHelix.Perp();
656 TVector3 displacementFromHelix = -perpendicularToHelix * newD0;
657 TVector3 testPosition = positionOnHelix + displacementFromHelix;
659 TVector3 expectedPerigee = -displacementFromHelix;
661 TEST_CONTEXT(
"Failed for chi " << chi);
662 TEST_CONTEXT(
"Failed for position on helix " << positionOnHelix);
663 TEST_CONTEXT(
"Failed for tangential to helix " << tangentialToHelix);
664 TEST_CONTEXT(
"Failed for perpendicular to helix " << perpendicularToHelix);
665 TEST_CONTEXT(
"Failed for test position " << testPosition);
667 float arcLength2D = helix.passiveMoveBy(testPosition);
669 ASSERT_NEAR(expectedArcLength2D, arcLength2D, absError);
671 ASSERT_ALL_NEAR(expectedPerigee, helix.getPerigee(), absError);
680 TEST_F(HelixTest, calcPassiveMoveByJacobian_identity)
682 TVector3 center(0.0, 1.0, 0.0);
686 Helix helix = helixFromCenter(center, radius, tanLambda);
688 TMatrixD noMoveJacobian = helix.calcPassiveMoveByJacobian(TVector3(0.0, 0.0, 0.0));
690 EXPECT_NEAR(1.0, noMoveJacobian(0, 0), 1e-7);
691 EXPECT_NEAR(0.0, noMoveJacobian(0, 1), 1e-7);
692 EXPECT_NEAR(0.0, noMoveJacobian(0, 2), 1e-7);
694 EXPECT_NEAR(0.0, noMoveJacobian(1, 0), 1e-7);
695 EXPECT_NEAR(1.0, noMoveJacobian(1, 1), 1e-7);
696 EXPECT_NEAR(0.0, noMoveJacobian(1, 2), 1e-7);
698 EXPECT_NEAR(0.0, noMoveJacobian(2, 0), 1e-7);
699 EXPECT_NEAR(0.0, noMoveJacobian(2, 1), 1e-7);
700 EXPECT_NEAR(1.0, noMoveJacobian(2, 2), 1e-7);
703 TEST_F(HelixTest, calcPassiveMoveByJacobian_orthogonalToPhi0)
705 TVector3 center(0.0, 1.0, 0.0);
709 Helix helix = helixFromCenter(center, radius, tanLambda);
711 TMatrixD moveByOneJacobian = helix.calcPassiveMoveByJacobian(TVector3(0.0, -1.0, 0.0));
713 EXPECT_NEAR(1.0, moveByOneJacobian(iD0, iD0), 1e-7);
714 EXPECT_NEAR(0.0, moveByOneJacobian(iD0, iPhi0), 1e-7);
715 EXPECT_NEAR(0.0, moveByOneJacobian(iD0, iOmega), 1e-7);
717 EXPECT_NEAR(0.0, moveByOneJacobian(iPhi0, iD0), 1e-7);
718 EXPECT_NEAR(1.0 / 2.0, moveByOneJacobian(iPhi0, iPhi0), 1e-7);
719 EXPECT_NEAR(0.0, moveByOneJacobian(iPhi0, iOmega), 1e-7);
721 EXPECT_NEAR(0.0, moveByOneJacobian(iOmega, iD0), 1e-7);
722 EXPECT_NEAR(0.0, moveByOneJacobian(iOmega, iPhi0), 1e-7);
723 EXPECT_NEAR(1.0, moveByOneJacobian(iOmega, iOmega), 1e-7);
727 TEST_F(HelixTest, calcPassiveMoveByJacobian_parallelToPhi0_compare_opposite_transformation)
729 TVector3 center(0.0, 1.0, 0.0);
733 Helix helix = helixFromCenter(center, radius, tanLambda);
735 TMatrixD moveByPlusTwoXJacobian = helix.calcPassiveMoveByJacobian(TVector3(2.0, 0.0, 0.0));
736 TMatrixD moveByMinusTwoXJacobian = helix.calcPassiveMoveByJacobian(TVector3(-2.0, 0.0, 0.0));
738 EXPECT_NEAR(1.0, moveByPlusTwoXJacobian(iOmega, iOmega), 1e-7);
739 EXPECT_NEAR(0.0, moveByPlusTwoXJacobian(iOmega, iPhi0), 1e-7);
740 EXPECT_NEAR(0.0, moveByPlusTwoXJacobian(iOmega, iD0), 1e-7);
742 EXPECT_NEAR(1.0, moveByMinusTwoXJacobian(iOmega, iOmega), 1e-7);
743 EXPECT_NEAR(0.0, moveByMinusTwoXJacobian(iOmega, iPhi0), 1e-7);
744 EXPECT_NEAR(0.0, moveByMinusTwoXJacobian(iOmega, iD0), 1e-7);
747 EXPECT_NEAR(moveByMinusTwoXJacobian(iD0, iOmega), moveByPlusTwoXJacobian(iD0, iOmega), 1e-7);
748 EXPECT_NEAR(moveByMinusTwoXJacobian(iPhi0, iPhi0), moveByPlusTwoXJacobian(iPhi0, iPhi0), 1e-7);
749 EXPECT_NEAR(moveByMinusTwoXJacobian(iD0, iD0), moveByPlusTwoXJacobian(iD0, iD0), 1e-7);
752 EXPECT_NEAR(moveByMinusTwoXJacobian(iD0, iPhi0), -moveByPlusTwoXJacobian(iD0, iPhi0), 1e-7);
753 EXPECT_NEAR(moveByMinusTwoXJacobian(iPhi0, iD0), -moveByPlusTwoXJacobian(iPhi0, iD0), 1e-7);
754 EXPECT_NEAR(moveByMinusTwoXJacobian(iPhi0, iOmega), -moveByPlusTwoXJacobian(iPhi0, iOmega), 1e-7);
757 EXPECT_NEGATIVE(moveByPlusTwoXJacobian(iD0, iPhi0));
758 EXPECT_POSITIVE(moveByPlusTwoXJacobian(iD0, iOmega));
759 EXPECT_POSITIVE(moveByPlusTwoXJacobian(iPhi0, iD0));
760 EXPECT_NEGATIVE(moveByPlusTwoXJacobian(iPhi0, iOmega));
762 EXPECT_NEGATIVE(moveByPlusTwoXJacobian(iD0, iD0) - 1);
763 EXPECT_NEGATIVE(moveByPlusTwoXJacobian(iPhi0, iPhi0) - 1);
766 TEST_F(HelixTest, calcPassiveMoveByJacobian_roundtrip)
768 for (
const TVector3& by : bys) {
769 TVector3 center(0.0, 1.0, 0.0);
773 Helix helix = helixFromCenter(center, radius, tanLambda);
775 TMatrixD jacobian = helix.calcPassiveMoveByJacobian(by);
776 helix.passiveMoveBy(by);
778 TMatrixD reverseJacobian = helix.calcPassiveMoveByJacobian(-by);
780 TMatrixD unitMatrix = reverseJacobian * jacobian;
782 for (
int i = 0; i < 5; ++i) {
783 for (
int j = 0; j < 5; ++j) {
786 EXPECT_NEAR(1, unitMatrix(i, j), 1e-7);
789 EXPECT_NEAR(0, unitMatrix(i, j), 1e-7);
798 TEST_F(HelixTest, calcPassiveMoveByJacobian_parallelToPhi0)
800 TVector3 center(1.0, 0.0, 0.0);
804 Helix helix = helixFromCenter(center, radius, tanLambda);
806 TMatrixD moveByTwoYJacobian = helix.calcPassiveMoveByJacobian(TVector3(0.0, 2.0, 0.0));
809 double deltaParallel = 2;
815 double lambda = 1.0 / (5.0 + 3.0 * sqrt(5.0));
816 double mu = sqrt(5.0) / 10.0;
819 EXPECT_NEAR(1.0, moveByTwoYJacobian(iOmega, iOmega), 1e-7);
820 EXPECT_NEAR(0.0, moveByTwoYJacobian(iOmega, iPhi0), 1e-7);
821 EXPECT_NEAR(0.0, moveByTwoYJacobian(iOmega, iD0), 1e-7);
823 EXPECT_NEAR(2.0 / 5.0, moveByTwoYJacobian(iPhi0, iOmega), 1e-7);
824 EXPECT_NEAR(1.0 / 5.0, moveByTwoYJacobian(iPhi0, iPhi0), 1e-7);
825 EXPECT_NEAR(-2.0 / 5.0, moveByTwoYJacobian(iPhi0, iD0), 1e-7);
827 EXPECT_NEAR(mu * zeta - A * lambda, moveByTwoYJacobian(iD0, iOmega), 1e-7);
828 EXPECT_NEAR(2.0 * mu * u * deltaParallel, moveByTwoYJacobian(iD0, iPhi0), 1e-7);
829 EXPECT_NEAR(2.0 * mu * nu, moveByTwoYJacobian(iD0, iD0), 1e-7);
834 TEST_F(HelixTest, calcPassiveMoveByJacobian_consistency_of_expansion)
836 TVector3 center(0.0, 1.0, 0.0);
840 Helix helix = helixFromCenter(center, radius, tanLambda);
842 for (
const TVector3& by : bys) {
843 TMatrixD jacobian = helix.calcPassiveMoveByJacobian(by);
844 TMatrixD jacobianWithExpansion = helix.calcPassiveMoveByJacobian(by, 10000);
846 EXPECT_NEAR(1.0, jacobian(iOmega, iOmega), 1e-7);
847 EXPECT_NEAR(0.0, jacobian(iOmega, iPhi0), 1e-7);
848 EXPECT_NEAR(0.0, jacobian(iOmega, iD0), 1e-7);
850 EXPECT_NEAR(1.0, jacobianWithExpansion(iOmega, iOmega), 1e-7);
851 EXPECT_NEAR(0.0, jacobianWithExpansion(iOmega, iPhi0), 1e-7);
852 EXPECT_NEAR(0.0, jacobianWithExpansion(iOmega, iD0), 1e-7);
854 EXPECT_NEAR(jacobianWithExpansion(iD0, iD0), jacobian(iD0, iD0), 1e-7);
855 EXPECT_NEAR(jacobianWithExpansion(iD0, iPhi0), jacobian(iD0, iPhi0), 1e-7);
856 EXPECT_NEAR(jacobianWithExpansion(iD0, iOmega), jacobian(iD0, iOmega), 1e-7);
858 EXPECT_NEAR(jacobianWithExpansion(iPhi0, iD0), jacobian(iPhi0, iD0), 1e-7);
859 EXPECT_NEAR(jacobianWithExpansion(iPhi0, iPhi0), jacobian(iPhi0, iPhi0), 1e-7);
860 EXPECT_NEAR(jacobianWithExpansion(iPhi0, iOmega), jacobian(iPhi0, iOmega), 1e-7);
864 TEST_F(HelixTest, realExample)
867 std::vector<double> helixParams(5);
868 helixParams[0] = 3.82384;
869 helixParams[1] = std::remainder(3.575595, 2 * M_PI);
870 helixParams[2] = 0.00530726;
871 helixParams[3] = -0.000317335;
872 helixParams[4] = 1.34536;
874 TVector3 momentum(-0.768755, -0.356297, 1.13994);
875 TVector3 position(-1.60794, 3.46933, -0.000317335);
881 const double bZ = 1.5;
884 EXPECT_NEAR(0.0, momentum.XYvector() * position.XYvector(), 1e-6);
886 Helix helix(helixParams[0], helixParams[1], helixParams[2], helixParams[3], helixParams[4]);
888 EXPECT_ALL_NEAR(momentum, helix.getMomentum(bZ), 1e-5);
889 EXPECT_ALL_NEAR(position, helix.getPerigee(), 1e-5);
893 TEST_F(HelixTest, omegaForUnitMomentum)
895 TVector3 expectedMomentum(1.0, 0.0, 0.0);
896 TVector3 expectedPosition(0.0, 1.0, 0.0);
897 const double expectedCharge = 1;
898 const double bZ = 1.5;
900 Helix helix(expectedPosition, expectedMomentum, expectedCharge, bZ);
902 double expectedOmega = 0.0044968868700000003;
903 EXPECT_ALL_NEAR(expectedOmega, helix.getOmega(), 1e-7);