ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EigenStepper.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EigenStepper.hpp
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 #pragma once
10 
11 #include <cmath>
12 #include <functional>
13 #include <limits>
14 
24 #include "Acts/Utilities/Units.hpp"
25 
26 namespace Acts {
27 
28 using namespace Acts::UnitLiterals;
29 
42 template <typename bfield_t,
43  typename extensionlist_t = StepperExtensionList<DefaultExtension>,
44  typename auctioneer_t = detail::VoidAuctioneer>
45 class EigenStepper {
46  public:
50  using BoundState = std::tuple<BoundParameters, Jacobian, double>;
51  using CurvilinearState = std::tuple<CurvilinearParameters, Jacobian, double>;
52  using BField = bfield_t;
53 
58  struct State {
60  State() = delete;
61 
72  template <typename parameters_t>
73  explicit State(std::reference_wrapper<const GeometryContext> gctx,
74  std::reference_wrapper<const MagneticFieldContext> mctx,
75  const parameters_t& par, NavigationDirection ndir = forward,
76  double ssize = std::numeric_limits<double>::max(),
77  double stolerance = s_onSurfaceTolerance)
78  : pos(par.position()),
79  dir(par.momentum().normalized()),
80  p(par.momentum().norm()),
81  q(par.charge()),
82  t(par.time()),
83  navDir(ndir),
84  stepSize(ndir * std::abs(ssize)),
85  tolerance(stolerance),
86  fieldCache(mctx),
87  geoContext(gctx) {
88  // Init the jacobian matrix if needed
89  if (par.covariance()) {
90  // Get the reference surface for navigation
91  const auto& surface = par.referenceSurface();
92  // set the covariance transport flag to true and copy
93  covTransport = true;
94  cov = BoundSymMatrix(*par.covariance());
95  surface.initJacobianToGlobal(gctx, jacToGlobal, pos, dir,
96  par.parameters());
97  }
98  }
99 
101  Vector3D pos = Vector3D(0., 0., 0.);
102 
104  Vector3D dir = Vector3D(1., 0., 0.);
105 
107  double p = 0.;
108 
110  double q = 1.;
111 
113  double t = 0.;
114 
117 
119  Jacobian jacobian = Jacobian::Identity();
120 
122  BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
123 
125  FreeMatrix jacTransport = FreeMatrix::Identity();
126 
128  FreeVector derivative = FreeVector::Zero();
129 
132  bool covTransport = false;
133  Covariance cov = Covariance::Zero();
134 
136  double pathAccumulated = 0.;
137 
140 
142  double previousStepSize = 0.;
143 
145  double tolerance = s_onSurfaceTolerance;
146 
150  typename BField::Cache fieldCache;
151 
153  std::reference_wrapper<const GeometryContext> geoContext;
154 
156  extensionlist_t extension;
157 
159  auctioneer_t auctioneer;
160 
162  struct {
164  Vector3D B_first, B_middle, B_last;
166  Vector3D k1, k2, k3, k4;
168  std::array<double, 4> kQoP;
169  } stepData;
170  };
171 
174 
181  Vector3D getField(State& state, const Vector3D& pos) const {
182  // get the field from the cell
183  return m_bField.getField(pos, state.fieldCache);
184  }
185 
189  Vector3D position(const State& state) const { return state.pos; }
190 
194  Vector3D direction(const State& state) const { return state.dir; }
195 
199  double momentum(const State& state) const { return state.p; }
200 
204  double charge(const State& state) const { return state.q; }
205 
209  double time(const State& state) const { return state.t; }
210 
219  Intersection::Status updateSurfaceStatus(State& state, const Surface& surface,
220  const BoundaryCheck& bcheck) const {
221  return detail::updateSingleSurfaceStatus<EigenStepper>(*this, state,
222  surface, bcheck);
223  }
224 
235  template <typename object_intersection_t>
236  void updateStepSize(State& state, const object_intersection_t& oIntersection,
237  bool release = true) const {
238  detail::updateSingleStepSize<EigenStepper>(state, oIntersection, release);
239  }
240 
246  void setStepSize(State& state, double stepSize,
248  state.previousStepSize = state.stepSize;
249  state.stepSize.update(stepSize, stype, true);
250  }
251 
255  void releaseStepSize(State& state) const {
256  state.stepSize.release(ConstrainedStep::actor);
257  }
258 
262  std::string outputStepSize(const State& state) const {
263  return state.stepSize.toString();
264  }
265 
269  double overstepLimit(const State& /*state*/) const {
270  // A dynamic overstep limit could sit here
271  return -m_overstepLimit;
272  }
273 
288  BoundState boundState(State& state, const Surface& surface) const;
289 
301  CurvilinearState curvilinearState(State& state) const;
302 
307  void update(State& state, const BoundParameters& pars) const;
308 
315  void update(State& state, const Vector3D& uposition,
316  const Vector3D& udirection, double up, double time) const;
317 
325  void covarianceTransport(State& state) const;
326 
336  void covarianceTransport(State& state, const Surface& surface) const;
337 
348  template <typename propagator_state_t>
349  Result<double> step(propagator_state_t& state) const;
350 
351  private:
354 
356  double m_overstepLimit = 100_um;
357 };
358 } // namespace Acts
359