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 <boost/range/irange.hpp>
12#include <framework/utilities/TestHelpers.h>
13#include <gtest/gtest.h>
14
15using boost::irange;
16
17using namespace std;
18using namespace Belle2;
19
20namespace {
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::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.