 |
ACTS
Experiment-independent tracking
|
|
ACTS
Experiment-independent tracking
|
Multi-component fitter for electrons with non-Gaussian energy loss.
The Kalman Filter (KF) is the optimal estimator only as long as every involved distribution is Gaussian. For electrons this breaks down: their dominant energy loss is bremsstrahlung, whose Bethe–Heitler distribution is strongly non-Gaussian and heavily tailed. Feeding that through the single-Gaussian pointwise update of the KF (see the energy-loss note) biases the momentum estimate and wrecks its error estimate.
The Gaussian Sum Filter (GSF) [5] addresses this by modelling the track state as a weighted mixture of Gaussians instead of a single one,
\[ p(\vec x) = \sum_i^{N_c} w_i\, \mathcal N(\vec x \mid \vec\mu_i, \mathbf\Sigma_i), \qquad \sum_i^{N_c} w_i = 1 , \]
and running, in effect, one Kalman filter per component. Each component carries a weight, a bound parameter vector and a bound covariance:
On a material surface the GSF replaces the KF's deterministic ionization loss and Landau straggling on \(q/p\) with an explicit mixture model of the Bremsstrahlung loss. The Bethe–Heitler probability density of the energy retention \(z = E_f/E_i\) depends only on the traversed thickness in radiation lengths \(x/X_0\),
\[ f(z) = \frac{[-\ln z]^{c-1}}{\Gamma(c)}, \quad 0 \le z \le 1, \qquad c = \frac{x/X_0}{\ln 2} . \]
This is approximated by a 1D Gaussian mixture in \(z\), \(f(z) \approx \sum_n^{N_{bh}} \pi_n\,\mathcal N(z \mid \mu_{z,n}, \sigma_{z,n})\). For a single component the natural choice keeps the first two moments of the exact distribution [8],
\[ \mu_z = e^{-t}, \qquad \sigma_z^2 = 3^{-c} - 4^{-c}, \qquad t = x/X_0 , \]
which is what Acts::BetheHeitlerApproxSingleCmp evaluates:
A single Gaussian reflects the true, tailed distribution very poorly (see the figure below), so in practice a multi-component approximation is used. Because the mixture cannot be derived in closed form, its weights, means and variances are pre-fitted (minimising either the Kullback–Leibler divergence or the CDF distance to \(f(z)\)) and stored as polynomials in \(x/X_0\) so they can be interpolated at run time. Any approximation is accessed through the abstract interface:
mixture() writes \(N_{bh}\) one-dimensional components (weight, mean, variance in \(z\)) into a caller-provided span:
The concrete Acts::PolynomialBetheHeitlerApprox implements the polynomial form; the default parametrisation shipped in the source (Acts::makeDefaultBetheHeitlerApprox) and the reference JSON configuration used by the examples (betheHeitler_geantSim_cdf_nC6_O5.json) both use a six-component, fifth-order CDF fit, split into a low- and a high-thickness range at \(x/X_0 = 0.1\) [8]. When a surface exceeds the valid \(x/X_0\) range the fitter counts the occurrence and emits a warning.
Applying the loss convolves every track-state component with every Bethe–Heitler component, so a mixture of \(N_c\) components becomes \(N_c \cdot N_{bh}\). In the backward pass an effective energy gain is applied instead.
The GSF actor (Acts::detail::Gsf::GsfActor) drives the fit as a propagator actor. When the multi-stepper reports the state on a surface, the actor executes the following, in this order (see also the figure below):
\[ w_{k|k}^i \;\propto\; w_{k|k-1}^i \, \mathcal N\!\bigl(m_k \mid \mathbf H_k \vec x_{k|k-1}^i,\; \mathbf V_k + \mathbf H_k \mathbf\Sigma_{k|k-1}^i \mathbf H_k^{\mathsf T}\bigr), \]
i.e. components incompatible with the measurement are exponentially suppressed. On a passive surface a no-measurement update is done instead (which may flag a hole). Concretely, the weight is multiplied by \(\sqrt{1/\det R}\,\exp(-\tfrac12\chi^2)\), with the smallest \(\chi^2\) over the components factored out for numerical stability and the weights normalised afterwards:Left unchecked, the component count would grow by a factor \(N_{bh}\) per material surface. To keep it bounded, a reducer is invoked after each convolution to bring the mixture down to min(stepper.maxComponents, maxComponents) components. The reducer is a delegate, so it can be swapped out:
Two production reducers are provided (a third, …Naive, is a reference/benchmark baseline):
The pairwise distance driving that greedy merge is the symmetric KL divergence, restricted to the \(q/p\) dimension:
Mixture merging (or component merging) collapses a mixture into a single parameter vector and covariance — as opposed to the mixture reduction above, which brings a mixture down to a smaller mixture. Several steps need it: storing an intermediate state, and producing the final fitted parameters. The method is selectable via the Acts::ComponentMergeMethod enum:
Merging must respect cyclic bound coordinates. Which coordinates are cyclic depends on the surface type, encoded as compile-time angle descriptions (note that on a cylinder the local \(R\phi\) coordinate is cyclic, scaled by the radius):
The mean itself is then formed with complex-phase arithmetic — each cyclic coordinate is mapped onto the unit circle, averaged as a complex number, and converted back with std::arg — so that angles wrap correctly:
Each component must be transported individually, so the GSF runs on the Acts::MultiEigenStepperLoop rather than the single-component stepper. The navigator, however, must see a single trajectory. The stepper therefore presents a reduced representation to the navigation, configurable through the reducer type; the default is the highest-weight component (Acts::MaxWeightReducerLoop, with Acts::MaxMomentumReducerLoop as an alternative), which keeps the navigation stream close to the bulk of the mixture.
Determining when the whole multi-component state has "reached" a surface is handled by Acts::MultiStepperSurfaceReached, which by default treats the state as on-surface once its average is within tolerance. This guards against a pathology described in [8] — low-momentum components approaching a cylinder on a straight-line intersection can spiral indefinitely while always reporting reachable. A step limit that engages once the first component lands on the surface (stepLimitAfterFirstComponentOnSurface, default 50) forces the remaining stragglers to unreachable and removes them, after which the weights are renormalised:
Acts::GaussianSumFitter is constructed from a propagator, a shared Bethe–Heitler approximation and a logger, and exposes two fit overloads: one for the standard Acts::Navigator, and one taking an explicit surface sequence for use with the Acts::DirectNavigator — the latter is the re-fitting entry point used in the electron workflow above.
A fit runs a forward pass from the start parameters, then a backward pass that starts from the last measurement with its covariance inflated by reverseFilteringCovarianceScaling (default 100) and targets the reference surface. Measurement surfaces that were seen going forward but not on the way back are flagged as outliers. The multi-component state is merged (Mixture merging) into the single set of parameters that downstream algorithms expect; the full final mixture can optionally be attached to the track.
The knobs on Acts::GsfOptions trade physics performance against runtime. The values below summarise the scan in [8]; the ACTS example chain uses 12 components, KL-distance reduction, eMaxWeight merging and a weight cutoff of \(10^{-4}\).
| Option | Effect | Guidance [8] |
|---|---|---|
| maxComponents | mixture size after each reduction | runtime grows \(\approx\) quadratically; physics plateaus beyond \(\sim 12\) (library default 4, example default 12) |
| weightCutoff | discard components below this weight | \(10^{-4}\) is a good default; \(0.1\) is too aggressive (fit failures spike) |
| mixtureReducer | reduction algorithm | KL-distance clearly beats the weight cut at modest extra cost |
| componentMergeMethod | mixture → single estimate | eMaxWeight avoids the \(q/p\) bias seen with eMean |
| Bethe–Heitler approx | mixture model of the loss | 6-component CDF polynomials, split at \(x/X_0=0.1\) |
| reverseFilteringCovarianceScaling | covariance inflation for the backward pass | default 100 (not tuned for all setups) |
| disableAllMaterialHandling | switch off convolution and scattering | debugging only |
The payoff: against the KF, the 12-component GSF turns a heavily one-sided \(q/p\) residual into a near-symmetric one and shrinks its width, while the KF's \(q/p\) pull — its error estimate — is badly distorted by the non-Gaussian loss [8]. A single-component GSF (equivalent to a KF using the Bethe–Heitler mean and variance) is visibly biased, which is what motivates the mixture in the first place.
The per-component Bethe–Heitler application of Bethe–Heitler energy loss as a mixture is not exclusive to the fitter. The Acts::CombinatorialKalmanFilter can optionally run in a bremsstrahlung-recovery mode that reuses the same machinery to find electron tracks that a single-component filter would otherwise lose to a large, non-Gaussian energy loss.
The mode is selected purely by the stepper type. When the CKF is built over a multi-component stepper (Acts::MultiEigenStepperLoop) an IsMultiStepper trait is true and the filter compiles in a multi-component path via if constexpr; the plain single-component filter therefore carries no runtime cost, and the extra per-actor state is elided entirely with [[no_unique_address]]. On each material surface the track state is convoluted with the Bethe–Heitler mixture through the shared Acts::detail::Gsf::applyBetheHeitler, the mixture is reduced (through the same mixtureReducer delegate, now also part of the CKF extensions) and merged back to a single representation before the measurement update — the GSF's surface algorithm of The algorithm on a surface, embedded in the combinatorial search.
The multi-component knobs (maxComponents, weightCutoff, mergeMethod, betheHeitlerApprox) live on Acts::BremCombinatorialKalmanFilterOptions. The filter's Options alias resolves to that type only for a multi-stepper, so a single-component configuration cannot even name the multi-component parameters. In the examples, electron-hypothesis seeds are routed to a brem-enabled finder built over a MultiStepperLoop, using KL-distance reduction and the default Bethe–Heitler approximation, while all other seeds use the plain finder.
The per-surface algorithm itself lives in the internal Acts::detail::Gsf code; the snippets above are extracted directly from the corresponding headers.
1.15.0