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