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