EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DenseEnvironmentExtension.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DenseEnvironmentExtension.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 
11 // Workaround for building on clang+libstdc++
13 
19 
20 #include <functional>
21 
22 namespace Acts {
23 
31  double currentMomentum = 0.;
33  double initialMomentum = 0.;
38  std::array<double, 4> dLdl;
40  std::array<double, 4> qop;
42  std::array<double, 4> dPds;
44  double dgdqopValue = 0.;
46  double g = 0.;
48  std::array<double, 4> tKi;
50  std::array<double, 4> Lambdappi;
52  std::array<double, 4> energy;
53 
55  DenseEnvironmentExtension() = default;
56 
63  template <typename propagator_state_t, typename stepper_t>
64  int bid(const propagator_state_t& state, const stepper_t& stepper) const {
65  // Check for valid particle properties
66  if (stepper.charge(state.stepping) == 0. || state.options.mass == 0. ||
67  stepper.momentum(state.stepping) < state.options.momentumCutOff) {
68  return 0;
69  }
70 
71  // Check existence of a volume with material
72  if (!state.navigation.currentVolume ||
73  !state.navigation.currentVolume->volumeMaterial()) {
74  return 0;
75  }
76  return 2;
77  }
78 
92  template <typename propagator_state_t, typename stepper_t>
93  bool k(const propagator_state_t& state, const stepper_t& stepper,
94  Vector3D& knew, const Vector3D& bField, std::array<double, 4>& kQoP,
95  const int i = 0, const double h = 0.,
96  const Vector3D& kprev = Vector3D()) {
97  // i = 0 is used for setup and evaluation of k
98  if (i == 0) {
99  // Set up container for energy loss
100  auto volumeMaterial = state.navigation.currentVolume->volumeMaterial();
101  Vector3D position = stepper.position(state.stepping);
102  material = (volumeMaterial->material(position));
103  initialMomentum = stepper.momentum(state.stepping);
105  qop[0] = stepper.charge(state.stepping) / initialMomentum;
106  initializeEnergyLoss(state);
107  // Evaluate k
108  knew = qop[0] * stepper.direction(state.stepping).cross(bField);
109  // Evaluate k for the time propagation
110  Lambdappi[0] =
111  -qop[0] * qop[0] * qop[0] * g * energy[0] /
112  (stepper.charge(state.stepping) * stepper.charge(state.stepping));
113  //~ tKi[0] = std::hypot(1, state.options.mass / initialMomentum);
114  tKi[0] = std::hypot(1, state.options.mass * qop[0]);
115  kQoP[0] = Lambdappi[0];
116  } else {
117  // Update parameters and check for momentum condition
118  updateEnergyLoss(state.options.mass, h, state.stepping, stepper, i);
119  if (currentMomentum < state.options.momentumCutOff) {
120  return false;
121  }
122  // Evaluate k
123  knew = qop[i] *
124  (stepper.direction(state.stepping) + h * kprev).cross(bField);
125  // Evaluate k_i for the time propagation
126  double qopNew = qop[0] + h * Lambdappi[i - 1];
127  Lambdappi[i] =
128  -qopNew * qopNew * qopNew * g * energy[i] /
129  (stepper.charge(state.stepping) * stepper.charge(state.stepping) *
131  tKi[i] = std::hypot(1, state.options.mass * qopNew);
132  kQoP[i] = Lambdappi[i];
133  }
134  return true;
135  }
136 
147  template <typename propagator_state_t, typename stepper_t>
148  bool finalize(propagator_state_t& state, const stepper_t& stepper,
149  const double h) const {
150  // Evaluate the new momentum
151  double newMomentum =
152  stepper.momentum(state.stepping) +
153  (h / 6.) * (dPds[0] + 2. * (dPds[1] + dPds[2]) + dPds[3]);
154 
155  // Break propagation if momentum becomes below cut-off
156  if (newMomentum < state.options.momentumCutOff) {
157  return false;
158  }
159 
160  // Add derivative dlambda/ds = Lambda''
161  state.stepping.derivative(7) =
162  -std::sqrt(state.options.mass * state.options.mass +
163  newMomentum * newMomentum) *
164  g / (newMomentum * newMomentum * newMomentum);
165 
166  // Update momentum
167  state.stepping.p = newMomentum;
168  // Add derivative dt/ds = 1/(beta * c) = sqrt(m^2 * p^{-2} + c^{-2})
169  state.stepping.derivative(3) =
170  std::hypot(1, state.options.mass / newMomentum);
171  // Update time
172  state.stepping.t += (h / 6.) * (tKi[0] + 2. * (tKi[1] + tKi[2]) + tKi[3]);
173 
174  return true;
175  }
176 
190  template <typename propagator_state_t, typename stepper_t>
191  bool finalize(propagator_state_t& state, const stepper_t& stepper,
192  const double h, FreeMatrix& D) const {
193  return finalize(state, stepper, h) && transportMatrix(state, stepper, h, D);
194  }
195 
196  private:
205  template <typename propagator_state_t, typename stepper_t>
206  bool transportMatrix(propagator_state_t& state, const stepper_t& stepper,
207  const double h, FreeMatrix& D) const {
227 
228  auto& sd = state.stepping.stepData;
229  auto dir = stepper.direction(state.stepping);
230 
231  D = FreeMatrix::Identity();
232  const double half_h = h * 0.5;
233 
234  // This sets the reference to the sub matrices
235  // dFdx is already initialised as (3x3) zero
236  auto dFdT = D.block<3, 3>(0, 4);
237  auto dFdL = D.block<3, 1>(0, 7);
238  // dGdx is already initialised as (3x3) identity
239  auto dGdT = D.block<3, 3>(4, 4);
240  auto dGdL = D.block<3, 1>(4, 7);
241 
246 
247  Vector3D dk1dL = Vector3D::Zero();
248  Vector3D dk2dL = Vector3D::Zero();
249  Vector3D dk3dL = Vector3D::Zero();
250  Vector3D dk4dL = Vector3D::Zero();
251 
253  std::array<double, 4> jdL;
254 
255  // Evaluation of the rightmost column without the last term.
256  jdL[0] = dLdl[0];
257  dk1dL = dir.cross(sd.B_first);
258 
259  jdL[1] = dLdl[1] * (1. + half_h * jdL[0]);
260  dk2dL = (1. + half_h * jdL[0]) * (dir + half_h * sd.k1).cross(sd.B_middle) +
261  qop[1] * half_h * dk1dL.cross(sd.B_middle);
262 
263  jdL[2] = dLdl[2] * (1. + half_h * jdL[1]);
264  dk3dL = (1. + half_h * jdL[1]) * (dir + half_h * sd.k2).cross(sd.B_middle) +
265  qop[2] * half_h * dk2dL.cross(sd.B_middle);
266 
267  jdL[3] = dLdl[3] * (1. + h * jdL[2]);
268  dk4dL = (1. + h * jdL[2]) * (dir + h * sd.k3).cross(sd.B_last) +
269  qop[3] * h * dk3dL.cross(sd.B_last);
270 
271  dk1dT(0, 1) = sd.B_first.z();
272  dk1dT(0, 2) = -sd.B_first.y();
273  dk1dT(1, 0) = -sd.B_first.z();
274  dk1dT(1, 2) = sd.B_first.x();
275  dk1dT(2, 0) = sd.B_first.y();
276  dk1dT(2, 1) = -sd.B_first.x();
277  dk1dT *= qop[0];
278 
279  dk2dT += half_h * dk1dT;
280  dk2dT = qop[1] * VectorHelpers::cross(dk2dT, sd.B_middle);
281 
282  dk3dT += half_h * dk2dT;
283  dk3dT = qop[2] * VectorHelpers::cross(dk3dT, sd.B_middle);
284 
285  dk4dT += h * dk3dT;
286  dk4dT = qop[3] * VectorHelpers::cross(dk4dT, sd.B_last);
287 
288  dFdT.setIdentity();
289  dFdT += h / 6. * (dk1dT + dk2dT + dk3dT);
290  dFdT *= h;
291 
292  dFdL = h * h / 6. * (dk1dL + dk2dL + dk3dL);
293 
294  dGdT += h / 6. * (dk1dT + 2. * (dk2dT + dk3dT) + dk4dT);
295 
296  dGdL = h / 6. * (dk1dL + 2. * (dk2dL + dk3dL) + dk4dL);
297 
298  // Evaluation of the dLambda''/dlambda term
299  D(7, 7) += (h / 6.) * (jdL[0] + 2. * (jdL[1] + jdL[2]) + jdL[3]);
300 
301  // The following comment lines refer to the application of the time being
302  // treated as a position. Since t and qop are treated independently for now,
303  // this just serves as entry point for building their relation
304  //~ double dtpp1dl = -state.options.mass * state.options.mass * qop[0] *
305  //~ qop[0] *
306  //~ (3. * g + qop[0] * dgdqop(energy[0], state.options.mass,
307  //~ state.options.absPdgCode,
308  //~ state.options.meanEnergyLoss));
309 
310  double dtp1dl = qop[0] * state.options.mass * state.options.mass /
311  std::hypot(1, qop[0] * state.options.mass);
312  double qopNew = qop[0] + half_h * Lambdappi[0];
313 
314  //~ double dtpp2dl = -state.options.mass * state.options.mass * qopNew *
315  //~ qopNew *
316  //~ (3. * g * (1. + half_h * jdL[0]) +
317  //~ qopNew * dgdqop(energy[1], state.options.mass,
318  //~ state.options.absPdgCode,
319  //~ state.options.meanEnergyLoss));
320 
321  double dtp2dl = qopNew * state.options.mass * state.options.mass /
322  std::hypot(1, qopNew * state.options.mass);
323  qopNew = qop[0] + half_h * Lambdappi[1];
324 
325  //~ double dtpp3dl = -state.options.mass * state.options.mass * qopNew *
326  //~ qopNew *
327  //~ (3. * g * (1. + half_h * jdL[1]) +
328  //~ qopNew * dgdqop(energy[2], state.options.mass,
329  //~ state.options.absPdgCode,
330  //~ state.options.meanEnergyLoss));
331 
332  double dtp3dl = qopNew * state.options.mass * state.options.mass /
333  std::hypot(1, qopNew * state.options.mass);
334  qopNew = qop[0] + half_h * Lambdappi[2];
335  double dtp4dl = qopNew * state.options.mass * state.options.mass /
336  std::hypot(1, qopNew * state.options.mass);
337 
338  //~ D(3, 7) = h * state.options.mass * state.options.mass * qop[0] /
339  //~ std::hypot(1., state.options.mass * qop[0])
340  //~ + h * h / 6. * (dtpp1dl + dtpp2dl + dtpp3dl);
341 
342  D(3, 7) = (h / 6.) * (dtp1dl + 2. * (dtp2dl + dtp3dl) + dtp4dl);
343  return true;
344  }
345 
351  template <typename propagator_state_t>
352  void initializeEnergyLoss(const propagator_state_t& state) {
353  energy[0] = std::hypot(initialMomentum, state.options.mass);
354  // use unit length as thickness to compute the energy loss per unit length
355  Acts::MaterialSlab slab(material, 1);
356  // Use the same energy loss throughout the step.
357  if (state.options.meanEnergyLoss) {
358  g = -computeEnergyLossMean(slab, state.options.absPdgCode,
359  state.options.mass, qop[0]);
360  } else {
361  // TODO using the unit path length is not quite right since the most
362  // probably energy loss is not independent from the path length.
363  g = -computeEnergyLossMode(slab, state.options.absPdgCode,
364  state.options.mass, qop[0]);
365  }
366  // Change of the momentum per path length
367  // dPds = dPdE * dEds
368  dPds[0] = g * energy[0] / initialMomentum;
369  if (state.stepping.covTransport) {
370  // Calculate the change of the energy loss per path length and
371  // inverse momentum
372  if (state.options.includeGgradient) {
373  if (state.options.meanEnergyLoss) {
375  slab, state.options.absPdgCode, state.options.mass, qop[0]);
376  } else {
377  // TODO path length dependence; see above
379  slab, state.options.absPdgCode, state.options.mass, qop[0]);
380  }
381  }
382  // Calculate term for later error propagation
383  dLdl[0] = (-qop[0] * qop[0] * g * energy[0] *
384  (3. - (initialMomentum * initialMomentum) /
385  (energy[0] * energy[0])) -
386  qop[0] * qop[0] * qop[0] * energy[0] * dgdqopValue);
387  }
388  }
389 
398  template <typename stepper_state_t, typename stepper_t>
399  void updateEnergyLoss(const double mass, const double h,
400  const stepper_state_t& state, const stepper_t& stepper,
401  const int i) {
402  // Update parameters related to a changed momentum
403  currentMomentum = initialMomentum + h * dPds[i - 1];
404  energy[i] = std::sqrt(currentMomentum * currentMomentum + mass * mass);
405  dPds[i] = g * energy[i] / currentMomentum;
406  qop[i] = stepper.charge(state) / currentMomentum;
407  // Calculate term for later error propagation
408  if (state.covTransport) {
409  dLdl[i] = (-qop[i] * qop[i] * g * energy[i] *
410  (3. - (currentMomentum * currentMomentum) /
411  (energy[i] * energy[i])) -
412  qop[i] * qop[i] * qop[i] * energy[i] * dgdqopValue);
413  }
414  }
415 };
416 
417 template <typename action_list_t = ActionList<>,
418  typename aborter_list_t = AbortList<>>
420  : public PropagatorOptions<action_list_t, aborter_list_t> {
424  dspo) = default;
425 
431  std::reference_wrapper<const GeometryContext> gctx,
432  std::reference_wrapper<const MagneticFieldContext> mctx,
433  LoggerWrapper logger_)
434  : PropagatorOptions<action_list_t, aborter_list_t>(gctx, mctx, logger_) {}
435 
437  bool meanEnergyLoss = true;
438 
440  bool includeGgradient = true;
441 
443  double momentumCutOff = 0.;
444 
450  template <typename extended_aborter_list_t>
452  extended_aborter_list_t aborters) const {
454  eoptions(this->geoContext, this->magFieldContext, this->logger);
455  // Copy the options over
456  eoptions.direction = this->direction;
457  eoptions.absPdgCode = this->absPdgCode;
458  eoptions.mass = this->mass;
459  eoptions.maxSteps = this->maxSteps;
460  eoptions.maxStepSize = this->maxStepSize;
461  eoptions.targetTolerance = this->targetTolerance;
462  eoptions.pathLimit = this->pathLimit;
463  eoptions.loopProtection = this->loopProtection;
464  eoptions.loopFraction = this->loopFraction;
465 
466  // Stepper options
467  eoptions.tolerance = this->tolerance;
468  eoptions.stepSizeCutOff = this->stepSizeCutOff;
469  // Action / abort list
470  eoptions.actionList = this->actionList;
471  eoptions.abortList = std::move(aborters);
472  // Copy dense environment specific parameters
473  eoptions.meanEnergyLoss = meanEnergyLoss;
475  eoptions.momentumCutOff = momentumCutOff;
476  // And return the options
477  return eoptions;
478  }
479 };
480 
481 } // namespace Acts