EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ExtrapolatorTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ExtrapolatorTests.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2018-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 
13 #include "Acts/EventData/TrackParameters.hpp"
14 #include "Acts/Geometry/GeometryContext.hpp"
15 #include "Acts/MagneticField/ConstantBField.hpp"
16 #include "Acts/MagneticField/MagneticFieldContext.hpp"
17 #include "Acts/Material/Material.hpp"
18 #include "Acts/Propagator/ActionList.hpp"
19 #include "Acts/Propagator/EigenStepper.hpp"
20 #include "Acts/Propagator/MaterialInteractor.hpp"
21 #include "Acts/Propagator/Navigator.hpp"
22 #include "Acts/Propagator/Propagator.hpp"
23 #include "Acts/Propagator/SurfaceCollector.hpp"
24 #include "Acts/Surfaces/CylinderSurface.hpp"
25 #include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp"
26 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
27 #include "Acts/Utilities/Definitions.hpp"
28 #include "Acts/Utilities/Units.hpp"
29 
30 namespace bdata = boost::unit_test::data;
31 namespace tt = boost::test_tools;
32 using namespace Acts::UnitLiterals;
33 
34 namespace Acts {
35 namespace Test {
36 
37 // Create a test context
38 GeometryContext tgContext = GeometryContext();
39 MagneticFieldContext mfContext = MagneticFieldContext();
40 
41 // Global definitions
42 // The path limit abort
43 using path_limit = PathLimitReached;
44 
45 std::vector<std::unique_ptr<const Surface>> stepState;
46 
47 CylindricalTrackingGeometry cGeometry(tgContext);
48 auto tGeometry = cGeometry();
49 
50 // Get the navigator and provide the TrackingGeometry
51 Navigator navigator(tGeometry);
52 
53 using BFieldType = ConstantBField;
54 using EigenStepperType = EigenStepper<BFieldType>;
55 using EigenPropagatorType = Propagator<EigenStepperType, Navigator>;
56 using Covariance = BoundSymMatrix;
57 
58 BFieldType bField(0, 0, 2_T);
59 EigenStepperType estepper(bField);
60 EigenPropagatorType epropagator(std::move(estepper), std::move(navigator));
61 
62 const int ntests = 100;
63 
64 // A plane selector for the SurfaceCollector
65 struct PlaneSelector {
68  bool operator()(const Surface& sf) const {
69  return (sf.type() == Surface::Plane);
70  }
71 };
72 
73 // This test case checks that no segmentation fault appears
74 // - simple extrapolation test
75 BOOST_DATA_TEST_CASE(
76  test_extrapolation_,
77  bdata::random((bdata::seed = 0,
78  bdata::distribution =
79  std::uniform_real_distribution<>(0.4_GeV, 10_GeV))) ^
80  bdata::random((bdata::seed = 1,
81  bdata::distribution =
82  std::uniform_real_distribution<>(-M_PI, M_PI))) ^
83  bdata::random((bdata::seed = 2,
84  bdata::distribution =
85  std::uniform_real_distribution<>(1.0, M_PI - 1.0))) ^
86  bdata::random(
87  (bdata::seed = 3,
88  bdata::distribution = std::uniform_int_distribution<>(0, 1))) ^
89  bdata::random(
90  (bdata::seed = 4,
91  bdata::distribution = std::uniform_int_distribution<>(0, 100))) ^
92  bdata::xrange(ntests),
93  pT, phi, theta, charge, time, index) {
94  double p = pT / sin(theta);
95  double q = -1 + 2 * charge;
96  (void)index;
97 
98  // define start parameters
100  Covariance cov;
101  // take some major correlations (off-diagonals)
102  cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
103  0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
104  0, 0;
105  CurvilinearTrackParameters start(Vector4D(0, 0, 0, time), phi, theta, p, q,
106  cov);
107 
108  PropagatorOptions<> options(tgContext, mfContext, getDummyLogger());
109  options.maxStepSize = 10_cm;
110  options.pathLimit = 25_cm;
111 
112  BOOST_CHECK(epropagator.propagate(start, options).value().endParameters !=
113  nullptr);
114 }
115 
116 // This test case checks that no segmentation fault appears
117 // - this tests the collection of surfaces
118 BOOST_DATA_TEST_CASE(
119  test_surface_collection_,
120  bdata::random((bdata::seed = 10,
121  bdata::distribution =
122  std::uniform_real_distribution<>(0.4_GeV, 10_GeV))) ^
123  bdata::random((bdata::seed = 11,
124  bdata::distribution =
125  std::uniform_real_distribution<>(-M_PI, M_PI))) ^
126  bdata::random((bdata::seed = 12,
127  bdata::distribution =
128  std::uniform_real_distribution<>(1.0, M_PI - 1.0))) ^
129  bdata::random(
130  (bdata::seed = 13,
131  bdata::distribution = std::uniform_int_distribution<>(0, 1))) ^
132  bdata::random(
133  (bdata::seed = 14,
134  bdata::distribution = std::uniform_int_distribution<>(0, 100))) ^
135  bdata::xrange(ntests),
136  pT, phi, theta, charge, time, index) {
137  double p = pT / sin(theta);
138  double q = -1 + 2 * charge;
139  (void)index;
140 
141  // define start parameters
143  Covariance cov;
144  // take some major correlations (off-diagonals)
145  cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
146  0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
147  0, 0;
148  CurvilinearTrackParameters start(Vector4D(0, 0, 0, time), phi, theta, p, q,
149  cov);
150 
151  // A PlaneSelector for the SurfaceCollector
152  using PlaneCollector = SurfaceCollector<PlaneSelector>;
153 
154  PropagatorOptions<ActionList<PlaneCollector>> options(tgContext, mfContext,
155  getDummyLogger());
156 
157  options.maxStepSize = 10_cm;
158  options.pathLimit = 25_cm;
159 
160  const auto& result = epropagator.propagate(start, options).value();
161  auto collector_result = result.get<PlaneCollector::result_type>();
162 
163  // step through the surfaces and go step by step
164  PropagatorOptions<> optionsEmpty(tgContext, mfContext, getDummyLogger());
165 
166  optionsEmpty.maxStepSize = 25_cm;
167  // Try propagation from start to each surface
168  for (const auto& colsf : collector_result.collected) {
169  const auto& csurface = colsf.surface;
170  // Avoid going to the same surface
171  // @todo: decide on strategy and write unit test for this
172  if (csurface == &start.referenceSurface()) {
173  continue;
174  }
175  // Extrapolate & check
176  const auto& cresult = epropagator.propagate(start, *csurface, optionsEmpty)
177  .value()
178  .endParameters;
179  BOOST_CHECK(cresult != nullptr);
180  }
181 }
182 
183 // This test case checks that no segmentation fault appears
184 // - this tests the collection of surfaces
185 BOOST_DATA_TEST_CASE(
186  test_material_interactor_,
187  bdata::random((bdata::seed = 20,
188  bdata::distribution =
189  std::uniform_real_distribution<>(0.4_GeV, 10_GeV))) ^
190  bdata::random((bdata::seed = 21,
191  bdata::distribution =
192  std::uniform_real_distribution<>(-M_PI, M_PI))) ^
193  bdata::random((bdata::seed = 22,
194  bdata::distribution =
195  std::uniform_real_distribution<>(1.0, M_PI - 1.0))) ^
196  bdata::random(
197  (bdata::seed = 23,
198  bdata::distribution = std::uniform_int_distribution<>(0, 1))) ^
199  bdata::random(
200  (bdata::seed = 24,
201  bdata::distribution = std::uniform_int_distribution<>(0, 100))) ^
202  bdata::xrange(ntests),
203  pT, phi, theta, charge, time, index) {
204  double p = pT / sin(theta);
205  double q = -1 + 2 * charge;
206  (void)index;
207 
208  // define start parameters
210  Covariance cov;
211  // take some major correlations (off-diagonals)
212  cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
213  0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
214  0, 0;
215  CurvilinearTrackParameters start(Vector4D(0, 0, 0, time), phi, theta, p, q,
216  cov);
217 
218  PropagatorOptions<ActionList<MaterialInteractor>> options(
219  tgContext, mfContext, getDummyLogger());
220  options.maxStepSize = 25_cm;
221  options.pathLimit = 25_cm;
222 
223  const auto& result = epropagator.propagate(start, options).value();
224  if (result.endParameters) {
225  // test that you actually lost some energy
226  BOOST_CHECK_LT(result.endParameters->absoluteMomentum(),
227  start.absoluteMomentum());
228  }
229 }
230 
231 // This test case checks that no segmentation fault appears
232 // - this tests the loop protection
233 BOOST_DATA_TEST_CASE(
234  loop_protection_test,
235  bdata::random((bdata::seed = 20,
236  bdata::distribution =
237  std::uniform_real_distribution<>(0.1_GeV, 0.5_GeV))) ^
238  bdata::random((bdata::seed = 21,
239  bdata::distribution =
240  std::uniform_real_distribution<>(-M_PI, M_PI))) ^
241  bdata::random((bdata::seed = 22,
242  bdata::distribution =
243  std::uniform_real_distribution<>(1.0, M_PI - 1.0))) ^
244  bdata::random(
245  (bdata::seed = 23,
246  bdata::distribution = std::uniform_int_distribution<>(0, 1))) ^
247  bdata::random(
248  (bdata::seed = 24,
249  bdata::distribution = std::uniform_int_distribution<>(0, 100))) ^
250  bdata::xrange(ntests),
251  pT, phi, theta, charge, time, index) {
252  double p = pT / sin(theta);
253  double q = -1 + 2 * charge;
254  (void)index;
255 
256  // define start parameters
258  Covariance cov;
259  // take some major correlations (off-diagonals)
260  cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
261  0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
262  0, 0;
263  CurvilinearTrackParameters start(Vector4D(0, 0, 0, time), phi, theta, p, q,
264  cov);
265 
266  // Action list and abort list
267  PropagatorOptions<ActionList<MaterialInteractor>> options(
268  tgContext, mfContext, getDummyLogger());
269  options.maxStepSize = 25_cm;
270  options.pathLimit = 1500_mm;
271 
272  const auto& status = epropagator.propagate(start, options).value();
273  // this test assumes state.options.loopFraction = 0.5
274  // maximum momentum allowed
275  double pmax = options.pathLimit *
276  bField.getField(start.position(tgContext)).norm() / M_PI;
277  if (p < pmax) {
278  BOOST_CHECK_LT(status.pathLength, options.pathLimit);
279  } else {
280  CHECK_CLOSE_REL(status.pathLength, options.pathLimit, 1e-3);
281  }
282 }
283 
284 } // namespace Test
285 } // namespace Acts