Aligning tracking detectors with a global chi-squared fit
A detector can measure a hit precisely and still reconstruct a biased track if its geometry is wrong. Alignment estimates sensor positions from many tracks, while allowing each track's own parameters to change. This notebook develops the global fit and demonstrates why reference constraints are essential.
Prerequisites: linear least squares, covariance matrices and NumPy.
Install projects/requirements.txt and run all cells in order.
The worked example is synthetic: six planar sensors and straight tracks in two
dimensions. It illustrates the mathematics, not ATLAS performance.
From ATLAS geometry to a small model
ATLAS aligns its Inner Detector by fitting track-hit residuals at several levels, from large structures to individual sensors. Geometry changes with time, and track-parameter biases require independent checks beyond small residuals. The ATLAS Run-2 alignment paper describes the full system. The ATLAS explanation of detector alignment provides an accessible description of the residual-based approach and movements during data-taking.
Our example keeps only one transverse shift $\alpha_j$ for each plane at known longitudinal position $z_j$. Track $t$ has intercept $a_t$ and slope $b_t$: $$m_{tj}=a_t+b_tz_j+\alpha_j+\varepsilon_{tj},\qquad \varepsilon_{tj}\sim\mathcal N(0,\sigma^2).$$ The sign convention is explicit: a positive plane offset adds to the recorded coordinate, so the corrected hit is $m_{tj}-\widehat\alpha_j$.
Show or hide code
import numpy as np
import matplotlib.pyplot as plt
rng = np.random.default_rng(2026)
z = np.linspace(0., 5., 6) # metres
sigma = 0.015 # metres: illustrative, not a detector specification
true_offsets = np.array([0., 0.04, -0.03, 0.02, 0.01, 0.])
n_tracks = 120
true_tracks = np.column_stack([rng.normal(0, 0.3, n_tracks), rng.normal(0, 0.08, n_tracks)])
X = np.column_stack([np.ones_like(z), z])
hits = true_tracks @ X.T + true_offsets + rng.normal(0, sigma, (n_tracks, len(z)))
fig, ax = plt.subplots(figsize=(8, 4))
for j, zj in enumerate(z):
ax.axvline(zj, color='grey', alpha=0.3)
for t in range(4):
ax.plot(z, true_tracks[t] @ X.T, '--', alpha=0.6)
ax.scatter(z, hits[t], s=20)
ax.set(xlabel='Longitudinal z [m]', ylabel='Recorded transverse coordinate [m]', title='Four tracks crossing six sensor planes')
plt.tight_layout(); plt.show()
Why tracks and geometry must be fitted together
With residual $r_{tj}=m_{tj}-a_t-b_tz_j-\alpha_j$, $$\chi^2=\sum_{t,j}r_{tj}^2/\sigma^2.$$ A change in a plane position changes the best-fit parameters of every crossing track. Fitting tracks once and treating their parameters as exact loses this correlation. Linearising a general detector model gives block normal equations $$\begin{pmatrix}A&B\\B^T&C\end{pmatrix} \begin{pmatrix}\Delta\alpha\\\Delta p\end{pmatrix} =\begin{pmatrix}g_\alpha\\g_p\end{pmatrix}.$$ Here $\alpha$ denotes shared geometry and $p$ all local track parameters. Eliminating the local parameters yields the Schur complement $$(A-BC^{-1}B^T)\Delta\alpha=g_\alpha-BC^{-1}g_p.$$ The straight-line toy is already linear, so a single solve suffices. Nonlinear geometry generally requires relinearisation, refitting and convergence checks.
Gauge freedom and weak modes
If every offset is changed by $c+dz_j$, then changing every track to $a_t-c,b_t-d$ leaves each predicted hit unchanged. These are two exact null modes: a global translation and a tilt. More tracks cannot resolve them. We choose a coordinate frame by fixing the first and last plane offsets to zero. The simulated truth uses the same convention. Real alignment can instead use survey information or other external constraints, with their uncertainties. Weak modes are near-degeneracies rather than exact null modes; a small residual alone does not establish an unbiased geometry.
Show or hide code
c, d = 0.07, -0.01
shifted_offsets = true_offsets + c + d*z
shifted_tracks = true_tracks - np.array([c, d])
np.testing.assert_allclose(true_tracks @ X.T + true_offsets,
shifted_tracks @ X.T + shifted_offsets, atol=1e-15)
# Eliminating each track leaves a projector orthogonal to [1, z].
projector = np.eye(len(z)) - X @ np.linalg.solve(X.T @ X, X.T)
np.testing.assert_allclose(projector @ X, np.zeros_like(X), atol=1e-14)
print('Unconstrained geometry rank:', np.linalg.matrix_rank(projector), 'of', len(z))
assert np.linalg.matrix_rank(projector) == len(z)-2
Solve the constrained alignment
Let $D$ insert the four free shifts between the two reference planes. With equal independent measurement variances, eliminating the tracks gives $$S=\frac{N}{\sigma^2}D^TPD,\quad g=\frac{1}{\sigma^2}D^TP\sum_t m_t,\quad P=I-X(X^TX)^{-1}X^T.$$ We solve $S\widehat\alpha=g$. The free-parameter covariance is $S^{-1}$ under this linear Gaussian model with known measurement variance and exact references. It excludes any uncertainty in the assumed reference-plane positions.
Show or hide code
D = np.eye(len(z))[:, 1:-1]
S = n_tracks * D.T @ projector @ D / sigma**2
g = D.T @ projector @ hits.sum(axis=0) / sigma**2
free_offsets = np.linalg.solve(S, g)
free_covariance = np.linalg.solve(S, np.eye(len(z)-2))
fitted_offsets = D @ free_offsets
fitted_tracks = np.linalg.lstsq(X, (hits-fitted_offsets).T, rcond=None)[0].T
residuals = hits - fitted_tracks @ X.T - fitted_offsets
chi_squared = np.sum((residuals/sigma)**2)
dof = hits.size - (2*n_tracks + len(z)-2)
print('Plane truth [mm] estimate [mm] uncertainty [mm]')
for j in range(len(z)):
error = 0. if j in (0, len(z)-1) else np.sqrt(free_covariance[j-1, j-1])
print(f'{j:3d} {1000*true_offsets[j]:10.3f} {1000*fitted_offsets[j]:14.3f} {1000*error:16.3f}')
print(f'chi2 / dof = {chi_squared:.1f} / {dof} = {chi_squared/dof:.3f}')
assert np.max(np.abs((free_offsets-true_offsets[1:-1])/np.sqrt(np.diag(free_covariance)))) < 4
assert 0.7 < chi_squared/dof < 1.3
Independent check: one joint least-squares fit
We now build the full design matrix for all track and geometry parameters. This is less efficient for a large detector, but it independently checks the Schur-complement signs, parameter ordering and marginal geometry covariance.
Show or hide code
joint = np.zeros((hits.size, 2*n_tracks+len(z)-2))
for t in range(n_tracks):
rows = slice(t*len(z), (t+1)*len(z))
joint[rows, 2*t:2*t+2] = X
joint[rows, 2*n_tracks:] = D
joint_fit = np.linalg.lstsq(joint/sigma, hits.ravel()/sigma, rcond=None)[0]
joint_cov = np.linalg.solve(joint.T @ joint/sigma**2, np.eye(joint.shape[1]))
np.testing.assert_allclose(joint_fit[2*n_tracks:], free_offsets, atol=1e-12)
np.testing.assert_allclose(joint_fit[:2*n_tracks].reshape(-1, 2), fitted_tracks, atol=1e-12)
np.testing.assert_allclose(joint_cov[2*n_tracks:, 2*n_tracks:], free_covariance, atol=1e-12)
print('Joint fit and eliminated-track fit agree, including geometry covariance.')
Show or hide code
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(z, true_offsets*1000, 'o--', label='True shifts')
errors = np.r_[0., np.sqrt(np.diag(free_covariance)), 0.]
ax.errorbar(z, fitted_offsets*1000, yerr=errors*1000, fmt='s', capsize=4, label='Fit ±1 sigma')
ax.set(xlabel='z [m]', ylabel='Plane shift [mm]', title='Geometry in the chosen reference frame')
ax.legend(); plt.tight_layout(); plt.show()
Validate on independent tracks
New tracks were not used to determine the geometry. Fit their intercepts and slopes before and after correcting the hits, and compare the residual bias at each plane. Track fitting removes two components of every six-hit residual vector, so its spread is not directly the single-hit resolution $\sigma$. A real experiment also checks physics observables, different track populations, time stability and distortions that may leave residuals nearly unchanged.
Show or hide code
validation_tracks = np.column_stack([rng.normal(0, 0.3, 200), rng.normal(0, 0.08, 200)])
validation_hits = validation_tracks @ X.T + true_offsets + rng.normal(0, sigma, (200, len(z)))
before = validation_hits @ projector
corrected = validation_hits-fitted_offsets
after = corrected @ projector
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(z, before.mean(axis=0)*1000, 'o-', label='Before')
axes[0].plot(z, after.mean(axis=0)*1000, 's-', label='After')
axes[0].axhline(0, color='grey', lw=1)
axes[0].set(xlabel='z [m]', ylabel='Mean residual [mm]', title='Held-out track residuals')
axes[1].hist(before.ravel()*1000, bins=40, alpha=0.5, label='Before')
axes[1].hist(after.ravel()*1000, bins=40, alpha=0.5, label='After')
axes[1].set(xlabel='Residual [mm]', ylabel='Hits', title='Held-out residual distribution')
for ax in axes: ax.legend()
plt.tight_layout(); plt.show()
assert np.linalg.norm(after.mean(axis=0)) < np.linalg.norm(before.mean(axis=0))
Exercise: can more tracks determine the absolute origin?
Try increasing the number of tracks while freeing all six plane shifts. The two null modes remain. Fixing references is a coordinate choice; setting an arbitrary regularisation penalty instead silently introduces information. The following answer scaffold exhibits the null directions directly.
Show or hide code
np.testing.assert_allclose(projector @ np.ones(len(z)), 0., atol=1e-14)
np.testing.assert_allclose(projector @ z, 0., atol=1e-14)
print('Translation and tilt remain unconstrained, regardless of sample size.')
What the toy leaves out
An actual tracking detector needs curved track propagation in a magnetic field, multiple scattering and correlated measurement errors, rotations and deformations, nonuniform track acceptance, outlier handling and geometry changes with time. A hierarchy can first fit large assemblies, then smaller structures. External constraints and independent bias measurements are essential where tracks alone cannot determine a deformation. This toy uses exact straight-line propagation, no scattering and a single time-independent geometry.
Interactive experiment: recover detector shifts
Vary the number of tracks, measurement noise and injected shift, then compare the recovered offsets with their uncertainties. The reference constraints stay in place: increasing the sample size cannot resolve the translation and tilt that were unconstrained in the derivation.
The controls use a separate, seeded synthetic example and run in your browser
in the exported HTML article. Some notebook viewers do not execute embedded
JavaScript; run scripts/check_notebooks.py --html and open the article preview
to use the controls.
Show or hide code
from pathlib import Path
import sys
interactive_project = next(p / "projects" for p in [Path.cwd(), *Path.cwd().parents]
if (p / "projects" / "src").is_dir())
sys.path.insert(0, str(interactive_project / "src"))
from notebook_interactive import display_interactive
display_interactive("alignment")
Conclusion
Alignment requires fitting shared geometry while allowing each track's own parameters to change. In this six-plane example, eliminating the local track parameters gives the same offsets and marginal geometry covariance as the full joint fit, to the checked numerical tolerance. The recovered free-plane shifts are within about one reported standard deviation of the simulated values, with uncertainties of roughly $1.7$–$1.8$ mm. The fitted $\chi^2/\mathrm{DOF}=1.03$ is consistent with the noise scale used to generate this sample.
Those uncertainties apply within the chosen reference frame. The unconstrained geometry has rank four rather than six: a translation and a tilt can be absorbed by the track intercepts and slopes. More tracks cannot determine these two exact null modes. Fixing the first and last planes makes the solution identifiable; the quoted covariance treats those references as exact and excludes uncertainty in their positions.
Correcting a separate set of tracks reduces the mean residual bias, providing a check beyond the alignment sample. Small residuals still cannot establish an unbiased absolute geometry, because some distortions can be absorbed by fitted tracks. The useful workflow is to identify unconstrained modes, state the reference information explicitly, and validate on independent data and observables. This straight-line toy demonstrates that workflow; it does not establish the accuracy achievable in a full detector.
References
ATLAS Collaboration, Alignment of the ATLAS Inner Detector in Run-2, Eur. Phys. J. C 80 (2020) 1194, doi:10.1140/epjc/s10052-020-08700-6. The context above follows that paper; all numerical examples and figures here are generated by the explicitly stated synthetic model.