Track Reconstruction in Drift Chambers

This page explains how to reconstruct a particle track from position measurements in the individual planes of a drift chamber (DC).

Track Parameters and Measured Coordinates

Let $(x, y)$ denote the reconstructed particle position at $z=0$, and let $a=\dv{x}{z}$ and $b=\dv{y}{z}$ denote the track slopes. We group these four parameters into a vector $\bm{\theta}$:

$$\bm{\theta} = \mqty( x \\ y \\ a \\ b )$$

At the $i$-th plane, located at $z=z_i$, the track passes through $(x_i, y_i)$, where

$$x_i = x + a z_i, \quad y_i = y + b z_i .$$

Each DC plane measures position along a single axis. Let $m_i$ denote the coordinate measured in the $i$-th plane, and let $\bm{m}$ be the vector of these measurements. If the measurement axis makes an angle $\phi_i$ with the $x$ axis, then $m_i$ is the projection of $(x_i, y_i)$ onto that axis:

$$m_i = x_i \cos \phi_i + y_i \sin \phi_i = (x + a z_i) \cos \phi_i + (y + b z_i) \sin \phi_i$$

Note that $\phi_i$ describes the measurement axis, which is perpendicular to the wire, rather than the wire direction itself. The figure below uses a right-handed coordinate system, with $+y$ pointing upward and $+z$ into the screen. The gray line represents a sense wire, and the blue arrow shows the measured coordinate $m_i$ along the axis perpendicular to it.

In practice, a DC records the ID of the wire that registered a hit and the drift time. We convert the drift time to a drift distance $d_i$ and denote the wire's position along the measurement axis by $w_i$. The two possible particle positions are then

$$m_i = w_i + s_i d_i, \qquad s_i \in \{-1,+1\}$$

The sign $s_i$ indicates which side of the wire the particle passed on. The drift time alone does not determine $s_i$, leaving two possible values of $m_i$ for each plane. This is the left-right ambiguity discussed later.

Using a matrix $H$, we can express the measurements as

$$\bm{m} = \mqty( m_0 \\ m_1 \\ \vdots \\ m_{n-1} ) = \mqty( \cos \phi_0 & \sin \phi_0 & z_0 \cos \phi_0 & z_0 \sin \phi_0 \\ \cos \phi_1 & \sin \phi_1 & z_1 \cos \phi_1 & z_1 \sin \phi_1 \\ \vdots & \vdots & \vdots & \vdots\\ \cos \phi_{n-1} & \sin \phi_{n-1} & z_{n-1} \cos \phi_{n-1} & z_{n-1} \sin \phi_{n-1} \\ ) \mqty( x \\ y \\ a \\ b ) = H \bm{\theta}$$

Least-Squares Track Fitting

Because a DC has finite measurement resolution, the relation above does not hold exactly for the measured coordinates. We therefore find the least-squares estimate $\hat{\bm{\theta}}$. We first assume that all planes have the same resolution and give each measurement equal weight. The sum of squared residuals (SSR) is

$$\mathrm{SSR} = \sum_i \qty(m_i - \sum_p H_{ip} \theta_p)^2$$

and the least-squares solution $\hat{\bm{\theta}}$ satisfies

$$\eval{\pdv{\mathrm{SSR}}{\theta_q}}_{\bm{\theta}=\hat{\bm{\theta}}} = -2 \sum_i H_{iq}\qty( m_i - \sum_p H_{ip} \hat{\theta}_p ) = 0$$ $$\Rightarrow \sum_i H_{iq} m_i = \sum_p \sum_i H_{iq} H_{ip} \hat{\theta}_p$$ $$\Rightarrow \qty(H^\mathsf{T} \bm{m})_q = \sum_p \qty(H^\mathsf{T} H)_{qp} \hat{\theta}_p$$

This gives the normal equations:

$$H^\mathsf{T} \bm{m} = H^\mathsf{T} H \hat{\bm{\theta}} .$$

For the rest of this discussion, we write $N = H^\mathsf{T} H$ and $\bm{g} = H^\mathsf{T} \bm{m}$. The matrix $N$ is real and symmetric:

$$\begin{align*} N &= H^\mathsf{T} H \\ &= \mqty( \cos \phi_0 & \cos \phi_1 & \cdots & \cos \phi_{n-1} \\ \sin \phi_0 & \sin \phi_1 & \cdots & \sin \phi_{n-1} \\ z_0 \cos \phi_0 & z_1 \cos \phi_1 & \cdots & z_{n-1} \cos \phi_{n-1} \\ z_0 \sin \phi_0 & z_1 \sin \phi_1 & \cdots & z_{n-1} \sin \phi_{n-1} ) \mqty( \cos \phi_0 & \sin \phi_0 & z_0 \cos \phi_0 & z_0 \sin \phi_0 \\ \cos \phi_1 & \sin \phi_1 & z_1 \cos \phi_1 & z_1 \sin \phi_1 \\ \vdots & \vdots & \vdots & \vdots\\ \cos \phi_{n-1} & \sin \phi_{n-1} & z_{n-1} \cos \phi_{n-1} & z_{n-1} \sin \phi_{n-1} \\ )\\ &= \mqty( \sum_i (\cos \phi_i)^2 & \sum_i \cos \phi_i \sin \phi_i & \sum_i z_i (\cos \phi_i)^2 & \sum_i z_i \cos \phi_i \sin \phi_i \\ \sum_i \cos \phi_i \sin \phi_i & \sum_i (\sin \phi_i)^2 & \sum_i z_i \cos \phi_i \sin \phi_i & \sum_i z_i (\sin \phi_i)^2 \\ \sum_i z_i (\cos \phi_i)^2 & \sum_i z_i \cos \phi_i \sin \phi_i & \sum_i (z_i \cos \phi_i)^2 & \sum_i z_i^2 \cos \phi_i \sin \phi_i \\ \sum_i z_i \cos \phi_i \sin \phi_i & \sum_i z_i (\sin \phi_i)^2 & \sum_i z_i^2 \cos \phi_i \sin \phi_i & \sum_i (z_i \sin \phi_i)^2 ) \end{align*}$$

If $\mathrm{rank}(H)=4$, then $N$ is positive definite. Therefore, $N \hat{\bm{\theta}} = \bm{g}$ can be solved using Cholesky decomposition.

Planes with different resolutions

Let $\sigma_i^2$ be the measurement variance in the $i$-th plane, and define $W=\operatorname{diag}(1/\sigma_i^2)$. The weighted least-squares normal equations are then

$$H^\mathsf{T} W H\hat{\bm{\theta}} = H^\mathsf{T}W\bm{m} .$$

If all $\sigma_i$ are equal, $W$ is a constant multiple of the identity matrix, and these equations reduce to the equal-weight form used here.

note

An $n \times n$ real symmetric matrix $M$ is positive definite if the following holds for any nonzero vector $x \in \mathbb{R}^n$:

$$\bm{x}^\mathsf{T} M \bm{x} > 0$$

It is called positive semidefinite if

$$\bm{x}^\mathsf{T} M \bm{x} \geq 0 .$$

In this case, since $N = H^\mathsf{T}H$,

$$\bm{x}^\mathsf{T} N \bm{x} = \bm{x}^\mathsf{T} H^\mathsf{T} H \bm{x} = (H\bm{x})^\mathsf{T} (H\bm{x}) = \| H \bm{x} \|^2 \geq 0$$

so $N$ is positive semidefinite. If the plane configuration allows the track parameters $(x,y,a,b)$ to be determined uniquely, then $\mathrm{rank}(H)=4$ and

$$\| H\bm{x} \|^2 > 0$$

for every nonzero $\bm{x}$. Thus, $N$ is real, symmetric, and positive definite. If there are too few hits or the plane configuration does not provide enough independent measurements, then $\mathrm{rank}(H)<4$. In that case, $N$ is singular and positive semidefinite, and standard Cholesky decomposition cannot yield a unique set of track parameters.

Solving with Cholesky Decomposition

Cholesky decomposition factors a positive-definite Hermitian matrix into the product of a lower triangular matrix $L$ and its conjugate transpose $L^*$. It can be regarded as a special case of LU decomposition. Here, $N$ is real and symmetric, so the conjugate transpose reduces to the transpose:

$$N = L L^\mathsf{T} = \mqty( l_{00} & 0 & 0 & 0 \\ l_{10} & l_{11} & 0 & 0 \\ l_{20} & l_{21} & l_{22} & 0 \\ l_{30} & l_{31} & l_{32} & l_{33} ) \mqty( l_{00} & l_{10} & l_{20} & l_{30} \\ 0 & l_{11} & l_{21} & l_{31} \\ 0 & 0 & l_{22} & l_{32} \\ 0 & 0 & 0 & l_{33} )$$

Cholesky decomposition lets us solve $N \hat{\bm{\theta}} = \bm{g}$ efficiently without explicitly computing a matrix inverse. We will describe how to compute $L$ later. For now, assume that $L$ is known and use it to solve the system.

Define an intermediate vector $\bm{u} = L^\mathsf{T} \hat{\bm{\theta}}$. Then the equation to solve becomes

$$N \hat{\bm{\theta}} = L L^\mathsf{T} \hat{\bm{\theta}} = L \bm{u} = \bm{g}$$ $$\Rightarrow \mqty( l_{00} & 0 & 0 & 0 \\ l_{10} & l_{11} & 0 & 0 \\ l_{20} & l_{21} & l_{22} & 0 \\ l_{30} & l_{31} & l_{32} & l_{33} ) \mqty( u_0 \\ u_1 \\ u_2 \\ u_3 ) = \mqty( g_0 \\ g_1 \\ g_2 \\ g_3 )$$

Because $L$ is lower triangular, we can solve for the components of $\bm{u}$ in order, starting with the first row:

$$l_{00} u_0 = g_0 \quad \Rightarrow \quad u_0 = \frac{g_0}{l_{00}}$$ $$l_{10} u_{0} + l_{11} u_{1} = g_1 \quad \Rightarrow \quad u_1 = \frac{g_1 - l_{10} u_{0}}{l_{11}}$$ $$l_{20} u_{0} + l_{21} u_{1} + l_{22} u_{2} = g_2 \quad \Rightarrow \quad u_2 = \frac{g_2 - l_{21} u_1 - l_{20} u_{0}}{l_{22}}$$ $$l_{30} u_{0} + l_{31} u_{1} + l_{32} u_{2} + l_{33} u_{3} = g_3 \quad \Rightarrow \quad u_3 = \frac{g_3 - l_{32} u_{2} - l_{31} u_1 - l_{30} u_{0}}{l_{33}}$$

In general,

$$u_i = \frac{ g_i - \sum_{j<i} l_{ij} u_j }{ l_{ii} }$$

This is called forward substitution.

With $\bm{u}$ known, we next solve $L^{\mathsf{T}} \hat{\bm{\theta}} = \bm{u}$:

$$\mqty( l_{00} & l_{10} & l_{20} & l_{30} \\ 0 & l_{11} & l_{21} & l_{31} \\ 0 & 0 & l_{22} & l_{32} \\ 0 & 0 & 0 & l_{33} ) \mqty( \hat{\theta}_0 \\ \hat{\theta}_1 \\ \hat{\theta}_2 \\ \hat{\theta}_3 ) = \mqty( u_0 \\ u_1 \\ u_2 \\ u_3 )$$

This time, we start with the last row and work upward:

$$l_{33} \hat{\theta}_3 = u_3 \quad \Rightarrow \quad \hat{\theta}_3 = \frac{u_3}{l_{33}}$$ $$l_{22} \hat{\theta}_2 + l_{32} \hat{\theta}_3 = u_2 \quad \Rightarrow \quad \hat{\theta}_2 = \frac{u_2 - l_{32} \hat{\theta}_3}{l_{22}}$$ $$l_{11} \hat{\theta}_1 + l_{21} \hat{\theta}_2 + l_{31} \hat{\theta}_3 = u_1 \quad \Rightarrow \quad \hat{\theta}_1 = \frac{u_1 - l_{21} \hat{\theta}_2 - l_{31} \hat{\theta}_3}{l_{11}}$$ $$l_{00} \hat{\theta}_0 + l_{10} \hat{\theta}_1 + l_{20} \hat{\theta}_2 + l_{30} \hat{\theta}_3 = u_0 \quad \Rightarrow \quad \hat{\theta}_0 = \frac{u_0 - l_{10} \hat{\theta}_1 - l_{20} \hat{\theta}_2 - l_{30} \hat{\theta}_3}{l_{00}}$$

In general,

$$\hat{\theta}_i = \frac{ u_i - \sum_{j>i} l_{ji} \hat{\theta}_j }{ l_{ii} }.$$

This is called back substitution.

Together, these two steps give the least-squares solution $\hat{\bm{\theta}}$.

Faster Fitting with a Transformation Matrix $G$

We could repeat the procedure above for every fit, but we can reduce the work by precomputing a transformation matrix. Recall that the system to solve is

$$N\hat{\bm{\theta}} = H^\mathsf{T} \bm{m}.$$

Solving this for $\hat{\bm{\theta}}$ gives

$$\hat{\bm{\theta}} = N^{-1} H^\mathsf{T} \bm{m} = G \bm{m}.$$

Here, $G = N^{-1} H^\mathsf{T}$. Once $G$ has been computed, each new measurement vector $\bm{m}$ requires only a matrix-vector multiplication to obtain $\hat{\bm{\theta}}$. We can compute $G$ using the Cholesky factor $L$.

Since $N = LL^\mathsf{T}$,

$$NG = LL^\mathsf{T} G = H^\mathsf{T}$$ $$\begin{align*} \Rightarrow \mqty( l_{00} & 0 & 0 & 0 \\ l_{10} & l_{11} & 0 & 0 \\ l_{20} & l_{21} & l_{22} & 0 \\ l_{30} & l_{31} & l_{32} & l_{33} ) \mqty( l_{00} & l_{10} & l_{20} & l_{30} \\ 0 & l_{11} & l_{21} & l_{31} \\ 0 & 0 & l_{22} & l_{32} \\ 0 & 0 & 0 & l_{33} ) &\mqty( G_{00} & G_{01} & \cdots & G_{0 n-1} \\ G_{10} & G_{11} & \cdots & G_{1 n-1} \\ G_{20} & G_{21} & \cdots & G_{2 n-1} \\ G_{30} & G_{31} & \cdots & G_{3 n-1} ) \\ &= \mqty( \cos \phi_0 & \cos \phi_1 & \cdots & \cos \phi_{n-1} \\ \sin \phi_0 & \sin \phi_1 & \cdots & \sin \phi_{n-1} \\ z_0 \cos \phi_0 & z_1 \cos \phi_1 & \cdots & z_{n-1} \cos \phi_{n-1} \\ z_0 \sin \phi_0 & z_1 \sin \phi_1 & \cdots & z_{n-1} \sin \phi_{n-1} ) \end{align*}$$

We can therefore compute $G$ by applying forward and back substitution to each column.

Define $U$ as

$$L^\mathsf{T} G = U = (\bm{u}_0, \bm{u}_1, \cdots, \bm{u}_{n-1}) .$$

Since $LU=H^\mathsf{T}$, we obtain $U$ by solving for each column using forward substitution:

$$\mqty( l_{00} & 0 & 0 & 0 \\ l_{10} & l_{11} & 0 & 0 \\ l_{20} & l_{21} & l_{22} & 0 \\ l_{30} & l_{31} & l_{32} & l_{33} ) (\bm{u}_0, \bm{u}_1, \cdots, \bm{u}_{n-1}) = H^\mathsf{T}$$ $$U_{ij} = \begin{cases} \frac{ H_{j0} }{ L_{00} } & (i=0) \\ \frac{ H_{ji} - \sum_{k=0}^{i-1} L_{ik} U_{kj} }{ L_{ii} } & (i>0) \end{cases}$$

Next, we obtain $G$ by back substitution:

$$L^\mathsf{T} G = U$$ $$\mqty( l_{00} & l_{10} & l_{20} & l_{30} \\ 0 & l_{11} & l_{21} & l_{31} \\ 0 & 0 & l_{22} & l_{32} \\ 0 & 0 & 0 & l_{33} ) \mqty( G_{00} & G_{01} & \cdots & G_{0 n-1} \\ G_{10} & G_{11} & \cdots & G_{1 n-1} \\ G_{20} & G_{21} & \cdots & G_{2 n-1} \\ G_{30} & G_{31} & \cdots & G_{3 n-1} ) = (\bm{u}_0, \bm{u}_1, \cdots, \bm{u}_{n-1})$$ $$G_{ij} = \begin{cases} \frac{ U_{ij} - \sum_{k=i+1}^{3} L_{ki} G_{kj} }{ L_{ii} } & (i<3) \\ \frac{ U_{3j} }{ L_{33} } & (i=3) \end{cases}$$

With $G$ computed, we obtain the fitted track parameters for a measurement vector $\bm{m}$ as

$$\hat{\bm{\theta}} = G \bm{m}$$

Cholesky Decomposition Algorithm

There are several algorithms for Cholesky decomposition. Here we use the Cholesky-Crout method, which computes $L$ one column at a time.

$$\begin{align*} N &= L L^\mathsf{T} \\ &= \mqty( l_{00} & 0 & 0 & 0 \\ l_{10} & l_{11} & 0 & 0 \\ l_{20} & l_{21} & l_{22} & 0 \\ l_{30} & l_{31} & l_{32} & l_{33} ) \mqty( l_{00} & l_{10} & l_{20} & l_{30} \\ 0 & l_{11} & l_{21} & l_{31} \\ 0 & 0 & l_{22} & l_{32} \\ 0 & 0 & 0 & l_{33} ) \\ &= \mqty( l_{00}^2 & l_{00}l_{10} & l_{00}l_{20} & l_{00}l_{30} \\ l_{00}l_{10} & l_{10}^2+l_{11}^2 & l_{10}l_{20}+l_{11} l_{21} & l_{10}l_{30}+l_{11}l_{31} \\ l_{00}l_{20} & l_{10}l_{20}+l_{11}l_{21} & l_{20}^2+l_{21}^2+l_{22}^2 & l_{20}l_{30}+l_{21}l_{31}+l_{22}l_{32} \\ l_{00}l_{30} & l_{10}l_{30}+l_{11}l_{31} & l_{20}l_{30}+l_{21}l_{31}+l_{22}l_{32} & l_{30}^2+l_{31}^2+l_{32}^2+l_{33}^2 ) \end{align*}$$

Start with the first column. The $(0,0)$ entry gives

$$N_{00} = l_{00}^2 \quad \Rightarrow \quad l_{00} = \sqrt{N_{00}} .$$

With $l_{00}$ known, we obtain $l_{10}$, $l_{20}$, and $l_{30}$ from the $(1,0)$, $(2,0)$, and $(3,0)$ entries:

$$N_{10} = l_{00} l_{10}, \quad N_{20} = l_{00} l_{20}, \quad N_{30} = l_{00} l_{30} .$$ $$\Rightarrow l_{10} = \frac{N_{10}}{\sqrt{N_{00}}}, \quad l_{20} = \frac{N_{20}}{\sqrt{N_{00}}}, \quad l_{30} = \frac{N_{30}}{\sqrt{N_{00}}}.$$

This gives the first column:

$$L = \mqty( \sqrt{N_{00}} & 0 & 0 & 0\\ \frac{N_{10}}{\sqrt{N_{00}}} & ? & 0 & 0\\ \frac{N_{20}}{\sqrt{N_{00}}} & ? & ? & 0\\ \frac{N_{30}}{\sqrt{N_{00}}} & ? & ? & ? )$$

Next, move to the second column. The $(1,1)$ entry gives $l_{11}$:

$$N_{11} = l_{10}^2 + l_{11}^2$$ $$l_{11} = \sqrt{N_{11} - l_{10}^2}$$

Similarly, the $(2,1)$ and $(3,1)$ entries give $l_{21}$ and $l_{31}$:

$$N_{21} = l_{10}l_{20}+l_{11}l_{21}, \quad N_{31} = l_{10}l_{30}+l_{11}l_{31}.$$ $$l_{21} = \frac{N_{21} - l_{10}l_{20}}{l_{11}}, \quad l_{31} = \frac{N_{31} - l_{10}l_{30}}{l_{11}} .$$

Repeating this procedure gives $L$.

In general, we first compute the diagonal entry in column $j$:

$$L_{jj} = \begin{cases} \sqrt{ N_{00} } & (j=0) \\ \sqrt{ N_{jj} - \sum_{k=0}^{j-1}L_{jk}^2 } & (j>0) \end{cases}$$

We then compute the entries below it:

$$L_{ij} = \begin{cases} \frac{ N_{ij} }{ L_{jj} } & (j = 0) \\ \frac{ N_{ij} - \sum_{k=0}^{j-1} L_{ik}L_{jk} }{ L_{jj} } & (j > 0) \end{cases}$$

Left-Right Ambiguity

The fitting procedure above assumes that $\bm{m}$ is known. In practice, the drift time gives only the drift distance, leaving two possible positions on opposite sides of the wire. To help resolve this left-right ambiguity, multiwire drift chambers (MWDCs) often use planes whose wire positions are offset by half a cell relative to one another. We now consider how to select the correct combination of left-right assignments.

The simplest approach is to fit a track for every possible combination of left-right assignments and choose the one with the smallest SSR. This is the method used in the existing Artemis code. The pseudocode below illustrates the procedure, assuming that one hit has already been selected in each plane and only the left-right assignments remain unknown. For $n$ planes, this requires evaluating the SSR for all $2^n$ combinations, making it a bottleneck in DC tracking. Ways to speed up this step will be covered separately.

// A left shift by n bits gives 2^n combinations.
for (unsigned int i = 0; i < (1U << n); ++i) {
    for (int p = 0; p != n; ++p) {
        // Bit p selects the left-right sign for plane p.
        int lr = ((i >> p) & 1) ? 1 : -1;
        m[p] = cell_size * (wire_id[p] - center) + drift_length[p] * lr;
    }
    for (int q = 0; q != 4; ++q) {
        theta[q] = inner_product(G[q], m); // theta_hat = Gm
    }
    const double SSR = calcSSR(m, theta);
    if (SSR < minSSR) {
        minSSR = SSR;
        for (int q = 0; q != 4; ++q) {
            best_theta[q] = theta[q];
        }
    }
}