Combination of various data analysis techniques for efficient track reconstruction in very high multiplicity events

Abstract. A novel combination of established data analysis techniques for reconstructing charged-particles in high energy collisions is proposed. It uses all information available in a collision event while keeping competing choices open as long as possible. Suitable track candidates are selected by transforming measured hits to a binned, threeor fourdimensional, track parameter space. It is accomplished by the use of templates taking advantage of the translational and rotational symmetries of the detectors. Track candidates and their corresponding hits, the nodes, form a usually highly connected network, a bipartite graph, where we allow for multiple hit to track assignments, edges. In order to get a manageable problem, the graph is cut into very many minigraphs by removing a few of its vulnerable components, edges and nodes. Finally the hits are distributed among the track candidates by exploring a deterministic decision tree. A depth-limited search is performed maximizing the number of hits on tracks, and minimizing the sum of track-fit χ2. Simplified but realistic models of LHC silicon trackers including the relevant physics processes are used to test and study the performance (efficiency, purity, timing) of the proposed method in the case of single or many simultaneous proton-proton collisions (high pileup), and for single heavy-ion collisions at the highest available energies.


Introduction
Traditional methods of track reconstruction can be scaled to work in high multiplicity events, namely in many simultaneous collisions (pileup) of elementary particles [1,2] and in high multiplicity single heavy-ion collisions.Nevertheless the performances are not optimal, efficiency and purity are reduced, especially at low momentum.That is why present data taking conditions and further luminosity and energy upgrades of high energy particle colliders, as well as those of detector systems, call for new ideas.
Image transformation methods and neural networks [3] are often used in gas detectors (time projection chambers [4,5] and transition radiation trackers [6,7]).In the case of silicon trackers the combinatorial track finding methods employed for trajectory building mostly use local information [8,9].They start with a trajectory seed and build a trajectory by extending the seed through the detector layers, picking up compatible hits.In the case of very many compatible hits the number of concurrently built trajectory candidates must be limited.Only some of the best candidates are kept which biases the final result.In this sense decisions are made too early.In addition, trajectories are mostly treated separately, there is no interaction between their assigned hits.
In this study a combination of established data analysis techniques for charged-particle reconstruction is proposed.It uses all information available in an event while keeping competing choices open as long as possible.
2 Silicon detectors at particle colliders At currently operating particle colliders the interaction region is very narrow (of the order of 50 μm) in both x and y (transverse) directions, while in z (longitudinal or beam) direction it is long, with a characteristic size of about 10 cm [10].For single heavy-ion collisions the z 0 position of the primary interaction (vertex) is estimated with good precision using the copiously produced high p T particles, thanks to the small pointing uncertainty of their tracks, reconstructed with traditional methods.In the case of single or multiple pp collisions no such information on the locations of the interaction vertices exists.
The trajectory of a primary particle is primarily determined by its initial position (0, 0, z 0 ) and parameters (q, η, p T , φ 0 ) of its initial momentum at creation.Here q is the charge, η = − ln tan(θ 0 /2) is the pseudorapidity, p T is the transverse momentum, φ 0 and θ 0 are the azimuthal and polar angle of the initial momentum vector in spherical coordinates.In a large volume solenoid the magnetic field near the center of the detector is rather homogeneous and points in the z direction.Hence in small volumes the trajectories of charged particles can be approximated by piecewise helices.For practical purposes a primary particle is parametrized by (k T , sinh η, φ 0 , z 0 ) in the following, where k T = q/R is the signed curvature of the projection of its trajectory on the transverse (bending) plane, R is its radius.If the particle is singly charged the curvature is connected to p T as p T = eB z R, where e is the electric charge of a proton.
Silicon trackers generally consist of several concentric cylindrical layers.Those close to the nominal interaction point are equipped with tiny pixel sensors, while others contain long strip sensors.Some strip layers are double-sided, they are located very close to each other two by two, and have a small relative tilt angle.The main characteristics of the inner barrel silicon detectors of the studied experimental setups are given in Table 1.
The trajectory of the primary particle intersects the concentric cylindrical layers and leaves hits behind in the silicon (Fig. 1).In the simplified case, when the magnetic field is homogeneous and if  the detector material and its physical effects are neglected, the position of those hits could be precisely determined by simple equations.The physical effects of detector material changes this overly simple picture.

Physical effects
When a long lived charged particle propagates through material the most important effects which alter its momentum vector are multiple scattering and energy loss.The distribution of multiple Coulomb scattering is roughly Gaussian [11], the standard deviation of the planar scattering angle is where p, βc, and z are the momentum, velocity, and charge of the particle in electron charge units, and x/X 0 is the thickness of the scattering material in radiation lengths.Momentum and energy is lost during traversal of sensitive detector layers and support structures.To a good approximation the most probable energy loss Δ p , and the full width of the energy loss distribution at half maximum Γ Δ [12] are where ξ = K 2 z 2 Z A ρ x β 2 is the Landau parameter; K = 4πN A r 2 e m e c 2 ; m is the mass of the particle; Z, A, I, and ρ are the mass number, atomic number, excitation energy, and the density of the material, respectively [11].The density correction δ is neglected.

Hit clusters
An incoming charged particle loses energy in the sensitive detector elements by producing electronhole pairs.The neighboring channels collecting a charge above a given threshold are grouped to form a cluster, a reconstructed hit.The size (dimensions) of the cluster depends on the angle of incidence of the particle: bigger angles will result in larger clusters.Due to the large fluctuations in energy loss, the measured w rφ and w z dimensions of the clusters will differ from the expected ones.In order to model these effects in the simulation, the cluster dimensions are varied by one unit in both directions for pixels, and up and down by two units in rφ direction for strips.In addition, if the size of a pixel cluster is at least two units in both directions, the layout of pixels with charge deposit and its location relative to the nominal interaction point usually indicate the sign of the electric charge of the particle.This way a pixel cluster is characterized by the measured widths w rφ and w z (dimensions of its rectangular envelope), and the charge, which can be ±1 or left unknown.A strip cluster has only one such quantity, its w rφ width.

Particle tracking
The Kalman filter is widely used in particle physics experiments for charged track and vertex finding and fitting and provides a coherent framework to handle known physical effects and measurement uncertainties [13].It is equivalent to a global linear least-squares fit which takes into account all correlations coming from process noise.It is the optimum solution since it minimizes the mean square estimation error.
The state vector x = (k, θ, ψ, rφ, z) is five dimensional, consisting of the signed inverse momentum, local polar angle, local azimuthal angle, global azimuthal position, and global longitudinal position, respectively.The propagation function f (x) from layer to layer is calculated analytically using a helix model.Multiple scattering and energy loss in tracker layers is implemented with their Gaussian approximations shown in Eqs. ( 1)-( 2).The propagation matrix F = ∂ f /∂x is obtained by numerical derivation.The measurement vector for pixel hits m = (rφ, z) is two dimensional, for strip hits m = (rφ) it is one dimensional.The covariance of the process noise is , where σ k = kσ Δ /β, σ θ = σ ψ = θ 0 and F a = ∂ f /∂x a is a vector.Multiple scattering contributes both to the variation of θ and ψ, while energy loss affects only k.The covariance of measurement noise V is a function of σ 2 rφ and σ 2 z , with slightly differing form for pixels and strips, where this latter also incorporates the value of the tilt angle.
Simulated particles are tracked while they are in the volume of the tracker detector, that is, trajectories looping in the magnetic field are properly followed.In the case of pattern recognition and track reconstruction only inside out propagation is considered.

Pattern recognition
Our goal is to collect as much information as possible about potential track candidates, based on the location and shape of the measured hits in an event.To accomplish this, the position of each hit is transformed to a four-dimensional (k T , sinh η, φ 0 , z 0 ) accumulator space of track parameters.The accumulator space is not continuous but binned.(Bins are consecutive, adjacent, non-overlapping equal size intervals of a variable.)Ranges and the optimized number of bins corresponding to track parameters are shown in Table 2.The transformation is a variant of the well known Hough transform [14].In the absence of physical effects (Sec.2.1) the image of a point-like (rφ, z) hit would be a well-defined two-dimensional manifold in that space, while the image of a section-shaped strip hit would be a three-dimensional manifold.

Preparation of templates
The detector models studied here have translational symmetry in longitudinal (z) and rotational symmetry in azimuthal (rφ) direction.The φ − φ 0 angular difference primarily depends on k T , while the z − z 0 longitudinal difference is mostly a function of sinh η (Fig. 2).These difference distributions further depend on the shape of the hit cluster (Sec.2.2).For particles with a given (k T , sinh η) and mass, the (φ − φ 0 , z − z 0 ) values on a given detector layer populate a small rectangular area.The dimensions of that area result from the binning of the track parameters.With help of numerous simulated particles we determine the populated area with help of local linear approximations around the center of each (k T , sinh η) bin.In practice the central values and the Jacobian is deduced for each bin.The set of these values will be referred to as templates in the following.In order to have uniform coverage in all bins, the distribution of simulated particles is chosen to be constant in k T , sinh η, and φ 0 .To limit fluctuations, normal distributed random variables, used in the simulation of physics processes (Sec.2.1), are limited to values within 3.5 standard deviations (only about 0.05% lies outside this range).Altogether 2 × 10 6 pions are generated.
The prepared templates are used in two ways.During the early stage of image transformation they provide a (φ 0 , z 0 ) accumulator area to increment for each (pixel) hit, in the case of a given (k T , sinh η) bin.Later they are used to specify a search rectangle on the (rφ, z) plane of a (strip) layer for a given (k T , sinh η, φ 0 , z 0 ) bin.
Although detector models with only barrel silicon detectors are studied here, the above considerations can be adopted to other geometries, such as disks perpendicular to the beam axis.In that case, the translational symmetry would be lost and the templates would become more complex by introducing another dimension, namely the relative z position of the primary interaction with respect to the longitudinal coordinate of the disk.

Image transformation
The transformation of spatial information to track space proceeds as described in the following.First hits on the three innermost layers are dealt with, containing exclusively pixel hits.For each a hit all potential (k T , sinh η) accumulator bins are examined and the corresponding possible (φ 0 , z 0 ) values, enveloped by rectangles, are determined.Bins within such (k T , sinh η, φ 0 , z 0 ) area are incremented.Since we look for tracks with hits on all three innermost layers, only those accumulator bins are kept which gathered votes from all three layers.
The subsequent layers usually contain strip hits.The method used for the innermost layers would not be efficient here because there are far too many accumulator bins to handle.To this end, the kept accumulator bins are examined, corresponding to proto-tracks with three counts, obtained in the previous step.Using the bin coordinates (k T , sinh η, φ 0 , z 0 ) we look for compatible hits by determining a search rectangle on the (rφ, z) plane for each layer.For quick access, and in order to facilitate hit selection, strip hits are in advance partitioned on an equidistant grid using their (rφ, z) coordinates.
The search for compatible hits proceeds outwards.It is advantageous since the process can be abandoned if some layers provided no compatible hits while they are still reachable according to the curvature range of the examined bin.In other words, the number of layers with compatible hits should not be very different from the number of reachable layers.(We can allow for a few layers without hits.)

Trajectory building
Track candidates are built using hits collected in given accumulator bins.Trajectory propagation and track fitting is performed by the extended classical Kalman filter [13] including prediction, filtering, and smoothing, with pion mass assumption.The initial state vector is estimated by placing a helix to the innermost two hits and using the beam-line as a constraint by adding a zeroth point with value rφ = 0, and with an uncertainty of σ rφ ≈ 50 μm.In the case of an off-centered beam, σ rφ can be increased to properly contain the interaction region in the transverse plane.
Trajectory building normally proceeds from inside out and considers all hit combinations recursively by forming branches.At each layer there are usually multiple hits to add to the existing trajectory.The number of hits to be considered is especially large for the inner strip layers that would exponentially increase the number of trajectory branches.
A powerful solution for this problem is the effective detection of unsuitable (outlier) hits in the busy detector layers.For this purpose, trajectory building starts with the estimate of the initial state vector at the zeroth point, which is then propagated through the first three layers.At that point we x Figure 3. Hits (red boxes and line sections) belonging to a given bin in the accumulator space and trajectories propagated to the outermost detector layer (blue dashed curves).The beamline is indicated with the gray line.
already have a good knowledge about the parameters of the track that is being reconstructed.Instead of going to the next (strip) layer, the trajectory is at once propagated to one of the potential hits in the outermost detector layer (Fig. 3).During propagation the physical effects of the crossed layers are duly taken into account, but information about their hits is not used.The outliers in the intermittent omitted layers are detected next in the smoothing step of the deficient trajectories using the smoothed residual [13].Those (strip) hits in the upper 0.5% tail of the corresponding χ 2 distribution are discarded.
Once the list of compatible strip hits is narrowed down, full trajectory building with the selected hits is performed again.During trajectory building we can allow for a few (one or two) missing hits.It may mean no hit at all or the case when χ 2 is too large for a given number of degrees of freedom (ndf).If there are too many missing hits, the process is abandoned.In order not to lose a noticeable amount of track candidates, but also to have a good selection power, trajectories in the upper 0.5% tail of the corresponding χ 2 distribution are discarded and not developed further.

Optimal distribution of hits among tracks
In the end we have a set of track candidates with the somewhat unusual property that temporarily several track candidates may share some hits.Our goal is to resolve these ambiguities, hit confusion, by optimally allocating the hits among tracks, since all hits must be assigned to not more than one track.The hits and track candidates, and their relations, are best represented by a bipartite graph G.The nodes of G are from two disjoint sets: hits and track candidates, such that each edge connects a hit to a track candidate (Fig. 4).Our goal is to assign all hits to at most one track, while keeping an eye on the total merit of all tracks ( χ 2 ) in the event.
First those track candidates are privileged which are very likely real.To this end, we extract a subgraph selecting track candidates which have at least three hits not requested by other candidates, that is, at least three leaf hit nodes.(A leaf hit node is connected to exactly one track candidate node.)This subgraph is disconnected (Sec.n = 9), and their corresponding hits.Thanks to the high number of hits required, these track candidates are likely real.Selected tracks are again stored, nodes and edges removed from G as above.Then the process is restarted with the subgraph of track candidates with n − 1 hits, iterating down to the subgraph of track candidates with three hits.

Disconnecting a subgraph
The graphs encountered are usually highly connected.If we would try to allocate the hits to tracks one by one, the number of trials needed would explode exponentially with increasing number of nodes.In order to reduce the complexity of the problem the graphs should be partitioned into several small pieces.This can be accomplished by finding vulnerable components in the graph whose removal disconnects the graph.Such vulnerable elements are special edges (bridges) and special nodes (articulation points) whose deletion increases the number of connected components of the graph.Of course this way some tracks would lose a hit, but that is only a small price to pay.
Bridges and articulation points can be found in linear time with help of graph traversal techniques.The depth-first search is an algorithm for traversing and searching a graph.One starts at some arbitrary node as the root and explores as far as possible along each branch before backtracking.The nodes of bridges and the articulation points are found by requiring that their children nodes do not have a backedge.In such a "disconnecting" step bridges and articulation points are removed.process new vulnerable elements may come to light, hence the disconnecting steps are repeated until no new such elements are found.
In addition, the resulted graph is further partitioned into disjoint graphs.This task is best accomplished by the flood fill method embedded into the above detailed traversal technique.The output of the disconnecting step is a large set of disjoint minigraphs.
In high-pileup pp collision events there are usually several thousand track candidates.Their corresponding bipartite graph G and its subgraphs contain several hundred bridges and up to 50 articulation points.Once the subgraphs are disconnected, we get couple of thousand minigraphs.

Solving a minigraph
A minigraph usually has several hits (between 2 and 7 hits for the detector models studied here) with identical role: they are connected to the same set of track candidates.In the interest of reducing complexity hits with identical role are treated jointly, the set of such hits is "contracted" (Fig. 5).
The number of remaining nodes is usually small, the contracted hits can be distributed among tracks by building and solving a decision tree.The process is similar to exploring decision trees of deterministic, strategy board games, such as chess and go, including their horizon problem (limited search depth).What is different here is that our process is a single-player one.The optimal hit-totracks assignments are chosen in the following recursive method: 1. First the most important, highest ranked available (contracted) hit is located.Rank is calculated as the product of the number of similar hits (the number of hits the contracted hit represents) and the number of edges the hit node has (the number of associated track candidates).Such a definition gives preference to contracted hits which represent many particle hits and have a central role in the graph.decreases but stays in the 80-90% range.The fake track rate starts at the permille level and stays under a percent.

Summary
A combination of established data analysis techniques for charged-particle reconstruction was presented.It follows a global approach and uses all information available in a collision event.It employs image transformation based on precomputed templates taking advantage of the translational and ro-  tational symmetries of the detectors.Track candidates and their corresponding hits usually form a highly connected network, a graph.The graph is partitioned into very many minigraphs by removing a few of its vulnerable components, edges and nodes.The hits of a subgraph are distributed among the track candidates by solving a deterministic decision tree.
Tests using simplified computer models of LHC silicon trackers show that efficiency and purity of track reconstruction are excellent and the timing of the proposed method is reasonable, both in simultaneous proton-proton collisions (high pileup), and in single heavy-ion collisions at the highest available energies.

Figure 1 .
Figure 1.Left: Pixel hits (red open squares) and strip hits (red solid sections) of a single inelastic pp collision from simulation of Exp B. The location of the primary interaction is plotted with a green circle, the beam-line is indicated by a gray straight line.Charged particle trajectories (blue dashed curves) are also plotted.Right: Event with 40 simultaneous inelastic pp collisions from simulation of Exp C.

Figure 2 . 5 Connecting
Figure 2. Left: Distributions of φ − φ 0 differences of the hit and track azimuth angles as a function of k T for some selected detector layers.Right: Distributions of z − z 0 differences of the hit and track longitudinal coordinates as a function of sinh η for some selected detector layers.Values from simulation of Exp C are plotted with various symbols, while the oversimplified expectations φ − φ 0 ≈ − arcsin(r/2 • k T ) and z − z 0 ≈ 2 sinh η • arcsin(r/2 • k T )/k T at k T = 0.01 cm −1 are shown with the curves.

Figure 4 .
Figure 4.A small fraction of the bipartite graph G of hits (ellipses) and track candidates (blue diamonds) for a multiple pp collision event (40 simultaneous inelastic pp collisions).Directed arrows, graph edges, show potential hits-to-track candidate assignments.

During the 1 EPJFigure 5 .
Figure 5. Example minigraphs obtained after the removal of all bridges and articulation points from the bipartite graph G of hits and track candidates, in the case of a multiple pp collision event (40 simultaneous inelastic pp collisions).The colors of contracted hits (ellipses) refer the number of hits with identical role (red -6 or more, orange -4 or 5, green -3, white -1 or 2).Filled track candidates (blue diamonds) indicate true tracks, while the others (open diamonds) show candidates where one or more hits are not in place.

Figure 6 .
Figure 6.Performance of the proposed algorithm as a function of transverse momentum p T of the changed particles.Efficiency (top) and fake track rate (bottom) are plotted for Exp A (left) and Exp C (right).The values are separately given for single pp (open green triangles), multiple (40) pp (open red circles), and semi-central PbPb (open blue boxes) collisions.The horizontal scale is linear in the region 0-1 GeV/c, while it is logarithmic for 1-10 GeV/c.

Figure 7 .
Figure 7. Performance of the proposed algorithm as a function of the number of simultaneous inelastic pp collisions.Efficiency (open red circles, top left), running times (open blue diamonds, top right), and fake track rate (open green triangles, bottom) values are plotted.Lines are drawn to guide the eye.

Table 1 .
The main characteristics of the inner barrel silicon detectors of the studied experimental setups.

Table 2 .
Ranges and the optimized number of bins (working point) corresponding to track parameters.The value of p T,min is 0.1 GeV/c.