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 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 fouth 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}) \]

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.

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

class CylMapField¶

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

The field vector is stored as a cartesian \((x,y,z)\) value 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.

Field data input and options¶

JSON input for the field setup corresponds to the uniform field input celeritas::inp::UniformField and the rz-map field input:

struct RZMapFieldInput

Input data for an magnetic R-Z vector field stored on an R-Z grid.

The magnetic field is discretized at nodes on an R-Z grid, and each point the field vector is approximated by a 2-D vector in R-Z. The input units of this field are in NATIVE UNITS (cm/gauss when CGS). An optional _units field in the input can specify whether the input is in SI or CGS units, with allowable values of “si”, “cgs”, or “clhep”. The native CLHEP unit strength is 1000*tesla.

The field values are all indexed with R having stride 1: [Z][R]

Todo:: Use C indexing instead of Fortran? Or rename to ZR field?

Public Functions

inline explicit operator bool() const: Whether all data are assigned and valid.

Public Members

double min_z = {}: Lower z coordinate [len].

double max_z = {}: Last z coordinate [len].

double min_r = {}: Lower r coordinate [len].

double max_r = {}: Last r coordinate [len].

std::vector<double> field_z: Flattened Z field component [bfield].

std::vector<double> field_r: Flattened R field component [bfield].

as well as fully cartesian or cylindrical input:

struct CartMapFieldInput

Input data for a magnetic X-Y-Z vector field stored on an X-Y-Z grid.

The magnetic field is discretized at nodes on an X-Y-Z grid, and at each point the field vector is approximated by a 3-D vector in X-Y-Z. The input units of this field are in NATIVE UNITS (cm/gauss when CGS).

The field values are all indexed with Z having stride 3, for the 3-dimensional vector at that position, Y having stride (num_grid_z * 3), and X having stride (num_grid_y * num_grid_z * 3): [X][Y][Z][3]

Public Functions

inline explicit operator bool() const: Whether all data are assigned and valid.

Public Members

real_type min_x: Minimum X grid point [len].

real_type max_x: Maximum X grid point [len].

size_type num_x: Number of X grid points.

real_type min_y: Minimum Y grid point [len].

real_type max_y: Maximum Y grid point [len].

size_type num_y: Number of Y grid points.

real_type min_z: Minimum Z grid point [len].

real_type max_z: Maximum Z grid point [len].

size_type num_z: Number of Z grid points.

std::vector<real_type> field: Flattened X-Y-Z field component [bfield].

struct CylMapFieldInput

Input data for a magnetic R-Phi-Z vector field stored on an R-Phi-Z grid.

The magnetic field is discretized at nodes on an R-Phi-Z grid, and at each point the field vector is approximated by a 3-D vector in R-Phi-Z. The input units of this field are in NATIVE UNITS (cm/gauss when CGS). An optional _units field in the input can specify whether the input is in SI or CGS units, with allowable values of “si”, “cgs”, or “clhep”. The native CLHEP unit strength is 1000*tesla.

The field values are all indexed with Z having stride 1, Phi having stride (num_grid_z), and R having stride (num_grid_phi * num_grid_z): [R][Phi][Z]

Public Functions

inline explicit operator bool() const: Whether all data are assigned and valid.

Public Members

std::vector<real_type> grid_r: R grid points [len].

std::vector<RealTurn> grid_phi: Phi grid points [AU].

std::vector<real_type> grid_z: Z grid points [len].

std::vector<real_type> field: Flattened R-Phi-Z field component [bfield].

The field driver options are not yet a stable part of the API:

struct FieldDriverOptions

Configuration options for field propagation and substepping.

TODO: replace epsilon_rel_max with 1/epsilon_rel_max^2 TODO: replace safety with step_shrink_mul (or something to indicate that it’s a multiplicative factor for reducing the step, not anything with geometry) TODO: remove errcon TODO: for some of these we could probably use single-precision

Public Functions

inline explicit operator bool() const: Whether all data are assigned and valid.

Public Members

real_type minimum_step = 1.0e-5 * units::millimeter: The minimum length of the field step.

real_type delta_chord = 0.25 * units::millimeter: The maximum sagitta of each substep (“miss distance”)

real_type delta_intersection = 1.0e-4 * units::millimeter: Accuracy of intersection of the boundary crossing.

real_type epsilon_step = 1.0e-5: Discretization error tolerance for each field substep.

real_type epsilon_rel_max = 1.0e-3: Targeted discretization error for “integrate step”.

real_type errcon = 1.0e-4: UNUSED: Targeted discretization error for “one good step”.

real_type pgrow = -0.20: Exponent to increase a step size.

real_type pshrink = -0.25: Exponent to decrease a step size.

real_type safety = 0.9: Scale factor for the predicted step size.

real_type max_stepping_increase = 5: Largest allowable relative increase a step size.

real_type max_stepping_decrease = 0.1: Smallest allowable relative decrease in step size.

short int max_nsteps = 100: Maximum number of integrations (or trials)

short int max_substeps = 10: Maximum number of substeps in the field propagator.

Public Static Attributes

static constexpr real_type initial_step_tol = 1e-6: Initial step tolerance.

static constexpr real_type dchord_tol = 1e-5 * units::millimeter: Chord distance fudge factor.

static constexpr real_type min_chord_shrink = 0.5: Lowest allowable scaling factor when searching for a chord.