1 Introduction

The production of hadronic jets is a common feature of particle collisions. Jets are widely studied, as they can be used to test the standard model and measure its parameters, they can signal new physics, and provide important background for new physics searches as well.

Jets are defined through jet clustering algorithms: they take final state particles as an input and combine them according to their prescription into larger objects, what we then call jets. The algorithms have a set of resolution parameters, which defines the jet structure: fixing the values of jet algorithm parameters determines what happens in each step of the clusterization, what particles get combined into jets eventually. Although jet algorithms are required in all kind of jet analysis, one particularly important observable is the so-called jet rate. The jet rate, as a function of its parameters, directly connects to the clustering algorithm, as it provides useful information about how the number of jets depends on the choice of parameters.

Jet rate measures the relative production rate of n-jets compared to all hadronic events. It is given by the ratio of the n-jet cross section \(\sigma _{{\mathrm{n}\text {-}\mathrm{jet}}}\) and the total hadronic cross section \(\sigma _{\mathrm {tot}}\) at center-of-mass energy \(Q^2\):

$$\begin{aligned} R_n(\mathbf {a}) = \frac{\sigma _{{\mathrm{n}\text {-}\mathrm{jet}}}(\mathbf {a})}{\sigma _{\mathrm {tot}}}, \end{aligned}$$
(1)

where \(\mathbf {a}\) denotes the set of jet resolution parameters characteristic to a given jet algorithm. Jet rates are mostly studied as a function of one or more of their resolution parameters. This means clustering the same set of momenta with a wide range of chosen values of the jet resolution parameters. This set of momenta might represent a point in the phase space of final state particles or a physical event, but the actual representation is not important in the scope of the paper. Thereby we use the umbrella term ‘partonic event’.

Since repeated clusterization is usually computationally inefficient, in practice one tries to exploit the properties of the algorithm to enhance performance in computations. This is also important on the theory side, since making higher order predictions in perturbation theory typically requires the generation of millions of phase space points, which all need to be clustered individually. Although the bulk of computational cost is coming from the calculation of amplitudes and subtraction terms, slow clusterization might add a non-negligible time contribution as well.

In the case of \(e^+e^-\) colliders the most common jet clustering algorithm is the \(k_{\bot }\) (or Durham) algorithm [1], which has auspicious properties to do computations efficiently, illustrated in the next section. Today, in the LHC era, the commonly used jet algorithm is the anti-\(k_{\bot }\) algorithm [2]. Though jet rate studies similar to \(k_{\bot }\) ones are not prevalent, they are also hampered by the lack of properties that the \(k_{\bot }\) has. Hence computations can be slowed down significantly due to clusterization only.

In this paper we present a reformulation of the anti-\(k_{\bot }\) algorithm equivalent to the original, which makes it possible to define transition values in a similar fashion to the \(k_{\bot }\) algorithm. Furthermore we show how these transition values can be computed and used to speed up calculations. Our method can be used for both hadron and \(e^+e^-\) colliders.

2 The \(k_{\bot }\) algorithm

We start with a short review of the \(k_{\bot }\) algorithm, and discuss how it is used in calculations in practice. The algorithm depends on a single jet resolution parameter \(y_{\mathrm {cut}}\) and the distance measure is defined as

$$\begin{aligned} y_{ij} = \frac{2 {\mathrm {min}}(E_i^2,E_j^2) (1- \cos \theta _{ij})}{Q^2} . \end{aligned}$$
(2)

\(E_i\) and \(E_j\) denote the energy of particle i and j respectively, while \(\theta _{ij}\) labels the angle between the three-momenta \(\mathbf {p}_i\) and \(\mathbf {p}_j\). During clusterization we compute \(y_{ij}\) for each pair of particles and find the smallest one \(y_{kl}={\mathrm {min}} \, y_{ij}\). If \(y_{kl} < y_{\mathrm {cut}}\) holds we combine particles k and l, then start the procedure again with the new list of objects. Otherwise we stop the clusterization and the resulting objects are considered jets.

In the case of the \(k_{\bot }\) algorithm one can uniquely define transition values. Transition values \(y_{i-1 \leftarrow i}\) are certain values of \(y_{\mathrm {cut}}\), where the number of jets changes from i into \(i-1\) for a given final state configuration. The distribution of the transition value behaves as an event shape observable. Using the \(k_{\bot }\) algorithm every transition value \(y_{i-1 \leftarrow i}\) can be computed performing the clusterization only once independently of \(y_{\mathrm {cut}}\), such that in every clusterization step the smallest \(y_{kl}\) value provides the corresponding \(y_{i-1 \leftarrow i}\) transition value. We repeat the steps until all particles are clustered into two jets. When jets are defined through the \(k_{\bot }\) algorithm the number of jets is a monotonically decreasing function of \(y_{\mathrm {cut}}\) for every possible partonic event.

These two properties of the algorithm described previously make possible to connect the \({\mathrm {d}} \sigma /{\mathrm {d}} y_{i-1 \leftarrow i}\) differential distributions and the \(\sigma _{{\mathrm{n}\text {-}\mathrm{jet}}}(y_{\mathrm {cut}})\) cross section. For example the three-jet cross section can be computed as

$$\begin{aligned} \sigma _{3-jet}(y_{\mathrm {cut}}) = \int _{y_{\mathrm {cut}}}^1 {\mathrm {d}}y_{2 \leftarrow 3} \frac{{\mathrm {d}}\sigma }{{\mathrm {d}}y_{2 \leftarrow 3}} - \int _{y_{\mathrm {cut}}}^1 {\mathrm {d}}y_{3 \leftarrow 4} \frac{{\mathrm {d}}\sigma }{{\mathrm {d}}y_{3 \leftarrow 4}} .\nonumber \\ \end{aligned}$$
(3)

The meaning of the two terms are the following: the three-jet cross section for a chosen \(y_{\mathrm {cut}}\) gets contributions from the \({\mathrm {d}}\sigma /{\mathrm {d}}y_{2 \leftarrow 3}\) differential cross section for every \(y_{2 \leftarrow 3}\) value which is greater than \(y_{\mathrm {cut}}\). This gives the first term in Eq. (3). However the resulting quantity in itself would include all events with \(y_{3 \leftarrow 4} \in [0,1]\). Events with \(y_{3 \leftarrow 4} \in [0,y_{\mathrm {cut}}]\) are indeed events which cluster into three-jets, however for events with \(y_{3 \leftarrow 4} \in [y_{\mathrm {cut}},1]\) clustering stops at four-jets. Thus we need to subtract the integrated \({\mathrm {d}}\sigma /{\mathrm {d}}y_{3 \leftarrow 4}\) distribution to get the correct three-jet cross section, which gives the second term. Equation (3) provides a very useful relation to speed up numerical calculations. One has to perform the clusterization only once per partonic event, calculate the differential cross sections, then do a simple integration with the desired \(y_{\mathrm {cut}}\) according to the formula to obtain the n-jet cross section.

3 The anti-\(k_{\bot }\) algorithm

Now we turn our interest towards the anti-\(k_{\bot }\) algorithm and discuss its shortcomings in computational time compared to the \(k_{\bot }\) algorithm. The anti-\(k_{\bot }\) algorithm uses two different measures: a two-particle measure \(d_{ij}\) and a beam jet measure \(d_{iB}\). They are defined as

$$\begin{aligned} d_{ij}= & {} {\mathrm {min}}(k_{\bot , i}^{2p},k_{\bot , j}^{2p}) \frac{(y_i - y_j)^2 + (\phi _i - \phi _j)^2}{R^2} , \nonumber \\ d_{iB}= & {} k_{\bot , i}^{2p} , \end{aligned}$$
(4)

for hadron colliders, where \(k_{\bot , i}\), \(y_i\) and \(\phi _i\) denote the transverse momentum, rapidity and azimuth of particle i respectively. In the case of \(e^+e^-\) colliders we have

$$\begin{aligned} d_{ij}= & {} {\mathrm {min}}(E_i^{2p},E_j^{2p}) \frac{(1-\cos \theta _{ij})}{1 - \cos R} , \nonumber \\ d_{iB}= & {} E_i^{2p} , \end{aligned}$$
(5)

with the notation being identical to the one introduced in the previous section. Choosing \(p=-1,0,1\) we obtain the anti-\(k_{\bot }\) [2], the Cambridge/Aachen [3] and the inclusive \(k_{\bot }\) [1] algorithms respectively. Collectively they are named as the general inclusive \(k_{\bot }\) algorithm.

The anti-\(k_{\bot }\) algorithm has two jet resolution parameters: R and \(E_{\mathrm {cut}}\). During clustering we calculate \(d_{iB}\) for every particle i and \(d_{ij}\) for every particle pair i, j. If \(d_{kl}\) is the smallest measure, we combine particle k, l, but if \(d_{kB}\) is the smallest one, particle k is considered a jet candidate, and we remove it from the list of objects. We repeat these steps until every particle becomes part of a jet candidate. Finally we apply energy cut(s), and every jet candidate with \(E_i > E_{\mathrm {cut}}\) is a resolved jet.

The anti-\(k_{\bot }\) algorithm has characteristics and properties, which makes it preferable for experimental use [2], for example cone-like jet shapes. However the algorithm has certain other properties, which unfortunately make computational shortcuts like Eq. (3) absent, therefore making clusterization more expensive in the study of the jet rate observable. This is due to the fact that in general the number of jets is not a monotonic function of \(R^2\) or \(1-\cos R\) as it can be seen in Fig. 1. The reason is partially the presence of the additional \(E_{\mathrm {cut}}\) parameter. Although we obtain more and more jet candidates when we increase the spatial resolution, many of them would not survive the last cut on the energy. Furthermore the presence of the beam jet measure, \(d_{iB}\) prevents the same definition of jet transition values as in the case of the \(k_{\bot }\) algorithm.

This would leave us in an unfortunate situation where clustering would need to be done for each different choice of the jet resolution parameters, in particular when we vary R. We note that this still can be an issue, even if one uses the improved version of the anti-\(k_{\bot }\) algorithm [4], which scales \(\mathcal {O}(N \log N)\) compared to the \(\mathcal {O}(N^3)\) cost of the original formulation. The choice of the efficient method depends on the number of histogram bins and the number of partons to be clustered.

Fortunately we can still define jet transition values, which can be used in calculations.

4 Transition values

We start with an equivalent reformulation of the anti-\(k_{\bot }\) algorithm, which is more suitable to define and find transition values. First we combine the two measures \(d_{ij}\) and \(d_{iB}\) the following way

$$\begin{aligned} y_{ijk} \equiv y_{\mathrm {cut}}\frac{{\mathrm {min}}_{i,j} \, d_{ij}}{{\mathrm {min}}_k \, d_{kB}} , \end{aligned}$$
(6)

where we define \(y_{\mathrm {cut}}\equiv R^2\) and \(y_{\mathrm {cut}}\equiv 1 - \cos R\) for hadron and \(e^+e^-\) colliders respectively. Note that \(y_{ijk}\) is independent of \(y_{\mathrm {cut}}\).

Now clusterization is done as it follows: first we calculate \(y_{ijk}\). If \(y_{ijk} < y_{\mathrm {cut}}\), we combine particle i, j; otherwise we consider particle k to be a jet candidate and remove it from the list. We repeat the procedure until the list is empty. Finally we apply the energy cuts on our jet candidates.

The clustering procedure is now similar to the \(k_{\bot }\) algorithm, hence we can define jet transition values in a similar fashion. We call \(y_{t} \equiv y_{\mathrm {cut}}\) a transition value when the clustered particle configuration changes. It is important to notice that it does not necessarily imply a change in the number of jet candidates. Two different \(y_{\mathrm {cut}}\) values can result in the same number of jet candidates, but these candidates may differ in their momenta configuration. As an illustration let us consider the following: we have 4 partons such that they can be separated into two hemispheres and we cluster them into 3 jets. The parton in the first hemisphere is the hardest and is widely separated in angle from the other three, we label it by H. In the second hemisphere one parton is soft and two are hard, labeled as s, h and \(h'\) respectively, with the following angle separation \(1 - \cos \theta _{s h} \sim 1 - \cos \theta _{s h'} > 1 - \cos \theta _{h h'}\). If we choose \(y_{\mathrm {cut}}\) such that \(1 - \cos \theta _{s h}> y_{\mathrm {cut}}> 1 - \cos \theta _{h h'}\), then in the clustering process we first remove the soft parton s from the list, then combine partons h and \(h'\) together and finally obtain \((s),\, (h h'),\, (H)\) as jet candidates. In an other case with \(y_{\mathrm {cut}}> 1 - \cos \theta _{s h}\) we first combine partons s and h, and obtain \((sh), \, (h'),\, (H)\) as jet candidates. Both configurations have the same number of jet candidates, but the way each parton is associated to a jet is different, hence they are in two different regions separated by a transition value. This behavior is due to the presence of the beam jet measure \(d_{iB}\). Then the final number of resolved jets depends on the chosen value of \(E_{\mathrm {cut}}\) as well.

Using this definition is convenient in practice. The transition values must be calculated only once, then one can apply as many different energy cuts as wanted without repeating the clusterization again. Nevertheless the calculation of \(y_{t}\) values is not straightforward. For the \(k_{\bot }\) algorithm the sequence of clustering is independent of \(y_{\mathrm {cut}}\) and relevant information can be fully retrieved for any \(y_{\mathrm {cut}}\) value from one complete clusterization. In contrast, the clusterization sequence of the anti-\(k_{\bot }\) algorithm depends on the actual choice of \(y_{\mathrm {cut}}\), due to the presence of the two different distance measures.

It was shown that in the Cambridge algorithm one faces a similar problem, but transition values can still be found systematically [5]. Here we can adopt the method of Ref. [5] as well to find transition values for the anti-\(k_{\bot }\) algorithm in the following way:

  1. 1.

    First set an initial value for \(y_{\mathrm {ini}}\) and set \(y_{\mathrm {cut}}= y_{\mathrm {ini}}\).

  2. 2.

    If \(y_{\mathrm {cut}}\) is less than some preset lower limit \(y_{\mathrm {stop}}\), stop the algorithm.

  3. 3.

    Perform clusterization with the chosen \(y_{\mathrm {cut}}\), and find the maximum value of \(y_{ijk}\) during the process.

  4. 4.

    Store the transition value \(y_{t} = y_{ijk}^{max}\) and apply energy cuts to obtain the corresponding number of jets.

  5. 5.

    Set \(y_{\mathrm {cut}}= y_{ijk}^{max}\) and go to Step 2.

Clusterization between two transition values is completely determined, choosing two different \(y_{\mathrm {cut}}\) in this set will lead to the same jet configuration. This leads to an improvement in speed in the calculation of jet rates. We can fill histograms more easily between two transition values, we do not have to consider each bin separately and perform clusterization over and over again.

It is worth to mention that the method is independent of the definition of \(d_{ij}\) and \(d_{iB}\) and also independent of how \(d_{ij}\) and \(d_{iB}\) are calculated, given that the maximum value of \(y_{ijk}\) can be obtained during clusterization. Therefore the transition method can be used both in the hadron and \(e^+e^-\) collider version of the anti-\(k_{\bot }\) algorithm and in fact for any version of the general inclusive \(k_{\bot }\) algorithm, that being the original \(\mathcal {O}(N^3)\) or the improved \(\mathcal {O}(N \log N)\) version.

Fig. 1
figure 1

The number of jets as function of \(y_{\mathrm {cut}}\) obtained from clustering a randomly generated partonic event with 10 particles in two different ways. The two approaches provide identical results and the non-monotonic behavior of function is also visible

Full size image

On Fig. 1 we show the number of jets as a function of \(y_{\mathrm {cut}}\). We used a randomly generated partonic event with 10 particles in the final state at \(\sqrt{Q^2} = 100\) GeV center-of-mass energy. \(E_{\mathrm {cut}}\) was chosen 8 GeV. The 10 particle configuration was clustered with the \(e^+e^-\) version of the anti-\(k_{\bot }\) algorithm using both approaches: bin-by-bin with 30 \(y_{\mathrm {cut}}\) values denoted by blue dots and via the transition values method denoted by the red line. Both methods produce identical results, but while the bin-by-bin method required 30 repeated full clusterizations, the red curve was reproduced from 14 transition values. Figure 1 also illustrates the general non-monotonic behavior of the number of jets as a function of \(y_{\mathrm {cut}}\).

5 Performance

Finally we explore the performance of the new method compared to the traditional approach. We employ three different methods: we name the method computing the number of jets over a wide range of \(y_{\mathrm {cut}}\) through transition values as transition, while the bin-by-bin version is dubbed as direct. Both the transition and the direct method are based on the original formulation of anti-\(k_{\bot }\), which scales as \(\mathcal {O}(N^3)\). As third we include the FJcore version of the anti-\(k_{\bot }\) algorithm from the FastJet package [6], and perform bin-by-bin clustering with it. We call this method fjcore. We note that the FJcore package provides only a \(\mathcal {O}(N^2)\) scaling in contrast to the full FastJet version scaling as \(\mathcal {O}(N \log N)\). In exchange FJcore is easier to integrate and is more likely to be used in cases, where the whole apparatus of FastJet is not required. We implemented the transition and the direct methods in a Fortran90 program and included the FJcore algorithm through the provided wrapper. For simplicity we chose the \(e^+e^-\) collider version of the anti-\(k_{\bot }\) algorithm. Using RAMBO [7] we generated 1000 partonic events with 5, 10, 15 and 20 particles in the final state, and clustered them with all three methods. We checked that the three methods give the same results, as it is illustrated in Fig. 1. To perform clusterization with the direct and fjcore methods we selected 30, 60, 90 and 120 bins for \(y_{\mathrm {cut}}\), the first number of bins being closer to experimental setups, while the last one is more typical for theoretical predictions. In the transition method \(y_{\mathrm {ini}}\) was always set to the largest value of \(y_{\mathrm {cut}}\) of the histogram, while \(y_{\mathrm {stop}}\) was chosen to be the smallest. This way we ensured that the range of search for transition values coincides with the range of the histograms.

We summarize our results in Table 1.

Table 1 Required time to perform clusterization of 1000 partonic events using the three different methods. Time values are shown is seconds. Various number of particles and bins were used to illustrate performance behavior
Full size table

The computations were performed on a simple everyday laptop. We emphasize that our numbers in Table 1 are shown just to illustrate the behavior of the new method compared to the usual one, it is not an exhaustive study on performance. For example fluctuations in computational time were not taken into account. Nevertheless Table 1 still provides useful information about how the transition method performs.

As we can see the timing of the direct and the fjcore methods scale with the number of bins, as one would expect it. The numbers indicate a linear relation. The transition method depends non-linearly on the number of particles, as more particles introduce more and more possible final jet configurations, hence more transition values to compute. This method also depends on the range of \(y_{\mathrm {cut}}\) values. Although a large number of particles would mean plenty of transition values, many of them could fall outside of the range of interest, hence they would be not computed in the end. For large number of partons the fjcore method is the best, however there is a turnover at 10 partons, where the transition method starts to take over and for 5 partons it clearly outperforms the other two methods, an order of magnitude speed up can be achieved. It is even more obvious when the number of bins is larger. Interestingly at low multiplicities fjcore is the slowest, which is probably due to the complexity of the algorithm.

Table 1 clearly shows that the transition method can be used to improve the speed of clustering in the calculation of fixed order parton level distributions, like jet rates. The calculation of fixed order predictions typically involve only a small number of strongly interacting final state particles, but a large number of bins in order to produce smooth histogram curves. In addition, millions of phase space points are generated, which all require clusterization, therefore faster methods are preferred.

We note that according to Table 1, the transition method is always faster than the direct method, when the number of bins is large and the multiplicity is moderate or small. As mentioned earlier both methods employ a ‘naive’ \(\mathcal {O}(N^3)\) clusterization. Our new method does not depend on whether the clusterization is done via a ‘naive’ \(\mathcal {O}(N^3)\) or an improved \(\mathcal {O}(N \log N)\) algorithm, given that the value of \(y_{ijk}^{max}\) can be tracked during repeated clustering. Hence we expect the integration of the transition method into the FastJet framework to be possible.

6 Summary

In this paper we defined transition values for the anti-\(k_{\bot }\) algorithm and we presented a way to compute them. The knowledge of these values can speed up computations, which involve large number of variations of the \(y_{\mathrm {cut}}\) jet parameter. Our simple performance test shows that the new method could be applied best to improve performance significantly in the calculation of fixed order predictions for jet rates with the anti-\(k_{\bot }\) algorithm, which might regain interest in the upcoming precision era and future electron-positron colliders. The definition of transition values could also serve as starting point for the development of new observables. Our method can be used both for the hadron and the \(e^+e^-\) collider version of the anti-\(k_{\bot }\) algorithm, in fact for any version of the general inclusive \(k_{\bot }\) algorithm. Furthermore the new transition method can be combined either with the ‘naive’ \(\mathcal {O}(N^3)\) or the improved \(\mathcal {O}(N \log N)\) clusterization method, hence it is compatible with the FastJet framework.