Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KalmanExtrapolatorTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KalmanExtrapolatorTests.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/unit_test.hpp>
10 
28 
29 #include <algorithm>
30 #include <array>
31 #include <map>
32 #include <memory>
33 #include <optional>
34 #include <tuple>
35 #include <utility>
36 #include <vector>
37 
38 namespace Acts {
39 class Logger;
40 struct EndOfWorldReached;
41 } // namespace Acts
42 
43 using namespace Acts::UnitLiterals;
44 
45 namespace Acts {
46 namespace Test {
47 
48 using Jacobian = BoundMatrix;
50 
51 // Create a test context
57 struct StepWiseActor {
59  struct this_result {
60  std::vector<Jacobian> jacobians = {};
61  std::vector<double> paths = {};
62 
63  double fullPath = 0.;
64 
65  bool finalized = false;
66  };
69 
80  template <typename propagator_state_t, typename stepper_t,
81  typename navigator_t>
82  void operator()(propagator_state_t& state, const stepper_t& stepper,
83  const navigator_t& navigator, result_type& result,
84  const Logger& /*logger*/) const {
85  // Listen to the surface and create bound state where necessary
86  auto surface = navigator.currentSurface(state.navigation);
87  if (surface and surface->associatedDetectorElement()) {
88  // Create a bound state and log the jacobian
89  auto boundState = stepper.boundState(state.stepping, *surface).value();
90  result.jacobians.push_back(std::move(std::get<Jacobian>(boundState)));
91  result.paths.push_back(std::get<double>(boundState));
92  }
93  // Also store the jacobian and full path
94  if ((navigator.navigationBreak(state.navigation) or
95  navigator.targetReached(state.navigation)) and
96  not result.finalized) {
97  // Set the last stepping parameter
98  result.paths.push_back(state.stepping.pathAccumulated);
99  // Set the full parameter
100  result.fullPath = state.stepping.pathAccumulated;
101  // Remember that you finalized this
102  result.finalized = true;
103  }
104  }
105 };
106 
110 BOOST_AUTO_TEST_CASE(kalman_extrapolator) {
111  // Build detector, take the cubic detector
113  auto detector = cGeometry();
114 
115  // The Navigator through the detector geometry
117  cfg.resolvePassive = false;
118  cfg.resolveMaterial = true;
119  cfg.resolveSensitive = true;
121 
122  // Configure propagation with deactivated B-field
123  auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 0.));
124  using Stepper = EigenStepper<>;
125  Stepper stepper(bField);
127  Propagator propagator(stepper, navigator);
128 
129  // Set initial parameters for the particle track
130  Covariance cov;
131  cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
132  0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
133  0, 0;
134  // The start parameters
135  CurvilinearTrackParameters start(Vector4(-3_m, 0, 0, 42_ns), 0_degree,
136  90_degree, 1_e / 1_GeV, cov,
138 
139  // Create the ActionList and AbortList
140  using StepWiseResult = StepWiseActor::result_type;
141  using StepWiseActors = ActionList<StepWiseActor>;
142  using Aborters = AbortList<EndOfWorldReached>;
143 
144  // Create some options
145  using StepWiseOptions = PropagatorOptions<StepWiseActors, Aborters>;
146  StepWiseOptions swOptions(tgContext, mfContext);
147 
148  using PlainActors = ActionList<>;
149  using PlainOptions = PropagatorOptions<PlainActors, Aborters>;
150  PlainOptions pOptions(tgContext, mfContext);
151 
152  // Run the standard propagation
153  const auto& pResult = propagator.propagate(start, pOptions).value();
154  // Let's get the end parameters and jacobian matrix
155  const auto& pJacobian = *(pResult.transportJacobian);
156 
157  // Run the stepwise propagation
158  const auto& swResult = propagator.propagate(start, swOptions).value();
159  auto swJacobianTest = swResult.template get<StepWiseResult>();
160 
161  // (1) Path length test
162  double accPath = 0.;
163  auto swPaths = swJacobianTest.paths;
164  // Sum up the step-wise paths, they are total though
165  for (auto cpath = swPaths.rbegin(); cpath != swPaths.rend(); ++cpath) {
166  if (cpath != swPaths.rend() - 1) {
167  accPath += (*cpath) - (*(cpath + 1));
168  continue;
169  }
170  accPath += (*cpath) - 0.;
171  }
172  CHECK_CLOSE_REL(swJacobianTest.fullPath, accPath, 1e-6);
173 
174  // (2) Jacobian test
175  Jacobian accJacobian = Jacobian::Identity();
176  // The stepwise jacobians
177  auto swJacobians = swJacobianTest.jacobians;
178  // The last-step jacobian, needed for the full jacobian transport
179  const auto& swlJacobian = *(swResult.transportJacobian);
180 
181  // Build up the step-wise jacobians
182  for (auto& j : swJacobians) {
183  accJacobian = j * accJacobian;
184  }
185  accJacobian = swlJacobian * accJacobian;
186  CHECK_CLOSE_OR_SMALL(pJacobian, accJacobian, 1e-6, 1e-9);
187 }
188 
189 } // namespace Test
190 } // namespace Acts