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 <Math/VectorUtil.h>
11#include <Math/Vector2D.h>
12#include <Math/Vector3D.h>
13#include <TMatrixD.h>
14#include <TRandom3.h>
15
16#include <framework/utilities/TestHelpers.h>
17
18#include <gtest/gtest.h>
19
20using namespace std;
21using namespace Belle2;
22using namespace HelixParameterIndex;
23
25std::ostream& operator<<(::std::ostream& output, const ROOT::Math::XYZVector& vector3)
26{
27 return output
28 << "ROOT::Math::XYZVector("
29 << vector3.X() << ", "
30 << vector3.Y() << ", "
31 << vector3.Z() << ")";
32}
33
34namespace Belle2::TestHelpers {
35
40 template<>
41 bool allNear<Belle2::Helix>(const Belle2::Helix& expected,
42 const Belle2::Helix& actual,
43 double tolerance)
44 {
45 bool d0Near = fabs(expected.getD0() - actual.getD0()) < tolerance;
46 bool phi0Near = angleNear(expected.getPhi0(), actual.getPhi0(), tolerance);
47 bool omegaNear = fabs(expected.getOmega() - actual.getOmega()) < tolerance;
48 bool z0Near = fabs(expected.getZ0() - actual.getZ0()) < tolerance;
49 bool tanLambdaNear = fabs(expected.getTanLambda() - actual.getTanLambda()) < tolerance;
50
51 return d0Near and phi0Near and omegaNear and z0Near and tanLambdaNear;
52 }
53}
54
55namespace {
56 Belle2::Helix helixFromCenter(const ROOT::Math::XYZVector& center, const float& radius, const float& tanLambda)
57 {
58 double omega = 1 / radius;
59 double phi0 = center.Phi() + copysign(M_PI / 2.0, radius);
60 double d0 = copysign(center.Rho(), radius) - radius;
61 double z0 = center.Z();
62
63 return Belle2::Helix(d0, phi0, omega, z0, tanLambda);
64 }
65
67 vector<float> linspace(const float& start, const float& end, const int n)
68 {
69 std::vector<float> result(n);
70 result[0] = start;
71 result[n - 1] = end;
72
73 for (int i = 1; i < n - 1; ++i) {
74 float start_weight = (float)(n - 1 - i) / (n - 1);
75 float end_weight = 1 - start_weight;
76 result[i] = start * start_weight + end * end_weight;
77 }
78
79 return result;
80 }
81
83 class HelixTest : public ::testing::Test {
84
85 protected:
86 // Common level precision for all tests.
87 double absError = 1e-6;
88 double nominalBz = 1.5;
89
90 std::vector<float> omegas { -1, 0, 1};
91 //std::vector<float> omegas {1};
92 std::vector<float> phi0s = linspace(-M_PI, M_PI, 11);
93 std::vector<float> d0s { -0.5, -0.2, 0, 0.2, 0.5};
94 //std::vector<float> d0s {0.5};
95 std::vector<float> chis = linspace(-5 * M_PI / 6, 5 * M_PI / 6, 11);
96
97 std::vector<ROOT::Math::XYZVector> bys = {
98 ROOT::Math::XYZVector(-2.0, 0.5, 0.0),
99 ROOT::Math::XYZVector(-1.0, 0.5, 0.0),
100 ROOT::Math::XYZVector(0.0, 0.5, 0.0),
101 ROOT::Math::XYZVector(1.0, 0.5, 0.0),
102 ROOT::Math::XYZVector(2.0, 0.5, 0.0),
103
104 ROOT::Math::XYZVector(-2.0, 0.0, 0.0),
105 ROOT::Math::XYZVector(-1.0, 0.0, 0.0),
106 ROOT::Math::XYZVector(0.0, 0.0, 0.0),
107 ROOT::Math::XYZVector(1.0, 0.0, 0.0),
108 ROOT::Math::XYZVector(2.0, 0.0, 0.0),
109
110 ROOT::Math::XYZVector(-2.0, -1.0, 0.0),
111 ROOT::Math::XYZVector(-1.0, -1.0, 0.0),
112 ROOT::Math::XYZVector(0.0, -1.0, 0.0),
113 ROOT::Math::XYZVector(1.0, -1.0, 0.0),
114 ROOT::Math::XYZVector(2.0, -1.0, 0.0),
115 };
116
117 };
118
120 TEST_F(HelixTest, Getters)
121 {
122 TRandom3 generator;
123 unsigned int nCases = 1;
124 double bField = nominalBz;
125
126 for (unsigned int i = 0; i < nCases; ++i) {
127
128 short int charge = generator.Uniform(-1, 1) > 0 ? 1 : -1;
129
130 // Generate a random put orthogonal pair of vectors in the r-phi plane
131 ROOT::Math::XYVector d(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
132 ROOT::Math::XYVector pt(generator.Uniform(-1, 1), generator.Uniform(-1, 1));
133 d.SetCoordinates(d.X(), -d.X() * pt.X() / pt.Y());
134
135 // Add a random z component
136 ROOT::Math::XYZVector position(d.X(), d.Y(), generator.Uniform(-1, 1));
137 ROOT::Math::XYZVector momentum(pt.X(), pt.Y(), generator.Uniform(-1, 1));
138
139 // Set up class for testing
140 Helix helix(position, momentum, charge, bField);
141
142 // Test all vector elements
143 EXPECT_ALL_NEAR(position, helix.getPerigee(), absError);
144 EXPECT_ALL_NEAR(momentum, helix.getMomentum(bField), absError);
145
146 // Test getter for transverse momentum
147 EXPECT_NEAR(momentum.Rho(), helix.getTransverseMomentum(bField), absError);
148
149 // Test getter of kappa
150 EXPECT_NEAR(charge / momentum.Rho(), helix.getKappa(bField), absError);
151
152 // Test getter of d0
153 // Check absolute value of d0
154 EXPECT_NEAR(position.Rho(), fabs(helix.getD0()), absError);
155
156 // Check sign of d0
157 EXPECT_SAME_SIGN(position.Cross(momentum).Z(), helix.getD0());
158
159 // Test getter of phi0
160 EXPECT_NEAR(momentum.Phi(), helix.getPhi0(), absError);
161
162 // Test getter of d0
163 EXPECT_NEAR(position.Z(), helix.getZ0(), absError);
164
165 // Test getter of tan lambda
166 EXPECT_NEAR(1 / tan(momentum.Theta()), helix.getTanLambda(), absError);
167
168 // Test getter of cot theta
169 EXPECT_NEAR(1 / tan(momentum.Theta()), helix.getCotTheta(), absError);
170
171 // Test other variables
172 EXPECT_EQ(charge, helix.getChargeSign());
173
174 }
175 } // Testcases for getters
176
177 TEST_F(HelixTest, SignOfD0)
178 {
179 // This tests the assumption that the sign of d0 is given by the sign of position x momentum as a two dimensional cross product.
180
181 const ROOT::Math::XYZVector position(1, 0, 0);
182 const ROOT::Math::XYZVector momentum(0, 1, 0);
183 const ROOT::Math::XYZVector oppositeMomentum(0, -1, 0);
184 const float charge = 1;
185 const float bField = nominalBz;
186
187 Helix helix(position, momentum, charge, bField);
188 EXPECT_NEAR(1, helix.getD0(), absError);
189
190 // D0 does not change with the charge
191 Helix helix2(position, momentum, -charge, bField);
192 EXPECT_NEAR(1, helix2.getD0(), absError);
193
194 // But with reversal of momentum
195 Helix oppositeMomentumHelix(position, oppositeMomentum, charge, bField);
196 EXPECT_NEAR(-1, oppositeMomentumHelix.getD0(), absError);
197
198 Helix oppositeMomentumHelix2(position, oppositeMomentum, -charge, bField);
199 EXPECT_NEAR(-1, oppositeMomentumHelix2.getD0(), absError);
200
201 }
202
203 TEST_F(HelixTest, Center)
204 {
211 ROOT::Math::XYZVector center(0.0, -2.0, 0.0);
212 float radius = -1;
213 // Keep it flat
214 float tanLambda = 0;
215
216 Helix helix = helixFromCenter(center, radius, tanLambda);
217
218 ROOT::Math::XYZVector actualCenter = helix.getPerigee() * ((1 / helix.getOmega() + helix.getD0()) / helix.getD0());
219 EXPECT_NEAR(center.X(), actualCenter.X(), absError);
220 EXPECT_NEAR(center.Y(), actualCenter.Y(), absError);
221 }
222
223 TEST_F(HelixTest, Explicit)
224 {
231 ROOT::Math::XYZVector center(0.0, -2.0, 0.0);
232 float radius = -1;
233 // Keep it flat
234 float tanLambda = 0;
235
236 Helix helix = helixFromCenter(center, radius, tanLambda);
237 EXPECT_NEAR(-1, helix.getD0(), absError);
238 EXPECT_ANGLE_NEAR(-M_PI, helix.getPhi0(), absError);
239 EXPECT_NEAR(-1, helix.getOmega(), absError);
240
241 // Positions on the helix
242 {
243 // Start point
244 float arcLength2D = 0;
245 ROOT::Math::XYZVector position = helix.getPositionAtArcLength2D(arcLength2D);
246 ROOT::Math::XYZVector tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
247
248 EXPECT_ALL_NEAR(ROOT::Math::XYZVector(0.0, -1.0, 0.0), position, absError);
249 EXPECT_ANGLE_NEAR(-M_PI, tangential.Phi(), absError);
250 }
251
252 {
253 float arcLength2D = M_PI / 2;
254 ROOT::Math::XYZVector position = helix.getPositionAtArcLength2D(arcLength2D);
255 ROOT::Math::XYZVector tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
256 EXPECT_ALL_NEAR(ROOT::Math::XYZVector(-1.0, -2.0, 0.0), position, absError);
257 EXPECT_ANGLE_NEAR(-M_PI / 2, tangential.Phi(), absError);
258 }
259
260 {
261 float arcLength2D = M_PI;
262 ROOT::Math::XYZVector position = helix.getPositionAtArcLength2D(arcLength2D);
263 ROOT::Math::XYZVector tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
264 EXPECT_ALL_NEAR(ROOT::Math::XYZVector(0.0, -3.0, 0.0), position, absError);
265 EXPECT_ANGLE_NEAR(0, tangential.Phi(), absError);
266 }
267
268 {
269 float arcLength2D = 3 * M_PI / 2 ;
270 ROOT::Math::XYZVector position = helix.getPositionAtArcLength2D(arcLength2D);
271 ROOT::Math::XYZVector tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
272 EXPECT_ALL_NEAR(ROOT::Math::XYZVector(1.0, -2.0, 0.0), position, absError);
273 EXPECT_ANGLE_NEAR(M_PI / 2, tangential.Phi(), absError);
274 }
275 }
276
277 TEST_F(HelixTest, Tangential)
278 {
279 float z0 = 0;
280 float tanLambda = 2;
281
282 for (const float d0 : d0s) {
283 for (const float phi0 : phi0s) {
284 for (const float omega : omegas) {
285
286 Helix helix(d0, phi0, omega, z0, tanLambda);
287 TEST_CONTEXT("Failed for " << helix);
288
289 ROOT::Math::XYZVector tangentialAtPerigee = helix.getUnitTangentialAtArcLength2D(0.0);
290
291 EXPECT_ANGLE_NEAR(phi0, tangentialAtPerigee.Phi(), absError);
292 EXPECT_FLOAT_EQ(1.0, tangentialAtPerigee.R());
293 EXPECT_FLOAT_EQ(tanLambda, 1 / tan(tangentialAtPerigee.Theta()));
294
295 for (const float chi : chis) {
296
297 if (omega == 0) {
298 // Use chi as the arc length in the straight line case
299 float arcLength2D = chi;
300
301 // Tangential vector shall not change along the line
302 ROOT::Math::XYZVector tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
303 EXPECT_ALL_NEAR(tangentialAtPerigee, tangential, absError);
304
305 } else {
306 float arcLength2D = -chi / omega;
307 ROOT::Math::XYZVector tangential = helix.getUnitTangentialAtArcLength2D(arcLength2D);
308
309 float actualChi = ROOT::Math::VectorUtil::DeltaPhi(tangentialAtPerigee, tangential);
310 EXPECT_ANGLE_NEAR(chi, actualChi, absError);
311 EXPECT_FLOAT_EQ(tangentialAtPerigee.Theta(), tangential.Theta());
312 EXPECT_FLOAT_EQ(1, tangential.R());
313
314 }
315
316 }
317 }
318 }
319 }
320 }
321
322
323 TEST_F(HelixTest, MomentumExtrapolation)
324 {
325 float z0 = 0;
326 float tanLambda = -2;
327
328 for (const float d0 : d0s) {
329 for (const float phi0 : phi0s) {
330 for (const float omega : omegas) {
331 if (omega != 0) {
332
333 Helix helix(d0, phi0, omega, z0, tanLambda);
334 ROOT::Math::XYZVector momentumAtPerigee = helix.getMomentum(nominalBz);
335 for (const float chi : chis) {
336
337 float arcLength2D = -chi / omega;
338 ROOT::Math::XYZVector extrapolatedMomentum = helix.getMomentumAtArcLength2D(arcLength2D, nominalBz);
339
340 float actualChi = ROOT::Math::VectorUtil::DeltaPhi(momentumAtPerigee, extrapolatedMomentum);
341 EXPECT_ANGLE_NEAR(chi, actualChi, absError);
342 EXPECT_FLOAT_EQ(momentumAtPerigee.Theta(), extrapolatedMomentum.Theta());
343 EXPECT_FLOAT_EQ(momentumAtPerigee.R(), extrapolatedMomentum.R());
344 }
345 }
346 }
347 }
348 }
349 }
350
351
352
353 TEST_F(HelixTest, Extrapolation)
354 {
355
356 // z coordinates do not matter for this test.
357 float z0 = 0;
358 float tanLambda = 2;
359
360 for (const float d0 : d0s) {
361 for (const float phi0 : phi0s) {
362 for (const float omega : omegas) {
363
364 Helix helix(d0, phi0, omega, z0, tanLambda);
365 ROOT::Math::XYZVector perigee = helix.getPerigee();
366
367 ROOT::Math::XYZVector tangentialAtPerigee = helix.getUnitTangentialAtArcLength2D(0.0);
368 TEST_CONTEXT("Failed for " << helix);
369
370 //continue;
371
372 for (const float chi : chis) {
373 TEST_CONTEXT("Failed for chi = " << chi);
374
375 // In the cases where omega is 0 (straight line case) chi become undefined.
376 // Use chi sample as transverse travel distance instead.
377 float expectedArcLength2D = omega != 0 ? -chi / omega : chi;
378 ROOT::Math::XYZVector pointOnHelix = helix.getPositionAtArcLength2D(expectedArcLength2D);
379
380 float cylindricalR = pointOnHelix.Rho();
381 float arcLength2D = helix.getArcLength2DAtCylindricalR(cylindricalR);
382
383 // Only the absolute value is returned.
384 EXPECT_NEAR(fabs(expectedArcLength2D), arcLength2D, absError);
385
386 // Also check it the extrapolation lies in the forward direction.
387 ROOT::Math::XYZVector secantVector = pointOnHelix - perigee;
388
389 if (expectedArcLength2D == 0) {
390 EXPECT_NEAR(0, secantVector.R(), absError);
391 } else {
392 ROOT::Math::XYVector secantVectorXY(secantVector.X(), secantVector.Y());
393
394 ROOT::Math::XYVector tangentialXY(tangentialAtPerigee.X(), tangentialAtPerigee.Y());
395 float coalignment = secantVectorXY.Dot(tangentialXY);
396
397 bool extrapolationIsForward = coalignment > 0;
398 bool expectedIsForward = expectedArcLength2D > 0;
399 EXPECT_EQ(expectedIsForward, extrapolationIsForward);
400 }
401 }
402 }
403 }
404 }
405 } // end TEST_F
406
407
408 TEST_F(HelixTest, ExtrapolationToNormalPlane)
409 {
410 {
411 ROOT::Math::XYZVector center(0.0, 2.0, 0.0);
412 float radius = -1;
413 // Keep it flat
414 float tanLambda = 0;
415 Helix helix = helixFromCenter(center, radius, tanLambda);
416
417 // Plane coordinate
418 double x = 0.0;
419 double y = 3.0;
420 double nx = -1.0;
421 double ny = 1.0;
422 double arcLength2D = helix.getArcLength2DAtNormalPlane(x, y, nx, ny);
423 EXPECT_NEAR(-M_PI / 2, arcLength2D, absError);
424 }
425
426 {
427 ROOT::Math::XYZVector center(0.0, 2.0, 0.0);
428 float radius = -1;
429 // Keep it flat
430 float tanLambda = 0;
431 Helix helix = helixFromCenter(center, radius, tanLambda);
432
433 // Plane coordinate
434 double x = 0.0;
435 double y = 3.0;
436 double nx = 1.0;
437 double ny = 1.0;
438 double arcLength2D = helix.getArcLength2DAtNormalPlane(x, y, nx, ny);
439 EXPECT_NEAR(M_PI / 2, arcLength2D, absError);
440 }
441
442 {
443 ROOT::Math::XYZVector center(0.0, 2.0, 0.0);
444 float radius = 1;
445 // Keep it flat
446 float tanLambda = 0;
447 Helix helix = helixFromCenter(center, radius, tanLambda);
448
449 // Plane coordinate
450 double x = 0.0;
451 double y = 3.0;
452 double nx = -1.0;
453 double ny = 1.0;
454 double arcLength2D = helix.getArcLength2DAtNormalPlane(x, y, nx, ny);
455 EXPECT_NEAR(M_PI / 2, arcLength2D, absError);
456 }
457
458 {
459 ROOT::Math::XYZVector center(0.0, 2.0, 0.0);
460 float radius = 1;
461 // Keep it flat
462 float tanLambda = 0;
463 Helix helix = helixFromCenter(center, radius, tanLambda);
464
465 // Plane coordinate
466 double x = 0.0;
467 double y = 3.0;
468 double nx = 1.0;
469 double ny = 1.0;
470 double arcLength2D = helix.getArcLength2DAtNormalPlane(x, y, nx, ny);
471 EXPECT_NEAR(-M_PI / 2, arcLength2D, absError);
472 }
473
474
475
476 }
477
478 TEST_F(HelixTest, PerigeeExtrapolateRoundTrip)
479 {
480 float z0 = 0;
481 float tanLambda = -2;
482
483 for (const float d0 : d0s) {
484 for (const float phi0 : phi0s) {
485 for (const float omega : omegas) {
486
487 // Extrapolations involving the momentum only makes sense with finit momenta
488 if (omega != 0) {
489 Helix expectedHelix(d0, phi0, omega, z0, tanLambda);
490
491 for (const float chi : chis) {
492 float arcLength2D = -chi / omega;
493 ROOT::Math::XYZVector position = expectedHelix.getPositionAtArcLength2D(arcLength2D);
494 ROOT::Math::XYZVector momentum = expectedHelix.getMomentumAtArcLength2D(arcLength2D, nominalBz);
495 int chargeSign = expectedHelix.getChargeSign();
496
497 EXPECT_NEAR(tanLambda, 1 / tan(momentum.Theta()), absError);
498 EXPECT_ANGLE_NEAR(phi0 + chi, momentum.Phi(), absError);
499 EXPECT_NEAR(z0 + tanLambda * arcLength2D, position.Z(), absError);
500
501 //B2INFO("chi out " << chi);
502 Helix helix(position, momentum, chargeSign, nominalBz);
503
504 EXPECT_NEAR(expectedHelix.getOmega(), helix.getOmega(), absError);
505 EXPECT_ANGLE_NEAR(expectedHelix.getPhi0(), helix.getPhi0(), absError);
506 EXPECT_NEAR(expectedHelix.getD0(), helix.getD0(), absError);
507 EXPECT_NEAR(expectedHelix.getTanLambda(), helix.getTanLambda(), absError);
508 EXPECT_NEAR(expectedHelix.getZ0(), helix.getZ0(), absError);
509 }
510 }
511
512 }
513 }
514 }
515 }
516
517
518
519
520 /*
521 TEST_F(HelixTest, CalculateDrExplicit)
522 {
523 float tanLambda = 3;
524
525 ROOT::Math::XYZVector center(0.0, -2.0, 0.0);
526 float radius = -1;
527
528 Helix helix = helixFromCenter(center, radius, tanLambda);
529 EXPECT_NEAR(-1, helix.getD0(), absError);
530 EXPECT_ANGLE_NEAR(-M_PI, helix.getPhi0(), absError);
531 EXPECT_NEAR(-1, helix.getOmega(), absError);
532 {
533 ROOT::Math::XYZVector position(0.0, 0.0, 0.0);
534 float newD0 = helix.getDr(position);
535 EXPECT_NEAR(-1, newD0, absError);
536 }
537
538 {
539 ROOT::Math::XYZVector position(2.0, -2.0, 0.0);
540 float newD0 = helix.getDr(position);
541 EXPECT_NEAR(-1, newD0, absError);
542 }
543 {
544 ROOT::Math::XYZVector position(-2.0, -2.0, 0.0);
545 float newD0 = helix.getDr(position);
546 EXPECT_NEAR(-1, newD0, absError);
547 }
548
549 {
550 ROOT::Math::XYZVector position(1.0, -1.0, 0.0);
551 float newD0 = helix.getDr(position);
552 EXPECT_NEAR(-(sqrt(2) - 1) , newD0, absError);
553 }
554 }
555
556 TEST_F(HelixTest, CalculateDr)
557 {
558 float z0 = 2;
559 float tanLambda = 3;
560
561 for (const float phi0 : phi0s) {
562 for (const float omega : omegas) {
563 for (const float d0 : d0s) {
564 Helix helix(d0, phi0, omega, z0, tanLambda);
565 TEST_CONTEXT("Failed for " << helix);
566
567 EXPECT_NEAR(d0, helix.getDr(ROOT::Math::XYZVector(0.0, 0.0, 0.0)), absError);
568
569 for (const float chi : chis) {
570 for (const float newD0 : d0s) {
571 // In the line case use the chi value directly as the arc length
572
573 float arcLength2D = omega == 0 ? chi : -chi / omega;
574 ROOT::Math::XYZVector positionOnHelix = helix.getPositionAtArcLength2D(arcLength2D);
575
576 ROOT::Math::XYZVector tangentialToHelix = helix.getUnitTangentialAtArcLength2D(arcLength2D);
577
578 ROOT::Math::XYZVector perpendicularToHelix = tangentialToHelix;
579 perpendicularToHelix = ROOT::Math::VectorUtil::RotateZ(perpendicularToHelix, M_PI / 2.0);
580 // Normalize the xy part
581 perpendicularToHelix *= 1 / perpendicularToHelix.Perp();
582
583 ROOT::Math::XYZVector displacementFromHelix = perpendicularToHelix * newD0;
584 ROOT::Math::XYZVector testPosition = positionOnHelix + displacementFromHelix;
585
586 TEST_CONTEXT("Failed for chi " << chi);
587 TEST_CONTEXT("Failed for position on helix " << positionOnHelix);
588 TEST_CONTEXT("Failed for tangential to helix " << tangentialToHelix);
589 TEST_CONTEXT("Failed for perpendicular to helix " << perpendicularToHelix);
590 TEST_CONTEXT("Failed for test position " << testPosition);
591
592 float testDr = helix.getDr(testPosition);
593 EXPECT_NEAR(newD0, testDr, absError);
594 }
595 }
596 }
597 }
598 }
599 }
600 */
601
602 TEST_F(HelixTest, PassiveMoveExplicit)
603 {
604 ROOT::Math::XYZVector center(0.0, 1.0, 0.0);
605 float radius = -1;
606 float tanLambda = 3;
607
608 Helix helix = helixFromCenter(center, radius, tanLambda);
609
610 ASSERT_NEAR(0, helix.getD0(), absError);
611 ASSERT_ANGLE_NEAR(0, helix.getPhi0(), absError);
612 ASSERT_NEAR(-1, helix.getOmega(), absError);
613
614 // Save the untransformed Helix
615 Helix expectedHelix(helix);
616
617 // Vector by which the coordinate system should move.
618 // (To the right of the circle)
619 ROOT::Math::XYZVector by(1.0, 1.0, 0.0);
620
621 float arcLength2D = helix.passiveMoveBy(by);
622
623 // The right of the circle lies in the counterclockwise direction
624 // The forward direction is counterclockwise, so we expect to move forward.
625 ASSERT_NEAR(M_PI / 2, arcLength2D, absError);
626
627 ASSERT_NEAR(0, helix.getD0(), absError);
628 ASSERT_ANGLE_NEAR(M_PI / 2, helix.getPhi0(), absError);
629 ASSERT_NEAR(-1, helix.getOmega(), absError);
630
631 ASSERT_NEAR(3 * M_PI / 2, helix.getZ0(), absError);
632 ASSERT_NEAR(3, helix.getTanLambda(), absError);
633
634 // Now transform back to the original point
635 float arcLength2DBackward = helix.passiveMoveBy(-by);
636
637 ASSERT_NEAR(arcLength2D, -arcLength2DBackward, absError);
638 ASSERT_ALL_NEAR(expectedHelix, helix, absError);
639 }
640
641
642 TEST_F(HelixTest, PassiveMove)
643 {
644 float z0 = 2;
645 float tanLambda = 2;
646
647 for (const float phi0 : phi0s) {
648 for (const float omega : omegas) {
649 for (const float d0 : d0s) {
650 for (const float chi : chis) {
651 for (const float newD0 : d0s) {
652 Helix helix(d0, phi0, omega, z0, tanLambda);
653 TEST_CONTEXT("Failed for " << helix);
654
655 // In the line case use the chi value directly as the arc length
656
657 float expectedArcLength2D = omega == 0 ? chi : -chi / omega;
658 ROOT::Math::XYZVector positionOnHelix = helix.getPositionAtArcLength2D(expectedArcLength2D);
659 ROOT::Math::XYZVector tangentialToHelix = helix.getUnitTangentialAtArcLength2D(expectedArcLength2D);
660
661 ROOT::Math::XYZVector perpendicularToHelix = tangentialToHelix;
662 perpendicularToHelix = ROOT::Math::VectorUtil::RotateZ(perpendicularToHelix, M_PI / 2.0);
663 // Normalize the xy part
664 perpendicularToHelix *= 1 / perpendicularToHelix.Rho();
665
666 ROOT::Math::XYZVector displacementFromHelix = -perpendicularToHelix * newD0;
667 ROOT::Math::XYZVector testPosition = positionOnHelix + displacementFromHelix;
668
669 ROOT::Math::XYZVector expectedPerigee = -displacementFromHelix;
670
671 TEST_CONTEXT("Failed for chi " << chi);
672 TEST_CONTEXT("Failed for position on helix " << positionOnHelix);
673 TEST_CONTEXT("Failed for tangential to helix " << tangentialToHelix);
674 TEST_CONTEXT("Failed for perpendicular to helix " << perpendicularToHelix);
675 TEST_CONTEXT("Failed for test position " << testPosition);
676
677 float arcLength2D = helix.passiveMoveBy(testPosition);
678
679 ASSERT_NEAR(expectedArcLength2D, arcLength2D, absError);
680
681 ASSERT_ALL_NEAR(expectedPerigee, helix.getPerigee(), absError);
682
683 }
684 }
685 }
686 }
687 }
688 }
689
690 TEST_F(HelixTest, calcPassiveMoveByJacobian_identity)
691 {
692 ROOT::Math::XYZVector center(0.0, 1.0, 0.0);
693 float radius = -1;
694 float tanLambda = 3;
695
696 Helix helix = helixFromCenter(center, radius, tanLambda);
697
698 TMatrixD noMoveJacobian = helix.calcPassiveMoveByJacobian(ROOT::Math::XYZVector(0.0, 0.0, 0.0));
699
700 EXPECT_NEAR(1.0, noMoveJacobian(0, 0), 1e-7);
701 EXPECT_NEAR(0.0, noMoveJacobian(0, 1), 1e-7);
702 EXPECT_NEAR(0.0, noMoveJacobian(0, 2), 1e-7);
703
704 EXPECT_NEAR(0.0, noMoveJacobian(1, 0), 1e-7);
705 EXPECT_NEAR(1.0, noMoveJacobian(1, 1), 1e-7);
706 EXPECT_NEAR(0.0, noMoveJacobian(1, 2), 1e-7);
707
708 EXPECT_NEAR(0.0, noMoveJacobian(2, 0), 1e-7);
709 EXPECT_NEAR(0.0, noMoveJacobian(2, 1), 1e-7);
710 EXPECT_NEAR(1.0, noMoveJacobian(2, 2), 1e-7);
711 }
712
713 TEST_F(HelixTest, calcPassiveMoveByJacobian_orthogonalToPhi0)
714 {
715 ROOT::Math::XYZVector center(0.0, 1.0, 0.0);
716 float radius = -1;
717 float tanLambda = 3;
718
719 Helix helix = helixFromCenter(center, radius, tanLambda);
720
721 TMatrixD moveByOneJacobian = helix.calcPassiveMoveByJacobian(ROOT::Math::XYZVector(0.0, -1.0, 0.0));
722
723 EXPECT_NEAR(1.0, moveByOneJacobian(iD0, iD0), 1e-7);
724 EXPECT_NEAR(0.0, moveByOneJacobian(iD0, iPhi0), 1e-7);
725 EXPECT_NEAR(0.0, moveByOneJacobian(iD0, iOmega), 1e-7);
726
727 EXPECT_NEAR(0.0, moveByOneJacobian(iPhi0, iD0), 1e-7);
728 EXPECT_NEAR(1.0 / 2.0, moveByOneJacobian(iPhi0, iPhi0), 1e-7);
729 EXPECT_NEAR(0.0, moveByOneJacobian(iPhi0, iOmega), 1e-7);
730
731 EXPECT_NEAR(0.0, moveByOneJacobian(iOmega, iD0), 1e-7);
732 EXPECT_NEAR(0.0, moveByOneJacobian(iOmega, iPhi0), 1e-7);
733 EXPECT_NEAR(1.0, moveByOneJacobian(iOmega, iOmega), 1e-7);
734 }
735
736
737 TEST_F(HelixTest, calcPassiveMoveByJacobian_parallelToPhi0_compare_opposite_transformation)
738 {
739 ROOT::Math::XYZVector center(0.0, 1.0, 0.0);
740 float radius = -1;
741 float tanLambda = 3;
742
743 Helix helix = helixFromCenter(center, radius, tanLambda);
744
745 TMatrixD moveByPlusTwoXJacobian = helix.calcPassiveMoveByJacobian(ROOT::Math::XYZVector(2.0, 0.0, 0.0));
746 TMatrixD moveByMinusTwoXJacobian = helix.calcPassiveMoveByJacobian(ROOT::Math::XYZVector(-2.0, 0.0, 0.0));
747
748 EXPECT_NEAR(1.0, moveByPlusTwoXJacobian(iOmega, iOmega), 1e-7);
749 EXPECT_NEAR(0.0, moveByPlusTwoXJacobian(iOmega, iPhi0), 1e-7);
750 EXPECT_NEAR(0.0, moveByPlusTwoXJacobian(iOmega, iD0), 1e-7);
751
752 EXPECT_NEAR(1.0, moveByMinusTwoXJacobian(iOmega, iOmega), 1e-7);
753 EXPECT_NEAR(0.0, moveByMinusTwoXJacobian(iOmega, iPhi0), 1e-7);
754 EXPECT_NEAR(0.0, moveByMinusTwoXJacobian(iOmega, iD0), 1e-7);
755
756 // Symmetric effects
757 EXPECT_NEAR(moveByMinusTwoXJacobian(iD0, iOmega), moveByPlusTwoXJacobian(iD0, iOmega), 1e-7);
758 EXPECT_NEAR(moveByMinusTwoXJacobian(iPhi0, iPhi0), moveByPlusTwoXJacobian(iPhi0, iPhi0), 1e-7);
759 EXPECT_NEAR(moveByMinusTwoXJacobian(iD0, iD0), moveByPlusTwoXJacobian(iD0, iD0), 1e-7);
760
761 // Asymmetric effects
762 EXPECT_NEAR(moveByMinusTwoXJacobian(iD0, iPhi0), -moveByPlusTwoXJacobian(iD0, iPhi0), 1e-7);
763 EXPECT_NEAR(moveByMinusTwoXJacobian(iPhi0, iD0), -moveByPlusTwoXJacobian(iPhi0, iD0), 1e-7);
764 EXPECT_NEAR(moveByMinusTwoXJacobian(iPhi0, iOmega), -moveByPlusTwoXJacobian(iPhi0, iOmega), 1e-7);
765
766 // Signs
767 EXPECT_NEGATIVE(moveByPlusTwoXJacobian(iD0, iPhi0));
768 EXPECT_POSITIVE(moveByPlusTwoXJacobian(iD0, iOmega));
769 EXPECT_POSITIVE(moveByPlusTwoXJacobian(iPhi0, iD0));
770 EXPECT_NEGATIVE(moveByPlusTwoXJacobian(iPhi0, iOmega));
771
772 EXPECT_NEGATIVE(moveByPlusTwoXJacobian(iD0, iD0) - 1);
773 EXPECT_NEGATIVE(moveByPlusTwoXJacobian(iPhi0, iPhi0) - 1);
774 }
775
776 TEST_F(HelixTest, calcPassiveMoveByJacobian_roundtrip)
777 {
778 for (const auto& by : bys) {
779 ROOT::Math::XYZVector center(0.0, 1.0, 0.0);
780 float radius = -1;
781 float tanLambda = 3;
782
783 Helix helix = helixFromCenter(center, radius, tanLambda);
784
785 TMatrixD jacobian = helix.calcPassiveMoveByJacobian(by);
786 helix.passiveMoveBy(by);
787
788 TMatrixD reverseJacobian = helix.calcPassiveMoveByJacobian(-by);
789
790 TMatrixD unitMatrix = reverseJacobian * jacobian;
791
792 for (int i = 0; i < 5; ++i) {
793 for (int j = 0; j < 5; ++j) {
794 if (i == j) {
795 // Diagonal is one.
796 EXPECT_NEAR(1, unitMatrix(i, j), 1e-7);
797 } else {
798 // Off diagonal is zero.
799 EXPECT_NEAR(0, unitMatrix(i, j), 1e-7);
800 }
801 }
802 }
803 }
804
805 }
806
807
808 TEST_F(HelixTest, calcPassiveMoveByJacobian_parallelToPhi0)
809 {
810 ROOT::Math::XYZVector center(1.0, 0.0, 0.0);
811 float radius = -1;
812 float tanLambda = 3;
813
814 Helix helix = helixFromCenter(center, radius, tanLambda);
815
816 TMatrixD moveByTwoYJacobian = helix.calcPassiveMoveByJacobian(ROOT::Math::XYZVector(0.0, 2.0, 0.0));
817
818 // Hand caluclated intermediate quantities;
819 double deltaParallel = 2;
820 double A = 4;
821 double u = 1;
822
823 double nu = 1;
824 //double xi = 1.0 / 5.0;
825 double lambda = 1.0 / (5.0 + 3.0 * sqrt(5.0));
826 double mu = sqrt(5.0) / 10.0;
827 double zeta = 4;
828
829 EXPECT_NEAR(1.0, moveByTwoYJacobian(iOmega, iOmega), 1e-7);
830 EXPECT_NEAR(0.0, moveByTwoYJacobian(iOmega, iPhi0), 1e-7);
831 EXPECT_NEAR(0.0, moveByTwoYJacobian(iOmega, iD0), 1e-7);
832
833 EXPECT_NEAR(2.0 / 5.0, moveByTwoYJacobian(iPhi0, iOmega), 1e-7);
834 EXPECT_NEAR(1.0 / 5.0, moveByTwoYJacobian(iPhi0, iPhi0), 1e-7);
835 EXPECT_NEAR(-2.0 / 5.0, moveByTwoYJacobian(iPhi0, iD0), 1e-7);
836
837 EXPECT_NEAR(mu * zeta - A * lambda, moveByTwoYJacobian(iD0, iOmega), 1e-7);
838 EXPECT_NEAR(2.0 * mu * u * deltaParallel, moveByTwoYJacobian(iD0, iPhi0), 1e-7);
839 EXPECT_NEAR(2.0 * mu * nu, moveByTwoYJacobian(iD0, iD0), 1e-7);
840
841 }
842
843
844 TEST_F(HelixTest, calcPassiveMoveByJacobian_consistency_of_expansion)
845 {
846 ROOT::Math::XYZVector center(0.0, 1.0, 0.0);
847 float radius = -1;
848 float tanLambda = 3;
849
850 Helix helix = helixFromCenter(center, radius, tanLambda);
851
852 for (const auto& by : bys) {
853 TMatrixD jacobian = helix.calcPassiveMoveByJacobian(by);
854 TMatrixD jacobianWithExpansion = helix.calcPassiveMoveByJacobian(by, 10000);
855
856 EXPECT_NEAR(1.0, jacobian(iOmega, iOmega), 1e-7);
857 EXPECT_NEAR(0.0, jacobian(iOmega, iPhi0), 1e-7);
858 EXPECT_NEAR(0.0, jacobian(iOmega, iD0), 1e-7);
859
860 EXPECT_NEAR(1.0, jacobianWithExpansion(iOmega, iOmega), 1e-7);
861 EXPECT_NEAR(0.0, jacobianWithExpansion(iOmega, iPhi0), 1e-7);
862 EXPECT_NEAR(0.0, jacobianWithExpansion(iOmega, iD0), 1e-7);
863
864 EXPECT_NEAR(jacobianWithExpansion(iD0, iD0), jacobian(iD0, iD0), 1e-7);
865 EXPECT_NEAR(jacobianWithExpansion(iD0, iPhi0), jacobian(iD0, iPhi0), 1e-7);
866 EXPECT_NEAR(jacobianWithExpansion(iD0, iOmega), jacobian(iD0, iOmega), 1e-7);
867
868 EXPECT_NEAR(jacobianWithExpansion(iPhi0, iD0), jacobian(iPhi0, iD0), 1e-7);
869 EXPECT_NEAR(jacobianWithExpansion(iPhi0, iPhi0), jacobian(iPhi0, iPhi0), 1e-7);
870 EXPECT_NEAR(jacobianWithExpansion(iPhi0, iOmega), jacobian(iPhi0, iOmega), 1e-7);
871 }
872 }
873
874 TEST_F(HelixTest, realExample)
875 {
876 // Example from Belle I data.
877 std::vector<double> helixParams(5);
878 helixParams[0] = 3.82384;
879 helixParams[1] = std::remainder(3.575595, 2 * M_PI);
880 helixParams[2] = 0.00530726;
881 helixParams[3] = -0.000317335;
882 helixParams[4] = 1.34536;
883
884 ROOT::Math::XYZVector momentum(-0.768755, -0.356297, 1.13994);
885 ROOT::Math::XYZVector position(-1.60794, 3.46933, -0.000317335);
886
887 // Note: The helix parameters already have small mismatches that can be fixed as follows
888 // helixParams[1] = std::atan2(static_cast<double>(momentum.Y()), static_cast<double>(momentum.X()));
889 // helixParams[0] = static_cast<double>(position.Rho());
890
891 const double bZ = 1.5;
892
893 // Test if the cartesian coordinates are at the perigee
894 EXPECT_NEAR(0.0, momentum.X() * position.X() + momentum.Y() * position.Y(), 1e-6);
895
896 Helix helix(helixParams[0], helixParams[1], helixParams[2], helixParams[3], helixParams[4]);
897
898 EXPECT_ALL_NEAR(momentum, helix.getMomentum(bZ), 1e-5);
899 EXPECT_ALL_NEAR(position, helix.getPerigee(), 1e-5);
900 }
901
902
903 TEST_F(HelixTest, omegaForUnitMomentum)
904 {
905 ROOT::Math::XYZVector expectedMomentum(1.0, 0.0, 0.0);
906 ROOT::Math::XYZVector expectedPosition(0.0, 1.0, 0.0);
907 const double expectedCharge = 1;
908 const double bZ = 1.5;
909
910 Helix helix(expectedPosition, expectedMomentum, expectedCharge, bZ);
911
912 double expectedOmega = 0.0044968868700000003; // Omega for one GeV
913 EXPECT_ALL_NEAR(expectedOmega, helix.getOmega(), 1e-7);
914 }
915}
Belle2::CDC::Helix
Helix parameter class.
Definition: Helix.h:48
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::EvtPDLUtil::charge
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44
Belle2::TestHelpers
Some utilities to help with writing unit tests.
Definition: Helix.cc:34
Belle2::TestHelpers::angleNear
bool angleNear(double expected, double actual, double tolerance)
Predicate checking that two angular values are close to each other modulus a 2 * PI difference.
Definition: TestHelpers.cc:71
Belle2::TestHelpers::allNear< Belle2::Helix >
bool allNear< Belle2::Helix >(const Belle2::Helix &expected, const Belle2::Helix &actual, double tolerance)
Predicate checking that all five components of the Helix are close by a maximal error of absError.
Definition: Helix.cc:41
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
std
STL namespace.