 |
ACTS
Experiment-independent tracking
|
|
ACTS
Experiment-independent tracking
|
How the Acts::SympyStepper kernels are derived with sympy, and the ATLAS-form arrangement they use.
The Acts::SympyStepper does not contain a hand-written Runge-Kutta step. Its two inner kernels — the vacuum one and the dense one, each emitted with and without covariance transport — come out of Core/src/Propagator/generate_sympy_stepper.py at build time, which derives them symbolically with sympy and prints them as C++.
A Runge-Kutta stepper has to produce two things per step: the new track state, and the Jacobian of that state with respect to the old one. The second is simply the derivative of the first, but written by hand the two drift apart — a term added to the value path is easy to forget in the derivative path, and the resulting bug is a covariance that is quietly a few percent wrong rather than a test that fails.
Deriving the kernel symbolically removes that failure mode. The generator states the equations of motion once, applies the Runge-Kutta recursion to them, and lets sympy differentiate the expressions that are actually evaluated. Structural zeros, constant folding and — where it pays — common-subexpression elimination then happen on the symbolic form, before a single line of C++ exists.
The same machinery generates the bound/free Jacobian helpers (generate_sympy_jac.py) and the covariance transport (generate_sympy_cov.py); the printer and the symbolic helpers they share live in the codegen package at the repository root, which is also used by the detray backend.
cmake/ActsCodegen.cmake runs each generator through uv into a throw-away environment pinned by codegen/requirements.txt, so a build never depends on the ambient Python. codegen/manifest.json is the single place that records which generator produces which header; releases ship the generated files pre-built, and a build fully covered by them never invokes Python at all.
With \(\vec r\) the position, \(\vec T\) the unit direction, \(\lambda = q/p\), and \(s\) the path length, the vacuum kernel integrates
\[ \frac{d\vec r}{ds} = \vec T, \qquad \frac{d\vec T}{ds} = \lambda\, \vec T \times \vec B, \qquad \frac{dt}{ds} = \frac{E}{p} = \sqrt{1 + \frac{m^2}{p^2}}, \qquad \frac{d\lambda}{ds} = 0 . \]
This is a second-order system in \(\vec r\), which is what makes the Runge-Kutta-Nyström form [2] the natural integrator. Stages 2 and 3 share the midpoint, so the field is sampled at three points per step — at \(\vec r\), at the midpoint, and near the end point — rather than four.
The dense kernel adds continuous energy loss, so \(\lambda\) — and with it \(dt/ds\) — evolve along the step:
\[ \frac{d\lambda}{ds} = \frac{\lambda^2}{q}\,\frac{dt}{ds}\,g , \qquad \frac{d}{ds}\!\left(\frac{dt}{ds}\right) = \frac{m^2\lambda^3}{q^3}\, g , \]
with \(g\) the mean energy loss per unit path length [9]. Because \(\lambda\) now enters every stage, the dense kernel cannot use the pre-scaled arrangement below and builds the step Jacobian explicitly instead.
The vacuum kernel does not evaluate the plain Runge-Kutta slopes \(\vec k_i\). It follows the arrangement of the ATLAS RungeKuttaPropagator — which Acts::AtlasStepper transcribes directly — and carries a half-step bend vector at each of the three field samples:
\[ \vec H_i = \frac{h\lambda}{2}\,\vec B_i . \]
Every stage slope then comes out already scaled by \(h/2\), as a bare cross product, and neither \(h\) nor \(\lambda\) appears again anywhere in the recursion:
\[\begin{aligned} \texttt{kick1} &= \vec T \times \vec H_1 &&= \tfrac{h}{2}\vec k_1 \\ \texttt{dir2} &= \vec T + \texttt{kick1} &&= \vec T + \tfrac{h}{2}\vec k_1 \\ \texttt{dir\_half\_sum} &= \vec T + \texttt{dir2} &&= 2\vec T + \tfrac{h}{2}\vec k_1 \\ \texttt{dir3} &= \vec T + \texttt{dir2} \times \vec H_2 &&= \vec T + \tfrac{h}{2}\vec k_2 \\ \texttt{dir4} &= \vec T + \texttt{dir3} \times \vec H_2 &&= \vec T + \tfrac{h}{2}\vec k_3 \\ \texttt{dir\_end} &= 2\,\texttt{dir4} - \vec T &&= \vec T + h\,\vec k_3 \\ \texttt{kick4} &= \texttt{dir\_end} \times \vec H_3 &&= \tfrac{h}{2}\vec k_4 \end{aligned} \]
The two intermediate field sample points, the step result and the embedded error estimate all fall out of these:
\[\begin{aligned} \vec r_2 &= \vec r + \tfrac{h}{4}\,\texttt{dir\_half\_sum}, \qquad \vec r_3 = \vec r + h\,\texttt{dir4}, \\ \vec r' &= \vec r + \tfrac{h}{3}\left(\texttt{dir2} + \texttt{dir3} + \texttt{dir4}\right), \\ 3\,\vec T'_{\text{un}} &= \texttt{kick1} + 2\,\texttt{dir3} + \texttt{dir\_end} + \texttt{kick4}, \\ \varepsilon &= 2|h| \left\|\,\texttt{dir\_half\_sum} + \texttt{kick4} - \texttt{dir3} - \texttt{dir4}\,\right\|_1 = h^2 \left\|\vec k_1 - \vec k_2 - \vec k_3 + \vec k_4\right\|_1 . \end{aligned} \]
The kernel is handed the bound-to-free Jacobian \(M\) and updates it in place, never forming the 8×8 free-to-free step Jacobian \(D\). Each live column is pushed through the same recursion as the state, which is why the tangent block mirrors the value block line for line — as ATLAS' d2A/d3A/d4A block mirrors its A0..A6; see Naming for the two sets of names.
Rows are the eight free parameters \(x_i\) now, columns the six bound parameters \(b_j\) at the start surface, \(M_{ij} = \partial x_i / \partial b_j\):
| \(l_0\) | \(l_1\) | \(\phi\) | \(\theta\) | \(\lambda\) | \(t\) | |
|---|---|---|---|---|---|---|
| \(\vec r\) (3 rows) | hold | hold | step | step | step | · |
| \(t\) | · | · | · | · | step | 1 |
| \(\vec T\) (3 rows) | · | · | step | step | step | · |
| \(\lambda\) | · | · | · | · | dense | · |
· is a structural zero, 1 a constant one, hold an entry no step writes, step one every step writes, dense one only a dense step writes. The generator declares the sparsity in this one table and reads its index sets back off it.
\(\lambda\) enters the recursion only through the bend vector \(\vec H = (h\lambda/2)\,\vec B\), so the \(\lambda\) column is the only one with a term from the field's own \(\lambda\) dependence: one 3-vector at each of the four stages that use a bend vector. Stored plainly, each of the four carries a factor \(M_{\lambda\lambda}\), and the column has to be scaled by \(\lambda\) into the recursion and unscaled out of it.
Storing it differentiated by \(\log|\lambda|\) of the current \(\lambda\), with the \(\lambda\) row kept plain, removes both:
\[ M_{i\lambda} \;\equiv\; \frac{\partial x_i}{\partial \log|\lambda|} \;=\; \lambda \, \frac{\partial x_i / \partial \lambda_0} {\partial \lambda / \partial \lambda_0} , \qquad i < 7 . \]
This is exact because nothing but \(\lambda_0\) can change \(\lambda\), so \(M_{\lambda\lambda}\) is the whole chain rule from \(\partial/\partial\lambda_0\) to \(\partial/\partial\lambda\): the division applies it, the factor carries on to the log, and the plain row inverts both.
It is cheaper because \(\vec H\) is homogeneous of degree one in \(\lambda\), so \(\partial\vec H/\partial\log|\lambda| = \vec H\). Each stage's field term is then that stage's bend-linear part, with no \(M_{\lambda\lambda}\) factor and no scaling around the recursion: nineteen multiplications and a division fewer per step, 395 floating-point operations with covariance transport instead of 415, against ATLAS' 394.
That identity is the one place the recursion departs from the plain chain rule, so the generator forms the chain-rule product as well and checks the two agree (Derivation.check_same).
\(\lambda\) and \(M_{\lambda\lambda}\) are constant across a vacuum step, so it carries the scaled form into itself; the dense kernel moves \(\lambda\) and converts explicitly. The stepper state holds the scaled form, and detail::sympy::toScaledBoundToFree and its inverse convert where the covariance engine wants the plain Jacobian. The convention is singular at \(\lambda = 0\), where the plain column already is.
ATLAS' pVector[40] block is the \(M_{\lambda\lambda} = 1\) case: without dense material the row stays one, and the block stays permanently scaled by \(\lambda\).
The kernels were originally transcribed with ATLAS' variable names, which are positional rather than descriptive and are not defined in any published note — they are Athena source convention. The generator now uses names that say what the quantity is. The correspondence, for anyone reading the two side by side:
| generated name | ATLAS name | quantity |
|---|---|---|
| pos, dir, qop | R, A, P[7] (CM) | position, unit direction, \(q/p\) |
| half_h_qop | PS2 ( \(=\) Pi \(\cdot h\), with Pi \(=\lambda/2\)) | \(h\lambda/2\) |
| bend1, bend2, bend3 | H0, H1, H2 | \((h\lambda/2)\vec B_i\) at the three field samples |
| kick1 | A0, B0, C0 | \(\tfrac{h}{2}\vec k_1\) |
| dir2 | A2, B2, C2 | \(\vec T + \tfrac{h}{2}\vec k_1\) |
| dir_half_sum | A1, B1, C1 | \(2\vec T + \tfrac{h}{2}\vec k_1\) |
| dir3 | A3, B3, C3 | \(\vec T + \tfrac{h}{2}\vec k_2\) |
| dir4 | A4, B4, C4 | \(\vec T + \tfrac{h}{2}\vec k_3\) |
| dir_end | A5, B5, C5 | \(\vec T + h\vec k_3\) |
| kick4 | A6, B6, C6 | \(\tfrac{h}{2}\vec k_4\) |
| new_dir_x3 | — | \(3\times\) the unnormalised new direction |
| h_third, h_quarter, two_over_h | S3, S4, Sl | \(h/3\), \(h/4\), \(2/h\) |
| dphi_*, dtheta_*, dqop_* | d2A*, d3A*, d4A* | the tangent of the correspondingly named stage, one set per live column |
| M | pVector[8..55] | bound-to-free Jacobian, column major |
| dEds, dEds1..4 | — | energy loss per unit path, per stage (Acts::AtlasStepper has no material) |
Two things the rename fixes rather than preserves. ATLAS numbers its bend vectors H0..H2 against field samples B1..B3; here both are numbered 1..3. And the A-family index is not a stage index — A1 is used after A2 — whereas kick1/dir2/dir3/dir4/dir_end/kick4 say which stage each belongs to.
The Runge-Kutta-Nyström track model and the semi-analytic transport of the derivatives alongside the trajectory come from [2]. Adaptive Runge-Kutta-Nyström step control and error estimation in ATLAS are studied in [10], and the extrapolation package the RungeKuttaPropagator belongs to is described in [12].
None of these notes define the variable names the Athena implementation uses — they are source convention, and the table above is the only written record of what they mean.
1.15.0