Optimal jet radius in kinematic dijet reconstruction

Obtaining a good momentum reconstruction of a jet is a compromise between taking it large enough to catch the perturbative final-state radiation and small enough to avoid too much contamination from the underlying event and initial-state radiation. In this paper, we compute analytically the optimal jet radius for dijet reconstructions and study its scale dependence. We also compare our results with previous Monte-Carlo studies.


Introduction
With the start of the LHC as a main source of motivation, jet physics has seen a tremendous development over the last couple of years. A significant amount of effort has been placed in trying to define jets in an optimal way for new physics searches. Just to mention a few examples, this includes the introduction of new jet algorithms [1,2], new techniques to clean the Underlying Event (UE) contribution to the jets [3,4,5,6], or a large series of sub-jet techniques aimed at tagging boosted objects [3,5,7,8,9,10,11,12].
An aspect that has been slightly less investigated is the values that one should use for the parameters inherent to a jet definition. A noticeable example is the case of the "radius" of a jet, often denoted by R, common to all the standard jet definitions used in hadronic collisions. It has been noticed [4] that if one wants to optimise the kinematic reconstruction of dijet events at the LHC, it is mandatory to adapt the radius of the jet when varying the scale of the process and the type of partons involved in the reconstruction.
With the appearance of new techniques, new parameters are added and one might expect the corresponding multi-dimensional optimisation to become more and more involved. There is therefore a risk that a poor choice of parameters for the definition of jets counteracts the positive effect of these new refinements.
In this paper we thus want to take an alternative, complementary, approach and, instead of adding a new procedure with its own parameters, study how it is possible to fix one parameter in an optimal way given the properties of the events we want to reconstruct. In other words, we shall remove one free parameter by determining it analytically from the scales in the process we shall look at. Practically, we shall focus on the determination of the most natural parameter, the radius R, for the reconstruction of a massive colourneutral object decaying into two jets, as done within a Monte Carlo approach in [4]. In that case, we want to choose R large enough to catch the perturbative QCD radiation, but not too large to avoid an excessive contamination from the UE. Our approach thus goes along the same line as [13], except that we shall extend the discussion to include extra features requested by the fact that we are looking at a peaked distribution and wish to obtain a more precise determination of R. Contrarily to the case of [14], where an attempt was made to adapt dynamically the size of the jet with its hardness by replacing the size parameter with a dimensionful scale varying with the process under consideration, we shall determine a unique value for R from the hard scale of the event and the scale of the UE, i.e. not introducing any alternative parameter 1 .
The paper is organised as follow: in Section 2 we describe the process we will study as well as the details on how we proceed with the event analysis. We also describe in that Section our strategy for the analytic computation of the optimal radius R best . We then proceed with the computation of R best itself. Section 3 concentrates on the situation at the partonic level, where only perturbative radiation has to be taken into account; while in Section 4 we add to the picture the contamination due to the UE. In both cases, we shall compare our results with what we obtain from Monte-Carlo studies, as in [4]. Finally, in Section 5, we shall discuss the extraction of the optimal jet radius for situations where we perform an underlying event background subtraction.
2 Dijet decay of a massive object As in [4], we shall study the hard processes qq → Z ′ → qq and gg → H → gg, where the (fictitious) Z ′ and H are made very narrow. The main advantage of such simple processes is that one can easily study the definition of jets at a given scale, by varying the mass M of the colourless resonance. It also allows one to study the differences between quark and gluon jets, and to test the validity of our computations with Monte-Carlo simulations.
The procedure for the event analysis is straightforward: we cluster the events with a given jet definition, select the two highest-p t jets and reconstruct the heavy object from them. In Monte-Carlo studies, the events were generated either with Pythia [15] (tune DWT) or with Herwig 2 [16] with the default Jimmy tune for the UE. We have required that the jets satisfy p t ≥ 10 GeV and |y| ≤ 5, and we have further imposed that the two hardest jets, the ones used to reconstruct the heavy object, are close enough in rapidity, |∆y| ≤ 1 (so that the transverse momentum of the jets remains close to M/2).
The quantification of the performance of a given jet definition is also borrowed from [4]: we define Q w f =z as the width of the smallest mass window that contains a fraction f = z of the reconstructed objects 3 . With that definition, a smaller Q w f =z would correspond to a narrower peak and thus a better reconstruction quality. Finding the optimal R, R best , is therefore equivalent to finding the minimum of Q w f =z , seen as a function of R. Our task in this paper is to perform an analytic computation of the spectrum of the reconstructed mass peak dP/dm rec , or dP/dδm with δm = m rec −M, the difference between the reconstructed and the nominal mass. By integrating that spectrum, we can compute the probability P (δm 1 , δm 2 ) for the difference δm to be between δm 1 and δm 2 . If we write the quality measure as Q w f =z = q 2 − q 1 , where q 1 and q 2 are the lower and upper end of the mass window defining Q w f =z , it is a straightforward exercise to show that it can be computed from the analytic spectrum by solving dP dδm δm=q 1 = dP dδm δm=q 2 and P (q 1 , q 2 ) = z.
Repeating this for different values of R, one can then compute R best , the value of R for which the quality measure is minimal. We shall use z = 0.25. As far as the analytic computation itself is concerned, there are three physical contributions that will affect the width of the reconstructed peak: final-state radiation, initial-state radiation 4 and the UE. The first two are perturbative; while final-state radiation leads to losses by emissions out of the jet, initial-state radiation adds extra, unwanted, radiation inside the jet. Already at the purely partonic level, there is thus an optimal radius. At this stage, the mass spectrum will depend on the mass of the reconstructed object, M, and whether it decays into quarks or gluons since, perturbatively, gluons radiate more than quarks.
The third contamination, the underlying event, is softer. Like the initial-state radiation, it tends to move the reconstructed mass towards larger values by clustering soft particles into the jet. The amount of contamination involves another scale: the density of UE (per unit area) ρ. The most important source of dispersion in the mass spectrum comes from the fact that ρ varies from an event to another.
In the next two Sections, we first concentrate on the purely perturbative behaviour, then include the effect of the UE. But before getting our hands dirty, there is one additional point that needs to be discussed: the fact that we should also expect some effect due to hadronisation. Actually, for the purposes of this paper, hadronisation plays a limited role. In theory [13], it leads to a loss that increases like 1/R at small R, to be compared with a 3 In [4], Q w f =z was defined as a fraction of the generated objects rather than a fraction of the reconstructed ones. Though the former is practically more reliable, choosing the later does not affect R best in practice and simplifies considerably the analytic computation. 4 Since the produced massive object is colourless, there is no interferences between gluons emitted from the initial-state and final-state partons.
log(1/R) for the perturbative component, and so should dominate the small-R behaviour. However, the normalisation is such that for practical values of R the perturbative contribution dominates, and the effect of hadronisation on the value of R best remains negligible.
A final comment concerns the jet algorithm that we shall consider: we will mostly focus on the anti-k t algorithm [2] as its behaviour is simpler. In Section 4.3 we shall discuss the case of the Cambridge/Aachen [17] and Cambridge/Aachen with filtering [3,4] algorithms.

Perturbative QCD spectrum
The part of the spectrum that is the most naturally computed is the pure perturbative QCD contribution. This is equivalent to concentrating on partonic events. In that case, two types of contribution can affect the reconstruction of the jet momentum: final-state radiation outside of the jet which leads to an underestimation of the jet momentum, and in-jet initial-state radiation yielding an over-estimation of the momentum.
In this Section, we shall thus compute the pure partonic spectrum. We do this at the one-gluon-emission (OGE) level. Since we are interested in the behaviour of the spectrum close to the mass peak, we are dominated by the 1/z singularity directly coming from the soft divergence of QCD. In that region, to get an integrable spectrum around the mass peak, it is also important to resum the (α s log(1/z)) n contributions to all orders, which we shall do assuming a Sudakov-like exponentiation, neglecting non-global logarithms [18,19,20].
For the clarity of the computation, we first deal with the soft approximation, where we only keep the 1/z part of the gluon emission. We then introduce the sub-leading corrections coming from the PDF in the initial state which can have non-negligible effects, especially in the case of gluon resonances at large mass.

Soft-gluon emission approximation
Because of longitudinal boost invariance, we can assume that the heavy object is produced at y = 0. Furthermore, because of our cut |∆y| ≤ 1 on the dijets, we shall assume that the decay products are at the same rapidity. The correction due to the finiteness of ∆y can be computed but has very little effect (at most a few percent). Since the calculation for a nonzero ∆y cut becomes rather technical, we shall keep things simple and assume ∆y = 0. At the lowest order of perturbation theory, the incoming (p 1,2 ) and outgoing (k 1,2 ) partons thus have the following momenta We shall consider the emission of an extra soft gluon: with z ≪ 1. The probability for such an emission to happen is given by the antenna formula dP for an emission from the initial state, and a similar expression with p i replaced by k i for final-state radiation. The pre-factor C R is a colour factor which should be equal to C F for an emission from a quark line and C A for an emission from a gluon line. Note that in the case of a single soft emission, the gluon is caught if and only if it is at a geometric distance smaller than R from one of the two "leading" partons in the final state 5 .
If we denote by δm the difference, m rec − M, between the reconstructed mass and the nominal one, initial-state radiation will lead to an overestimation of the jet mass with δm = zM/2 cosh(y) whenever the additional gluon is radiated inside the jet, while finalstate radiation outside the jet will lead to an underestimation of the mass with δm = −zM/2 cosh(y).
Integrating over the rapidity and azimuth of the emitted gluon and using (2), we find for the initial-state radiation spectrum, and for final-state radiation, with capturing the geometric factors and the R dependence. One recognises the typical expected behaviour: the contamination due to initial-state radiation is proportional to the jet area, i.e. to πR 2 , while the logarithmic behaviour at small R of the final-state radiation is the trace of the collinear divergence in QCD. In both cases, the 1/δm behaviour corresponds to the soft divergence of QCD.
Note that the coefficient of the final-state radiation A f (R) has been expanded in series of R. The logarithmic behaviour at small R corresponds to the leading collinear term as computed in [13] but we have kept enough terms in the expansion to have a correct description over the whole R range. In particular, the numerator in the argument of the logarithm will play a mandatory role to allow the best radius to go above 1, which is the case in the purely perturbative spectrum and, more generally, for large-mass gluon resonances.
It is also important to discuss the choice of scale µ for α s . When using equation (2), the argument of the coupling is, for initial and final-state radiation respectively, Assuming that emissions close to the edge of the jet are dominant, we will use the following simplified approximation: As we are interested in the behaviour of the mass spectrum in the vicinity of the peak, the last thing we need to consider is the virtual corrections and the higher-order terms that diverge as α n s log n (1/z) in that region. In principle, this would depend on the jet algorithm under consideration and an exact computation would involve non-global logarithms [18,19,20]. In practice, we shall adopt a simpler approach and assume that the one-gluon-emission result exponentiates. This is equivalent to the introduction of a Sudakov factor where the integration is performed up to δm max = M/2, corresponding to z = 1 at y = 0. This is an arbitrary choice though a different one would only introduce corrections beyond the level of precision we are working at. Using equations (3), (4), (7) and (8) we find with We note that as long as K i,f is smaller than 1, the spectrum goes to 0 when δm tends to Λ i,f . Alternatively, we could regularise the Landau pole in the expression for α s but we have noticed no significant effect in our conclusions when doing so.  Figure 1: Coefficients 2α a C R A i,f (R)/π 2 , governing the behaviour of the initial or final-state QCD spectrum as a function of R for different masses. The left plot covers the case of quark jets for which we use C R = C F = 4/3, and the right plot is for gluon jets for which Together with the spectrum, we will also need the total probability for the reconstructed mass to be within a certain range, i.e. the integrated spectrum. Denoting by P (δm 1 , δm 2 ) the probability for m rec − M to be between δm 1 and δm 2 , we easily find corresponding to the spectrum (9). Note that in this computation, we have assume that the spectrum is cut at the Landau pole, corresponding to |δm| = Λ i,f . Before proceeding with a more involved computation, it is important to comment on the limit of validity of this perturbative computation. From equations (3), (4) and (5), one might expect higher-order corrections to be important when 2α s C R A i,f (R)/π 2 becomes of order 1. Fixing the scale in the coupling to M/2 for the sake of the argument, we have plotted that quantity on Fig. 1 for various cases of interest. In the case of quark-jets, both K i and K f remain smaller than 1. The case of gluon jets is more interesting; we see that the coefficient for final-state radiation goes above one at small R. In the case of a small-mass resonance, this can happen for values of R as large as 0.5. As a rule of thumb, we can thus say that, in the case of a small-mass resonance decaying to gluon jets, our perturbative computation would only be valid for R 0.5. It is particularly interesting to relate that to the observation made in [4] that in this precise case of small-mass resonances decaying into gluons, we observed a disappearance of the peak in the reconstructed mass spectrum for small values of R. Though this was happening at slightly smaller values of R (around R ≈ 0.4), the present computation provides an explanation of this effect: when 2α a C R A i,f (R)/π 2 grows above 1, the peak in the final-state radiation spectrum (9) is no longer around z = 0 and the mass peak disappears. We shall come back to this point later in Section 4.2.
Note also that, from Fig. 1, one can have a feeling of what the optimal radius will be for parton-level events. Indeed, one can expect that, for R = R best , the spectrum will be nearly symmetric around the nominal mass (i.e. around z = 0), that is K i (R best ) = K f (R best ). For both quarks and gluons and for all masses, this means that one should expect R best ≈ 1.1. We will see later on that this is indeed relatively close to what we shall obtain but that PDF effects, computed in the next Section, induce some departure from that situation.

PDF effects
In order to get a better description of the variety of processes we consider over the whole kinematic range, we shall see that it is important to take into account the effect of the PDF on initial-state radiation. This addresses the fact that the emission of an initial-state gluon imposes to take the parton distribution functions at a larger momentum fraction x. This should therefore introduce an additional suppression, more important at large mass and for gluon-jets than for quark jets.
Let us consider the production of a heavy particle at a fixed rapidity y. The leadingorder cross section is easily obtained: where the index a represents the colour of the colliding partons, f a (x) is the PDF implicitly taken at a scale M, and s is the centre-of-mass energy of the collision. This expression also comes with the kinematic constraint |y| ≤ log(M/ √ s).
We should now add to that process an extra parton emitted at a rapidity η. For our purpose of computing the initial-state radiation contamination to a jet, we can safely make the simplifying assumption that the rapidity of the gluon is the same as the rapidity of the "leading" parton, i.e., η = y. By doing so, we decouple the effects of the PDF from the effects of the jet clustering. The latter will be reinserted later on by multiplying our results by the geometric factor A i (R). We thus have In that expression, ξ is the longitudinal fraction of x 1 carried by the parton entering the collision satisfying 6 1 − ξ ≪ 1. If we rewrite it in terms of the longitudinal momentum of the emitted parton (measured w.r.t. the momentum of the beam) ζ = (1 − ξ)x 1 , one gets after a bit of algebra Keeping the kinematic conventions introduced in the previous section, the longitudinal fraction ζ can be rewritten in terms of the transverse momentum k t = zM/2 ≈ δm of the parton as which finally implies with the kinematic constraint In the limit of small δm, i.e. neglecting the δm offset in the PDF, one gets which is exactly the result we have obtained in the soft limit up to the geometric factor 7 2A i (R) that we will need to reintroduce at the end of the computation.
To simplify the discussion, we shall assume that the PDFs take the form The coefficients N a and β a will depend on the type of parton and the scale M but we shall assume a unique value for λ. At this level, we could proceed by integrating over the rapidity y. This is feasible analytically with the choice of PDF (17) but it leads to a Gauss hypergeometric function in dP (0) i,f /dδm for which the Sudakov factor (8) cannot be computed analytically. We will therefore carry on with the computation at a fixed rapidity y.
Taking into account the emissions from both partons entering the collisions, normalising to the LO cross-section and reinserting the geometric factor, one obtains where we have introduced κ = M/ √ s and τ = 1 + δm/M, and it is understood that the terms in the squared brackets are restricted to the appropriate phase space i.e. κτ e ±y ≤ 1, respectively for each term.
In order to compute the Sudakov form factor analytically, we shall make the approximation that the exponents β a are integers (in practice, we shall use β q = 3 and β g = 5 for (anti-)quark and gluon-jets respectively), and expand the small-x power behaviour to first order in δm: the fact that β a is an integer allows one to see (18) as a polynomial in δm.
We can then easily compute the Sudakov factor using the following fundamental integrals: where Ei(x) is the exponential integral. With these simplifying assumptions, the final perturbative spectrum for initial-state radiation can be written as with the integrated probability and the coefficients Note that to obtain (28), we have to integrate δm up to M/2, which is strictly valid only when |y| ≤ log(κ/2), but similar results can easily be derived at forward rapidities. For the analytic computation, we show both the soft approximation (red curve) and the full case including PDF effects (blue curve). The left plot is for a qq dijet system at a nominal mass of 100 GeV while the right plot is for the gluonic case at 2 TeV. In both cases, the anti-k t algorithm with a radius of 1 has been used for the clustering.

Comparison with Monte Carlo
We want to conclude this Section by a comparison of the predictions we obtain for the optimal radius from our analytical studies and from Monte-Carlo approaches. When dealing with the analytic result, we shall use the initial and final-state results computed in the previous sections in the two different situations for which we have analytic results: the soft limit with a running-coupling approximation, eqs. (9) and (11), and the "full" spectrum where we also include the PDF effects in the initial-state radiation, eq. (23). The total perturbative QCD spectrum is the convolution of the initial and final-state pieces which, in practice, will be carried numerically. We shall use Λ QCD = 200 MeV and N f = 5 to compute the running of the QCD coupling, and, for the PDF effects in the initial-state-radiation spectrum, we shall adopt λ = 0.3, β q = 3 and β g = 5, and assume y = 0, where we have checked that most of the massive objects in the Monte Carlo sample were produced. The quality measure is computed from the spectrum using (1).
We start our comparison directly with the spectrum for the reconstructed dijet mass, as shown on Fig. 2 for quark jets, M = 100 GeV and gluon jets, M = 2 TeV. Generally speaking, the agreement between our analytic computation and Herwig is good, even very good for the 100 GeV quark-jet case. Note that this agreement is always better in the region of the peak, where the soft-gluon emissions approximation we are working with is supposed to hold. Similarly the small differences in the tail at small and large m rec are beyond the scope of our approximation. In both the quark and gluon-jet cases, we see that the inclusion of PDF effects, reducing the initial-state radiation, significantly improves Next, we turn to the computation of the quality measure, presented in Fig. 3. After properly taking into account the PDF effects, our computation agrees in shape with the Monte-Carlo results. For the normalisation, we are closer to the Herwig simulation in the case of quark jets, while we better agree with Pythia simulations for gluon jets. Since our analytic computation is identical in both cases, it it somehow hard to find an explanation for those differences.
Finally, let us concentrate on the values of R best extracted from the previous examples. This is shown on Fig. 4 for quark jets (left panel) and gluon jets (right panel). Despite the apparent normalisation issue between Pythia and Herwig observed on Fig. 3, the fact that they agree on the shape of the quality means that they yield very similar extracted optimal radii. If we turn to the analytic computation, we first see that in the soft approximation, the optimal radius is independent on the process under consideration with R best ≈ 0.96 regardless of the mass or the parton type 8 . The reduction of the initial-state radiation due to the PDF effects has a sizeable effect on the extracted optimal radius in such a way that, at the end of the day, our analytic extraction of R best in the "full" model reproduces very well the Monte-Carlo simulations with an error 0.1.

Description of the Underlying Event
One of the most important effects when clustering jets at hadronic colliders is the presence of the underlying event (UE) coming from multiple interaction and beam remnants. This tends to produce a fairly uniform soft radiation of a few GeV per unit area that contaminates the jets and affects their momentum reconstruction. In this Section, we therefore study the effects of this UE on the mass spectrum we have computed from perturbative QCD in the previous Section. If we work within the assumption that the UE is uniform in rapidity and azimuth, the average p t contamination to each jet will be A jet ρ where A jet is the (active) area of the jet [22] and ρ the UE density per unit area. For a given event, this contamination is affected by two kinds of fluctuations: first, the fluctuations of the background itself which go like A jet σ, with σ the background fluctuations per unit area. Then, keeping in mind the possibility to discuss our results for different jet algorithms than the anti-k t algorithm, the active area of a jet is in general known to fluctuate, with an average value proportional to πR 2 that increases with p t in a different way for each jet algorithm.
For a given background density ρ, let us denote by dP i /dp t,k the probability distribution for a jet k to receive a contamination p t,k from the UE. For a jet of area A jet = πR 2 a k , this distribution has an average of πR 2 a k ρ and a dispersion √ πR 2 a k σ. We therefore also have to introduce dP a /da k the distribution, of averageā and dispersion σ a of the area divided by πR 2 , of the jet k.
If the two jets were perfect circles of radius R, one could easily compute that the offset on the reconstructed mass would be 9 δm = 2πRI 1 (R)(ρ 1 + ρ 2 ), up to corrections proportional to δm 2 /M, where I 1 (x) is the modified Bessel function of the first kind. For generic jets we shall then use where p t,1 and p t,2 are the transverse momenta of the background contamination of our dijets.
On top of that, when we consider an event sample, the values of ρ and σ will fluctuate from one event to another. For the event-by-event fluctuations of the background densities, we shall assume a distribution dP e /dρ, of average ρ and dispersion σ ρ . For the intra-event fluctuations (per unit area) σ we shall simply assume that they remain proportional 10 to ρ, i.e., µ = σ/ρ = σ/ρ . Given the event-by-event ρ distribution dP e /dρ, the intra-event fluctuations, dP i /dp t,k and the area distribution dP a /da k , the UE spectrum can be inferred from (29) dP UE dδm = dρ dP e dρ da 1 da 2 dP a da 1 dP a da 2 dp t,1 dp t,2 dP i dp t,1 dP i dp t,2 δ δm − 2I 1 (R) R (p t,1 + p t,2 ) .
(30) Without knowing exactly that distribution, we can compute its average and dispersion in terms of the properties of the fundamental distributions: The average is consistent with what one would naively expect, and the dispersion is a combination of the three separate sources of dispersion we are including: event-by-event fluctuations in ρ (the first term), area fluctuations (the second term) and intra-event fluctuations (the third term). As expected, the first two are proportional to the area of the jet i.e. to πR 2 , while the intra-event fluctuations are proportional to its square root.
So, instead of modelling the various distributions separately and perform the integrations in (30), we shall directly model the total UE spectrum and adjust its parameters to reproduce the total average and dispersion. More specifically, we shall assume The parameters δm 0 and α can be adjusted to reproduce the average (31) and dispersion (32) of the UE distribution. One finds This choice has the advantage over a Gaussian distribution that it remains positive.

Computation and comparison with Monte Carlo
Our strategy is rather straightforward: as we did when combining the initial-state and final-state radiation spectra, we shall convolute the perturbative QCD spectrum, discussed in Section 3.3, with the UE spectrum (33). For each mass and parton-type, we wish to compare our results with the Monte-Carlo studies [4]. For the perturbative QCD part of the spectrum, we use the convolution already discussed in Section 3.3 and consider the same models: the pure soft approximation and the spectrum including PDF effects in the initial-state. The result is itself convoluted with the underlying-event distribution for which we still have to determine the parameters appearing in eqs. (31) and (32). The area propertiesā and σ a have been calculated, together with their scaling violations in [22]. The properties of the background are computed directly from the Monte-Carlo event sample using FastJet [24] and the method advocated in [25] and [23]: each event is first clustered with the Cambridge/Aachen algorithm with R = 0.6; ρ and σ are then estimated from the jets that are within a rapidity strip of 1 unit around the hardest dijet system, excluding the two hardest jets in the event. The average properties ρ , µ and σ ρ can then be easily extracted from the complete event sample. For completeness, we give a table of the values we have obtained in Appendix A.
With that method, we are able to compute the mass spectrum, and thus the quality measure, at a given mass and with a given parton type and UE characteristics. We then extract R best by searching the minimum of the quality measure by steps of 0.01 in R.
As for the parton-level case, we shall compare our computations with the Monte Carlo simulations from Pythia (tune DWT) and Herwig (with default Jimmy tune) as the details of the Underlying Event can differ. We shall again proceed in three steps: first consider the histograms, then the quality measures and finally the optimal radius. For simplicity, we shall focus on the anti-k t algorithm for whichā = 1 and σ 2 a = 0. Let us start with the histograms, Figure 5, where we illustrate our description of the spectrum for the quark dijets at small mass and gluon dijets at large mass. In both cases, our analytic histogram correctly reproduces the behaviour observed with Herwig. We note  however that our analytic computation slightly underestimating the UE in the peak region in the case of quark jets at small nominal mass 11 M. If we then move to the case of the quality measure, Figure 6, we see that, given the differences we have already observed at the parton level (see Figure 3), the addition of the underlying event gives a reasonable description of the Monte-Carlo results, especially the shape of the quality at large R. Note that only the curves corresponding to the model including the PDF effects are shown on Figure 6 for clarity reasons. Finally, the extracted optimal radius is presented on Figures 7 and 8 for Pythia and Herwig simulations respectively, as a function of the mass for both quark and gluon jets. Generally speaking we obtain a good description of the Monte-Carlo results. We see once again that the inclusion of the PDF effects significantly improves the description of the optimal radius for gluon jets at large mass, as expected. We slightly underestimate the optimal radius in the case of gluon jets at large mass, but we have to notice that this is the regime on which it matters the least as in that region the minimum in the quality measure is rather flat.
We thus see that our complete analytic model, the solid (blue) curve on Figures 7 and 8, corresponding to our analytic pQCD computation including PDF effects, eq. (23), convoluted with our model of the Underlying Event, eq. (33), gives a good extraction of the optimal radius R best .

Computation for other algorithms
Now that we have checked the agreement between our analytic computations and Monte-Carlo simulations for the simpler case of the anti-k t algorithm, let us briefly discuss other possible jet algorithms. In this section, we shall apply our analytic computations to the case of the Cambridge/Aachen (C/A) as well as to the Cambridge/Aachen algorithm supplemented by a filter (C/A+filt). While the former usually leads to quality measurements close to what is obtained in the anti-k t case, the latter produces narrower peaks with an optimal radius slightly larger than the one obtained from C/A without filtering. For simplicity, we shall follow what has been done in [4] and use R filt = R/2 and n filt = 2 for the parameters of the filter, i.e., recluster each jet with C/A and a radius of half the original one and keep the two hardest subjets, discarding the softer ones.
The analytic computation goes exactly as in the case of the anti-k t algorithm. The only thing that will change in our approach is the inclusion of the area contribution in the UE distribution. The expressions for the area are computed analytically, as in [22], up to  order α s :ā where the coefficient A 0 , d, Σ 0 and s are computed in perturbation theory and have been calculated 12 in [22], while the scales Q 0 are of non-perturbative origin and have been determined in [22] from a fit to Herwig simulations. In our case, we shall use p t = M/2 together with the following values for the coefficients 13 :  Figure 9, the final result of the extracted optimal radius as a function of the mass for the three algorithms under consideration, compared to Pythia simulations. We see that, in agreement with the Monte Carlo, the differences between anti-k t and Cambridge/Aachen are very small and thus our analytic approach reproduces the expected behaviour correctly. For the case of the C/A algorithm with filtering, we obtain a good description in the case of quark dijet systems but systematically underestimate the optimal radius, by about 0.1-0.2, in the case of gluon jets. We believe that this is due to the fact that we should also introduce some dependence on the algorithm at the level of the perturbative spectrum. In our soft-gluon approximation all the algorithm will have the same initial and final-state radiation spectra. However, with a more complete treatment, e.g. including non-global logarithms as in [19], differences will start appearing. Another effect comes from the combination of the initial and final-state radiation spectra: the simple convolution we have used is reasonable in the case of the anti-k t algorithm where the shape of the jet is not affected by the radiation, but more subtle phenomena will appear for more complex algorithms, especially in the filtered case where, in some geometric regions, the capture of ISR gluons can be affected by internal gluon radiation in the jet.
All these effects are beyond the reach of the present analysis. As expected from the previous discussion, discrepancies are more pronounced in the case of gluon jets at large mass, where perturbative radiation dominated the spectrum. Our computation nevertheless gives reasonable results at the end of the day, given also that for gluon dijets at larger mass, the minimum in the quality measure is less pronounced than in other situations.
Finally, let us mention that the case of the k t algorithm [21] is also very close to what we obtain for the anti-k t and C/A algorithm, hence well reproduced by our calculation. For the SISCone algorithm, differences in the perturbative spectrum also appear when the radiated gluon is not soft, we therefore postpone that case for a further study.

Comparison with earlier approaches
To some extent, the present calculation can be seen as an extended version of what has been done in [13], where the authors compute the average p t shift of a jet due to final-state radiation, hadronisation and the UE. The approach we pursue in this paper is however much more complete in the sense that we obtain the full spectrum instead of just the average value 14 . As a consequence, we can extract the width of the peak really as the dispersion in the spectrum, while in [13] the width is approximated by i.e. the total dispersion is approximated by the sum of each independent average squared. The individual contributions are found to be [13] δp t pert = α s π log(R)L a p t , gluon Pythia Figure 10: Best radius as a function of the mass for our approach compared to the previous work of [13]. The solid (blue) line corresponds to our present computation in this paper while the dotted (red) line shows the result obtained in Ref. [13]. The Pythia results, the (green) triangles, are shown for comparison. The left (right) panel shows the case of quark (gluon) jets.
with 2C F A(µ I ) = 0.5 GeV for µ I = 2 GeV and The optimal radius is then obtained by minimising (40) w.r.t. R.
In practice, we have used p t = M/2, N f = 5, Λ QCD = 200 MeV and computed the QCD coupling at the scale Rp t . Instead of using a fixed value for Λ UE as in [13], we have used the fact that it is related to the average background density, ρ , through (see Section 4.1) and used the values extracted from the Monte Carlo event samples and already used in Section 4.1 for ρ . The extracted value of R best is shown on Fig. 10 in the case of Pythia simulations, for the anti-k t algorithm. The results are compared with the Monte Carlo results and with our complete analytic calculation. We see that the simpler approach of [13] manages to extract the correct value for R best at small mass but fails to reproduce the observed behaviour at larger mass, and this for both quark and gluon jets. The reason for this failure is rather hard to pinpoint: it can come, e.g., from the crudeness of the approximation (40), the lack of an initial-state radiation contribution and the corresponding PDF effects, or the description of the underlying event in terms of the average density ρ instead of its fluctuations σ ρ .

Description
The last point we wish to discuss is the situation where we perform background subtraction using jet areas [22], as advocated in [25]. The idea behind UE subtraction is to get rid of the shift due to the UE density contaminating the mass reconstruction. Since subtraction is performed on each single event, one hopes to get rid of the fluctuations of the background density across the events, hence obtaining a narrower peak. In this respect, all we should be left with is the background fluctuations inside an event which have a linear dependence on R rather than the quadratic one computed in the unsubtracted case (32). This is a much smaller effect and one thus expect the quality measure to be better (i.e. smaller) than in the unsubtracted case as well as the optimal radius to be larger.
In practice, one would thus naively convolute the perturbative spectrum computed in Section 3 with a subtracted UE spectrum taking into account these intra-event fluctuations. Doing that would actually lead to a quality much better than what is observed in Monte-Carlo studies and to a systematic over-estimation of R best .
The additional effect that is mandatory to take into account is the fact that subtraction is performed using an estimated value for the background density ρ est which can differ from the exact value of ρ for that event. This results in some smearing effect left over from the unsubtracted spectrum. Because of this, the computation of the subtracted UE spectrum turns out to be more involved than the corresponding unsubtracted case computed in Section 4.1, mostly because one also would have to consider the distribution probability for ρ est .
The properties of the estimated value of ρ have actually recently been discussed at length in [23], where we learn that the estimation of the background density using a medianbased approach suffers from two sources of bias: one of soft origin, related to the fact that the median of a set of pure-UE jets may differ from its average; and a contamination from radiation coming from the hard jets. If we call ρ and σ the average UE density and its fluctuations for one given event, and define δ = (ρ est − ρ)/ρ as the relative error in the estimation of ρ, the average and dispersion of δ are [23] where, as in Section 4.1 we used µ = σ/ρ, and the coefficients are 15 and with In the previous set of equations, c J ≃ 2.04 is a pure number extracted in [23], R ρ is the jet radius used to estimate the background properties, A tot is the total area of the range we are using to estimate the background, Q 0 is a soft cut-off that we fixed to 1 GeV, and Q is the hard scale that we naturally took equal to M/2. These properties constrain the distribution of δ, that we shall denote by dP δ /dδ, as a function of µ. Since the intra-event fluctuations will play a substantial role in this UEsubtracted situation, we shall also consider a distribution for them. For simplicity, we shall parametrise the distribution in terms of µ rather than in terms of σ as background estimations are naturally expressed in that variable. We therefore introduce dP µ /dµ, of average µ and dispersion σ 2 µ . In the course of the computation below, we shall encounter higher cumulants for this distribution. To keep things simple, we shall work with the Gaussian assumption 16 i.e. µ 3 = µ ( µ 2 + 3σ 2 µ ), µ 4 = µ 4 + 6 µ 2 σ 2 µ + 3σ 4 µ . The values of µ and σ 2 µ can again be extracted from the Monte-Carlo samples. With this at hand, we proceed, as in Section 4.1, by computing the mass shift for a given configuration and deducing its average and dispersion. For given values of ρ, µ, ρ est , jet areas a k and jet UE contaminations p t,k , the offset in the reconstructed mass after subtraction is where the last term is the subtracted piece. After a slightly lengthy computation we reach 17 Once we have these expressions, the last step is to parametrise the subtracted UE spectrum. Since subtraction can produce negative values for δm, we cannot use (33).
The technique we shall adopt is to keep the same distribution using the unsubtracted average together with the dispersion we have just computed, and then to shift the complete distribution in order to reproduce the subtracted average δm. In practice, this gives where the parameters δm 0 , δm 1 and α are given by with δm and σ 2 δm given by eqs. (50) and (51), and δm unsub by eq. (31). Finally, it is interesting to notice that in the case of a perfect estimation of the background density ρ est = ρ, i.e. d 1 = d 2 = S d = 0, one recovers i.e. no average shift and a (squared) dispersion that is proportional to the area (linearly) and to σ 2 .

Comparison with Monte-Carlo simulations
As for the previous cases we have analysed, we shall perform a comparison of our analytic predictions with Monte Carlo simulations for the UE-subtracted case. The method goes exactly as before, so we just highlight the main steps and present our results for the optimal radius R best and stay with the anti-k t algorithm.
For the analytic part of the computation, the perturbative QCD spectrum is convoluted with the subtracted UE spectrum calculated in the previous Section. Concerning the parameters of the latter, all background properties, ρ , σ 2 ρ , µ and σ 2 µ are extracted directly from the Monte-Carlo data sample and the estimation of the background density was performed using the Cambridge/Aachen algorithm with a radius of 0.6. We have kept all the jets within one unit of rapidity around the hardest dijet, excluding the two hardest jets in the event, giving a total area A tot = 4π. It is interesting to notice that when background subtraction is performed, the final results will slightly depend on the details of how ρ is estimated, through the coefficients d 1 , d 2 and S d in our analytic approach. Extending the range, i.e. increasing its area A tot , would decrease the bias in the estimation of ρ, but in that case, there would be a risk that the background would not be uniform over the range (e.g. vary with rapidity). Following [23], we could therefore also try to optimise the choice of range and jet definition used for estimating ρ.
The computed value of R best is shown on Fig. 11

Conclusions and perspectives
Throughout this paper we have investigated analytically the optimisation of the radius R one uses with a jet algorithm to reconstruct dijet resonances. We have already learnt [4] that carefully choosing the jet radius can lead to significant improvement of the reconstructed signal and we find it valuable to supplement previous Monte-Carlo studies with an analytic understanding of these dijet reconstructions. Our approach has been to calculate the reconstructed dijet mass spectrum. This was done by convoluting fundamental pieces: the perturbative QCD radiation in the initial and final state, and the contamination due to the Underlying Event. For the perturbative part, we have worked in the soft-gluon emission approximation, a choice motivated by the fact that we are mostly interested in the behaviour in the vicinity of the mass peak. In the case of the initial-state radiation spectrum, we have also included PDF effects that can have a non-negligible impact, especially in the case of gluon jets and at large dijet mass. For the UE description, we have modelled it in such a way as to reproduce the observed properties of the background. From the reconstructed mass spectrum, we can compute a quality measure and extract an optimal radius, following the method outlined in [4].
Our analytic results have been considered and compared with Monte-Carlo simulations in various cases. We have started with parton-level events, where we can test our analytic computation of the perturbative spectrum; then we have focused on full hadronic events including the Underlying Event; and, finally, we have considered the situation with a jetarea-based background subtraction [25]. We have mostly considered the case of the anti-k t algorithm but have also considered other options.
Generally speaking, we can say that we obtain a good description of the Monte-Carlo results, especially when including the PDF effects in initial-state radiation. The optimal radius R best computed analytically is in most cases less than 0.1 away from the Monte-Carlo expectation. The comparison at a more basic level, the quality measure or even the mass spectrum, also shows a good agreement.
Most of the deviations between our analytic results and Monte Carlo simulations can probably be traced back to what we observe at the partonic level. Though our analytic computation reproduces nicely the Monte-Carlo simulations, we have observed some deviations in Section 3.3. It is interesting to notice that the differences between our description and the Monte-Carlo results is usually of the same order as the differences obtained between Pythia and Herwig. These differences propagate to the situation where the UE is included.
The one case where our analytic results deviate a bit more from the Monte-Carlo simulations is the case of the Cambridge/Aachen algorithm supplemented by a filter (see Section 4.3), for events including the UE, while our description remains good for the anti-k t and the Cambridge/Aachen algorithm without any filter. We think that a better description goes beyond the reach of the present analysis, e.g. already at the partonic level, a simple convolution between the initial and final-state spectra may be too simplistic when filtering is applied.
Let us also notice that our description of the UE allows us to treat successfully both the case where the UE is subtracted and the case where it is not. Our approach includes all the effects that contribute to the dispersion of the reconstructed mass: mostly background fluctuations from one event to another, intra-event fluctuations and area fluctuations. It is interesting to notice that in the case of subtracted UE, it was also mandatory to take into account the fact that the estimated background can differ from its "true" value. This leads to a significant additional contribution to the dispersion. Using the analytic computations from [23] led to good results for the extracted R best .
If one wishes to improve the present description, additional effects have to be included in the analytic computation: relaxing the soft-gluon-emission approximation, obtaining sub-leading logarithmic corrections, performing a better resummation of the soft virtual corrections than a simple Sudakov factor, e.g. using non-global logarithms, and including hadronisation corrections. However, we have seen that the description we can achieve with our simple approach is already satisfactory without having to introduce these corrections. It is nevertheless important to remember their existence as they might become relevant in more specific situations, like the case of filtering where we have observed that we start to reach the limits of our simple approach. Some effects, like subleading logarithmic corrections, may even be simpler to compute analytically than to implement in a Monte Carlo.
Another point that we have not investigated and that might be interesting in a future study is the case of the rapidity dependence of R best . As we have seen in Section 3.2, PDF effects introduce a dependence on the rapidity at which the heavy object is produced [6] and it might be interesting to see how R best varies with that rapidity.
Also, all the computations we have performed were for a complete event-set. A next major step would be to optimise the radius on an event-by-event basis. This is much more complicated as the previous approach would then only give access to an average dispersion for a single event and the additional dispersion arising when combining the whole set of events would be much harder to handle. We therefore also leave that step to future studies.
The question of knowing which jet algorithm should be preferred has already triggered many discussions both on the theoretical and on the experimental side, and more complex jet definitions introducing more parameters have often been proposed to improve jet reconstruction. The present work, together with [19], can be seen as opening a new direction of improvement: analytically optimising the parameters in the jet definitions. Considering a single parameter, the jet radius, for a simple process, dijet reconstructions, is a first step in that direction. Gaining analytic control on the parameters entering jet definitions can be extremely valuable. This is especially true when new parameters are introduced: being able to predict them can avoid often costly determinations based on Monte-Carlo simulations or rough empirical choices. The application of the techniques developed in this paper to other parameters and more general processes would therefore be very interesting.