$$\text{Var}\begin{bmatrix}\begin{pmatrix}\Delta x\\ \Delta \mathbf{y}\\ \Delta t\_{\mathbf{x}}\\ \Delta t\_{\mathbf{y}}\end{pmatrix}\end{bmatrix} = \sigma\_{0}^{2}(1+t\_{\mathbf{x}}^{2}+t\_{\mathbf{y}}^{2})\begin{pmatrix}(1+t\_{\mathbf{x}}^{2})D^{3}/3 & t\_{\mathbf{x}}t\_{\mathbf{y}}D^{3}/3 & (1+t\_{\mathbf{x}}^{2})D^{2}/2 & t\_{\mathbf{x}}t\_{\mathbf{y}}D^{2}/2\\ t\_{\mathbf{x}}t\_{\mathbf{y}}D^{3}/3 & (1+t\_{\mathbf{y}}^{2})D^{3}/3 & t\_{\mathbf{x}}t\_{\mathbf{y}}D^{2}/2 & (1+t\_{\mathbf{y}}^{2})D^{2}/2\\ (1+t\_{\mathbf{x}}^{2})D^{2}/2 & t\_{\mathbf{x}}t\_{\mathbf{y}}D^{2}/2 & (1+t\_{\mathbf{x}}^{2})D & t\_{\mathbf{x}}t\_{\mathbf{y}}D\\ t\_{\mathbf{x}}t\_{\mathbf{y}}D^{2}/2 & (1+t\_{\mathbf{y}}^{2})D^{2}/2 & t\_{\mathbf{x}}t\_{\mathbf{y}}D & (1+t\_{\mathbf{y}}^{2})D\end{pmatrix},$$

where *D* = *(*1 + *t* 2 *x* + *t* 2 *y )*1*/*2*d* is the effective thickness crossed. If the direction is represented by *θ* and *ϕ*, the covariance matrix can be computed via the transformation

$$t\_{\chi} = \tan(\theta)\cos(\varphi), \quad t\_{\mathcal{Y}} = \tan(\theta)\sin(\varphi),$$

and linear error propagation with the Jacobian

$$T = \frac{\partial(\theta, \varphi)}{\partial(t\_{\mathcal{X}}, t\_{\mathcal{Y}})} = \begin{pmatrix} \frac{t\_{\mathcal{X}}}{\sqrt{t\_{\mathcal{X}}^2 + t\_{\mathcal{Y}}^2 (1 + t\_{\mathcal{X}}^2 + t\_{\mathcal{Y}}^2)}} & \frac{t\_{\mathcal{Y}}}{\sqrt{t\_{\mathcal{X}}^2 + t\_{\mathcal{Y}}^2 (1 + t\_{\mathcal{X}}^2 + t\_{\mathcal{Y}}^2)}}\\ -\frac{t\_{\mathcal{Y}}}{t\_{\mathcal{X}}^2 + t\_{\mathcal{Y}}^2} & \frac{t\_{\mathcal{X}}}{t\_{\mathcal{X}}^2 + t\_{\mathcal{Y}}^2} \end{pmatrix} = \begin{pmatrix} \frac{\cos(\varphi)}{1 + \tan^2(\theta)} & \frac{\sin(\varphi)}{1 + \tan^2(\theta)}\\ -\frac{\sin(\varphi)}{\tan(\theta)} & \frac{\cos(\varphi)}{\tan(\theta)} \end{pmatrix}.$$

For analogous formulas in cylindrical coordinates, see [46].

If the curvature of the track cannot be neglected, the simplest approach is a stepwise integration of the equation of motion, assuming the validity of a helical track model within each step and considering each such step as a thin scatterer [39, 47].

#### Energy Loss

For particles other than electrons the energy loss in material is almost exclusively due to scattering on electrons. The momentum correction *p* in a material layer of thickness *d* is calculated by integrating the Bethe-Bloch formula [37]:

$$\Delta p = \int\_0^d \frac{\mathrm{d}p}{\mathrm{d}x} \mathrm{d}x = \int\_0^d \frac{1}{\beta} \frac{\mathrm{d}E}{\mathrm{d}x} \mathrm{d}x = \int\_0^d \frac{\mathrm{K}}{\beta^3} \left[ \ln \frac{2m\_ec^2 \beta^2 \gamma^2}{\langle I \rangle} - \beta^2 \right] \mathrm{d}x,\tag{13.4}$$

where *K* is a constant depending on the material, *m*e is the electron mass, *I* is the average ionization potential of the material, and *β* = *v/c* and *γ* = *E/mc*2 are the usual kinematic parameters. The ratio *I /Z* is about 20 eV for hydrogen and helium, between 12 and 16 eV for light nuclei, and around 10 eV for heavy nuclei [44]. For practical purposes, the differential energy loss d*E/*d*x* is a function only of *β*. For small *β*, it decreases like 1*/β*2. It has a minimum, the position of which drops with increasing *Z* from *βγ* ≈ 3*.*5 (carbon) to *βγ* ≈ 3 (lead). In terms of momentum, the minimum is at *p* = *βγ mc* and thus depends on the mass of the particle. This dependency is used for particle identification. The energy loss at the minimum can be parameterized for *Z* ≥ 6 by [44]:

$$(\mathrm{d}E/\mathrm{d}\mathbf{x})\_{\mathrm{min}} = (2.35 - 0.64 \ln \mathrm{l} \,\mathrm{0} \,\mathrm{Z}) \,\mathrm{MeV} \,\mathrm{g}^{-1} \mathrm{cm}^2.$$

From this the constant *K* in Eq. (13.4) can be calculated. For large *βγ* the energy loss increases like ln*(βγ )*; this is called the relativistic rise. For momenta in the vicinity of the minimum d*E/*d*x* can be considered as constant, giving *p* ≈ *(*d*E/*d*x)*min · *dρ/β*, *ρ* being the density of the material.

#### Bremsstrahlung

For an electron (or positron) passing through matter the most significant contribution to energy loss is bremsstrahlung, the emission of photons in the electric field of an atomic nucleus. In the Bethe–Heitler model [48] the relative energy loss is distributed independently of the energy. Let *d* be the path length in the material in units of radiation length, and *z* the fraction of energy remaining after the material is traversed. Then the distribution of *z* is given by the following probability density function:

$$f(z) = \frac{(-\ln z)^{c-1}}{\Gamma(c)}, \quad 0 \le z \le 1,$$

where *(x)* is Euler's gamma function and *c* = *d/*ln 2. For high energy electrons *p* ≈ *E*, so the momentum correction is *p* ≈ *p(z* − 1*)*. The first two moments of *p* are

$$E(\Delta p) = p(2^{-c} - 1), \quad \text{var}(\Delta p) = p^2(3^{-c} - 4^{-c}).$$

The moments can be used for a Gaussian representation of bremsstrahlung as an additional process noise in the Kalman filter (see Sect. <IP_ADDRESS>). As this is a very crude approximation, more sophisticated methods have been developed that take into account the actual shape of the distribution. One of them is the Gaussiansum filter [49, 50], see Sect. <IP_ADDRESS>. A computationally less intensive approach is described in [51].

#### **<IP_ADDRESS> Estimation Methods**

The main task of the track fit is to estimate the values of a set of parameters describing the state of a particle somewhere in the detector, often at a reference surface close to the interaction vertex. The information from the measurements created by the particle while traversing the tracking detector should be processed in an optimal manner. If the track model is truly linear, i.e., if the measurements are strictly linear functions of the track parameters, and all stochastic disturbances entering the estimation procedure are Gaussian, the linear least-squares method is the optimal one [37]. Since track parameter propagation in general is a nonlinear procedure, strict linearity holds very rarely in practice. The relation between the track parameter vector *q*0 at a reference surface and the measurement vector *mk* at a detector layer *k* is a function *dk* given by

$$
\mathfrak{m}\_k = d\_k(\mathfrak{q}\_0) + \mathfrak{p}\_k,
$$

where *γ k* is a noise term containing the measurement error of *mk* and all multiple scattering in front of *mk*. The function *dk* is a composition of the measurement model function *hk* and the track propagator functions *f i*|*i*−1 (see Sect. <IP_ADDRESS>):

$$d\_k = h\_k \circ f\_{k|k-1} \circ \cdots \circ f\_{2|1} \circ f\_{1|0}.$$

For the linear least-squares method *dk* has to be linearized around some expansion point, providing the Jacobian *Dk* of each *dk*:

$$D\_k = H\_k F\_{k|k-1} \cdots F\_{2|1} F\_{1|0},$$

with *Hk* from Eq. (13.3). The covariance matrix of *γ k* is obtained by linear error propagation:

$$\text{var}(\mathbf{y}\_k) = \mathbf{V}\_k + \mathbf{H}\_k (\mathbf{F}\_{k|1} \mathbf{Q}\_1 \mathbf{F}\_{k|1} ^\mathrm{T} + \dots + \mathbf{F}\_{k|k-1} \mathbf{Q}\_{k-1} \mathbf{F}\_{k|k-1} ^\mathrm{T} + \mathbf{Q}\_k) \mathbf{H}\_k ^\mathrm{T},$$

where *V k* is the covariance matrix of the measurement error *εk* of *mk*, and *Qj* is the covariance matrix of multiple scattering after layer *j* −1 up to and including layer *j* . The part of *Qj* originating from scattering between the layers has to be transported to layer *j* by the appropriate Jacobian. Because of the cumulative effect of multiple scattering *γ i* and *γ k* are correlated. If *i<k*, the covariance is given by

$$\text{cov}(\boldsymbol{\mathcal{Y}}\_{l}, \boldsymbol{\mathcal{Y}}\_{k}) = H\_{l}(\boldsymbol{F}\_{l|1}\boldsymbol{\mathcal{Q}}\_{1}\boldsymbol{F}\_{k|1}^{\mathrm{T}} + \cdots + \boldsymbol{F}\_{l|l-1}\boldsymbol{\mathcal{Q}}\_{l-1}\boldsymbol{F}\_{k|l-1}^{\mathrm{T}} + \boldsymbol{\mathcal{Q}}\_{l}\boldsymbol{F}\_{k|l}^{\mathrm{T}})\boldsymbol{H}\_{k}^{\mathrm{T}}.$$

The observations *mk*, the functions *dk*, their Jacobians *Dk*, and the noise *γ k* are now collected in single vectors and a matrix:

$$m = \begin{pmatrix} m\_1 \\ \cdot \\ \cdot \\ m\_n \end{pmatrix}, \quad d = \begin{pmatrix} d\_1 \\ \cdot \\ \cdot \\ d\_n \end{pmatrix}, \quad D = \begin{pmatrix} D\_1 \\ \cdot \\ \cdot \\ \cdot \\ D\_n \end{pmatrix}, \quad \nu = \begin{pmatrix} \nu\_1 \\ \cdot \\ \cdot \\ \cdot \\ \mathbf{\varPsi}\_n \end{pmatrix}.$$

where *n* is the total number of measurement layers. This gives the following model:

$$
\mathfrak{m} = d(q\_0) + \mathfrak{y},
$$

which now can be linearized into

$$
\mathfrak{m} = Dq\_0 + \mathfrak{c} + \mathfrak{y},
$$

where *c* is a constant vector. The global least-squares estimate of *q*0 is given by

$$
\tilde{q}\_0 = (\boldsymbol{\mathcal{D}}^{\mathsf{T}} \boldsymbol{G} \boldsymbol{D})^{-1} \boldsymbol{D}^{\mathsf{T}} \boldsymbol{G} \, (\boldsymbol{\mathfrak{m}} - \boldsymbol{c}),
$$

where *V* = *G*−1 is the non-diagonal covariance matrix of *γ* . The quality of the initial expansion point can be monitored by using the obtained estimate as a new expansion point, and the state vector estimate can hence be re-calculated. Such a procedure is repeated until convergence, defined by a suitable stopping criterion.

If the track model is a circle and multiple scattering and energy loss can be neglected, the estimation can be simplified substantially. Explicit estimators are given in [52] for the center and radius of the circle, and in [53] for the curvature, the direction and the distance from a fixed point. Other algorithms are based on conformal mapping in the plane [2] or on a mapping to the Riemann sphere [54– 56].

If there is strong multiple scattering, the estimated track can be quite far away from the real track. In order to follow the actual track more closely, two projected scattering angles can be explicitly estimated at each detector layer or at a set of virtual breakpoints inside a continuous scatterer [45, 57]. The breakpoint method, also known as General Broken Lines [58], and the global least-squares method are equivalent, as far as the estimate of the state vector *q*0 is concerned [59].

If the number of measurements or the number of breakpoints is substantial, the computational cost of these methods can be high due to the necessity of inverting large matrices during the estimation procedure. The Kalman filter, a recursive formulation of the least-squares method, requires the inversion of only small matrices and exhibits the same attractive feature as the breakpoint method of following the actual track quite closely [26, 60].

As mentioned earlier, the Kalman filter proceeds by alternating prediction and update steps. The prediction step is the propagation of the track parameter vector from one detector layer containing a measurement to the next,

$$\mathfrak{q}\_{k|k-1} = f\_{k|k-1}(\mathfrak{q}\_{k-1|k-1}),$$

and the associated covariance matrix,

$$\mathbf{C}\_{k|k-1} = F\_{k|k-1} \mathbf{C}\_{k-1|k-1} F\_{k|k-1} \mathbf{T} + \mathbf{Q}\_k.$$

The update step is the correction of the predicted state vector due to the information from the measurement in layer *k*:

$$
\mathfrak{q}\_{k|k} = \mathfrak{q}\_{k|k-1} + \mathbf{K}\_k \left[ \mathfrak{m}\_k - \mathfrak{h}\_k(\mathfrak{q}\_{k|k-1}) \right],
$$

where the gain matrix *Kk* is given by

$$\mathbf{K}\_k = \mathbf{C}\_{k|k-1} \mathbf{H}\_k^\top \left(\mathbf{V}\_k + \mathbf{H}\_k \mathbf{C}\_{k|k-1} \mathbf{H}\_k^\top \right)^{-1}.$$

The update of the covariance matrix is given by

$$\mathbf{C}\_{k|k} = (I - K\_k H\_k) \mathbf{C}\_{k|k-1}.$$

The information filter is a mathematically equivalent, but numerically more stable formulation of the Kalman filter. In the information filter, the update of the state vector reads

$$\boldsymbol{\mathfrak{q}}\_{k|k} = \mathbf{C}\_{k|k} \left[ \left( \mathbf{C}\_{k|k-1} \right)^{-1} \boldsymbol{\mathfrak{q}}\_{k|k-1} + \boldsymbol{H}\_k^{\mathrm{T}} \boldsymbol{V}\_k^{-1} \boldsymbol{\mathfrak{m}}\_k \right],$$

whereas the update of the covariance matrix is given by

$$\mathcal{C}\_{k|k} = \left[ \left( \mathcal{C}\_{k|k-1} \right)^{-1} + \boldsymbol{H}\_k^\top \boldsymbol{V}\_k^{-1} \boldsymbol{H}\_k \right]^{-1}.$$

The implementation of the Kalman filter requires the computation of the Jacobians *Fk*|*k*−1 and *H k*. A compilation of analytical formulas for two important cases (fixedtarget configuration and solenoidal configuration) is given in [61].

Full information of the track parameters at the end of the track is obtained when all *n* measurements in the track candidate have been processed by the filter. The full information can be propagated back to all previous estimates by another iterative procedure, the Kalman smoother. A step of the smoother from layer *k* + 1 to layer

$$
\mathfrak{q}\_{k|n} = \mathfrak{q}\_{k|k} + A\_k (\mathfrak{q}\_{k+1|n} - \mathfrak{q}\_{k+1|k}),
$$

where the smoother gain matrix is given by

$$A\_k = \mathbf{C}\_{k|k} \mathbf{F}\_{k+1|k} ^\mathrm{T} (\mathbf{C}\_{k+1|k})^{-1} \mathrm{.}$$

The smoothed covariance matrix is

$$\mathbf{C}\_{k|n} = \mathbf{C}\_{k|k} - A\_k (\mathbf{C}\_{k+1|k} - \mathbf{C}\_{k+1|n}) A\_k^{\mathrm{T}}.$$

The smoother can also be realized by combining two filters running in opposite directions: a forward filter from *m*1 to *mn* and a backward filter from *mn* to *m*1. The smoothed states are the weighted mean of the predicted states of one filter and the updated states of the other filter. This approach is numerically more stable than the gain matrix formulation of the smoother.

#### **<IP_ADDRESS> Track Quality and Robust Estimation**

Robust estimators are insensitive to outliers, i.e., measurements that are biased or do not originate from the particle creating the majority of the hits in a track candidate. Some estimators are inherently robust by construction; other estimators can be made robust by finding and discarding outliers.

In the Kalman filter, the residual of the measurement in layer *k* with respect to the updated state vector is

$$r\_{k|k} = \mathfrak{m}\_k - \mathfrak{h}\_k(\mathfrak{q}\_{k|k}),$$

and the covariance matrix of this residual is

$$\mathcal{R}\_{k|k} = \mathcal{V}\_k - H\_k \mathcal{C}\_{k|k} H^\top.$$

The chi-square increment in layer *k* is

$$
\chi\_{k,+}^2 = \mathbf{r}\_{k|k}{}^\mathrm{T} \mathbf{R}\_{k|k}^{-1} \mathbf{r}\_{k|k},
$$

and the total chi-square of the track is found by summing up the chi-square increments for all measurements in the track candidate. The total chi-square is used to evaluate the quality of the track candidate. A too large value of this test statistic indicates that one or more of the measurements of the track candidate do not originate from the particle creating the majority of the measurements. Such measurements are called outliers.

An outlier rejection procedure can make use of the chi-squares of the measurements with respect to the smoothed predictions, i.e., a weighted mean of the predicted states of a forward and a backward Kalman filter. The measurement with the largest value of the chi-square is removed, and the total chi-square is again calculated. This procedure is repeated until the value of the total chi-square falls below a defined threshold.

In the presence of a potentially large fraction of outliers in a track candidate, the sequential outlier rejection procedure outlined above might become unstable, because the smoothed predictions may themselves be biased by outliers. An alternative approach is the Gaussian-sum filter [62]. This algorithm is based on the assumption that the probability distribution of the measurement error can be modeled as a two-component Gaussian mixture, where a narrow component represents the hypothesis that the measurement is real and a wider component represents the hypothesis that the measurement is an outlier. It takes the form of a set of Kalman filters running in parallel, each Kalman filter representing a specific hypothesis of a subset of the measurements that should be classified as outliers. A weight attached to each Kalman filter can be interpreted as the probability of correctness of the hypothesis. In the end, the Kalman filter with the largest weight or a weighted mean of the different filters can be taken as the final estimate.

The Gaussian-sum filter can also be used to deal with a mixture model of the process noise, i.e., the stochastic disturbance of the track because of interactions with the detector material [63]. In the case of bremsstrahlung, a successful application to the reconstruction of electrons is described in [50].

For the treatment of outliers, the Gaussian-sum filter has two disadvantages. First, it may create a large number of Kalman filters running in parallel because of poor knowledge of the track parameters in the early stages of the filter, making the approach expensive in terms of computing time. Second, an explicit outlier model is required. A faster and even more robust alternative is the Deterministic Annealing Filter [64]. This filter is an iterated Kalman filter with annealing, which assigns small weights to measurements far away from the track. A temperature parameter is introduced, facilitating convergence to the globally optimal solution. The iterations start at a high temperature, continue with a gradual lowering of the temperature and converge at the nominal value of the temperature. The procedure is easily generalized to the situation of several measurements being present in the same detector layer. In this case the measurements compete for inclusion in the track. As opposed to a standard outlier rejection approach, the assignment of measurements is soft. This means that several measurements in the same detector layer might contribute to the final estimate of the track parameters, each with a weight equal to the assignment probability. A further generalization is the multi-track filter, where several tracks are allowed to compete for compatible hits in all detector layers [65]. For an experimental application, see [66].

#### **<IP_ADDRESS> Jet Reconstruction**

Jets are bundles of collimated hadrons, reflecting hard scattering processes at the parton level. In order to carry out detailed comparisons between parton-level predictions and hadron-level observations a well-defined "jet finder" is required. In the jet finding information from both the tracking devices and the calorimeters is used.

Jet finding can be understood as finding clusters in the set of reconstructed tracks, including neutral tracks. As in the case of vertex finding (see Sect. 13.2.2), various types of clustering methods have been proposed and investigated. The performance strongly depends on the underlying physics, and usually a jet finder is optimized for specific physics requirements. For instance, the widely used *k*⊥ clustering algorithm comes in several versions, for instance one for *e*+*e*− collisions [67], and one for hadron-hadron collisions [68].

Hierarchical cluster algorithms offer a large variety of jet finders, differing mainly by the definition of the measure of distance between objects (tracks and jets), but sometimes also by the order in which the objects are combined. Some examples of agglomerative clustering algorithms are described and studied in [69]. Table 13.1 gives a summary of the distance measures used. The names refer to the ones used in [69]. *Ei* is the energy of cluster *i*, *pi* is its momentum, *θij* is the opening angle between the momentum vectors of the two clusters, and *E*vis is the visible energy.

A divisive hierarchical clustering algorithm is described in [76]. It is based on the following measure of distance between two tracks:

$$d\_{lj} = \frac{\theta\_{lj}^2}{p\_i p\_j},$$

**Table 13.1** Some distance measures used for agglomerative jet finding with respective references

but can be generalized to any other measure of distance. The method first constructs a minimum spanning tree [77] in the edge-weighted graph connecting all particles with each other and then proceeds to cut the tree along its longest edges. The procedure stops when the longest remaining edge is shorter than a fixed multiple of the median of all edge lengths.

Several non-hierarchical cluster algorithms have been proposed as well. Some of them employ general unsupervised learning methods, such as deterministic annealing [78] or *k*-means [79]. Others are specially designed for jet finding, for instance the cone algorithm described in [72]. It is an iterating procedure which constructs jets out of seeds. In contrast to the hierarchical clustering method the jets may overlap and a unique assignment has to be forced at the end. A modified cone algorithm suitable for the much larger multiplicity of heavy-ion collisions is proposed in [80]. A specialized jet finder for the reconstruction of hadronic *τ* -decays is described in [81].

### *13.1.4 Detector Alignment***1**

Alignment is the general term used in experimental high energy physics to refer to the process of obtaining and applying corrections to the nominal setup of a given experiment. These corrections are typically related to geometrical displacements of devices with a spatial resolution, in contrast to calibrations, where the corrections are usually extracted from pedestal or reference measurements to compensate for offsets in scalar measurements. Misalignment compromises tracking and vertex finding [82] and thus directly affects physics measurements such as momentum and invariant mass resolutions, or the efficiency of b-tagging algorithms. There are various possibilities for the treatment of alignment corrections, ranging from simple translations and rotations, equivalent to those of a rigid body, to more complex deformations, like sags or twists.

To this end experiments typically use several independent strategies [83]. For testing the long-term stability or the alignment of sub-detectors with respect to each other, very often so-called hardware alignment is utilized, where special reference markers are measured directly e.g. via optical systems or photogrammetry. However, these techniques reach only a limited precision in the range of several tens to hundreds of microns. If the intrinsic resolution of a tracking device is smaller, an improved resolution can only be obtained with track-based alignment, where the information from recorded particle tracks is used to obtain the alignment parameters [83, 84]. For various examples of the track-based alignment methods used in experiments since the LEP era, see [85–100].

1The section on detector alignment was contributed by E. Widl (Institute of High Energy Physics, Vienna; now at Austrian Institute of Technology).

#### **<IP_ADDRESS> General Overview**

The basis of all track-based alignment algorithms is an extended track model *d*, where the measurements *m* depend not only on the true track-parameters *q*0, but also on a set of alignment parameters *p*0 that describe the effects of sufficiently small deviations from the ideal geometry:

$$m = d(q\_0, p\_0) + \mathcal{Y}, \quad \text{cov}(\mathcal{Y}) = V.$$

The stochastic term *γ* , which describes the intrinsic resolution of the tracking devices and the effects of multiple scattering, is dealt with via its covariance matrix *V* . Since typically high momentum particles are used, energy-loss effects can be assumed to be deterministic and hence directly taken care of in the track model *d* itself.

With an initial guess *q*˘ for the track parameters and *p*˘ for the alignment parameters, this model allows to define residuals that are functions of the unknowns *q* and *p*:

$$r(q, p) = \mathfrak{m} - d(q, p) \approx \mathfrak{m} - d - D\_{\mathfrak{q}} \Delta q - D\_{\mathfrak{p}} \Delta p \tag{13.5}$$

with

$$\check{d} = d(\check{q}, \check{p}), \quad \Delta q = q - \check{q}, \quad \Delta p = p - \check{p},$$

$$D\_q = \left. \frac{\partial d}{\partial q} \right|\_{\check{q}, \check{p}}, \quad D\_p = \left. \frac{\partial d}{\partial p} \right|\_{\check{q}, \check{p}}.$$

The goal of a track-based alignment algorithm is to determine *p* from the residuals *r*, by minimizing the quadratic form *χ*2 = *r*T*V* −1*r*, using a sufficiently large set of recorded tracks. The methods used are quite diverse, but can be grouped into two categories: biased and unbiased algorithms.

Biased algorithms initially ignore the fact that the initial guess of the track parameters *q*0 is in general biased by the factual misalignment. In other words, by setting *q* = ˘*q* for every track, the residuals become a function of *p* alone, i.e., *r(q, p)* → *r(p)*. In general, the influence of the biased track information has to be compensated by iterating several times over the track sample, where at each iteration step the previously determined parameters are applied to the track reconstruction.

Unbiased algorithms on the other hand, minimize the residuals or the normalized residuals, respectively, estimating at the same time the track parameters. The problem with such an approach is the resulting huge number of parameters. In the presence of *N* alignment parameters and a sample of *M* tracks with *m* track parameters each, a total of *N* + *m* · *M* parameters have to be dealt with. While the value of *N* depends on the experimental setup, and *m* usually equals 5, the number of tracks *M* has always to be of considerable size to acquire reasonable statistics. On the other hand, unbiased algorithms usually do not require iterations, with the possible exception of problems like non-linearities or rejection of outliers.

Besides the differences between various algorithms it should be noted that the final result of any track-based alignment is always limited by the tracks used. Basic quality cuts, like the selection of high momentum tracks to minimize the influence of multiple scattering or cuts on the minimum number of hits, have a strong influence on the convergence. More subtle is the effect of an unbalanced mixture of tracks or the complete absence of some types of tracks, such as tracks from collisions and cosmic events or tracks taken with and without a magnetic field. This is due to the fact that any kind of tracks has several unconstrained degrees of freedom, usually referred to as weak modes, weakly defined modes or *χ*2-invariant modes. As an example, typical weak modes for straight tracks are shears but not bends, and vice versa for curved tracks. Combining the information of both kinds of tracks is therefore a reasonable strategy to avoid these deformations in the final result. The most obvious weak mode is a translation or rotation of the entire tracking device, which can be only fixed with some kind of reference frame, be it an external system or by definition. This, however, is less severe and sometimes even not considered at all, as it does not affect the internal alignment of the tracking device.

Once a set of alignment parameters is calculated, it should always be validated [83, chapter 11]. Apart from checking the improvement of the residuals, several physics measurements can be utilized, especially to probe for remaining weak modes. Known charge, forward-backward or *ϕ*-symmetries of distinct physics processes can be used. Distributions of the signed curvature or the signed transverse impact parameter are also sensitive observables.

#### **<IP_ADDRESS> Examples of Alignment Algorithms**

Some modern experiments deploy large tracking devices that require a large number of alignment parameters, of the order of 105. In such a case the computation of parameters by using straightforward recipes might become unreasonably slow or cause numerical problems. The two algorithms presented in this section are examples of how to cope with such challenging circumstances.

#### The HIP Algorithm

The HIP algorithm [101] is a straightforward and easy-to-implement biased alignment algorithm. It computes the alignment parameters for each alignable object separately. Only when iterating on the track sample a certain kind of indirect feedback between the alignable objects is established due to the track refit.

Since only individual alignable objects are regarded, Eq. (13.5) can be partitioned. This is simply done by evaluating the corresponding expressions for each alignable object *i* together with its associated parameters *pi*:

$$r\_l(\mathbf{p}\_l) = \mathbf{m}\_l - \mathbf{d}\_l(\check{\mathbf{q}}, \mathbf{p}\_l) \approx \mathbf{m}\_l - \mathbf{d}\_l - \mathbf{D}\_{pl}\Delta \mathbf{p}\_l$$

*.*

$$\check{\mathbf{d}}\_{l} = \mathbf{d}\_{l}(\check{\mathbf{q}}, \check{\mathbf{p}}\_{l}), \quad \Delta \mathbf{p}\_{l} = \mathbf{p}\_{l} - \check{\mathbf{p}}\_{l}, \quad \mathbf{D}\_{pl} = \frac{\partial \mathbf{d}\_{l}}{\partial \mathbf{p}\_{l}}\Big|\_{\check{\mathbf{q}}\_{l}, \check{\mathbf{p}}\_{l}}$$

The result is determined by minimizing the normalized squared residuals from a given set of tracks, again for each alignable object separately. The formal solution is given by

$$\Delta p\_{l} = \left(\sum\_{\text{tracks}} \mathbf{D}\_{pl} \prescript{\text{T}}{}{\mathbf{V}}\_{l}^{-1} \mathbf{D}\_{pl}\right)^{-1} \left(\sum\_{\text{tracks}} \mathbf{D}\_{pl} \prescript{\text{T}}{}{\mathbf{V}}\_{l}^{-1} \mathbf{r}\_{l}(\check{\mathbf{p}}\_{l})\right)^{-1}$$

The Millepede Algorithm

The Millepede algorithm [102] is an unbiased algorithm that minimizes the sum of the squared residuals of all tracks at once. To this end a system of linear equations, equivalent to the formal solution of an ordinary *χ*2-fit, is solved. However, to achieve this within a reasonable amount of time, only the solution for the alignment parameters is computed, while the computation of the improved track parameters is skipped. This is possible because of the special structure of the system: Firstly, the coefficient matrix is symmetric and, mostly due to the independence of the individual tracks, relatively sparse. Secondly, only the alignment parameters are common parameters for all track measurements, while the specific track parameters are only relevant for each corresponding track. Due to the latter, the solutions for the alignment and track parameters are only coupled via coefficient matrices of the form

$$\mathbf{G} = \mathbf{D}\_p \prescript{\mathbf{T}}{}{V}^{-1} \mathbf{D}\_q \prescript{}{}{}{\mathbf{D}}\_q \prescript{}{}{\mathbf{q}}$$

To set up the reduced system of equations, for each track the following information has to be extracted:

$$
\Gamma = \mathcal{D}\_q^{\,\,\,\,\,\,} V^{-1} \mathcal{D}\_q, \quad \mathcal{B} = \mathcal{D}\_q^{\,\,\,\,\,\,} V^{-1} \left( \mathfrak{m} - \check{\mathfrak{d}} - \mathcal{D}\_p \,\, \Delta \, p' \right) \,\, \, \,\,
$$

Here *p* = *p* − ˘*p* may already include an estimate *p* on the actual alignment. Then compute

$$\Delta \mathbf{C} = \mathbf{D}\_p^\top \mathbf{V}^{-1} \mathbf{D}\_p - \mathbf{G}^\top \mathbf{I}^{-1} \mathbf{G}^\top, \quad \Delta \mathbf{g} = \mathbf{D}\_p^\top \mathbf{V}^{-1} \left( \mathbf{m} - \check{\mathbf{d}} - \mathbf{D}\_p \Delta \mathbf{p}' + \mathbf{D}\_q \boldsymbol{\Gamma}^{-1} \boldsymbol{\mathfrak{f}} \right).$$

with

Note the expression −−1*β* instead of *q*. These are all necessary terms, including implicitly the full information from all track parameters. The complete system of equations to determine the alignment parameters then reads

$$
\mathbf{C} \,\Delta p = -\mathbf{g},
$$

with

$$\mathcal{C} = \sum\_{\text{tracks}} \mathbf{A} \mathcal{C}, \quad \mathbf{g} = \sum\_{\text{tracks}} \mathbf{A} \mathbf{g}.$$

The solution by matrix inversion is only feasible if the number of parameters is fairly small (*N* ≤ 103). The matrix *C* is usually relatively sparse, so that less time-consuming and more reliable methods can be used, such as the GMRES algorithm [103].

It is also possible to introduce constraints into the solution, which allows to align on various hierarchical levels at once. When aligning for instance on module- and layer-level at the same time, these constraints can remove redundant degrees of freedom by forcing the average movement of all modules within one layer to zero.

Millepede is a well-tested algorithm. To use it efficiently, some knowledge of its inner workings is of advantage. The HIP algorithm is simpler to implement, but less suitable for very large setups than Millepede. Another unbiased algorithm is the Kalman Alignment Algorithm [104]. It is a sequential method, derived from the Kalman filter (see also [105]).

### *13.1.5 Momentum Resolution*

The momentum resolution that can be achieved by a tracking detector is determined by the magnetic field, the arrangement and precision of the tracking detectors, and the amount of material crossed by the particle. Simple approximate formulas can be obtained for two cases:

The units are the same as in Sect. <IP_ADDRESS>: momentum in GeV/*c*, length in meters, and magnetic field in Tesla.

#### **<IP_ADDRESS> Two-arm Spectrometer**

We assume that the trajectory of the particle is parallel to the *z* axis and that *By* is the only significant component of the magnetic field. The angle of deflection is then given by [37]

$$\alpha \approx -\frac{kq}{p} \int\_{L} B\_{\text{y}} \,\mathrm{d}z = -\frac{kq}{p} \bar{B}\_{\text{y}} L\_{\text{y}}$$

where *L* is the length of the magnet, *B*¯*y* is the average value of the field along the trajectory, *p* is the momentum, and *q,k* are as in Eq. (13.2). Assuming that |*q*| = 1, linear error propagation gives

$$\frac{\sigma(p)}{p} = \frac{p\sigma(\alpha)}{k|\bar{B}\_{\underline{\chi}}|L}.$$

Assume that each arm consists of *m* identical position detectors spread over a length *l*, and that the standard deviation of the measurement error of *x* is equal to *δ*. The best angular resolution is obtained if in each arm half of the detectors is placed at each end of the arm. Neglecting all multiple scattering, it is equal to

$$
\sigma(\alpha) = \frac{2\delta}{l\sqrt{m/2}}.
$$

The relative momentum resolution due to measurement errors is therefore

$$\frac{\sigma\_{\rm me}(p)}{p} = \frac{2p\delta}{l\sqrt{m/2}k|\bar{B}\_\circ|L}.$$

Although this arrangement optimizes the precision in terms of geometry, it offers little redundancy for track finding and should be used only in setups with trivial pattern recognition requirements, for instance in the forward direction of fixed target experiments.

At low energies, multiple scattering can no longer be neglected. Whereas *σ (p)/p* arising from position measurement errors only is proportional to *p*, the term *σ*ms*(p)/p* arising from multiple scattering is proportional to 1*/(β*|*B*¯*y*|*L)*, which is large for small *β* and constant for high momenta (*β* ≈ 1). Under the same assumptions about the detector positions as above, the following formula is obtained:

$$\frac{\sigma\_{\rm ms}(p)}{p} = \frac{0.0136}{\beta k |\bar{B}\_{\rm y}| L} \left(\frac{md}{X\_0}\right)^{1/2}.$$

where *d/X*0 is the thickness of the detectors in units of radiation length. The total resolution is obtained by adding the corresponding variances and taking the square root,

$$\frac{\sigma(p)}{p} = \frac{\sigma\_{\rm me}(p)}{p} \oplus \frac{\sigma\_{\rm ms}(p)}{p} = ap \oplus b,$$

with *a* and *b* depending on the detector and the magnetic field.

#### **<IP_ADDRESS> Cylindrical Spectrometer**

Assume that there are *m* cylindrical detectors immersed in a homogeneous magnetic field *Bz* parallel to the *z* axis. The projection of the track on the *x*–*y* plane is a circle with curvature *κ*. For high momentum the circle can be approximated by a parabola, the detector cylinders can be approximated by planes, and multiple scattering can be neglected. For this case closed formulas for the joint covariance matrix of *κ* and the tangent *tϕ* = tan *ϕ* of the initial track direction *ϕ* can be given [106, 107]. For equidistant detectors, uniform resolution *δ* and *tϕ* = 0 it is given by

$$\operatorname{Cov} \begin{pmatrix} t\_{\varphi} \\ \kappa \end{pmatrix} = \frac{\delta^2}{m(m+1)(m+2)} \begin{pmatrix} \frac{12\left(m-1\right)(2m-1)(8m-11)}{L^2(m-2)} - \frac{360\left(m-1\right)^3}{L^3(m-2)} \\\\ -\frac{360\left(m-1\right)^3}{L^3(m-2)} & \frac{720\left(m-1\right)^3}{L^4(m-2)} \end{pmatrix},$$

where *L* is now the track length in the *x*–*y* projection. *L* is approximately equal to the radial distance between the innermost and the outermost detector. As *κ* = *kBz/p*T,

$$\frac{\sigma\_{\rm me}(p\_{\rm T})}{p\_{\rm T}} = \frac{p\_{\rm T}}{k|B\_{\rm \mathbb{Z}}|L} \frac{\delta}{L} \left[ \frac{720 \, (m-1)^3}{(m-2)m(m+1)(m+2)} \right]^{1/2}.$$

There is a high negative correlation between 1*/p*T and the direction tangent *tϕ*. For large *m*, the asymptotic values are

$$\begin{split} \frac{\sigma\_{\rm me}(p\_{\rm T})}{p\_{\rm T}} &= \frac{p\_{\rm T}}{k|B\_{\rm \bar{z}}|L} \frac{\delta}{L} \left[ \frac{720}{m+4} \right]^{1/2}, \ \sigma(t\_{\psi}) = \frac{\delta}{L} \left[ \frac{192}{m+3.875} \right]^{1/2}, \\ \text{cov}(t\_{\psi}, 1/p\_{\rm T}) &= -\frac{\sqrt{15}}{4} = -0.968. \end{split}$$

More general closed formulas for *tϕ* = 0 are given in [107].

*,*

**Table 13.2** Values of *Cm* in Eq. (13.6)

If one half of the detectors is placed at the center of the track and one quarter at either end, the variance of the curvature is minimal, and the covariance matrix reads

$$\operatorname{Cov} \begin{pmatrix} t\_{\phi} \\ \kappa \end{pmatrix} = \delta^2 \begin{pmatrix} \frac{72}{L^2 m} & -\frac{128}{L^3 m} \\ -\frac{128}{L^3 m} & \frac{256}{L^4 m} \end{pmatrix},$$

which is considerably smaller than in the equidistant case. This arrangement, however, is not particularly well suited for track finding and moreover difficult to realize.

The contribution of multiple scattering to the transverse momentum resolution can be approximated by

$$\frac{\sigma\_{\rm ms}(p\_{\rm T})}{p\_{\rm T}} = C\_m \cdot \frac{s}{\beta k \, |B\_\varepsilon| L} \left(\frac{md}{X\_0 \cos \lambda}\right)^{1/2},\tag{13.6}$$

*,*

where *d/X*0 is the thickness of the detectors in units of radiation length, *λ* = *π/*2−*θ* is the dip angle of the track, *s* = 0*.*0136*(*1 + 0*.*038 ln*(d/X*0*))*, *k* is as in Eq. (13.2), and *Cm* is a factor depending on *m*. Values of *Cm* for small *m*, obtained by the program described in [108], are given in Table 13.2. Note that the values are different from the ones given in [106]. In a time projection chamber *md/X*0 has to be replaced by *L/X*0, where *X*0 is the radiation length of the gas. The factor cos *λ* in the denominator accounts for the actual amount of matter traversed by a track with dip angle *λ*. Approximate formulas for the best possible resolution including multiple scattering can be found in [109].

The total transverse momentum resolution is calculated by quadratic addition,

$$\frac{\sigma(p\_{\rm T})}{p\_{\rm T}} = \frac{\sigma\_{\rm me}(p\_{\rm T})}{p\_{\rm T}} \oplus \frac{\sigma\_{\rm ms}(p\_{\rm T})}{p\_{\rm T}}$$

which can be written in the form

$$\frac{\sigma(p\_{\mathrm{T}})}{p\_{\mathrm{T}}} = \frac{a \, p\_{\mathrm{T}}}{\sqrt{m+4}} \oplus \frac{b\sqrt{m}}{\sqrt{\cos \lambda}}.$$

This shows that an optimal *m* exists for every *p*T and *λ* if the projected track length *L* is kept fixed. Overinstrumentation will deteriorate the resolution for low momenta unless additional measurements can be included without increasing the amount of matter to be traversed.

In order to calculate the error of the momentum *p* = *p*T*/* cos *λ* the error in *λ* must be taken into account:

$$
\sigma^2(p) = \sigma^2(p\_\Gamma) / \cos^2 \lambda + \sigma^2(\lambda) \, p\_\Gamma^{-2} \sin^2 \lambda / \cos^4 \lambda,
$$

the correlation between *p*T and *λ* being negligible in practice. Because of *σ (p)/p* = *pσ (*1*/p)* it follows that:

$$\frac{\sigma(p)}{p} = \frac{\sigma(p\_{\rm T})}{p\_{\rm T}} \oplus \sigma(\lambda)\tan\lambda...$$

With the exception of very low momenta the track can be approximated by a straight line in the *r*–*z* projection, where *r* = *(x*2 + *y*2*)*1*/*2. For *m* equidistant detectors and uniform resolution *δ*, the variance of the direction tangent *tλ* = tan *λ* due to the measurement errors is given by [106]:

$$
\sigma\_{\rm me}^2(t\_\lambda) = \frac{\delta^2}{L^2} \frac{12(m-1)}{m(m+1)} \frac{1}{\cos^4 \lambda}.
$$

If the measurement error in *z* is very small, the variance of *tλ* is dominated by multiple scattering. For equidistant layers of uniform thickness *d* an approximate formula can be given. Remarkably, it does not depend on the number of layers:

$$
\sigma\_{\rm ms}^2(t\_\lambda) \approx \frac{s^2}{p^2} \frac{d}{X\_0 \cos \lambda} \frac{1}{\cos^4 \lambda} = \frac{s^2}{p\Gamma^2} \frac{d}{X\_0 \cos \lambda} \frac{1}{\cos^2 \lambda},
$$

with *s* = 0*.*0136*(*1 + 0*.*038 ln*(d/X*0*))*.

In the design and optimization phase of the detector a precise evaluation of the resolution of all track parameters is mandatory. There are several software packages that allow a fast track simulation plus reconstruction in a general detector setup, for instance [110] (in FORTRAN), [111] (in Matlab/Octave), or [108] (in Java).

### **13.2 Vertex Reconstruction**

### *13.2.1 Introduction*

Vertex reconstruction is the task of finding and estimating the production point of a set of particles. The pattern recognition algorithms and statistical estimation methods involved are in many respects similar to the ones used in track reconstruction. For an overview of vertex reconstruction algorithms used in past or active experiments see for instance [112–116].

In practice it is useful to distinguish between several types of vertices:

Vertex reconstruction frequently proceeds in several steps:

Vertex finding can be accomplished in many different ways. A few of them will be described in Sect. 13.2.2. The vertex fit takes a vertex candidate and estimates the vertex location from the estimated track parameters of the outgoing particles (Sect. 13.2.3). As a rule, only charged particles are used, but sometimes also neutral particles contribute to the vertex fit. In the test stage (Sect. <IP_ADDRESS>) outliers are identified, i.e., particles that apparently do not belong to the estimated vertex. As this can lead to a different assignment of particles to vertices, it can be considered as a method of vertex finding. If outliers are expected, the estimation procedure should be robust so that the estimated vertex is not significantly biased by the outliers (Sect. <IP_ADDRESS>). Kinematic constraints (Sect. 13.2.4) are usually imposed via Lagrange multipliers. By repeating the kinematic fit under various mass hypotheses of the mother and/or daughter particles the most likely mass assignment can be found out.

### *13.2.2 Vertex Finding*

Vertex finding is the process of dividing the reconstructed tracks in an event into classes such that presumably all tracks in a class are produced at the same vertex. The primary vertex in an event is usually easy to find, especially if prior information about its location is available (beam profile, target position). On the other hand, secondary decay vertices of short-lived decays are hard to find, as some of the decay products may also be compatible with the primary vertex. Vertex finding methods can be roughly divided in three main types: generic clustering algorithms, topological methods, and iterated estimators. The latter can be considered as a special divisive clustering method.

#### **<IP_ADDRESS> Clustering Methods**

As mentioned above in the context of jet finding (see Sect. <IP_ADDRESS>), clustering methods are based on a distance matrix or a similarity matrix of the objects to be classified. A cluster is then a group with small distances (large similarities) inside the group and large distances (small similarities) to objects outside the group. The distance measure reflects only the geometry of the tracks.

Various clustering methods have been evaluated in the context of vertex finding, of both the hierarchical and the non-hierarchical type [117]. Hierarchical clustering can be agglomerative or divisive. In agglomerative clustering each track starts out as a single cluster. Clusters are merged iteratively on the basis of a distance measure. The shortest distance in space between two tracks is peculiar insofar as it does not satisfy the triangle inequality: if tracks *a* and *b* are close, and tracks *b* and *c* are close, it does not follow that tracks *a* and *c* are close as well. The distance between two clusters of tracks should therefore be defined as the maximum of the individual pairwise distances, known as complete linkage in the clustering literature. Alternatively, the distance between two clusters can be the distance between the two vertices fitted from the clusters. Divisive clustering starts out with a single cluster containing all tracks. Further division of this cluster can be based on repeated vertex estimation with outlier identification (see Sect. <IP_ADDRESS>). Examples of nonhierarchical clustering methods used in vertex finding are vector quantization, the *k*-means algorithm and deterministic annealing [113].

#### **<IP_ADDRESS> Topological Methods**

A very general topological vertex finder was proposed in [118]. It is related to the Radon transform, which is a continuous version of the Hough transform used for track finding (Sect. <IP_ADDRESS>). The search for vertices is based on a function *V (v)* which quantifies the probability of a vertex at location *v*. For each track a Gaussian probability tube *fi(v)* is constructed. The function *V (v)* is defined taking into account that the value of *fi(v)* must be significant for at least two tracks:

$$V(\boldsymbol{\upsilon}) = \sum\_{i=0}^{n} f\_i(\boldsymbol{\upsilon}) - \frac{\sum\_{i=0}^{n} f\_i^2(\boldsymbol{\upsilon})}{\sum\_{i=0}^{n} f\_i(\boldsymbol{\upsilon})}$$

Due to the second term on the right-hand side, *V (v)* ≈ 0 in regions where *fi(v)* is significant for only one track. The form of *V (v)* can be modified to fold in known physics information about probable vertex locations. For instance, *V (v)* can be augmented by a further function *f*0*(v)* describing the location and spread of the interaction point. In addition, *V (v)* may be modified by a factor dependent on the angular location of the point *v*.

Vertex finding amounts to finding local maxima of the function *V (v)*. The search starts at the calculated maxima of the products *fi(v)fj (v)* for all track pairs. For each of these points the nearest maximum of *V (v)* is found. These maxima are clustered together to form candidate vertex regions. The final association of the tracks to the vertex candidates can be done on the basis of the respective *χ*2 contributions or by an adaptive fit (see Sect. <IP_ADDRESS>). In [119] the topological vertex finder was augmented by a procedure based on the concept of the minimum spanning tree of a graph.

#### **<IP_ADDRESS> Iterated Estimators**

Vertex finding can also be accomplished by iterated vertex fits (see Sect. 13.2.3). The procedure can be summarized in the following way:

The iteration stops when no vertex with at least two tracks can be successfully fitted. Step 2 might itself be iterative, especially if the vertex fit is not robust, so that the incompatible tracks have to be removed sequentially. Iterative vertex finders based on a least-squares fit (Sect. <IP_ADDRESS>) and an adaptive fit (Sect. <IP_ADDRESS>) are implemented in the RAVE toolbox [120, 121].

### *13.2.3 Vertex Fitting*

The input to the vertex fit is a vertex candidate, i.e., a set of estimated track parameters { ˜*q*1*,..., q*˜ *n*} located at one or more reference surfaces, along with their covariance matrices {*C*1*,..., Cn*}. For instance, in the primary vertex fit in a collider experiment the reference surface may be the beam tube. If possible, the reference surface(s) should be chosen such that multiple scattering between the vertex and the location of the track parameters is negligible.

The parameters to be fitted are the vertex position *v* and the track momenta *pi* at the vertex. The functional dependence of the track parameters on the vertex parameters requires a track model, which depends on the shape of the magnetic field in the vicinity of the vertex. If the field is homogeneous, the track model is a helix; if the field is zero, the track model is a straight line. In other cases the track model may have to be computed numerically (see Sect. <IP_ADDRESS>).

#### **<IP_ADDRESS> Least-Squares Methods**

The conventional approach to estimating the vertex position is the minimization of some quadratic objective function, yielding a least-squares estimate. There are two main flavors of least-squares estimation in vertex fitting, constrained and unconstrained minimization. In the first case the vertex constraint is introduced into the objective function via a Lagrange multiplier, in the second case the constraint is implicit in the track model.

As an example, consider a vertex fit with *n* straight tracks. The *n* straight tracks originating from the common vertex *v* = *(x*v*, y*v*, z*v*)*T can be represented by *n* straight lines with parameters *λi*:

$$\mathbf{x} = \mathbf{x}\_{\mathrm{V}} + \lambda\_{l} a\_{l}, \quad \mathbf{y} = \mathbf{y}\_{\mathrm{V}} + \lambda\_{l} b\_{l}, \quad z = z\_{\mathrm{V}} + \lambda\_{l}, \quad i = 1, \ldots, n, \mathbf{y}$$

where *ai* and *bi* are the direction tangents at the vertex. At the reference surface *z* = *z*ref track *i* is specified by its parameter vector *qi* = *(xi, yi, ai, bi)*T, consisting of the intersection point *(xi, yi)* and the two direction tangents *(ai, bi)*. The track fit delivers estimates *q*˜*i* and information matrices *Gi* for *i* = 1*,...,n*. In the constrained problem the sum of the squared residuals

$$M(\mathbf{q}\_1, \dots, \mathbf{q}\_n) = \sum\_{l=1}^n \mathbf{e}\_l^\top \mathbf{G}\_l \mathbf{e}\_l, \quad \mathbf{e}\_l = \tilde{\mathbf{q}}\_l - \mathbf{q}\_l,\tag{13.7}$$

must be minimized under the 2*n* nonlinear constraints

$$\mathbf{x}\_{\rm V} = \mathbf{x}\_{l} + a\_{l}(\mathbf{z} - z\_{\rm ref}), \quad \mathbf{y}\_{\rm V} = \mathbf{y}\_{l} + b\_{l}(\mathbf{z} - z\_{\rm ref}), \quad i = 1, \ldots, n$$

There are 4*n* + 3 unknowns, 4*n* observations and 2*n* constraints, giving 4*n* + 2*n* − *(*4*n* + 3*)* = 2*n* − 3 degrees of freedom. The resulting track parameters *q*¯*i* fit best, in the least-squares sense, to the track fit estimates *q*˜*i* and at the same time have a common vertex. For the solution of the constrained vertex fit see Sect. <IP_ADDRESS>.

In the example, the constraints can be rewritten as

$$\mathbf{x}\_{l} = \mathbf{x}\_{\mathrm{V}} + (z\_{\mathrm{ref}} - z)a\_{l}, \quad \mathbf{y}\_{l} = \mathbf{y}\_{\mathrm{V}} + (z\_{\mathrm{ref}} - z)b\_{l}, \quad i = 1, \ldots, n. \tag{13.8}$$

Insertion of Eq. (13.8) into Eq. (13.7) gives the objective function of the unconstrained nonlinear least-squares problem:

$$M(\mathfrak{v}, a\_{\mathbb{I}}, b\_{\mathbb{I}}, \dots, a\_n, b\_n) = \sum\_{i=1}^n \mathbf{e}\_i^{\mathbb{T}} \mathbf{G}\_i \mathbf{e}\_i, \quad \mathbf{e}\_i = \tilde{\mathbf{q}}\_i - \mathbf{q}\_i.$$

There are now 4*n* observations and 2*n* + 3 unknown parameters, namely the vertex position and the track directions at the vertex, giving again 4*n* − *(*2*n* + 3*)* = 2*n* − 3 degrees of freedom.

A generalization of this simple case to helix tracks can be found in [122–124]. In the general case, the unconstrained problem can be formulated in terms of the unknown vertex position *v* and the unknown track momentum vectors *pi* at the vertex [26, 59]. The measurement equation reads

$$\mathbf{q}\_{l} = \hbar(\mathbf{v}, \mathbf{p}\_{l}), i = 1, \ldots, n,\tag{13.9}$$

where the function *hi* incorporates the track model in the magnetic field. The objective function is equal to

$$M(\mathfrak{v}, \mathfrak{p}\_1, \dots, \mathfrak{p}\_n) = \sum\_{i=1}^n \mathfrak{e}\_i^\mathrm{T} \mathbf{G}\_i \mathbf{e}\_i, \quad \mathfrak{e}\_i = \tilde{\mathfrak{q}}\_i - \mathfrak{q}\_i.$$

Minimization of the objective function can proceed in several ways. For a detailed exposition of non-linear least-squares estimation see e.g. [125].

#### Gauss-Newton Method

Assume that there are approximate values *v*˘ and *p*˘ *i* for all *i*. Then Eq. (13.9) can be approximated by an affine function:

$$\mathbf{q}\_{l} \approx \mathbf{h}\_{l}(\check{\mathbf{v}}, \check{\mathbf{p}}\_{l}) + A\_{l}(\mathbf{v} - \check{\mathbf{v}}) + \mathbf{B}\_{l}(\mathbf{p}\_{l} - \check{\mathbf{p}}\_{l}) = \mathbf{c}\_{l} + A\_{l}\mathbf{v} + \mathbf{B}\_{l}\mathbf{p}\_{l},$$

with

$$A\_{l} = \frac{\partial h\_{l}(\mathbf{v}, \mathbf{p}\_{l})}{\partial \mathbf{v}}\Big|\_{\breve{\mathbf{v}}, \breve{\mathbf{p}}\_{l}}, \quad \mathbf{B}\_{l} = \left.\frac{\partial h\_{l}(\mathbf{v}, \mathbf{p}\_{l})}{\partial \mathbf{p}\_{l}}\right|\_{\breve{\mathbf{v}}, \breve{\mathbf{p}}\_{l}}, \quad \mathbf{c}\_{l} = h\_{l}(\breve{\mathbf{v}}, \breve{\mathbf{p}}\_{l}) - A\_{l}\breve{\mathbf{v}} - \mathbf{B}\_{l}\breve{\mathbf{p}}\_{l}.$$

The objective function then reads

$$M(\mathfrak{v}, \mathfrak{p}\_1, \dots, \mathfrak{p}\_n) = \sum\_{l=1}^n (\tilde{\mathfrak{q}}\_l - \mathfrak{c}\_l - A\_l \mathfrak{v} - \mathfrak{B}\_l \mathfrak{p}\_l)^\mathrm{T} \mathbf{G}\_l (\tilde{\mathfrak{q}}\_l - \mathfrak{c}\_l - A\_l \mathfrak{v} - \mathfrak{B}\_l \mathfrak{p}\_l) .$$

As *M* is now quadratic in the unknown parameters, the minimum can be computed explicitly. The estimated vertex position and its covariance matrix are given by

$$\tilde{\boldsymbol{\nu}}\_{n} = \mathbf{C}\_{n} \sum\_{l=1}^{n} \mathbf{A}\_{l}^{\top} \mathbf{G}\_{l}^{\,{}^{B}} (\tilde{\mathbf{q}}\_{l} - \mathbf{c}\_{l}), \quad \text{Var}(\tilde{\boldsymbol{\nu}}\_{n}) = \mathbf{C}\_{n} = \left(\sum\_{l=1}^{n} \mathbf{A}\_{l}^{\top} \mathbf{G}\_{l}^{\,{}^{B}} \mathbf{A}\_{l}\right)^{-1}, \tag{13.10}$$

#### 13 Pattern Recognition and Reconstruction 589

with

$$\mathbf{G}\_{l}^{\;B} = \mathbf{G}\_{l} - \mathbf{G}\_{l}\mathbf{B}\_{l}\mathbf{W}\_{l}\mathbf{B}\_{l}^{\;T}\mathbf{G}\_{l}, \quad \mathbf{W}\_{l} = (\mathbf{B}\_{l}^{\;T}\mathbf{G}\_{l}\mathbf{B}\_{l})^{-1}.$$

In general, the procedure has to be iterated. The measurement equation is expanded at the new estimate, and the estimate is recomputed until convergence is obtained. The formulas required for the implementation of two important cases, fixed-target configuration and solenoidal configuration, are given in [61].

Once *v*˜*n* is known, the track momenta and the full covariance matrix can be computed:

$$\begin{aligned} \tilde{p}\_{l}^{n} &= \mathbf{W}\_{l} \mathbf{B}\_{l}^{\top} \mathbf{G}\_{l} (\tilde{\mathbf{q}}\_{l} - \mathbf{c}\_{l} - A\_{l} \tilde{\mathbf{v}}\_{n}), \\ \text{Var}(\tilde{p}\_{l}^{n}) &= \mathbf{D}\_{l}^{n} = \mathbf{W}\_{l} + \mathbf{W}\_{l} \mathbf{B}\_{l}^{\top} \mathbf{G}\_{l} A\_{l} \mathbf{C}\_{n} A\_{l}^{\top} \mathbf{G}\_{l} \mathbf{B}\_{l} \mathbf{W}\_{l}, \\ \text{Cov}(\tilde{\mathbf{p}}\_{l}^{n}, \tilde{\mathbf{v}}\_{n}) &= \mathbf{E}\_{l}^{n} = -\mathbf{W}\_{l} \mathbf{B}\_{l}^{\top} \mathbf{G}\_{l} A\_{l} \mathbf{C}\_{n}. \end{aligned} \tag{13.11}$$

The estimates can also be computed recursively, resulting in an extended Kalman filter [25, 26, 59]:

$$\begin{aligned} \tilde{\boldsymbol{\nu}}\_{i} &= \mathbf{C}\_{i}[\mathbf{C}\_{i-1}^{-1}\tilde{\boldsymbol{\nu}}\_{i-1} + \mathbf{A}\_{i}^{\mathrm{T}}\mathbf{G}\_{i}^{\mathrm{B}}(\tilde{\mathbf{q}}\_{i} - \mathbf{c}\_{i})], & \mathbf{C}\_{i} &= (\mathbf{C}\_{i-1}^{-1} + \mathbf{A}\_{i}^{\mathrm{T}}\mathbf{G}\_{i}\mathbf{A}\_{i})^{-1} \\ \tilde{\boldsymbol{p}}\_{i} &= \mathbf{W}\_{i}\mathbf{B}\_{i}^{\mathrm{T}}\mathbf{G}\_{i}(\tilde{\mathbf{q}}\_{i} - \mathbf{c}\_{i} - \mathbf{A}\_{i}\tilde{\mathbf{v}}\_{i}), & \mathbf{D}\_{i} &= \mathbf{W}\_{i} + \mathbf{W}\_{i}\mathbf{B}\_{i}^{\mathrm{T}}\mathbf{G}\_{i}\mathbf{A}\_{i}\mathbf{C}\_{i}\mathbf{A}\_{i}^{\mathrm{T}}\mathbf{G}\_{i}\mathbf{B}\_{i}\mathbf{W}\_{i}, \\ \tilde{\mathbf{E}}\_{i} &= -\mathbf{W}\_{i}\mathbf{B}\_{i}^{\mathrm{T}}\mathbf{G}\_{i}\mathbf{A}\_{i}\mathbf{C}\_{i}. \end{aligned}$$

The associated smoother is tantamount to recomputing the track momenta using the last vertex estimate *v*˜*n*, i.e., Eq. (13.11).

#### Newton–Raphson Method

This method uses a local quadratic approximation to the objective function. In order to simplify the notation we introduce *α* = *(v, p*1*,..., pn)*T, *q*˜ = *(q*˜ 1*,..., q*˜ *n)*T and *h* = *(h*1*,..., hn)*T. Then the objective function can be written as

$$M(\mathfrak{a}) = \left[\tilde{q} - h(\mathfrak{a})\right]^\mathrm{T} \mathbf{G}[\tilde{q} - h(\mathfrak{a})], \quad \mathbf{G} = \mathrm{diag}(\mathbf{G}\_1, \dots, \mathbf{G}\_n).$$

If *α*˘ is an appropriate expansion point, *M(α)* is approximated by

$$M(\mathfrak{a}) \approx M(\check{\mathfrak{a}}) + \mathfrak{g}^{\mathrm{T}}(\mathfrak{a} - \check{\mathfrak{a}}) + \frac{1}{2}(\mathfrak{a} - \check{\mathfrak{a}})^{\mathrm{T}}\mathfrak{Q}(\mathfrak{a} - \check{\mathfrak{a}}),$$

where

$$\mathbf{g} = \frac{\partial M}{\partial \mathbf{a}} = -2H^T G [\tilde{\mathbf{q}} - \hbar(\breve{\mathbf{a}})], \quad \mathbf{\mathfrak{Q}} = \frac{\partial^2 M}{\partial \mathbf{a} \partial \mathbf{a}^T} = 2H^T G H - 2 \frac{\partial H^T}{\partial \mathbf{a}^T} \mathbf{G} [\tilde{\mathbf{q}} - \hbar(\breve{\mathbf{a}})]$$

are the gradient and the Hessian of *M*, respectively, evaluated at *α*˘, and *H* is the Jacobian of the track model *h(α)*. If  is positive definite, *M* has a minimum when its gradient is zero, leading to

$$
\ddot{\mathfrak{a}} = \check{\mathfrak{a}} - \mathfrak{Q}^{-1} \mathfrak{g} \cdot \mathfrak{g}
$$

If the second term of the Hessian is set to zero, the Gauss–Newton method is recovered. Clearly, the Newton–Raphson method is more complex, but it gives some additional information about the problem. In particular, a Hessian that is not positive definite indicates that the expansion point is too far from the true global minimum.

#### Levenberg–Marquardt Method

In this method the matrix *H*T*GH* is inflated by a diagonal matrix *kI* . As a consequence, the direction of the parameter update is intermediate between the direction of the Gauss–Newton step (*k* = 0) and the direction of steepest descent (*k*→∞). An example of a vertex fit with the Levenberg–Marquardt method is given in [122].

#### Fast Vertex Fits

The estimated track parameters *q*˜*i* are frequently given at the innermost detector surface or at the beam tube. If the *q*˜*i* are propagated to the vicinity of the presumed vertex, the vertex estimation can be speeded up by applying some approximations.

The "perigee" parametrization for helical tracks was introduced in [123], with a correction in [124]. The track is parameterized around the point of closest approach (the perigee point *v*P) of the helix to the *z*-axis. The variation of transverse errors along the track is neglected in the vicinity of the perigee, and the track direction and curvature at the vertex is considered to be constant. The approximate objective function of the vertex fit can then be written entirely in terms of the perigee points:

$$M(\boldsymbol{\upsilon}) = \sum\_{i=1}^{n} (\boldsymbol{\upsilon}\_{i}^{\text{P}} - \boldsymbol{\upsilon})^{\text{T}} T\_{i} (\boldsymbol{\upsilon}\_{i}^{\text{P}} - \boldsymbol{\upsilon}),\tag{13.12}$$

where *T i* is a weight matrix of rank 2. The vertex estimate is then

$$
\tilde{\boldsymbol{v}} = \left(\sum\_{i=1}^n \boldsymbol{\mathcal{T}}\_i\right)^{-1} \left(\sum\_{i=1}^n \boldsymbol{\mathcal{T}}\_i \boldsymbol{\mathfrak{v}}\_i^{\boldsymbol{\mathcal{P}}}\right).
$$

The Jacobians required to compute the *T i* are spelled out in [123, 124].

A further simplification was proposed in [126]. In the vicinity of the vertex the track is approximated by a straight line. The estimated track parameters are transformed to a coordinate system the *x*-axis of which is parallel to the track. The vertex is then estimated by minimizing the sum of the weighted transverse distances of the tracks to the vertex. The resulting objective function has the same form as in Eq. (13.12), again with weight matrices of rank 2. The estimate is exact for straight tracks.

A different type of a fast vertex fitting algorithm is described in [127]. It is based on approximating the tracks by straight lines both in the *x*–*y* plane and in the *x*–*z* plane. In either projection, the lines representing the tracks are Houghtransformed to points in the dual plane of line parameters. The vertex coordinates are then obtained by a weighted linear least-squares fit in the dual plane.

#### Adding Prior Information

If the vertex to be fitted is the primary vertex, there may be prior information about the vertex position from the beam profile in a collider experiment or the target location in a fixed target experiment. The prior information usually comes in the form of a position *v*0 plus a covariance matrix *C*0. The objective function is then augmented by an additional term

$$(\boldsymbol{\upsilon}\_0 - \boldsymbol{\upsilon})^\mathsf{T} \boldsymbol{C}\_0^{-1} (\boldsymbol{\upsilon}\_0 - \boldsymbol{\upsilon}) .$$

For instance, the Gauss–Newton estimate Eq. (13.10) is modified in the following way:

$$\tilde{\boldsymbol{\nu}}\_{n} = \mathbf{C}\_{n} \left[ \mathbf{C}\_{0}^{-1} \boldsymbol{\nu}\_{0} + \sum\_{i=1}^{n} \mathbf{A}\_{i}^{\top} \mathbf{G}\_{i}^{\,{}^{B}} (\tilde{\boldsymbol{q}}\_{i} - \mathbf{c}\_{i}) \right], \\ \text{Var}(\tilde{\boldsymbol{\nu}}\_{n}) = \mathbf{C}\_{n} = \left( \mathbf{C}\_{0}^{-1} + \sum\_{i=1}^{n} \mathbf{A}\_{i}^{\top} \mathbf{G}\_{i}^{\,{}^{B}} \mathbf{A}\_{i} \right)^{-1}.$$

Similar modifications apply to the Newton–Raphson estimate and the fast vertex fits.

#### **<IP_ADDRESS> Vertex Quality and Outlier Removal**

Some tracks used in the vertex fit may be outliers in the sense that they do not actually belong to the vertex. Also, the estimated track parameters may be distorted by outliers or distorted hits in the track fit. Both types of outliers distort the vertex estimate and need to be identified.

In the case of Gaussian errors and a linear model the contribution of each track to the minimum value of the objective function is distributed according to a *χ*2 distribution with two degrees of freedom. The contribution *χ*2 *i* of track *i* has to be computed relative to the vertex estimated without track *i*. For instance, in the Gauss–Newton algorithm:

$$\chi\_i^2 = \mathbf{r}\_i^{n\mathbf{T}} \mathbf{G}\_l \mathbf{r}\_i^n + (\tilde{\boldsymbol{\nu}}\_n - \tilde{\boldsymbol{\nu}}\_n^{-i})^\mathbf{T} (\mathbf{C}\_n^{-i})^{-1} (\tilde{\boldsymbol{\nu}}\_n - \tilde{\boldsymbol{\nu}}\_n^{-i}),$$

where *rn i* = ˜*qi* − *ci* − *Aiv*˜*n* − *Bip*˜ *n i* is the residual of track *i* and *v*˜ −*i n* is the vertex estimate with track *i* removed:

$$
\tilde{\boldsymbol{\sigma}}\_{n}^{-l} = \mathbf{C}\_{n}^{-l} \left[ \mathbf{C}\_{n}^{-1} \tilde{\mathbf{v}}\_{n} - \mathbf{A}\_{l}^{\mathrm{T}} \mathbf{G}\_{l}^{\mathrm{B}} (\tilde{\mathbf{q}}\_{l} - \mathbf{c}\_{l}) \right], \quad \mathbf{C}\_{n}^{-l} = \left( \mathbf{C}\_{n}^{-1} - \mathbf{A}\_{l}^{\mathrm{T}} \mathbf{G}\_{l}^{\mathrm{B}} \mathbf{A}\_{l} \right)^{-1}.
$$

Analogous but somewhat simpler formulas hold for the fast vertex fits.

The test statistic *χ*2 *i* can be computed for all *i*, and the track with the largest *χ*2 *i* is a candidate for removal. This procedure can be repeated until all *χ*2 *i* are below the cut. Even if there is only a single outlier, all *χ*2 *i* are no longer *χ*2-distributed and the power of the test is impaired. This loss of power can be compensated by robust estimation of the vertex.

#### **<IP_ADDRESS> Robust and Adaptive Estimators**

Robust estimators are less influenced or not influenced at all by outlying observations. This can be achieved by downweighting outliers or by excluding them from the estimate. For example, in the case of a one-dimensional location estimate, the Mestimator [128] downweights outliers, whereas the LMS (least median of squares) estimator [129] uses only one half of the sample (the one spanning the shortest interval) and ignores the other one.

Robust estimators tend to be statistically less efficient and computationally more expensive than least-squares estimators. On the other hand, estimation and outlier detection are performed in parallel, whereas a least-squares estimator has to be recomputed after an outlier has been identified and removed.

One of the earliest proposals for a robust vertex fit is in [130]. The method is an M-estimator with Huber's *ψ*-function [131]. It is implemented as a re-weighted least-squares estimator. The initial vertex estimate is a plain least-squares estimate. Then, for each track, the residuals are rotated to the eigensystem of the covariance matrix of the track, and weight factors are computed according to

$$w\_{l} = \frac{\psi(r\_{l}/\sigma\_{l})}{r\_{l}/\sigma\_{l}} = \begin{cases} 1, & |r\_{l}| \le c\sigma\_{l}, \\ c\sigma\_{l}/|r\_{l}|, & |r\_{l}| > c\sigma\_{l}, \end{cases}$$

where *ri* is one of the residuals in the rotated frame, *σi* is the standard deviation in the rotated frame, and *c* is the robustness constant, usually chosen between 1 and 3. The weight factors are applied and the estimate is recomputed. The entire procedure is iterated until convergence.

A different kind of re-weighted least-squares estimator is proposed in [132]. The weights are computed according to Tukey's bi-square function [128]:

$$w\_{l} = \begin{cases} \left(1 - \frac{r\_{l}^{2}/\sigma\_{l}^{2}}{c^{2}}\right)^{2}, & |r\_{l}| \le c\sigma\_{l}, \\ 0, & \text{otherwise}, \end{cases}$$

where *r*2 *i* is the squared residual of track *i* with respect to the vertex, *σ*2 *i* is its variance, and *c* is again the robustness constant. The estimator is now equivalent to a redescending M-estimator, and consequently less sensitive to outliers than Huber's M-estimator.

The combination of a redescending M-estimator with the concept of deterministic annealing [133] leads to the adaptive method of vertex fitting [113, 117, 134, 135]. The concept of the adaptive vertex fit is derived from the Deterministic Annealing Filter [64] (see Sect. <IP_ADDRESS>). The weights are computed according to

$$w\_l = \frac{\exp\left(-\chi\_l^2/2T\right)}{\exp\left(-\chi\_l^2/2T\right) + \exp\left(-\chi\_{\text{cut}}^2/2T\right)},$$

where *χ*2 *i* is the *χ*2-contribution of track *i*, *χ*2 cut is a cutoff value, and *T* is a temperature parameter. The computation of the redescending M-estimator can be interpreted as an EM (expectation–maximization) algorithm [136, 137]. Alternatively it can be viewed as the minimization of the energy function of an elastic arm algorithm [18, 19]. If annealing is employed, the iteration starts at high *T* . The temperature is then gradually decreased. At low *T* the weights approach either zero or one. The final weights can be used for classification of the tracks as inliers or outliers. A comparison of the adaptive method with other robust estimators can be found in [138]. The adaptive estimator has been extended to a multi-vertex estimator fitting several vertices simultaneously, including competition of all vertices for all tracks [139].

Iterated re-weighted least-squares estimators require a good starting point, in order to ensure convergence to the global minimum and to minimize the number of iterations required. In many cases a standard least-squares estimate is sufficient. In the presence of a large number of outliers also the starting point should be estimated robustly, preferably by an estimator with a high breakdown point [129]. Several such initial estimators have been proposed and studied in [117].

The M-estimators and the adaptive estimator presented above do not presuppose an explicit outlier model. If it is possible to describe the outliers by a Gaussian mixture model, estimation of the vertex can be carried out by the Gaussian-sum filter [140].

Several of the estimators described here are implemented in RAVE, a detectorindependent toolkit for reconstruction of interaction vertices [120, 121].

### *13.2.4 Kinematic Fitting*

Kinematic fitting imposes physical constraints on the particles participating in an interaction and thereby improves the measured track momenta and positions. At the same time hypotheses about the interaction and the participating particles can be tested.

#### **<IP_ADDRESS> Lagrange Multiplier Method**

The most commonly used method of imposing constraints on the measured tracks is by way of Lagrange multipliers [141]. Let *q*˜ = *(q*˜ 1*,..., q*˜ *n)*T be the unconstrained estimated parameters of a set of *n* tracks, along with their joint information matrix *G* = diag*(G*1*,..., Gn)* = *V* −1. The *r* functions describing the constraints can be written as *g(q)* = **0**. Taylor expansion around a suitable point *q*˘ yields the linearized equation

$$
\check{\mathbf{g}} + D(\boldsymbol{q} - \check{\boldsymbol{q}}) = \mathbf{0},
$$

where *D* is the Jacobian of *g* with respect to *q*, evaluated at *q*˘, and *g*˘ = *g(q*˘*)*. The obvious expansion point is *q*˘ = ˜*q*. The constrained track parameters *q*¯*i* are obtained by minimizing the objective function

$$M(\boldsymbol{q}, \boldsymbol{\lambda}) = (\boldsymbol{q} - \boldsymbol{\tilde{q}})^{\mathrm{T}} \mathbf{G}(\boldsymbol{q} - \boldsymbol{\tilde{q}}) + 2\boldsymbol{\lambda}^{\mathrm{T}} \left[ \boldsymbol{\check{\mathbf{g}}} + \mathbf{D}(\boldsymbol{q} - \boldsymbol{\tilde{q}}) \right]$$

with respect to *q* and *λ*. *λ* is a vector of *r* unknowns, the Lagrange multipliers. The solution is

$$
\bar{\boldsymbol{q}} = \tilde{\boldsymbol{q}} - \boldsymbol{V} \boldsymbol{D}^{\mathrm{T}} \bar{\boldsymbol{\lambda}}, \quad \text{with} \quad \bar{\boldsymbol{\lambda}} = \boldsymbol{G}\_{\boldsymbol{D}} \left[ \check{\boldsymbol{\mathfrak{g}}} + \boldsymbol{D} (\tilde{\boldsymbol{\mathfrak{q}}} - \check{\boldsymbol{\mathfrak{q}}}) \right] \quad \text{and} \quad \boldsymbol{G}\_{\boldsymbol{D}} = (\boldsymbol{D} \boldsymbol{V} \boldsymbol{D}^{\mathrm{T}})^{-1} \dots
$$

The covariance matrix *V*¯ and the *χ*2 statistic are given by

$$
\bar{\mathbf{V}} = \mathbf{V} - \mathbf{V}\mathbf{D}^{\mathrm{T}}\mathbf{G}\_{\mathrm{D}}\mathbf{D}\mathbf{V}, \quad \chi^{2} = \bar{\lambda}^{\mathrm{T}}\mathbf{G}\_{\mathrm{D}}^{-1}\bar{\lambda} = \bar{\lambda}^{\mathrm{T}}\left[\check{\mathbf{g}} + \mathbf{D}(\tilde{\mathbf{q}} - \check{\mathbf{q}})\right].
$$

If required, the constraint function *g* can be re-expanded at the new point *q*˘ = ¯*q*, and the constrained track parameters can be recomputed.

The Jacobian *D* depends both on the parametrization of the tracks and on the type of constraint to be imposed. For kinematic constraints it is often convenient to choose a parametrization that uses physically meaningful quantities. In [142] it is proposed to use the four-momentum and a point in space, i.e., *q* = *(px, py , pz, E, x, y, z)*. With this parametrization the following constraints can be formulated in a straightforward manner (for further examples see [142]).

#### 13 Pattern Recognition and Reconstruction 595

1. Invariant mass constraint. The equation that constrains a track to have an invariant mass *mc* is

$$E^2 - p\_x^2 - p\_y^2 - p\_z^2 - m\_c^2 = 0.$$

Expanding at *q*˘ = *(p*˘*x, p*˘*y , p*˘*z, E,*˘ *x,*˘ *y,*˘ *z)*˘ yields

$$\mathbf{D} = \left(-2\breve{p}\_{\mathbf{x}} - 2\breve{p}\_{\mathbf{y}} - 2\breve{p}\_{\mathbf{z}}\ 2\breve{E}\ 0\ 0\ 0\right), \quad \breve{\mathbf{g}} = \breve{E}^2 - \breve{p}\_{\mathbf{x}}^2 - \breve{p}\_{\mathbf{y}}^2 - \breve{p}\_{\mathbf{z}}^2 - m\_c^2.$$

2. Total energy constraint. The equation that constrains a track to have a total energy *Ec* is

$$E - E\_c = 0.$$

It follows that

$$D = \begin{pmatrix} 0 \ 0 \ 0 \ 1 \ 0 \ 0 \ 0 \end{pmatrix}, \quad \check{\mathbf{g}} = \check{E} - E\_c.$$

3. Total momentum constraint. The equation that constrains a track to have a total momentum *pc* is

$$
\sqrt{p\_\chi^2 + p\_\chi^2 + p\_\varepsilon^2} - p\_c = 0.
$$

Expanding at *q*˘ = *(p*˘*x, p*˘*y , p*˘*z, E,*˘ *x,*˘ *y,*˘ *z)*˘ yields

$$\mathcal{D} = \left( \frac{\check{p}\_{\times}}{\check{p}} \, \frac{\check{p}\_{\circ}}{\check{p}} \, \frac{\check{p}\_{\circ}}{\check{p}} \, 0 \, 0 \, 0 \, 0 \right), \quad \check{\mathbf{g}} = \sqrt{\check{p}\_{\times}^{2} + \check{p}\_{\circ}^{2} + \check{p}\_{\circ}^{2}} - p\_c.$$

#### **<IP_ADDRESS> Vertex Constraint**

If a vertex constraint is added to the kinematic constraints, the constraint functions depend on the unknown vertex position *v* and are extended to *g(q, v)* = **0**. Taylor expansion around a suitable point *(q*˘*, v*˘*)* yields the linearized equation

$$
\check{\mathbf{g}} + \mathbf{D}(\boldsymbol{q} - \check{\boldsymbol{q}}) + E(\boldsymbol{v} - \check{\boldsymbol{v}}) = \mathbf{0},
$$

where *E* is the Jacobian of *g* with respect to *v*, evaluated at *v*˘, and *g*˘ = *g(q*˘*, v*˘*)*. It is assumed that there is prior information about the vertex position, represented by the position *v*˜ and the covariance matrix *C*. The position *v*˜ can be used as the expansion point *v*˘.

The constrained track parameters *q*¯*i* and the estimated vertex position *v*¯ are obtained by minimizing the objective function

$$M(q, \boldsymbol{\upsilon}, \lambda) = (q - \tilde{q})^T G (q - \tilde{q}) + (\boldsymbol{\upsilon} - \tilde{\boldsymbol{\upsilon}})^T \boldsymbol{C}^{-1} (\boldsymbol{\upsilon} - \tilde{\boldsymbol{\upsilon}}) + 2\lambda^T \left[ \check{\mathbf{g}} + D(q - \tilde{q}) + E(\boldsymbol{\upsilon} - \tilde{\boldsymbol{\upsilon}}) \right]$$

with respect to *q*, *v*, and *λ*. The solution is

$$
\bar{\lambda} = W \left[ \check{\mathfrak{g}} + D(\tilde{\mathfrak{q}} - \check{\mathfrak{q}}) + E(\tilde{\mathfrak{v}} - \check{\mathfrak{v}}) \right], \quad \bar{\mathfrak{v}} = \tilde{\mathfrak{v}} - C E^{\mathrm{T}} \bar{\lambda}, \quad \bar{\mathfrak{q}} = \tilde{\mathfrak{q}} - V D^{\mathrm{T}} \bar{\lambda},
$$

with *W* = *DV D*T + *ECE*T−1 . The covariance matrices are

$$\text{Var}(\bar{\mathfrak{v}}) = \mathcal{C} - \mathcal{C}E^{\mathrm{T}}WEC, \quad \text{Var}(\bar{\mathfrak{q}}) = V - VD^{\mathrm{T}}WDV, \quad \text{Cov}(\bar{\mathfrak{q}}, \bar{\mathfrak{v}}) = -VD^{\mathrm{T}}WEC.$$

The *χ*2 statistic is

$$\chi^2 = \bar{\lambda}^T W^{-1} \bar{\lambda} = \bar{\lambda}^T \left[ \check{\mathbf{g}} + \mathcal{D}(\tilde{\mathbf{q}} - \check{\mathbf{q}}) + E(\tilde{\mathbf{v}} - \check{\mathbf{v}}) \right],$$

with *r* degrees of freedom, where *r* is the number of constraint functions. If the vertex constraint is the only constraint imposed on the tracks, the *χ*2 has 2*n* degrees of freedom. If there is no prior information about the vertex, the prior vertex position is assigned an infinitely large covariance matrix, and *W* is replaced by its limiting value:

$$W = \lim\_{C \to \infty} \left( DVD^\mathrm{T} + ECE^\mathrm{T} \right)^{-1} = \mathbf{G}\_D - \mathbf{G}\_D \mathbf{E} \left( E^\mathrm{T} \mathbf{G}\_D \mathbf{E} \right)^{-1} \mathbf{E}^\mathrm{T} \mathbf{G}\_D.$$

The number of degrees of freedom is reduced to 2*n* − 3.

### **13.3 Track Reconstruction in the LHC Experiments**

### *13.3.1 ALICE*

ALICE [143] is the experiment at the LHC that is devoted to the physics of high energy ion collisions. Its main goal is to investigate the physics of strongly interacting matter and the quark-gluon plasma at extreme values of energy density and temperature in nucleus-nucleus collisions. Among the four experiments at the LHC, ALICE is equipped with the largest number of subdetectors in order to face the reconstruction complexity of ion physics events. In particular, three subdetectors focus on measuring the passage of charged particles using the bending power of the magnetic field. They are assembled in a cylindrical fashion: the Inner Tracking System (ITS) with six planes of high-resolution silicon pixel, drift, and strip detectors, the cylindrical Time-Projection Chamber (TPC) and the Transition Radiation Detector (TRD). The principal functions of the ITS are the identification and reconstruction of secondary vertices, the track reconstruction of low-*p*T particles and the improvement of the impact parameter and momentum resolution. The TPC is the most important tracking sub-detector. Thanks to its time information, it can provide an efficient and robust tracking also in a very high multiplicity environments (in the order of 10,000 charged particles). Finally, the TRD is also used for tracking in the central region and for improving the *p*T resolution at high momentum.

The first step in the track reconstruction in ALICE is the clusterization, which is performed separately for each of the three subdetectors [144]. Tracking then proceeds by determining the preliminary interaction vertex using tracklets defined as lines built with pairs of clusters in the first two layers of the ITS. The preliminary interaction vertex is thus found as a space point to which a maximum number of tracklets converge. In the next step, track finding and fitting is performed in three stages using a inward-outward-inward strategy:

**Fig. 13.4** ITS–TPC matching efficiency vs. *p*T for data and Monte Carlo for pp (left) for Pb-Pb (right) collisions in the ALICE experiment [144]

similar way, the tracks reaching the time-of-flight (TOF) detector are matched to TOF clusters.

• At the final stage of the track reconstruction, all tracks in both ITS and TPC subdetectors are propagated inwards and refitted one last time to determine the final estimate of the track position, direction, inverse curvature, and its associated covariance matrix.

The final interaction vertex is then re-determined using the all tracks reconstructed in TPC and ITS. The precise vertex fit is performed using track weighting to suppress the contribution of any remaining outliers. For data-taking conditions where a high pileup rate is expected, a more robust version of vertex finding inspired by the algorithm described in [132] is used. It is based on iterative vertex finding and fitting using Tukey bisquare weights to suppress outliers. The algorithm stops when no more vertices are identified in the scan along the beam direction. Once the tracks and the interaction vertex have been found, a search for photon conversions and decays of strange hadrons such as *K*0 *S* and0 concludes the central-barrel tracking procedure.

### *13.3.2 ATLAS*

ATLAS [145] is the largest of the four LHC experiments, measuring 25 m in diameter and 44 m in length. Its magnet system is composed of a Central Solenoid Magnet with a 2 T field, a Barrel Toroid and an Endcap Toroids with 4 T each. The Inner Detector (ID) is very compact and highly sensitive in order to measure accurately the decay products of each collision. It consists of three different systems of sensors immersed in the solenoid magnetic field: the Pixel Detector, the Semiconductor Tracker (SCT), and the Transition Radiation Tracker (TRT). The Pixel Detector is situated closest to the interaction point and has the highest granularity with about 80 million readout channels. The intrinsic spatial resolution of the Pixel Detector sensors is 10μm in *r*–*φ* and 115μm in *z*. The SCT is a silicon microstrip detector surrounding the Pixel Detector. It provides eight measurements per track with an overall resolution of 16μm in *r*–*φ* and 580μm in *z*. In the outermost region, the TRT is placed. It is a light-weight detector composed of proportional gas counters (70% Xe, 27% CO2 and 3% O2 straws) embedded in a radiator material and its operational drift radius accuracy is about 130μm. The TRT contributes both to the track pattern recognition stage, featuring typically around 30 hits per track, and to particle identification.

The basic concepts of the ATLAS track reconstruction are described in [146, 147]. The tracking in the ID consists of two principal sequences: an initial inside-out tracking, and a subsequent outside-in tracking. Inside-out tracking starts with space point formation in the silicon part of the ID. Using the space points, track seeds are generated with or without a constraint on the longitudinal vertex position. The seeds are then followed through the SCT by a combinatorial Kalman filter/smoother. After ambiguity solving, the remaining track candidates are extended into the TRT. Outside-in tracking first finds track segments in the TRT, using a Hough transform of the straw centers. A Kalman filter/smoother using also the drift times builds the final track segments. These track segments are then extrapolated back into the SCT and the Pixel Detector. Muons are reconstructed in the ID like any other charged particles; for the standalone reconstruction of muons in the muon system and the combined reconstruction, see [148].

Based on the experience gained in Run 1, several improvements to track reconstruction were made for Run 2 [149]. For example, the tracking was adapted to the new insertable B-layer (IBL) [150], and track reconstruction in dense environments (TIDE) was optimized [151]. This included an artificial neural network based approach to identify pixel clusters created by multiple charged particles. The effect of these two developments is shown in Fig. 13.5. In Fig. 13.5a the transverse impact parameter as a function of track momentum resolution is shown for data taken in 2015 at 13 TeV with the inclusion of the IBL information and for data in 2012 at 8 TeV without the IBL. The data in 2015 was collected with a minimum bias trigger. The data in 2012 is derived from a mixture of jet, tau and missing *E*T triggers [150]. Figure 13.5b shows the improvement of the track reconstruction efficiency in the jet core due to the TIDE optimization [151].

ATLAS track reconstruction efficiency as a function of pseudorapidity and transverse momentum with simulated data at a center-of-mass energy of 13 TeV is shown in Fig. 13.6 [152].

**Fig. 13.5** (**a**) Upper panel: unfolded transverse impact parameter resolution measured from data in 2015 at 13 TeV with the Inner Detector including the IBL, as a function of track *p*T for values of 0*.*0 *<η<* 0*.*2, compared to that measured from data in 2012 at 8 TeV [150]; lower panel: ratio of the resolution in 2015 over the resolution in 2012. (**b**) Improvement of the track reconstruction efficiency due to the TIDE optimization, as a function of the angular distance of the particle from the jet axis. The track selection is explained in [151]

**Fig. 13.6** The track reconstruction efficiency (**a**) as a function of pseudorapidity and (**b**) as a function of transverse momentum, as predicted by Pythia 8 A2 simulation. The statistical uncertainties are shown as black lines, the total uncertainties as green shaded areas [152]

### *13.3.3 CMS*

CMS [153], together with ATLAS, is one of the two general-purpose experiments at the LHC. Its main distinguishing feature is a 3.8 T superconducting solenoid. With a length of 13 m and a diameter of 6 m, it provides a high bending power to precisely measure the momentum of charged particles. The solenoid magnetic field lines run parallel to the beam direction in the central region, where the tracking system is placed. The tracking system is designed to provide a precise and efficient measurement of particle trajectories using position-sensitive detectors. The CMS tracker is a silicon-based system [154]. It splits into two parts, the Pixel Tracker and the Strip Tracker and covers a pseudorapidity range up to |*η*| = 2*.*5. The Pixel Tracker is the innermost CMS detector sub-system and is composed of 66 million silicon pixels with dimensions 100 × 250 × 250 μm, covering a total area of about 1 m2. In the barrel layers the magnetic field induces a Lorentz angle which increases charge sharing between neighbouring pixels. Charge sharing in conjunction with analog readout allows to achieve 10 μm position resolution for the *(r, φ)* coordinate and 15 μm in the *z* direction. The pixel detectors in the forward direction are tilted at an angle of 20◦ to induce charge sharing which allows to achieve 15 μm and 20 μm resolution respectively. This resolution is not only necessary for a precise track reconstruction, but also for the determination of both the vertices produced in the primary interaction and the decay vertices of short-lived particles.

The Strip Tracker constitutes the outer part of the tracking system. Its basic building blocks are silicon strip modules. Each module is equipped with one or two silicon sensors and a so-called Front-End hybrid containing readout electronics. In total, the CMS silicon strip tracker has 9.3 million strips and covers 198 m2 of active silicon area. The resolution in *(r, φ)* is 30μm in all layers. The inner layers of the strip tracker are equipped with double-sided sensors, one side of which is rotated by a stereo angle of 100 mrad, achieving a resolution along the *z* coordinate of about 230 μm and allowing the reconstruction of the hit position in 3-D. In the outer layers the sensors are single-sided, and the *z* resolution can be approximated by the strip length over √12, or about 15 mm. In order to maintain excellent tracking performance until the Long Shutdown 3 of LHC, the Pixel Tracker was replaced in the year-end technical stop of 2016*/*2017 with a new Pixel Tracker composed of four barrel layers and six forward disks providing four-hit pixel coverage up to |*η*| = 2*.*5. After the Long Shutdown 3, the High Luminosity phase of the LHC (HL-LHC) is scheduled where the accelerator will provide an unprecedented instantaneous luminosity of 5−7*.*5×1034 cm−2 s−1. In [155] the new CMS silicon tracker and its tracking and vertexing performance for different event types, pileup scenarios and detector geometries are presented.

The CMS track reconstruction algorithm is based on an iterative approach [156]. The main idea is to search for easier-to-find tracks first, to mask the hits associated to the found tracks, and to proceed to the next iteration. In this way the combinatorial problem is reduced, and the search for more difficult classes of tracks is simplified. Moreover, this approach introduces the possibility of developing special iterations that can improve track reconstruction in high-density environments such as jets, or to use the information from other subsystems such as muon chambers and calorimeters. In each iteration, the Combinatorial Track Finder is run. It can be divided into four different steps:

As an example, the twelve tracking iterations foreseen for 2017 data taking is listed in Table 13.3 [157]. The main difference between iterations is the configuration of the seed generation and the target tracks.

Figure 13.7 shows the tracking efficiency, using a standard sample of tt events simulated with √*s* = 13 TeV with different superimposed pileup conditions. The contribution of different iterations for 2017 track reconstruction is also shown as a function of the *p*T of the simulated particle. It can be seen how iterations targeting low-*p*T tracks are more efficient in the region between 100 and 500 MeV.

**Table 13.3** List of different tracking iterations used after the Pixel Tracker upgrade with the corresponding seeding configuration used and target tracks [157]

**Fig. 13.7** Track reconstruction efficiency as a function of simulated track pseudorapidity for 2017 tracker at different pileup conditions (left) and cumulative contributions to the overall tracking performance from the twelve iterations in 2017 track reconstruction shown as a function of the simulated track *p*T (right) [157]. The 2017 tracking reconstruction includes the Cellular Automaton-based Hit Chain-Maker (CA) seeding [158]

Figure 13.8 shows the muon tracking efficiency and the corresponding ratios between real and simulated data for 2016 collisions data coming from the Z resonance using the tag and probe method. The measured track efficiency as a function of |*η*| is found to be between 99.5% and 100% for the collection including all tracks. It degrades, however, with increasing number of primary vertices.

**Fig. 13.8** Data (black dots) and simulation (rectangles) tracking efficiency and respective ratio for muons coming from the Z decay as a function of the absolute pseudorapidity of the probe muon (left) and the number of primary vertices (right). The data are based on an integrated luminosity of 36 fb−1 [155]

### *13.3.4 LHCb*

As its name indicates, LHCb [159] focuses on physics involving bottom quarks and investigates CP violation phenomena. These studies require the measurement of the rare decays of Bd, Bs, and D mesons which are produced with a large crosssection at the LHC. Given the fact that b hadrons are predominantly produced in the forward or backward cone, the LHCb experiment is a single-arm spectrometer in contrast to the other three experiments. In order to exploit this large number of b hadrons, it requires a robust and flexible trigger and a data acquisition that allows high bandwidth data taking and provides powerful online data processing. Furthermore, superior vertex and momentum resolution are crucial to study the rapidly oscillating Bs − Bs meson system. LHCb is thus equipped with the highly sophisticated silicon microstrip detector close to the interaction point, the Vertex Locator (VELO). It can be moved to a distance of only 7 mm from the proton beams and measures the position of the primary vertices and the impact parameters of the track with extremely high precision. A further silicon microstrip detector, the Tracker Turicensis (TT) is placed before the dipole magnet. Its task is to improve the momentum resolution of reconstructed tracks and reject pairs of tracks that in reality belong to the same particle. The magnet is placed behind the TT. It bends the flight path of the particles in the *x* − *z* plane and therefore allows the determination of their momenta. The tracking system is completed by the T stations (T1-T2-T3), which, together with the information from the VELO, determine the momentum and flight direction of the particles. The T stations are composed of silicon microstrip sensors close to the beam pipe and by straw tubes in the outer regions.

Track reconstruction uses hits in the VELO, TT and T stations. Depending on which detectors are crossed, different track types are defined [160, 161]:

Reconstruction of long tracks starts in the VELO. There are two complementary algorithms to add information from the downstream tracking stations to these VELO tracks. The first one combines the VELO tracks with information from the T stations. The second one combines the VELO tracks with track segments found after the magnet in the T stations, using a standalone track finding algorithm. The candidate tracks found by each algorithm are then combined, removing duplicates, to form the final set of long tracks used for analysis. Finally, hits in the TT consistent with the extrapolated trajectories of each track are added to improve their momentum determination.

Downstream tracks are found starting with T tracks, extrapolating them through the magnetic field and searching for corresponding hits in the TT. Upstream tracks are found by extrapolating VELO tracks to the TT where matching hits are then added in a procedure similar to that used by the downstream tracking. At least three TT hits are required to be present by these algorithms.

The found tracks are fitted using a Kalman filter, taking into account multiple scattering and energy loss due to ionisation. The *χ*2-statistic of the fit is used to determine the quality of the reconstructed track. If two or more tracks have many hits in common, only the one with most hits is kept.

The track reconstruction efficiency for the 2012 and the 2015 data as a function of the momentum can be seen in Fig. 13.9 [162]. The results of the two periods are compatible.

### **13.4 Conclusion**

An overview of current methods in track and vertex reconstruction and alignment has been presented. Many of them have been developed in response to the requirements of the current experimental program at the Large Hadron Collider. The most difficult challenges are:

Every experiment has to meet these challenges on its own terms. The outlines of the solutions found by the four major LHC experiments are described in Sect. 13.3 and, in more detail, in the references given there. In addition, the repertory of the methods discussed in this contribution can certainly not lay claim to completeness. We have tried to select widely applicable methods, thereby neglecting by necessity many experiment specific adaptations, improvements and innovations, for which we again refer to the references.

### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 14 Distributed Computing**

**Manuel Delfino**

Distributed computing is an established discipline in computer science and engineering. It has evolved over the past 40 years to be one of the most important methodologies for implementing the data processing services needed by almost every activity in society. Distributed computing is still evolving rather rapidly, with major innovations introduced every few years. The aim of this chapter is to introduce the reader to the basic concepts of distributed computing in the context of particle physics, allowing a better understanding of what is behind the scenes when using distributed systems, and providing starting points in order to seek further information on the subject.

### **14.1 Usage by Particle Physics**

Distributed computing, i.e. the coherent use of many computers to accomplish a given task, is extensively used in particle physics. It is particularly suited to simulating and analyzing data from particle collisions, where processing the data of one collision is largely independent from processing data of all other collisions and hence very little communication between the different computers is necessary. This type of distributed computing is called "High Throughput Computing" (HTC),1 because there is a very large number of quasi-independent tasks to accomplish,

M. Delfino (-)

1Livny, Miron, et al. "Mechanisms for high throughput computing." SPEEDUP journal 11.1 (1997): 36–40.

Departament de Física, Facultat de Cièncias, Universitat Autònoma de Barcelona, Barcelona, Spain e-mail<EMAIL_ADDRESS>
and the performance perceived by the user is the rate of completion of tasks, or "throughput". This type of computing is complimentary to "High Performance Computing" (HPC), i.e. the execution of a single task at maximum speed on a classical supercomputer. The particular case of event or collision processing requiring little inter-processor communication is an example of "loosely coupled parallel computing" or even "embarrassingly parallel computing".2,3 Certain kinds of accelerator simulations, particularly those that use ray-tracing techniques, can also be executed as distributed computing tasks. Multiple computer systems can also be interconnected in "tightly-bound" configurations via low-latency networks,4 essentially yielding a supercomputer. Tightly-bound systems are used in particle physics, for example, for solving equations from theories using numerical methods (e.g. lattice-gauge theories); they are mentioned here for completeness and will not be covered further.

From the decade of the 2010s, the boundary between HTC and HPC has started to blur, for a number of disparate reasons:

*Multi-core Applications* Detectors have become more complex and have more channels. This impacts the memory footprint of reconstruction, simulation and analysis programs, which has grown substantially. In parallel, processors have become multi-core5 with shared RAM.6 The growth in number of cores per processor has been faster than the drop in price of RAM chips, resulting in an effective memory shortage. This means that running an independent copy of the operating system and the application in each core is not economical. The way out is to implement the safe execution of parallel threads7 of a single copy of the program on each multi-core processor. This introduces a number of dependencies, for example competition for RAM and for input–output services, that break the loose coupling and make the applications behave more like HPC programs.

2Wilkinson, Barry, and Michael Allen. Parallel programming. Pearson India, 2004.

3Birrittella, Mark S., et al. "Intel® Omni-path architecture: Enabling scalable, high performance fabrics." High-Performance Interconnects (HOTI), 2015 IEEE 23rd Annual Symposium on. IEEE, 2015.

4Pfister, Gregory F. "An introduction to the infiniband architecture." High Performance Mass Storage and Parallel I/O 42 (2001): 617–632.

5A core corresponds to the Central Processing Unit (CPU) of a classical computer, which is able to execute a single stream of instructions. Technology allows a growing number of ever smaller transistors to be placed on a single chip. However, profiting from these large number of transistors in a single core design would require prohibitively complex processor designs and impossible to achieve clock speeds. The alternative is to populate the chips with many copies of the same processor core, which work independently except for sharing external connections, leading to multi-core processors. Gepner, Pawel, and Michal Filip Kowalik. "Multi-core processors: New way to achieve high system performance." Parallel Computing in Electrical Engineering, 2006. PAR ELEC 2006. International Symposium on. IEEE, 2006.

6Random Access Memory, the external solid state memory used by a processor.

7A thread is a stream of processing instructions coming from a shared program image which has its own private data instances in processor hardware registers, instruction and data stack and caches, and RAM.

*Supercomputers as High Core-Count Clusters* The original design of supercomputers was based on a rather small number of the fastest single core processors available at the time, interconnected with custom very low latency links. In addition, the design usually included a much larger amount of RAM than in standard computers. The design of supercomputers has changed completely over the last 20 years.8 The main driver of this change is the difficulty in building processors with ever shorter clock cycles. The so-called "clock speed wall" has been hit, and practically all processors operate within a narrow range around a few GHz. Hence, the only way to make a supercomputer faster is through parallelism. Modern supercomputers are in fact huge clusters of relatively standard multi-core processors. Low latency links continue to be used and their hardware costs have become much lower. In addition, "latency hiding" techniques sometimes allow the use of standard networks, such as Ethernet. Finally, the amount of RAM per core does not differ much from standard computers. Hence, executing an HTC workload on a supercomputer targeted for HPC is no longer wasteful. At worst, only the low latency interconnect will be underutilized.

*Increase of Workloads with Low Input–Output* Particle physics has traditionally been a heavy user of HTC systems because they were the least expensive architecture for executing detector track reconstruction and analysis, applications with a relatively low CPU to input–output ratio. The needs for CPU for reconstruction and analysis have grown, however, as detectors have become much more complex. In parallel, the precision needed in the simulations has vastly increased. These high CPU, low input–output applications currently represent the largest computing demand of a modern particle physics detector. Hence, the global workload profile has moved closer to HPC in the last decades.

### **14.2 Functional Decomposition of a Distributed Computing Environment**

It is useful to introduce a functional decomposition, or reference framework, to discuss distributed computing systems. The decomposition that has dominated particle physics data processing since the 1990s, is the "SHIFT Model", introduced in the early 1990s by Robertson and collaborators.9 It will be described here with some updates of terminology. The basic assumption in the SHIFT framework is that a Local Area Network (LAN) can be built with enough capacity and flexibility so that the rest of the elements in the distributed environment can

8Xie, Xianghui, et al. "Evolution of supercomputers." Frontiers of Computer Science in China 4.4 (2010): 428–436.

9J.P. Baud, et al. "SHIFT, the Scalable Heterogeneous Integrated Facility for HEP Computing", Proc. Conference on Computing in the High Energy Physics, CHEP91, Tsukuba, Japan, Universal Academic Press.

communicate freely amongst themselves. Furthermore, the environment is loosely coupled and is therefore not sensitive to relatively long round-trip time for network messaging (large network latency). Connected to the network we have the following elements:

The elements of the functional decomposition presented, illustrated in Fig. 14.1, together with the aforementioned powerful network and the maintenance of coherent security10 domains, constitute a framework which is sufficient to analyze the most commonly used distributed computing environments.

The computing industry places a lot of emphasis on how storage is attached to the computers: Direct-Attached Storage (DAS) does not use a network, attaching storage hardware via an internal communication bus of a computer; Storage Area Network (SAN) uses a short-distance network (originally Fiberchannel, with an increasing tendency towards Ethernet) to couple storage hardware to computers; and Network Attached Storage (NAS) which uses computers hosting storage devices and connected to a LAN (or WAN configured to behave like a LAN) to present storage services to the other computers on the network. These distinctions are important in environments where few computers are used. In the case of particle physics, DAS and SAN are used to construct disk servers, which in turn are exposed almost exclusively in NAS mode.

10A more precise term is "Authentication and Authorization Infrastructure", abbreviated as AAI or AA.

### **14.3 Data Processing Clusters**

The most common way to deploy a distributed computing environment is a "cluster".11 A cluster consists of a number of CPU, Disk and Control servers connected on a LAN12 and sharing status information through local Information servers while maintaining a *single, coherent, security domain*. A cluster may also implement Tape servers and Remote Data servers.

Historically, the first complete and reliable commercial implementation of a cluster was the VAXcluster13 (or VMScluster), implemented by Digital Equipment Corp. starting in 1983 on VAX computers and workstations using the VMS operating system. Today, essentially all clusters used by particle physics are implemented using the Linux operating system14 as a basis, and adding a number of additional packages for handling authentication and storage.

11Unfortunately, after many decades of using the term "cluster" as in the present text, some commercial firms have taken to using the terms "Grid" or "Cloud" to mean "cluster". This is considered by the author as incorrect and confusing.

12There are also implementations using several LANs transparently interconnected via WANs, but maintaining a single security domain. This configuration is often used in business environments to gain fault-tolerance and disaster recovery. Problems related to performance limitations of the WAN, as well as inter-institutional security management complications, limit the use of this configuration in particle physics data processing. Interconnection of clusters via a Grid or a federated Cloud is used instead.

13Kronenberg, Nancy P., Henry M. Levy, and William D. Strecker. "VAXcluster: a closely-coupled distributed system." ACM Transactions on Computer Systems (TOCS) 4.2 (1986): 130–146.

14Linux (or more precisely GNU/Linux) is mentioned specifically as it has become the dominant Unix-like operating system, and has become better known than the original Unix operating system. Most of the remarks in this chapter, however, are equally applicable to other Unix-like systems, notably macOS and FreeBSD.

### *14.3.1 Authentication and User Identification*

The machinery behind authentication and user identification is often poorly understood by users, which leads to usage and security problems. It is important to gain an understanding of the basic concepts, especially as more complex distributed environments, such as Grids and Clouds (see sections below), have come into widespread use.

A number of Information servers are deployed to integrate individual computers to form a cluster.15 The most crucial ones are the ones that generate a single, coherent, security domain.16 Security under Linux and most other operating systems is based on the concept of an "account" which, ideally, is used exclusively by a single trusted individual. An account is represented externally by an alphanumeric string (the "username") and internally by a numeric code (the *userid* or *uid*). Although modern operating systems allow very long usernames, local restrictions may apply in order to maintain compatibility with older system software and utilities (the most common restriction is a maximum length of eight characters). Each account is associated to specific directives (called "rights" or "privileges") to allow or disallow access to operating system services, and can also be used to control access to files and other resources directly or via Access Control Lists ("ACL"). Since people often work in teams, it is useful to associate an account with a "group" named by an alphanumeric string and represented internally by a numeric code (the *groupid* or *gid*). Certain rights and, very importantly, file access control can thus be quickly managed for all accounts belonging to a group.

The process of authentication involves an exchange of credentials which establishes the identity of the user wanting to access a system which results in the creation of a process running under the corresponding *uid*. Authentication has been traditionally accomplished by a "password", an alphanumeric string which, ideally, is known only to the individual owner of an account. Password related issues are major contributors to computer security problems. Practically all computers are nowadays connected directly or indirectly to the Internet. An inherent weakness of a cluster (which is more than compensated by the gained functionality) is that having a single security domain means that gaining access to any component of the cluster grants access to all of it. The components of clusters directly connected to the Internet are under constant attacks, the most common being automated attempts to guess passwords (usernames are easily obtained from public information sources, such as Web pages listing email addresses). Attacks from within the cluster must also be considered, especially in clusters with many accounts. A

15The term Information server is used in a broad manner. For example, Domain Name Service (DNS) servers which translate alphabetic Internet addresses to numerical, and Network Time servers which ensure all cluster elements have their clocks synchronized, are considered Information servers.

16The actual implementation is a highly technical matter, using NIS, ldap, Microsoft Active Directory or other secure database sharing schemes.

password which is easily guessed is denominated "weak". Simple measures exist to avoid weak passwords:17 The password should be as long as possible and contain a mix of numbers and upper and lower case letters. Many installations require users to change their passwords periodically in order to improve computer security. A very dangerous practice which should be completely avoided, but which is unfortunately common practice in the particle physics community, is the "sharing" of personal accounts (and hence, their passwords) or the use of "service" accounts with passwords known by dozens of persons, written in documents and blackboards or even posted on Web sites. Modern operating system features make these practices completely unnecessary, as they can be configured so that each user can first authenticate with their personal account and then gain access to a common environment to perform the tasks required.

Portable devices ("laptops" and "smartphones") can be deployed so that they are elements of a cluster and therefore part of the security domain. Since the portable computer must remain useable when temporarily disconnected from the network, the operating system will keep a local copy (or "cache"18) of the security data, including the passwords. This means that a stolen portable computer is a computer security threat. Care should be taken in configuring the "suspended" or "hibernated" modes in portable computers to ask for a password when they are turned back on.

Another major issue with authentication based on usernames and passwords is that each cluster (and each online service such as electronic mail, social networks, and document and file repositories), being a separate security environment, requires a separate username and password and has its own policies for requiring password changes. Users are becoming overwhelmed in keeping track of their usernames and passwords and their reaction is often to use weak passwords or keep a list of passwords in a file which can be stolen. In order to reduce username/password proliferation, many organizations are linking their clusters and online services to shared authentication servers, thus providing a "single-sign-on" environment.

In the 2000s, organizations, most notably the CERN Large Hadron Collider experiments, started to use an alternative authentication and authorization frame-

17Users should consult their local security rules for specific recommendations on password choices.

18The term "cache" is used in computer architecture to describe a local copy of a limited amount of data which is normally stored elsewhere. The introduction of this local copy serves to increase the performance and efficiency of use of the element connected to it, by smoothing out peaks and valleys in storing and retrieving information, or avoiding retrieves altogether if the same information is needed repetitively. The cache is transparent to the element accessing the information; if the data needed is in cache, it will be delivered from there (a cache "hit") whereas if the data is not in cache, the normal access path to the normal storage place will be used (a cache "miss"). The key to successful implementation of a cache is the right choice of algorithm to choose which data to copy and keep in the cache. Examples are high-speed memory caches in CPU chips which hold local copies of instructions and data, and memory caches in disk controllers which hold local copies of disk blocks.

work based on digital certificates.19 This was one of the pillars for constructing the Worldwide LHC Computing Grid (WLCG), later generalized to the European Grid Infrastructure (EGI), the U.S. Open Science Grid (OSG) and others. The certificates used are based on the X.509 standard, supported by all major Web servers and browsers. These schemes were very successful and are still being used, allowing tens of thousands of computers to provide coherent services to thousands of users. They are, however, extremely difficult to manage and expensive to operate. Hence, the tendency is to abandon these schemes in user-facing services, keeping them mostly for machine-to-machine services.

Industry and academia are putting a lot of effort on the emerging area of "distributed" authentication and authorization (AA) schemes.20 The ultimate aim is to simplify the AA process and make it easier to use, enabling "single-sign-on" to a wide variety of resources, a functionality that users have come to expect from their experience with social networks. Unfortunately, progress has been rather slow and solutions proposed by industry and academia are often incompatible. Part of the problem is that there are two competing standards: SAML21 and OpenID.22

The main problem, however, is that the usage by large collaborations, such as those in particle physics, ideally requires the simultaneous use of many different AA sources (something known as Identity Federation). The largest Identity Federation currently deployed is the *eduGAIN* federation,23 where the AA sources come from the universities and research institutes who employ the users. *eduGAIN* is modeled after the successful *eduroam* system24 used to grant worldwide access to academic WiFi networks worldwide. Solving the general access problem is much more difficult, however, as it requires the maintenance of many more attributes for each user. For example, each university and research institute would have to include and maintain in their databases which experiment or project each employee is participating in, something which is not practical. In order to solve this issue, hybrid schemes are being worked on, where the authentication would come from *eduGAIN* but the authorization information would come from an attribute server managed and operated by a specific experiment or project. This still leaves the problem of reliably operating a service build from thousands of independently managed AA servers with varying degrees of service quality.

19Thompson, Mary R., Abdelilah Essiari, and Srilekha Mudumbai. "Certificate-based authorization policy in a PKI environment." ACM Transactions on Information and System Security (TISSEC) 6.4 (2003): 566–588.

20These schemes centralize AA information and then make it available to a distributed set of heterogeneous resources, hence the term "distributed".

21Rosenberg, Jonathan B., and David L. Remy. Securing web services with WS-security: Demystifying WS-security, WS-policy, SAML, XML signature, and XML encryption. Sams, 2004. 22https://openid.net/developers/specs/

23López, D. "eduGAIN: Federation interoperation by design." TERENA Networking Conference. 2006.

24López, Gabriel, et al. "A proposal for extending the eduroam infrastructure with authorization mechanisms." Computer Standards & Interfaces 30.6 (2008): 418–423.

Some interesting alternatives are starting to emerge, which avoid Identity Federation. One approach is to identify a single reliable provider for AA, for example the OpenID service from ORCID25 combined with an appropriate attribute server. Another approach is for a project or experiment to deploy a single AA scheme available to services *for that project* around the world. This can be accomplished in a simple and economical manner by re-using the project's personnel database and exposing it through the Internet in a secure manner using an *ldap* server.26 These type of implementations of single-sign-on are ideal for service providers that serve a single project. On the other hand, they do shift the diversity problem to the service providers, which must locally adjust crucial attributes such as *uid* and *gid* in order to avoid duplications in multi-project or multi-experiment environments.

### *14.3.2 Processing and Storage*

Two configurations of CPU servers are usually deployed:

Several configurations of Disk servers are usually deployed:

• "Network File servers", based on industry-standard or widely used protocols, such as NFS,29 smb230 or http with WebDAV,31 are used to provide "home" and "project" directories holding text and binary files of relatively small size. These directories appear in the CPU servers as "mount points" or "network shares" that largely behave as a local, conventional disk resource. The POSIX input–output

25Haak, Laurel L., et al. "ORCID: a system to uniquely identify researchers." Learned Publishing 25.4 (2012): 259–264.

26Johner, Heinz, et al. Understanding LDAP. Vol. 6. IBM, 1998.

27See for example Yoo, Andy B., Morris A. Jette, and Mark Grondona. "Slurm: Simple linux utility for resource management." Workshop on Job Scheduling Strategies for Parallel Processing. Springer, Berlin, Heidelberg, 2003.

28Fajardo, E. M., et al. "How much higher can HTCondor fly?" Journal of Physics: Conference Series. Vol. 664. No. 6. IOP Publishing, 2015.

29Shepler, Spencer, et al. Network file system (NFS) version 4 protocol. No. RFC 3530. 2003.

30French, Steven M., and Samba Team. "A New Network File System is Born: Comparison of SMB2, CIFS and NFS." Linux Symposium. sn, 2007.

31Goland, Yaron, et al. HTTP Extensions for Distributed Authoring – WEBDAV. No. RFC 2518. 1999.

standard32 is the most often used definition of this behavior. The implementations are far from trivial, but fortunately they are widely used and hence available as packages fully supported by operating systems, or even as part of the operating system. These protocols support operations such as simultaneous access by multiple users in read/write/append mode as well as record-level locking. This forces the deployment in CPU servers of some Control and Information server functionality, with information being constantly exchanged between all nodes accessing a mount point in order to keep file system coherency and distribute lock information. In order to fulfill performance expectations, these Control and Information functionalities must often be implemented as drivers, often executing inside the operating system kernel. The price of this is the risk of all CPU servers becoming inoperative (or "hung") if a network file server has a fault or the status information becomes corrupted. Some Network File servers, such as the CernVM-FS33 system used for software installation, serve very specialized purposes. They implement only part of the functionalities, for example serving files only in read-only mode, with the benefit of having a smaller impact on the operating system environment.

• "Database servers", are disk servers34 which organize their data using database products such as Oracle, MySQL or PostgreSQL. They are mostly used to store and manage detector configuration and calibration data, though on occasion they can store high-level analysis data such as the "Event Tags" used by some experiments.35 CPU servers read these data through the network by issuing database queries and retrieving the results. Database servers must be configured in accordance to the type of most-often used queries (a process called "tuning") in order to achieve the required performance. Because the overall cost of purchasing and maintaining a Database server is high, situations often arise where the deployed capacity is insufficient to serve all the CPU servers unless specific caches or buffers36 are deployed. An emerging alternative type of database

32Gallmeister, Bill. POSIX. 4 Programmers Guide: Programming for the real world. "O'Reilly Media, Inc.", 1995.

33Blomer, Jakob, et al. "Status and future perspectives of CernVM-FS." Journal of Physics: Conference Series. Vol. 396. No. 5. IOP Publishing, 2012.

34It may seem odd to classify a database server as a Disk server. It is considered correct in this context, as their usage pattern is unusual in that the record update frequency is very low compared to the record reading frequency. Note that the technical database servers used to implement Data and Tape server features are considered Information servers.

35See for example Cranshaw, Jack, et al. "Event selection services in ATLAS." Journal of Physics: Conference Series. Vol. 219. No. 4. IOP Publishing, 2010.

36The term "buffer" is used in computer architecture to mean "temporary storage location". Buffers are used to accumulate information and then transmit it as a block, or to receive information as a block and then distribute it. The difference between a buffer and a cache is that in the case of buffers all the information is temporarily stored, whereas a cache implements an algorithm to store only the most relevant information. In addition, whereas a cache is always transparent, a buffer may be transparent or it may allow explicit buffer manipulation by the application.

servers are those based on "Big Data" tools, such as NoSQL databases37 with underlying *Hadoop* storage management systems.38

• "Data servers", based on products such as dCache,39 GPFS,40 Lustre,41 EOS42 and DPM,43 are used to hold the large-size binary files which contain physics data. They use a set of specific Control and Information Servers to present to the other elements of the cluster, in particular to the CPU servers, an interface approximating a conventional (POSIX-compliant) disk resource, but with some important differences in order to achieve the large capacities and high performance needed for particle physics applications (and data-intensive applications in other fields). Much of the optimization is possible because of the particular usage pattern of particle physics data, which is essentially written once and then read many times.44 The key point is to separate45 the "namespace" information (directory tree or file folder structures, file attributes such as name, creation date, etc.), the "file storage" information (physical location of the blocks that comprise a file, whether it is open for read, write or append, etc.) and the "file data" information (the actual data blocks that comprise the file). This separation makes it possible to serve the data blocks from a large number of Data servers with a combined capacity far exceeding what can be provided by a conventional file system. In addition, the data blocks of a given file may be stored in multiple servers in order to enhance read performance or provide high-availability (for reading).46 The cost of this separation (apart from the complexity of the multiple

37Han, Jing, et al. "Survey on NoSQL database." Pervasive computing and applications (ICPCA), 2011 6th international conference on. IEEE, 2011.

38Shvachko, Konstantin, et al. "The hadoop distributed file system." Mass storage systems and technologies (MSST), 2010 IEEE 26th symposium on. Ieee, 2010.

39Millar, Paul, et al. "Storage for advanced scientific use-cases and beyond." Parallel, Distributed and Network-based Processing (PDP), 2018 26th Euromicro International Conference on. IEEE, 2018.

40Schmuck, Frank B., and Roger L. Haskin. "GPFS: A Shared-Disk File System for Large Computing Clusters." FAST. Vol. 2. No. 19. 2002.

41Schwan, Philip. "Lustre: Building a file system for 1000-node clusters." Proceedings of the 2003 Linux symposium. Vol. 2003. 2003.

42Peters, A. J., E. A. Sindrilaru, and G. Adde. "EOS as the present and future solution for data storage at CERN." Journal of Physics: Conference Series. Vol. 664. No. 4. IOP Publishing, 2015.

43Alvarez, Alejandro, et al. "DPM: future proof storage." Journal of Physics: Conference Series. Vol. 396. No. 3. IOP Publishing, 2012.

44Aside from the write-once/read-many access pattern, particle physics datasets almost never have data appended to them nor have their already existing records updated. For an extreme contrast, consider the access pattern of bank account datasets, which constantly have data appended and updated.

45Conventional file systems, whether local or networked, tightly bind the management of the three types of information described in order to optimize file open/close operations and implement sharing of files in read/write/append mode.

46This is conceptually similar to the "striping" of files across many disks performed by RAID controllers on DAS devices.