Drift Chamberのトラッキング解析

Drift chamber (DC) の各面で得られた測定座標から, 粒子の飛跡を求める方法を説明する.

トラックパラメータと測定座標

$z=0$ におけるトラッキング結果の粒子の位置を $(x, y)$ , そのslopeを $a=\dv{x}{z}$ , $b=\dv{y}{z}$ と書く. これらをまとめてパラメータベクトル $\bm{\theta}$ と表す.

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

$i$ 番目の面の $z$ 座標を $z_i$ と表すと, $z_i$ での粒子の飛跡位置 $(x_i, y_i)$ は次のように表せる.

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

DCの各面はある軸方向しか位置が分からない. $i$ 番目の面の軸での測定座標を $m_i$ と表し, これらをまとめて $\bm{m}$ と表す. 各面の軸の角度を $\phi_i$ ( $x$ 軸を基準に取って測る) とすると, $m_i$ は $(x_i, y_i)$ のこの軸への射影であるから, 次のように表せる.

$$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$$

ここで $\phi_i$ はwireの方向ではなく, その面が位置を測定する軸, すなわちwireに垂直な方向の角度である点に注意せよ. 以下の図では, $+y$ を上向き, $+z$ を画面奥向きに取り, 右手系を採用している. 灰色の線はsense wire, 青色の矢印はsense wireに垂直な方向の測定座標 $m_i$ を表す.

実際のDCでは, hitしたwireのIDとdrift timeが測定される. drift timeをdrift length $d_i$ に変換し, hitしたwireの測定軸上の位置を $w_i$ とすると, 粒子の通過位置の候補は

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

と表せる. $s_i$ は粒子がwireのどちら側を通過したかを表す符号である. drift 時間だけからは $s_i$ を決められないため, 各面について2つの $m_i$ が候補となる. この曖昧さが, 後で扱うLeft-Right ambiguityである.

$\bm{m}$ は行列 $H$ を使って次のように表せる.

$$\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}$$

最小二乗法によるトラックフィット

DCは有限な分解能を持つため, これが厳密に成り立つことはなく, 最小二乗解 $\hat{\bm{\theta}}$ を探すというのが現実の操作になる. ここでは, まず各面の測定分解能が同じであると仮定し, 等しい重みでfitする. 残差の二乗和(SSR)は

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

と表せ, 最小二乗解 $\hat{\bm{\theta}}$ は次の条件を満たす.

$$\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$$

したがって, 最小二乗解を得るために解くべき方程式は

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

となる. 以降, $N = H^\mathsf{T} H$ , $\bm{g} = H^\mathsf{T} \bm{m}$ と表すことにする. ところで $N$ は次のような実対称行列となる.

$$\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*}$$

さらに, $\mathrm{rank}(H)=4$ならば$N$は正定値であるから, Cholesky分解を利用して $N \hat{\bm{\theta}} = \bm{g}$を解くことができる.

各面の分解能が異なる場合

各面の測定分散を$\sigma_i^2$とし, $W=\operatorname{diag}(1/\sigma_i^2)$とすると, 重み付き最小二乗法で解くべき正規方程式は

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

となる. すべての$\sigma_i$が同じ場合は, 定数倍を除いて$W=I$とみなせるため, ここで扱う等重みの式に帰着する.

注記

$n \times n$ 実対称行列 $M$ が正定値であるとは, 任意の非零ベクトル $x \in \mathbb{R}^n$ に対して次が成り立つことを意味する.

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

また, 次が成り立つとき, $M$ が半正定値であるという.

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

今の場合, $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$$

となり, 半正定値であることがわかる. 加えて, track parameter $(x,y,a,b)$を一意に決定できる面構成では $\mathrm{rank}(H)=4$となるから,

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

が成り立ち, $N$は実対称正定値行列である. 一方, hit数や面構成が不十分で $\mathrm{rank}(H)<4$となる場合, $N$は特異な半正定値行列となり, 通常のCholesky分解で track parameterを一意に求めることはできない.

Cholesky分解を使った解法

Cholesky分解とは正定値エルミート行列を下三角行列 $L$ と その共役転置 $L^*$ との積に分解することであり, LU分解の特殊な場合である. 今の場合, $N$ は実対称行列なので共役転置はただの転置となり, 次のように分解できる.

$$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分解を使うことのメリットは逆行列を計算せずに, $N \hat{\bm{\theta}} = \bm{g}$ を高速に解けることにある. Cholesky分解の方法については後で取り扱うこととして, まずはCholesky分解できたとして $N \hat{\bm{\theta}} = \bm{g}$ を解く.

まず中間ベクトルとして $\bm{u} = L^\mathsf{T} \hat{\bm{\theta}}$ とおくと, 解くべき方程式は次のような形になる.

$$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 )$$

$L$ の下三角行列であるという性質がここで効力を発揮して, 上から順番に解いていくことができる.

$$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}}$$

一般には, 次のように書ける.

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

これを前進代入と呼ぶ.

$\bm{u}$ が求まったので, 今度は $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 )$$

今度は下から順番に解いていけば良い.

$$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}}$$

一般には

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

これを後退代入と呼ぶ.

こうして最小二乗解 $\hat{\bm{\theta}}$ を得ることができる.

変換行列 $G$ を用いた高速化

以上の手順を毎回行っても良いが, もう少し工夫するとより計算量を減らすことができる. もともと解くべき方程式は次であった.

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

これを $\hat{\bm{\theta}}$ について解くと,

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

となる. ここで, $G = N^{-1} H^\mathsf{T}$ とおいた. このような $G$ を求めておけば, $\bm{m}$ が得られたときに通常の行列計算ですぐに $\hat{\bm{\theta}}$ を得ることができる. Cholesky分解で得た $L$ を使ってこの $G$ を求める.

$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*}$$

となる. この形を見ればわかるように前進代入, 後退代入を各列に対して行うことで $G$ を得ることができる.

次のように $U$ をおく.

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

すると $LU=H^\mathsf{T}$ であるから, 各列に対して前進代入を行うことで $U$ を得ることができる.

$$\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}$$

次に 後退代入で $G$ を求める.

$$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}$$

こうして $G$ を求めることができる. $\bm{m}$ が得られれば, $G$ をこれに乗ずることで, $\hat{\bm{\theta}}$ を得ることができる.

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

Cholesky分解のアルゴリズム

Cholesky分解のアルゴリズムについて述べる. Cholesky分解にはいくつかの方法が存在するが, ここでは Cholesky-Crout法 (列指向) と呼ばれる方法を紹介する.

$$\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*}$$

1列目から取り掛かる. まず $(0,0)$ 成分は

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

$l_{00}$が分かったので, $(1,0)$, $(2,0)$, $(3,0)$ 成分から $l_{10}$, $l_{20}$, $l_{30}$ を得ることができる.

$$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}}}.$$

こうして, 1列目が解けた.

$$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}}} & ? & ? & ? )$$

次に, 2列目に移る. $(1,1)$成分から $l_{11}$ を求めることができる.

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

そして, 同じように $(2,1)$, $(3,1)$ 成分から $l_{21}$, $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}}$$

これを繰り返すことで $L$ を得ることができる.

一般には $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}$$

で求めた後に, 対角成分より下の要素を

$$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

上記の手続きは $\bm{m}$ がわかったものとして解説しているが, 実際にはドリフトタイムからはドリフト長しか分からず, wireの左右どちらを通ったか分からない. このLRの組み合わせを識別しやすくするために, MWDCでは半セルずらした面構成を採用している場合が多い. ここではこのLRを解く手法について解説する.

この左右を解く最もシンプルな方法は全ての面に対してLR両方のパターンを総当たりで計算し, SSRが最小となる組み合わせを採用するというものである. Artemisの既存コードではこの手法が使われている. 適当に模擬コードを書いてみると以下のような感じになる. ここでは, 各面で使用するhitは一つに決まっており, LRだけが未決定であるとする. これをみてわかるように全ての面のLRのパターンに対してSSRを $2^n$ 回計算しており, これがDCのトラッキング上のボトルネックとなっている部分である. これを工夫して高速化する手法についてはまた今度取り扱う.

// 2進数の1をn bit左shiftすると値は2^nになる.
for (unsigned int i = 0; i < (1U << n); ++i) {
    for (int p = 0; p != n; ++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); // θ̂ = Gm
    }
    const double SSR = calcSSR(m, theta);
    if (SSR < minSSR) {
        minSSR = SSR;
        for (int q = 0; q != 4; ++q) {
            best_theta[q] = theta[q];
        }
    }
}