Propagation and magnetic field¶

The propagation interface is built on top of the geometry to allow both curved and straight-line movement. Field propagation is based on a composition of:

Field: Maps a point in space and time to a field vector.
Equation of motion: Calculates the path derivative of position and momentum given their current state and the templated field.
Integrator: Numerically integrates a new position/momentum state given the start, path derivative, and step length.
Substepper: Integrate a path segment that satisfies certain truncation error conditions, solving for the required segment length.
Propagator: Given a maximum physics step, advance the geometry state and momentum along the field lines, satisfying constraints (see field driver options) for the maximum geometry error.

Propagation¶

template<class GTV> class LinearPropagator¶: Propagate (move) a particle in a straight line.

template<class SubstepperT, class GTV> class FieldPropagator¶

Propagate a charged particle in a field.

For a given initial state (position, momentum), it propagates a charged particle along a curved trajectory up to an interaction length proposed by a chosen physics process for the step, possibly integrating sub-steps by an adaptive step control with a required accuracy of tracking in a field. It updates the final state (position, momentum, boundary) along with the step actually taken. If the final position is outside the current volume, it returns a geometry limited step and the state at the intersection between the curve trajectory and the first volume boundary using an iterative step control method within a tolerance error imposed on the closest distance between two positions by the field stepper and the linear projection to the volume boundary.

Note

This follows similar methods to Geant4’s G4PropagatorInField class.

template<class IntegratorT> class FieldSubstepper¶

Advance the field state by a single substep based on user tolerances.

The substep length is based on the radius of curvature for the step, ensuring that the “miss distance” (sagitta, the distance between the straight-line arc and the furthest point) is less than the delta_chord option. This target length is reduced into sub-substeps if necessary to meet a targeted relative error epsilon_rel_max based on the position and momentum update.

This iteratively reduces the given step length until the sagitta is no more than delta_chord . The sagitta is calculated as the projection of the mid-step point onto the line between the start and end-step points.

Each iteration reduces the step length by a factor of no more than min_chord_shrink , but is based on an approximate “exact” correction factor if the chord length is very small and the curve is circular. The sagitta h is related to the chord length s and radius of curvature r with the trig expression:

\[ r - h = r \cos \frac{s}{2r} \]

For small chord lengths or a large radius, we expand \( \cos \theta \sim 1 - \frac{\theta^2}{2} \), giving a radius of curvature

\[ r = \frac{s^2}{8h} \; . \]

Given a trial step (chord length) s and resulting sagitta of h, the exact step needed to give a chord length of \( \epsilon = {} \) delta_chord is

\[ s' = s \sqrt{\frac{\epsilon}{h}} \,. \]

Note

This class is based on G4ChordFinder and G4MagIntegratorDriver.

Field integration¶

template<class EquationT> class DormandPrinceIntegrator¶

Integrate the field ODEs using Dormand-Prince RK5(4)7M.

The algorithm, RK5(4)7M and the coefficients have been adapted from J. R. Dormand and P. J. Prince, “A family of embedded Runge-Kutta formulae” Journal Computational and Applied Mathematics, volume 6, no 1 (1980) and the coefficients to locate the mid point are taken from Shampine [1986] .

For a given ordinary differential equation, \(\difd{y}{x} = f(x, y)\), the fifth order solution, \( y_{n+1} \), an embedded fourth order solution, \( y^{*}_{n+1} \), and the error estimate as difference between them are as follows,

\[\begin{split} \begin{aligned} y_{n+1} &= y_n + h \sum_{n=1}^{6} b_i k_i + O(h^6) \\ y^{*}_{n+1} &= y_n + h \sum_{n=1}^{7} b*_i k_i + O(h^5) \\ y_{error} &= y_{n+1} - y^{*}_{n+1} = \sum_{n=1}^{7} (b^{*}_i - b_i) k_i \end{aligned} \end{split}\]

where \(h\) is the step to advance and \(k_i\) is the right hand side of the function at \(x_n + h c_i\), and the coefficients (The Butcher table) for Dormand-Prince RK5(4)7M are

  c_i | a_ij
  0   |
  1/5 | 1/5
  3/10| 3/40        9/40
  4/5 | 44/45      -56/15      32/9
  8/9 | 19372/6561 -25360/2187 64448/6561 -212/729
  1   | 9017/3168  -355/33     46732/5247  49/176  -5103/18656
  1   | 35/384      0          500/1113    125/192 -2187/6784    11/84
 ----------------------------------------------------------------------------
  b*_i| 35/384      0          500/1113    125/192 -2187/6784    11/84    0
  b_i | 5179/57600  0          7571/16695  393/640 -92097/339200 187/2100 1/40

The result at the mid point is computed

\[ y_{n+1/2} = y_n + (h/2) \sum_{n=1}^{7} c^{*}_i k_i \]

with the coefficients \(c^{*}\) taken from Shampine.

Note

The analogous class in Geant4 is G4DormandPrince745 , which inherits from G4MagIntegratorStepper.

template<class EquationT> class RungeKuttaIntegrator¶

Integrate the field ODEs using the 4th order classical Runge-Kutta method.

This method estimates the updated state from an initial state and evaluates the truncation error, with fourth-order accuracy based on description in Numerical Recipes in C, The Art of Scientific Computing, Sec. 16.2, Adaptive Stepsize Control for Runge-Kutta.

For a magnetic field equation f along a charged particle trajectory with state y, which includes position and momentum but may also include time and spin. For N-variables (i = 1, … N), the right hand side of the equation

\[ \difd{y_{i}}{s} = f_i (s, y_{i}) \]

and the fourth order Runge-Kutta solution for a given step size, h is

\[ y_{n+1} - y_{n} = h f(x_n, y_n) = \frac{h}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) \]

which is the average slope at four different points. The truncation error is the difference of the final states of one full step (y1) and two half steps (y2)

\[ \Delta = y_2 - y_1, y(x+h) = y_2 + \frac{\Delta}{15} + \mathrm{O}(h^{6}) \]

template<class EquationT> class ZHelixIntegrator¶

Analytically step along a helical path for a uniform Z magnetic field.

Given a uniform magnetic field along the z axis, \(B = (0, 0, B_z)\), the motion of a charged particle is described by a helix trajectory. For this algorithm, the radius of the helix, \(R = \frac{m v}{q B_z}\) and the helicity, defined as \( -\sgn(q B_z)\), are evaluated through the right hand side of the ODE equation where \( q \) is the charge of the particle.

The midpoint and endpoint states are calculated analytically.

Magnetic field types¶

class UniformField¶

A uniform field.

The values of the field should be in native units. For magnetic fields, this unit is gauss for the CGS system.

class RZMapField¶

Evaluate the value of magnetic field based on a volume-based RZ field map.

Values outside the grid are returned as zero.

class CartMapField¶: Interpolate a magnetic field vector on an x/y/z grid.

Warning

Accessing values outside the grid clamps to boundary values. This behavior differs from other field maps, where values outside the map are assumed zero.

class CylMapField¶

Interpolate a magnetic field vector on an r/phi/z grid.

The field vector is stored as cylindrical \((r,\phi,z)\) components on the cylindrical mesh grid points, and trilinear interpolation is performed within each grid cell. The value outside the grid is zero.

Currently the grid requires a full \(2\pi\) azimuthal grid.