Jet reconstruction at high-energy lepton colliders

In this paper we study the performance in $e^+e^-$ collisions of classical $e^+e^-$ jet reconstruction algorithms, longitudinally invariant algorithms and the recently-proposed Valencia algorithm. The study includes a comparison of perturbative and non-perturbative jet energy corrections and the response under realistic background conditions. Several algorithms are benchmarked with a detailed detector simulation at $\sqrt{s}= 3$~\tev. We find that the classical $e^+e^-$ algorithms, with or without beam jets, have the best response, but are inadequate in environments with non-negligible background. The Valencia algorithm and longitudinally invariant $k_t$ algorithms have a much more robust performance, with a slight advantage for the former.


Introduction
The next large collider facility could be a high-energy electron-positron collider. Linear e + e − colliders are the best tool to explore the energy range from several 100 GeV to a few TeV. The Technical Design Report of the International Linear Collider (ILC [1,2]) project envisages a programme of precision Higgs and top physics at centre-of-mass energies √ s = 250 GeV, 500 GeV and, after an energy upgrade, 1 TeV. The Compact Linear Collider (CLIC [3]) scheme has been shown to reach accelerating gradients that extend the e + e − programme into the multi-TeV regime. CLIC envisages a first phase at √ s = 380 GeV, followed by its high-energy programme of multi-TeV operation, with stages at 1.5 TeV and 3 TeV [4]. A large circular machine, as envisaged by the FCCee [5] and CEPC [6] projects, could provide high luminosity at 250 GeV. FCCee may reach the top quark pair production threshold [5]. On a longer time scale a muon collider [7] also has the potential to reach the multi-TeV regime.
Measurements of hadronic final states are a key ingredient of the programme of any next-generation lepton collider. Excellent jet reconstruction is essential to characterize the couplings of the Higgs boson and top quark at the sub-percent level. To distinguish hadronic W and Z boson decays a jet energy resolution of approximately 3% is required. The linear collider detector concepts [8,9] achieve this performance with highly granular calorimeters [10,11] and particle-flow algorithms [12]. Excellent jet clustering is needed to benefit fully from their potential.
The increase in energy comes with a number of challenges for jet reconstruction. Compared to LEP and SLC, high-energy machines produce an abundance of multi-jet final states, final states with multiple energy scales (in associated production), forward-peaked processes and highly boosted objects. Backgrounds such as γγ → hadrons production are increasingly important at high energy [13]. The classical e + e − algorithms cannot cope with this environment, in particular with the beam-induced background [3,14,15]. A critical evaluation of jet reconstruction at lepton colliders is therefore mandatory.
In this paper we study the performance of jet reconstruction in multi-TeV e + e − collisions in detail. We benchmark the performance of a number of sequential recombination algorithms: the classical e + e − k t algorithm [16] used by the LEP experiments and SLD, the longitudinally invariant k t algorithm [17,18] used at hadron colliders and a generalization of the Valencia algorithm [19], a robust e + e − algorithm.
Several aspects of jet reconstruction performance are studied in simulated events. In a particlelevel study we estimate the size of perturbative and non-perturbative corrections to the jet energy. We establish their dependence on the process and the centre-of-mass energy, and on the parameters of the jet algorithms. These particle-level simulations are also used to study the impact of energy deposited on the signal event by background processes. We conclude by using a realistic simulation of the CLIC detector and particle flow reconstruction to study the jet reconstruction performance in top quark pair production and di-Higgs boson production.
The layout of this paper is as follows. In Section 2 we review the challenges of jet reconstruction at high-energy lepton colliders. In Section 3 the jet reconstruction algorithms are introduced. In Section 4 the perturbative and non-perturbative corrections are studied. In Section 5 we study the response of the algorithms to the signal jet and the background. Section 6 presents the results of a full-simulation study of a few benchmark processes. In Section 7 we discuss possible directions for future work and in Section 8 the most important findings of this work are summarized.
2 Challenges for jet reconstruction at high-energy e + e − colliders The experimental environment at previous lepton colliders, such as LEP and SLC, was very benign compared to that at hadron colliders. While this remains true at future high-energy lepton machines, jet reconstruction faces a number of new challenges.

Multi-jet final states
Future lepton colliders offer the possibility to study 2 → 4, 2 → 6, and even 2 → 8 processes. The dominant branching ratios of the W , Z and Higgs bosons are from hadronic decays. Final states with four jets, most notably e + e − → Zh, play a key role in the physics programme of any future electronpositron collider At high energy final states with six jets (e.g. e + e − → tt), or even eight or ten jets (e.g. e + e − → tth), become important. Imperfect clustering of final-state particles can affect the reconstruction of hadronic decays in an important way.
The impact of incorrect assignments of final-state particles to jets is illustrated with an example. We consider the Higgsstrahlung process with hadronic decays of Zand Higgs bosons, where reliable reconstruction of Z and h-boson candidates is the key to a precise measurement of the Higgs boson couplings [20]. The hard scattering process e + e − → Zh is simulated with the MadGraph_aMC@NLO [21] package and the Z → qq and h → bb decays and subsequent hadronization with Pythia8 [22]. No beam energy spread, initial-state-radiation, or modelling of background or the detector are included. Stable particles are clustered into jets with the Durham algorithm (exclusive clustering with N = 4).
In Figure 1 the invariant mass distribution of the Higgs boson candidates is shown for three centreof-mass energies. We find that the distribution has a non-zero width, even in this relatively perfect simulation. The finite resolution is purely due to imperfect clustering of final-state particles into jets. The effect of confusion in jet clustering is less pronounced at higher centre-of-mass energy, as the greater boost of the Zand Higgs bosons leads to a cleaner separation of the jets.  The most important challenge stems from the larger jet multiplicity. In events with only two jets (i.e. e + e − → Zh events with Z → νν and h → bb) the clustering contribution to the mass resolution is negligible. In the ILC and CLIC analyses of di-Higgs boson [23] and tth production [24,25] jet clustering is found to be the limiting factor for the Higgs mass resolution.
Finally, we note that as the centre-of-mass energy increases, t-channel processes become more important. The most obvious example is vector-boson fusion production of the Higgs boson. The final-state products at high energy are strongly forward-peaked [26] and special care is needed to ensure robust jet reconstruction performance over the full polar angle coverage of the experiment. Di-Higgs boson production through vector boson fusion (e + e − → ννhh) presents a double challenge of high jet multiplicity and forward jets. We therefore take this analysis at √ s = 3 TeV as a benchmark.

Jet substructure
The production of very energetic gauge bosons, Higgs bosons and top quarks with hadronic decays, collectively denoted as boosted objects, poses a challenge to the experiments at a future high-energy collider. Whenever the energy of the decaying object exceeds its mass m significantly the highly collimated decay products typically cannot be resolved. In such cases, an analysis of the internal structure of jets at a scale 2 R < 2m/p T is performed to identify the boosted object [27,28]. At high-energy linear colliders, the separation of boosted W -and Z-bosons and the reconstruction of boosted top quarks challenge the detector and reconstruction algorithms. The highly granular calorimeters [10,11] of the Linear Collider detector concepts and the particle flow paradigm are eminently suited for substructure analyses. We analyze the large-R jet mass resolution in top quark pair production at √ s = 3 TeV, as a first exploration of the jet substructure performance of experiments at future lepton colliders.

Beam-induced background
While the environment at lepton colliders remains much more benign than the pile-up conditions of high-energy hadron colliders, several background sources cannot be ignored in the detector design and evaluation of the performance of the linear collider experiments. The most relevant background source for jet reconstruction at linear e + e − colliders is γγ → hadrons production [3]: photons emitted from the incoming electron and positron beams (bremsstrahlung and beamstrahlung) collide and produce minijets of hadrons 3 . In high-energy colliders the probability to produce a mini-jet event in a given bunch crossing is of the order of one.
To evaluate the detector performance, ILC and CLIC detector concepts superpose a number of γγ → hadrons background events on the hard scattering process. The distribution of the particles formed in γγ collisions is forward-peaked, with approximately constant density per unit of rapidity over the instrumented region (a feature also present in the pile-up due to minimum-bias events in proton-proton collisions). For CLIC at 3 TeV approximately 90% of the energy is deposited in the endcap calorimeters and only 10% in the central (barrel) regions of the experiment [3].
The impact of the background on the performance depends on the bunch structure of the accelerator and the read-out speed of the detector systems. In particular for machines based on radio-frequency cavities operated at room temperature the bunch spacing can be very small; CLIC envisages a bunch spacing of 0.5 ns. Background processes deposit 19 TeV in the detectors during a complete bunch train of 312 consecutive bunch crossings at √ s = 3 TeV [3]. A selection based on the time stamp and transverse momentum of the reconstructed particles reduces this background to approximately 100 GeV on each reconstructed physics event. The relatively large bunch spacing of approximately 500 ns at the ILC allows the detector to distinguish individual bunch crossings. In our full-simulation study simulated γγ → hadrons events are overlaid on the signal events. 2 At hadron colliders the scale R is usually expressed in terms of the ∆R distance between two objects, defined as ∆R = (∆φ ) 2 + (∆η) 2 , where φ and η are the azimuthal angle and pseudo-rapidity 3 Other sources, in particular incoherent e + e − pair production due to beamstrahlung photons, have a non-negligible impact on the design of the innermost detector elements, but can be ignored in the study of the jet reconstruction performance. A detailed discussion is found in Reference [29].

Initial State Radiation
In e + e − annihilation processes the system formed by the final state products is, to first approximation, produced at rest in the laboratory. Initial state radiation (ISR) photons emitted by the incoming electrons and positrons changes this picture somewhat. To estimate the magnitude of the boost introduced by ISR we generate events using a parton-level calculation 4 in MadGraph_aMC@NLO [21]. For each 2 → 2 process e + e − → XY we include also the 2 → 3 process e + e − → XY γ. In Fig. 2 the fraction of the energy carried by the XY system is shown for several 2 → 2 processes and for several centre-of-mass energies. For most s-channel processes, the ISR photon energy spectrum falls off very rapidly and the visible energy distribution displays a sharp peak at 1.  Figure 2: The fraction of the visible energy (the energy carried by all final state products except the photon) in several processes, where a photon radiated off the initial or final-state particles may accompany the final state products. The distributions correspond to pair production of light fermions in association with a photon (e + e − → ff (γ)), W -boson and top-quark pair production (e + e − → W + W − (γ), e + e − → tt(γ)) and the Higgsstrahlung process e + e − → Zh(γ). The centre-of-mass energy is indicated on the figure for each process.
The boost of the system along the z-axis due to ISR remains relatively small close to the production threshold: for e + e − → Zh(γ) at √ s = 250 GeV and e + e − → tt(γ) at 500 GeV β z = v z /c is smaller than 0.1 in over 95% and 90% of the events, respectively 5 . For processes with a cross section that grows with √ s the peak is even narrower. The role of ISR is only significant for radiative return to the Z in the process e + e − → ff (γ), where f is any fermion with mass less than half that of the Z boson.
At linear colliders with very narrow beams the luminosity spectrum displays a sizeable tail towards lower centre-of-mass energy [30]. Beam energy spread and beamstrahlung may therefore lead to a pronounced boost of the visible final state objects in a small fraction of events. This effect is included in the full simulation study of Section 6. 4 The effects of beam energy spread and beamstrahlung are not included. 5 Compared to hadron colliders this boost is very small indeed: for di-jet production at the LHC β z of the di-jet system is very close to 1 and even a massive system such as a top quark pair acquires a typical β z = 0.5.

Jet reconstruction algorithms
In this section the jet algorithms considered in this paper are introduced. We discuss the most important differences and their implications for the performance. Three classes of sequential recombination algorithms are considered here: • The classical e + e − algorithms [16] and their generalization with beam jets [31] • The longitudinally invariant algorithms developed for hadron colliders [17,18] • The Valencia algorithm proposed in a previous publication [19], which is further generalized here.

The VLC algorithm
In Ref. [19] a robust jet reconstruction algorithm was proposed, that maintains a Durham-like distance criterion based on energy and polar angle. It achieves a background resilience that can compete with the longitudinally invariant k t algorithm. Here, we further generalize the definition of the algorithm.
The VLC algorithm has the following inter-particle distance: where R is the radius or resolution parameter. For β =1 the distance is given by the transverse momentum squared of the softer of the two particles relative to the harder one, as in the Durham algorithm. The beam distance of the algorithm is: where θ iB is the angle with respect to the beam axis, i.e. the polar angle. The two parameters β and γ allow independent control of the clustering order and the background resilience 6 . The β and γ parameters are real numbers that can take any value. For β = γ = 1 the expression for the beam distance simplifies to d iB = E 2 sin 2 θ iB = p 2 Ti . We discuss the impact of different choices in Sec 3.3.
This new version of the algorithm fulfils the standard IR-safety tests of the FastJet team. The VLC algorithm is available as a plug-in for the FastJet [31,32] package. The code can be obtained from the "contrib" area [33].

Comparison of the distance criteria
The generalized distance criteria for three families of algorithms are summarized in Tab. 1.
For all algorithms the clustering order can be modified by an appropriate choice of the n in the exponent of the energy (or p T ) in the inter-particle distance (β in the VLC algorithm). The Durham (or k t ) algorithm, with n = 1, clusters pairs of particles starting with those that are soft and collinear (i.e. the inverse of the virtuality-ordered emission during the parton shower). Choosing n = 0 yields the Cambridge/Aachen algorithm, that has a purely angular distance criterion. The anti-k t algorithm has n = −1.
In the generalized algorithms the area of jets is limited. Any particle with a beam distance [35] smaller than the distance to any other particle is associated with the beam jet, and therefore not considered part of the visible final state. The radius parameter R governs the relative size of the interparticle and 6 The first version of the algorithm [19] had a single parameter β . Equation 1 furthermore differs by a factor two from the inter-particle distance of Ref. [19]. To distinguish the two algorithms we refer to the more general expression as the VLC algorithm, while the name Valencia is reserved for the setting β = γ that recovers the first proposal (adjusting R by factor √ 2). Table 1: Summary of the distance criteria used in sequential recombination algorithms. Generalized inter-particle and beam distances are given for three main classes of algorithms: the classical e + e − algorithms (comprising a version with beam jets of the Cambridge [34] and Durham algorithms), the longitudinally invariant algorithms used at hadron colliders, which comprise longitudinally invariant k t , Cambridge-Aachen and anti-k t and the robust e + e − algorithms introduced in Section 3.1.
beam distance and thus determines the size of the jet. This modification renders jet reconstruction very resilient to backgrounds.
The generalized e + e − algorithm and the VLC algorithm have virtually the same inter-particle distance. However, the radius parameter R is redefined: the inter-particle distance d i j denominator is R 2 instead of 1 − cos R. The hadron collider algorithms replace the particle energy E i and angle θ i j with quantities that are invariant under boosts along the beam axis, the transverse momentum p Ti and the dis- where φ is the azimuthal angle in the usual cylindrical coordinates and y denotes the rapidity.  Figure 3: The area or footprint of jets reconstructed with R = 0.5 with the three major families of sequential recombination algorithms. The two shaded areas in each column correspond to a jet in the central detector (θ = π/2) and to a forward jet (θ = 7π/8). The jet axis is indicated with a cross.
Detailed studies [3,14] show that the the longitudinally invariant k t algorithm is much more resilient to the γγ → hadrons background than the classical and generalized e + e − algorithms.
For each of the algorithms the catchment areas of a single central and forward jet with n =1 and R = 0.5 are indicated in Figure 3. The footprint of the central jet (at θ = π/2) is approximately circular for all algorithms. The area of the jet in the forward detector (at θ = 7π/8) shrinks considerably for the longitudinally invariant algorithms and the VLC algorithm. The reduced exposure in this region, where backgrounds are most pronounced, is the crucial feature for the enhanced resilience of these algorithms.
An analytical understanding of this property can be obtained by considering two test particles with energies E i and E j and separated by a fixed angle Ω i j . For the generalized e + e − algorithms, both the distance d i j between the two particles and the ratio d i j /d iB of the inter-particle distance and the beam distance are independent of polar and azimuthal angle. For the longitudinally invariant algorithms the ratio d i j /d iB increases (while the inter-particle distance decreases) as the two-particle system is rotated into the forward region. Finally, for the VLC algorithm the ratio d i j /d iB increases as 1/ sin 2γ θ in the forward region, with a slope that is similar to that of the longitudinally invariant algorithms for γ = 1. The distance d i j is constant, as in classical e + e − algorithms.
A closer comparison of the shape of the footprint of the longitudinally invariant algorithms and the VLC algorithm show that, given identical jet axes, the former extend further into the forward region. This causes a slight difference in background resilience of both classes of algorithms.

Interpolation between algorithms
The two parameters β and γ of the VLC algorithm allow one to tailor the algorithm to a specific application. As these parameters are real numbers, one can interpolate smoothly between different clustering schemes.
The β -parameter that exponentiates the energy in inter-particle and beam distance governs the clustering order (similar to the exponent n in the generalized k t algorithm). For β =1 clustering starts with soft, collinear radiation. Choosing β =0 yields purely angular clustering, while β = −1 corresponds to clustering starting from hard, collinear radiation. These integer choices of β correspond to k t , Cambridge/Aachen and anti-k t clustering, respectively. Non-integer values of β interpolate smoothly between these three schemes.
The parameter γ in the exponent of the beam distance of the VLC algorithm provides a handle to control the shrinking of the jet catchment area in the forward regions of the experiment. After setting the R-parameter to the optimal value for central jets, the area of forward jets can be tuned by the choice of γ to ensure the required background resilience.
We have seen that γ = 1 yields forward jets with a similar size of those of the longitudinally invariant algorithms for hadron colliders 7 . Values of γ greater than 1 further enhance the rise of the d i j d iB ratio in the forward region, causing the jet footprint to shrink faster. Values between 0 and 1 yield a slower decrease of the area when the polar angle goes to 0 or π.
For γ = 0, d iB = E 2β i and we retrieve the generalized e + e − algorithms with constant angular opening: the generalized Cambridge algorithm [34] for β = 0 and generalized k t or Durham [16] for β = 1. Choosing β = -1 yields an e + e − variant of the anti-k t algorithm [36]. A schematic overview of the algorithms in (β , γ) space is given in Figure 4. 7 In both algorithms the ratio d i j /d iB for two test particles separated by a constant angle depends on the polar angle of the system. In the VLC algorithm the ratio is proportional to sin −2γ θ . In longitudinally invariant algorithms the ratio follows the evolution of the pseudo-rapidity and grows approximately as η 2 = log 2 tan θ /2 in the forward region. Both gives rise to qualitative the same behaviour.

Jet energy corrections
Before we turn to a detailed simulation including overlaid backgrounds and a model for the detector response, we study the perturbative and non-perturbative jet energy corrections of the algorithms. Both types of corrections are closely connected to the jet area [37]. In this Section we quantify their impact, following the analysis of Ref. [37]. This first exploration of the stability of the algorithms should be extended in future work to quantify the impact of next-to-leading correction, as performed for instance in Ref. [38]. Also the robustness of the conclusions for a variety of different sets of parameters (tunes) of the Monte Carlo simulation merits further study.

Monte Carlo setup
The Monte Carlo simulation chain uses the MadGraph5_aMC@NLO package [21] to generate the matrix elements of the hard scattering 2 → 2 event. Several processes are studied, but results in this Section focus on e + e − → qq at √ s = 250 GeV and e + e − → tt with fully hadronic top decays at √ s = 3 TeV. The four-vectors of the outgoing quarks are fed into Pythia 8.180 [22], with the default tune to LEP data, that performs the simulation of top quark and W boson decays, the parton shower and hadronization. No detector simulation is performed and initial state radiation and beam energy spread are not included in the simulation. Particles or partons from the Pythia event record are clustered using FastJet 3.0.6 [31] exclusive clustering with N = 2. The default ("E-scheme") recombination algorithm is used to merge (pseudo-) jets.

Definition of response and resolution
The jet energy and mass distributions often display substantial non-Gaussian tails and the choice of robust estimators has non-trivial implications. To estimate the centre of the distribution (i.e. the response) the mean and median are used. In some cases we present both, to give an indication for the skewness of the distribution. The width of the distribution (resolution) is estimated using the inter-quantile range IQR 34 , that measures half the width of the interval centered on the median that contains 68% of all jets. We also use RMS 90 , the root-mean-square of the values after discarding 5% outliers in both the low and high tails of the distribution.

Perturbative corrections
Following Reference [37] we estimate the total energy correction by comparing the parton from the hard scatter to the jet of stable particles. For jets of finite size this correction is dominated by energy that leaks out of the jet. We indeed find that the distribution of the difference of parton and jet energy is asymmetric, with a long tail towards negative corrections, where the parton energy is larger than the energy captured in the jet. This energy leakage is most pronounced for jets with a small radius parameter, as expected.
In Figure 5 the average (dashed line) and median (continuous line) relative energy correction are presented. The left plot corresponds to e + e − → qq collisions at relatively low energy ( √ s = 250 GeV), while the right plot corresponds to e + e − → tt at √ s = 3 TeV. At a quantitative level the results show some dependence on the process, centre-of-mass energy and the generator tune for which they are obtained, but qualitatively the same pattern emerges in all cases. The energy correction decreases as the catchment area of the jet increases. The energy corrections for the generalized e + e − algorithm vanish relatively rapidly, with the median correction reaching sub-% level for R ∼ 1. The VLC and longitudinally invariant k t algorithm show much slower convergence towards zero correction. This is entirely due to jets close to the beam axis. For central jets the three classes of algorithms yield identical results (within the statistical accuracy). The VLC and k t algorithms have similar footprints and, indeed, very similar energy corrections.
The clustering order (as controlled by n in the generalized algorithm and by β in the VLC algorithm has a minor impact on the energy corrections. The (inclusive) Cambridge/Aachen algorithm and anti-k t algorithm give similar results to the k t variants of the same algorithm shown here. The largest part of the jet energy correction due to the finite size is amenable to perturbative calculations. A small residual correction is related to the hadronization and must be extracted from (or tuned to) data. The non-perturbative energy correction is estimated as the difference between the energy of the parton-level jet, clustering all partons before hadronization, and the jet reconstructed from stable finalstate particles. The difference in energy between the parton-level and particle-level jet is typically small, but the distribution is offset from 0 and has a long asymmetric tail. Mean and median are again different and even have opposite signs.

Non-perturbative corrections
The dependence of this correction on R is shown in Figure 6. The non-perturbative part is very small compared to the total correction. It is well below 1% at √ s = 250 GeV, for any value of R studied here. For high-energy collisions the correction is well below the per mille level. The generalized e + e − algorithm again has the best convergence, while for both VLC and longitudinally invariant k t the median or mean remain sizeable even for R = 1.5.

Jet mass corrections
The previous discussion has focussed on the jet energy response. Corrections to other jet properties may also be important. Here, we study the corrections to the jet mass 8 , which can be taken as a proxy for the substructure of the jet. The non-perturbative jet mass correction is defined (analogously to the nonperturbative energy correction) as the difference between the masses of the parton-level and particle-level jet.
The dependence on the radius parameter R is shown in Figure 7. The non-perturbative contribution to the jet mass is quite large. The relative correction can be several tens of % at low energy up to R = 1. It drops to a few % for √ s = 3 TeV. In this case, the algorithms with the e + e − inter-particle distance (generalized Durham and VLC) converge slightly faster than the longitudinally invariant algorithms.

Particle-level results
In this Section the response of several algorithms is studied on simulated e + e − → tt events at √ s = 3 TeV. Clustering is exclusive, with N = 2. Both highly boosted top quarks are reconstructed as a single, large-R jet. We gain insight in the impact of the background by superposing randomly distributed background on the signal events. To this end the Monte Carlo setup described in Section 4.1 is extended with a simple mechanism to superpose a random energy flow on the signal event.

Jet energy response without background
Before we study the impact of the background, the response of the jet algorithms to the signal event is estimated. The energy response is determined as the median reconstructed energy. The mass distribution has a sharp peak at the top quark mass and a long tail towards larger masses. The mass response is therefore estimated as the mean reconstructed jet mass. In both cases the response of the Durham algorithmwhich clusters all final state particles into the jets -is taken as a reference. The reconstructed energy is divided by 1.5 TeV, the reconstructed jet mass by the average jet mass of ∼ 370 GeV.
In Figure 8 the energy and mass response is shown as a function of polar angle for three algorithms: the generalized e + e − k t algorithm (black), the longitudinally invariant k t algorithm (blue dashed) and VLC with β = γ = 1 (red). The R-parameter is set to 1.5 for all three algorithms. The generalized e + e − k t algorithm recovers over 99.9% of the top quark energy for R =1.5, independent of the jet polar angle. The shrinking jet areas in the forward region of the longitudinally invariant k t and VLC algorithms lead to a slightly smaller response for | cos θ | > 0. 6 k t is more pronounced.
The mass response of all three algorithms is substantially lower than for the Durham algorithm. The generalized e + e − k t algorithm has a flat response at nearly 80%. The VLC and longitudinally invariant k t algorithms display the same pattern as for the energy response: VLC starts off with a lower response in the central region, but the response is much flatter versus polar angle.

Jet energy response with background
To gain insight in the performance in a more realistic environment with background, we overlay two hundred 1 GeV particles on each signal event. The background distribution is strongly peaked in the forward direction following an exponential distribution peaked at θ = 0, an approximation to the γγ → hadrons background in energy-frontier electron-positron colliders (a more realistic simulation of this background follows in Section 6.1).
The two event displays in Figure 9 provide a zoom image of the θ − φ plane for a single event. The location of the jet axis is indicated as a red circle. The approximate catchment area of both jets is shown in grey. The green squares represent particles from the top decay that are associated with each jet, the blue squares to background particles clustered into the jet. Both algorithms find a very similar jet axis, centered on the high-energy core of the jet. However, the algorithms have quite distinctive footprints. The longitudinally invariant algorithms expose a larger area in the forward region, which renders it more vulnerable to background in this region.
A quantitative view is obtained by comparing the energy and mass of jets obtained when clustering the same events with and without background particles. The bias (the average difference) in the jet energy and jet mass is shown in Figure 10. The background leads to a significant bias for forward jets reconstructed with the longitudinally invariant k t algorithm. The VLC algorithm, on the other hand, is only affected in the very forward region and the bias is much less pronounced.
The jet mass is known to be quite sensitive to soft and diffuse radiation, with the contribution scaling as the third power of the jet area [39]. We indeed find that the mass is strongly affected. A comparison of the jet reconstruction performance of the same process in a fully realistic environment is presented in

Results from full simulation
The performance of the different algorithms is compared in full-simulation samples. We choose two benchmark scenarios with fully hadronic final states that challenge jet reconstruction: di-Higgs boson production (with h → bb, i.e. a final state with four b-jets) and tt production. Both analyses are performed at √ s =3 TeV, with the CLIC_ILD detector and a realistic background of γγ →hadrons.

Monte Carlo setup
The studies in this Section are performed on CLIC 3 TeV Monte Carlo samples. Events are generated with WHIZARD [40] (version 1.95). The response of the CLIC_ILD detector [8] is simulated with GEANT4 [41]. Multi-peripheral γγ →hadrons events are generated with Pythia and superposed as pileup on the signal events. At CLIC bunches are spaced by 0.5 ns and detector systems are expected to integrate the background of a number of subsequent bunch crossings. In this study, the background corresponding to 60 bunch crossings is overlaid. In the event reconstruction, the information of the tracking system and the calorimeters is combined to form particle-flow objects (PFO) with the Pandora [42] algorithm. Timing cuts on PFOs reduce the background level, with a very small impact on the signal energy flow. The nominal (or default) selection of Ref. [3,14] reduces the 19 TeV of energy deposited in the calorimeters by the entire bunch train to approximately 200 GeV superposed on a reconstructed events. A more stringent set of cuts, referred to as the tight selection in Ref. [3], reduces the background energy by another factor of two. Both scenarios are studied in the following.
The event simulation and reconstruction of the large data samples used in this study was performed using the ILCDIRAC [43,44] grid production tools.

Higgs pair production
The study of Higgs boson pair production is crucial to assess the strength of the Higgs self-coupling. The analysis is very challenging at both hadron and lepton colliders due to the very small cross section. At an e + e − collider, the significance of this signal is enhanced at large centre-of-mass energy, as the production rate in the vector-boson-fusion channel e + e − → ννhh grows strongly with centre-of-mass energy. In this section we focus on events where both Higgs bosons decay to hadrons, through the dominant h → bb decay of the SM Higgs boson. This final state can be isolated [45] provided the four jets are reconstructed with excellent energy resolution.
The challenge of this measurement lies in the fact that both Higgs bosons are typically emitted at small polar angle [26]. The most frequently observed topology has both Higgs bosons emitted in opposite directions: one in the forward direction and the other in the backward direction. At 3 TeV (at least) one of the Higgs bosons is emitted with | cos θ | > 0.9 in approximately 85% of events. In this area of the detector the background level due to γγ → hadrons production is most prominent.
Despite the large centre-of-mass energy, the Higgs bosons are produced with rather moderate energy: in 3 TeV collisions the most probable energy of the Higgs bosons is approximately 200 GeV, with a long, scarcely populated tail extending to 1.5 TeV. The modest Higgs boost is sufficient for the b-quarks to continue in the same hemisphere as their parent Higgs boson, but it is rarely large enough for the Higgs boson to form a single jet.
We perform exlusive jet reconstruction with N jets = 4. Higgs boson candidates are reconstructed by pairing two out of the four jets. The combination is retained that yields the best di-jet masses (i.e. that minimizes χ 2 = (m i j − m h ) 2 + (m kl − m h ) 2 , where m i j and m kl are the the masses of the two di-jet systems and m h = 126 GeV is the nominal Higgs boson mass used in the simulation).
The distribution of the reconstructed mass of both di-jet systems forming the Higgs boson candidates is shown in Figure 11. The results of two algorithms are shown, both with the radius parameter R set to 1.3. The red line denotes the result of the VLC algorithm with β = γ = 1, the black line that of the longitudinally invariant k t algorithm. Numerical results of the centre and width of the reconstructed di-jet mass distribution are presented in Table 2. The response of both algorithms is found to agree to within 0.5%, for all methods to estimate the central value of the distribution. The Higgs mass resolution The dependence of the IQR 34 resolution on the radius parameter and the parameter γ of the VLC algorithm is shown in Fig. 12. The best mass resolution is obtained for large values of R in both algorithms. The choice of R ∼ 1.3 is close to optimal for both algorithms. Variation of the γ parameter, that controls the evolution of the VLC jet area in the forward region, leads to a shift of the optimal value of R. With γ < 1 the jet area is reduced at a slightly slower rate and the best resolution is obtained for smaller R. Choosing γ > 1 the jet area shrinks more rapidly and a larger R is required to capture the complete energy flow.

Top quark pair production
The second benchmark we analyze is pair production of boosted top quarks in multi-TeV operation of the CLIC e + e − collider. At these energies the top quark decay products are so collimated that hadronic top quarks can be reconstructed as a single large-R top-jet (R ∼1). Only the fully hadronic final state e + e − → tt → bbqq q q is considered. Events where either the top or anti-top quark is emitted in the forward or backward direction(| cos θ | > 0.7) are discarded to avoid the incomplete acceptance in that region. To cope with the increased background at 3 TeV the tight PFO selection is applied. Jets are reconstructed with exclusive (N =2) clustering with R =1.2. For comparison the same algorithm is also Table 2: Response and resolution of the di-jet mass distributions obtained with the k t and VLC algorithms in the left panel of Figure 11. The columns list the median di-jet mass, the IQR 34   run on all stable Monte Carlo particles. These include neutrinos, but not the particles from γγ → hadrons background.
The jet energy is a fairly good measure of the top quark energy. The correction to the top quark energy is typically 3.5% and the energy resolution is typically 8%. To measure the performance we compare the jet reconstructed from particle flow objects with the jet found by the same algorithm on the stable particles from the signal event (i.e. excluding the γγ → hadrons). The jet energy residual is defined as the difference of the energy of detector-level and particle-level jets. The distribution is shown in Fig. 13. The response is measured as the median of the residual distribution. The resolution is measured as the IQR 34 .  Figure 14: The reconstructed jet mass distribution for fully hadronic decays of tt events at a 3 TeV CLIC. No backgrounds are added in the left plot. In the right plot 60 bunch crossings of γγ → hadrons background are overlaid on the signal and particle flow objects are selected using the tight selection. Table 3: The bias and resolution of the energy and mass measurements of reconstructed top jets in top quark pair production with fully hadronic top quark decay at a centre-of-mass energy of 3 TeV. Results are presented for the median response and two estimates of the resolution: the 34% inter-quantile range (IQR 34 ) and the RMS of central 90% of jets (RMS 90 ). All results are obtained by comparing the jet energy reconstructed from particle flow objects to the jet of stable MC particles from the signal event.
The performance of the classical e + e − algorithm is such that the figures-of-merit cannot be estimated reliably under nominal background conditions (indicated by "-" entries in the Table) Quantitative results are presented in Table 3. The RMS 90 is also presented to facilitate comparison to other studies. In the absence of background, all algorithms reconstruct the energy of the jet quite precisely, with a bias of less than 1% and a resolution of 2-4%. The performance of the classical e + e − algorithms is degraded as soon as the γγ → hadrons background with tight PFO selection is added. The VLC and longitudinally invariant k t algorithms show very little performance degradation even with the nominal PFO selection.
The jet invariant mass is much more sensitive to soft background contamination [37,39]. The jet mass distributions are presented in Figure 14. In the left panel, which corresponds to tt events without background overlay, all algorithms are seen to reconstruct a narrow peak close to the top quark mass. The long tail toward large mass is due to radiation off the top quark and its decay products and is also  Figure 15: The bias and resolution of the reconstructed jet mass versus radius parameter R of two jet algorithms. The jets are reconstructed in fully hadronic top quark decays in tt events at a center-of-mass energy of 3 TeV. The jets reconstruced from particle-flow objects are compared to jets reconstructed from all stable particles from the signal event. The black curves correspond to the results obtained without background, the blue dashed and red curves to 60 bunch crossings of γγ →hadrons background overlaid on the signal, with the tight PFO selection.
present in the jets reconstructed from stable MC particles. The plots in the right panel show a severe degradation when the γγ → hadrons background with tight PFO selection is added, most noticeably for the Durham algorihtm. The bias and resolution of the jet mass is shown as a function of radius parameter in Fig. 15. A quantitative summary is presented in the second part of Table 3. The bias on the jet mass without background is sub-% for most algorithms. The resolution of the VLC and longitudinally invariant k t algorihms is significantly better than that of the classical e + e − algorithms. The 4.1% resolution is a testimony to the potential of highly granular calorimeters and particle flow reconstruction for jet substructure measurements.
The γγ → hadrons background has a profound effect on the performance. The performance of the classical algorithms is clearly inadequate, with a strong bias and a severe degradation, even with the tight PFO selection. The VLC and longitudinally invariant k t algorihms are much less affected, as expected from the smaller exposed area. The VLC algorithm is found to be more resilient than the longitudinally invariant k t , confirming the result anticipated at the particle level in Section 5.

Discussion
In this Section we discuss the implications for lower-energy colliders and identify several topics that merit further study.

Implications for lower-energy lepton colliders
In this paper we have focused on CLIC operation at √ s = 3 TeV, arguably the most challenging environment that lepton colliders might face in the next decades. We have chosen this environment because subtle differences in jet definitions lead to significant differences in performance. This has helped us to gain a deeper understanding of the intricacies of jet reconstruction at lepton colliders and to establish solid conclusions about the resilience of the different algorithms.
These findings are by no means limited to the CLIC environment. Subtle but significant differences in performance are expected also for the ILC at 250 GeV or 500 GeV and at circular colliders.

Further R&D on jet algorithms
The set of jet algorithms studied in this paper is by no means exhaustive. This study does not address several recent proposals, such as the XCone algorithm [46] or the global jet clustering proposed by Georgi [47].
A broad range of new techniques developed for the LHC have so far remained unexplored. This is particularly true for a set of tools that has proven extremely powerful in pile-up mitigation and correction in ATLAS and CMS.
Jet grooming (the collective name for (mass-drop) filtering [48], pruning [49] and trimming [50]) effectively reduces the exposed jet area to several small regions with large energy flow. This provides an effective means of capturing a large fraction of the jet energy while reducing the impact of soft contamination. Tests of the trimming algorithm in the CLIC environment yield very good results, improving the jet mass resolution significantly.
Techniques to correct for the effect of pile-up based on an event-by-event measurement of the pileup activity [51] are quite successful at the LHC. Subtraction at the constituent level with dynamical thresholds [52][53][54] is under active development. An adaptation to the environment at lepton colliders, with a very sparse background energy flow, may prove useful.

Conclusions
We have studied the jet reconstruction performance of several sequential reconstruction algorithms at high-energy lepton colliders. In addition to the classical e + e − algorithms we include a version of the same algorithms with beam jets. We also study the performance of the longitudinally invariant k t algorithm and a new e + e − algorithm, called VLC [19], which are expected to be more resilient to the impact of backgrounds.
The study is based on detailed Monte Carlo simulation. For two benchmark processes we use a full simulation of the linear collider detector concepts, including the relevant background processes.
The perturbative energy corrections of all algorithms with finite size jets are sizeable for small values of the radius parameter (10-15% for R =0.5) and decrease to 1-5% for R = 1.5. This result is approximately independent of the process and centre-of-mass energy. Convergence with R is faster for the generalized e + e − algorithm than for the longitudinally invariant k t algorithm and VLC, that expose a smaller area to forward jets.
The non-perturbative hadronization correction represents a very small part of the total correction. Its relevance decreases with increasing centre-of-mass energy: it is less than 1% at √ s = 250 GeV for all algorithms and less than a per mille at √ s = 3 TeV. The generalized e + e − algorithm again converges fastest. We have estimated the non-perturbative correction also for the jet invariant mass: these corrections are much larger than the non-perturbative energy corrections and remain of the order of a few % even for R = 1.5 and √ s = 3 TeV. The forward-peaked γγ →hadrons background at future high-energy linear lepton colilders is one of the most important factors in the jet reconstruction performance. It has motivated ILC and CLIC to abandon classical, inclusive algorithms in favour of algorithms with a finite jet size. Algorithms that expose a reduced solid angle in the forward region of the detector, such as longitudinally invariant k t or VLC are more robust. A particle-level study shows that these two algorithms have a different response, with VLC showing a lower response, but one that is more stable versus polar angle. VLC is found to be less susceptible to background.
We present two studies in full simulation, namely di-Higgs production and top quark pair production at √ s = 3 TeV, which present a combination of a relatively harsh background level, high jet multiplicity and forward jets. In both cases the classical e + e − algorithms offer an inadequate performance. The same is true for the generalized version with beam jets. VLC provides significantly better mass resolution for the Higgs study and considerably better jet mass reconstruction than the longitudinally invariant k t algorithm.
Jet clustering is key technique for many analyses of multi-jet final states at future high-energy electron-positron colliders. This study shows that a considerable increase in performance can be obtained by a careful choice of the clustering algorithm. We recommend, therefore, that studies into the physics potential of future e + e − colliders carefully optimize the choice of the jet reconstruction algorithm and its parameters. We also encourage further work on robust algorithms for e + e − collisions.