Paul Avery CBX 92–39 August 20, 1992

Applied Fitting Theory V Track Fitting Using the Kalman Filter

Section 1: Introduction Section 2: How Tracks are Measured Section 3: The Standard Track Fit Section 4: Problems with the Standard Fit Section 5: Track Fitting with the Kalman Filter Section 6: Startingthe Kalman Fit

Appendix 1: Multiple Scattering Formulas Appendix 2: Energy Loss Formulas Appendix 3: Stepping Through Gases

1. Introduction

I present in this note a relatively new approach to track fitting that uses an iterative algorithm to correctly account for the effects of multiple scattering and energy loss along the track trajectory. This technique, known generally as a Kalman filter, was first applied to track fitting by P. Billoir^{1} and is now being used by several CERN experiments. Some useful discussions have recently started appearing in the literature^{2,3}. I implemented it in the FTMONTE package some time ago but it is not used in our standard data reduction yet. My intention is that this note will serve as a mathematical reference for an installation into DUET sometime later this year.

My discussion here focuses exclusively on track fitting rather than track finding, although it is not hard to figure out how the filter can be used to select hits as well as fit them. I have also left out a nice discussion of how the method can be applied to vertex finding and fitting. I will discuss that topic in a future note, but I will say in passing that I have implemented a very powerful set of vertex fitting routines in KNLIB based on the filter.

An extensive discussion of the motion of charged particles in magnetic fields can be found in CBX 92–45^{4}. Appendices 1 and 2 provide a useful set of formulas for dealing with multiple scattering and energy loss when fitting tracks with the Kalman filter.

2. How Tracks are Measured

Charged particle tracks in high energy physics experiments are measured through a two step process. First, detector measurements (e.g., drift chambers, straws, MWPCs, silicon strips, etc.) are put through a pattern finding algorithm to select a subset that seems consistent with belonging to a single track. This set of measurements is then fit statistically through a maximum likelihood method to determine the most probable set of track parameters consistent with the measurements. Errors in these parameters are estimated from the measurement uncertainties and other factors, as discussed below. Measurements may be further eliminated during the track fit if they are found to be of sufficiently poor quality or not consistent with belonging to the track. Because of its well understood properties, the least squares algorithm is most commonly used to fit the track parameters and estimate the parameter errors.

In a constant magnetic field oriented along the z axis (such as found in CLEO) a charged particle moves along a helical trajectory, a path consisting of circular motion in the x − y plane combined with a constant velocity along the z axis. Unfortunately, several detector related effects complicate the trajectory and therefore the fitting process: (1) inhomogeneous magnetic fields, (2) energy loss and (3) multiple scattering. I consider only dE/dx and multiple scattering corrections in this note.

Energy loss from particle interactions with detector walls and gases causes the particle to gradually decrease its radius of curvature, resulting in a spiralling-like effect for very soft tracks. The energy loss has a deterministic (average) component that can be calculated via the Bethe–Bloch dE/dx formula and a stochastic (fluctuating) piece caused by the large Landau tail of the dE/dx distribution. Correcting the mean energy loss is straightforward and is implemented at the DST level in CLEO, not in the DUET fit (see CBX 92-40 for a discussion^{5}). This “after the fact” correction is approximate but works pretty well since most of the energy loss occurs before the outer drift chamber where the momentum is mainly determined. It works less well for soft tracks where the momentum determination depends more equally on measurements in the inner tracking chambers and the original track fit attempted to fit the whole track to a single helix.

As a charged particle traverses a medium it is deflected by many small angle scatters, mostly due to coulomb interactions with atomic electrons. The scattering occurs in the two perpendicular planes to the direction of motion and for small angles has a gaussian distribution in each plane (uncorrelated), hence it is known as multiple coulomb scattering or just multiple scattering. For larger angles, one should also consider the hard scattering component which adds a non-gaussian tail to the multiple scattering distribution, but we won’t worry about this effect here.

3. The Standard Track Fit

To appreciate how the Kalman fitting algorithm works, let’s compare it to the method used in DUET which uses a general least squares fit. The helix is specified by the 5 parameters α = (c, ϕ_{0}, D, λ, z_{0}) where c is 1/2 the inverse curvature of the track, ϕ_{0} is the ϕ coordinate of the momentum at the r − ϕ point of closest approach to the origin, D is the signed distance of closest approach to the origin in r −ϕ, λ = cot θ and z_{0} is the z position at the distance of closest approach to the origin. The equations of motion of tracks using these parameters are described in CBX 92–45. The general least squares algorithm was discussed in detail in CBX 91–72^{6} and is summarized as follows using matrix notation. Suppose we want to fit a set of n measurements y = (y_{1}, y_{2}, … y_{n}) to a set of m parameters α (m = 5 here) through the relation y_{l} = f_{l}(α) for 1 ≤ l ≤ n. In track fitting, the meaning of y_{l} depends on the type of detector, i.e., it can be the distance of closest approach to a wire (drift chambers and straws), the z coordinate (cathode strip) or a distance from a metal strip in a plane (silicon strip detector).

Since the f_{l}(α) are nonlinear we expand them about an approximate solution α = α_{A} (corresponding to a first guess of the track parameters): y_{l} = f_{l}(α_{A}) + (∂f_{l}/∂α_{i})(α_{i} − α_{A i}). This linearization permits us to define the χ^{2} statistic

χ^{2} = ∑_{l}(y_{l} − (α_{A}) − A_{li}(α_{l} − α_{A l}))^{2}/σ_{l}^{2}
= (y − f(α_{A}) − A(α − α_{A}))^{T} V_{y}^{−1} (y − f(α_{A}) − A(α_{A} − α))
≡ (Δy − A(α − α_{A}))^{T} V_{y}^{−1} (Δy − A(α − α_{A}))

where Δy = y − f(α_{A}), A_{li} = ∂f_{l}(α)/∂α_{i}_{α_{A}} is a constant matrix of derivatives, and

V_{y}^{−1} = ( \begin{matrix}1/σ_{1}^{2} & 0 & ⋯ & 0 \\ 0 & 1/σ_{2}^{2} & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ⋯ & 1/σ_{n}^{2} \end{matrix} )

is the inverse of the covariance matrix of the measurements. It is diagonal for independent measurements. Since the χ^{2} measures how much the measurements “miss” the function, the solution we want is that which minimizes χ^{2}, i.e., ∂χ^{2}/∂α_{i} = 0. The solution was found to be α = α_{A} + V_{A}A^{T} V_{y}^{−1}Δy with covariance matrix V_{α}, where V_{α} = V_{A} = (A^{T} V_{y}^{−1}A)^{−1}. In theory one should iterate this procedure to guarantee that the final track is close to the χ^{2} minimum.

4. Problems with the Standard Fit

The standard fit simultaneously fits all the measurements and so can be regarded as a global method. When the measurement errors σ_{l} are independent of each other it executes in a time proportional to n, the number of measurements. The situation, however, changes drastically when we introduce multiple scattering because a single scattering event affects the drift distance measurements at all subsequent points. The measurement covariance matrix is now V_{y} = V_{meas} + V_{MS}, where V_{meas} is the original contribution from independent measurement errors and V_{MS} is the correlated component due to the presence of the multiple scattering event and contains off diagonal terms (see CBX 91–74^{7} for the calculation of V_{MS}). The problem now is obtaining V_{y}^{−1} when V_{y} is no longer diagonal: the inversion of the covariance matrix is prohibitively time consuming when the number of measurements n is large (execution time is O(n^{3})). In the CLEO central detector, for instance, a track can have as many as 71 measurements. Inverting a 71 × 71 matrix at least once per track to account for multiple scattering is enormously expensive in computer time and clearly unacceptable.

There is an extensive literature about the multiple scattering problem in high energy physics. DUET partially accounts for multiple scattering by including a kink in the r − ϕ plane at the entrance to the DR drift chamber. The rest of it is included by a method described in CBX 91–74^{7} in which multiple scattering is ignored in the fit but its effects on the covariance matrix is accounted for by propagation of errors. The implementation does not work well for soft tracks which bend significantly and in any case is non-optimal because the multiple scattering information is not used in the determination of the track parameters.

The standard fit also does not account for energy loss as the track propagates. Implementing dE/dx corrections cannot be done by linear expansion because soft tracks lose large amounts of energy in detector walls. In any case, the equations become very difficult to write down.

5. Track Fitting with the Kalman Filter

The Kalman filter method, in contrast to the standard algorithm, is an iterative procedure that is, in some sense, intuitively obvious (in fact, parts of the idea have been used for years in track fitting). The procedure iteratively traces the track from its outermost point to the origin, picking up measurements and accounting for multiple scattering and dE/dx losses along the way. Unlike global methods which fit all the measurements to a single set of track parameters, the filter causes the track to “follow the measurements” through the detector.

Let the n hits on the track trajectory be divided among r regions, labelled 1, 2, … r starting from the outermost one, where the regions are separated from one another by material boundaries. A particularly simple case is the one where every hit lies in a separate region because, as we shall see, no matrix inversion need be done. Assume for the time being that the parameters α_{1} and their covariance matrix V_{α 1} have somehow been determined in region 1 (this avoids the tricky subject of starting the fit, a topic that I defer till later). The track parameters in region 2 are obtained by minimizing a χ^{2} which includes contributions from the measurements in region 2 and the projected track from region 1. The projection must include corrections to the track parameters from the expected energy loss in the material (α_{1} → α′_{1}) and to the track covariance matrix to account for multiple scattering and fluctuations in the energy loss (V_{α 1} → V′_{α 1}). The formulas for these corrections can be found in Appendices 1 and 2.

The change in χ^{2} for region 2 can be written as a sum of two terms, the first accounting for the measurements in region 2 and the second allowing for the pull of the track parameters from their previously fitted values (after dE/dx corrections), e.g.:

Δχ^{2} = (Δy_{2} − A_{2}(α_{2} − α′_{1}))^{T} V_{meas 2}^{−1}(Δy_{2} − A_{2}(α_{2} − α′_{1})) + (α_{2} − α′_{1})^{T} V′_{α 1}^{−1}(α_{2} − α′_{1})

where A_{2 li} = ∂f_{2 l}(α)/∂α_{i} are the derivatives of calculated drift distance in region 2 with respect to the five parameters, α′_{1} are the track parameters from region 1 corrected for energy loss and V′_{α 1} = V_{α 1} + V_{α MS} + V_{α energy} is the covariance matrix of the track from region 1 corrected by the contributions from multiple scattering and energy fluctuations. As before, Δy_{2} = y_{2} −f(α′_{1}) represents the vector of differences between the measurements and the predictions based on the parameters before fitting. The solution, obtained by setting ∂χ^{2}/∂α_{i} to zero, is

α_{2} = α′_{1} + V_{A g}A_{2}^{T} V_{meas 2}^{−1}Δy_{2},
V_{α 2} = V_{A g},
Δχ^{2} = (Δy_{2})^{T} V_{meas 2}^{−1}Δy_{2} − (Δy_{2})^{T} V_{meas 2}^{−1}A_{2}V_{A g}A_{2}^{T} V_{meas 2}^{−1}Δy_{2},

where

V_{A g} = (A_{2}^{T} V_{meas 2}^{−1}A_{2} + V′_{α 1}^{−1})^{−1} = V′_{α 1}(1 + A_{2}^{T} V_{meas 2}^{−1}A_{2}V′_{α 1})^{−1}.

The second form of V_{A g} is preferred because it avoids the calculation of V′_{α 1}^{−1}.

This procedure is repeated for the subsequent regions, each iteration improving the track parameters and their covariance matrix, and terminates when the beginning of the track is reached (usually when it breaks through the beam pipe). The net effect of the Kalman filter is to continuously update the track parameters and their covariance matrix.

If we include the measurements a few at a time the fitting process can be speeded up by avoiding the 5 × 5 matrix inversion in each iteration of the above formulas. This follows from applying the Woodbury matrix inversion theorem^{6} to the calculation of V_{A g} in the above equations. The theorem states that if A is an invertible n × n matrix and U and V are n × p matrices (with p < n), then

(A + UV^{T})^{−1} = A^{−1} − A^{−1}U(1 + V^{T} A^{−1}U)^{−1}V^{T} A^{−1}

where (1 + V^{T} A^{−1}U) is a p × p matrix. Here n = 5 and p is the number of hits being added from a region. The formula for V_{A g} now becomes

V_{A g} = V_{α 1} − V_{α 1}A_{2}^{T} (V_{meas 2} + A_{2}V_{α 1}A_{2}^{T} )^{−1} A_{2}V_{α 1}

We must be careful here to include the multiple scattering terms in the equation for V_{meas 2} as worked out in reference 7. Of course, if a large number of hits are included at one time the inversion of the measurement covariance matrix becomes prohibitively expensive in computer time. Appendix 3 shows a faster way of adding the hits that relies on a careful choice of step size in the gas.

For the special case p = 1, i.e., hits added one at a time, the entire process becomes very fast since the matrix inversion reduces to a simple arithmetic inverse. The expressions for α_{2} and its covariance matrix V_{α 2} now become

α_{2} = α′_{1} + V_{A g}A_{2}^{T} \frac{δy_{i}}{σ_{i}^{2}} ,
V_{α 2} = V_{A g} = V_{α 1} − \frac{V_{α 1}A_{2}^{T} A_{2}V_{α 1}}{σ_{i}^{2} + A_{2}V_{α 1}A_{2}^{T}} .

These equations assume that the track has been traced as closely to the point of closest approach as possible so that further multiple scattering and energy corrections need not be done.

There are several interesting points about the Kalman filter that should be noted:

1. The Kalman filter method uses all the information and cannot, if used correctly, give poorer track parameters by adding more measurements. For example, the contribution from the back half of a curler would be small because multiple scattering errors would reduce the significance of the parameters determined in this region. Including these hits would not cause the track parameters to worsen, unlike many other track fitting methods.

2. Because the algorithm traces the track backwards, the parameters on the outer part of the track are much more poorly determined than the ones on the inner part. Although this is good for physics reasons (we want the parameters at the production point), it means that the projection of the track to the outer region is unreliable. To do the projection properly, the track must be refit by tracing it outward from the production point to the outermost part of the drift chamber. It might be necessary to store these parameters in ROAR, but probably not the covariance matrix. A study should be performed to determine how accurately we need to carry out the projection.

3. The Kalman filter is mathematically equivalent to the global fit when there is no multiple scattering or energy loss. That this is true can be seen from the first form of the solution for α_{2} above.

4. Every region does not have to contain measurements. If a particular region lacks measurements, the track is simply projected through the region to the next one.

5. The natural iterativeness of this procedure allows tracks to be fit in pieces. For example, tracks can be found and fit in the drift chambers and then projected to the silicon with the full covariance matrix. The silicon hits can be added to the fit trivially without having to do the entire fit again. Many experiements already exploit this technique.

6. The Kalman filter technique is sensitive, especially at low momentum, to the particle hypothesis used in the fit (e, μ, π, K or p). For best results, we will probably want to fit each fit 5 times and store all the results, although we will need to devise a clever coding mechanism at the ROAR level to avoid storing excessive numbers of bytes (especially for the covariance matrices). An added benefit of separate hypothesis fitting is that we can adjust the cell drift times using the correct time of flight, rather than always using the π hypothesis. It should also be possible to perform some crude particle ID by comparing the pattern of measurement residuals.

7. The Kalman filter is well suited to object oriented programming (OOP) techniques. In OOP language, a region object takes an incoming track object, extrapolates the track across its boundaries, modifies the track parameters by adding its own measurements (if any), and extrapolates the track to the next boundary, where the next region takes over.

6. Starting the Kalman Fit

I have avoided until now the issue of how to start the Kalman filter fit. The problem is that in the outer region there are only a few measurements and the track parameters determined there are very poorly determined, causing the extrapolation into the next region to be unreliable. Moreover, the track parameters in this region become extremely sensitive to the vagaries of these few hits and, even worse, there may not be enough measurements to start the fit, especially if the hits lack z information.

There seems to be little published information on how to solve this problem. However, I have come up with a technique that might be robust enough work in data. The stratagem, which is used in FTMONTE, goes as follows. The track is first fit using all the hits in the outer drift chamber. I then multiply the 5 × 5 track covariance by a scale factor in the range of 200 to 1000. The track parameters with the new covariance matrix becomes the initial track that I use to start the fit in region 1 (I pretend that the track enters region 1 with these parameters). The scale factor must be small enough so that the track has some initial “stiffness” so that it cannot easily be pulled away from its starting value, but it must be large so that the contribution of the drift chamber hits is not double counted in the final fit.

I have not carried out any systematic studies on how to start the fit in an optimal way. Currently only the diagonal part of the covariance matrix are used and I scale the 5 terms by a different amount in FTMONTE (I use 200, 1000, 1000, 500, 500 for the 5 parameters). It might also be useful to choose the outermost region to include more hits so that the track is reasonably stabilized after the first region. These issues need to be explored more fully.

Appendix 1

Multiple scattering formulas Multiple scattering causes the particle to scatter in the two planes perpendicular to its

path without losing energy. The distribution of each angle is approximately gaussian with a standard deviation given by

σ_{θ} = \frac{0.0141}{pβ}\sqrt{\frac{L}{X_{R}}} ≡ \sqrt{HL}

where p is the particle momentum in GeV/c, β is its velocity, L is the length of the path and X_{R} is the radiation length of the material. The track also changes position in the two planes with a standard deviation in each plane given by

σ_{x} = \frac{0.0141}{pβ}L\sqrt{\frac{L}{3X_{R}}} ≡ L\sqrt{\frac{HL}{3}}

The position and angle changes in each plane are correlated with a correlation coefficient ρ = \sqrt{3}/2 ≃ 0.866. The scatterings in each plane are independent of one another.

Suppose we want to trace the track to a radius r. In the absence of multiple scattering the track would end up at (x, y, z) (or (r, ϕ, θ) in spherical coordinates) and its direction would be described by the spherical angles (ϕ_{t}, θ_{t}). Multiple scattering causes the track to scatter in the plane oriented perpendicular to its path. The unit vectors defining the plane are ϕ̂_{t} = (− sin ϕ_{t}, cos ϕ_{t}, 0) and θ̂_{t} = (cos θ_{t} cos ϕ_{t}, cos θ sin ϕ_{t}, − sin θ_{t}).

The parameters will be changed by an amount

δα = M ( \begin{matrix}sin θ_{t}δϕ_{t} \\ δθ_{t} \\ δx_{ϕ_{t}}/L \\ δx_{θ_{t}}/L \end{matrix} )

where M_{i1} = ∂α_{i}/ sin θ_{t}∂ϕ_{t}, M_{i2} = ∂α_{i}/∂θ_{t}, M_{i3} = L∂α_{i}/∂x_{ϕ_{t}} and M_{i4} = L∂α_{i}/∂x_{θ_{t}}, are the derivatives of the track parameters with respect to the scattering quantities and x_{ϕ_{t}} and x_{θ_{t}} are distances along the unit vectors ϕ̂_{t} and θ̂_{t}, respectively. Note that the perpendicular plane is in general not tangent at radius r to the sphere centered at 0 because of the bending of the track.

Using the correlations, the contribution to the track covariance matrix due to multiple scattering is

V_{α MS} ≡ ⟨δαδα^{T} ⟩ = (HL)M ( \begin{matrix}1 & 0 & 1/2 & 0 \\ 0 & 1 & 0 & 1/2 \\ 1/2 & 0 & 1/3 & 0 \\ 0 & 1/2 & 0 & 1/3 \end{matrix} ) M^{T} .

The derivative matrix can be computed using the notation of CBX 92–45 with r^{2} = x^{2} + y^{2} and T = \sqrt{p_{⊥}^{2} − 2a(xp_{y} − yp_{x}) + a^{2}r^{2}} = p_{⊥} + aD = p_{⊥}(1 + 2cD). In the interests of fitting the results on a single page we write M as the direct sum of two 5 × 2 matrices, e.g., M = M_{angle} ⊕ M_{position}, where M_{angle} and M_{position} account for the angular and positional multiple scattering effects, respectively. A straightforward but tedious calculation yields

M_{angle} = ( \begin{matrix}0 & −\frac{1}{2}aλ/p_{⊥} \\ (p/p_{⊥})(p_{⊥}^{2} − axp_{y} + ayp_{x})/T^{2} & aλ(xp_{x} + yp_{y})/T^{2} \\ −(p/p_{⊥})(xp_{x} + yp_{y})/T & −(aλ/2)(r^{2} − D^{2})/T \\ 0 & −(p/p_{⊥})^{2} \\ λp(xp_{y} − yp_{x} − ar^{2})/T^{2} & s_{⊥} + λ^{2}p_{⊥}(xp_{x} + yp_{y})/T^{2} \end{matrix} )

= ( \begin{matrix}0 & −\frac{1}{2}aλ/p_{⊥} \\ pp_{⊥}(1 − ρr sin α)/T^{2} & aλp_{⊥}r cos α/T^{2} \\ −pr cos α/T & −(aλ/2)(r^{2} − D^{2})/T \\ 0 & −(p/p_{⊥})^{2} \\ −λpp_{⊥}(D + c(D^{2} + r^{2})/T^{2} & s_{⊥} + λ^{2}p_{⊥}^{2}r cos α/T^{2} \end{matrix} )

M_{position} = ( \begin{matrix}0 & 0 \\ ρp(xp_{x} + yp_{y})/T^{2} & −ρ cos θ(p_{⊥}^{2} − axp_{y} + ayp_{x})/T^{2} \\ (p/p_{⊥}^{2})(p_{⊥}^{2} − axp_{y} + ayp_{x})/T & ρ cos θ(xp_{x} + yp_{y})/T \\ 0 & 0 \\ ρλp(xp_{x} + yp_{y})/T^{2} & [−T^{2} sin θ + λ cos θ(p_{⊥}^{2} − axp_{y} + ayp_{x})]/T^{2} \end{matrix} )

= ( \begin{matrix}0 & 0 \\ apr cos α/T^{2} & −ap_{z} sin θ(1 − ρr sin α)/T^{2} \\ p(1 − ρr sin α)/T & ar cos θ cos α/T \\ 0 & 0 \\ aλpr cos α/T^{2} & sin θ[−T^{2} + p_{z}^{2}(1 − ρr sin α)]/T^{2} \end{matrix} )

where sin θ = 1/\sqrt{1 + λ^{2}} and s_{⊥} is calculated from

\begin{aligned}T sin ρs_{⊥} & = ρ(xp_{x} + yp_{y}), \\ T cos ρs_{⊥} & = p_{⊥} − ρ(xp_{y} − yp_{x}). \end{aligned}

α is the angle between the radius vector and the track direction. Recall from CBX 92–45 that

\begin{aligned}sin α & = \frac{xp_{y} − yp_{x}}{rp_{⊥}} = rc − \frac{D}{r}(1 + cD), \\ cos α & = \frac{xp_{x} + yp_{y}}{rp_{⊥}} = ϵ\sqrt{(1 − D^{2}/r^{2})[(1 + cD)^{2} − r^{2}c^{2}]}, \end{aligned}

where ϵ = +1(−1) refers to the outgoing (incoming track). It is equal to the sign of xp_{x}+yp_{y}.

Appendix 2

Energy loss formulas There is a longer discussion of dE/dx corrections in CBX 92–40.^{5} The corrected momen-

tum p′ is calculated from the initial momentum p by p′ = p + Δp, where Δp is the energy loss correction calculated by integrating the Bethe-Bloch formula

Δp = ∫_{0}^{t} \frac{dp}{dx}dx = ∫_{0}^{t} \frac{1}{β} \frac{dE}{dx} dx = ∫_{0}^{t} \frac{A}{β^{3}} [ln (\frac{2m_{e}c^{2}β^{2}γ^{2}}{I_{0}} ) − β^{2}] dx

where t is the thickness of the material to be traversed, A is a constant that depends on the material, Z is the atomic number, m_{e} is the electron mass, I_{0} is the ionization potential and β and γ are the usual relativistic parameters. We use the approximation I_{0} = Z^{0.9}eV, and take Z = 9 for representative media. A is determined from our knowledge of (dE/dx)_{min} at minimum ionization, the quantity listed in the Particle Data Book. The minimum occurs at βγ = 3.40 and gives A = (dE/dx)_{min}/11.528.

It is easier to write the differential equation in terms of τ = p/m and the normalized thickness x′ = x/t:

\begin{aligned}\frac{dτ}{dx′} & = \frac{ΔE_{min}}{m} F (τ) \\ F (τ) & = \frac{(1 + τ^{2})^{3/2}}{11.528τ^{3}} [9.0872 + 2 ln(τ) − \frac{τ^{2}}{1 + τ^{2}} ] \end{aligned}

where ΔE_{min} = (dE/dx)_{min}t is the average energy loss suffered by a minimum ionizing particle passing through a distance t of material (including path length corrections from non-normal tracks). The change in momentum can then be calculated from the integral

Δp = mΔτ = ΔE_{min} ∫_{0}^{1} F (τ)dx′

For moderate velocities (β > 0.7), F (τ) is approximately constant giving Δp ≃ ΔE_{min}F (τ). For low velocities, where F (τ) is changing rapidly, it is probably necessary to integrate the dτ/dx′ differential equation numerically using a method such as the Runge–Kutte technique, as discussed in CBX 92–40.^{5} The integration is aided by transforming to a new variable ξ = τ^{4} which removes the 1/τ^{3} term from the differential equation and makes the numerical solution far more reliable. The new equation using this variable is

\begin{aligned}\frac{dξ}{dx′} & = \frac{ΔE_{min}}{m} F′(ξ) \\ F′(ξ) & = \frac{(1 + \sqrt{ξ})^{3/2}}{4 · 11.528} [9.0872 + \frac{1}{2} ln(ξ) − \frac{\sqrt{ξ}}{1 + \sqrt{ξ}} ] \end{aligned}

Once the new momentum has been determined, the formulas in Section 1.4 in CBX 92–45 can be used to calculate the new helix parameters.

Appendix 3

Stepping Through Gases If several hits are to be fit simultaneously in a region, it is important to add the multiple

scattering corrections to the measurement covariance matrix V_{meas} using formulas found in thelast section of CBX 91–74.^{7} Since these corrections add off diagonal terms to V_{meas}, computing its inverse is time consuming if a moderately large number of measurements are added at one time. In this section I show that it is possible to ignore these corrections during the fit providing that we choose a step size small enough. After the fit is complete, one should add the multiple scattering contribution due to the material in the region to the new track covariance matrix using the formulas in Appendix 1.

The criterion I use is that the expected deviation of the track over the whole region from multiple scattering must be less than a certain amount, both in the r − ϕ and r − z planes. Let Δ_{ϕ} and Δ_{z} be the maximum average deviation that will be tolerated in r − ϕ and along z, respectively. The RMS deviation in each plane from multiple scattering over a length L in material of radiation length X_{R} is (0.0141L/pβ)\sqrt{L/3X_{R}}, where the step ΔL must satisfy

\frac{0.0141ΔL}{pβ}\sqrt{\frac{ΔL}{3X_{R}}} < Δ_{ϕ},
\frac{0.0141ΔL}{pβ sin θ}\sqrt{\frac{ΔL}{3X_{R}}} < Δ_{z},

where 1/ sin θ is a geometrical effect accounting for the error in z at fixed radius for steeply falling tracks. Rearranging, we find that ΔL must satisfy both

ΔL < Δ_{ϕ}^{max} = (\frac{3Δ_{ϕ}^{2}p^{2}β^{2}X_{R}}{0.0141^{2}} )^{1/3}
ΔL < Δ_{z}^{max} = (\frac{3Δ_{z}^{2}p^{2}β^{2}X_{R} sin^{2} θ}{0.0141^{2}})^{1/3}

Note that the step size is much more sensitive to the momentum (∼ (pβ)^{2/3}) than to the radiation length of the gas (∼ X_{R}^{1/3}). The following table shows maximum step lengths for 50–50 Argon-Ethane (X_{R} = 166 m) as a function of momentum (assuming sin θ = 0.8) with maximum allowed scattering Δ_{ϕ} = 50μm and Δ_{z} = 100μm.

Maximum size step size (in cm) in 50–50 Argon–Ethane for π, K, p

based on Δ_{ϕ} = 50μm (Δ_{z} = 100μm)

p π K p
Δ_{ϕ}^{max} (Δ_{z}^{max}) Δ_{ϕ}^{max} (Δ_{z}^{max}) Δ_{ϕ}^{max} (Δ_{z}^{max})
0.05 1.21 (1.66) 0.54 (0.74) 0.35 (0.48)
0.10 2.77 (3.79) 1.35 (1.85) 0.89 (1.22)
0.25 6.68 (9.14) 4.31 (5.89) 2.96 (4.05)
0.50 11.3 (15.5) 9.25 (12.7) 7.02 (9.61)
1.00 18.3 (25.1) 17.1 (23.5) 14.9 (20.4)
1.50 24.1 (33.0) 23.3 (31.9) 21.6 (29.6)
2.00 29.2 (40.0) 28.7 (39.2) 27.4 (37.5)
3.00 38.3 (52.4) 38.0 (52.0) 37.2 (50.8)
4.00 46.4 (63.5) 46.2 (63.2) 45.6 (62.4)
5.00 53.9 (73.5) 53.7 (73.5) 53.3 (72.9)

References

1. P. Billor, Track Fitting with Multiple Scattering: A New Method, NIM 225, 352-366 (1984).

2. M. Regler and R. Fruwirth, Reconstruction of Charged Particle Trajectories, Plenum Publishing Company (1989).

3. R. Bock, H. Grote, D. Notz and M. Regler, Data Analysis Techniques for High-Energy Physics Experiments, Cambridge University Press (1990).

4. P. Avery, Applied Fitting Theory IV: Formulas for Track Fitting, CBX 92–45.

5. P. Avery, The Charged Particle Energy Correction in CLEO, CBX 92–40.

6. P. Avery, Applied Fitting Theory I: General Least Squares Theory, CBX 91–72.

7. P. Avery, Applied Fitting Theory III: Non-Optimal Least Squares Fitting and Multiple Scattering, CBX 91–74.