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