Belle II Software  release-08-01-10
UncertainHelix.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 
9 #include <framework/dataobjects/UncertainHelix.h>
10 #include <TVectorD.h>
11 #include <boost/range/irange.hpp>
12 #include <framework/utilities/TestHelpers.h>
13 #include <gtest/gtest.h>
14 
15 using boost::irange;
16 
17 using namespace std;
18 using namespace Belle2;
19 
20 namespace {
21 
22  TEST(UncertainHelixTest, CartesianCovarianceRoundtripWithPerigeeOnOrigin)
23  {
24  ROOT::Math::XYZVector expectedMomentum(1.0, 0.0, 0.0);
25  ROOT::Math::XYZVector expectedPosition(0.0, 0.0, 0.0);
26  const double expectedCharge = -1;
27  const double bZ = 2;
28  const double pValue = 0.5;
29 
30  TMatrixDSym expectedCov6(6);
31  expectedCov6.Zero();
32  expectedCov6(0, 0) = 0; // There cannot be a covariance in the x direction since the direction along the track has no constrained.
33  expectedCov6(1, 1) = 2;
34  expectedCov6(2, 2) = 3;
35  expectedCov6(3, 3) = 4;
36  expectedCov6(4, 4) = 5;
37  expectedCov6(5, 5) = 6;
38 
39  UncertainHelix uncertainHelix(expectedPosition,
40  expectedMomentum,
41  expectedCharge,
42  bZ,
43  expectedCov6,
44  pValue);
45 
46  const ROOT::Math::XYZVector position = uncertainHelix.getPerigee();
47  const ROOT::Math::XYZVector momentum = uncertainHelix.getMomentum(bZ);
48  const double charge = uncertainHelix.getChargeSign();
49  const TMatrixDSym cov6 = uncertainHelix.getCartesianCovariance(bZ);
50 
51  EXPECT_NEAR(expectedPosition.X(), position.X(), 1e-7);
52  EXPECT_NEAR(expectedPosition.Y(), position.Y(), 1e-7);
53  EXPECT_NEAR(expectedPosition.Z(), position.Z(), 1e-7);
54 
55  EXPECT_NEAR(expectedMomentum.X(), momentum.X(), 1e-7);
56  EXPECT_NEAR(expectedMomentum.Y(), momentum.Y(), 1e-7);
57  EXPECT_NEAR(expectedMomentum.Z(), momentum.Z(), 1e-7);
58 
59  EXPECT_EQ(expectedCharge, charge);
60 
61  for (int i : irange(0, 6)) {
62  for (int j : irange(0, 6)) {
63  EXPECT_NEAR(expectedCov6(i, j), cov6(i, j), 1e-7);
64  }
65  }
66  }
67 
68 
69  TEST(UncertainHelixTest, CartesianCovarianceRoundtripAtPerigeeBesidesOrigin)
70  {
71  ROOT::Math::XYZVector expectedMomentum(1.0, 0.0, 0.0);
72  ROOT::Math::XYZVector expectedPosition(0.0, 1.0, 0.0);
73  const double expectedCharge = -1;
74  const double bZ = 2;
75  const double pValue = 0.5;
76 
77  TMatrixDSym expectedCov6(6);
78  expectedCov6.Zero();
79  expectedCov6(0, 0) = 0; // There cannot be a covariance in the x direction since the direction along the track has no constraint.
80  for (int i : irange(1, 6)) {
81  for (int j : irange(1, 6)) {
82  // Assign some values to the covariances to check the correct propagation of signs.
83  expectedCov6(i, j) = 1;
84  }
85  }
86 
87  expectedCov6(1, 1) = 2;
88  expectedCov6(2, 2) = 3;
89  expectedCov6(3, 3) = 4;
90  expectedCov6(4, 4) = 5;
91  expectedCov6(5, 5) = 6;
92 
93  UncertainHelix uncertainHelix(expectedPosition,
94  expectedMomentum,
95  expectedCharge,
96  bZ,
97  expectedCov6,
98  pValue);
99 
100  const ROOT::Math::XYZVector position = uncertainHelix.getPerigee();
101  const ROOT::Math::XYZVector momentum = uncertainHelix.getMomentum(bZ);
102  const double charge = uncertainHelix.getChargeSign();
103  const TMatrixDSym cov6 = uncertainHelix.getCartesianCovariance(bZ);
104 
105  EXPECT_NEAR(expectedPosition.X(), position.X(), 1e-7);
106  EXPECT_NEAR(expectedPosition.Y(), position.Y(), 1e-7);
107  EXPECT_NEAR(expectedPosition.Z(), position.Z(), 1e-7);
108 
109  EXPECT_NEAR(expectedMomentum.X(), momentum.X(), 1e-7);
110  EXPECT_NEAR(expectedMomentum.Y(), momentum.Y(), 1e-7);
111  EXPECT_NEAR(expectedMomentum.Z(), momentum.Z(), 1e-7);
112 
113  EXPECT_EQ(expectedCharge, charge);
114 
115  // Test the covariance matrix for the right null space
116  TVectorD nullSpace(6);
117  const double factor =
118  (expectedPosition.X() * expectedMomentum.Y() - expectedPosition.Y() * expectedMomentum.X()) / expectedMomentum.Perp2();
119  nullSpace[0] = expectedMomentum.X();
120  nullSpace[1] = expectedMomentum.Y();
121  nullSpace[2] = 0;
122  nullSpace[3] = factor * expectedMomentum.Y();
123  nullSpace[4] = - factor * expectedMomentum.X();
124  nullSpace[5] = 0;
125 
126  TVectorD nullVector = cov6 * nullSpace;
127  for (int j : irange(0, 6)) {
128  EXPECT_NEAR(0, nullVector(j), 1e-7);
129  }
130 
131  // Make a second round trip to assert the covariance
132  // is stable after one adjustment of the null space
133  UncertainHelix uncertainHelix2(position,
134  momentum,
135  charge,
136  bZ,
137  cov6,
138  pValue);
139 
140  const TMatrixDSym cov6_2 = uncertainHelix2.getCartesianCovariance(bZ);
141  for (int i : irange(0, 6)) {
142  for (int j : irange(0, 6)) {
143  EXPECT_NEAR(cov6(i, j), cov6_2(i, j), 1e-7);
144  }
145  }
146  }
147 
148  TEST(UncertainHelixTest, CartesianCovarianceRoundtripAtPerigeeBesidesOriginIndividualCoordinates)
149  {
150  ROOT::Math::XYZVector expectedMomentum(1.0, 0.0, 0.0);
151  ROOT::Math::XYZVector expectedPosition(0.0, 1.0, 0.0);
152  const double expectedCharge = -1;
153  const double bZ = 2;
154  const double pValue = 0.5;
155 
156  for (int i : irange(1, 6)) {
157 
158  TMatrixDSym expectedCov6(6);
159  expectedCov6.Zero();
160  expectedCov6(0, 0) = 0; // There cannot be a covariance in the x direction since the direction along the track has no constraint.
161 
162  expectedCov6(i, i) = 1;
163 
164  UncertainHelix uncertainHelix(expectedPosition,
165  expectedMomentum,
166  expectedCharge,
167  bZ,
168  expectedCov6,
169  pValue);
170 
171  const ROOT::Math::XYZVector position = uncertainHelix.getPerigee();
172  const ROOT::Math::XYZVector momentum = uncertainHelix.getMomentum(bZ);
173  const double charge = uncertainHelix.getChargeSign();
174  const TMatrixDSym cov6 = uncertainHelix.getCartesianCovariance(bZ);
175 
176  EXPECT_NEAR(expectedPosition.X(), position.X(), 1e-7);
177  EXPECT_NEAR(expectedPosition.Y(), position.Y(), 1e-7);
178  EXPECT_NEAR(expectedPosition.Z(), position.Z(), 1e-7);
179 
180  EXPECT_NEAR(expectedMomentum.X(), momentum.X(), 1e-7);
181  EXPECT_NEAR(expectedMomentum.Y(), momentum.Y(), 1e-7);
182  EXPECT_NEAR(expectedMomentum.Z(), momentum.Z(), 1e-7);
183 
184  EXPECT_EQ(expectedCharge, charge);
185 
186  // Test the covariance matrix for the right null space
187  TVectorD nullSpace(6);
188  const double factor =
189  (expectedPosition.X() * expectedMomentum.Y() - expectedPosition.Y() * expectedMomentum.X()) / expectedMomentum.Perp2();
190  nullSpace[0] = expectedMomentum.X();
191  nullSpace[1] = expectedMomentum.Y();
192  nullSpace[2] = 0;
193  nullSpace[3] = factor * expectedMomentum.Y();
194  nullSpace[4] = - factor * expectedMomentum.X();
195  nullSpace[5] = 0;
196 
197  TVectorD nullVector = cov6 * nullSpace;
198  for (int j : irange(0, 6)) {
199  EXPECT_NEAR(0, nullVector(j), 1e-7);
200  }
201 
202  // Make a second round trip to assert the covariance
203  // is stable after one adjustment of the null space
204  UncertainHelix uncertainHelix2(position,
205  momentum,
206  charge,
207  bZ,
208  cov6,
209  pValue);
210 
211  const TMatrixDSym cov6_2 = uncertainHelix2.getCartesianCovariance(bZ);
212  for (int j : irange(0, 6)) {
213  for (int k : irange(0, 6)) {
214  EXPECT_NEAR(cov6(j, k), cov6_2(j, k), 1e-7);
215  }
216  }
217  }
218  }
219 
220  TEST(UncertainHelixTest, PerigeeCovarianceRoundtripAtPerigeeBesidesOriginIndividualCoordinates)
221  {
222  double d0 = 1;
223  double phi0 = 0;
224  double omega = -0.005;
225  double tanLambda = 1;
226  double z0 = 0;
227 
228  const double bZ = 2;
229  const double pValue = 0.5;
230 
231  for (int i : irange(0, 5)) {
232 
233  TMatrixDSym expectedCov5(5);
234  expectedCov5.Zero();
235  expectedCov5(i, i) = 1;
236 
237  UncertainHelix uncertainHelix(d0, phi0, omega, z0, tanLambda, expectedCov5, pValue);
238 
239  const ROOT::Math::XYZVector position = uncertainHelix.getPerigee();
240  const ROOT::Math::XYZVector momentum = uncertainHelix.getMomentum(bZ);
241  const double charge = uncertainHelix.getChargeSign();
242  const TMatrixDSym cov6 = uncertainHelix.getCartesianCovariance(bZ);
243 
244  UncertainHelix uncertainHelix2(position,
245  momentum,
246  charge,
247  bZ,
248  cov6,
249  pValue);
250 
251  EXPECT_NEAR(d0, uncertainHelix2.getD0(), 1e-7);
252  EXPECT_NEAR(phi0, uncertainHelix2.getPhi0(), 1e-7);
253  EXPECT_NEAR(omega, uncertainHelix2.getOmega(), 1e-7);
254 
255  EXPECT_NEAR(z0, uncertainHelix2.getZ0(), 1e-7);
256  EXPECT_NEAR(tanLambda, uncertainHelix2.getTanLambda(), 1e-7);
257 
258  const TMatrixDSym& cov5 = uncertainHelix2.getCovariance();
259  EXPECT_EQ(-1, charge);
260 
261  for (int j : irange(0, 5)) {
262  for (int k : irange(0, 5)) {
263  EXPECT_NEAR(expectedCov5(j, k), cov5(j, k), 1e-7);
264  }
265  }
266  }
267  }
268 
269  TEST(UncertainHelixTest, RealPerigeeCovarianceRoundTrip)
270  {
271  const double bZ = 1.5;
272  const double pValue = 0.5;
273 
274  std::vector<float> helixParams(5);
275  std::vector<float> helixCovariance(15);
276 
277  helixParams[0] = 3.82384;
278  helixParams[1] = std::remainder(3.57559, 2 * M_PI);
279  helixParams[2] = 0.00530726;
280  helixParams[3] = -0.000317335;
281  helixParams[4] = 1.34536;
282 
283  helixCovariance[0] = 0.000648818;
284  helixCovariance[1] = 6.357e-05;
285  helixCovariance[2] = 2.04344e-07;
286  helixCovariance[3] = -0.00101879;
287  helixCovariance[4] = 2.29407e-05;
288 
289  helixCovariance[5] = 1.02934e-05;
290  helixCovariance[6] = 1.39609e-08;
291  helixCovariance[7] = 9.57685e-07;
292  helixCovariance[8] = 1.19863e-06;
293 
294  helixCovariance[9] = 4.96093e-10;
295  helixCovariance[10] = -8.41612e-07;
296  helixCovariance[11] = 6.50684e-08;
297 
298  helixCovariance[12] = 0.0317385;
299  helixCovariance[13] = -0.00065476;
300 
301  helixCovariance[14] = 3.6521e-05;
302 
303  TMatrixDSym expectedCov5(5);
304  unsigned int counter = 0;
305  for (unsigned int i = 0; i < 5; ++i) {
306  for (unsigned int j = i; j < 5; ++j) {
307  expectedCov5(i, j) = expectedCov5(j, i) = helixCovariance[counter];
308  ++counter;
309  }
310  }
311 
312  UncertainHelix uncertainHelix(helixParams[0],
313  helixParams[1],
314  helixParams[2],
315  helixParams[3],
316  helixParams[4],
317  expectedCov5,
318  pValue);
319 
320  const ROOT::Math::XYZVector position = uncertainHelix.getPerigee();
321  const ROOT::Math::XYZVector momentum = uncertainHelix.getMomentum(bZ);
322  const double charge = uncertainHelix.getChargeSign();
323  const TMatrixDSym cov6 = uncertainHelix.getCartesianCovariance(bZ);
324 
325  UncertainHelix uncertainHelix2(position,
326  momentum,
327  charge,
328  bZ,
329  cov6,
330  pValue);
331 
332  EXPECT_NEAR(helixParams[0], uncertainHelix2.getD0(), 1e-7);
333  EXPECT_NEAR(helixParams[1], uncertainHelix2.getPhi0(), 1e-7);
334  EXPECT_NEAR(helixParams[2], uncertainHelix2.getOmega(), 1e-7);
335 
336  EXPECT_NEAR(helixParams[3], uncertainHelix2.getZ0(), 1e-7);
337  EXPECT_NEAR(helixParams[4], uncertainHelix2.getTanLambda(), 1e-7);
338 
339  const TMatrixDSym& cov5 = uncertainHelix2.getCovariance();
340  EXPECT_EQ(1, charge);
341 
342  for (int i : boost::irange(0, 5)) {
343  for (int j : boost::irange(0, 5)) {
344  EXPECT_NEAR(expectedCov5(i, j), cov5(i, j), 1e-7);
345  }
346  }
347  }
348 
349 
350  TEST(UncertainHelixTest, Explicit_D0)
351  {
352  double d0 = 1;
353  double phi0 = 0;
354  double omega = -0.005;
355  double tanLambda = 1;
356  double z0 = 0;
357 
358  const double bZ = 2;
359  const double pValue = 0.5;
360 
361  TMatrixDSym expectedCov5(5);
362  expectedCov5.Zero();
363  expectedCov5(0, 0) = 1;
364 
365  TMatrixDSym expectedCov6(6);
366  expectedCov6.Zero();
367  expectedCov6(1, 1) = 1;
368 
369  UncertainHelix uncertainHelix(d0, phi0, omega, z0, tanLambda, expectedCov5, pValue);
370 
371  const TMatrixDSym cov6 = uncertainHelix.getCartesianCovariance(bZ);
372 
373  for (int i : irange(0, 6)) {
374  for (int j : irange(0, 6)) {
375  TEST_CONTEXT("Failed for index (" << i << ", " << j << ")");
376  EXPECT_NEAR(expectedCov6(i, j), cov6(i, j), 1e-7);
377  }
378  }
379  }
380 
381  TEST(UncertainHelixTest, Explicit_Phi0)
382  {
383  double d0 = 0;
384  double phi0 = 0;
385  double omega = -0.006;
386  double tanLambda = 0;
387  double z0 = 0;
388 
389  const double bZ = 1.5;
390  const double pValue = 0.5;
391 
392  TMatrixDSym expectedCov5(5);
393  expectedCov5.Zero();
394  expectedCov5(1, 1) = 1;
395 
396  UncertainHelix uncertainHelix(d0, phi0, omega, z0, tanLambda, expectedCov5, pValue);
397 
398  const ROOT::Math::XYZVector position = uncertainHelix.getPerigee();
399  const ROOT::Math::XYZVector momentum = uncertainHelix.getMomentum(bZ);
400  const double charge = uncertainHelix.getChargeSign();
401  const TMatrixDSym cov6 = uncertainHelix.getCartesianCovariance(bZ);
402 
403  double absMom2D = uncertainHelix.getTransverseMomentum(1.5);
404 
405  TMatrixDSym expectedCov6(6);
406  expectedCov6.Zero();
407  expectedCov6(4, 4) = absMom2D * absMom2D;
408 
409  for (int i : irange(0, 6)) {
410  for (int j : irange(0, 6)) {
411  TEST_CONTEXT("Failed for index (" << i << ", " << j << ")");
412  EXPECT_NEAR(expectedCov6(i, j), cov6(i, j), 1e-7);
413  }
414  }
415 
416  UncertainHelix uncertainHelix2(position,
417  momentum,
418  charge,
419  bZ,
420  cov6,
421  pValue);
422 
423  EXPECT_NEAR(d0, uncertainHelix2.getD0(), 1e-7);
424  EXPECT_NEAR(phi0, uncertainHelix2.getPhi0(), 1e-7);
425  EXPECT_NEAR(omega, uncertainHelix2.getOmega(), 1e-7);
426 
427  EXPECT_NEAR(z0, uncertainHelix2.getZ0(), 1e-7);
428  EXPECT_NEAR(tanLambda, uncertainHelix2.getTanLambda(), 1e-7);
429 
430  const TMatrixDSym& cov5 = uncertainHelix2.getCovariance();
431  EXPECT_EQ(-1, charge);
432 
433  for (int i : irange(0, 5)) {
434  for (int j : irange(0, 5)) {
435  TEST_CONTEXT("Failed for index (" << i << ", " << j << ")");
436  EXPECT_NEAR(expectedCov5(i, j), cov5(i, j), 1e-7);
437  }
438  }
439  }
440 
441  TEST(UncertainHelixTest, MoveInCartesianVersusMoveInPerigeeCoordinates)
442  {
443  // Move of the coordinate system
444  ROOT::Math::XYZVector moveBy(0.5, -1, 0.5);
445 
446  double d0 = 1;
447  double phi0 = 0;
448  double omega = -0.005;
449  double tanLambda = 1;
450  double z0 = 0;
451 
452  const double bZ = 2;
453  const double pValue = 0.5;
454 
455  TMatrixDSym initialCov5(5);
456  initialCov5.Zero();
457  for (int i : irange(0, 5)) {
458  for (int j : irange(0, 5)) {
459  initialCov5(i, j) = i + j;
460  }
461  }
462 
463  UncertainHelix uncertainHelix(d0, phi0, omega, z0, tanLambda, initialCov5, pValue);
464 
465  const ROOT::Math::XYZVector position = uncertainHelix.getPerigee();
466  const ROOT::Math::XYZVector momentum = uncertainHelix.getMomentum(bZ);
467  const double charge = uncertainHelix.getChargeSign();
468  const TMatrixDSym cov6 = uncertainHelix.getCartesianCovariance(bZ);
469 
470  // Execute the move in cartesian coordinates
471  const ROOT::Math::XYZVector movedPosition = position - moveBy;
472  UncertainHelix movedUncertainHelix(movedPosition,
473  momentum,
474  charge,
475  bZ,
476  cov6,
477  pValue);
478 
479  const TMatrixDSym& cov5 = movedUncertainHelix.getCovariance();
480 
481  // Execute the move in perigee coordinates
482  UncertainHelix expectedMovedUncertainHelix = uncertainHelix;
483  expectedMovedUncertainHelix.passiveMoveBy(moveBy);
484  const TMatrixDSym expectedCov5 = expectedMovedUncertainHelix.getCovariance();
485 
486  // Test equivalence of both transformation
487  EXPECT_NEAR(expectedMovedUncertainHelix.getD0(), movedUncertainHelix.getD0(), 1e-7);
488  EXPECT_NEAR(expectedMovedUncertainHelix.getPhi0(), movedUncertainHelix.getPhi0(), 1e-7);
489  EXPECT_NEAR(expectedMovedUncertainHelix.getOmega(), movedUncertainHelix.getOmega(), 1e-7);
490  EXPECT_NEAR(expectedMovedUncertainHelix.getZ0(), movedUncertainHelix.getZ0(), 1e-7);
491  EXPECT_NEAR(expectedMovedUncertainHelix.getTanLambda(), movedUncertainHelix.getTanLambda(), 1e-7);
492 
493  for (int i : irange(0, 5)) {
494  for (int j : irange(0, 5)) {
495  TEST_CONTEXT("Failed for index (" << i << ", " << j << ")");
496  EXPECT_NEAR(expectedCov5(i, j), cov5(i, j), 1e-7);
497  }
498  }
499  }
500 }
Belle2::UncertainHelix
This class represents an ideal helix in perigee parameterization including the covariance matrix of t...
Definition: UncertainHelix.h:36
Belle2::UncertainHelix::passiveMoveBy
double passiveMoveBy(const ROOT::Math::XYZVector &by)
Moves origin of the coordinate system (passive transformation) by the given vector.
Definition: UncertainHelix.h:117
Belle2::UncertainHelix::getCovariance
const TMatrixDSym & getCovariance() const
Getter for covariance matrix of perigee parameters in matrix form.
Definition: UncertainHelix.h:100
Belle2::TEST
TEST(TestgetDetectorRegion, TestgetDetectorRegion)
Test Constructors.
Definition: utilityFunctions.cc:25
Belle2::EvtPDLUtil::charge
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17