Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EigenStepper.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EigenStepper.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019-2022 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 
12 
13 template <typename E, typename A>
15  std::shared_ptr<const MagneticFieldProvider> bField, double overstepLimit)
16  : m_bField(std::move(bField)), m_overstepLimit(overstepLimit) {}
17 
18 template <typename E, typename A>
20  std::reference_wrapper<const GeometryContext> gctx,
21  std::reference_wrapper<const MagneticFieldContext> mctx,
22  const BoundTrackParameters& par, double ssize) const -> State {
23  return State{gctx, m_bField->makeCache(mctx), par, ssize};
24 }
25 
26 template <typename E, typename A>
28  const BoundVector& boundParams,
29  const BoundSquareMatrix& cov,
30  const Surface& surface,
31  const double stepSize) const {
32  // Update the stepping state
33  update(state,
35  boundParams),
36  boundParams, cov, surface);
37  state.stepSize = ConstrainedStep(stepSize);
38  state.pathAccumulated = 0.;
39 
40  // Reinitialize the stepping jacobian
41  state.jacToGlobal =
42  surface.boundToFreeJacobian(state.geoContext, boundParams);
43  state.jacobian = BoundMatrix::Identity();
44  state.jacTransport = FreeMatrix::Identity();
45  state.derivative = FreeVector::Zero();
46 }
47 
48 template <typename E, typename A>
50  State& state, const Surface& surface, bool transportCov,
51  const FreeToBoundCorrection& freeToBoundCorrection) const
53  return detail::boundState(
54  state.geoContext, state.cov, state.jacobian, state.jacTransport,
55  state.derivative, state.jacToGlobal, state.pars, state.particleHypothesis,
56  state.covTransport && transportCov, state.pathAccumulated, surface,
57  freeToBoundCorrection);
58 }
59 
60 template <typename E, typename A>
62  bool transportCov) const
63  -> CurvilinearState {
65  state.cov, state.jacobian, state.jacTransport, state.derivative,
66  state.jacToGlobal, state.pars, state.particleHypothesis,
67  state.covTransport && transportCov, state.pathAccumulated);
68 }
69 
70 template <typename E, typename A>
72  const FreeVector& freeParams,
73  const BoundVector& boundParams,
74  const Covariance& covariance,
75  const Surface& surface) const {
76  state.pars = freeParams;
77  state.cov = covariance;
78  state.jacToGlobal =
79  surface.boundToFreeJacobian(state.geoContext, boundParams);
80 }
81 
82 template <typename E, typename A>
84  const Vector3& udirection, double qOverP,
85  double time) const {
86  state.pars.template segment<3>(eFreePos0) = uposition;
87  state.pars.template segment<3>(eFreeDir0) = udirection;
88  state.pars[eFreeTime] = time;
89  state.pars[eFreeQOverP] = qOverP;
90 }
91 
92 template <typename E, typename A>
94  State& state) const {
96  state.jacTransport, state.derivative,
97  state.jacToGlobal, direction(state));
98 }
99 
100 template <typename E, typename A>
102  State& state, const Surface& surface,
103  const FreeToBoundCorrection& freeToBoundCorrection) const {
105  state.geoContext.get(), state.cov, state.jacobian, state.jacTransport,
106  state.derivative, state.jacToGlobal, state.pars, surface,
107  freeToBoundCorrection);
108 }
109 
110 template <typename E, typename A>
111 template <typename propagator_state_t, typename navigator_t>
113  propagator_state_t& state, const navigator_t& navigator) const {
114  using namespace UnitLiterals;
115 
116  // Runge-Kutta integrator state
117  auto& sd = state.stepping.stepData;
118  double error_estimate = 0.;
119  double h2 = 0, half_h = 0;
120 
121  auto pos = position(state.stepping);
122  auto dir = direction(state.stepping);
123 
124  // First Runge-Kutta point (at current position)
125  auto fieldRes = getField(state.stepping, pos);
126  if (!fieldRes.ok()) {
127  return fieldRes.error();
128  }
129  sd.B_first = *fieldRes;
130  if (!state.stepping.extension.validExtensionForStep(state, *this,
131  navigator) ||
132  !state.stepping.extension.k1(state, *this, navigator, sd.k1, sd.B_first,
133  sd.kQoP)) {
134  return 0.;
135  }
136 
137  // The following functor starts to perform a Runge-Kutta step of a certain
138  // size, going up to the point where it can return an estimate of the local
139  // integration error. The results are stated in the local variables above,
140  // allowing integration to continue once the error is deemed satisfactory
141  const auto tryRungeKuttaStep = [&](const double h) -> Result<bool> {
142  // helpers because bool and std::error_code are ambiguous
143  constexpr auto success = &Result<bool>::success;
144  constexpr auto failure = &Result<bool>::failure;
145 
146  // State the square and half of the step size
147  h2 = h * h;
148  half_h = h * 0.5;
149 
150  // Second Runge-Kutta point
151  const Vector3 pos1 = pos + half_h * dir + h2 * 0.125 * sd.k1;
152  auto field = getField(state.stepping, pos1);
153  if (!field.ok()) {
154  return failure(field.error());
155  }
156  sd.B_middle = *field;
157 
158  if (!state.stepping.extension.k2(state, *this, navigator, sd.k2,
159  sd.B_middle, sd.kQoP, half_h, sd.k1)) {
160  return success(false);
161  }
162 
163  // Third Runge-Kutta point
164  if (!state.stepping.extension.k3(state, *this, navigator, sd.k3,
165  sd.B_middle, sd.kQoP, half_h, sd.k2)) {
166  return success(false);
167  }
168 
169  // Last Runge-Kutta point
170  const Vector3 pos2 = pos + h * dir + h2 * 0.5 * sd.k3;
171  field = getField(state.stepping, pos2);
172  if (!field.ok()) {
173  return failure(field.error());
174  }
175  sd.B_last = *field;
176  if (!state.stepping.extension.k4(state, *this, navigator, sd.k4, sd.B_last,
177  sd.kQoP, h, sd.k3)) {
178  return success(false);
179  }
180 
181  // Compute and check the local integration error estimate
182  error_estimate =
183  h2 * ((sd.k1 - sd.k2 - sd.k3 + sd.k4).template lpNorm<1>() +
184  std::abs(sd.kQoP[0] - sd.kQoP[1] - sd.kQoP[2] + sd.kQoP[3]));
185  error_estimate = std::max(error_estimate, 1e-20);
186 
187  return success(error_estimate <= state.options.tolerance);
188  };
189 
190  const double initialH =
191  state.stepping.stepSize.value() * state.options.direction;
192  double h = initialH;
193  size_t nStepTrials = 0;
194  // Select and adjust the appropriate Runge-Kutta step size as given
195  // ATL-SOFT-PUB-2009-001
196  while (true) {
197  auto res = tryRungeKuttaStep(h);
198  if (!res.ok()) {
199  return res.error();
200  }
201  if (!!res.value()) {
202  break;
203  }
204 
205  const double stepSizeScaling =
206  std::min(std::max(0.25f, std::sqrt(std::sqrt(static_cast<float>(
207  state.options.tolerance /
208  std::abs(2. * error_estimate))))),
209  4.0f);
210  h *= stepSizeScaling;
211 
212  // If step size becomes too small the particle remains at the initial
213  // place
214  if (std::abs(h) < std::abs(state.options.stepSizeCutOff)) {
215  // Not moving due to too low momentum needs an aborter
216  return EigenStepperError::StepSizeStalled;
217  }
218 
219  // If the parameter is off track too much or given stepSize is not
220  // appropriate
221  if (nStepTrials > state.options.maxRungeKuttaStepTrials) {
222  // Too many trials, have to abort
223  return EigenStepperError::StepSizeAdjustmentFailed;
224  }
225  nStepTrials++;
226  }
227 
228  // When doing error propagation, update the associated Jacobian matrix
229  if (state.stepping.covTransport) {
230  // The step transport matrix in global coordinates
231  FreeMatrix D;
232  if (!state.stepping.extension.finalize(state, *this, navigator, h, D)) {
233  return EigenStepperError::StepInvalid;
234  }
235 
236  // for moment, only update the transport part
237  state.stepping.jacTransport = D * state.stepping.jacTransport;
238  } else {
239  if (!state.stepping.extension.finalize(state, *this, navigator, h)) {
240  return EigenStepperError::StepInvalid;
241  }
242  }
243 
244  // Update the track parameters according to the equations of motion
245  state.stepping.pars.template segment<3>(eFreePos0) +=
246  h * dir + h2 / 6. * (sd.k1 + sd.k2 + sd.k3);
247  state.stepping.pars.template segment<3>(eFreeDir0) +=
248  h / 6. * (sd.k1 + 2. * (sd.k2 + sd.k3) + sd.k4);
249  (state.stepping.pars.template segment<3>(eFreeDir0)).normalize();
250 
251  if (state.stepping.covTransport) {
252  state.stepping.derivative.template head<3>() =
253  state.stepping.pars.template segment<3>(eFreeDir0);
254  state.stepping.derivative.template segment<3>(4) = sd.k4;
255  }
256  state.stepping.pathAccumulated += h;
257  const double stepSizeScaling = std::min(
258  std::max(0.25f,
259  std::sqrt(std::sqrt(static_cast<float>(
260  state.options.tolerance / std::abs(error_estimate))))),
261  4.0f);
262  const double nextAccuracy = std::abs(h * stepSizeScaling);
263  const double previousAccuracy = std::abs(state.stepping.stepSize.accuracy());
264  const double initialStepLength = std::abs(initialH);
265  if (nextAccuracy < initialStepLength || nextAccuracy > previousAccuracy) {
266  state.stepping.stepSize.setAccuracy(nextAccuracy);
267  }
268  state.stepping.stepSize.nStepTrials = nStepTrials;
269 
270  return h;
271 }
272 
273 template <typename E, typename A>
275  state.jacobian = BoundMatrix::Identity();
276 }