Paul Avery CBX 91 − 73 October 18, 1991

Applied Fitting Theory II Determining Systematic Effects by Fitting

Section 1: Introduction Section 2: Finding Systematic Effects Using Unconstrained Fits Section 3: Finding Systematic Effects Using Constrained Fits Section 4: Tracking Example Section 5: Limitations of This Algorithm

1. Introduction

This is the second paper in my thrilling series on fitting theory. In this note I will concentrate on determining systematic contributions to measurement error. These systematic effects can degrade the resolution of the parameters beyond what is reported in the statistical error. Examples of effects that are commonly found in track fitting are misalignment of detector positions and orientations, nonuniformity of the magnetic field, noncircularity of a drift chamber layer, unaccounted bowing in the endplates (which affects stereo angles), extra wire sag, twisting of the endplates relative to each other, and so forth.

I present here a method which allows one to determine these systematic effects provided that one can parametrize the effect of any given parameter on the measurement. The heart of the technique lies in performing many thousands of fits (constrained or unconstrained), each of which depends on parameters local to that fit together with the unknown parameters describing the systematic effect. One naively would expect this procedure to be hopeless because the χ^{2} minimization would lead to thousands of equations to be solved simultaneously. It turns out, however, that one can remove the parameters describing the many distinct fits and only determine the small number of remaining systematic unknowns. Better still, one can accumulate statistics steadily, without any intermediate matrix inversion, in order to reduce the errors in these parameters.

Note: I make no claim to be the only person ever to have invented this algorithm. My treatment is very general, however, since it permits the calculation of the full covariance matrix of the systematic parameters. I also include the interesting case where the individual fits were performed with constraints.

2. Finding Systematic Effects Using Unconstrained Fits

In the first paper in the series (CBX 91−72), I gave a general overview of fitting theory using matrix notation. I showed that if the dependence of the measurements on the parameters can be linearized, then the solution which is best from the point of view of minimum variance, can be obtained rather easily, even in the presence of constraints, by minimizing the χ^{2} function,

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}))^{t}V_{y}^{−1} (y − f(α_{A}) − A(α_{A} − α))
≡ (Δy − Aη)^{t}V_{y}^{−1} (Δy − Aη)

where Δy = y−f(α_{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}.

Now suppose we have many hundreds, or even thousands, of such fits and that each fit l depends on parameters η_{l} which are local to it as well as some unknown parameters v (s of them) which affect all fits (note: when I use the word “fit” I refer to any complete fit, for example a single track or even a set of tracks within an event). Our goal is to find v. For this case the χ^{2} can be written as the sum

χ^{2} = ∑_{l}(Δy_{l} − A_{l}η_{l} − B_{l}v)^{t}V_{y l}^{−1}(Δy_{l} − A_{l}η_{l} − B_{l}v),

where I have assumed that the parameters v are small enough that their effect on the measurements y_{l} can be expressed linearly as δy_{l} = B_{l}v, with B_{l} a n×s matrix. Minimizing this total χ^{2} function with respect to the variables η and v we obtain the set of equations

η_{l} :  A_{l}^{t}V_{y}l^{−1}(Δy_{l} − A_{l}η_{l} − B_{l}v) = 0,
v :  ∑_{l} B_{l}^{t}V_{y}l^{−1}(Δy_{l} − A_{l}η_{l} − B_{l}v) = 0.

Although the number of equations is huge, straightforward algebra gives a surprisingly simple solution for v and its covariance matrix V_{v} in terms of two matrices, s and S:

\begin{aligned}v & = S^{−1}s \\ V_{v} & = S^{−1} \end{aligned}

where s and S are defined as

s = ∑_{l} B_{l}^{t}V_{y l}^{−1}(Δy_{l} − A_{l}η_{l0}) ≡ ∑_{l} B_{l}^{t}V_{y l}^{−1}r_{l}
S = ∑_{l}(B_{l}^{t}V_{y l}^{−1}B_{l} − U_{l}^{t}V_{y l}^{−1}U_{l})
η_{l0} = V_{A l}A_{l}^{t}V_{y l}^{−1}Δy_{l}
η_{l} = η_{l0} − V_{A l}U_{l}v.

The auxiliary matrices V_{A l} and U_{l} are defined as V_{A l} = (A_{l}^{t}V_{y l}^{−1}A_{l})^{−1} and U_{l} = A_{l}^{t}V_{y l}^{−1}B_{l}. Note that the equation for s involves the term r_{l} = Δy_{l} − A_{l}η_{l0}, which is the just the vector of residuals for each individual fit l, ignoring v.

The end result is that the parameters v can be obtained by performing a large number of individual fits (ignoring the effect of v) and then accumulating the results of each fit in two matrices, s and S. Recall that the word “accumulating” refers to a sum over complete fits. In high energy physics language this normally means a sum over individual tracks, but it could even be a sum over events if all the tracks within an event were fit together. The sequence of steps that should be followed is shown below.

1. For each fit in which the data values y_{l} and coefficients A_{l} are given, compute Δy = y − f(α_{A}), and V_{A l}, solve for η_{0l} as shown above and compute the residuals r_{l}. This step assumes nothing about the effect of the v parameters.

2. Compute B_{l} and U_{l} and accumulate over some sample of fits the matrices s = ∑_{l} B_{l}^{t}V_{y l}^{−1}r_{l} and S = ∑_{l}(B_{l}^{t}V_{y l}^{−1}B_{l} − U_{l}^{t}V_{y l}^{−1}U_{l}).

3. After accumulating the above matrices for awhile, occasionally invert S to get V_{v} (remember that V_{v} = S^{−1}) to see if the errors on v are small enough.

4. If the errors on v are not small enough, go to step 1 and accumulate more data. If enough data has been taken, solve for v = S^{−1}s.

Notice that you will not be able to compute the corrected values of η_{l} during the accumulation since they depend on the value of v. Usually this is not important since what is most often desired is v and its covariance matrix V_{v}. Of course, if you want the best parameters η_{l} for each individual fit l, you can always go back and reanalyze the data using the now known values of v via the correction η_{l} = η_{l0} − V_{A l}U_{l}v.

3. Finding Systematic Effects Using Constrained Fits

It is fairly easy to extend the results of the previous section to the case where the parameters v describing the systematic effects must be determined using fits which were performed with constraints. This situation occasionally arises in high energy physics when individual tracks carry insufficient information to determine a systematic parameter but a set of tracks satisfying a global constraint does. For example, if the lower and upper hemispheres of a drift chamber were displaced relative to one another, the shift would be in principle unobservable in individual tracks because each track would lie completely within a single hemisphere. However, the shift would have an effect if one constrained the dimuons in e^{+}e^{−} → μ^{+}μ^{−} events to meet in the middle.

Let the r constraint equations for the l^{th} fit be written H_{l}(α_{l}) = 0. Expanding about α_{A} as before, we obtain the linearized equations H_{l}(α_{A l}) + (α_{l} − α_{A l})∂H_{l}(α_{A l})/∂α_{l} ≡ D_{l}η_{l} + d_{l} = 0, using obvious notation, where D_{l} is a r × m matrix and d_{l} is a vector of length r. Note that the number of constraints r can vary from fit to fit. The full χ^{2} equation can now be written, using the method of Lagrange multipliers, as

χ^{2} = ∑_{l}(Δy_{l} − A_{l}η_{l} − B_{l}v)^{t}V_{y l}^{−1}(Δy_{l} − A_{l}η_{l} − B_{l}v) + ∑_{l} 2λ_{l}^{t}(D_{l}η_{l} + d_{l}),

which, when minimized with respect to the parameters η_{l}, v and λ_{l}, yields the equations

η_{l} :  − A_{l}^{t}V_{y l}^{−1}(Δy_{l} − A_{l}η_{l} − B_{l}v) + D_{l}^{t}λ_{l} = 0
v :  ∑_{l} B_{l}^{t}V_{y l}^{−1}(Δy_{l} − A_{l}η_{l} − B_{l}v) = 0
λ_{l} :  D_{l}η_{l} + d_{l} = 0

Again, straightforward matrix algebra yields the solution for v and its covariance matrix V_{v} in terms of the accumulating matrices s and S:

\begin{aligned}v & = S^{−1}s, \\ V_{v} & = S^{−1}, \end{aligned}

where s and S are given by

s = ∑_{l} B_{l}^{t}V_{y l}^{−1}(Δy_{l} − A_{l}η′_{l0}) ≡ ∑_{l} B_{l}^{t}V_{y l}^{−1}r_{l}
S = ∑_{l}(B_{l}^{t}V_{y l}^{−1}B_{l} − U_{l}^{t}V_{A l}U_{l} + W_{l}^{t}V_{D l}W_{l})
η′_{l0} = η_{l0} − V_{A l}D_{l}^{t}λ_{l0}
η_{l0} = V_{A l}A_{l}^{t}V_{y l}^{−1}Δy_{l}
λ_{l0} = V_{D l}(D_{l}η_{l0} + d_{l})
λ_{l} = λ_{l0} − V_{D l}W_{l}v
η_{l} = η_{l0} − V_{A l}D_{l}^{t}λ_{l} − V_{A l}U_{l}v.

The auxiliary matrices are given by V_{A l} = (A_{l}^{t}V_{y l}^{−1}A_{l})^{−1}, U_{l} = A_{l}^{t}V_{y l}^{−1}B_{l}, V_{D l} = (D_{l}V_{A l}D_{l}^{t})^{−1} and W_{l} = D_{l}V_{A l}U_{l}. Note that η_{l0} and η′_{l0} are, respectively, the unconstrained and constrained solutions for the individual fits obtained without knowledge of the systematic parameters v.

The solution should be obtained using the following sequence of steps.

1. For each fit in which the data values y_{l} and coefficients A_{l} are given, compute Δy = y − f(α_{A}) and V_{A l} and solve for the unconstrained parameters η_{0l} as shown above. The effects of v are ignored here.

2. Determine V_{D l}, solve for λ_{l0} and the constrained parameters η′_{l0} and compute the residuals r_{l}. The effects of v still play no role here.

3. Compute B_{l}, U_{l}, W_{l} and accumulate over some sample of fits the two matrices s = ∑_{l} B_{l}^{t}V_{y l}^{−1}r_{l} and S = ∑_{l}(B_{l}^{t}V_{y l}^{−1}B_{l} − U_{l}^{t}V_{y l}^{−1}U_{l} + W_{l}^{t}V_{D l}W_{l}).

4. After accumulating the above matrices for awhile, occasionally invert S to get V_{v} (remember that V_{v} = S^{−1}) to see if the errors on v are small enough.

5. If the errors on v are not small enough, go to step 1 and accumulate more data. If enough data has been taken, solve for v = S^{−1}s.

As mentioned in the last section, you will not be able to compute the corrected values of η_{l} during the accumulation since they depend on the value of v. However, after determining v you can go back and reanalyze the data.

4. Tracking example

This example is derived from a real application in the CLEO central tracking detector. The central detector consisted at one time of three separate pieces: (1) a 3 layer microvertex chamber consisting of straws; (2) a 10 layer gaseous vertex detector; and (3) a 51 layer drift chamber. The large drift chamber defined the coordinate system and the other two chambers had no z information available on the anode wires. Since the accuracy of chamber placement was far inferior to the internal machining of each chamber, 10 positioning parameters, 5 for each chamber, had to be determined:

1. Offset in x

2. Offset in y

3. Rotation in the x − z plane

4. Rotation in the y − z plane

5. Rotation in the x − y plane

The z offset could not be determined because there was no z information on the anode wires of the inner chambers.

The quantity that is measured by all the drift chambers is the drift distance, which is the distance of closest approach of the track to the wire. Let ϕ_{k} be the azimuthal position at which the track crosses the radial position of the wire and α_{k} be the crossing angle in the r − ϕ plane (a perfectly radial track would have α_{k} = 0). Let v_{1−5} be the positioning parameters for the first chamber and v_{6−10} be the parameters for the second chamber. If we call δd_{k} the change in the drift distance for layer k due to v_{i}, then we can write the contribution to the drift distance as

δd_{k} = ∑_{i} B_{ki}v_{i}.

For 1 ≤ k ≤ 3 (the first chamber), we can determine B_{ki} to be

B_{k1} = sin ϕ_{k} cos α_{k}
B_{k2} = − cos ϕ_{k} cos α_{k}
B_{k3} = z_{k} sin ϕ_{k} cos α_{k}
B_{k4} = −z_{k} cos ϕ_{k} cos α_{k}
B_{k5} = −r_{k} cos α_{k}
B_{ki} = 0 i ≥ 6,

where r_{k} and z_{k} are, respectively, the radius of the layer and z position of the hit (as determined by an initial, approximate fit). For 4 ≤ k ≤ 13 (the second chamber), B_{ki} is

B_{k6} = sin ϕ_{k} cos α_{k}
B_{k7} = − cos ϕ_{k} cos α_{k}
B_{k8} = z_{k} sin ϕ_{k} cos α_{k}
B_{k9} = −z_{k} cos ϕ_{k} cos α_{k}
B_{k10} = −r_{k} cos α_{k}
B_{ki} = 0 i ≤ 5,

For 14 ≤ k ≤ 64 (the third chamber), B_{ki} = 0.

The procedure for finding these 10 constants is straightforward. For each track in the data sample, we perform a fit using all three chambers to determine the track parameters η_{0} and their covariance matrix V_{A}. Each track is then traced through the detector to compute the values of r_{k} (the residual), z_{k}, r_{k}, ϕ_{k} and α_{k} for each layer. Finally, we calculate the matrix B_{ki} and accumulate the two matrices s and S over the data sample. After sufficient data has been gathered, we compute the 10 constants v = S^{−1}s and their covariance matrix V_{v} = S^{−1}. Note that this approach requires only a single pass over the data. By comparison, the conventional approach in which highly selected data (i.e., tracks lying close to the x or y axes) are used to isolate each individual parameter requires several iterations to converge properly.

5. Limitations of This Algorithm

I should comment before closing this paper that determining drift chamber positioning constants with this χ^{2} technique can lead to problems if you try to push it too far in terms of accuracy, especially if the constants are highly correlated. The reason is that once you get to the few micron level, many new effects arise that can partially mimic the ones you are trying to determine. These effects can systematically skew the fit depending on the precise set of tracks you pick. In CLEO, for example, there is a correction to the drift measurement due to pulse height saturation that is a function of θ, the dip angle. Unfortunately, θ is highly correlated with the z value of the track at any layer, causing problems in determining the tilts of chambers since the effect of tilt on drift distance is linear in z.

Another example from CLEO is illustrative. Layer by layer measurements of drift distance have an E × B correction that has the net effect of making the time–distance (t − d) relation asymmetric, i.e., left is different from right. This effect interferes with the ability to precisely determine the x − y rotation angle of one chamber relative to another. If it were symmetric, then the fact that there would be roughly an equal number of tracks on the left or right of any cell would cause any uncertainty in the t−d relation to cancel when determining the rotation angle. A possible solution is to avoid using the innermost or outermost layers (since these have the largest asymmetry) and only use drift distances well within the cell, so that E × B effects are minimized.

Avoidance of systematic bias can be tricky. For example, if you plan to use Bhabha events (e^{+}e^{−} → e^{+}e^{−}) to measure chamber positions and orientations, be aware of the fact that Bhabha events are primarily forward scattered, meaning that positive charged tracks tend to have a particular sign of p_{z} and negative tracks another. This can cause great problems in determing chamber positions and orientations if there are systematic effects which are charge dependent. Suppose, for example, that the drift chamber endplates are slightly rotated with respect to one another. This causes the curvature of tracks to be slightly shifted systematically, resulting in a change of momentum. If positive tracks along the +z direction get an increase in momentum, then they are reduced in momentum along the −z direction. The opposite behavior holds for negative tracks. Now, if an equal number of positive and negative tracks are used in each direction (for instance by using e^{+}e^{−} → μ^{+}μ^{−} events) then the effect cancels. It does not cancel for Bhabha events because of the aforementioned asymmetry.