Jet reconstruction in hadronic collisions by Gaussian filtering
Abstract
A new algorithm for jet finding in hadronic collisions is presented. The algorithm, based on a Gaussian filter in , is specifically intended for use in heavy ion collisions and/or for detectors with limited acceptance. The performance of the algorithm is compared to two conventional algorithms, a seedless cone algorithm and a algorithm, for Pythia simulated dijet events in collisions with . The Gaussian filter is found to perform as well as, and in some instances better than, the conventional algorithms.
pacs:
13.87.aI Introduction
The observation of jet quenching in heavy ion collisions at RHIC Adcox et al. (2001) and the problems inherent in the use of single or dihadron observables to quantify the energy loss of hardscattered partons Adler et al. (2006) together provide a strong motivation for complete jet measurements in heavy ion collisions. However, at RHIC energies, the underlying heavy ion event makes application of conventional jet algorithms difficult. Commonly used jet algorithms employ either fixedsize cones Huth et al. (1990), augmented by split/merging procedures, or iterative clustering (“”) techniques Catani et al. (1993); Ellis and Soper (1993) to determine the solid angle coverage of the jets. The kinematic parameters of the jets are obtained by summing over the particles or calorimeter elements within the jet with a constant (i.e. angleindependent) weight. Such a flat weighting may not be optimal in the presence of a fluctuating background because a typical jet has most of its energy concentrated near the center of the jet and only a small fraction of the energy in the periphery of the jet. In contrast, an angular weighting that enhances the center of the jet compared to the periphery would provide an improved signal to background in the measurement of the jet energy.
Backgrounds from the underlying event in hadron–hadron collisions or heavy ion collisions or from pileup at high luminosity colliders can also distort the jet finding algorithms themselves. High energy particles from the background can prevent or modify the convergence of the meanshift iteration in cone algorithms or can modify the order of combination of fragments/protojets in algorithms. An algorithm that finds maxima in the the angular distribution of fragment (transverse) energies using a weighting function or filter that is strongly peaked has the potential to find the jet position with less sensitivity to the presence of background. Such an algorithm also has the advantage of reducing the impact of the limited acceptance of certain heavy ion detectors used at RHIC and the LHC.
We describe in this paper a new algorithm for finding jets and extracting jet kinematic quantities that uses a linear, Gaussian filter applied in the space of pseudorapidity and azimuthal angle, . Jets are found as local maxima of the filter output. Because the filter is strongly peaked, the algorithm is expected to improve the reconstruction of jet positions and jet energies in the presence of background and to reduce the impact of restricted aperture on jet measurements. This paper takes the first step in exploring the basic characteristics of the Gaussian filter jet finder by studying its performance on Pythia simulated events and comparing the performance to two conventional algorithms, the SISCone Salam and Soyez (2007) seedless cone algorithm and an implementation Butterworth et al. (2003) of the algorithm.
We note that an angular weighting of final state particles has been used previously to define energy flow variables Berger et al. (2003), and the algorithm described here is inspired by that work. Though we focus here on jet finding and energy estimation, as with energy flow variables, the filtered event shape contains more information than just the locations and energies of reconstructed jets, A previous use of lowpass filtering by Donati et al. Donati et al. (1983) should also be recognized. However, that work only used the filter to bootstrap other clustering algorithms, while in this work, filtering is used for complete jet reconstruction.
Ii Algorithm
Instead of the traditional picture of discrete final state particles used by traditional jet reconstruction algorithms, we would like to consider a set of particles of generating the event density
(1) 
Defining the jet reconstruction over a density has the ability to compensate for a nonlocalized background, which is useful to apply the jet reconstruction algorithm to heavy ion collisions and background in high luminosity hadronic colliders.
The jet reconstruction procedure can be expressed as finding the discrete set of all jets with the transverse momentum , pseudorapidity and azimuth ,
(2) 
with the filtered density being the linear–circular convolution
on the cylindrical topology.
When implemented solely with an iterative maximum finding, (2) would have implicitly the issue of proper initialization, not unlike the seeding dilemma of the cone algorithm. However, discrete digital filtering provide an efficient mean to calculate for a large number of sample points. Thus one could find every possible maximum, i.e. in the analogy of the cone algorithm, achieving “seedlessnes” by means of sufficiently sample the entire range with a number of seeds , with being the characteristic separation of close jet pairs, and the pseudorapidity range (typically that of the collider experiment calorimetry).
Therefore we implement (2) by a multistage algorithm, consisting of the following steps:

Accumulate a rectangularly binned density of the event. This can be thought as a histogram, or in term of the density distribution (1) and modulo a constant normalization, to approximate by a discretized for each , with (and denoting the floor function) being the discretized pseudorapidity, and analogously for .

Apply the discrete realization of the filter on the binned density to obtain the initial approximation for the discrete bins or pixels. This filter could be either implemented in the position space using e.g. a finite impulse response (FIR) or infinite impulse response (IIR) version of the filter, or realized in the Fourier space.

The local maxima are localized by comparing the filtered density of each discrete pixel against that of the surrounding pixels.

The set of locally maximal pixel centers is used to initialize a suitable local optimization algorithm, operating on the unbinned density (1), to obtain the true positions.
A jet definition like (2) is inherently collinearly and infrared safe, as is a noniteratively defined quantity not involving any thresholds or cutoffs, and therefore like event shape variables such as thrust, is neither sensitive to infinitely soft radiation, nor collinear splitting. When implemented properly, the maximum finding only depends on the collinearly and infrared safe .
For the filtering kernel, we would like to propose a bivariate Gaussian distribution function with the normalization ,
(3) 
with (and denoting the ceiling function) being the angularly reduced azimuth. The parameter determines the radial scale of the jet reconstruction, analogous to the fixedcone radius . Figure 1 illustrates the Gaussian filter applied to a Pythia event, and a comparison of the filter reconstructed jet distribution functions for jets that the algorithm and SISCone reconstruct to , , and is shown in Figure 2.
The mean shift iteration (such as the cone algorithm) is closely related to a bound maximization with respect to the convolution with a “shadow function” Cheng (1995); Fashing and Tomasi (2005). This has the consequence that for the pure jet direction finding, the filtering based algorithm with subsequent maximization presented here has an exact correspondence in the picture of a weighted mean shift iteration, and therefore can be regarded as a generalized form of the cone algorithm. The correspondence also leads naturally to the definition of jets as the local maxima in (2). For the cone algorithm, the shadow function kernel is . Comparing this with the functional form of (3) suggest for the limit of narrowly focused jet, the angular behavior of the Gaussian filter and the cone algorithm (without split/merge) would correspond with the choice .
There are several reasons to prefer a Gaussian density distribution instead of e.g. the parabolic shadow function implied by the cone algorithm. A Gaussian density does not possess the undesirable feature of maximum creation (or in the language of a cone iteration, stable midpoint axes). And when a slow varying background is present, the foreground to background ratio is accounted for both in the reconstructed jet energy and the directional jet finding, since the Gaussian weighting is present in both and its gradient. As a rapidly decreasing test function, the Gaussian density possesses the mathematical property necessary to regulate the singular directional distribution of an outgoing high fragment in (1).
Gaussian filtering has been used extensively in fields such as computer vision, and several fast approximations are known beside the direct and Fourier space convolution. We use a partial fraction expanded form Hale (2006) of the IIR approximation originally described by Young & van Vliet Young and van Vliet (1995); van Vliet et al. (1998); Young et al. (2002), which compared to a fast Fourier transform (FFT) realization requires no padding along the pseudorapidity axis, but at least two filtering passes to simulate a circular response in azimuth. It should be noted that multipass filtering does not exactly reproduce in the the angularly truncated Gaussian distribution as in (3), but rather the sum of its circular frequency components, and the filtering kernel therefore in fact takes the form
with being the Jacobi theta function Abramowitz and Stegun (1964). The difference is however negligible for all meaningful choices of , especially at the initial, approximate stage of filtering.
As input to the filter we consider final state particles within an acceptance of . The discrete filtering operates on a grid. For the continuous maximization, we use the Newton optimization algorithm with analytically calculated gradient and Hessian, and for robustness, modified to handle of indefinite Hessian by spectral decomposition.
Iii Simulation and event selection
In order to assess the performance of jet reconstruction by filtering, we compared the reconstructed jets against the and cone algorithms at the ideal detector level, by using event generator truth particles from Pythia 6.4.16 Sjöstrand et al. (2006).
For all comparisons we present here, we consider the collision system at , which as of Run8 accounts for approximately of the hadronic RHIC operating mode. The algorithm used is the KtJet Butterworth et al. (2003) implementation in the inclusive clustering mode, with the inclusive stopping parameter and the longitudinally invariant QCD distance scheme Catani et al. (1993). We further compare against SISCone Salam and Soyez (2007), a variant of the cone algorithm, with the overlap threshold parameter and the cone radius , close to the typically used and in asymptotic angular correspondence to the Gaussian filter with for narrow jets. Each of the comparison presented here are based on (triggered) Pythia events.
Two form of triggering were used. For the event multiplicity, dijet and three jet balance (Figures 3, 5, 6, 7), we are comparing two algorithms at a certain scale of hard scattering. Here we triggered on the Pythia hard scattering with a centered, wide window, i.e. with , and obtained the differential cross sections per unit trigger window at the triggered . The same applies to the triggered spectrum that illustrates the sensitivity to soft QCD background (Figure 4). For the comparison of filter reconstructed against the reconstructed of the and cone algorithms (Figure 2), we used a large trigger window with to avoid biasing the jet, and the results are shown as normalized probability density functions.
To test the dijet balance of the filtering algorithm, we arrange the reconstructed jets in descending , i.e. . Events with a prominent dijet structure can be selected by requiring , with as the relative gap, which we chose to be . This approach avoids sensitivity to the individual energy scales of the algorithm. In the pairwise comparision, the selection is applied on both jet algorithms in the logical “and” sense to ensure symmetric event selection (i.e. both algorithms agree that a clear dijet structure is present).
An analogous selection for the three jet events can be made by requiring (with the same choice as with the dijet event selection). If the pairwise angle among the three leading jets are , the resulting three jet opening angle is defined as . Since three algorithms are compared in Figure 7, we used the logical “or” sense (i.e. only one algorithm detecting a clear three jet structure is sufficient).
The inbalance expressed by the three momentum sum scales with the scale the jet algorithm reconstructs to. Therefore the comparison is made against the algorithm with the momentum sum convoluted by the distribution function, i.e. if one would “simulate” the filter behavior by apply the distribution function event by event to the algorithm.
While not explicitly shown, we also checked our results against herwig 5.6.10 Corcella et al. (2001), and found no algorithm specific, dissimilar behaviors, beside differences among the event generators that are reproduced by all algorithms.
Iv Discussion
The prominent feature of jet reconstruction by filtering is the effect of the angular weight. It has the effect that the reconstructed jet is shifted for wide angle fragmentation, while becoming less noticeable for sharply focused, high jets (Figure 2). This accounts for the performance of the filter algorithm in the rejection of the soft background, which usually only contributes significantly to the event when summed over large angles. This inherent discrimination against the background provide us with a good start position to apply the formalism presented here to stronger background levels, including those usually found in heavy ion collisions. And while this effect complicates the determination of jet , the parameter dependent shift of jet and the subsequent need for calibration is present with all jet reconstruction algorithms.
To show the insensitivity to the background, we calculated both the jet multiplicity and jet specturm at initial hard scattering. In term of the jet multiplicity (Figure 3), the filtering reproduces more accurately an initial hard scattering like event shape, which is dominated by two and three jet events.
The spectra (Figure 4) shows a two component mixture, consisting of a high jet cross section, low exponential distribution resulting from the soft QCD background, and a peak from the triggered, initial hard scattering. Here both the algorithm and SISCone in fact produce a large amount of “jets” from background particles, thus creating the large jet multiplicity visible in Figure 3, while the Gaussian filter suppresses the background approximately by a factor of 2 at an initial hard scattering of , thus producing a significantly more isolated and prominent peak for the trigger. This effect is consistent over a range of filter sizes and , and becomes even stronger with increasing filter size or at lower . While SISCone e.g. fails to resolve the hard scattering cross section peak at from the background for any choice of , this is possible with the filter algorithm and , albeit not at .
In dijet events, the Gaussian filter reproduces the angular behavior that previous jet reconstruction algorithms exhibit (Figures 5), and a slight difference in the three momentum sum is visible, with the filter having a slightly weaker exponential tail (Figure 6). The similarities are expected, since with adequate jet reconstruction performance, these variables should be dominated by the physical processes rather than the jet definitions. There is also a slight difference in the three jet opening angle when compared to the algorithm (7), and difference to the algorithm is smaller than the algorithm to SISCone.
V Conclusion
In this paper, we presented a new algorithm for jet reconstruction by filtering, that beside essential pQCD properties of collinear and infrared safety, provides robustness against soft background. Applied on event generator final state particles, we showed that its performance in key areas matches the algorithm and SISCone, such as for the dijet balance and three jet opening angle, which are given by the physical jet production and jet fragmentation and not significantly modified by the jet definition. While the scale differs from the reconstructed scale of the and cone algorithms, we note that an energy calibration is usually required for all algorithms for the extraction of the accurate parton energy. We also demonstrated, in this paper for collisions, its superior ability in discriminating jets from the background, thus rejecting the associated, spurious jets contaminating events at low . We hope to apply the algorithm presented here to most of the collision systems available at RHIC and LHC, thus providing an unified jet definition suitable across a large range of species, collision energies and jet that so far was not attainable with the traditional jet reconstruction algorithms.
References
 Adcox et al. (2001) K. Adcox et al., Phys. Rev. Lett. 88, 022301 (2001).
 Adler et al. (2006) S. S. Adler et al., Phys. Rev. Lett. 97, 052301 (2006).
 Huth et al. (1990) J. E. Huth et al., in Research directions for the decade: proceedings of the 1990 Summer Study on High Energy Physics (1990), pp. 134–136.
 Catani et al. (1993) S. Catani, Y. L. Dokshitzer, M. H. Seymour, and B. R. Webber, Nucl. Phys. B 406, 187 (1993).
 Ellis and Soper (1993) S. D. Ellis and D. E. Soper, Phys. Rev. D 48, 3160 (1993).
 Salam and Soyez (2007) G. P. Salam and G. Soyez, J. High Energy Phys. 2007, 086 (2007).
 Butterworth et al. (2003) J. M. Butterworth, J. P. Couchman, B. E. Cox, and B. M. Waugh, Comput. Phys. Commun. 153, 85 (2003).
 Berger et al. (2003) C. F. Berger, T. Kúcs, and G. Sterman, Phys. Rev. D 68, 014012 (2003).
 Donati et al. (1983) A. Donati, R. Odorico, and V. Roberto, Z. Phys. C 20, 9 (1983).
 Cheng (1995) Y.z. Cheng, IEEE Trans. Pattern Anal. Mach. Intell. 17, 790 (1995).
 Fashing and Tomasi (2005) M. Fashing and C. Tomasi, IEEE Trans. Pattern Anal. Mach. Intell. 27, 471 (2005).
 Hale (2006) D. Hale, Tech. Rep. 546, Center for Wave Phenomena, Colorado School of Mines (2006).
 Young and van Vliet (1995) I. T. Young and L. J. van Vliet, Signal Processing 44, 139 (1995).
 van Vliet et al. (1998) L. J. van Vliet, I. T. Young, and P. W. Verbeek, in Proceedings of the 14th International Conference on Pattern Recognition (IEEE Computer Society, Los Alamitos, CA, 1998), vol. 1, pp. 509–514.
 Young et al. (2002) I. T. Young, L. J. van Vliet, and M. van Ginkel, IEEE Trans. Signal Processing 50, 2798 (2002).
 Abramowitz and Stegun (1964) M. Abramowitz and I. A. Stegun, eds., Handbook of Mathematical Functions (U.S. Government Printing Office, Washington, DC, 1964), p. 576, no. 55 in Applied Mathematics Series.
 Sjöstrand et al. (2006) T. Sjöstrand, S. Mrenna, and P. Skands, J. High Energy Phys. 2006, 026 (2006).
 Corcella et al. (2001) G. Corcella, I. G. Knowles, G. Marchesini, S. Moretti, K. Odagiri, P. Richardson, M. H. Seymour, and B. R. Webber, J. High Energy Phys. 2001, 010 (2001).