26 template <
typename propagator_state_t,
typename stepper_t>
27 int bid(
const propagator_state_t& ,
28 const stepper_t& )
const {
46 template <
typename propagator_state_t,
typename stepper_t>
47 bool k(
const propagator_state_t& state,
const stepper_t&
stepper,
49 const int i = 0,
const double h = 0.,
52 stepper.charge(state.stepping) / stepper.momentum(state.stepping);
55 knew = qop * stepper.direction(state.stepping).cross(bField);
56 kQoP = {0., 0., 0., 0.};
59 qop * (stepper.direction(state.stepping) +
h * kprev).
cross(bField);
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 {
92 template <
typename propagator_state_t,
typename stepper_t>
93 bool finalize(propagator_state_t& state,
const stepper_t& stepper,
107 template <
typename propagator_state_t,
typename stepper_t>
109 const double h)
const {
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;
130 template <
typename propagator_state_t,
typename stepper_t>
152 auto& sd = state.stepping.stepData;
153 auto dir = stepper.direction(state.stepping);
155 stepper.charge(state.stepping) / stepper.momentum(state.stepping);
157 D = FreeMatrix::Identity();
159 double half_h = h * 0.5;
162 auto dFdT = D.block<3, 3>(0, 4);
163 auto dFdL = D.block<3, 1>(0, 7);
165 auto dGdT = D.block<3, 3>(4, 4);
166 auto dGdL = D.block<3, 1>(4, 7);
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);
185 (
dir + h * sd.k3).
cross(sd.B_last) + qop * h * dk3dL.cross(sd.B_last);
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();
195 dk2dT += half_h * dk1dT;
198 dk3dT += half_h * dk2dT;
205 dFdT += h / 6. * (dk1dT + dk2dT + dk3dT);
208 dFdL = (h *
h) / 6. * (dk1dL + dk2dL + dk3dL);
210 dGdT += h / 6. * (dk1dT + 2. * (dk2dT + dk3dT) + dk4dT);
212 dGdL = h / 6. * (dk1dL + 2. * (dk2dL + dk3dL) + dk4dL);
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)));