Paul Avery CBX 91 − 74 October 18, 1991

Applied Fitting Theory III

Non-Optimal Least Squares Fitting and Multiple Scattering

Section 1: The Problem

Section 2: Non-Optimal Least Squares Fitting Section 3: Non-Optimal Fitting with Constraints Section 4: Simple Application to Multiple Scattering Section 5: Summary of Multiple Scattering Formulas

1. The Problem

In the first paper in this series, I gave a general overview of least squares fitting theory and showed that if the dependence of the measurements on the parameters can be linearized, the parameters which minimize the χ^{2} function can be obtained rather easily, even in the presence of constraints.

To recapitulate the argument, suppose we want to fit a set of n measurements y to a set of m parameters α through the relation y_{l} = f_{l}(α) for 1 ≤ l ≤ n. If the f_{l}(α) are nonlinear we can expand them about an approximate solution α = α_{A}: y_{l} = f_{l}(α_{A})+(∂f_{l}/∂α_{i})(α_{i} −α_{A i}). This linearization permits us to define the χ^{2} statistic as

χ^{2} = (y − f_{A}(α_{A}) − A(α − α_{A}))^{t}V_{y}^{−1} (y − f_{A}(α_{A}) − A(α_{A} − α))
≡ (Δy − Aη)^{t}V_{y}^{−1} (Δy − Aη)

where Δy = y − f_{A}(α_{A}), A_{li} = ∂f_{l}(α)/∂α_{i}_{α_{A}} is a constant matrix, η = α − α_{A} is the new vector of unknowns, and V_{y}^{−1} is the inverse of the covariance matrix of the 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 η = V_{A}A^{t}V_{y}^{−1}Δy with covariance matrix V_{η}, where V_{η} = V_{A} = (A^{t}V_{y}^{−1}A)^{−1}. The Gauss-Markov theorem states that the parameters obtained by this procedure are both unbiased and have minimum uncertainties, i.e., they are the best parameters that can be determined by any method. So far, so good.

Under normal circumstances the measurment errors σ_{l} are independent of each other, i.e., the weight matrix V_{y}^{−1} can be written in diagonal form:

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

which is obtained by inversion of the diagonal covariance matrix V_{y ij} = ⟨δy_{i}δy_{j}⟩ = σ_{i}^{2}δ_{ij}. Thus it is a trivial matter to fit an arbitrary number of measurements.

There are situations, however, where the data measurements y_{l} are not independent of each other, so that movement of a particular measurement will cause the simultaneous adjustment of other measurements. An example of this is multiple scattering, in which a single scattering event at a particular point affects the drift distance measurements at all subsequent points. In this case the deviation of the measurements from the original track will have an independent component due to random measurement errors and a correlated component due to the presence of the multiple scattering event. The correlations among the data points can be computed and put into the measurement covariance matrix V_{y}. The problem now is obtaining V_{y}^{−1}: the inversion of the covariance matrix can be prohibitively time consuming when the number of measurements n is large. 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 just not very practical, unless one has a lot of CPU time to waste.

2. Non-Optimal Least Squares Fitting

Several attempts have been made in the literature to account for multiple scattering effects in track fitting. The methods employed range from the mundane to the elegant, but a common goal links all of them: finding a set of track parameters whose errors are as close as possible to the optimal values calculated in Section 1, but without the high cost in CPU time which that technique entails. I call this the “optimal” approach.

For most tracks inclusion of multiple scattering effects in the fitting algorithm does not improve the parameters of the fit significantly, except for very slow tracks. A large amount of programming and CPU time is invested in making marginal improvements in the track parameters. In fact, what is many times desired is not the best possible parameters, but parameters whose errors are well understood. Knowledge of these errors is important when they are used as input to kinematic fitting or when a sensitive lifetime measurement depends critically on uncertainties in the vertex resolution. In both cases, lack of understanding of the parameter errors and their correlations can give erroneous physics results, especially when those results are based on a quoted number of standard deviations.

I start by noting that the procedure whereby one looks for parameters which minimize the χ^{2} function is merely one possible estimation scheme — although it is the optimal one — and that others are possible. Assume, for instance, that an approximation for V_{y}^{−1} exists which we call Ṽ_{y}^{−1}. Ṽ_{y}^{−1} will almost certainly be diagonal, but this assumption is not necessary in what is to follow. We define a modified form of the χ^{2} function

χ̃^{2} = (y − f_{A} − Aη)^{t}Ṽ_{y}^{−1}(y − f_{A} − Aη)
≡ (Δy − Aη)^{t}Ṽ_{y}^{−1}(Δy − Aη)

The parameters obtained by minimizing χ̃^{2} are given by η = Ṽ_{A}A^{t}Ṽ_{y}^{−1}Δy, where Ṽ_{A} = (A^{t}Ṽ_{y}^{−1}A)^{−1}. The covariance matrix must be calculated from the definition

V_{η} ≡ ⟨δηδη^{t}⟩ = Ṽ_{A}A^{t}Ṽ_{y}^{−1}⟨δyδy^{t}⟩Ṽ_{y}^{−1}AṼ_{A} = Ṽ_{A}A^{t}Ṽ_{y}^{−1}V_{y}Ṽ_{y}^{−1}AṼ_{A}.

This equation collapses to V_{η} = V_{A} when Ṽ_{y}^{−1} = V_{y}^{−1}.

Notice what has been accomplished here: I have produced an estimate of the fitting parameters and their associated covariance matrix without inverting the measurement covariance matrix V_{y}. Furthermore, the estimate is unbiased (the parameters converge to the true values if the experiment is repeated indefinitely) and the covariance matrix — and hence the errors — are known exactly. The price paid for this estimation scheme is that the parameter errors are no longer the best possible; that would require the minimization of the original χ^{2} function. However, if Ṽ_{y}^{−1} is a reasonable representation of V_{y}^{−1} then the “non-optimal” estimate may in fact be close to optimal.

3. Non-Optimal Fitting with Constraints

The effect of constraints can easily be taken into account by the standard Lagrange Multiplier method where the “initial” unconstrained values and their covariance matrix are obtained by the non-optimal fitting procedure. This procedure works exactly the same as if the parameters were obtained by the normal optimal least squares technique. The reason is because constraints act only to move the fitted parameters away from their unconstrained values and the relative motion of each of the parameters is governed by the covariance matrix of the fit, which has absorbed all the information regarding details of the measurement errors and correlations. This point was discussed in Section 5 of the first paper in this series, CBX 91–72.

4. Simple Application to Multiple Scattering

Particles moving through a detector suffer innumerable collisions which alter the trajectory of the particle by a stochastic process. In general, a particle traversing material of length Δx and a radiation length L_{R} will be deflected in any particular plane by an random angle θ whose r.m.s. value can be calculated from

⟨θ^{2}⟩ = (\frac{0.0141}{pβ})^{2} \frac{Δx}{L_{R}} = HΔx,

where p is the momentum in GeV/c and β is the velocity in units of c. Note that this is the width in each of two planes perpendicular to each other and lying along the flight path. I am only accounting for standard multiple scattering effects that can be described by Gaussian statistics; single scattering theory is too complicated to incorporate in least squares fitting, at least by me.

Now consider a particle moving in a solenoidal detector (having a magnetic field in the z direction) from a point near the center and assume that the drift distance to a wire is measured in every layer. If the particle has sufficient momentum so that there is not much bending, we can equate this situation to a simple geometry in which n planar drift chamber are arrayed along the x axis which is defined to be the direction of the particle in the r − ϕ plane (i.e., x is really a measurement of the radius r and y measures the drift distance). The covariance matrix of the measurements can be represented by the sum V_{y} = V_{y0} + V_{yD} + V_{yC}, where V_{y0} = diag(σ_{i}^{2}) is the component due to independent measurement errors, V_{yD} is due to multiple scattering at discrete surfaces and V_{yC} is a result of continuous multiple scattering in gas.

Let’s consider discrete scatterers first. If material of thickness t is present at x = x_{S}, the track will be deflected stochastically from its nominal flight path by an angle θ_{s}, causing the drift distance measurements at plane i to be changed by an amount δy_{i} = θ_{s}(x_{i} − x_{S}). The covariance matrix element (V_{yD})_{ij} can then be calculated from its definition, viz.

(V_{yD})_{ij} = ⟨δy_{i}δy_{j}⟩ = ⟨θ_{s}^{2}⟩(x_{i} − x_{S})(x_{j} − x_{S}) ≡ Ht(x_{i} − x_{S})(x_{j} − x_{S}),

and is 0 if x_{i} ≤ x_{S} or x_{j} ≤ x_{S}. Note that t should reflect the total path length of the track through the scatterer; this is expecially important for tracks having significant components in the z direction.

Continuous scattering is slightly more complicated. Assume the presence of gas in the region x_{G1} ≤ x ≤ x_{G2}, where L = x_{G2} − x_{G1}. If we divide up the path into small regions Δx the total scattering angle and deviation up to a point x in the gas volume is

\begin{aligned}θ & = ∑_{l=1}^{N} δθ_{l} \\ δy & = ∑_{l=1}^{N} θ_{l}Δx = ∑_{l=1}^{N}(N − l)δθ_{l}Δx, \end{aligned}

where Δx = (x − x_{G1})/N. Thus the covariance matrix element for points x_{i} and x_{j} (where x_{j} ≥ x_{i}) in the gas volume is

\begin{aligned}(V_{yG})_{ij} & = ∑_{l=1}^{N_{i}} ∑_{l′=1}^{N_{j}}(N_{i} − l)(N_{j} − l′)⟨δθ_{l}δθ_{l′}⟩Δx^{2} = ∫_{x_{G1}}^{x_{i}} H(x_{i} − x)(x_{j} − x)dx \\ & = H \frac{1}{2}(x_{i} − x_{G1})^{2}(x_{j} − x_{G1} − \frac{1}{3}(x_{i} − x_{G1})), \end{aligned}

where I have used ⟨δθ_{l}δθ_{l′}⟩ = HΔxδ_{ll′}.

If x_{i} and x_{j} both are beyond the gas region the above equations have to be modified to

\begin{aligned}θ & = ∑_{l=1}^{N_{G}} δθ_{l} \\ δy & = ∑_{l=1}^{N_{G}} θ_{l}Δx + θ(x − x_{G2}) = ∑_{l=1}^{N_{G}}(N_{G} − l)δθ_{l}Δx + θ(x − x_{G2}), \end{aligned}

where N_{G} = L/Δx. In this case the covariance matrix is

(V_{yD})_{ij} = HL \{\frac{1}{3}L^{2} + \frac{1}{2}L[(x_{i} − x_{G2}) + (x_{j} − x_{G2})] + (x_{j} − x_{G2})(x_{i} − x_{G2})\} .

If x_{j} is beyond the gas region and x_{i} is inside it the covariance matrix becomes

(V_{yD})_{ij} = H \frac{1}{2}(x_{i} − x_{G1})^{2}(L − \frac{1}{3}(x_{i} − x_{G1}) + (x_{j} − x_{G2})).

x_{i} and x_{j} are interchanged in the above equation if x_{i} is beyond the gas region and x_{j} is inside of it.

5. Summary of Multiple Scattering Formulas

A summary of the contributions to the multiple scattering matrix V_{y} = V_{y0}+V_{yD}+V_{yC} is shown below (remember that x_{j} and x_{i} are defined such that x_{j} ≥ x_{i}).

1. Measurement error:

(V_{y0})_{ij} = σ_{i}^{2}δ_{ij}

where σ_{i} is the measurement error of layer i.

2. Scattering from discrete material (thin) at x = x_{S} and thickness t (remember that H = (0.0141/pβ)^{2}/L_{R} is the square of the r.m.s. scattering angle per unit length in a plane and L_{R} is the radiation length):

(V_{yD})_{ij} = Ht(x_{i} − x_{S})(x_{j} − x_{S}) for x_{i}, x_{j} ≥ x_{S}

= 0 otherwise

3. Scattering in a gas region x_{G1} ≤ x ≤ x_{G2}, L = x_{G1} − x_{G2}:

(V_{yG})_{ij} = H \frac{1}{2}(x_{i} − x_{G1})^{2}(x_{j} − x_{G1} − \frac{1}{3}(x_{i} − x_{G1}))

for x_{G1} ≤ x_{i}, x_{j} ≤ x_{G2}

(V_{yD})_{ij} = H \frac{1}{2}(x_{i} − x_{G1})^{2}(L − \frac{1}{3}(x_{i} − x_{G1}) + (x_{j} − x_{G2}))

for x_{G1} ≤ x_{i} ≤ x_{G2}, x_{j} ≥ x_{G2}

(V_{yD})_{ij} = H \frac{1}{2}(x_{j} − x_{G1})^{2}(L − \frac{1}{3}(x_{j} − x_{G1}) + (x_{i} − x_{G2}))

for x_{G1} ≤ x_{j} ≤ x_{G2}, x_{i} ≥ x_{G2}

(V_{yD})_{ij} = HL \{\frac{1}{3}L2 + \frac{1}{2}L[(x_{i} − x_{G2}) + (x_{j} − x_{G2})] + (x_{j} − x_{G2})(x_{i} − x_{G2})\}

for x_{i}, x_{j} ≥ x_{G2} = 0 otherwise