Dispelling the N^3 myth for the Kt jet-finder

At high-energy colliders, jets of hadrons are the observable counterparts of the perturbative concepts of quarks and gluons. Good procedures for identifying jets are central to experimental analyses and comparisons with theory. The Kt family of successive recombination jet finders has been widely advocated because of its conceptual simplicity and flexibility and its unique ability to approximately reconstruct the partonic branching sequence in an event. Until now however, it had been believed that for an ensemble of N particles the algorithmic complexity of the Kt jet finder scaled as N^3, a severe issue in the high multiplicity environments of LHC and heavy-ion colliders. We here show that the computationally complex part of Kt jet-clustering can be reduced to two-dimensional nearest neighbour location for a dynamic set of points. Borrowing techniques developed for this extensively studied problem in computational geometry, Kt jet-finding can then be performed in N ln N time. Code based on these ideas is found to run faster than all other jet finders in current use.


Introduction
Partons (quarks and gluons), are the concepts that are central to discussions of the QCD aspects of high-energy collisions such as those at the Fermilab Tevatron and the future Large Hadron Collider (LHC) at CERN. Quarks and gluons, however, are not observable, and in their place one sees jets, collimated bunches of high-energy hadrons which are the result of the fragmentation and hadronisation of the original hard (high-energy) partons. Today's limited understanding of non-perturbative QCD is such that it is not currently possible to predict the exact patterns of hadrons produced. Instead one makes predictions in terms of quarks and gluons and relates these to observations in terms of hadron jets. Naively, jets are easily identified -one simply searches for bunches of collimated hadrons. However, to carry out accurate comparisons between parton-level predictions and hadron-level observations one needs a well-defined 'jet-finding' procedure. The jet-finder is applied both to perturbatively predicted partonic configurations and to observed hadronic configurations and one then directly compares distributions for the predicted partonic jets and the observed hadronic jets. Though partonic and hadronic jets are not equivalent, there is strong evidence (theoretical [1] and experimental [2]) that the comparison can be performed with controlled accuracy.
Cone jet-finders are the most frequently used at the Tevatron. They are based on identifying energy-flow into cones in (pseudo)rapidity η = − ln tan θ/2 and azimuth φ, together with various steps of iteration, merging and splitting of the cones to obtain the final jets. Cone jet-finders tend to be rather complex, different experiments have used different variants (some of them infrared unsafe), and it is often difficult to know exactly which jet-finder to use in theoretical comparisons.
In contrast, the cluster-type jet-finders, generally based on successive pair-wise recombination of particles, have simple definitions and are all infrared safe (for reviews see [12,13]). We shall focus here on the most widely used of them, the k t jet-finder [5], defined below. Among its physics advantages are (a) that it purposely mimics a walk backwards through the QCD branching sequence, which means that reconstructed jets naturally collect most of the particles radiated from an original hard parton, giving better particle mass measurements [14,15], general kinematic reconstruction [16] and gaps-between-jets identification [17] (of relevance to Higgs searches); and (b) it allows one to decompose a jet into constituent subjets, which is useful for identifying decay products of fast-moving heavy particles (see e.g. [18]) and various QCD studies. This has led to the widespread adoption of the k t jet-finder in the LEP (e + e − collisions) and HERA (ep) communities.
Despite its advantages, k t clustering has so far seen only limited study [19][20][21] at the Tevatron. The reasons for this are not entirely clear. One known drawback of the k t jet finder for high-multiplicity hadron-collider environments is its apparent algorithmic slowness: to cluster N particles into jets requires O (N 3 ) operations in current implementations [22]. For a typical event at the upcoming LHC, with an expected multiplicity of N = O (2000), this translates into a clustering time of O (10 s) of CPU time on a modern O (3 GHz) processor; this is considerable given that the clustering has to be repeated for millions of events. For a typical heavy-ion event at LHC, where N = O (50000), the clustering time would grow to an unsustainable O (10 5 s), i.e. more than one day! Even at the Tevatron, where the multiplicity is quite modest, the fact that noise may cause the number of active calorimeter cells to be far larger than the number of particles has led to the use of a complex (and physically questionable) preclustering procedure prior to running the k t jet finder, so as to reduce the effective value of N to something that is manageable [20].
The slowness of the k t jet-finder has been one of the motivating factors behind proposals for alternative jet-finders [9,10]. Here we will show that the k t jet-finder can in fact be formulated in an algorithmically fast (N ln N) manner. A C++ implementation of this (and a related N 2 ) algorithm 1 will be shown to run orders of magnitude faster than currently available implementations, making it feasible (and easy) to use the k t jet finder for efficiently studying high-multiplicity events.

The k t jet-finder
The k t jet finder, in the longitudinally invariant formulation suitable for hadron colliders, is defined as follows.

For each pair of particles
where k ti , η i and φ i are the transverse momentum, rapidity and azimuth of particle i; for each parton i also work out the beam distance If d min is a d ij merge particles i and j into a single particle, summing their four-momenta (alternative recombination schemes are possible); if it is a d iB then declare particle i to be a final jet and remove it from the list.
3. Repeat from step 1 until no particles are left.
There exist extensions of this basic procedure, (a) where d ij is rescaled relative to d iB by a user-chosen factor 1/R 2 ∼ 1 or (b) where clustering is stopped when all d ij , d iB are above a jet resolution threshold d cut . We here consider only the simplest version, as given above, but the arguments below are identical for the general case. Now we reconsider the above procedure, making explicit the computational overheads of the various steps as implemented in standard jet finding codes [ Step 2 dominates, requiring O ( To obtain a better algorithm we isolate the geometrical aspects of the problem, with the help of the following observation. Lemma: If i, j form the smallest d ij , and k ti < k tj , then R ij < R iℓ for all ℓ = j, i.e. j is the geometrical nearest neighbour of particle i.
Proof: Suppose the Lemma is wrong and that there exists a particle ℓ such that This means that if we can identify each particle's geometrical nearest neighbour (in terms of the geometrical R ij distance), then we need not construct a size- We can therefore write the following algorithm: 1. For each particle i establish its nearest neighbour G i and construct the arrays of the d iG i and d iB .
2. Find the minimal value d min of the d iG i , d iB .
3. Merge or remove the particles corresponding to d min as appropriate. We note, though, that three steps of this algorithm -initial nearest neighbour identification, finding d min at each iteration, and updating the nearest neighbour information at each iteration -bear close resemblance to problems studied in the computer science literature and for which efficient solutions are known: • Given an ensemble of vertices in a plane (specified by the η i and φ i of the particles), to find the nearest neighbour of each vertex one can use a structure known as a Voronoi diagram [23]  Therefore both the geometrical and minimum-finding aspects of the k t jet-finder can be related to known problems whose solutions require O (N ln N) operations.

Timings
The program FastJet 5 codes concrete implementations of the N 2 and N ln N algorithms described above. It has been written in C++ and for the N ln N case makes use of a number of pre-existing components. Construction of a Voronoi diagram is a sufficiently common task (useful in areas of science ranging from astronomy to zoology) that several codes are publicly available. Of these, the only one that we found that also straightforwardly allows the addition and removal of points from a pre-constructed Voronoi diagram, was the fortunately one can show that on average, the number of particles that had i as a nearest neighbour is of O (1). One also needs to establish if any particles acquire the newly created particle ℓ as their nearest neighbour -this can be done in O (N ) time by comparing each particle's current nearest neighbour distance with its distance to ℓ. 5 Available from http://www.lpthe.jussieu.fr/~salam/fastjet. Computational Geometry Algorithms Library (CGAL) [26], in particular its triangulation components [27]. 6 For the binary tree we made use of a (red-black) balanced tree. 7 6 One issue relates to the fact that we need nearest-neighbour location on a cylinder (η-φ space) whereas CGAL works on the Euclidean plane. This problem can solved by making mirror copies of a small (∼ 1/ √ N ) fraction of the vertices near the 0 − 2π border. 7 Balanced trees are the basis of the map and set classes in the C++ Standard Template Library. Figure 2 shows the running times for the two algorithms in FastJet as well as for KtJet, a standard implementation [22] of the N 3 algorithm. Our "N 2 algorithm" actually departs slightly from exact N 2 behaviour owing to certain further optimisations carried out. 8 The scaling with N of the Voronoi-based algorithm has been verified to go as N ln N, as expected. It is the fastest algorithm only for N 10 4 , owing to a large coefficient in front of its N ln N behaviour, mostly associated with the computational geometry tasks. This situation could conceivably be improved in the future by optimisations of the CGAL package or by replacing it with a dedicated implementation of the construction and updating of the Voronoi diagram.
The better of the N 2 and N ln N algorithms (which can be selected based on the value of N) runs at least an order of magnitude faster than the N 3 algorithm for all values of N shown, vastly more at large N. Figure 3 compares the running time of our combined N 2 -N ln N FastJet implementation of the k t jet-finder with other jet-finders whose code is publicly available. One sees that it runs at least an order of magnitude faster than all others (except the almost IR unsafe JetClu).

Perspectives
Since the FastJet algorithm is functionally equivalent to the standard N 3 algorithms used for the k t jet finder till now, the results of the clustering are of course identical to those of other implementations. Howewer, its enhanced speed opens up new ways of using k t clustering in the analysis of high-multiplicity events.
Historically, one apparent drawback of k t -type jets with respect to cone-type jets in hadron-hadron collisions was considered to be the larger fluctuations of the areas of the jets defined by the clustering procedure. Such fluctuations would seem to make it more difficult to subtract, from the hard event, the energy coming from the non-perturbative underlying event and from any additional minimum bias interactions taking place in the same bunch crossing (pileup).
However, the fluctuations become irrelevant if one can easily estimate the area of each individual jet. This can be done on an event-by-event basis, as follows: because of the infrared safety of the k t jet-finder algorithm, one can add a large number of extremely soft particles ("ghosts") to the event without modifying the properties of the hard jets. After clustering, each jet will contain a subset of the ghosts, and if the ghosts had a uniform density in η and φ, then the number of ghosts in a given jet will be a measure of its area. In practice we have found that the use of O (10 4 ) ghost particles is necessary to obtain reliable area estimations. This definition for the area of a k t jet can of course be implemented with any coding of the jet-finder. It is however impractical, indeed nearly impossible, to deal is a widely-used cone-type jet-finder, however it is 'almost infrared unsafe', i.e. perturbative predictions have large logarithmic dependence on small parameters (e.g. seed threshold) [29,30]. MidPoint [29] is an infrared safe cone-type jet finder. For both we use code and parameters from CDF [31]. The Optimal Jet Finder [9] (OJF) has been run with Ω cut = 0.15 and a maximum of 8 jets, so as to produce a final state similar to that returned by the k t and cone jet-finders and to limit its run time.
with the required large number of ghost particles without a fast code. Preliminary studies have shown that with simple assumptions about the uniformity of the underlying event and pileup, one can readily determine its size and subtract it from the hard jets, leading to good determinations of kinematical quantities (e.g. invariant masses) in high-luminosity pp collisions, or of single inclusive jet distributions in Pb-Pb collisions at the LHC. Full results will be shown in [32].
Two more observations are worth making before closing this section. They will both be discussed in more detail in [32].
The first is that it can also be interesting to examine alternative definitions of jet areas. One option is to make use of the areas of the Voronoi cells of all the real particles belonging to a given jet. This definition avoids the need to cluster thousands of ghost particles together with the real ones. It instead rests on the geometrical properties of the event, and on the computational geometry component of the FastJet implementation.
The second observation is that there exist clustering-type jet-finders other than the k t jet-finder that share a large fraction of its features (including, of course, infrared safety), and the possibility of a fast implementation. One such example is the "Cambridge" jet-finder. It was originally formulated in the context of e + e − collisions in [33] and an inclusive version, adapted to hadron collisions, was given in [34]. We shall call this inclusive version the Cambridge/Aachen algorithm. It is defined in the same way as the k t jet-finder at the beginning of Section 2, but replacing the particle-particle distance by d ij = R 2 ij /R 2 , and the particle-beam distance by d iB = 1. We shall show in [32] that the Cambridge/Aachen jet-finder has smaller average areas than the k t jet-finder, making it perhaps even better suited for jet studies in high-multiplicity environments.

Conclusions
To conclude, we have identified an unexpected relation between clustering type jet-finders and a widely studied problem in computational geometry. The resulting reduction of the complexity of the k t jet-finding problem, from N 3 to N ln N, opens up the previously inconceivable option of using the k t jet-finder for the large values of N that arise when considering all cells of a finely segmented calorimeter and for heavy-ion events. For moderate N, the one or two orders of magnitude in speed that we gain with a related N 2 approach pave the way to much wider use of the k t jet finder for normal hadron-collider jet studies, especially at the LHC. More generally, the speed gains discussed in this paper also suggest novel ways of using the k t jet finder, which are the subject of ongoing investigation. One example, given in section 5, is the determination of jet areas, knowledge of which is crucial for optimal subtraction of pileup contamination in high luminosity environments.