ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DefaultExtension.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DefaultExtension.hpp
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 #pragma once
10 
12 
13 namespace Acts {
14 
19  DefaultExtension() = default;
20 
26  template <typename propagator_state_t, typename stepper_t>
27  int bid(const propagator_state_t& /*unused*/,
28  const stepper_t& /*unused*/) const {
29  return 1;
30  }
31 
46  template <typename propagator_state_t, typename stepper_t>
47  bool k(const propagator_state_t& state, const stepper_t& stepper,
48  Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP,
49  const int i = 0, const double h = 0.,
50  const Vector3D& kprev = Vector3D()) {
51  auto qop =
52  stepper.charge(state.stepping) / stepper.momentum(state.stepping);
53  // First step does not rely on previous data
54  if (i == 0) {
55  knew = qop * stepper.direction(state.stepping).cross(bField);
56  kQoP = {0., 0., 0., 0.};
57  } else {
58  knew =
59  qop * (stepper.direction(state.stepping) + h * kprev).cross(bField);
60  }
61  return true;
62  }
63 
74  template <typename propagator_state_t, typename stepper_t>
75  bool finalize(propagator_state_t& state, const stepper_t& stepper,
76  const double h) const {
77  propagateTime(state, stepper, h);
78  return true;
79  }
80 
92  template <typename propagator_state_t, typename stepper_t>
93  bool finalize(propagator_state_t& state, const stepper_t& stepper,
94  const double h, FreeMatrix& D) const {
95  propagateTime(state, stepper, h);
96  return transportMatrix(state, stepper, h, D);
97  }
98 
99  private:
107  template <typename propagator_state_t, typename stepper_t>
108  void propagateTime(propagator_state_t& state, const stepper_t& stepper,
109  const double h) const {
113  auto derivative =
114  std::hypot(1, state.options.mass / stepper.momentum(state.stepping));
115  state.stepping.t += h * derivative;
116  if (state.stepping.covTransport) {
117  state.stepping.derivative(3) = derivative;
118  }
119  }
120 
130  template <typename propagator_state_t, typename stepper_t>
131  bool transportMatrix(propagator_state_t& state, const stepper_t& stepper,
132  const double h, FreeMatrix& D) const {
151 
152  auto& sd = state.stepping.stepData;
153  auto dir = stepper.direction(state.stepping);
154  auto qop =
155  stepper.charge(state.stepping) / stepper.momentum(state.stepping);
156 
157  D = FreeMatrix::Identity();
158 
159  double half_h = h * 0.5;
160  // This sets the reference to the sub matrices
161  // dFdx is already initialised as (3x3) idendity
162  auto dFdT = D.block<3, 3>(0, 4);
163  auto dFdL = D.block<3, 1>(0, 7);
164  // dGdx is already initialised as (3x3) zero
165  auto dGdT = D.block<3, 3>(4, 4);
166  auto dGdL = D.block<3, 1>(4, 7);
167 
172 
177 
178  // For the case without energy loss
179  dk1dL = dir.cross(sd.B_first);
180  dk2dL = (dir + half_h * sd.k1).cross(sd.B_middle) +
181  qop * half_h * dk1dL.cross(sd.B_middle);
182  dk3dL = (dir + half_h * sd.k2).cross(sd.B_middle) +
183  qop * half_h * dk2dL.cross(sd.B_middle);
184  dk4dL =
185  (dir + h * sd.k3).cross(sd.B_last) + qop * h * dk3dL.cross(sd.B_last);
186 
187  dk1dT(0, 1) = sd.B_first.z();
188  dk1dT(0, 2) = -sd.B_first.y();
189  dk1dT(1, 0) = -sd.B_first.z();
190  dk1dT(1, 2) = sd.B_first.x();
191  dk1dT(2, 0) = sd.B_first.y();
192  dk1dT(2, 1) = -sd.B_first.x();
193  dk1dT *= qop;
194 
195  dk2dT += half_h * dk1dT;
196  dk2dT = qop * VectorHelpers::cross(dk2dT, sd.B_middle);
197 
198  dk3dT += half_h * dk2dT;
199  dk3dT = qop * VectorHelpers::cross(dk3dT, sd.B_middle);
200 
201  dk4dT += h * dk3dT;
202  dk4dT = qop * VectorHelpers::cross(dk4dT, sd.B_last);
203 
204  dFdT.setIdentity();
205  dFdT += h / 6. * (dk1dT + dk2dT + dk3dT);
206  dFdT *= h;
207 
208  dFdL = (h * h) / 6. * (dk1dL + dk2dL + dk3dL);
209 
210  dGdT += h / 6. * (dk1dT + 2. * (dk2dT + dk3dT) + dk4dT);
211 
212  dGdL = h / 6. * (dk1dL + 2. * (dk2dL + dk3dL) + dk4dL);
213 
214  D(3, 7) =
215  h * state.options.mass * state.options.mass *
216  stepper.charge(state.stepping) /
217  (stepper.momentum(state.stepping) *
218  std::hypot(1., state.options.mass / stepper.momentum(state.stepping)));
219  return true;
220  }
221 };
222 } // namespace Acts