Charged particle tracking without magnetic field: optimal measurement of track momentum by a Bayesian analysis of the multiple measurements of deflections due to multiple scattering

TL;DR

This paper develops an optimal Bayesian method for measuring particle track momentum in zero magnetic field environments, enhancing detector design and performance for neutrino and gamma-ray experiments.

Contribution

It introduces a Bayesian analysis approach to estimate track momentum from multiple scattering measurements, improving upon existing methods without magnetic fields.

Findings

Optimal estimator for track momentum derived

Enhanced precision in track parameter measurement

Guidance for designing advanced particle detectors

Abstract

We revisit the precision of the measurement of track parameters (position, angle) with optimal methods in the presence of detector resolution, multiple scattering and zero magnetic field. We then obtain an optimal estimator of the track momentum by a Bayesian analysis of the filtering innovations of a series of Kalman filters applied to the track. This work could pave the way to the development of autonomous high-performance gas time-projection chambers (TPC) or silicon wafer gamma-ray space telescopes and be a powerful guide in the optimisation of the design of the multi-kilo-ton liquid argon TPCs that are under development for neutrino studies.

Tables

Table 1. Table 1: Parameters of two trackers considered in the text.

	gas argon	liquid argon	silicon
	TPC	TPC	detector
$X_0$	2351.	14.0	9.4	$\mathrm{cm}$
$l$	0.1	0.3	1.0	$\mathrm{cm}$
$\Delta x$	$l$	$l$	0.0500	$\mathrm{cm}$
$\sigma$	0.1	0.1	0.0070	$\mathrm{cm}$
$L$	30.	1000.		$\mathrm{cm}$
$N$	300	3 333	56
$p_1$	0.112	1.739	0.239	$\mathrm{MeV}/c$	eq. (38)
$p_u$	1277.	10 614 042.	71 098.	$\mathrm{MeV}/c$	eq. (40)
$p_x$	2.2	149.	3 542.	$\mathrm{MeV}/c$	eq. (42)
$p_s$	0.024	5.4	16.6	$\mathrm{MeV}/c$	eq. (61)
$p_{\ell }$	352.	2 931 742.	19 638.	$\mathrm{MeV}/c$	eq. (62)

Equations

$${\theta }_0\approx \frac{p_0}{\beta p}\sqrt{\frac{\Delta x}{X_0}},$$

$$y=a+b\times x.$$

$$I_{n+1}=D^T{\left(I_n^{-1}+B\right)}^{-1}D+M,$$

$$D=\left[\begin{array}{cc}1& l\\ 0& 1\end{array}\right],$$

$$B=\left[\begin{array}{cc}0& 0\\ 0& sl\end{array}\right],$$

$$M=\left[\begin{array}{cc}ıl& 0\\ 0& 0\end{array}\right],$$

$$\left[\begin{array}{c}{A}_{n+1}\\ B_{n+1}\end{array}\right]=\left[\begin{array}{cc}{D}^T+M{D}^{-1}B& M{D}^{-1}\\ D^{-1}B& D^{-1}\end{array}\right]\left[\begin{array}{c}{A}_n\\ B_n\end{array}\right],$$

$$I_{n+1}=A_{n+1}{B}_{n+1}^{-1}.$$

$$\Phi \equiv \left[\begin{array}{cc}{D}^T+M{D}^{-1}B& M{D}^{-1}\\ D^{-1}B& D^{-1}\end{array}\right],$$

$$\left[\begin{array}{c}{A}_n\\ B_n\end{array}\right]={\Phi }^n\left[\begin{array}{c}{A}_0\\ B_0\end{array}\right].$$

$$\frac{1}{{\beta }^n}\left[\begin{array}{c}{A}_n\\ B_n\end{array}\right]=\frac{1}{{\beta }^n}{\Phi }^n\left[\begin{array}{c}{A}_0\\ B_0\end{array}\right]$$

$$\left[\begin{array}{c}{\tilde{A}}_n\\ {\tilde{B}}_n\end{array}\right]={\tilde{\Phi }}^n\left[\begin{array}{c}{\tilde{A}}_0\\ {\tilde{B}}_0\end{array}\right],$$

$${\Phi }^{T}J\Phi =J,$$

$$J=\left[\begin{array}{ccccc}0& 0& 1& 0& \\ 0& 0& 0& 1& \\ -1& 0& 0& 0& \\ 0& -1& 0& 0& \end{array}\right],$$

$$\left\{\alpha ,{\alpha }^{\ast },\frac{1}{\alpha },\frac{1}{{\alpha }^{\ast }}\right\},$$

$$\alpha (x)=\frac{1}{2}j{x}^2+\frac{1}{2}{\left(-x^4+4j{x}^2\right)}^{\frac{1}{2}}+1,$$

$$\lambda \equiv \frac{1}{\sqrt[4]{ıs}}.$$

$$u\equiv \frac{L}{\lambda },$$

$$I=D^T{\left(I^{-1}+B\right)}^{-1}D+M.$$

$$V=\left(\begin{array}[]{cc}\frac{4l^{3}s}{x^{3}\left(2x+\sqrt{x^{2}+4j}-\sqrt{x^{2}-4j}\right)}&\frac{sl^{2}\left(\sqrt{x^{2}+4j}+\sqrt{x^{2}-4j}\right)}{x^{2}\left(\sqrt{-x^{2}-4j}-\sqrt{4j-x^{2}}-2jx\right)}\\ -\frac{4l^{2}s}{x^{2}\left(\sqrt{-x^{2}-4j}-jx\right)\left(jx+\sqrt{4j-x^{2}}\right)}&\frac{2ls\left(\sqrt{x^{2}+4j}+\sqrt{x^{2}-4j}\right)}{\left(x+\sqrt{x^{2}+4j}\right)\left(jx+\sqrt{4j-x^{2}}\right)}\\ \end{array}\right).$$

$$V_{aa}=\frac{\sqrt{2}}{ı\lambda }\left[1-\frac{x}{\sqrt{2}}+\frac{3}{8}{x}^2-\frac{\sqrt{2}}{8}{x}^3+\frac{9}{128}{x}^4+O[x^5]\right],$$

$$V_{bb}=\frac{\sqrt{2}}{ı{\lambda }^3}\left[1-\frac{x}{\sqrt{2}}+\frac{1}{8}{x}^2+\frac{x^3}{128}-\frac{1}{1024}{x}^4+O[x^5]\right].$$

$$V_{aa}={\sigma }^2,V_{bb}=2{\left(\frac{\sigma }{l}\right)}^2:$$

$$V_{aa}$$

$$V_{bb}$$

$$\dot{I}(L)=D^{\prime T}I(L)+I(L)D^{\prime }-I(L)B^{\prime }I(L)+M^{\prime },$$

$$D^{\prime }=\left[\begin{array}{cc}0& 1\\ 0& 0\end{array}\right],B^{\prime }=\left[\begin{array}{cc}0& 0\\ 0& s\end{array}\right],M^{\prime }=\left[\begin{array}{cc}ı& 0\\ 0& 0\end{array}\right].$$

$$D^{\prime T}I(L)+I(L)D^{\prime }-I(L)B^{\prime }I(L)+M^{\prime }=0,$$

$$I^0$$

$$I^1$$

$$I^{k+1}$$

$$V_{aa}=\frac{4}{ı\lambda u}\left[1+\frac{u^4}{416}-\frac{127u^8}{15891876000}+O(u^{12})\right]$$

Full text

This work could pave the way to the development of autonomous high-performance gas time-projection chambers (TPC) or silicon wafer $\gamma$-ray space telescopes and be a powerful guide in the optimisation of the design of the multi-kilo-ton liquid argon TPCs that are under development for neutrino studies.

††journal: Nucl. Instrum. Meth. A

1 Introduction

1.1 $\gamma$-ray astronomy

A huge effort is in progress to design $\gamma$-ray telescopes able to bridge the sensitivity gap that extends between the upper end of the high-sensitivity energy range of past and present X-ray telescopes and the lower end of the high-sensitivity energy range of the Fermi-LAT telescope, that is, approximately $0.1-100\,\mathrm{M}\mathrm{eV}$.

On the low-energy side of the gap, tracking of the electron issued from the first Compton scattering of an incident photon enables a major improvement of the precision of the reconstruction of the direction of the incident photon ([1] and references therein) that induces an impressive improvement of the true-photon-background rejection and therefore of the point-like-source sensitivity. A serious limitation of that ETCC (electron tracking Compton camera) scheme arises though, as the effective area undergoes a sharp drop for photon energies above $0.5\,\mathrm{M}\mathrm{eV}$, due to the fact that the recoil electron can exit on the side and escape energy measurement [1]: electron momentum measurement inside the time projection chamber (TPC) itself is highly desirable.

On the high-energy side of the gap, novel approaches improve the sensitivity by improving the single-photon angular resolution by using converters having a lower-$Z$ than that of the tungsten plates of the EGRET / Fermi-LAT series. Using a series of silicon wafer active targets placed at a distance of each other, at the same time the material in which the photon converts and in which the tracks are tracked, enables an improvement of $\approx$ a factor of three in the angular resolution at $100\,\mathrm{M}\mathrm{eV}$ with respect to the LAT [2, 3, 4, 5, 6, 7, 8, 9] at the cost of a lower average active target density. Similar values of the angular resolution are achieved using a high-spatial-resolution, homogeneous, high-density material such as an emulsion [10].

If the trend to lower densities is pushed to the use of a gaseous detector, the angular resolution with respect to the LAT can increase up to a factor of ten at $100\,\mathrm{M}\mathrm{eV}$ [11] and the single-track angular resolution is so good that the azimuthal angle of the $e^{+}e^{-}$ pair can be measured with a good enough precision to enable the measurement of the linear polarization fraction of the incident radiation [12, 13, 14]. Gas detectors enable the detection of low-energy photons close to the pair-creation threshold where most of the statistics lie for cosmic sources (Fig. 1 of [15]), something which is critical for polarimetry.

Astrophysicists also need to measure the energy of incoming photons and therefore the momentum of the conversion electron(s). This can be achieved using a number of techniques.

In a calorimeter, the total energy of the photon is absorbed and measured.

• 2.

In a magnetic spectrometer, the trajectory of a particle with electric charge $q$ and momentum $p$ in a magnetic field $B$ is curved with a curvature radius $\rho=p/(qB)$: from a measurement of $\rho$, one obtains a measurement of $p$ and in the end of the photon energy $E$.

• 3.

In a transition radiation detector (TRD), the energy of the radiation emitted in the forward direction by a charged particle at the interface between two media that have different refraction indices is proportional to the Lorentz factor $\gamma$ of the particle, enabling a direct measurement. The low number of emitted photons per track per interface has lead to the development of multi-foil systems that suffer destructive interference at high energies. Appropriate configurations have showed saturation values larger than $\gamma\approx 10^{4}$, which implies that a measurement can be done up to a photon energy of $\approx 10\,\mathrm{G}\mathrm{eV}$ [16].

The low-density active targets that have been considered above can provide a large effective area telescope only with a large volume (${\cal O}(\mathrm{m}^{3})$) and therefore the mass of the additional device used for energy measurement is a serious issue onboard a space mission. In this document we first address the performance of the track momentum measurement from measurements of the angular deflections of charged tracks due to multiple scattering during the propagation in the tracker itself.

1.2 Large noble-liquid TPCs for neutrino physics

Neutrino oscillation is a well established phenomenon and several experiments are being prepared with the goals:

To test the occurrence of CP violation in the neutral lepton sector, i.e. to measure the only free complex phase $\delta$ of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix with enough precision to determine its non-compatibility with zero,

• 2.

To determine unambiguously the 3 neutrino mass ordering, i.e. to solve the sign ambiguity of the square mass difference $\Delta m^{2}_{31}$.

Not only the (vacuum propagation) phase term that involves $\delta$ changes sign upon $\nu\leftrightarrow\overline{\nu}$ exchange, but the term that describes the interaction with matter changes sign too as our Earth contains much more electrons than positrons. “In the few-GeV energy range, the asymmetry from the matter effect increases with baseline as the neutrinos pass through more matter, therefore an experiment with a longer baseline [is] more sensitive to the neutrino mass hierarchy. For baselines longer than $\approx$ 1200 km, the degeneracy between the asymmetries from matter and CP-violation effects can be resolved” [17]. Large distances imply low fluxes, that is, huge detectors and, to match the $\sin{(\Delta mL/4E_{\nu})}$ oscillation function, high-energy neutrinos. So we should be prepared to measure the momentum of high-momentum muons in huge non-magnetised detectors such as liquid argon (lAr) TPCs.

The DUNE experiment expects to be able to measure muon momenta with a relative precision of $\approx 18\%$ [18], based on a past ICARUS work [19]. They “anticipate that the resolution will deteriorate for higher-energy muons because they scatter less”, though. Given the $\mbox{d}E/\mbox{d}x$ of $0.2\,\mathrm{G}\mathrm{eV}/\mathrm{m}$ of minimum ionising particles in lAr, a typical $6\,\mathrm{G}\mathrm{eV}/c$ muon produces a long track: it should be interesting to study to what extent an optimal analysis of the thousands of measurements per track, at their $\approx 3\,\mathrm{m}\mathrm{m}$ sampling pitch, can do better.

1.3 Track momentum measurement from multiple scattering

The measurement of track momentum using multiple scattering was pioneered by Molière [20] and has been used since, in particular in the context of emulsion detectors (recent accounts can be found in [21, 22]).

In a practical detector consisting of $N$ detection layers, the precision of the deflection measurement and therefore of the momentum measurement is affected by the precision, $\sigma$, of the measurement of the position of the track when crossing each layer: the combined square deflection angle summed up over the whole track length therefore includes contributions from both the scattering angle and the detector precision. Bernard has optimized the longitudinal “cell” length over which each deflection angle is measured [11] and obtains a value of the relative momentum precision $\sigma_{p}/p$ that scales as $p^{1/3}$, but the fact that the track position precision can improve when the cell length is extended and several measurements are combined was not taken into account in [11]. **In the present document we study an optimal method of momentum measurement with a tracker that has a finite (non zero) precision. **

In Section 2 we revisit optimal tracking methods in a context where the momentum of the particle is known. This allows us to present concepts and notations that are used later in the paper. We also extend the results published in the past by the use of more powerful methods.

The optimal precision of track measurements obtained in Sec. 2 can be obtained by performing the fit with a Kalman filter (KF), a tool that was imported in our field by Frühwirth [23]. In section 3 we give a brief summary of Kalman filter tracking in a Bayesian formalism. In magnetic spectrometers, the particle momentum takes part both in the particle state vector through the curvature of the trajectory and in the magnitude of multiple scattering. The precision of the magnetic measurement is most often so good that the momentum can rightfully be considered as being perfectly known in the expression of the multiple scattering. In our case of a zero magnetic field, it is not the case. A Kalman filter is the optimal estimate for linear system models with additive white noise, such is the case for multiple scattering (process noise) and detector precision (measurement noise), but at the condition that the optimal Kalman gain be used in the expression, that is, that the track momentum be known. In section 4, we use the Bayesian method developed by Matisko and Havlena [24] to obtain an optimal estimator of the process noise covariance, and therefore of the track momentum 111Attempts of estimation of track momenta based on the use of a Kalman filter have been performed in the past, with little success. The un-validated un-characterized study of Ref. [25], for example, shows a poor relative resolution of $\sigma_{p}/p=30-40\%$ and that does not vary with the true particle momentum between $50\,\mathrm{M}\mathrm{eV}/c$ and $2\,\mathrm{G}\mathrm{eV}/c$, which is a bad symptom.. We implement this method and characterize its performance on Monte Carlo (MC) simulated tracks. We check that the momentum measurement is unbiased within uncertainties. We obtain a heuristic analytical expression of the relative momentum uncertainty.

Numerical examples are given for a homogeneous gas detector such as an argon TPC and for a silicon-wafer detector:

TPC gas, argon, 5 bar, $\sigma=l=0.1\,\mathrm{c}\mathrm{m}$, $L=30\,\mathrm{c}\mathrm{m}$ [12];

• 2.

Silicon detector $N=56$, $\Delta x=500\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$-thick wafers spaced by $l=1\,\mathrm{c}\mathrm{m}$, with a single point precision of $\sigma=70\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ [8].

In this work a number of approximations are done: only the Gaussian core of the multiple-scattering angle distribution is considered and the non-Gaussian tails due to large-angle single scatters are neglected. The small logarithmic correction term in the expression of the RMS multiple scattering angle, $\theta_{0}$, is neglected

[TABLE]

where $p_{0}=13.6\,\mathrm{M}\mathrm{eV}/c$ is the “multiple-scattering constant”, $\Delta x$ is the matter thickness through which the particle propagates and $X_{0}$ is its radiation length (Eqs. (33.14), (33.15), (33.17) of [26]). In the case of a homogeneous detector, the thickness of the scatterer is equal to the length of the longitudinal sampling, $l=\Delta x$. We assume relativistic particles ($\beta\approx 1$) without loss of generality. Only the first-order term (angle deflection) of multiple scattering is taken into account which is legitimate for the thin detectors considered here; the 2nd-order transverse displacement (eq. (33.19) of [26]) is neglected. Continuous ($\mbox{d}E/\mbox{d}x$) and discrete (BremsStrahlung radiation) energy losses are also neglected. In TPCs in which the signal is sampled, most often the electronics applies a shaping of the pulse before digitisation, that creates a short scale longitudinal correlation between successive measurements that we neglect too. Also the limitations of pattern recognition, that is, in the case of $\gamma$-ray telescopes, of the assignment of each hit to one of two close tracks, are not addressed.

Note that in the two main parts of this work (section 2 and sections 3-4) we have made our best to follow the notations of Refs. [27, 28] and of [24], respectively, and that they turn out to differ to some extent.

2 Tracking

An optimal tracking makes use of the full $N\times N$ covariance matrix of the $N$ measurements, including multiple scattering (correlation terms). This is most often impractical in modern trackers that provide a huge number of measurements for each track. The first successful attempt to perform a recursive determination of the covariance matrix was achieved by Billoir [27]. He considered the paraxial propagation of a charged track along the $x$ axis inside a magnetic field oriented along $z$: close to the particle origin, the trajectory is a straight line in the $(x,z)$ plane, and a parabola osculatrix to the true circle in the $(x,y)$ plane. As we examine here the case of a magnetic-field-free tracker, the propagation (in the $(x,z)$ and in the $(x,y)$ planes) is approximated by straight lines (already using Innes notations [28] but assuming $B=0$):

[TABLE]

Astronomers obviously have a special interest in the slope $b$, that is, in the paraxial direction of the track at the conversion vertex.

The $(a,b)$ correlation matrix is named $V$ and the information matrix, $I\equiv V^{-1}$. Billoir develops a recursive method in which the fit propagates along the track, adding the information gain (measurement) and loss (scattering) at each layer. He obtains the information matrix at layer $n+1$ from the information matrix at layer $n$ [27, 28]:

[TABLE]

where $D$ is the drift matrix that propagates the track from layer $n$ to layer $n+1$,

[TABLE]

$l$ is the layer spacing. $B$ is the scattering matrix,

[TABLE]

where $s\equiv\left(\displaystyle\frac{p_{0}}{p}\right)^{2}\displaystyle\frac{\Delta x}{lX_{0}}$ is the average multiple-scattering angle variance per unit track length, $\theta_{0}^{2}=s\times l$. $M$ is the measurement matrix,

[TABLE]

where $\imath\equiv\displaystyle\frac{N+5}{L\sigma^{2}}\approx\displaystyle\frac{1}{l\sigma^{2}}$ is the information density per unit track length, $L=N\times l$ is the full detector thickness.

Billoir considers the two particular cases of “scatters at one point” (detector layers separated by an empty space) that we name here a segmented detector and “uniformly distributed scattering” (homogeneous detector) [27]. These concepts are defined more precisely below, following Innes [28].

2.1 Segmented detector

Expressing $I_{n}$ as $I_{n}=A_{n}B_{n}^{-1}$ ([29], page 149) we obtain

[TABLE]

and

[TABLE]

Noting

[TABLE]

we obtain

[TABLE]

$A_{n}$ and $B_{n}$ are obtained from the eigenvalues of $\Phi$. The covariance matrix becomes $V_{n}=B_{n}A_{n}^{-1}$. We initialise the recurrence with $A_{0}=I_{0}=0$, $B_{0}=1$. If $\begin{bmatrix}A_{n}\\ B_{n}\end{bmatrix}$ is a solution of eq. (8), then for any $\beta>0$,

[TABLE]

is a solution too. Noting $\tilde{A}_{n}=\displaystyle\frac{A_{n}}{\beta^{n}}$, $\tilde{B}_{n}=\displaystyle\frac{B_{n}}{\beta^{n}}$, $\tilde{\Phi}=\displaystyle\frac{\Phi}{\beta}$, we obtain

[TABLE]

with $V_{n}=\tilde{B}_{n}\tilde{A}_{n}^{-1}$. $\Phi$ is found to satisfy

[TABLE]

with

[TABLE]

from which $\Phi$ is a symplectic matrix. General theorems enable a classification of $\Phi$ eigenvalues into two “invert” and “conjugate” blocks, respectively (eq. (10) of [30])

[TABLE]

where “$*$” denotes complex conjugation. We choose $\alpha$ to have a norm larger than unity, $|\alpha|>1$. We obtain [31]

[TABLE]

where $j$ is the imaginary unit and $x\equiv\displaystyle\frac{l}{\lambda}$ is the detector longitudinal sampling normalized to the detector scattering length at momentum $p$ [28]:

[TABLE]

An exploration of the consequences of a variation of the initialisation of the recurrence parameters shows that the system converges to the same solution $\alpha(x)$ regardless of the values of $A_{0},B_{0}$. $B_{0}\neq 0$ is needed so that $I_{0}$ is defined. With $B_{0}=1$, $I_{0}=A_{0}=0$ simply assumes that no a priori information is known about the track.

We study the convergence of the covariance matrix while the Billoir mechanism is in progress along the track (increasing $n$) by setting $\beta=|\alpha|$, that is, $\tilde{\Phi}=\displaystyle\frac{1}{|\alpha|}\Phi$. $\tilde{\Phi}$ has two eigenvalues with modulus unity and two eigenvalues with modulus $\displaystyle\frac{1}{|\alpha|^{2}}$. The unity-modulus eigenvalues could be a major nuisance in the behaviour of $I_{n}$ as a function of $n$, but when applying the Billoir mechanism we observe that for some reason the amplitude the so-induced oscillating terms is zero. The convergence behaviour is then driven by the two other eigenvalues, that is, by terms proportional to $\displaystyle\frac{1}{|\alpha|^{2n}}$. That exponential convergence is illustrated in Fig. 1 that shows the value of the detector thickness normalized to the detector scattering length, $u$ [28],

[TABLE]

for which $\displaystyle\frac{1}{|\alpha|^{2n}}<10^{-4}$, as a function of $x$.

Note that the homogeneousness parameter $x$ and the thickness parameter $u$ have a similar dependence on track momentum $p$, as $u=x\times N$.

2.1.1 Segmented detector: Thick detector limit

The asymptotic expression at high $n$, that is, at high $u$ (thick detector) is reached after the Billoir mechanism (eq. (3)) has converged: we obtain the discrete Riccati equation:

[TABLE]

When the geometric, the multiple scattering and the measurement properties of the detector are uniform (at least piecewise) the dynamics of the particle is described by a time-invariant system and eq. (19) is referred to as the “algebraic” Riccati equation (DARE). Equation (19) has four solutions, but the fact that the asymptotically stable solution must be positive definite (Theorem 2.2 of [30]) leaves us with only one.

Segmented detector: Exact solution

We obtain [31]

[TABLE]

Even though it is not explicit from eq. (20), $V$ is found to be a real matrix, which is decent for a covariance matrix.

Segmented thick detector: Small $x$ behaviour: Homogeneous detector limit

The Taylor expansion close to $x=0$ is found to be

[TABLE]

[TABLE]

These expressions are similar222We have detected a misprint, though: the factor $-5/8$ in their expressions of $V_{aa}$ is here corrected to $3/8$. to what was found by Billoir ([27], p364, no magnetic field) and Innes ([28], eq. (8), magnetic field). The Taylor expansion is found to converge for $x\lesssim 2$ (Fig. 2).

Segmented thick detector: Large $x$ behaviour: Coarse segmentation limit

The asymptotic behaviour of the coarsely instrumented detector (high $x$) is presented in Fig. 3.

For $1/x=0$ we obtain $\imath lV_{aa}=1$ and $\imath l^{3}V_{bb}=2$, that is the obvious

[TABLE]

the scattering is so intense that the intercept (angle) measurement is based on the first (two first) layer(s), respectively. A thick coarse detector can be defined by $1/x<0.5$, that is, $l>2\lambda$ (Fig. 3). The $\displaystyle\frac{1}{x}$ Taylor expansion is:

[TABLE]

2.2 Homogeneous Detector

A homogeneous detector is described having $l$ tend to [math] while $s$ and $\imath$ are kept constants. Fig. 4 shows that for all values of $u$, the intercept and angle variances become very close to the homogeneous limit $(x=0)$ for $x\lesssim 0.2$.

From the discrete evolution equation, eq. (3), and denoting $I_{n}=I(nl)=I(L)$, we obtain

[TABLE]

where the dot denotes the derivation with respect to $L$ and with

[TABLE]

After convergence (thick detector), we obtain the continuous algebraic Riccati equation (CARE):

[TABLE]

Homogeneous Detector: Small $u$ behaviour

We first use Innes’ method to compute an approximate solution. Attempting a Taylor expansion in $u$, $I(u)=\sum I^{k}u^{k}$, we obtain:

[TABLE]

In our case ($B=0$, no curvature) we obtain:

[TABLE]

[TABLE]

For $B\neq 0$ and a fit with curvature, we obtain the results of Innes (eq. (9) of [28]).

These Taylor expansions converge for $u\lesssim 3.5$ (no curvature, Fig. 5) and $u\lesssim 7.0$ (with curvature, [28]). The thin detector ($u=0$) value of $V_{bb}$ for $B=0$ is found to be smaller than that for $B\neq 0$ by a factor of 16, as was discussed in the corrigendum of [11]: in a fit with curvature, the correlation between the curvature and the angle at the end(s) of the track degrades the angular resolution badly; this lasts until $u\approx 1$ that is $L\approx\lambda$ (plot not shown), after which all is flooded by multiple scattering anyway.

Homogeneous Detector: Large $u$ behaviour: Thick detector limit

Searching for expressions that are valid at high $u$, we follow again Innes and search a solution of the continuous equation for $V$ that is similar to eq. (26) for $I$. Here the Taylor expansion is searched in $1/u$. Searching a solution parametrized as $V(u)=V_{0}+\displaystyle\frac{1}{u}V_{1}$, we obtain

[TABLE]

These values agree with that of eqs. (21), (22) for $x=0$. The term proportional to $1/u$ that was present in the case with curvature (eq. (11) of [28]) cancels here, which is related to the exponential convergence seen on eq. (37) (see also Fig. 7).

2.2.1 Homogeneous thick detector: Exact solution

We solve the continuous algebraic Riccati equation in a way similar to the discrete case [29]: Expressing

[TABLE]

and $I(L)=X(L)Y^{-1}(L)$, and taking $I(0)=0$, we obtain:

[TABLE]

$\Phi^{\prime}$ satisfies

[TABLE]

and is a hamiltonian matrix (section (4.8) of [29], [32]), which implies that $\exp\left[\Phi^{\prime}\right]$ is a symplectic matrix and therefore that $\Re\{\mathrm{Tr}(\Phi^{\prime})\}=0$. Furthermore all eigenvalues of $\Phi^{\prime}$ are found to be non singular [31]:

[TABLE]

Solving eq. (34) we obtain [31]

[TABLE]

Even though it is not explicit from eq. (37), $V$ is found to be a real matrix, which is decent for a covariance matrix. The convergence is driven by a term proportional to $e^{-2ue^{\frac{j\pi}{4}}}$, which implies a convergence in $e^{{-l\sqrt{2}}/\lambda}$.

2.2.2 Variation of the angle variance along the track

We have considered above the optimal measurement of the track parameters at the vertex, $z=0$. Here we examine the measurement at any position along the track. A track now consists of two segments (left and right), the fit of each of which provides an estimate with its own covariance matrix.

A first combination attempt is performed on the two variables ($a$ and $b$) separately. As is obvious, if the detector is thick on both sides, the two estimates (right, left) of a track parameter (say: the angle) have the same uncertainty on the plateau and their optimal combination provides a gain in RMS precision of a factor of $\sqrt{2}$. But that neglected the fact that in the combination, the other variable should be constrained to have the same value on both sides too. With a weighted variance estimation, a further gain of a factor of $\sqrt{2}$ is obtained on the plateau, that is a total gain of a factor of $2$ with respect to individual measurements (Fig. 6) as was observed experimentally by running a KF on simulated tracks (Fig. 18 of [12]).

2.3 Optimal tracking: wrap up

We now build on the results of the previous subsections to obtain expressions of the variances in terms of the detector parameters. We do so for segmented detectors. The expressions for continuous detectors can be obtained with $\Delta x=l$.

[TABLE]

With:

small $x$ (large $p$), homogeneous detector (continuous equation),

• 2.

large $x$ (small $p$), segmented detector (discrete equation)

and

small $u$ (large $p$), thin detector,

• 2.

large $u$ (small $p$), thick detector.

The variances are found to be asymptotically:

[TABLE]

where $p_{1}$ is a momentum that characterises the tracking angular-resolution properties of a detector affected by multiple scattering [12]

[TABLE]

The two $V_{aa}$ asymptotes cross for $u=u_{c,a}=2\sqrt{2}\approx 2.83$;

• 2.

The two $V_{bb}$ asymptotes cross for $u=u_{c,b}=(12/\sqrt{2})^{1/3}\approx 2.04$.

This, for a given detector, takes place for a value of the momentum $p_{u}$ for which

[TABLE]

that is,

[TABLE]

from which

[TABLE]

In short, a homogeneous detector is a thick detector, $u>u_{c}$, at low momentum, $p<p_{u}$ and a thin detector at higher momentum. In Table 1, we use $u_{c}=2.5$ to compute the value of $p_{u}$.

In the same way, a detector is a homogeneous detector, $x<x_{c}$ at high momentum, $p>p_{x}$, with

[TABLE]

and $x_{c}=0.2$ (Fig. 4). And similarly:

[TABLE]

Argon gas TPC. We see that $p_{u}>p>p_{x}$ for most of the [1 MeV - 1 GeV] momentum range that is the primary target of the high-performance $\gamma$-ray telescopes mentioned above: the telescope is both a homogeneous and a thick detector.

• 2.

Silicon detector. Here the telescope is a segmented and a thick detector for most of the momentum range.

Note that the equality $p_{u}=p_{x}$ holds for $N=u_{c}/x_{c}$, that is, $N=2.5/0.2=12.5$, so for most conceivable detectors, $p_{u}>p_{x}$, that is,

if $p>p_{u}$ then $p>p_{x}$, if a detector is thin for a given track, then it’s also homogeneous;

• 2.

if $p<p_{x}$ then $p<p_{u}$, if a detector is segmented for a given track, then it’s also thick.

3 Kalman filter

A Kalman filter is an estimator of the state of a linear dynamic system perturbed by Gaussian white noise using measurements that are linear functions of the system state and corrupted by additive Gaussian white noise [29]. The paraxial propagation of a high-momentum particle inside a detector is affected by angular deflections due to scattering on the charged particles (electrons, nuclei) present in the detector matter. Deflections undergone at different locations on the track are uncorrelated, “white process noise”, and are approximated to have a Gaussian distribution under the multiple-scattering approximation. Transverse position measurements are performed at several locations along the track. They are affected by an uncertainty that does not correlate from layer to layer, “white measurement noise”, and that most often can be approximated by a Gaussian distribution. Angular deflections and measurement uncertainty are not correlated. When the system is non-linear, such as for the propagation in a magnetic field, it is linearized locally, “extended Kalman filter” ([33] and references therein). In the case of most particle detectors, the geometric, the multiple scattering and the measurement properties of the detector are uniform (at least piecewise) so the dynamics of the particle is described by a time-invariant system.

Since the founding work by Frühwirth [23], KF tracking has been used largely in high-energy physics. We present here a short description of the elements that are used in the next section, in a Bayesian formulation. Denoting $\{z^{0}_{n}\}$ and $\{z^{m}_{n}\}$ the true and the measured positions of a particle at layer $n$, respectively, and $x_{n}$ the corresponding state vector, $\hat{x}_{n}=\mathbb{E}(x_{n}|z_{0}^{m},\cdots,z_{n}^{m})$ is the estimator of $x_{n}$ conditioned to $\{z_{0}^{m},\cdots,z_{n}^{m}\}$ and $x_{n}^{n-1}=\mathbb{E}(x_{n}|z_{0}^{m},\cdots,z_{n-1}^{m})$ is the prediction of $x_{n}$ given $\{z_{0}^{m},\cdots,z_{n-1}^{m}\}$. $x_{n}$ is obtained from $x_{n-1}$

[TABLE]

$u_{n}$ is the Gaussian-distributed deflection angle with variance $sl$. The covariance matrix of the state vectors is $P_{n}=\mathbb{E}((\hat{x}_{n}-x_{n})(\hat{x}_{n}-x_{n})^{T}|z_{0},\cdots,z_{n}^{m})$. The optimal estimator of $x_{n}$ is obtained from $\hat{x}_{n-1}$ and from the measurements $\{z_{0}^{m},\cdots,z_{n-1}^{m}\}$,

[TABLE]

with

[TABLE]

where $H=\begin{bmatrix}0&1\end{bmatrix}$ is the measurement matrix and $v_{n}$ is the measurement uncertainty which is Gaussian-distributed with variance $\sigma^{2}=\displaystyle\frac{1}{\imath l}$. The innovations are the difference between measurement and prediction,

[TABLE]

and their variance is

[TABLE]

The gain matrix of the filter is

[TABLE]

For the optimal value of the gain that minimizes the variance of the innovations, we obtain [29]

[TABLE]

Noting $Z^{n}=z_{0}^{m},\cdots,z_{n}^{m}$ the set of measurements up to layer $n$ and $p$ the probability density,

[TABLE]

As $z_{i}^{m}|Z^{i-1}$ is Gaussian distributed ${\mathcal{N}}(z_{i},S_{i})$:

[TABLE]

We have implemented such a KF tracking software. Figure 8 shows a couple of sanity-check validation plots, the RMS of the position residues (left) and of the innovation residues (right) as a function of the longitudinal position in the tracker, for a sample of $10^{6}$ $50\,\mathrm{M}\mathrm{eV}/c$ tracks in a silicon detector, compared to the RMS computed from their variance, $P_{n}$ and $S_{n}$, respectively.

4 Momentum Measurement

A Kalman filter is the optimal linear estimator of the state vector of a dynamical system at the condition that the model be an accurate description of the dynamics of the system and that the process and measurement noise covariance matrices be known, that is here, that the track momentum be known. The estimation of the noise covariance matrices of a dynamic system was pioneered by Mehra [34, 35] who studied and compared several methods:

A Bayesian method that is the root of that we use in this work.

• 2.

Maximum likelihood methods, if necessary of both the state vector of the system and the noise matrices at the same time.

Bayesian and maximum likelihood methods were deemed to be too CPU consuming for the time.

Covariance-matching techniques, making the innovation residuals consistent with their theoretical covariances; these methods were shown later to give biased estimates of the covariance matrices.

• 2.

Correlation methods, in particular based on the observation that when the KF gain $K$ is optimal, the innovations of the filtering process are white and Gaussian.

Mehra showed that the optimal gain $K$ can be determined uniquely, after which many efforts and publications have then been spent in determining the convergence of these methods and to which values of the process and measurement noise matrices they were, eventually, converging. In 2006, Odelson et al. [36] re-examined Mehra’s work, showed that the definite positiveness of the matrices was not assured; based on the fact that the autocorrelation of the innovation sequence is linearly dependent on the noise covariances [37], they developed an autocovariance least-square (ALS) method that provides unbiased estimates of the noise matrices and that includes a mechanism that enforces definite positiveness [36].

Kalman filters have already been used for momentum measurement in non-magnetic particle physics detectors in the past ([38, 39]. The trick is to augment the state vector with the parameter vector to ($x,y,\mbox{d}x/\mbox{d}z,\mbox{d}y/\mbox{d}z,1/p$) so that the KF performs their estimation simultaneously. However, this augmentation approach has been originally intended for estimating parameters in deterministic part of the model and its straightforward application for noise covariance matrices does not result in appropriate estimates ([40] and references therein).

In the case of charged particle tracking in a magnetic-field-free detector, the augmentation method was found to provide unbiased results though [38, 39], most likely as the tracking part and the deflection part of the filter behave as two separate filters, “only” linked by the joint uses of the track momentum, one for the process noise matrix, the other as part of the state vector. Also it enables the optimal treatment of energy loss and therefore it provides an improvement of about a factor of two with respect to the Molière method [39]. If would be interesting to examine to what extent they are efficient or even whether they are optimal. Note that in that scheme, the track has to be segmented to measure the track angle on each segment (in $\approx 19\,\mathrm{c}\mathrm{m}$ long segments that contain $\approx 57$ hits on average for[39], from which they obtain a relative resolution of 16 % on a sample of $4\,\mathrm{m}$ tracks with momenta ranging from 0.5 to $4.5\,\mathrm{G}\mathrm{eV}/c$). (See also [41]).

4.1 Single track momentum measurement: Bayesian method

Following Matisko and Havlena [24] we obtain 333 We assume that the detector spatial resolution $\sigma$ is known either from calibration on beam or from the analysis of high momentum tracks., from the measurements, the most probable value of $s$, and we extract from it an optimal estimator $\hat{p}$ of $p$. For an event $A$, defining $p_{n}(A)\equiv p(A|Z^{n})$, we have

[TABLE]

We name $s$-filter a KF with a gain matrix computed with a given value of $s$. We remember (eq. (49)) that the $s$-filter innovation probability density function (pdf) $\beta_{n}=p_{n-1}(z_{n}^{m}|s)$, is a normal, $\beta_{n}=\mathcal{N}(\nu_{n}(s),0,S_{n}(s))$, where $\nu_{n}(s)$ and $S_{n}(s)$ are computed during the filtering process. The $1/p_{n-1}(z_{n}^{m})$ factor does not vary with $s$ and is therefore neglected. We obtain

[TABLE]

The distribution of $p(s)$ for one simulated $50\,\mathrm{M}\mathrm{eV}/c$ track is shown in Fig. 9. The track momentum is then obtained from the value of $s$ that maximizes $p_{n}(s)$:

[TABLE]

From the full width half maximum (FWHM) of $p(s)$ we calculate $\mathrm{RMS}_{s}/s$ and $\mathrm{RMS}_{p}/p=(\mathrm{RMS}_{s}/s)/2$. The average value of $\mathrm{RMS}_{p}/p$ is found to be much smaller than the relative RMS $\sigma_{p}/p$ of the momentum measurements performed on a sample of tracks and shown in Figs. 10 and 11. Our interpretation is that the range of $s$ values compatible within uncertainties with the deflection sample of a given track and the fluctuation of the deflection sample from track to track are separate quantities. In addition, we observe that the single-track $p(s)$ width and the measured momentum for that track are weakly correlated, so the former would add little information to the measurement of the latter.

We were not able to obtain an analytical expression for the relative precision of the momentum measurement. Instead we performed a parametric variation study for a silicon detector with

$l=0.5,1.0,2.0\,\mathrm{c}\mathrm{m}$,

• 2.

$X_{0}=4.685,9.37,18.74\,\mathrm{c}\mathrm{m}$,

• 3.

$N=23,46,56,92$,

• 4.

$\sigma^{2}=2.5,5.0,10.0\times 10^{-5}\,\mathrm{c}\mathrm{m}^{2}$,

• 5.

$p=1\cdots 2048\,\mathrm{M}\mathrm{eV}/c$

and with $\Delta x=500\,\mathrm{\SIUnitSymbolMicro}\mathrm{m}$. A good representation of these data is obtained with the following expression:

[TABLE]

from which we obtain the obvious low-momentum asymptote

[TABLE]

• 2.

and the high-momentum asymptote

[TABLE]

Of particular interest is the momentum, $p_{s}$, above which $\sigma_{p}/p$ starts to depart from the low momentum asymptote,

[TABLE]

We define also the momentum, $p_{l}$, above which $\sigma_{p}/p$ is larger than unity, which means that the measurement becomes meaningless:

[TABLE]

The only thing that can be said then is that that track is a straight track within uncertainties, that is, with inverse momentum $1/p$ compatible with zero. These two momenta are characteristics of the ability to measure track momenta with a given detector and are related to each other,

[TABLE]

Finally, we obtain a simpler expression of the relative momentum resolution,

[TABLE]

The target relative precision of the DUNE project of 18 % is within reach for $10\,\mathrm{m}$ tracks up to a momentum of $17.1\,\mathrm{G}\mathrm{eV}/c$ with detector parameter values from Table 1 (Fig. 12).

4.2 Comparison with the cell-optimization result

For the continuous detector, $\Delta x=l$, eq. (60) becomes

[TABLE]

that we can compare to the cell-optimization expression (eq. (12) of [11]):

[TABLE]

with $C\equiv 5^{1/6}+5^{-5/6}\approx 1.57$. We see that the precisions are commensurate at small $N$ and that the present approach becomes more precise at larger $N$, within the high-momentum approximation, $p\gg p_{s}$.

4.3 Cramér-Rao Bounds

The Cramér-Rao bound is a lower bound on the variance of an estimator. If the variance of the estimator reaches the Cramér-Rao bound, it can be stated that the estimate is optimal. The Cramér-Rao criterion for an estimator $\hat{\theta}$ of a parameter $\theta$ obtained from measurements $Z^{N}$ is [42]:

[TABLE]

where $I$ is the Fischer information. If $\hat{\theta}$ is an unbiased estimator of $\theta$, then

[TABLE]

Following the recursive method of [42] we obtain

[TABLE]

that is, finally, the obvious

[TABLE]

No major insight obtained with the Cramér-Rao Bounds then.

4.4 Smoothing and Momentum Measurement

In this section 4 we have obtained an optimal estimator of a charged particle momentum based on the analysis of the filtering innovations of KFs with variable $s$ parameters. After filtering, a KF provides an optimal estimate of the state vector parameters (transverse position and angle) of the track at the end of the track. An optimal estimate all along the track can be obtained by an additional, backward, pass named smoothing [23]. One might consider a scheme for momentum measurement based on the smoothing innovations rather than on the filtering innovations in the hope that the performance would be even better.

Actually smoothing is equivalent to an optimal linear combination of two independent filterings performed in the direct and in the backward directions, respectively [43]. We have compared the values of the estimators of the particle momentum obtained in the two directions and found them to be equal for each track. Therefore no further improvement is to be expected with such a combination, nor with a measurement based on smoothing innovations.

5 Conclusion

We first reconsider tracking with multiple scattering and detector resolution in magnetic-field-free detectors under the assumption that the track momentum is known, using optimal methods. This is done under a number of approximations, including Gaussian-distributed multiple-scattering deflections and the absence of energy loss during propagation. The information matrix is updated recursively while the track proceeds through the detector: after this mechanism has converged, the information matrix is found to be a solution of a Riccati equation which is not surprising as this optimal estimation can be performed with a Kalman filter.

For segmented detectors (discrete Riccati equation) and homogeneous detectors (continuous Riccati equation), we obtain exact expressions of the variances of the intercept and of the angle from the solution of that equation (eqs. (20) and (37), respectively). We compare their Taylor expansions with expressions published in the past. Convergence (thick detector) takes place after a detector thickness $L\gtrsim 2.5\lambda$ for homogeneous detectors ($l\lesssim 0.2\lambda$), and for somewhat larger values of $L$ for segmented detectors (Fig. 1). $\lambda$ is the detector scattering length for track momentum $p$.

For a given track momentum, a homogeneous detector is defined as the small longitudinal sampling limit, $l\to 0$, $\imath$ and $s$ being kept constant. In practice a limit of $l/\lambda\lesssim 0.2$ is found (Fig. 4). In contrast with magnetic spectrometers, for which the large $L/\lambda$ Taylor expansion contains $1/(L/\lambda)^{n}$ terms, here ($\vec{B}=\vec{0}$), the expansion contains only exponential terms and convergence is therefore much faster (Fig. 7). For coarse segmented detectors for which $l/\lambda\gtrsim 2$, e.g. for $p\lesssim 35\,\mathrm{M}\mathrm{eV}/c$ for eASTROGAM [8] or AMEGO [7] a KF becomes useless as the angular resolution is determined mainly by the measurements in the two first wafers (Fig. 3).

We then obtain an optimal estimator of the track momentum by a Bayesian analysis of the filtering innovations of a series of Kalman filters applied to the track. A numerical characterisation of the method shows that for a given detector the method is reliable up to some limit momentum $p_{\ell}$ above which the relative precision $\sigma_{p}/p$ becomes larger than unity. For lower momentum tracks, $p\ll p_{\ell}$, the momentum estimation is found to be unbiased. We perform a parametric study of the estimator from which we extract a heuristic analytical description of the relative uncertainty of the momentum measurement (eq. (58)).

6 Acknowledgment

It is a pleasure to acknowledge the support of the French National Research Agency (ANR-13-BS05-0002).

[TABLE]

[TABLE]