Particle-level pileup subtraction for jets and jet shapes

We present an extension to the jet area-based pileup subtraction for both jet kinematics and jet shapes. A particle-level approach is explored whereby the jet constituents are corrected or removed using an extension of the methods currently being employed by the LHC experiments. Several jet shapes and nominal jet radii are used to assess the performance in simulated events with pileup levels equivalent to approximately 30 and 100 interactions per bunch crossing, which are characteristic of both the LHC Run I and Run II conditions. An improved performance in removing the pileup contributions is found when using the new subtraction method. The performance of the new procedure is also compared to other existing methods.


Introduction
The intense environment created by the high luminosity of the Large Hadron Collider (LHC) necessitates novel methods for isolating, mitigating, and, where possible, correcting for the contributions of multiple uncorrelated proton-proton interactions (pileup) to the measured hadronic final state. Pileup has a substantial impact on both jet kinematics and substructure, thereby degrading critical tools for identifying new physics via highly boosted hadronic decays of W , Z, and Higgs bosons, or top quarks.
The last few years has seen the development of several effective tools for pileup mitigation and removal. Simple subtraction techniques remove a constant offset from the measured transverse momentum that is proportional to the number of observed pileup events [1]. So-called grooming techniques such as filtering [2], pruning [3], and trimming [4], actively remove potential pileup constituents from jets. Other approaches, such as the jet cleansing method [5] or charged hadron subtraction [6], use tracking information to identify a given hadronic energy deposition with charged particles originating in pileup interactions.
Techniques that utilize event-by-event and jet-by-jet information to determine the extent of contamination from pileup provide a new approach to perform jet physics at very high luminosities. The area-based subtraction procedure [7] corrects the jet 4-momentum and it is extended to account for hadron masses in Ref. [8]. The shape-expansion method [8] provides general approach to correct jet shapes.
The extension of the area-based subtraction procedure that we propose here allows for a particle-by-particle approach to this concept. We find improved performance in removing the contributions due to pileup using this new procedure even for previously intractable jet shape observables (such as planar flow). Furthermore, this approach provides the possibility to perform pileup subtraction without explicit consideration of a specific jet algorithm, reducing the constraints and increasing the flexibility of the jet area-based subtraction procedure overall. Therefore, this method may be used also in heavy ion physics where the jet reconstruction is challenging due to sizable underlying event [9,10].

Subtraction algorithm
The novel feature of the approach described here is the local subtraction of pileup at the level of individual jet constituents. In contrast to the area-based subtraction and the shape-expansion method, the constituent-level subtraction is performed particle-by-particle, thereby correcting both the 4-momentum of the jet and its substructure, simultaneously. This is achieved by combining the kinematics of particles within a specific jet with the kinematics of soft "negative" particles that are added to balance the pileup contribution.
The basic ingredient of the particle-level subtraction is the pileup energy density estimation which is identical to that used in the shape-expansion method proposed in Ref. [8].
The contamination due to pileup is described in terms of the transverse momentum density ρ and mass density ρ m . The expected pileup deposition in a small region of ∆y∆φ is expressed by the 4-momentum where the pileup p T and mass densities, ρ and ρ m , are assumed to be weakly dependent on rapidity y and azimuth φ. In the shape-expansion method, all particles (or effective particles such as calorimeter towers [12] or clusters [13]) in the event are grouped into patches in order to estimate the densities used in Eq. (2.1). The patches are defined by jets reconstructed using the k t algorithm [14,15]. The transverse momentum, p Tpatch , and mass, m δpatch , of each patch is determined by summing over all particles within that patch: where p Ti and m i are the transverse momentum and mass of particle i, respectively. Each patch covers certain area A patch in the (y − φ) plane. The overall background p T and mass densities are estimated as although several modifications exist, including y-dependent ρ and ρ m [16]. The estimation of the background densities is the first step which needs to be followed by a scheme by which to subtract a specified amount of those densities. In our approach, massless particles with very low momentum are incorporated into the event such that they uniformly cover the y − φ plane with high density. These soft particles are referred to as ghosts and they are most commonly used to define the area of a jet [11] or to perform the shape-expansion correction. Each ghost covers a certain fixed area, A g , in the y − φ plane which is defined by the ghost number density (A g is its inverse). The 4-momentum of each particle or ghost is expressed as where m δ = m 2 + p 2 T − p T (in what follows, we will use superscript g to denote the kinematic variables of ghosts). After adding ghosts into the event, the jet clustering algorithm runs over all particles and ghosts delivering the same jets as in the case without the ghosts. Now, the jets contain except the real particles also ghosts which can be used to correct for the pileup contribution within each jet. Eq. (2.1) is translated into the 4-momentum of each ghost by identifying the transverse momentum p g T and mass m g δ with the amount of pileup within area A g : An iterative procedure is used to define the scheme for calculating the specified amount of transverse momentum and mass m δ to subtract from each jet constituent. For each pair of particle i and ghost k, a matching scheme is implemented using the distance measure, ∆R i,k , defined as (2.6) For complete generality, α is allowed to be any real number, but is taken to be zero in the studies performed here. The list of all distance measures, {∆R i,k }, is sorted from the lowest to the highest values. The pileup removal proceeds iteratively, starting from the particle-ghost pair with the lowest ∆R i,k . At each step, the momentum p T and mass m δ of each particle i and ghost k are modified as follows. (2.7) The azimuth and rapidity of the particles and ghosts remain unchanged. The iterative process is terminated when the end of the sorted list is reached. Alternatively, a threshold ∆R max can be introduced to stop the iterations when ∆R i,k > ∆R max . In principle, introducing the ∆R max threshold also guarantees that only ghosts neighbouring a given particle are used to correct the kinematics of that particle. Particles with zero transverse momentum after the iterative process are discarded and the 4-momentum of a given jet is recalculated following a desired recombination scheme (commonly, the 4-momentum recombination scheme is used [17]). It can happen that after the subtraction no real particle remains. This may be a signal that such a jet originates from pileup. The scalar subtraction in Eq. (2.7) is chosen instead of 4-momentum subtraction since it allows to reduce the local differences between the actual background deposit and its estimate from Eq. (2.1) 1 . The scalar subtraction also eliminates the occurrence of unphysical negative squared masses which may be present when subtracting 4-momenta. Furthermore, the scalar summation is also used in the calculation of the transverse momentum of patches in Eq. (2.2) and thereby it avoids the dependence of background densities in Eq. (2.3) on the shape and size of the patches.
The above described subtraction procedure -referred to as constituent subtraction -corrects both the 4-momentum of a jet as well as its substructure. The constituent subtraction works equally well when applied directly to Monte Carlo truth particles from simulation as when applied to a coarse pseudo-detector grid over which the energy from the truth particles is distributed (see Sec. 3.3). An important feature of the algorithm is that it preserves longitudinal invariance -an arbitrary jet after the correction and a subsequent boost in the direction of colliding beams has the same constituents as the same jet which is first boosted and then corrected. The algorithm can be easily extended to account for rapidity dependence of background densities ρ and ρ m in Eq. (2.5).
It is also straightforward to extend the method further to the whole event instead of correcting just particles within a jet. Jet finding can then be performed using the subtracted event. Global event shapes can also benefit from the correction in addition to individual jet observables. The performance and resolution of missing transverse energy calculated from calorimeter energy deposits may also improve, thereby enhancing several searches for physics beyond the Standard Model. Studies of whole-event constituent subtraction are left to future work, as well as testing directly within the experimental communities.
An important advantage of the constituent subtraction is the speed -it can be as much as twenty times faster compared to the shape-expansion method, depending on the type of the jet shape and the nominal jet radius. Furthermore, the shape-expansion correction must be determined for each jet shape in consideration, whereas the constituent subtraction approach provides a corrected set of constituents, from which any shape may be determined. Corrected constituents may also then be used as inputs to jet grooming and tagging algorithms, e.g. the top-quark tagging using the shower deconstruction method [18]. In comparison to the jet cleansing method or charged hadron subtraction used by the CMS experiment, the constituent subtraction does not require any knowledge about the connection of each charged particle with the signal vertex or pileup vertices, though such a knowledge might in principle be used to further enhance capabilities of the algorithm.
The constituent subtraction procedure has following free parameters: A g , ∆R max , and α. The basic recommended settings are: A g = 0.01, ∆R max → ∞, and α = 0. These settings were used in the performance studies presented in Sec. 3. The subtraction is stable with respect to varying A g . The variation of A g by a factor of two does not lead to a change in any of the studied quantities that would be significant with respect to the statistical uncertainty shown on plots in Sec. 3. Introducing a finite ∆R max may improve the performance of the correction and the speed of the algorithm when running over the full event. The configuration with α > 0 prefers to subtract the lower p T constituents which more often originate from background. This configuration will not be discussed in this paper but it appears to lead to an improvement in the correction of some of the jet shapes.
The software for the constituent subtraction will be implemented as a part of the FastJet Contrib project [19].

Performance of the subtraction
The constituent subtraction algorithm corrects both the jet kinematics (p T and mass) and the jet shapes. Studies of the algorithm performance for p T are discussed in Section 3.1 along with comparisons to the area-based subtraction which follows Ref. [8] where the pileup 4-momentum (2.1) is subtracted using the jet area 4-vector A µ as The performance of the subtraction applied to both the jet mass and several jet shapes is presented in Section 3.3. Comparisons of the constituent subtraction approach with the shape-expansion method [8] are also presented. The studies presented are performed using events generated with PYTHIA 8.180, tune 4C [20,21] but without any detector simulation. The effect of additional proton-proton collisions is emulated by using inclusive events (often referred to as "minimum bias" events) overlaid with the hard scattering interaction, which are also generated with PYTHIA 8.180. The CTEQ 5L, LO parton density functions [22], configured to simulate the LHC conditions at √ s = 8 TeV, are used for all event generation. Two processes are simulated without underlying event: di-jet events covering the p T range of 10-800 GeV and events with boosted top quarks from decay Z → tt of hypothetical boson Z with mass of 1.5 TeV. The performance of the subtraction is tested using jets clustered with the anti-k t algorithm [23] with the distance parameter R = 0.7 or R = 1.0 and jets clustered with Cambridge-Aachen (C/A) algorithm [24] with R = 1.2. These are representative jet definitions for both the ATLAS [25] and CMS [26] experimental collaborations. The number of pileup events, n P U , has a Poisson distribution with a mean n P U . Two pileup conditions are simulated, n P U = 30 and n P U = 100. These pileup configurations represent realistic conditions for LHC Run I and upcoming LHC Run II as well as for the high luminosity LHC running [27]. On average, the pileup contribution to the hard-scatter event can be destribed through mean value of transverse momentum densities, ρ , and pileup fluctuations characterised by standard deviation, σ[ρ]. For the used configuration n P U = 100, these quantities are ρ ≈ 75 GeV and σ[ρ] ≈ 13 GeV.
All jet finding and background estimation is performed using FastJet 3.0.6 [16,28]. The shape-expansion correction is performed using FastJet Contrib 1.003 [19]. The patches in Eq. (2.3) are obtained by clustering particles with the k t algorithm [14,15] with distance parameter R = 0.4. The non-negligible dependence of the background densities ρ and ρ m on rapidity impacts each of the corrections methods discussed below. Consequently, in order to focus the comparisons and performance evaluations, only patches with rapidity |y| < 2.0 are used in Eq. (2.3) and jets are required to be fairly central, with |η| < 2.0.

Jet kinematics
The ability of the subtraction to correctly recover the kinematics of the jet can be characterized in terms of following quantities: jet momentum response, jet momentum resolution, jet position resolution, and jet finding efficiency. These quantities are commonly used to evaluate the performance of the jet reconstruction, see e.g. Refs. [29][30][31][32]. In this section we evaluate these quantities both for the constituent subtraction and area-based subtraction.
The first quantity characterizing the basic performance of the subtraction is the jet momentum response. It can be defined as ∆p T /p orig is the original jet momentum with no pileup and p T is either the pileup corrected jet momentum or the momentum of uncorrected jet that is jet clustered in the presence of pileup with no subtraction. This quantity is also often referred to as the jet energy scale. In the optimal situation, the jet momentum response should be zero which means that, on average, the algorithm can reconstruct the same p T as with no pileup. Left panel of Fig. 1 shows the jet momentum response as a function of the number of pileup collisions, n P U , that is the size of the pileup. The jet momentum response of subtracted jets differs from zero by less then a 1%. The jet momentum response of subtracted jets is stable with respect to the pileup which is a crucial condition for the jet reconstruction. Without satisfying this condition any cut applied on the jet p T would lead to a choice of different subset of jets depending on the size of the pileup. Small deviation from zero of the jet momentum response is resulting from ignoring the rapidity dependence of the pileup density ρ and ρ m and from a small average bias of pileup densities by a presence of the hard-scatter event. Such small deviation can be easily corrected after the jet reconstruction by multiplying the jet momentum by a correction factor. The constituent subtraction performs equally well as the area-based subtraction.
The small deviation of jet momentum response from zero for corrected jets can be contrasted with the jet momentum response of uncorrected jets -in the case of low pileup scenario the momentum response is 10-20%, in the case of high pileup scenario the jet momentum response is around 40%. This means that in the high pileup case, there is on average approximately 40 GeV of the pileup background underneath each jet leading to a reconstruction of a typical 100 GeV jet as a 140 GeV jet if performing no subtraction.
The second quantity characterizing the basic performance of the subtraction is the jet momentum resolution defined as σ[∆p T ]/p orig The jet finding efficiency is defined as the number of original jets having a matching corrected (or uncorrected) jet divided by the number of original jets. The matching criterion is the distance in the η − φ plane between the original jet and corrected (or uncorrected) jet satisfying the condition ∆R = ∆η 2 + ∆φ 2 < 0.2. This quantity is plotted in the lower middle panel of Fig. 1. It shows that in the case of high pileup events it is difficult in principle to reconstruct the jets with p orig T < 50 GeV due to the presence of sizable fluctuating background. The jet efficiency is better for the constituent subtraction than for the area-based subtraction.
The last quantity evaluated is the jet position resolution which characterizes the ability to recover the original jet axis in pseudorapidity, η, or azimuth, φ. The jet position resolution in η and φ are similar. The jet position resolution in η, σ[∆η], is defined as the standard deviation of the difference between the original jet η position and η position of the corrected (or uncorrected) jet. The σ[∆η] is plotted in the right panel of Fig. 1 where a clear difference between jets with pileup contribution and corrected jets is present. Jets corrected by the area-based subtraction have slightly worse jet position resolution than jets corrected by the constituent subtraction.
Based on the analysis of the basic performance, we can conclude the constituent subtraction has a good ability to correct for the pileup background and to recover the original jet kinematics even in the presence of a sizable pileup. The constituent subtraction has generally similar performance as the area-based subtraction. A slightly better performance in terms of jet position resolution and jet efficiency may be attributed to the correction of the jet internal structure which is done by the constituent subtraction. The ability of the constituent subtraction to correct the jet internal structure is discussed in Section 3.3.

Jet shape definitions
The ability of the constituent subtraction method to recover the internal structure of jet has been tested by evaluating four jet shape variables that are discussed in the literature to be useful for analyzing the boosted objects or to perform jet tagging [33,34]. Here we briefly introduce these four jet shapes: • The jet mass which can also be used to identify the hadronic decays of boosted heavy particles [37,38].
• The N -subjettiness, τ N , defined as [35] where R is the distance parameter of the jet algorithm, p Tk is the transverse momentum of constituent k and ∆R ik is the distance between a subjet i and a constituent k.
The N subjets are defined by re-clustering the constituents of the jet with exclusive version of the k t algorithm [15] and requiring that exactly N subjets are found. Beside the N -subjettiness also the subjettiness ratio, τ M N = τ M /τ N , can be used to characterize the jet substructure. Typically, the three-to-two ratio, τ 32 = τ 3 /τ 2 , is used which provides a good discrimination between standard QCD jets and jets formed e.g. by boosted top quarks.
• The k t splitting scale, √ d 12 , defined as [36] where p 1 T and p 2 T are the transverse momenta of two subjets and ∆R 12 is the distance between these two subjets. The two subjets are found by going back one step in the clustering history of the k t algorithm. The variable √ d 12 can be used to distinguish heavy-particle decays, which tend to be reasonably symmetric when the decay is to like-mass particles, from the largely asymmetric splittings that originate from QCD radiation in light-quark or gluon jets.
• The longitudinally invariant version of the planar flow, Pf, defined as [8]: where λ 1 and λ 2 are the eigenvalues of 2 × 2 matrix: where α and β correspond to rapidity or azimuth.

Jet shape subtraction
The constituent subtraction method is tested on various combinations of signal samples, pileup conditions, clustering algorithms, and jet shapes defined in Sec. 3.2. The details of the configuration are provided at the beginning of Sec. 3. The constituent subtraction can recover the original jet shape with a good accuracy in all evaluated combinations. pileup, the distributions corrected by the constituent subtraction and the shape-expansion 3 methods. To quantify precisely the performance of the correction, two quantities have been evaluated for the differences between jet shape x and its original value without pileup x orig : the mean value of these differences ∆x = x − x orig and the standard deviation of these differences σ[∆x] = σ[x − x orig ] which represents the resolution. For each combination of configurations, the uncorrected distributions differ significantly from the corresponding original distribution, and have a significant dependence on n P U . A substantial improvement is achieved by the constituent subtraction. The mean difference ∆x does not exhibit the n P U dependence and it is always centered near zero after the subtraction. The resolution σ[∆x] is improved as well. The constituent subtraction method performs similarly or mostly better when compared to the shape-expansion method.
For any of the studied jet shapes, the shape-expansion method can lead to negative corrected jet shapes that are unphysical. To better visualize the contribution of such values, the first bin with negative jet shape in plots of Figs tion of the jet mass when the corrected energy is smaller than the magnitude of corrected momentum. Again, the negative bin represents the fraction of such jets. The fraction of unphysical jet shapes obtained from the shape-expansion method reaches up to ∼ 12% depending on the p T interval and the type of the jet shape. The constituent subtraction method has been tested also on the jets clustered with C/A algorithm which is often employed in various studies of the jet substructure and boosted objects [33]. The clustering in the C/A algorithm is based purely on the geometry and thus it leads to jet with a different jet area compared to the anti-k t algorithm [11]. The performance of the constituent subtraction method for C/A algorithm with distance parameter R = 1.2 is shown in Fig. 3. For this configuration, the impact of the pileup on jet shapes is much stronger compared to the configuration with the anti-k t algorithm. The constituent subtraction can recover the original distributions and it exhibits significantly better ability to subtract the pileup compared to the shape-expansion method.
Further, the constituent subtraction method has been tested on jets reconstructed in events run through a simple simulation of a segmented detector. In this simulation, the η−φ plane is divided into cells of size 0.1 × 0.1. Particles pointing to the same cell are combined into one new effective particle by summing their energies. The mass of the cell is set to zero and its η − φ position is set to the center of the cell. These new effective particles have the same properties as the calorimeter clusters or towers used in real experiments and they are a combination of the pileup and signal. The jet finding algorithm runs over these events delivering jets that are corrected in the same way as for events composed of standard particles. A typical example of the performance of the subtraction methods in case of this simulation is shown in Fig. 4. The constituent subtraction exhibits very similar performance as without the detector simulation which also applies for the shape-expansion method while the constituent subtraction again outperforms the shape-expansion method.
An important test is to evaluate the jet tagging performance. The splitting scale, √ d 12 , is used to tag the boosted top quarks from the Z decay and tagging efficiencies have been evaluated. The signal sample can be compared with the di-jet sample which in this case provides a reasonable estimate of the background for the Z → tt decay. The result is shown on Fig. 5. One can see that both the constituent subtraction and shape-expansion can achieve the same tagging efficiency as in the case of no pileup.
The above presented results demonstrate the stability and good performance of the constituent subtraction method.

Conclusions
We have introduced a new tool to correct for the pileup in high-luminosity LHC running that represents an extension and a simplification of the current state of the art. The constituent subtraction method operates at the level of the jet constituents and provides both a performance improvement and a simplification compared to existing methods: the precision of the reconstruction of jet shapes is improved as well as the speed of the correction itself.
The constituent subtraction method is tested by evaluating the pileup dependence, and other key metrics, of several jet shapes and jet kinematics using multiple jet definitions. Improvements are demonstrated in both the reduction of fluctuations in the resulting jet shapes as a function of pileup and the ability to remove the pileup dependence of the corrected quantities. It also provides a better jet position resolution and jet finding efficiency for reconstructing jets from boosted objects, which directly impacts the experimental sensitivity to new physics.
Since the correction proceeds without knowledge of a jet algorithm, a novel application of this approach would be to correct the whole event prior to jet finding. Not only would this have implications for the technical performance and computational resources currently devoted to evaluation the pileup subtraction for each jet individually, but it also has the potential to improve the determination of the missing transverse energy. Furthermore, the constituent subtraction approach can be used to correct for the underlying event in heavy ion collisions.
It is very important that these studies be verified by the experiments using both fully simulated data samples of signal and background events, as well as in situ studies using data. Given the excellent correspondence between the similar preliminary studies of earlier pileup correction methods and the experimental reality, we expect that the performance observed here is quantitatively representative of what can be achieved by the experimental collaborations. Nonetheless, any new approach must be vetted and tested thoroughly.