Global Decay Chain Vertex Fitting at B-Factories

TL;DR

This paper introduces a Kalman Filter-based vertex fitting method for B-factories that improves the reconstruction of complex decay chains, including neutral particles, enhancing analysis precision and background suppression.

Contribution

It presents a novel global decay chain vertex fitting technique using Kalman Filter, capable of handling neutral particles and providing comprehensive parameter estimates.

Findings

Enhanced vertex parameter extraction accuracy.

Effective suppression of combinatorial background.

Improved resolution for decay chain measurements.

Abstract

We present a particle vertex fitting method designed for B factories. The presented method uses a Kalman Filter to solve a least squares estimate to globally fit decay chains, as opposed to traditional methods that fit each vertex at a time. It allows for the extraction of particle momenta, energies, vertex positions and flight lengths, as well as the uncertainty estimates of these quantities. Furthermore, it allows for the precise extraction of vertex parameters in complex decay chains containing neutral final state particles, such as ${\gamma}$ or $K^0_L$ , which cannot properly be tracked due to limited spatial resolution of longitudinally segmented single-layer crystal calorimeters like the Belle II ECL. The presented technique can be used to suppress combinatorial background and improve resolutions on measured parameters. We present studies using Monte Carlo simulations of collisions in the Belle II experiment, where modes with neutrals are crucial to the physics analysis program.

Figures

Tables

Table 1. Table 1 : Definitions of the Belle II helix parametrisation and their dependencies.

$d_0$	The distance of closest approach to the $z$ axis (POCA) signed with the $z$ component of the angular momentum w.r.t. to the origin.
${\varphi }_0$	Azimuthal angle of the momentum at the POCA.
$\omega$	Scaled inverse of the track momentum.
$z_0$	The pivotal point is the perigee that is the POCA.
$\tan \lambda$	The angle of the momentum at the POCA with respect to the $x-y$ plane.
$\varphi =\mathrm{atan2}(p_y,p_x)$	Angle of the $x-y$ plane to the helix.
${\Delta }_{\parallel }=-x\cdot \cos \varphi -y\cdot \sin \varphi$	Quantities used to transport the coordinate system, such that $\{x,y,z\}$ points to the perigee.
${\Delta }_{⟂}=-y\cdot \cos \varphi +x\cdot \sin \varphi$
$A=2\cdot {\Delta }_{⟂}+\omega \cdot \left({\Delta }_{⟂}^2+{\Delta }_{\parallel }^2\right)$
$U=\sqrt{1+\omega \cdot A}$
$\lambda$	Angle between the $x-y$ plane and the helix.
$l=\mathrm{atan2}(\omega \cdot {\Delta }_{\parallel },1+\omega \cdot {\Delta }_{⟂})$	Arc length from the perigee to a point.
$p_x,p_y,p_z$	Momenta along the $x,y,z$ axes.
$q$	Charge of the particle.
$B_z$	Magnetic field strength in $z$-direction.
$a=B_z/c$	Magnetic field strength in the $z$ direction divided by the speed of light.
$p_t=\sqrt{p_x^2+p_y^2}$	Transverse momentum.
$p_{t0}=\sqrt{p_{x0}^2+p_{y0}^2}$	Pseudo transverse momentum.
$p_{x0}=p_x-a\cdot q\cdot y$	Pseudo momentum along the $x$ direction.
$p_{y0}=p_y+a\cdot q\cdot x$	Pseudo momentum along the $y$ direction.
$r^2=x^2+y^2$	Radius squared.
$\beta =1+\frac{p_{t0}}{p_t}$	Quantity used for reading convenience.

Table 2. Table 2 : Selection criteria used in the analyses in this section. All particles use the same criteria if appearing in the decay chains. The invariant mass is obtained by summing the particle’s daughter four-momenta before performing the fit. The particle ID for K + / π + superscript 𝐾 superscript 𝜋 K^{+}/\pi^{+} is defined as the likelihood ratio ℒ ​ ℛ ​ ( K + ​ ( π + ) ) = ℒ ​ ( K + ​ ( π + ) ) / ( ℒ ​ ( π + ) + ℒ ​ ( K + ) ) ℒ ℛ superscript 𝐾 superscript 𝜋 ℒ superscript 𝐾 superscript 𝜋 ℒ superscript 𝜋 ℒ superscript 𝐾 \mathcal{LR}(K^{+}(\pi^{+}))=\mathcal{L}(K^{+}(\pi^{+}))/(\mathcal{L}(\pi^{+})+\mathcal{L}(K^{+})) and the beam energy constrained mass m bc subscript 𝑚 bc m_{{\rm{bc}}} [ 10 ] .

particle	pre selection applied
$\gamma$	$E_{\gamma }>0.075\mathrm{GeV}$
${\pi }^0$	$0.145\mathrm{GeV}>m\left(\gamma \gamma \right)>0.125\mathrm{GeV}$
${\pi }^{+}$	$ℒℛ\left({\pi }^{+}\right)>0.5$
$K^{+}$	$ℒℛ\left(K^{+}\right)>0.5$
$D^0$	$1.9\mathrm{GeV}>m\left(K^{-}{\pi }^{+}{\pi }^0\right)>1.7\mathrm{GeV}$
$D^{\ast }$
$B^0$	$m_{\mathrm{bc}}=\sqrt{E_{\mathrm{beam}}^2/4-p^2}>5.27\mathrm{GeV}$

Equations

$$\mathbf{x}=\{x_1,y_1,z_1,{\theta }_1,p_{x1},p_{y1},p_{z1},E_1,...,x_n,y_n,z_n,{\theta }_n,p_{xn},p_{yn},p_{zn},E_n\},$$

$${\chi }_{\alpha }^2=$$

$$+$$

$${\mathbf{x}}_k^{\alpha }={\mathbf{x}}_{k-1}^{\alpha }-{\mathbf{K}}_k^{\alpha }{\mathbf{r}}_k^{\alpha },$$

$${\mathbf{r}}_k^{\alpha }={\mathbf{m}}_k-{\mathbf{h}}_k({\mathbf{x}}^{\alpha }).$$

$${\mathbf{H}}_k^{\alpha -1}=\frac{\partial {\mathbf{h}}_k({\mathbf{x}}^{\alpha -1})}{\partial {\mathbf{x}}^{\alpha -1}},$$

$${\mathbf{h}}_k({\mathbf{x}}^{\alpha })={\mathbf{h}}_k({\mathbf{x}}^{\alpha -1})-{\mathbf{H}}_k^{\alpha -1}\cdot ({\mathbf{x}}^{\alpha }-{\mathbf{x}}^{\alpha -1}).$$

$${\mathbf{r}}_k^{\alpha }={\mathbf{m}}_k-{\mathbf{h}}_k({\mathbf{x}}^{\alpha -1})+{\mathbf{H}}_k^{\alpha -1}\cdot ({\mathbf{x}}^{\alpha }-{\mathbf{x}}^{\alpha -1}).$$

$${\mathbf{K}}_k={\mathbf{C}}_{k-1}{\mathbf{H}}_k^T{\mathbf{R}}_k^{-1}.$$

$${\mathbf{C}}_k=$$

$$=$$

$${\mathbf{R}}_k={\mathbf{V}}_k+{\mathbf{H}}_k{\mathbf{C}}_{k-1}{\mathbf{H}}_k^T.$$

$$\frac{\partial\bm{h}}{\partial\bm{x}}=\bm{H}=\left(\begin{array}[]{ccccc}&\\ \mbox{\Large 0}&&\Large{\bm{\mathbbm{1}}}_{3}&&\mbox{\Large 0}\\ &&&&\\ \end{array}\right)\leavevmode\nobreak\ .$$

$${\mathbf{h}}_{\mathrm{track}}(\mathbf{x})=\left(\begin{array}{c}{d}_0\\ {\varphi }_0\\ \omega \\ z_0\\ \tan \lambda \end{array}\right)=\left(\begin{array}{c}A(1+U{)}^{-1}\\ atan2(p_y,p_x)-atan2(\omega \cdot {\Delta }_{\parallel },1+\omega \cdot {\Delta }_{\perp })\\ a\cdot q/p_t\\ z+l\cdot \tan \lambda \\ p_z/p_t\end{array}\right).$$

$${\mathbf{r}}_{track}^{\alpha }(\mathbf{x})=\left(\begin{array}{c}{d}_{0,m}-d_0\\ {\varphi }_{0,m}-{\varphi }_0\\ {\omega }_m-\omega \\ z_{0,m}-z_0\\ \tan {\lambda }_m-\tan \lambda \end{array}\right)+{\mathbf{H}}^{\alpha -1}\cdot ({\mathbf{x}}^{\alpha }-{\mathbf{x}}^{\alpha -1})\text{leavevmode}.$$

$$\bm{H}=\left(\begin{array}[]{ccccccccc}...&&&&...&&&&...\\ ...&&&&...&&&&...\\ ...&&\bm{A}&&...&&\bm{B}&&...\\ ...&&&&...&&&&...\\ ...&&&&...&&&&...\\ \end{array}\right)\leavevmode\nobreak\ .$$

$${\mathbf{A}}_{d_0,x}$$

$${\mathbf{A}}_{{\varphi }_0,x}$$

$${\mathbf{A}}_{z_0,x}$$

$${\mathbf{B}}_{d_0,p_x}=\begin{array}{rl}& -\frac{y((aq{)}^{2}r+2aq{p}_{y}x+2p_y^2\beta )}{p_{t0}{p}_t^2{\beta }^2}\\ & -\frac{p_x(2p_{y}x\beta +aq(y^2(-2+\beta )+x^2\beta )}{p_{t0}{p}_t^2{\beta }^2},\end{array}$$

$${\mathbf{B}}_{d_0,p_y}=\begin{array}{rl}& \frac{2p_x^{2}x\beta +2p_{x}y(p_y-aqx+p_y{p}_{t0}/p_t)}{p_{t0}{p}_t^2{\beta }^2}\\ & +\frac{aq(aqrx-p_y(x^2(-2+\beta )+y^2\beta ))}{p_{t0}{p}_t^2{\beta }^2},\end{array}$$

$${\mathbf{B}}_{\omega ,p_x}$$

$${\mathbf{B}}_{z_0,p_x}$$

$${\mathbf{B}}_{z_0,p_y}$$

$${\mathbf{B}}_{\tan \lambda ,p_x}$$

$$0={\mathbf{u}}_{\mathrm{mother}}+\mathbf{\delta }-\mathbf{m},$$

$${\mathbf{h}}_{\mathrm{photon}}(\mathbf{x})=\left(\begin{array}{c}{u}_x+\tau \cdot p_x\\ u_y+\tau \cdot p_y\\ u_z+\tau \cdot p_z\\ \sqrt{p_x^2+p_y^2+p_z^2}\end{array}\right)\mathrm{and}{\mathbf{m}}_{\mathrm{photon}}(\mathbf{x})=\left(\begin{array}{c}{m}_x\\ m_y\\ m_z\\ E_m\end{array}\right),$$

$${\mathbf{r}}_{\mathrm{photon}}^{\prime \alpha }(\mathbf{x})=\left(\begin{array}{c}(m_i-u_i)-(m_k-u_k)\frac{p_i}{p_k}\\ (m_j-u_j)-(m_k-u_k)\frac{p_j}{p_k}\\ E_m-\sqrt{p_i^2+p_j^2+p_k^2}\end{array}\right)+{\mathbf{H}}^{\alpha -1}\cdot ({\mathbf{x}}^{\alpha }-{\mathbf{x}}^{\alpha -1})\text{leavevmode},$$

$${\mathbf{A}}_{0,u_k}$$

$${\mathbf{A}}_{1,u_k}$$

$${\mathbf{B}}_{0,p_k}$$

$${\mathbf{B}}_{1,p_k}$$

$${\mathbf{B}}_{2,p_k}$$

$$\bm{H}=\left(\begin{array}[]{ccccccccc}...&&&&...&&&&...\\ ...&&\bm{A}&&...&&\bm{B}&&...\\ ...&&&&...&&&&...\\ \end{array}\right)\leavevmode\nobreak\ .$$

$${\mathbf{V}}^{\prime }=\mathbf{F}\mathbf{V}{\mathbf{F}}^T,$$

Full text

††journal: Elsevier Science

1 Introduction

Particle vertex fitting techniques are widely used in particle and nuclear physics. Beyond the suppression of background, applications range from the improvement of particle momentum resolution (under the assumption they originate from some vertex point), to the determination of the presence of intermediate particles and the precision determination of decay vertex positions. One can, for example, combine the measurements of two charged pion tracks originating from the decay of a $K^{0}_{\rm{S}}$ to extract the decay vertex position, flight length and four-vector and their uncertainties. By performing a kinematic fit, one obtains an improvement of the pion track momenta and can use the $\chi^{2}$ probability of the fit result to suppress background.

In order to construct more complex decay topologies, one usually combines cascades of these fits starting with long lived stable particles such as electrons, muons, pions and photons, forming intermediate resonances and finally the full decay of interest. For example, in the decay $B^{0}\to J/\Psi K^{0}_{\rm{S}}$, where $J/\Psi\to\mu^{+}\mu^{-}$ and $K^{0}_{\rm{S}}\to\pi^{+}\pi^{-}$, one would first fit the $J/\Psi$ and $K^{0}_{\rm{S}}$ candidates and then use these to construct the $B^{0}$ candidates, as depicted in Fig. 1a and Fig. 1b. However, this only works well if the final state particles are charged and leave traces in the tracking detectors. Neutral particles can not be tracked; they are only measured by their energy deposition in the calorimeter. Single layer crystal calorimeters do not offer directional information on where the particle associated to an energy deposition originated from. This means that the decay vertex can not be extracted from the fit. In order to obtain the momentum vector, it has to be assumed that the particle originates from the primary interaction point and travels directly into the cluster’s center of gravity. This can introduce a large bias on the momentum components. Consider, for example, the decay $B^{0}\to J/\Psi K^{0}_{\rm{S}}$, where the kaon instead decays to neutral final states, $K^{0}_{\rm{S}}\to\pi^{0}\pi^{0}$, and the pions decay, $\pi^{0}\to\gamma\gamma$, as displayed in Fig. 2. Weakly decaying intermediate particles such as $K^{0}_{\rm{S}}$ have relatively long flight lengths, up to $15\mathrm{\leavevmode\nobreak\ cm}$ for $K^{0}_{\rm{S}}$ at Belle II. Thus for a neutral particle in such a particle decay chain, the assumption that it originates from the primary $e^{+}e^{-}$ collision is not sufficient.

The method we present in this paper overcomes these issues by globally fitting the entire decay tree in a single fit, taking into account all intermediate particles, extracting all involved particle’s four-momenta, vertex positions, flight lengths and their covariance matrices, using a Kalman Filter as described in Ref. [1]. We use the software environment of Belle II and rely on the C**++**  template library EIGEN [2] for matrix operations, which provides a fast execution time for the fit algorithm. We furthermore present physics applications of the fitter with Belle II Monte Carlo samples.

1.1 The Belle II Experiment

The algorithm described in this paper was developed for the analysis software framework of Belle II. The Belle II experiment takes place at the asymmetric $e^{+}e^{-}$ collider, SuperKEKB. SuperKEKB provides a beam energy slightly above the mass of the $\Upsilon(4\rm{S})$ resonance $(10.58\mathrm{\leavevmode\nobreak\ GeV}/c^{2})$ at an instantaneous luminosity of $8\cdot 10^{35}\mathrm{\leavevmode\nobreak\ cm^{-2}s^{-1}}$. The $\Upsilon(4\rm{S})$ resonance decays into pairs of $B$-mesons just above production threshold, hence this type of experiment is referred to as a $B$-factory. The asymmetric beam energy gives the $B$-meson a relativistic boost along a direction close to the detector’s axis of symmetry, increasing its flight length in the lab frame, which makes it easier to study the time evolution of $B$ decays - a key observable in the study of CP symmetry violation.

The Belle II detector has a cylindrical structure designed to study the decays of $B$- and $D$-mesons, $\tau$-leptons and other processes produced in $e^{+}e^{-}$ collisions. Six layers of silicon vertex detectors (2 layers of silicon pixels (PXD), and 4 layers of double sided silicon detectors (SVD)) are located in the central volume of the detector, designed to accurately track the flight paths of charged particles. The following layers are, a central drift chamber (CDC) used to measure track trajectories within a solenoid magnetic field, Cherenkov light based particle identification devices surrounding the CDC in the barrel (TOP) and forward regions (ARICH), followed by the CsI(Tl) electromagnetic calorimeter (ECL). The outermost layers are composed of, a magnet solenoid, a $K_{\rm{L}}^{0}$ and muon detector system (KLM), which is also used as the flux return yoke of the magnet. The magnetic field is aligned along the detector’s axis of symmetry. For a more detailed description of the detector see Ref. [3].

Fig. 3 and Fig. 4 show an event display depicting simulated particles traversing the Belle II detector. In the decay $\Upsilon(4{\rm{S}})\to\bar{B}^{0}B^{0}$, one meson decays as $\bar{B}^{0}\to D\omega\pi\pi$, and the other as $B^{0}\to K_{\rm{S}}J/\Psi$, with $J/\Psi\to\mu^{+}\mu^{-}$ and $K^{0}_{\rm{S}}\to\pi^{0}\pi^{0}$ with $\pi^{0}\to\gamma\gamma$. Fig. 3, depicts the full detector geometry and Fig. 4 shows a close-up of the inner vertex detectors. In this example we show that the decay vertex of the $K_{\rm{S}}^{0}$ can be very displaced.

2 Extended Kalman Filter

Vertex fitting is a least squares minimization problem. The best least squares estimator (LSE) is the solution to this problem. The computational challenge in finding a LSE lies in matrix inversions, which naively scale as $\mathcal{O}(n^{3})$, where $n$ is the dimension of the matrix. In a naive approach, this is equal to the number of parameters extracted in the fit. A Kalman Filter is an iterative approach to find the LSE by defining a series of constraints (knowledge of parameters from measurements and symmetries) on a hypothesis (a set of particle parameters in this case). The hypothesis has different states during the stages of the filtering process. The filtering process is an iterative algorithm applying the constraints sequentially and updating the state with respect to the constraints. The sequence is repeated until convergence is reached or divergence is observed. In the case of convergence, the last state describes the best hypothesis for the parameters that can be found. We use an Extended Kalman Filter in the gain matrix formulation [4, 5, 1] in our algorithm. The state vector $\bm{x}$ holds the particle parameters to be optimized by the fit. The parameters depend on the type of particle. The most general parametrization takes the form

[TABLE]

where vertex coordinates of the $i$-th particle are denoted as $\{x_{i},y_{i},z_{i}\}$, its decay length is denoted as $\theta_{i}$ and the four momentum is $\{p_{x,i},p_{y,i},p_{z,i},E_{i}\}$, and $n$ is the number of particles in the fitted topology. This implies that $n\approx 8\times\text{number of particles}$. The problem is split into $k$ constraints given by measurements and other knowledge of parameters (see Sec. 3 for the definitions of the measurement $\bm{m}$ and hypothesis $\bm{h}$ for each of the different constraints). The minimized $\chi^{2}$ can be expressed as a weighted sum of $k$ sets of equations, which are the constraints. For a given iteration, $\alpha$, of the Kalman Filter, one can write

[TABLE]

using the measurement covariance matrix $\bm{V}_{k}^{-1}$ transported to the $\bm{m}-\bm{h}$ system as a weight, such that the hypothesis of each constraint $k$ depends on the outcome of the previous constraint.

Minimizing Eq. 3 yields a rule to find a new state $\bm{x}_{k}^{\alpha}$ for the current iteration and constraint, that is

[TABLE]

where $\bm{K}$ is a gain matrix, which will be defined below. However, we first define the residual $r$ as the distance between measurement and hypothesis

[TABLE]

The hypothesis $\bm{h}_{k}(\bm{x}^{\alpha})$ can be linearised around a reference state of the previous iteration $\bm{x}^{\alpha-1}$ using the Jacobian

[TABLE]

such that

[TABLE]

Eq. 5 thus becomes

[TABLE]

The gain matrix is calculated for every constraint in every iteration and is defined as

[TABLE]

The only matrix to be inverted is of the dimension of the constraint and not the dimension of the statespace. The current state’s covariance matrix is obtained via propagation of uncertainties

[TABLE]

The covariance matrix of the residual system can be found to be

[TABLE]

For constraints where no measurement is available, for example the mass constraint, we set $\bm{V}_{k}=0$.

Linearising the hypothesis around a reference state, see Eq. 7, makes the fits more stable than a normal Kalman Filter and reduces the overall number of signal candidates for which the fit diverges. In this formulation of the LSE problem, the update of the covariance matrix, Eq. 10, is the computational bottleneck, as it is a dense $m\times m$ matrix with the dimension of the state vector $m$ and must be calculated and filled $k$ times for each iteration.

3 Parametrising and constraining the decay chain

To parametrise the decay chain, we use a set of parameters that describe the properties of the particles. We perform a number of reductions on these properties to reduce the dimensionality of the problem. For final state particles, we only save the momenta, as they do not have decay vertices. For their production vertices we use the decay vertices of their respective composite particles, later referred to as “mothers”. The energy of each each final state particle is calculated using the momenta and it’s nominal mass hypotheses, taking the mass values provided by the particle data group (PDG) [6]. Intermediate particles are classified in two categories: particles that decay dominantly via the strong force and via the weak force. Strongly decaying particles with a lifetime of less than $10^{-14}\mathrm{\leavevmode\nobreak\ s}$, which corresponds to boosted flight lengths of the order of less than $1\mathrm{\leavevmode\nobreak\ \mu m}$, are treated as if they instantly decay (the detector has a vertex position resolution in the x-y plane for charged particles of $20-30\mathrm{\leavevmode\nobreak\ \mu m}$). The hypothesised quantities are their energy and momenta, while the production and decay vertices coincide and are taken from their mother particle’s decay vertex. For weakly decaying particles, we additionally measure a decay vertex and a flight length, defined as the distance between the production and decay vertices in three dimensions.

3.1 Parametrising the constraints

Constraints are defined by Eq. 7. The resulting Jacobians take the form of $m\times n$ matrices where $m$ is the dimension of the state vector and $n$ is the dimension of the respective constraint. Thus, only few of its elements are non-zero. For example, for a three dimensional point constraint $k$, the hypothesis of particle number $j$ with $\bm{h}_{j}=\{x_{j},y_{j},z_{j}\}$ and $\bm{x}$ as in Eq. 2, only the $j$-th diagonal block is non-zero

[TABLE]

The blocks filled with zero correspond to the parameters of particle $\bm{x}_{i}\neq\bm{x}_{j}$. We will omit the columns filled with zeros throughout this section, for brevity.

In the following we list the definitions of constraints that have been implemented in the Belle II software, based on the specific geometry of Belle II.

3.1.1 Reconstructed track

A track can be parametrised with a five parameter helix. In Belle II it was chosen to use a perigee-parametrised helix, such that the helix is defined at the perigee, the point of closest approach of the helix to the origin of the coordinate system. The corresponding transformations to transport a helix to that point are discussed in Ref. [7]. A description of the parameters can be found in Table 1, and a depiction of the helix is in Fig. 5. We parametrise tracks such that we can express the model’s dependence on Cartesian parameters as

[TABLE]

Where atan2 refers to the phi domain corrected inverse tangent function. We use the same parametrisation for the hypothesis and the measurement. We label the measurement quantities with the index $m$. The residuals of iteration $\alpha$ then become

[TABLE]

We define the Jacobian block $\bm{A}:=\partial\bm{h}/\partial\bm{x}$ as the derivatives with respect to the vertex position, and $\bm{B}:=\partial\bm{h}/\partial\bm{p}$ as the derivatives with respect to momentum. The positions of these blocks in the Jacobian depend on the topology fitted and the particle represented by the track. We choose to order the state vector hierarchically. This means that the decay vertex parameters come before its momentum, followed by the daughter particle’s parameters. The full Jacobian $\bm{H}$ then takes the following form

[TABLE]

For the non-zero elements of the Jacobian blocks, denoted by $\partial d_{0}/\partial x=\bm{A}_{d_{0},x}$, we derive the spatial components as

[TABLE]

and for the momenta

[TABLE]

and

[TABLE]

3.1.2 Reconstructed photon

For photons we measure the position of the calorimeter cluster and its energy and can infer the vertex parameters. The geometry, depicted in Fig. 6, gives

[TABLE]

substituting $\bm{\delta}=\tau\cdot\bm{p}$ and inserting the energy relation, we get

[TABLE]

where $\{u_{x},u_{y},u_{z}\}$ are the production vertex coordinates, $\{p_{x},p_{y},p_{z}\}$ are the parameters of the momentum vector pointing from the production vertex to the calorimeter cluster, $\{m_{{x}},m_{{y}},m_{{z}},E_{{m}}\}$ are the position and measured energy of the corresponding ECL cluster. The parameter $\tau$ is the decay time, which can be eliminated when writing down the residual to reduce the dimensionality of the equation system and avoid a trivial local minimum of $\bm{r}_{\gamma}$ at $\tau=0$ when taking $\{u_{x},u_{y},u_{z}\}=0$ as the starting point of the first iteration. Since the geometry of the detector is cylindrical, we can not simply eliminate any of the dimensions as this could introduce a pole in the residual equations. Thus, we sort the momenta and eliminate the dimension with the highest momentum. Such that we get a 3-dimensional equation system

[TABLE]

where the indices $i,j,k$ indicate the dimensions by order of increasing momentum $p_{k}\geq p_{i}\geq p_{j}$. We define $\bm{A}_{i,u_{k}}:=\partial h^{\prime}_{i}/\partial u_{k}$ and $\bm{B}_{i,p_{k}}:=\partial h^{\prime}_{i}/\partial p_{k}$ with the hypothesis of the reduced system $r^{\prime}$. Thus, the non-zero entries are

[TABLE]

The full Jacobian then takes the form

[TABLE]

We must transform the covariance matrix of the measurement into the reduced system. For that, we use

[TABLE]

with the transport matrix $\bm{F}=\partial r^{\prime}/\partial m$, which depends on the sorting of the momenta such that the non-zero entries are

[TABLE]

We do not parametrise this constraint in $p_{t}$ in order to keep the derivatives in Eq. 22 as computationally simple as possible.

3.1.3 Reconstructed $K^{0}_{L}$

We treat $K^{0}_{L}$ in the same way as photons, except that we use the nominal mass provided by the PDG in the energy calculation. The KLM detector is used for the cluster position measurement instead of the calorimeter. It cannot provide a precise energy measurement. Instead, it extrapolates the energy deposited by a particle as $E=c\cdot n$, where $n$ is the number of hit cells in the cluster and $c$ is a constant with the units $\mathrm{GeV}$. This approach makes the energy measurement for $K^{0}_{L}$ much less resolved than for photons.

3.1.4 Kinematic constraint

The kinematic constraint enforces four-momentum conservation, meaning it fits the four-momentum of the mother as the sum of the daughter momenta

[TABLE]

The Jacobian for a particle with $two$ daughters, which in this example are taken to be stable particles, can be defined

[TABLE]

with the blocks $\bm{A},\bm{B}$ as

[TABLE]

and

[TABLE]

Note that the energy row of $\bm{B}_{i}$ depends on how the particle is parametrised, composite particles, for example, are parametrised with an energy variable in the state vector, resulting in $\bm{B}=-\bm{\mathbbm{1}}_{4}$, while for stable particles Eq. 29 is used.

3.1.5 Geometric constraint

The geometric constraint fits the decay length parameter $\theta$ for composite particles, see Fig. 6. Accounting for the geometry we have

[TABLE]

Instead of directly extracting a flight vector $\bm{\Delta}$, we use the unit vector of the momentum as it is well constrained by the previously filtered kinematic constraints, substituting $\bm{\Delta}=\theta\cdot\bm{p}/|\bm{p}|$, allows for a more accurate estimation of $\theta$. Thus we define the residual as

[TABLE]

using

[TABLE]

and

[TABLE]

such that

[TABLE]

3.1.6 Mass constraint

The mass constraint requires a particle four-vector to be consistent with its nominal mass. We treat the particle as a measurement with infinite precision and use the mass value provided by PDG such that

[TABLE]

with

[TABLE]

3.1.7 Beam spot constraint

$B$-mesons have a short lifetime and decay very close the beam spot, therefore it is useful to constrain the reconstructed particles to a volume within that region. The constraint is implemented by considering the initial $e^{+}e^{-}$ collision as an abstract mother particle with the parameters and uncertainties of the beam spot, thus constraining the $production$ vertex position of its daughters to an area within the beam spot’s three dimensional uncertainty region using

[TABLE]

with the Jacobian

[TABLE]

Here $\bm{s}$ denotes the beam spot position vector and $\bm{h}_{\rm{production}}$ is the fitted production vertex of the particle. The beam spot is determined as part of the detector calibration, by averaging the $x$, $y$ and $z$ positions of the $e^{+}e^{-}$ collisions over many interactions. The $decay$ vertex can be constrained similarly.

3.1.8 Custom Origin constraint

Similar to the beam spot constraint, one can construct another geometric constraint by defining a custom vertex position and an associated uncertainty to be the origin of the decay chain. This can be very useful when the decay chain contains particles that can not be detected, for example, neutrinos or long lived dark matter particles. In these decay chains it is not useful to fit the full chain, because the missing particle’s four momentum makes the kinematic constraint of the mother particle an incorrect assumption. However, knowing that the particle originates from a $B$-meson decay, one can define a geometric constraint corresponding to the volume where $B$-mesons decay on average.

A beam energy constraint is not needed, as one can mass constrain the $\Upsilon(4S)$ particle, since all four-vectors in $e^{+}e^{-}\to\Upsilon(4\rm{S})$ are well known.

4 Fitting complex decay topologies and extracting parameters from the fit

In this section we present two use cases of the algorithm containing one and two $\pi^{0}$-mesons in a decay chain, respectively. There are numerous channels where the phase space is large enough to add $\pi^{0}$ to a vertex with two charged tracks, which makes this a very common structure in decay trees and therefore an interesting target to fit in a wide spectrum of analyses. We use Monte Carlo samples from the Belle II experiment and evaluate background rejection and signal resolution improvements, as well as the extraction of various other parameters. The selection criteria used for all analyses in this section are listed in Table 2. We use implicit charge conjugation, meaning that charged particles appearing in the decay chains imply including the opposite charges. We abbreviate decays $x\to yz$ as $x(yz)$.

4.1 Fitting decay chains containing a neutral particle

The first mode, with a single $\pi^{0}$, is $B^{0}\to D^{*+}({D}^{0}(K^{-}\pi^{+}\pi^{0})\pi^{+})\pi^{-}$. This channel is of interest as it is the most common $D^{0}$ decay modes. We use the fitter in this example to suppress background. If we reject combinations for which the fit failed - i.e. those with a p-value $\leq 0$ - we are able to reject about $20\%$ of the background while rejecting only about $4\%$ signal, as displayed in Fig. 7(a). Note that since some of the constraints are nonlinear, the $\chi^{2}$ which is used to calculate the p-values can be degenerate, where the effective number of degrees of freedom is smaller than the theoretical one used to calculate the p-values. As a result, the p-values can have non-flat shapes [9].

4.2 Fitting ill-defined decay vertices

We study the decay chain $B^{0}\to K^{0}_{\rm{S}}(\pi^{0}(\gamma\gamma)\pi^{0}(\gamma\gamma))J/\Psi(\mu^{+}\mu^{-})$. By performing the fit a large amount of the background can be removed by requiring that it passes the fit, see Fig. 8(a). The only well defined vertex in this chain is given by the $J/\Psi\to\mu^{+}\mu^{-}$ decay. The other vertices are ill-defined, the best assumption that can be made for the production vertex position of the four photons is that they originate from the interaction point. This assumption biases the reconstructed mass of the $K^{0}_{\rm{S}}$. Performing a fit with mass constrained $\pi^{0}$-mesons improves the extracted mass of the $K^{0}_{\rm{S}}$, so that after the fit, it is centred around the true value, as depicted in Fig. 8(b). It is then possible to further reject background outside the nominal mass window and improve the signal purity when analysing this channel.

4.3 Extracting the decay length of $D^{0}$-mesons from $D^{*+}$ decays using a beam-spot constraint

The geometric constraints allow for the extraction of all the production and decay vertices of all particles in the decay chain. This allows for the extraction of flight lengths and thus lifetimes of intermediate particles such as $D^{0}$-mesons. We perform this study on $B^{0}\to D^{*+}({D}^{0}(K^{-}\pi^{+})\pi^{+})\pi^{-}$ decays and extract the decay length of the $D^{0}$-meson. The results are depicted in Fig. 9. Since $D^{*+}$-mesons decay almost instantly, we will use the three dimensional distance between the $B^{0}$ and the $D^{0}$ decay vertices. To improve the resolution on the $B^{0}$ vertex we apply a beam-spot constraint.

4.4 Using a custom origin constraint to improve the $D^{*+}$ vertex resolution in $B^{0}\to D^{*+}({D}^{0}(K^{-}\pi^{+})\pi^{+})l^{-}\bar{\nu}$ decays

In some scenarios one can not fit the entire decay chain e.g. in semi-leptonic decays $B^{0}\to D^{*}\mu\nu_{\mu}$, where the neutrino escapes detection. Since the neutrino’s four vector is unknown it is not possible to use a kinematic constraint to fit the $B$-meson. To circumvent this problem, one can omit fitting the $B^{0}$, and instead account for its presence with a geometric constraint. In $B^{0}\to D^{*+}({D}^{0}(K^{-}\pi^{+})\pi^{+})l^{-}\bar{\nu}$ decays we use a custom origin constraint to constrain the $D^{*+}$ production vertices to a volume around the median $B$-meson decay vertex position obtained from Monte Carlo simulations. The results are depicted in Fig. 10. This technique can improve the resolution on the $D^{*+}$ vertex position.

5 Conclusion

We presented a versatile fitting tool tailored for the environment of $B$ factories with a cylindrical detector geometry. It can be used for various purposes, such as the extraction of particle production and decay vertices, decay lengths, particle four-momenta and rejection of backgrounds, as well as the extraction of the respective uncertainties for all parameters involved in three dimensions. This global fitting technique is particularly powerful in fitting and reducing background in modes that contain neutral particles, and will be very important for the Belle II physics analysis program. As a future extension of the fitter we are looking into the possibility to map the residuals found by the Kalman Filter to a scalar signal probability $f:\bm{r}\to p$. Such a probability may outperform methods using a p-value derived from the fit $\chi^{2}$ for the purpose of signal-background-separation.

6 Acknowledgements

This work was made possible by the Belle II collaboration and funding from ARC (Australia), ARRS (Slovenia), BMBF (Germany), EXC153 (Germany), HGF (Germany), INFN (Italy), MSMT (Czech Republic), NSERC (Canada) and U.S. DOE. We would like to thank Wouter Hulsbergen for his original paper, the discussion and sharing his code.