Belle II Software  release-08-01-10
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 
19 using namespace std;
20 using namespace Belle2;
21 using namespace HelixParameterIndex;
22 
24 std::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 
33 namespace 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 
54 namespace {
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::TEST_F
TEST_F(GlobalLabelTest, LargeNumberOfTimeDependentParameters)
Test large number of time-dep params for registration and retrieval.
Definition: globalLabel.cc:72
Belle2::operator<<
std::ostream & operator<<(std::ostream &output, const Helix &helix)
Output operator for debugging and the generation of unittest error messages.
Definition: Helix.cc:631
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::tan
double tan(double a)
tan for double
Definition: beamHelpers.h:31
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
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17