ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PropagatorTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PropagatorTests.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2016-2019 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/data/test_case.hpp>
10 #include <boost/test/tools/output_test_stream.hpp>
11 #include <boost/test/unit_test.hpp>
12 
25 #include "Acts/Utilities/Units.hpp"
26 
27 namespace bdata = boost::unit_test::data;
28 namespace tt = boost::test_tools;
29 using namespace Acts::UnitLiterals;
31 
32 namespace Acts {
33 namespace Test {
34 
35 // Create a test context
38 
40 
44  struct this_result {
45  double distance = std::numeric_limits<double>::max();
46  };
47 
49 
50  PerpendicularMeasure() = default;
51 
52  template <typename propagator_state_t, typename stepper_t>
53  void operator()(propagator_state_t& state, const stepper_t& stepper,
54  result_type& result) const {
55  result.distance = perp(stepper.position(state.stepping));
56  }
57 
58  template <typename propagator_state_t, typename stepper_t>
59  void operator()(propagator_state_t& /*unused*/,
60  const stepper_t& /*unused*/) const {}
61 };
62 
64 template <typename Surface>
66  // the surface to be intersected
67  const Surface* surface = nullptr;
68  // the tolerance for intersection
69  double tolerance = 1.e-5;
70 
72  struct this_result {
73  size_t surfaces_passed = 0;
74  double surface_passed_r = std::numeric_limits<double>::max();
75  };
76 
78 
79  SurfaceObserver() = default;
80 
81  template <typename propagator_state_t, typename stepper_t>
82  void operator()(propagator_state_t& state, const stepper_t& stepper,
83  result_type& result) const {
84  if (surface && !result.surfaces_passed) {
85  // calculate the distance to the surface
86  const double distance =
87  surface
88  ->intersectionEstimate(state.geoContext,
89  stepper.position(state.stepping),
90  stepper.direction(state.stepping), true)
91  .pathLength;
92  // Adjust the step size so that we cannot cross the target surface
93  state.stepping.stepSize.update(distance, ConstrainedStep::actor);
94  // return true if you fall below tolerance
95  if (std::abs(distance) <= tolerance) {
96  ++result.surfaces_passed;
97  result.surface_passed_r = perp(stepper.position(state.stepping));
98  // release the step size, will be re-adjusted
99  state.stepping.stepSize.release(ConstrainedStep::actor);
100  }
101  }
102  }
103 
104  template <typename propagator_state_t, typename stepper_t>
105  void operator()(propagator_state_t& /*unused*/,
106  const stepper_t& /*unused*/) const {}
107 };
108 
109 // Global definitions
110 // The path limit abort
112 
113 using BFieldType = ConstantBField;
116 
117 const double Bz = 2_T;
118 BFieldType bField(0, 0, Bz);
121 
122 auto mCylinder = std::make_shared<CylinderBounds>(10_mm, 1000_mm);
123 auto mSurface = Surface::makeShared<CylinderSurface>(nullptr, mCylinder);
124 auto cCylinder = std::make_shared<CylinderBounds>(150_mm, 1000_mm);
125 auto cSurface = Surface::makeShared<CylinderSurface>(nullptr, cCylinder);
126 
127 const int ntests = 5;
128 
129 // This tests the Options
130 BOOST_AUTO_TEST_CASE(PropagatorOptions_) {
131  using null_optionsType = PropagatorOptions<>;
132  null_optionsType null_options(tgContext, mfContext);
133  // todo write null options test
134 
135  using ActionListType = ActionList<PerpendicularMeasure>;
136  using AbortConditionsType = AbortList<>;
137 
139 
140  optionsType options(tgContext, mfContext);
141 }
142 
144  cylinder_passage_observer_,
145  bdata::random((bdata::seed = 0,
146  bdata::distribution =
147  std::uniform_real_distribution<>(0.4_GeV, 10_GeV))) ^
148  bdata::random((bdata::seed = 1,
149  bdata::distribution =
150  std::uniform_real_distribution<>(-M_PI, M_PI))) ^
151  bdata::random((bdata::seed = 2,
152  bdata::distribution =
153  std::uniform_real_distribution<>(1.0, M_PI - 1.0))) ^
154  bdata::random(
155  (bdata::seed = 3,
156  bdata::distribution = std::uniform_int_distribution<>(0, 1))) ^
157  bdata::random(
158  (bdata::seed = 4,
159  bdata::distribution = std::uniform_int_distribution<>(0, 100))) ^
160  bdata::xrange(ntests),
161  pT, phi, theta, charge, time, index) {
162  double dcharge = -1 + 2 * charge;
163  (void)index;
164 
165  using CylinderObserver = SurfaceObserver<CylinderSurface>;
166  using ActionListType = ActionList<CylinderObserver>;
167  using AbortConditionsType = AbortList<>;
168 
169  // setup propagation options
171  mfContext);
172 
173  options.pathLimit = 20_m;
174  options.maxStepSize = 1_cm;
175 
176  // set the surface to be passed by
177  options.actionList.get<CylinderObserver>().surface = mSurface.get();
178 
179  using so_result = typename CylinderObserver::result_type;
180 
181  // define start parameters
182  double x = 0;
183  double y = 0;
184  double z = 0;
185  double px = pT * cos(phi);
186  double py = pT * sin(phi);
187  double pz = pT / tan(theta);
188  double q = dcharge;
189  Vector3D pos(x, y, z);
190  Vector3D mom(px, py, pz);
191  CurvilinearParameters start(std::nullopt, pos, mom, q, time);
192  // propagate to the cylinder surface
193  const auto& result = epropagator.propagate(start, *cSurface, options).value();
194  auto& sor = result.get<so_result>();
195 
196  BOOST_CHECK_EQUAL(sor.surfaces_passed, 1u);
197  CHECK_CLOSE_ABS(sor.surface_passed_r, 10., 1e-5);
198 }
199 
201  curvilinear_additive_,
202  bdata::random((bdata::seed = 0,
203  bdata::distribution =
204  std::uniform_real_distribution<>(0.4_GeV, 10_GeV))) ^
205  bdata::random((bdata::seed = 1,
206  bdata::distribution =
207  std::uniform_real_distribution<>(-M_PI, M_PI))) ^
208  bdata::random((bdata::seed = 2,
209  bdata::distribution =
210  std::uniform_real_distribution<>(1.0, M_PI - 1.0))) ^
211  bdata::random(
212  (bdata::seed = 3,
213  bdata::distribution = std::uniform_int_distribution<>(0, 1))) ^
214  bdata::random(
215  (bdata::seed = 4,
216  bdata::distribution = std::uniform_int_distribution<>(0, 100))) ^
217  bdata::xrange(ntests),
218  pT, phi, theta, charge, time, index) {
219  double dcharge = -1 + 2 * charge;
220  (void)index;
221 
222  // setup propagation options - the tow step options
224  options_2s.pathLimit = 50_cm;
225  options_2s.maxStepSize = 1_cm;
226 
227  // define start parameters
228  double x = 0;
229  double y = 0;
230  double z = 0;
231  double px = pT * cos(phi);
232  double py = pT * sin(phi);
233  double pz = pT / tan(theta);
234  double q = dcharge;
235  Vector3D pos(x, y, z);
236  Vector3D mom(px, py, pz);
238  Covariance cov;
239  // take some major correlations (off-diagonals)
240  cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
241  0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
242  0, 0;
243  CurvilinearParameters start(cov, pos, mom, q, time);
244  // propagate to a path length of 100 with two steps of 50
245  const auto& mid_parameters =
246  epropagator.propagate(start, options_2s).value().endParameters;
247  const auto& end_parameters_2s =
248  epropagator.propagate(*mid_parameters, options_2s).value().endParameters;
249 
250  // setup propagation options - the one step options
252  options_1s.pathLimit = 100_cm;
253  options_1s.maxStepSize = 1_cm;
254  // propagate to a path length of 100 in one step
255  const auto& end_parameters_1s =
256  epropagator.propagate(start, options_1s).value().endParameters;
257 
258  // test that the propagation is additive
259  CHECK_CLOSE_REL(end_parameters_1s->position(), end_parameters_2s->position(),
260  0.001);
261 
262  const auto& cov_1s = *(end_parameters_1s->covariance());
263  const auto& cov_2s = *(end_parameters_2s->covariance());
264 
265  // CHECK_CLOSE_COVARIANCE(cov_1s, cov_2s, 0.001);
266  for (unsigned int i = 0; i < cov_1s.rows(); i++) {
267  for (unsigned int j = 0; j < cov_1s.cols(); j++) {
268  CHECK_CLOSE_OR_SMALL(cov_1s(i, j), cov_2s(i, j), 0.001, 1e-6);
269  }
270  }
271 }
272 
274  cylinder_additive_,
275  bdata::random((bdata::seed = 0,
276  bdata::distribution =
277  std::uniform_real_distribution<>(0.4_GeV, 10_GeV))) ^
278  bdata::random((bdata::seed = 1,
279  bdata::distribution =
280  std::uniform_real_distribution<>(-M_PI, M_PI))) ^
281  bdata::random((bdata::seed = 2,
282  bdata::distribution =
283  std::uniform_real_distribution<>(1.0, M_PI - 1.0))) ^
284  bdata::random(
285  (bdata::seed = 3,
286  bdata::distribution = std::uniform_int_distribution<>(0, 1))) ^
287  bdata::random(
288  (bdata::seed = 4,
289  bdata::distribution = std::uniform_int_distribution<>(0, 100))) ^
290  bdata::xrange(ntests),
291  pT, phi, theta, charge, time, index) {
292  double dcharge = -1 + 2 * charge;
293  (void)index;
294 
295  // setup propagation options - 2 setp options
297  options_2s.pathLimit = 10_m;
298  options_2s.maxStepSize = 1_cm;
299 
300  // define start parameters
301  double x = 0;
302  double y = 0;
303  double z = 0;
304  double px = pT * cos(phi);
305  double py = pT * sin(phi);
306  double pz = pT / tan(theta);
307  double q = dcharge;
308  Vector3D pos(x, y, z);
309  Vector3D mom(px, py, pz);
311  Covariance cov;
312  // take some major correlations (off-diagonals)
313  cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
314  0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
315  0, 0;
316  CurvilinearParameters start(cov, pos, mom, q, time);
317  // propagate to a final surface with one stop in between
318  const auto& mid_parameters =
319  epropagator.propagate(start, *mSurface, options_2s).value().endParameters;
320 
321  const auto& end_parameters_2s =
322  epropagator.propagate(*mid_parameters, *cSurface, options_2s)
323  .value()
324  .endParameters;
325 
326  // setup propagation options - one step options
328  options_1s.pathLimit = 10_m;
329  options_1s.maxStepSize = 1_cm;
330  // propagate to a final surface in one stop
331  const auto& end_parameters_1s =
332  epropagator.propagate(start, *cSurface, options_1s).value().endParameters;
333 
334  // test that the propagation is additive
335  CHECK_CLOSE_REL(end_parameters_1s->position(), end_parameters_2s->position(),
336  0.001);
337 
338  const auto& cov_1s = (*(end_parameters_1s->covariance()));
339  const auto& cov_2s = (*(end_parameters_2s->covariance()));
340 
341  // CHECK_CLOSE_COVARIANCE(cov_1s, cov_2s, 0.001);
342  for (unsigned int i = 0; i < cov_1s.rows(); i++) {
343  for (unsigned int j = 0; j < cov_1s.cols(); j++) {
344  CHECK_CLOSE_OR_SMALL(cov_1s(i, j), cov_2s(i, j), 0.001, 1e-6);
345  }
346  }
347 }
348 
349 } // namespace Test
350 } // namespace Acts