Kalman Track Fitting

Before any of the mathematics, the shape of the problem. A collision leaves the detector full of unlabelled hits — clusters of charge on silicon, each one telling you where something crossed, and nothing about what. Turning that into a list of measured particles happens in three distinct stages, and it is worth being precise about what each one produces, because in casual use all three get called “tracking”.

• Seeding. Find small groups of hits — almost always three — consistent with a single track coming from near the beamline. A seed is a hypothesis, together with a first, crude estimate of the track parameters and their covariance. Its only job is to give the next stage somewhere to start.
• Track finding. Carry that hypothesis outward and decide, layer by layer, which hits belong to it. The output is a set of hits — still not a measurement. This is section 6.
• Track fitting. Given the hits, extract the best estimate of the track parameters and, more importantly, the covariance the physics will actually consume. This is sections 1–5 and 7.

The figure runs all three on one simulated event in a nine-layer barrel tracker, seen down the beam axis. Every number in it is computed live, from the same equations the rest of the page derives.

Three things in that figure are worth carrying into the rest of the page.

The seed is bad, and that is fine. Three pixel hits over an 11 cm lever arm pin the momentum to a few percent at 1 GeV, and to about 15% at 5 GeV — the stiffer the track, the less a short lever arm can say about its curvature. The seed is not trying to measure anything. It is trying to be right about which hits, and to be wrong about the parameters by a knowable amount, because that amount is what sets the search road.

The road is the seed's error, extrapolated. It opens from ±0.2 mm at the fourth layer to ±24 mm at the ninth, and it pinches back every time a hit is absorbed. That fanning shape is not decoration — it is ±3√S, computed from the covariance, and section 6 is about nothing else.

Fitting is where the measurement happens. Finding hands over a list of hits; fitting turns that list into a number with an error bar, and it is worth roughly a factor of fifty in momentum resolution over the seed that started it. Everything from section 1 to section 7 is about doing that step properly — and about section 8's question, which is whether the error bar it hands back is telling the truth.

One last thing the figure hides. The finding stage and the fitting stage are the same filter. Finding runs it forward with branching, scoring candidate hits as it goes; fitting runs it over the hits that survived, then smooths backward. Once you have section 4's five equations, you have most of both.

The three-stage framing, and the figure it is drawn from, follow the ACTS documentation's Tracking in a nutshell.

You measure the same quantity twice, independently: m₁ with uncertainty σ₁, and m₂ with σ₂. The maximum-likelihood combination is the inverse-variance weighted mean, and its variance follows:

Variances do not add. Information — inverse variance — adds. That is the whole idea, and the rest is bookkeeping. Now rewrite the same equation as an update to the first measurement rather than a symmetric average of the two:

Algebraically nothing happened. Conceptually everything did. You now have a current best estimate; a new measurement arrives; you move toward it by a fraction K of the disagreement. That fraction is the Kalman gain, and it is nothing but a statement of relative trust:

• \(K\to 1\) when you knew almost nothing before (\(\sigma_1\gg\sigma_2\)). Take the new measurement at face value.
• \(K\to 0\) when the new measurement is much worse than what you have. Barely move.
• \(\hat\sigma^{2}=(1-K)\,\sigma_1^{2}<\sigma_1^{2}\) always. Absorbing a measurement can never make you less certain.

Drag the sliders below until the gain stops being a formula and starts being obvious.

A charged particle in a magnetic field follows a helix, which takes five parameters (a sixth, time, comes along if your detector measures it). What varies between experiments is only which five, and what surface they are referred to. The state vector x is those parameters; the covariance matrix C is their 5×5 error matrix. A track, to a fitter, is exactly this pair — nothing else.

Every one of those real parameterisations carries the momentum as q/p, the signed inverse momentum. This is not a convention, it is load-bearing, for three reasons:

1. It is what the detector measures. Curvature is proportional to q/p, so the measurement model is close to linear in this parameter and badly nonlinear in p.
2. Its errors are close to Gaussian. A symmetric Gaussian in q/p is a long asymmetric tail in p. The filter assumes Gaussians; give it the variable where that assumption is least of a lie.
3. Infinite momentum is an ordinary point. A very stiff track has q/p → 0, which is finite and perfectly well behaved. Parameterise in p and every straight track blows up your covariance matrix.

Strip away everything not needed to see the mechanism: work in the bending plane only, with planes perpendicular to z, and approximate the helix by its parabola (valid while the bend angle stays well under a radian). Three parameters remain:

The third parameter c is the curvature, and it is the momentum measurement in disguise. For a track of transverse momentum p_T in a field B,

So fitting c and quoting σ_c is fitting the momentum and quoting its resolution — the toy is small, but it is not a cartoon.

Given the state at plane k−1, the equation of motion tells you where the particle should be at plane k. Write that map as f_k and its Jacobian as F_k = ∂f/∂x. The prediction is then

Read a|b as “at plane a, using hits up to b”. In the toy, transport over a gap Δz is exact and linear:

Q_k is the process noise — the randomness the material itself adds between one plane and the next. Dominantly that is multiple Coulomb scattering, whose projected angle has RMS given by the Highland formula,

plus energy-loss straggling in the q/p element. Notice what Q does structurally: it lets the true state itself change between measurements. A least-squares fit forces one rigid helix through every hit; a Kalman filter fits a trajectory that is allowed to kink at each scatterer by an amount the material budget says is reasonable. That is why it is the natural fitter for a detector made of stuff.

In the toy, a thin scatterer kicks the slope and leaves the position alone — hence the single non-zero element of Q in equation (6). It does not touch the curvature: the field sets that, not the material.

A measurement m_k arrives with error matrix V_k. The projection matrix H_k says which combination of state parameters the sensor actually sees — one row for a silicon strip, two for a pixel. Then:

Set H = 1 and all the matrices to scalars and equations (10)–(12) collapse exactly onto equation (2). The gain is still “my variance divided by the total variance”; the update is still “move by K times the disagreement.”

The line that explains why anyone bothers is (10). The only matrix ever inverted is S, whose size is the dimension of a single measurement — 1×1 or 2×2. There is no large matrix inversion anywhere, and the cost of the whole fit is linear in the number of hits. A global least-squares fit of the same track with scattering handled properly requires inverting an N×N correlation matrix, because scattering at layer 3 correlates the residuals at layers 4 through N.

Below is the toy of section 2, fitted with the equations of section 4 — the same code, verified line by line against a NumPy reference. Step the filter forward with the arrows and watch three things.

• The gain starts at exactly 1 and decays. With three parameters and fewer than three hits, the filter has no opinion to defend: it takes each hit at face value.
• χ² stays at zero until the fourth hit. A three-parameter curve passes exactly through any three points, so the first three residuals are not tests of anything. This is exactly why \(\mathrm{ndf}=N-3\) — and why a real five-parameter track fit has \(\mathrm{ndf}=N-5\).
• The filtered estimate is always best at the plane you just left. Watch the curvature error in the lower panel: it falls monotonically and only reaches its final value at the last plane. Everywhere else, the forward filter is worse than it needs to be. Fixing that is the entire job of section 7.

You can also drag any hit vertically. Pull one badly off the track and watch what a single outlier does to the fit — and to every residual downstream of it.

Every hit so far arrived pre-assigned to your track. Nothing in a real event does that. A layer in a busy event carries thousands of clusters and no labels, and the question “which of these belongs to the track I am following?” has to be answered before the filter can absorb anything.

The answer is already in section 4. At layer k you hold a prediction and its covariance, and equation (9) gives S_k — the covariance of the residual before you know which hit will produce it. So any candidate hit m can be scored with the same expression the fit already uses,

The region of the surface satisfying (14) is the search window — also called the road, or the validation gate. In one dimension it is the interval \(|r| < n\sqrt{S}\) with \(n=\sqrt{\chi^{2}_{\max}}\). Everything else in this section is consequences of that one line.

Two things about S that decide whether your window is right

It is built from S, not from the prediction error. S = HCH^T + V carries both the uncertainty in where the track will be and the uncertainty in where the sensor puts the hit. Leaving V out is a tempting shortcut — the window is “about the extrapolation”, after all — and it gets worse precisely as the fit gets better, because HCH^T shrinks while V does not:

In more than one dimension it is an ellipse, not a box. On a pixel sensor S is 2×2 and its off-diagonal term is generally not zero — a track crossing a module at a shallow angle has correlated errors in the two local coordinates. A rectangle is still useful, but only as a pre-filter: it is cheap, it can be evaluated against a layer's binned hit index without touching a covariance matrix, and it cuts the list down before the real test runs. It has to circumscribe the ellipse, never crop it.

The window narrows as the fit tightens

Take the three innermost planes as given — that is what a seeding algorithm produces — and search outward from plane 3. The window shrinks at every step, and its composition changes as it does:

At the first search layer the window is 95% your ignorance about the track. By the last it is 71%, and heading toward the floor set by the sensor itself: however good the fit gets, the window can never shrink below V. A tracker whose windows have bottomed out at the hit resolution has extracted everything the geometry has to give.

The strategy, in the order production code runs it

1. Navigate. Ask the geometry which surfaces the trajectory crosses next, and in what order. This is a geometry problem, not a statistics one, and in a real detector it is most of the running time.
2. Propagate to the surface, transporting C with F and adding Q — equation (5).
3. Bound. Form S, take a rectangle that circumscribes the χ² ellipse, and query the layer's binned hit index for everything inside it.
4. Gate. Score that shortlist with equation (14) and keep what passes.
5. Branch. Every surviving hit spawns its own copy of the track state, updated with that hit. Spawn a hole branch as well — the candidate continuing with no hit on this layer — so that a dead module or a genuine inefficiency does not silently kill a real track.
6. Prune, or drown. Cap the branches per seed, cut on accumulated χ², cut on the number of holes, and drop candidates that share too many hits with a better one.
7. Repeat until the layers run out, keep the surviving candidate or candidates, and only then smooth — section 7.

Steps 1–4 and 7 are the fit you already know. Steps 5 and 6 are the whole difference between a fitter and a combinatorial Kalman filter, and they are where the CPU goes.

Choosing χ²_max

The gate is a straight trade. The true hit lands inside a symmetric \(n\sqrt{S}\) window with probability \(\mathrm{erf}(n/\sqrt{2})\), while the number of unrelated hits admitted grows linearly with n — and each one is a branch, which is a track candidate, which costs a full propagation at the next layer:

Efficiency saturates; the combinatorics do not. That asymmetry is why n = 3 is such a common default — it buys 99.73% of the true hits, and everything past it costs branches to purchase a fraction of a percent. Note also that this compounds: 99.73% per layer across five search layers is 98.7% of tracks keeping every hit, and a gate at n = 2 would leave you 79.2%.

Why one branch is not enough

That 99.73% is the probability the true hit lies inside the window given a correct prediction. It is not the probability that a search which always follows its single best candidate ends up on the right track, because the two are coupled: take the wrong hit once and the filtered state moves with it, so the next window is centred in the wrong place, so the next pick is likelier to be wrong again. Failures compound. Here is the greedy strategy — always absorb the smallest χ², never branch — measured over five search layers at n = 3:

On a clean layer the greedy search is essentially perfect. Add four tenths of an unrelated hit per centimetre — a quiet event by any collider standard — and one search in six already ends on the wrong hit or on nothing at all. None of that is the window's fault: the window is correct at every step, and the true hit is inside it 99.73% of the time it is given a fair chance. It is the cost of committing to one answer at a layer where the data do not yet distinguish between two. Carrying both and deciding later, once more layers have voted, is exactly what branching buys — and why no production tracker searches greedily.

The figure below runs the greedy search, so you can watch it fail: it follows the best χ² at every layer and marks the places where a real CKF would have forked instead.

Four things are worth doing to it. Turn the gate down to n = 1 and the true hit starts falling outside its own window: the track loses hits it genuinely produced, and the fit continues on whatever happened to stay inside. Turn the gate up to 6 and nothing is ever missed — but count the branches. Turn the unrelated-hit density up and the window does not change at all, yet the candidate count climbs as a product over layers; that product is what pruning exists to contain. Push the scattering slider and watch √S grow: at 15 mrad the window at layer 5 goes from 1183 to 2189 µm, and a fitter that had left Q out would still be searching the narrow one. Underneath all of it, watch the window narrow layer by layer — the search gets cheaper and sharper the further along the track you already are, which is why seeding well matters so much.

Where Q comes back

Since S = H(FCF^T + Q)H^T + V, material widens the window as well as the errors. Whether that matters at all is a single comparison — the scattering displacement across the gap, θ₀Δz, against √S:

So in a fine-grained silicon tracker with 10 cm gaps, dropping Q from the window costs a stiff track essentially nothing — 20 µm against a 1183 µm window is not a measurable effect, and you can convince yourself the term is optional. It stops being optional exactly where tracking is hard: soft tracks, thick material, and the long extrapolations across a gap between subdetectors, where θ₀Δz can exceed the extrapolation error outright. In figure 3 the window at plane 5 goes 1183 → 1247 → 1729 → 2189 µm as the scattering slider moves 0 → 3 → 10 → 15 mrad — it nearly doubles. A gate tuned on 10 GeV muons and shipped for 1 GeV pions is running windows too narrow in the one regime that needed them wide.

The filtered state at plane k has seen hits 1…k. Only the last plane has seen everything. But the parameters an analysis wants — impact parameter, direction and momentum at the innermost layer, where the track meets the vertex — live at the plane the filter knew least about. The forward pass alone is the wrong tool for the job you actually have.

The Rauch–Tung–Striebel smoother fixes this by walking back, carrying the information that arrived later:

Read (16) carefully: the correction is driven by x_k+1|n − x_k+1|k, the amount by which later hits moved the state away from what plane k predicted. If the prediction turned out to be right, the smoother leaves plane k alone. No measurement appears anywhere in (15)–(17) — the backward pass only redistributes information the forward pass already extracted, which is why you must store both x_k|k−1, C_k|k−1 and x_k|k, C_k|k at every plane on the way out — a real memory cost, and the reason a filter that only needs the final parameters can skip the smoother entirely.

Frühwirth's original paper gives an alternative, the two-filter smoother: run a second, independent filter inward from the last plane, then combine the outward and inward estimates at each plane with the weighted mean of equation (1). It costs more, but it hands you an estimate at plane k that is statistically independent of the hits on one side — convenient when you want a residual that excludes the hit under test.

From the default toy — 8 planes over 70 cm, 500 µm hits, no scattering. Two different stories, one for each kind of parameter:

Curvature is flat after smoothing. It has to be: c is a property of the whole trajectory, so once every hit has been seen, the answer is the same wherever you evaluate it. The filtered column, by contrast, is a record of the filter slowly learning — and if you read the momentum off the innermost plane of a filter-only fit, you get the seed's error, not the tracker's.

Position is not flat. The smoothed σ_y is best in the middle (241 µm) and worst at both ends (421 µm), because the interior planes are interpolations between hits and the end planes are extrapolations beyond them. This is the same reason a straight-line fit's error band has a waist. Note also that both end values equal the filtered value at the last plane — the smoothed fit run in either direction must agree at the point where all the information has arrived.

Track parameter uncertainties propagate into everything downstream — invariant mass resolution, impact-parameter significance for b-tagging, vertex χ², isolation. A fit that is 15% optimistic about its momentum error does not announce itself; it quietly biases a whole analysis. So you have to test C directly, and there are exactly two standard instruments.

Pulls

Generate tracks whose true parameters you know, fit them, and histogram

If the model is right this is a unit Gaussian, and the two ways it fails are diagnostic: a shifted mean means the fit is biased; a width above 1 means the errors are underestimated; below 1, overestimated or double-counted.

χ² per degree of freedom

Sum equation (13) over planes. With a diffuse seed, the degrees of freedom are (total measurement dimensions) − (number of track parameters): N−5 for a real 5-parameter fit, N−3 here.

Both instruments are below. The selector deliberately breaks the fit in three different ways — the point is to learn what each failure looks like, because you will meet all of them.

Units, and the condition number

Put q/p in MeV⁻¹ instead of GeV⁻¹ and that element of C changes by 10⁶; mix it with a position in millimetres and the diagonal of your covariance matrix spans a dozen orders of magnitude. Every subsequent inversion loses digits to that spread. Pick one unit system, write it at the top of the file, and never guess.

The seed, from both sides

C₀ has to be diffuse enough that the answer does not depend on it — and small enough that equation (17) does not cancel catastrophically. Both walls are real, and the second one is much closer than people expect. Here is this page's toy, with everything held fixed except the curvature element of the seed. The correct answer, from an unconstrained least-squares fit of the same eight hits, is \(\sigma_c=7.7152\times10^{-5}\,\mathrm{cm}^{-1}\):

Read the last two rows twice. Equation (17) subtracts two nearly equal large matrices; once the seed is wide enough, the digits that survive are noise. At \(\sigma_c=2\) that noise happens to be negative and you get a loud NaN. At 5 and 10 it happens to be positive, the covariance passes every positive-definiteness check you might write, and the fit calmly reports momentum errors twenty-two times too large. Positive-definite is not the same as correct. The only real test is a pull distribution — which is why section 8 exists.

The alternative, standard in production fitters, is to seed properly from the first few hits (a three-point circle fit) and inflate that covariance by a modest factor, so the dynamic range never gets large in the first place.

Keep the covariance a covariance

Symmetrise after every update — C ← (C+C^T)/2 costs nothing and stops asymmetry compounding. Use the Joseph form. Assert that the diagonal is positive rather than discovering it via a NaN. When you need real robustness, the square-root and UD-factorised filters propagate a factor of C instead of C, halving the dynamic range that has to fit in a double.

One bad hit ruins everything after it

A wrong hit — a neighbouring track's, a noise cluster — enters the fit at full weight and drags the state. Every prediction downstream is then wrong, so their residuals are wrong too, and the damage propagates forward. The cheap defence is a cut on the per-hit χ²_k of equation (13): reject a hit whose predicted residual is implausible before it is absorbed. Try it in figure 2 by dragging a hit far off the track. The principled version is the deterministic annealing filter, which gives every hit a weight between 0 and 1 rather than a binary verdict.

Direction matters, physically

Filtering inside-out and outside-in give the same smoothed answer in exact arithmetic, but they are not interchangeable in practice: energy loss has a sign that depends on which way you are propagating, and the linearisation errors differ. Choose deliberately, and be consistent between fitting and any subsequent extrapolation.

• Material done properly. Real Q needs a geometry that knows its own radiation lengths, and energy loss as a bias in q/p rather than just a variance. This is most of the engineering in a production fitter.
• When Gaussian is wrong. Electrons radiate, and Bethe–Heitler energy loss is nothing like a Gaussian. The Gaussian-sum filter carries a weighted mixture of Gaussians instead of one, running a bank of Kalman filters in parallel. The deterministic annealing filter handles contaminated hits the same way, with soft weights annealed over iterations.
• Vertices. Vertex fitting is the identical machinery with a different H: the state is the vertex position (plus track momenta), each track is a “measurement”, and tracks are added one at a time exactly as hits are here. Frühwirth's 1987 paper introduced both applications together for precisely this reason.

For code to read rather than write: ACTS is the modern, experiment-independent reference implementation, and its track-fitting documentation maps directly onto the vocabulary above — predicted, filtered, smoothed. GenFit is a smaller, more readable toolkit if you want to follow a fit end to end in an afternoon.