EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ParameterSetTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ParameterSetTests.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2016-2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <boost/test/unit_test.hpp>
10 
14 
15 #include <cmath>
16 #include <memory>
17 #include <random>
18 
19 using namespace Acts::detail;
20 
24 namespace Acts {
28 namespace Test {
30 namespace {
31 // tolerance used for floating point comparison in this translation unit
32 const double tol = 1e-6;
33 
34 double get_cyclic_value(double value, double min, double max) {
35  return value - (max - min) * std::floor((value - min) / (max - min));
36 }
37 
38 double get_cyclic_difference(double a, double b, double min, double max) {
39  const double period = max - min;
40  const double half_period = period / 2;
41  a = get_cyclic_value(a, min, max);
42  b = get_cyclic_value(b, min, max);
43  double raw_diff = a - b;
44  double diff =
45  (raw_diff > half_period)
46  ? raw_diff - period
47  : ((raw_diff < -half_period) ? period + raw_diff : raw_diff);
48  return diff;
49 }
50 
51 void check_residuals_for_bound_parameters() {
54  double theta_1 = 0.7 * M_PI;
55  double theta_2 = 0.4 * M_PI;
56  ActsVectorD<1> dTheta;
57  dTheta << (theta_1 - theta_2);
58 
59  // both parameters inside bounds, difference is positive
60  ParameterSet<BoundIndices, eBoundTheta> bound1(std::nullopt, theta_1);
61  ParameterSet<BoundIndices, eBoundTheta> bound2(std::nullopt, theta_2);
62  CHECK_CLOSE_REL(bound1.residual(bound2), dTheta, tol);
63 
64  // both parameters inside bound, difference negative
65  dTheta << (theta_2 - theta_1);
66  CHECK_CLOSE_REL(bound2.residual(bound1), dTheta, tol);
67 
68  // one parameter above upper bound, difference positive
69  theta_1 = max + 1;
70  bound1.setParameter<eBoundTheta>(theta_1);
71  dTheta << max - theta_2;
72  CHECK_CLOSE_REL(bound1.residual(bound2), dTheta, tol);
73 
74  // one parameter above upper bound, difference negative
75  dTheta << theta_2 - max;
76  CHECK_CLOSE_REL(bound2.residual(bound1), dTheta, tol);
77 
78  // one parameter below lower bound, difference positive
79  theta_1 = min - 1;
80  bound1.setParameter<eBoundTheta>(theta_1);
81  dTheta << theta_2 - min;
82  CHECK_CLOSE_REL(bound2.residual(bound1), dTheta, tol);
83 
84  // one parameter below lower bound, difference negative
85  dTheta << min - theta_2;
86  CHECK_CLOSE_REL(bound1.residual(bound2), dTheta, tol);
87 
88  // both parameters outside bounds, both below
89  theta_1 = min - 1;
90  theta_2 = min - 2;
91  bound1.setParameter<eBoundTheta>(theta_1);
92  bound2.setParameter<eBoundTheta>(theta_2);
93  CHECK_SMALL(bound1.residual(bound2), tol);
94 
95  // both parameters outside bounds, both above
96  theta_1 = max + 1;
97  theta_2 = max + 2;
98  bound1.setParameter<eBoundTheta>(theta_1);
99  bound2.setParameter<eBoundTheta>(theta_2);
100  CHECK_SMALL(bound1.residual(bound2), tol);
101 
102  // both parameters outside bounds, one above, one below
103  theta_1 = max + 1;
104  theta_2 = min - 2;
105  bound1.setParameter<eBoundTheta>(theta_1);
106  bound2.setParameter<eBoundTheta>(theta_2);
107  dTheta << max - min;
108  CHECK_CLOSE_REL(bound1.residual(bound2), dTheta, tol);
109  dTheta << min - max;
110  CHECK_CLOSE_REL(bound2.residual(bound1), dTheta, tol);
111 }
112 
113 void check_residuals_for_cyclic_parameters() {
116 
117  double phi_1 = 0.7 * M_PI;
118  double phi_2 = 0.4 * M_PI;
119  ActsVectorD<1> dPhi;
120  dPhi << (phi_1 - phi_2);
121 
122  ParameterSet<BoundIndices, eBoundPhi> cyclic1(std::nullopt, phi_1);
123  ParameterSet<BoundIndices, eBoundPhi> cyclic2(std::nullopt, phi_2);
124 
125  // no boundary crossing, difference is positive
126  CHECK_CLOSE_REL(cyclic1.residual(cyclic2), dPhi, tol);
127 
128  // no boundary crossing, difference is negative
129  CHECK_CLOSE_REL(cyclic2.residual(cyclic1), -dPhi, tol);
130 
131  // forward boundary crossing
132  phi_1 = -0.9 * M_PI;
133  cyclic1.setParameter<eBoundPhi>(phi_1);
134  dPhi << get_cyclic_difference(phi_1, phi_2, min, max);
135  CHECK_CLOSE_REL(cyclic1.residual(cyclic2), dPhi, tol);
136  CHECK_CLOSE_REL(cyclic2.residual(cyclic1), -dPhi, tol);
137 
138  // backward boundary crossing
139  phi_1 = 0.7 * M_PI;
140  phi_2 = -0.9 * M_PI;
141  cyclic1.setParameter<eBoundPhi>(phi_1);
142  cyclic2.setParameter<eBoundPhi>(phi_2);
143  dPhi << get_cyclic_difference(phi_1, phi_2, min, max);
144  CHECK_CLOSE_REL(cyclic1.residual(cyclic2), dPhi, tol);
145  CHECK_CLOSE_REL(cyclic2.residual(cyclic1), -dPhi, tol);
146 }
147 
148 void random_residual_tests() {
149  // random number generators
150  std::default_random_engine e;
151  std::uniform_real_distribution<float> uniform_dist(-1000, 300);
152 
153  const double theta_max = ParameterTraits<BoundIndices, eBoundTheta>::max;
154  const double theta_min = ParameterTraits<BoundIndices, eBoundTheta>::min;
155  const double phi_max = ParameterTraits<BoundIndices, eBoundPhi>::max;
156  const double phi_min = ParameterTraits<BoundIndices, eBoundPhi>::min;
157 
158  BoundVector parValues_1;
159  BoundVector parValues_2;
160  FullBoundParameterSet parSet_1(std::nullopt, parValues_1);
161  FullBoundParameterSet parSet_2(std::nullopt, parValues_2);
162  BoundVector residual;
163  const unsigned int toys = 1000;
164  for (unsigned int i = 0; i < toys; ++i) {
165  const double loc0_1 = uniform_dist(e);
166  const double loc1_1 = uniform_dist(e);
167  const double phi_1 = uniform_dist(e);
168  const double theta_1 = uniform_dist(e);
169  const double qop_1 = uniform_dist(e);
170  parValues_1 << loc0_1, loc1_1, phi_1, theta_1, qop_1, 0.;
171  parSet_1.setParameters(parValues_1);
172 
173  const double loc0_2 = uniform_dist(e);
174  const double loc1_2 = uniform_dist(e);
175  const double phi_2 = uniform_dist(e);
176  const double theta_2 = uniform_dist(e);
177  const double qop_2 = uniform_dist(e);
178  parValues_2 << loc0_2, loc1_2, phi_2, theta_2, qop_2, 0.;
179  parSet_2.setParameters(parValues_2);
180 
181  const double delta_loc0 = loc0_1 - loc0_2;
182  const double delta_loc1 = loc1_1 - loc1_2;
183  // for theta make sure that the difference calculation considers the
184  // restricted value range
185  const double delta_theta =
186  (theta_1 > theta_max ? theta_max
187  : (theta_1 < theta_min ? theta_min : theta_1)) -
188  (theta_2 > theta_max ? theta_max
189  : (theta_2 < theta_min ? theta_min : theta_2));
190  const double delta_qop = qop_1 - qop_2;
191  residual = parSet_1.residual(parSet_2);
192 
193  // local parameters are unbound -> check for usual difference
194  if (std::abs(residual(0) - delta_loc0) > tol) {
195  BOOST_CHECK(false);
196  break;
197  }
198  if (std::abs(residual(1) - delta_loc1) > tol) {
199  BOOST_CHECK(false);
200  break;
201  }
202  // phi is a cyclic parameter -> check that (unsigned) difference is not
203  // larger than half period
204  // check that corrected(corrected(phi_2) + residual) == corrected(phi_1)
205  if (std::abs(get_cyclic_value(
206  get_cyclic_value(phi_2, phi_min, phi_max) + residual(2),
207  phi_min, phi_max) -
208  get_cyclic_value(phi_1, phi_min, phi_max)) > tol or
209  std::abs(residual(2)) > (phi_max - phi_min) / 2) {
210  BOOST_CHECK(false);
211  break;
212  }
213  // theta is bound -> check that (unsigned) difference is not larger then
214  // allowed range, check corrected difference
215  if (std::abs(residual(3) - delta_theta) > tol or
216  std::abs(residual(3)) > (theta_max - theta_min)) {
217  BOOST_CHECK(false);
218  break;
219  }
220  // qop is unbound -> check usual difference
221  if (std::abs(residual(4) - delta_qop) > tol) {
222  BOOST_CHECK(false);
223  break;
224  }
225  }
226 }
227 
228 void free_random_residual_tests() {
229  // random number generators
230  std::default_random_engine e;
231  std::uniform_real_distribution<float> uniform_dist(-1000, 300);
232 
233  FreeVector parValues_1;
234  FreeVector parValues_2;
235  FullFreeParameterSet parSet_1(std::nullopt, parValues_1);
236  FullFreeParameterSet parSet_2(std::nullopt, parValues_2);
237  FreeVector residual;
238  const unsigned int toys = 1000;
239  for (unsigned int i = 0; i < toys; ++i) {
240  const double x1 = uniform_dist(e);
241  const double y1 = uniform_dist(e);
242  const double z1 = uniform_dist(e);
243  const double t1 = uniform_dist(e);
244  const double tx1 = uniform_dist(e);
245  const double ty1 = uniform_dist(e);
246  const double tz1 = uniform_dist(e);
247  const double qop1 = uniform_dist(e);
248  parValues_1 << x1, y1, z1, t1, tx1, ty1, tz1, qop1;
249  parSet_1.setParameters(parValues_1);
250 
251  const double x2 = uniform_dist(e);
252  const double y2 = uniform_dist(e);
253  const double z2 = uniform_dist(e);
254  const double t2 = uniform_dist(e);
255  const double tx2 = uniform_dist(e);
256  const double ty2 = uniform_dist(e);
257  const double tz2 = uniform_dist(e);
258  const double qop2 = uniform_dist(e);
259  parValues_2 << x2, y2, z2, t2, tx2, ty2, tz2, qop2;
260  parSet_2.setParameters(parValues_2);
261 
262  const double delta_x = x1 - x2;
263  const double delta_y = y1 - y2;
264  const double delta_z = z1 - z2;
265  const double delta_t = t1 - t2;
266  const double delta_tx = tx1 - tx2;
267  const double delta_ty = ty1 - ty2;
268  const double delta_tz = tz1 - tz2;
269  const double delta_qop = qop1 - qop2;
270  residual = parSet_1.residual(parSet_2);
271 
272  // local parameters are unbound -> check for usual difference
273  if (std::abs(residual(0) - delta_x) > tol) {
274  BOOST_CHECK(false);
275  break;
276  }
277  if (std::abs(residual(1) - delta_y) > tol) {
278  BOOST_CHECK(false);
279  break;
280  }
281  if (std::abs(residual(2) - delta_z) > tol) {
282  BOOST_CHECK(false);
283  break;
284  }
285  if (std::abs(residual(3) - delta_t) > tol) {
286  BOOST_CHECK(false);
287  break;
288  }
289  if (std::abs(residual(4) - delta_tx) > tol) {
290  BOOST_CHECK(false);
291  break;
292  }
293  if (std::abs(residual(5) - delta_ty) > tol) {
294  BOOST_CHECK(false);
295  break;
296  }
297  if (std::abs(residual(6) - delta_tz) > tol) {
298  BOOST_CHECK(false);
299  break;
300  }
301  if (std::abs(residual(7) - delta_qop) > tol) {
302  BOOST_CHECK(false);
303  break;
304  }
305  }
306 }
307 } // namespace
309 
323 BOOST_AUTO_TEST_CASE(parset_consistency_tests) {
324  // One-dimensional constructions
326  cov1D << 0.5;
327 
328  ActsVectorD<1> vec1D;
329  vec1D << 0.1;
330 
331  ParameterSet<BoundIndices, eBoundLoc0> parSet_with_cov1D_real(cov1D, 0.1);
332 
333  // This line does not compile
334  ParameterSet<BoundIndices, eBoundLoc0> parSet_with_cov1D_vec1D(cov1D, vec1D);
335 
336  // Two-dimensional constructions
338  cov2D << 0.5, 0., 0.8, 0.;
339 
340  ActsVectorD<2> vec2D;
341  vec2D << 0.1, 0.3;
342 
344  cov2D, 0.1, 0.3);
345 
347  cov2D, vec2D);
348 
349  // check template parameter based information
350  BOOST_CHECK(
352 
353  // covariance matrix
355  cov << 1, 0, 0, 0, 1.2, 0.2, 0, 0.2, 0.7;
356 
357  // parameter values
358  double loc0 = 0.5;
359  double loc1 = -0.2;
360  double phi = 0.3 * M_PI; // this should be within [-M_PI,M_PI) to avoid
361  // failed tests due to angle range corrections
362  Vector3D parValues(loc0, loc1, phi);
363 
364  // parameter set with covariance matrix, and three parameters as pack
366  cov, loc0, loc1, phi);
367 
368  // parameter set with covariance matrix, and a vector
369  ActsVectorD<3> vec;
370  vec << loc0, loc1, phi;
372  parSet_with_cov_vec(cov, vec);
373 
374  // check number and type of stored parameters
375  BOOST_CHECK(parSet_with_cov.size() == 3);
376  BOOST_CHECK(parSet_with_cov.contains<eBoundLoc0>());
377  BOOST_CHECK(parSet_with_cov.contains<eBoundLoc1>());
378  BOOST_CHECK(parSet_with_cov.contains<eBoundPhi>());
379  BOOST_CHECK(not parSet_with_cov.contains<eBoundTheta>());
380  BOOST_CHECK(not parSet_with_cov.contains<eBoundQOverP>());
381  BOOST_CHECK(not parSet_with_cov.contains<eBoundTime>());
382 
383  // check stored parameter values
384  BOOST_CHECK(parSet_with_cov.getParameter<eBoundLoc0>() == loc0);
385  BOOST_CHECK(parSet_with_cov.getParameter<eBoundLoc1>() == loc1);
386  BOOST_CHECK(parSet_with_cov.getParameter<eBoundPhi>() == phi);
387  BOOST_CHECK(parSet_with_cov.getParameters() == parValues);
388 
389  // check stored covariance
390  BOOST_CHECK(parSet_with_cov.getCovariance());
391  BOOST_CHECK(*parSet_with_cov.getCovariance() == cov);
392  BOOST_CHECK(parSet_with_cov.getUncertainty<eBoundLoc0>() == sqrt(cov(0, 0)));
393  BOOST_CHECK(parSet_with_cov.getUncertainty<eBoundLoc1>() == sqrt(cov(1, 1)));
394  BOOST_CHECK(parSet_with_cov.getUncertainty<eBoundPhi>() == sqrt(cov(2, 2)));
395 
396  // same parameter set without covariance matrix
398  parSet_without_cov(std::nullopt, parValues);
399 
400  BOOST_CHECK(!parSet_without_cov.getCovariance());
401  BOOST_CHECK(parSet_without_cov.getUncertainty<eBoundLoc0>() < 0);
402  BOOST_CHECK(parSet_without_cov.getUncertainty<eBoundLoc1>() < 0);
403  BOOST_CHECK(parSet_without_cov.getUncertainty<eBoundPhi>() < 0);
404  BOOST_CHECK(parSet_without_cov.getParameters() ==
405  parSet_with_cov.getParameters());
406 
407  // set new covariance matrix
408  parSet_without_cov.setCovariance(cov);
409 
410  BOOST_CHECK(parSet_without_cov.getCovariance());
411  BOOST_CHECK(*parSet_without_cov.getCovariance() == cov);
412 
413  // set new parameter values
414  double newLoc0 = 0.1;
415  double newLoc1 = 0.6;
416  double newPhi = -0.15 * M_PI;
417  parValues << newLoc0, newLoc1, newPhi;
418  parSet_with_cov.setParameter<eBoundLoc0>(newLoc0);
419  parSet_with_cov.setParameter<eBoundLoc1>(newLoc1);
420  parSet_with_cov.setParameter<eBoundPhi>(newPhi);
421 
422  BOOST_CHECK(parSet_with_cov.getParameter<eBoundLoc0>() == newLoc0);
423  BOOST_CHECK(parSet_with_cov.getParameter<eBoundLoc1>() == newLoc1);
424  BOOST_CHECK(parSet_with_cov.getParameter<eBoundPhi>() == newPhi);
425  BOOST_CHECK(parSet_with_cov.getParameters() == parValues);
426 }
427 
435 BOOST_AUTO_TEST_CASE(parset_copy_assignment_tests) {
436  // covariance matrix
438  cov << 1, 0, 0, 0, 1.2, 0.2, 0, 0.2, 0.7;
439 
440  // parameter values
441  double loc0 = 0.5;
442  double loc1 = -0.2;
443  double phi = 0.3 * M_PI; // this should be within [-M_PI,M_PI) to avoid
444  // failed tests due to angle range corrections
445  Vector3D firstParValues(loc0, loc1, phi);
446 
447  // parameter set with covariance matrix
449  cov, loc0, loc1, phi);
450 
451  // check copy constructor
453  BOOST_CHECK(first == copy);
454 
455  // check move constructor
457  std::move(copy));
458  BOOST_CHECK(first == moved);
459 
460  // check assignment operator
462  moved;
463  BOOST_CHECK(assigned == moved);
464 
466  std::nullopt, 0, 1.7, -0.15);
467  BOOST_CHECK(assigned != other);
468  assigned = other;
469  BOOST_CHECK(assigned == other);
470 
471  // check move assignment
472  BOOST_CHECK(first != assigned);
473  first =
475  BOOST_CHECK(first == assigned);
476 
477  // check swap method
479  loc1, phi);
481  std::nullopt, 2 * loc0, 2 * loc1, 2 * phi);
484 
485  BOOST_CHECK(lhs != rhs && lhs == lhs_copy && rhs == rhs_copy);
486 }
487 
493 BOOST_AUTO_TEST_CASE(parset_comparison_tests) {
494  // covariance matrix
496  cov << 1, 0, 0, 0, 1.2, 0.2, 0, 0.2, 0.7;
497 
498  // parameter values
499  double loc0 = 0.5;
500  double loc1 = -0.2;
501  double phi = 0.3 * M_PI; // this should be within [-M_PI,M_PI) to avoid
502  // failed tests due to angle range corrections
503 
504  // parameter set with covariance matrix
506  cov, loc0, loc1, phi);
508  std::nullopt, 2 * loc0, 2 * loc1, 2 * phi);
509 
510  // check self comparison
511  BOOST_CHECK(first == first);
512  BOOST_CHECK(not(first != first));
513 
514  // check mutual exclusivity
515  BOOST_CHECK(first != second);
516  BOOST_CHECK(not(first == second));
517  first = second;
518  BOOST_CHECK(first == second);
519 
520  // check that comparison fails for unequal parameter values
521  second.setParameter<eBoundLoc0>(3 * loc0);
522  BOOST_CHECK(first != second);
523  first = second;
524  BOOST_CHECK(first == second);
525 
526  second.setParameter<eBoundLoc1>(3 * loc1);
527  BOOST_CHECK(first != second);
528  first = second;
529  BOOST_CHECK(first == second);
530 
531  second.setParameter<eBoundPhi>(3 * phi);
532  BOOST_CHECK(first != second);
533  first = second;
534  BOOST_CHECK(first == second);
535 
536  // check that comparison fails for unequal covariance matrices
537  second.setCovariance(cov);
538  BOOST_CHECK(first != second);
539  first = second;
540  BOOST_CHECK(first == second);
541 
542  cov(0, 0) = 2 * cov(0, 0);
543  second.setCovariance(cov);
544  BOOST_CHECK(first != second);
545  first = second;
546  BOOST_CHECK(first == second);
547 }
548 
558 BOOST_AUTO_TEST_CASE(parset_projection_tests) {
559  // clang-format off
561  phi_proj << 0, 0, 1, 0, 0, 0;
562 
563  ActsMatrixD<2, eBoundSize> loc0_qop_proj;
564  loc0_qop_proj << 1, 0, 0, 0, 0, 0,
565  0, 0, 0, 0, 1, 0;
566 
567  ActsMatrixD<2, eBoundSize> loc1_theta_proj;
568  loc1_theta_proj << 0, 1, 0, 0, 0, 0,
569  0, 0, 0, 1, 0, 0;
570 
571  ActsMatrixD<3, eBoundSize> loc0_loc1_phi_proj;
572  loc0_loc1_phi_proj << 1, 0, 0, 0, 0, 0,
573  0, 1, 0, 0, 0, 0,
574  0, 0, 1, 0, 0, 0;
575 
576  ActsMatrixD<4, eBoundSize> loc0_phi_theta_qop_proj;
577  loc0_phi_theta_qop_proj << 1, 0, 0, 0, 0, 0,
578  0, 0, 1, 0, 0, 0,
579  0, 0, 0, 1, 0, 0,
580  0, 0, 0, 0, 1, 0;
581 
582  ActsMatrixD<5, eBoundSize> loc0_loc1_phi_theta_qop_proj;
583  loc0_loc1_phi_theta_qop_proj << 1, 0, 0, 0, 0, 0,
584  0, 1, 0, 0, 0, 0,
585  0, 0, 1, 0, 0, 0,
586  0, 0, 0, 1, 0, 0,
587  0, 0, 0, 0, 1, 0;
588 
590  loc0_loc1_phi_theta_qop_t_proj;
591  loc0_loc1_phi_theta_qop_t_proj << 1, 0, 0, 0, 0, 0,
592  0, 1, 0, 0, 0, 0,
593  0, 0, 1, 0, 0, 0,
594  0, 0, 0, 1, 0, 0,
595  0, 0, 0, 0, 1, 0,
596  0, 0, 0, 0, 0, 1;
597  // clang-format on
598 
599  BOOST_CHECK((ParameterSet<BoundIndices, eBoundPhi>::projector() == phi_proj));
600  BOOST_CHECK(
602  loc0_qop_proj));
603  BOOST_CHECK(
605  loc1_theta_proj));
607  eBoundPhi>::projector() == loc0_loc1_phi_proj));
608  BOOST_CHECK(
610  eBoundQOverP>::projector() == loc0_phi_theta_qop_proj));
612  eBoundTheta, eBoundQOverP>::projector() ==
613  loc0_loc1_phi_theta_qop_proj));
614  BOOST_CHECK(
616  eBoundTheta, eBoundQOverP, eBoundTime>::projector() ==
617  loc0_loc1_phi_theta_qop_t_proj));
618 }
619 
631 BOOST_AUTO_TEST_CASE(parset_residual_tests) {
632  // check unbound parameter type
633  const double large_number = 12443534120;
634  const double small_number = -924342675;
635  const double normal_number = 1.234;
637  std::nullopt, small_number, large_number, normal_number);
638  BOOST_CHECK(unbound.getParameter<eBoundLoc0>() == small_number);
639  BOOST_CHECK(unbound.getParameter<eBoundLoc1>() == large_number);
640  BOOST_CHECK(unbound.getParameter<eBoundQOverP>() == normal_number);
641 
642  // check bound parameter type
643  ParameterSet<BoundIndices, eBoundTheta> bound(std::nullopt, small_number);
644  BOOST_CHECK((bound.getParameter<eBoundTheta>() ==
646  bound.setParameter<eBoundTheta>(large_number);
647  BOOST_CHECK((bound.getParameter<eBoundTheta>() ==
649  bound.setParameter<eBoundTheta>(normal_number);
650  BOOST_CHECK((bound.getParameter<eBoundTheta>() == normal_number));
651 
652  // check cyclic parameter type
653  ParameterSet<BoundIndices, eBoundPhi> cyclic(std::nullopt, small_number);
654  // calculate expected results
657  // check that difference between original phi and stored phi is a multiple
658  // of the cyclic period
659  double multiple =
660  (cyclic.getParameter<eBoundPhi>() - small_number) / (max - min);
661  BOOST_CHECK(cyclic.getParameter<eBoundPhi>() >= min);
662  BOOST_CHECK(cyclic.getParameter<eBoundPhi>() < max);
663  BOOST_CHECK(std::abs(multiple - std::floor(multiple + 0.5)) < tol);
664 
665  cyclic.setParameter<eBoundPhi>(large_number);
666  multiple = (cyclic.getParameter<eBoundPhi>() - large_number) / (max - min);
667  BOOST_CHECK(cyclic.getParameter<eBoundPhi>() >= min);
668  BOOST_CHECK(cyclic.getParameter<eBoundPhi>() < max);
669  BOOST_CHECK(std::abs(multiple - std::floor(multiple + 0.5)) < tol);
670 
671  cyclic.setParameter<eBoundPhi>(normal_number);
672  multiple = (cyclic.getParameter<eBoundPhi>() - normal_number) / (max - min);
673  BOOST_CHECK(cyclic.getParameter<eBoundPhi>() >= min);
674  BOOST_CHECK(cyclic.getParameter<eBoundPhi>() < max);
675  BOOST_CHECK(std::abs(multiple - std::floor(multiple + 0.5)) < tol);
676 
677  // check residual calculation
678 
679  // input numbers
680  const double first_loc0 = 0.3;
681  const double first_phi = 0.9 * M_PI;
682  const double first_theta = 0.7 * M_PI;
683 
684  const double second_loc0 = 2.7;
685  const double second_phi = -0.9 * M_PI;
686  const double second_theta = 0.35 * M_PI;
687 
688  // expected results for residual second wrt first
689  const double delta_loc0 = second_loc0 - first_loc0;
690  const double delta_phi =
691  get_cyclic_difference(second_phi, first_phi, min, max);
692  const double delta_theta = second_theta - first_theta;
693  Vector3D residuals(delta_loc0, delta_phi, delta_theta);
694 
696  std::nullopt, first_loc0, first_phi, first_theta);
698  std::nullopt, second_loc0, second_phi, second_theta);
699  CHECK_CLOSE_REL(residuals, second.residual(first), 1e-6);
700 
701  // some more checks for bound variables
702  check_residuals_for_bound_parameters();
703 
704  // some more checks for cyclic variables
705  check_residuals_for_cyclic_parameters();
706 
707  // inspecific residual tests with random numbers
708  random_residual_tests();
709 }
710 
711 template <BoundIndices... params>
712 using ParSet = ParameterSet<BoundIndices, params...>;
713 
720 BOOST_AUTO_TEST_CASE(parset_parID_mapping) {
721  // check logic for type-based access
723  eBoundTime>::getIndex<eBoundLoc0>() == 0));
725  eBoundTime>::getIndex<eBoundLoc1>() == 1));
727  eBoundTime>::getIndex<eBoundPhi>() == 2));
729  eBoundTime>::getIndex<eBoundQOverP>() == 3));
731  eBoundTime>::getIndex<eBoundTime>() == 4));
732 
733  // check logic for index-based access
735  eBoundTime>::getParameterIndex<0>() == eBoundLoc0));
737  eBoundTime>::getParameterIndex<1>() == eBoundLoc1));
739  eBoundTime>::getParameterIndex<2>() == eBoundPhi));
741  eBoundTime>::getParameterIndex<3>() == eBoundQOverP));
743  eBoundTime>::getParameterIndex<4>() == eBoundTime));
744 
745  // check consistency
746  using FullSet = FullBoundParameterSet;
747  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<0>()>() == 0));
748  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<1>()>() == 1));
749  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<2>()>() == 2));
750  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<3>()>() == 3));
751  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<4>()>() == 4));
752  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<5>()>() == 5));
753 
754  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eBoundLoc0>()>() ==
755  eBoundLoc0));
756  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eBoundLoc1>()>() ==
757  eBoundLoc1));
758  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eBoundPhi>()>() ==
759  eBoundPhi));
760  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eBoundTheta>()>() ==
761  eBoundTheta));
762  BOOST_CHECK(
763  (FullSet::getParameterIndex<FullSet::getIndex<eBoundQOverP>()>() ==
764  eBoundQOverP));
765  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eBoundTime>()>() ==
766  eBoundTime));
767 
768  // consistency of types
769  BOOST_CHECK(
770  (std::is_same<std::remove_cv<decltype(
772  decltype(eBoundLoc0)>::value));
773  BOOST_CHECK((std::is_same<decltype(FullSet::getParameterIndex<0>()),
774  decltype(eBoundLoc0)>::value));
775 }
776 
790 BOOST_AUTO_TEST_CASE(free_parset_consistency_tests) {
791  // check template parameter based information
793 
794  // covariance matrix
796  cov << 1, 0, 0, 0, 1.2, 0.2, 0, 0.2, 0.7;
797 
798  // parameter values
799  double x = 0.5;
800  double y = -0.2;
801  double z = 0.3;
802  Vector3D parValues(x, y, z);
803 
804  // parameter set with covariance matrix
806  cov, x, y, z);
807 
808  // check number and type of stored parameters
809  BOOST_CHECK(parSet_with_cov.size() == 3);
810  BOOST_CHECK(parSet_with_cov.contains<eFreePos0>());
811  BOOST_CHECK(parSet_with_cov.contains<eFreePos1>());
812  BOOST_CHECK(parSet_with_cov.contains<eFreePos2>());
813  BOOST_CHECK(not parSet_with_cov.contains<eFreeTime>());
814  BOOST_CHECK(not parSet_with_cov.contains<eFreeDir0>());
815  BOOST_CHECK(not parSet_with_cov.contains<eFreeDir1>());
816  BOOST_CHECK(not parSet_with_cov.contains<eFreeDir2>());
817  BOOST_CHECK(not parSet_with_cov.contains<eFreeQOverP>());
818 
819  // check stored parameter values
820  BOOST_CHECK(parSet_with_cov.getParameter<eFreePos0>() == x);
821  BOOST_CHECK(parSet_with_cov.getParameter<eFreePos1>() == y);
822  BOOST_CHECK(parSet_with_cov.getParameter<eFreePos2>() == z);
823  BOOST_CHECK(parSet_with_cov.getParameters() == parValues);
824 
825  // check stored covariance
826  BOOST_CHECK(parSet_with_cov.getCovariance());
827  BOOST_CHECK(*parSet_with_cov.getCovariance() == cov);
828  BOOST_CHECK(parSet_with_cov.getUncertainty<eFreePos0>() == sqrt(cov(0, 0)));
829  BOOST_CHECK(parSet_with_cov.getUncertainty<eFreePos1>() == sqrt(cov(1, 1)));
830  BOOST_CHECK(parSet_with_cov.getUncertainty<eFreePos2>() == sqrt(cov(2, 2)));
831 
832  // same parameter set without covariance matrix
834  std::nullopt, parValues);
835 
836  BOOST_CHECK(!parSet_without_cov.getCovariance());
837  BOOST_CHECK(parSet_without_cov.getUncertainty<eFreePos0>() < 0);
838  BOOST_CHECK(parSet_without_cov.getUncertainty<eFreePos1>() < 0);
839  BOOST_CHECK(parSet_without_cov.getUncertainty<eFreePos2>() < 0);
840  BOOST_CHECK(parSet_without_cov.getParameters() ==
841  parSet_with_cov.getParameters());
842 
843  // set new covariance matrix
844  parSet_without_cov.setCovariance(cov);
845 
846  BOOST_CHECK(parSet_without_cov.getCovariance());
847  BOOST_CHECK(*parSet_without_cov.getCovariance() == cov);
848 
849  // set new parameter values
850  double newX = 0.1;
851  double newY = 0.6;
852  double newZ = -0.15;
853  parValues << newX, newY, newZ;
854  parSet_with_cov.setParameter<eFreePos0>(newX);
855  parSet_with_cov.setParameter<eFreePos1>(newY);
856  parSet_with_cov.setParameter<eFreePos2>(newZ);
857 
858  BOOST_CHECK(parSet_with_cov.getParameter<eFreePos0>() == newX);
859  BOOST_CHECK(parSet_with_cov.getParameter<eFreePos1>() == newY);
860  BOOST_CHECK(parSet_with_cov.getParameter<eFreePos2>() == newZ);
861  BOOST_CHECK(parSet_with_cov.getParameters() == parValues);
862 }
863 
871 BOOST_AUTO_TEST_CASE(free_parset_copy_assignment_tests) {
872  // covariance matrix
874  cov << 1, 0, 0, 0, 1.2, 0.2, 0, 0.2, 0.7;
875 
876  // parameter values
877  double x = 0.5;
878  double y = -0.2;
879  double z = 0.3;
880  Vector3D firstParValues(x, y, z);
881 
882  // parameter set with covariance matrix
884  z);
885 
886  // check copy constructor
888  BOOST_CHECK(first == copy);
889 
890  // check move constructor
892  std::move(copy));
893  BOOST_CHECK(first == moved);
894 
895  // check assignment operator
897  BOOST_CHECK(assigned == moved);
898 
900  std::nullopt, 0, 1.7, -0.15);
901  BOOST_CHECK(assigned != other);
902  assigned = other;
903  BOOST_CHECK(assigned == other);
904 
905  // check move assignment
906  BOOST_CHECK(first != assigned);
908  BOOST_CHECK(first == assigned);
909 
910  // check swap method
913  std::nullopt, 2 * x, 2 * y, 2 * z);
916 
917  BOOST_CHECK(lhs != rhs && lhs == lhs_copy && rhs == rhs_copy);
918 }
919 
925 BOOST_AUTO_TEST_CASE(free_parset_comparison_tests) {
926  // covariance matrix
928  cov << 1, 0, 0, 0, 1.2, 0.2, 0, 0.2, 0.7;
929 
930  // parameter values
931  double x = 0.5;
932  double y = -0.2;
933  double z = 0.3;
934 
935  // parameter set with covariance matrix
937  z);
939  std::nullopt, 2 * x, 2 * y, 2 * z);
940 
941  // check self comparison
942  BOOST_CHECK(first == first);
943  BOOST_CHECK(not(first != first));
944 
945  // check mutual exclusivity
946  BOOST_CHECK(first != second);
947  BOOST_CHECK(not(first == second));
948  first = second;
949  BOOST_CHECK(first == second);
950 
951  // check that comparison fails for unequal parameter values
952  second.setParameter<eFreePos0>(3 * x);
953  BOOST_CHECK(first != second);
954  first = second;
955  BOOST_CHECK(first == second);
956 
957  second.setParameter<eFreePos1>(3 * y);
958  BOOST_CHECK(first != second);
959  first = second;
960  BOOST_CHECK(first == second);
961 
962  second.setParameter<eFreePos2>(3 * z);
963  BOOST_CHECK(first != second);
964  first = second;
965  BOOST_CHECK(first == second);
966 
967  // check that comparison fails for unequal covariance matrices
968  second.setCovariance(cov);
969  BOOST_CHECK(first != second);
970  first = second;
971  BOOST_CHECK(first == second);
972 
973  cov(0, 0) = 2 * cov(0, 0);
974  second.setCovariance(cov);
975  BOOST_CHECK(first != second);
976  first = second;
977  BOOST_CHECK(first == second);
978 }
979 
989 BOOST_AUTO_TEST_CASE(free_parset_projection_tests) {
990  // clang-format off
992  z_proj << 0, 0, 1, 0, 0, 0, 0, 0;
993 
994  ActsMatrixD<2, eFreeSize> x_qop_proj;
995  x_qop_proj << 1, 0, 0, 0, 0, 0, 0, 0,
996  0, 0, 0, 0, 0, 0, 0, 1;
997 
998  ActsMatrixD<2, eFreeSize> y_tz_proj;
999  y_tz_proj << 0, 1, 0, 0, 0, 0, 0, 0,
1000  0, 0, 0, 0, 0, 0, 1, 0;
1001 
1002  ActsMatrixD<3, eFreeSize> x_y_z_proj;
1003  x_y_z_proj << 1, 0, 0, 0, 0, 0, 0, 0,
1004  0, 1, 0, 0, 0, 0, 0, 0,
1005  0, 0, 1, 0, 0, 0, 0, 0;
1006 
1007  ActsMatrixD<4, eFreeSize> x_z_tz_qop_proj;
1008  x_z_tz_qop_proj << 1, 0, 0, 0, 0, 0, 0, 0,
1009  0, 0, 1, 0, 0, 0, 0, 0,
1010  0, 0, 0, 0, 0, 0, 1, 0,
1011  0, 0, 0, 0, 0, 0, 0, 1;
1012 
1013  ActsMatrixD<5, eFreeSize> x_y_z_tz_qop_proj;
1014  x_y_z_tz_qop_proj << 1, 0, 0, 0, 0, 0, 0, 0,
1015  0, 1, 0, 0, 0, 0, 0, 0,
1016  0, 0, 1, 0, 0, 0, 0, 0,
1017  0, 0, 0, 0, 0, 0, 1, 0,
1018  0, 0, 0, 0, 0, 0, 0, 1;
1019 
1020  ActsMatrixD<6, eFreeSize> x_y_z_t_tz_qop_proj;
1021  x_y_z_t_tz_qop_proj << 1, 0, 0, 0, 0, 0, 0, 0,
1022  0, 1, 0, 0, 0, 0, 0, 0,
1023  0, 0, 1, 0, 0, 0, 0, 0,
1024  0, 0, 0, 1, 0, 0, 0, 0,
1025  0, 0, 0, 0, 0, 0, 1, 0,
1026  0, 0, 0, 0, 0, 0, 0, 1;
1027 
1028  ActsMatrixD<7, eFreeSize> x_y_z_t_ty_tz_qop_proj;
1029  x_y_z_t_ty_tz_qop_proj << 1, 0, 0, 0, 0, 0, 0, 0,
1030  0, 1, 0, 0, 0, 0, 0, 0,
1031  0, 0, 1, 0, 0, 0, 0, 0,
1032  0, 0, 0, 1, 0, 0, 0, 0,
1033  0, 0, 0, 0, 0, 1, 0, 0,
1034  0, 0, 0, 0, 0, 0, 1, 0,
1035  0, 0, 0, 0, 0, 0, 0, 1;
1036 
1038  x_y_z_t_tx_ty_tz_qop_proj;
1039  x_y_z_t_tx_ty_tz_qop_proj << 1, 0, 0, 0, 0, 0, 0, 0,
1040  0, 1, 0, 0, 0, 0, 0, 0,
1041  0, 0, 1, 0, 0, 0, 0, 0,
1042  0, 0, 0, 1, 0, 0, 0, 0,
1043  0, 0, 0, 0, 1, 0, 0, 0,
1044  0, 0, 0, 0, 0, 1, 0, 0,
1045  0, 0, 0, 0, 0, 0, 1, 0,
1046  0, 0, 0, 0, 0, 0, 0, 1;
1047  // clang-format on
1048 
1049  BOOST_CHECK((ParameterSet<FreeIndices, eFreePos2>::projector() == z_proj));
1051  x_qop_proj));
1053  y_tz_proj));
1054  BOOST_CHECK((
1056  x_y_z_proj));
1058  eFreeQOverP>::projector() == x_z_tz_qop_proj));
1059  BOOST_CHECK(
1061  eFreeQOverP>::projector() == x_y_z_tz_qop_proj));
1063  eFreeTime, eFreeDir2, eFreeQOverP>::projector() ==
1064  x_y_z_t_tz_qop_proj));
1065  BOOST_CHECK(
1067  eFreeDir1, eFreeDir2, eFreeQOverP>::projector() ==
1068  x_y_z_t_ty_tz_qop_proj));
1069  BOOST_CHECK((
1071  eFreeDir0, eFreeDir1, eFreeDir2, eFreeQOverP>::projector() ==
1072  x_y_z_t_tx_ty_tz_qop_proj));
1073 }
1074 
1085 BOOST_AUTO_TEST_CASE(free_parset_residual_tests) {
1086  // check unbound parameter type
1087  const double large_number = 12443534120;
1088  const double small_number = -924342675;
1089  const double normal_number = 0.1234;
1091  std::nullopt, small_number, large_number, normal_number);
1092  BOOST_CHECK(unbound.getParameter<eFreePos0>() == small_number);
1093  BOOST_CHECK(unbound.getParameter<eFreePos1>() == large_number);
1094  BOOST_CHECK(unbound.getParameter<eFreeQOverP>() == normal_number);
1095 
1096  // check residual calculation
1097  // input numbers
1098  const double first_x = 0.3;
1099  const double first_y = 0.9;
1100  const double first_z = 0.7;
1101 
1102  const double second_x = 2.7;
1103  const double second_y = -0.9;
1104  const double second_z = 0.35;
1105 
1106  // expected results for residual second wrt first
1107  const double delta_x = second_x - first_x;
1108  const double delta_y = second_y - first_y;
1109  const double delta_z = second_z - first_z;
1110  Vector3D residuals(delta_x, delta_y, delta_z);
1111 
1113  std::nullopt, first_x, first_y, first_z);
1115  std::nullopt, second_x, second_y, second_z);
1116  CHECK_CLOSE_REL(residuals, second.residual(first), 1e-6);
1117 
1118  // inspecific residual tests with random numbers
1119  free_random_residual_tests();
1120 }
1121 
1122 template <FreeIndices... params>
1124 
1131 BOOST_AUTO_TEST_CASE(free_parset_parID_mapping) {
1132  // check logic for type-based access
1133  BOOST_CHECK(
1135  eFreeDir1, eFreeDir2, eFreeQOverP>::getIndex<eFreePos0>() ==
1136  0));
1137  BOOST_CHECK(
1139  eFreeDir1, eFreeDir2, eFreeQOverP>::getIndex<eFreePos1>() ==
1140  1));
1141  BOOST_CHECK(
1143  eFreeDir1, eFreeDir2, eFreeQOverP>::getIndex<eFreePos2>() ==
1144  2));
1145  BOOST_CHECK(
1147  eFreeDir1, eFreeDir2, eFreeQOverP>::getIndex<eFreeTime>() ==
1148  3));
1149  BOOST_CHECK(
1151  eFreeDir1, eFreeDir2, eFreeQOverP>::getIndex<eFreeDir0>() ==
1152  4));
1153  BOOST_CHECK(
1155  eFreeDir1, eFreeDir2, eFreeQOverP>::getIndex<eFreeDir1>() ==
1156  5));
1157  BOOST_CHECK(
1159  eFreeDir1, eFreeDir2, eFreeQOverP>::getIndex<eFreeDir2>() ==
1160  6));
1161  BOOST_CHECK(
1163  eFreeDir1, eFreeDir2, eFreeQOverP>::getIndex<eFreeQOverP>() ==
1164  7));
1165 
1166  // check logic for index-based access
1167  BOOST_CHECK(
1169  eFreeDir1, eFreeDir2, eFreeQOverP>::getParameterIndex<0>() ==
1170  eFreePos0));
1171  BOOST_CHECK(
1173  eFreeDir1, eFreeDir2, eFreeQOverP>::getParameterIndex<1>() ==
1174  eFreePos1));
1175  BOOST_CHECK(
1177  eFreeDir1, eFreeDir2, eFreeQOverP>::getParameterIndex<2>() ==
1178  eFreePos2));
1179  BOOST_CHECK(
1181  eFreeDir1, eFreeDir2, eFreeQOverP>::getParameterIndex<3>() ==
1182  eFreeTime));
1183  BOOST_CHECK(
1185  eFreeDir1, eFreeDir2, eFreeQOverP>::getParameterIndex<4>() ==
1186  eFreeDir0));
1187  BOOST_CHECK(
1189  eFreeDir1, eFreeDir2, eFreeQOverP>::getParameterIndex<5>() ==
1190  eFreeDir1));
1191  BOOST_CHECK(
1193  eFreeDir1, eFreeDir2, eFreeQOverP>::getParameterIndex<6>() ==
1194  eFreeDir2));
1195  BOOST_CHECK(
1197  eFreeDir1, eFreeDir2, eFreeQOverP>::getParameterIndex<7>() ==
1198  eFreeQOverP));
1199 
1200  // check consistency
1201  using FullSet = FullFreeParameterSet;
1202  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<0>()>() == 0));
1203  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<1>()>() == 1));
1204  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<2>()>() == 2));
1205  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<3>()>() == 3));
1206  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<4>()>() == 4));
1207  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<5>()>() == 5));
1208  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<6>()>() == 6));
1209  BOOST_CHECK((FullSet::getIndex<FullSet::getParameterIndex<7>()>() == 7));
1210 
1211  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eFreePos0>()>() ==
1212  eFreePos0));
1213  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eFreePos1>()>() ==
1214  eFreePos1));
1215  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eFreePos2>()>() ==
1216  eFreePos2));
1217  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eFreeTime>()>() ==
1218  eFreeTime));
1219  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eFreeDir0>()>() ==
1220  eFreeDir0));
1221  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eFreeDir1>()>() ==
1222  eFreeDir1));
1223  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eFreeDir2>()>() ==
1224  eFreeDir2));
1225  BOOST_CHECK((FullSet::getParameterIndex<FullSet::getIndex<eFreeQOverP>()>() ==
1226  eFreeQOverP));
1227 
1228  // consistency of types
1229  BOOST_CHECK(
1230  (std::is_same<std::remove_cv<decltype(
1232  decltype(eFreePos0)>::value));
1233  BOOST_CHECK((std::is_same<decltype(FullSet::getParameterIndex<0>()),
1234  decltype(eFreePos0)>::value));
1235 }
1236 } // namespace Test
1237 } // namespace Acts