Phenomenology of event shapes at hadron colliders

We present results for matched distributions of a range of dijet event shapes at hadron colliders, combining next-to-leading logarithmic (NLL) accuracy in the resummation exponent, next-to-next-to leading logarithmic (NNLL) accuracy in its expansion and next-to-leading order (NLO) accuracy in a pure alpha_s expansion. This is the first time that such a matching has been carried out for hadronic final-state observables at hadron colliders. We compare our results to Monte Carlo predictions, with and without matching to multi-parton tree-level fixed-order calculations. These studies suggest that hadron-collider event shapes have significant scope for constraining both perturbative and non-perturbative aspects of hadron-collider QCD. The differences between various calculational methods also highlight the limits of relying on simultaneous variations of renormalisation and factorisation scale in making reliable estimates of uncertainties in QCD predictions. We also discuss the sensitivity of event shapes to the topology of multi-jet events, which are expected to appear in many New Physics scenarios.


Introduction
Event shapes measure the geometrical properties of the energy flow in QCD events and, notably, its deviation from that expected based on pure lowest order partonic predictions. Event shapes, as well as being among the first observables proposed to test QCD [1], have been inextricably tied with the progress of QCD. They have played a crucial role in the extraction of the strong coupling from properties of the final-state [2]. They have been essential in tuning the parton showers and non-perturbative components of Monte Carlo event generators [3,4,5,6] and have also provided a laboratory for developing and testing analytical insight into the hadronisation process (e.g. refs. [7] and the reviews [8,9]). From a technical point of view, the development of resummations and fixed-order calculations has benefited from comparisons of predictions for event-shape distributions obtained with both kinds of methods [10,11,12,13,14,15,16]. Additionally, they are one of the several tools that are used for classifying hadronic final states in new physics searches.
The majority of investigations of event shapes has been performed for e + e − colliders, with significant work also in DIS. A review of some of that work is given in [17]. In contrast, few dedicated studies have been performed on them at hadron colliders, with a handful of measurements at the Tevatron [18,19,20], a pure fixed order study [21], pure resummations in [22,23,24], and a recent experimental simulation study by CMS [25], as well as some investigation of the use of event shapes applied to jet contents for the identification of hadronic decays of boosted massive particles [26,27,28] (other approaches are reviewed e.g. in ref. [29]). The purpose of this article is to help bring our understanding of hadroncollider event-shape phenomenology closer to the level of sophistication that is standard in the e + e − and DIS cases, concentrating specifically on event shapes in hard QCD (dijet) events.
As is well known from the e + e − and DIS cases, accurate studies of event shapes involve the simultaneous use of two kinds of calculation. Fixed-order calculations provide expansions of event-shape distributions in powers of the strong coupling, α s . They are available up to next-leading-order (NLO) for hadron-collider event shapes, through the nlojet++ [21] program. When the event shape has a value v ≪ 1, for each power of α s in the distribution there can be up two powers of a large logarithm, L = ln 1/v, associated with soft and collinear enhancements. This compromises the convergence of the perturbative series. The enhanced terms can however be resummed to all orders, providing the dominant contribution to the distribution for v ≪ 1. Such resummed predictions tend to be carried out for the distribution integrated up to some value v, which generally has an exponentiated structure exp(Lg 1 (α s L) + g 2 (α s L) + α s g 3 (α s L) + . . .). 1 The Lg 1 (α s L) term gives leading logarithmic (LL) accuracy in the exponent, g 2 (α s L) is next-to-leading-logarithmic (NLL), etc. For suitable (recursively infrared and collinear safe, global) observables, the caesar program [31] calculates both Lg 1 (α s L) and g 2 (α s L). To obtain a reliable prediction for an event shape distribution, one must combine both types of calculations, via a "matching" procedure [32]. An appropriately performed matching of NLO fixed order and NLL exponentiated resummation allows one to ensure that in the expansion of the resummation one correctly accounts for all terms α n s L p , with 2n − 2 ≤ p ≤ 2n, which is NNLL in the expansion. While NLO+NLL with NNLL in the expansion is the state-of-the-art for generic e + e − and DIS event shapes, 2 matching of this kind had not so far been achieved for hadron-collider event shapes.
Our study here will give NLL+NLO predictions for event shape distributions both at Tevatron (pp, √ s = 1.96 TeV) and at the LHC (pp, √ s = 14 TeV). We start in Sec. 2 by recalling the definitions of three classes of global event shapes for hadron colliders, 3 addressing also the question of event shapes defined in terms of jets rather than particles. In Sec. 3 we describe the structure of the perturbative resummation as well as its matching to the NLO result. We also discuss possible general event-shape resummation issues associated with "super-leading" logarithms [39,40]. In Sec. 4 we present our results for matched distributions, paying particular attention to the issue of theoretical uncertainties. We also compare our results to those obtained with parton-level shower Monte Carlo event generators, in some cases also matched to exact multi-parton tree-level matrix elements. In Sec. 5 we briefly discuss the impact of non-perturbative effects, the hadronisation and the underlying event. Finally, switching to more phenomenological questions, in Sec. 6 we compare various event shapes' ability to distinguish characteristically different event topologies, and examine their robustness in such tasks, both with respect to parton showering and to the orientation of the final state event. Our results are summarised in Sec. 7. Some technical details are collected in Appendices A and B. Many further additional plots can be obtained from the URL [41].

Event-shape definitions at hadron colliders
In this article, we shall consider observables that measure the extent to which an event's energy flow departs from a dijet structure. The lowest-order contribution to a dijet event consists of just two incoming and two outgoing partons. Throughout the paper we refer to these QCD configurations as the "Born limit". Many of the event shapes studied here were presented for the first time in [23]. All share the property of continuous globalness [42,31], which is a necessary condition for being able to carry out a resummation to NLL accuracy without a leading-N C approximation, and which also contributes to the simplicity of caesar's generalised resummation approach (independently of the question of large-N C approximations). For an observable to be continuously global, it has to be sensitive to all emissions in an event (this is the requirement of globalness), and moreover it should have definite scaling properties with respect to secondary emission's transverse momenta (see sec. 3.1.1 for a mathematical formulation). The continuously global event shapes we propose fall into three main classes: observables that are directly global, others that are supplemented with "exponentially suppressed forward terms" and observables with "recoil terms".

Directly global observables
We first consider observables that are defined in terms of all hadrons in the event, therefore the name 'directly' global. The global transverse thrust is defined as where the sum runs over all particles q i in the final state, q ⊥i represents the two momentum components transverse to the beam, q ⊥i its modulus, and n T is the transverse vector that maximises the sum. The observable which is resummed is then τ ⊥,g ≡ 1 − T ⊥,g , which vanishes in the Born limit. The normalization of event shape observables to a hard transverse scale of the event is important because it reduces uncertainties associated with the experimental jet-energy scale, which partially cancel between numerator and denominator [43]. For most event shapes (except τ ⊥,g ) the choice of specific hard scale to which one normalises is arbitrary, and could for example also be the sum of the transverse momenta of the two hardest jets. The transverse thrust axis n T and the beam form the so-called event plane. One can then define a directly global thrust minor, which is a measure of the out-of-event-plane energy flow In close analogy with the e + e − case [44], one can formulate a transverse spherocity: where the minimisation is carried over all possible unit transverse 2-vectors n. 4 This variable ranges from 0 for pencil-like events, to a maximum of 1 for circularly symmetric events. An alternative observable, which makes use of a linearised version of the transverse momentum tensor (with direct analogy to the C and D parameters [45] used in e + e − ), is the F -parameter: where λ 1 ≥ λ 2 are the two eigenvalues of M lin . Related variables have been considered in the plane transverse to the thrust axis in e + e − (resummed for 3-jet events in [46]) and in the plane transverse to a jet in the context of boosted top-quark decays [26,28,27], where forms involving the determinant of M lin , e.g. 4λ 1 λ 2 /(λ 1 + λ 2 ) 2 = 4F/(1 + F ) 2 , have been used. There is a one-to-one mapping between different forms, and we have chosen eq. (2.4) because it gives clearer separation between different kinematic regions. Finally, we consider the exclusive variant of the k t -algorithm [47] (closely related to the inclusive variant [48] as adopted for Run II of the Tevatron [49] and expected to be used also at the LHC) 1. One defines, for all n final-state (pseudo)particles still in the event, and for each pair of final state particles where y i = 1 2 ln E i +p zi E i −p zi is the rapidity of particle i and φ i its azimuthal angle. The jet-radius parameter R sets the angular reach of the jet algorithm. Throughout this paper, we will take R = 0.7.
2. One determines the minimum over k and l of the d kl and the d kB and calls it d (n) .
If the smallest value is d iB then particle q i is included in the beam and eliminated from the final state particles. If the smallest value is d ij then particles q i and q j are recombined into a pseudoparticle (jet). A number of recombination procedures exist. We adopt the E-scheme, in which the particle four-momenta are simply added together, q ij = q i + q j . (2.7) 3. The procedure is repeated until only 3 pseudoparticles are left in the final state.
The observable we resum is the directly global three-jet resolution parameter where P ⊥ is defined by further clustering the event until only two jets remain and taking P ⊥ as the sum of the two jet transverse momenta, While directly global event-shapes are defined in terms of all particles in the event, experimental measurements can be carried out only up to some given pseudorapidity η max (η max ∼ 3.5 at the Tevatron and η max ∼ 5 at the LHC). However, as long as the eventshape's value v is not too small [22], v > v min , one can safely neglect the contribution of hadrons beyond the rapidity cut. The value v min up to which a NLL resummation is valid is observable specific. In particular, it depends on the behaviour of each event shape under a soft and collinear emission, as derived in [22,23]. Further discussion is given in appendix B.

Observables with exponentially suppressed forward terms
One way to address the difficulty in performing measurements near the beam is to define event-shapes using only particles in a central region and to add a term sensitive to emissions along the beam direction, so as to render them global, but with an exponential suppression in the forward or backward directions. We define the central region C by requiring that the rapidity of particles in C satisfies |η i | < η c = y j,max + δη, where y j,max specifies the rapidity region in which the two highest p t jets should lie, and δη is a rapidity buffer around the jets of size ∼ 1.
Given the central region C, we introduce the mean transverse-energy weighted rapidity η C of this region, and define the exponentially suppressed (boost-invariant) forward term as We can then define non-global variants of the event-shapes defined in sec. 2.1 by restricting the sums to just the central region. For example, we have a central transverse thrust, and a central three-jet resolution threshold, y 23,C defined by the algorithm of section 2.1 applied only to the final state particles in C (but maintaining the "beam" distance, eq. (2.5)). Finally, we define "exponentially suppressed" variants of the event-shapes of sec. 2.1 by adding to the central event-shapes a power of EC which makes the event-shape continuously global [23]. We obtain the exponentially suppressed transverse thrust, thrust minor and three-jet resolution, 14) T m,E = T m,C + EC , Additionally, one can consider event-shapes which are more naturally defined using particles only in a restricted region, like jet-masses and broadenings. Given a central transverse thrust axis n T,C , one can separate the central region C into an up part C U consisting of all particles in C with p ⊥ · n T,C > 0 and a down part C D with p ⊥ · n T,C < 0 respectively. One then defines, in analogy with e + e − [50], the normalised squared invariant masses of the two regions from which one can obtain a (non-global) central sum of masses and heavy-mass, 18) and the corresponding global event-shapes, the exponentially-suppressed sum of masses and heavy-mass With the same division into up and down regions as for the jet masses, one can define jet broadenings. To do so in a boost-invariant manner, one first introduces rapidities and azimuthal angles of axes for the up and down regions, 20) and defines broadenings for the two regions, from which one can obtain central total and wide-jet broadenings, Adding the forward term one obtains the global exponentially-suppressed total and wide-jet broadenings, We note that an observable that effectively has exponentially suppressed forward behaviour has also been studied in [24].

Observables with recoil term (indirectly global observables)
Because of transverse momentum conservation, if radiation is emitted in the forward region C, recoil effects will cause the vector sum of the transverse momenta in the complementary, central region C to be non-vanishing. It is then possible to exploit this effect to make observables (continuously) global, despite the fact that only a central subset of particles in the event effectively enters the definition of the event shapes. To do so, we add to the central event-shapes a suitable power of a recoil term, the two-dimensional vector sum of the transverse momenta in C, We obtain than recoil enhanced transverse thrust, thrust minor, three-jet resolution, sum-and heavy-jet masses, total and wide broadenings

Particles versus jets as inputs
The event shapes discussed so far have all been defined in terms of the particles in the event. The experiments don't measure particles directly. They do, however, have methods such as the combination of information from electromagnetic and hadronic calorimeters into "Topoclusters" (ATLAS [51]) and, with tracking, "particle flow" (CMS [52]), that provide inputs to jet algorithms that are quite close to particles. These same inputs would probably also be well suited to event-shape studies.
In uses of event shapes to cut on event topology in beyond-standard-model searches, as well as in the study of ref. [25], it is not particles but instead jets that have been used as inputs. The jets are usually defined through an angular resolution parameter R (as in eq. (2.6)) and a transverse momentum cutoff, which we will denote p t0 . One of the interests of using jets is that the p t0 cutoff eliminates much of the contamination from the underlying event, which can easily contribute O (100 GeV) of transverse momentum to the rapidity region covered by LHC detectors.
From the point of view of resummation, the use of jets as inputs poses two main problems. One comes from the presence of the new scale p t0 in the problem: in terms of the parameters a and b 1,2 defining the event-shape's sensitivity to radiation along the incoming legs (cf. table 1 and section 3.1.1), this new scale causes separate regions of event shape value to each involve different logarithmic structure: for cases with b 1,2 > 0 the potentially different regions are v ≫ (p t0 /Q) a , (p t0 /Q) a ≫ v ≫ (p t0 /Q) a+b 1,2 and v ≪ (p t0 /Q) a+b 1,2 . The first of these regions may be within the scope of caesar if p t0 /Q is sufficiently small.
A second problem is that of globalness. Emissions collinear to any outgoing hard parton will be clustered together with its emitter to form a jet. The observable's sensitivity to these emissions will then depend on the jet recombination scheme. In the E-scheme, the current default at the Tevatron and LHC, the jet four-momentum is constructed by simply adding the four momenta of its constituents. Therefore, all observables defined using transverse momenta will get no sensitivity to emissions inside each of the two hard jets, and will therefore be non-global. This statement is true for any recombination scheme that adds three-momenta vectorially. For variables with sensitivity to longitudinal degrees of freedom, globalness can only be assessed on a case by case basis. For instance if one considers any global version of the total and heavy-jet mass (with exponentially suppressed or recoil term), in the E-scheme the mass of each central hard jet will enter the hemisphere central jet mass ρ X,C in eq. (2.17). Therefore one obtains the same result for the central component of the event shape as would have been obtained using hadrons as inputs (modulo the fact that the jet clustering may affect which particles are considered central).
One alternative to the use of jets as inputs, in order to avoid the globalness issue, is the following: use as inputs the particles that are inside the two hardest jets, together with all the remaining jet momenta. Note that this does not eliminate the issue of the extra scales related to p t0 , though it does maintain the reduced sensitivity to underlying event that comes from the use of jets.

Structure of the perturbative calculation
Typically one wishes to consider event shapes only for events that are sufficiently hard, requiring for example at least one jet above some minimum transverse momentum threshold p t,min and in some central rapidity region. We will denote this kind of hardness selection cut by a function H(q 1 , . . . , q N ) of the N particles in the event; H(q 1 , . . . , q N ) is equal to 1 for events that pass the cuts and 0 otherwise. One can then define the cross section for events that pass the cuts, where dσ N /dΦ N is the differential cross section for producing N particles in some configuration Φ N . One can determine σ perturbatively as long as H corresponds to an infrared and collinear (IRC) safe selection procedure. One also defines the partial integrated cross section Σ(v) for events that pass the cut and for which additionally the event shape observable V (q 1 , . . . , q N ) is smaller than some value v, The differential normalised distribution for the event shape is then given by Perturbatively we will write σ and Σ(v) as expansions in the number of powers of the coupling that they contain, where σ 0 is the leading order (LO) result, σ 1 the NLO result, etc.; σ i is proportional to α 2+i s . We have chosen not to extract the powers of α s from the σ i coefficients, because the scale of α s may depend on the kinematics of the events over which one has integrated.
The expansion for Σ(v) is similar with the property that Σ 0 (v) ≡ σ 0 because the observable vanishes at O (α 2 s ). Σ 1 (v) looks like a NLO term in eq. (3.5), but it is usually determined from the LO α 3 s term for the differential cross section of v, The quantityΣ 2 (v) is similarly determined from the NLO term of the differential cross section of v. In the following we shall never use explicitly Σ 2 , since σ 2 , the NNLO correction to the dijet cross section, has yet to be calculated and since its effect would lead to terms that are beyond our accuracy in differential distributions.

Resummation
Resummations are relevant in the region of small v, where logarithmically enhanced contributions of soft and collinear origin, as large as (α s ln 2 v) n , appear at all orders in the integrated cross section Σ(v), thus making fixed-order predictions unreliable. There is a large class of observables for which one can write a common "master" resummation formula, as was done in [31], in order to sum such terms to all orders in α s . In this section we will first examine what the class of observables is, and then review the broad structure of the resummation.

Prerequisites for resummation with caesar
In order for an observable to be resummed within the caesar framework, its functional behaviour in the presence of an arbitrary number of soft and/or collinear emissions has to satisfy a number of conditions. These have been extensively discussed in [31], and are checked automatically by caesar given a computer subroutine that computes the value of the an observable given a set of four-momenta. The conditions are: 1. a specific functional form for the observable's dependence V ({p}, k) on the momentum of a single soft emission k, collinear to one of the hard "Born" partons ("legs") in the event: where {p} denote the Born momenta (including recoil effects) and k is the soft collinear emission; k (ℓ) t and η (ℓ) denote respectively its transverse momentum and rapidity, as measured with respect to the Born parton ('leg') labelled ℓ; φ is the azimuthal angle of the emission with respect to a suitably defined event plane (when relevant); g(φ) can be any function for which dφ ln g(φ) is well defined; Q is a hard scale of the problem (taken here to be the sum of the transverse momenta of the two hardest jets).
2. continuous globalness [11,53], a requirement on the observable's single-emission scaling properties in every region of the phase space. First, all the a ℓ have to be equal, a 1 = a 2 = . . . ≡ a, and the d ℓ have to be all non-zero. Second, the observable's scaling at the boundaries of the soft collinear region has to be consistent with eq. (3.7), i.e. in the soft large-angle region we require V ({p}, k) ∼ k a t for a fixed angle of k, whilst for hard emission collinear to leg ℓ we must have V ({p}, k) ∼ k a+b ℓ t at fixed energy for k.
3. recursive infrared and collinear safety, a subtle mathematical condition (see [31] for its precise formulation) concerning the observable's scaling in the presence of multiple soft/collinear emissions. Table 1 summarises the values of the coefficients a ℓ and b ℓ for the event shapes presented in section 2. We stress that central observables, like the central transverse thrust eq. (2.12), defined using only hadron momenta in a selected rapidity interval, tend to have d 1,2 = 0, and therefore be non-global. [but not true for τ ⊥,g .] The exponentially-suppressed term EC in eq. (2.11) or the recoil term R ⊥,C in eq. (2.24), as explained in sections 2.2 and 2.3, are added to central event shapes precisely so as to make them global. The different powers of EC and R ⊥,C that appear in the definition of these modified event shapes (see for instance eqs. (2.14) and (2.16)) are chosen so as to ensure their continuous globalness. This can be seen by observing that, for each event shape, the coefficients a ℓ corresponding to different legs are equal. This above discussion holds for most observables but there may be exceptions. For example, the central variant of the thrust-minor T m,C is actually a global observable because of an indirect sensitivity to non-central emissions due to recoil. This is the reason why the b 1,2 coefficients for T m,E are not those that usually appear for observables with exponentially-suppressed components, but are rather those typical of a (linear) recoil term. Recursive infrared and collinear (rIRC) safety, a detailed discussion of which is beyond the scope of the present paper, is trivially satisfied for all observables that we discuss here. Table 1: Table of event shapes being considered here and the powers of their parametric sensitivity to the transverse momentum (a ℓ ) and collinear angle (b ℓ ) of an emission along incoming (a 1,2 , b 1,2 ) and outgoing (a 3,4 , b 3,4 ) hard partons.

NLL resummation structure
For global observables, in events with v ≪ 1, it is possible, unambiguously, to associate the event kinematics with that of a 2 → 2 (Born) event. This is because the requirement v ≪ 1 forces all radiation to be either soft or collinear. At perturbative level it is also possible to unambiguously attribute a partonic subprocess to the event, for example qq → qq (doing so requires a flavour infrared and collinear safe procedure, as in [54], but the result is independent of the choice of the procedure). Here we will use B to label the event's 2 → 2 kinematics and δ to label its 2 → 2 flavour structure. Then we can write Σ(v) as a sum over partonic subprocesses and an integral over Born configurations that pass the hard event cuts, The ambiguities in such a decomposition of Σ(v) are suppressed by powers of v.
For observables that satisfy the properties of the previous section, the result of ref. [31] is that we can write where dσ (δ) 0 (v)/dB is the LO cross section, differential in the Born configuration, separated into subprocesses, and understood to have been evaluated with a factorisation scale µ F ∼
The order-by-order expansion of f (δ) B (v) involves terms of the form α n s L 2n . It is because of the property of "exponentiation" (a consequence of rIRC safety and of coherence 5 ) that one can write it in the form eq. (3.10), whose leading-logarithmic (LL) contribution in the exponent, Lg 2,B (α s L, µ R , µ F ) resums "next-to-leading logarithmic" (NLL) terms in the exponent, α n s L n , also referred to sometimes as singlelogarithmic terms.
The LL function Lg ℓ is the colour charge (C F or C A ) of hard parton ℓ for the hard-scattering subprocess δ, λ = α s β 0 L and β 0 = (11C A −4T R n f )/(12π). Since the coefficients a ≡ a 1 = a 2 = . . . and b ℓ do not depend on the particular momentum configuration B of the hard partons, g (δ) 1 (α s L) is also independent of B. Its dependence on the subprocess arises only through the colour charges of the incoming and outgoing partons.
The NLL function g (δ) 2,B (α s L, µ R , µ F ) can be decomposed into three types of terms [31] g (δ) The term g 2s (α s L) accounts for NLL corrections associated with the event kinematics, the particular values of the d ℓ and g ℓ (φ) coefficients in eq. (3.7), the choice of renormalisation scale and scheme used in α s in g 1 (α s L), as well as the non-trivial colour evolution of large-angle soft virtual gluon resummation [56,57,58,59,60]. The second term on the right-hand side (RHS) of eq. (3.12) involves parton distribution functions (PDF) for the parton flavours in the initial state of the given subprocess, q for each leg that depends on the Born kinematics. This term arises because the PDFs in dσ (δ) 0 /dB in eq. (3.9) were evaluated at a factorisation scale µ F ∼ Q. The presence of a PDF at scale µ F ∼ Q implies that one integrates over all possible incoming collinear emissions, up to k t ∼ Q. However the requirement that the event shape be small, V (k) v, restricts collinear emissions to have k t v 1/(a+b ℓ ) Q. Thus the PDFs should actually be evaluated at a factorisation scale ∼ v 1/(a+b ℓ ) Q ∼ v 1/(a+b ℓ ) µ F (as occurs also in Drell-Yan transverse-momentum resummations [55,61,62]). The ratio of PDFs in eq. (3.12) serves to replace q ℓ , µ F ) as used in the Born cross section with a PDF at the correct factorisation scale.
The third term on the RHS of eq. (3.12) accounts for the NLL corrections associated with the presence of multiple soft and collinear emissions, when each has V (k) ∼ v and they are all widely separated in rapidity. It is a function of and is known analytically for some observables (e.g. τ ⊥,g ), while in all other cases caesar can compute it numerically via a suitable Monte Carlo procedure. F (δ) (R ′ ) sometimes depends on the underlying scattering channel δ, but not (for the observables studied here) on the hard momentum configuration. The behaviour of F (R ′ ) with increasing R ′ (decreasing v) is a characteristic feature of each event shape. It depends on whether multiple emissions tend to increase or decrease the value of the event shape. In the first case, for a fixed value v, F (R ′ ) has to account for an extra suppression of emissions so as to keep the event shape's value less than v, i.e. F (R ′ ) < 1. For the special case of V (k 1 , k 2 , . . .) = max(V (k 1 ), V (k 2 ), . . .) then F (R ′ ) ≡ 1 (for example the y 23 jet resolution threshold [12] in the e + e − Cambridge jet algorithm [63]). Conversely if the contributions of multiple emissions tend to cancel, the function F (R ′ ) has to compensate the excessive suppression given by the LL function Lg 1 (α s L), therefore F (R ′ ) > 1.
This last case appears most dramatically when it is a cancellation between multiple emissions, and not a direct veto on real emissions, that is the dominant effect that keeps the event shape small. In this case the LL function Lg 1 (α s L) (whose functional form depends only on the effect of single emission) no longer accounts for the dominant contribution to the distribution. Furthermore, no NLL function such as F (R ′ ) can fully compensate for this. This inconsistency reveals itself through a divergence of F (R ′ ) at a given critical value R ′ c , which can be inferred from considerations on the cancellation mechanism, as explained in refs. [23,31]. Such a divergence is present, for example, for T m,C and T m,E and occurs at R ′ c = C T /(C 1 + C 2 ), where C 1 and C 2 are the colour charges of incoming partons and C T the total colour charge of the hard parton system. It will prevent us from obtaining sensible NLL resummed results for these observables. In the case of recoil observables there is also a divergence, but at larger R ′ c (smaller v), for example at R ′ c = 2C T /(C 1 + C 2 ) for T m,R and B T /W,R . The effect on the corresponding differential distributions will be discussed when presenting matched results.

NNLL Σ accuracy
As well as discussing the LL, NLL, etc. accuracy of resummation in the exponent of eq. (3.10), one can also discuss the accuracy in the order by order expansion of Σ itself.
In this way of counting logarithms, "LL Σ " terms involve powers α n s L 2n , NLL Σ involve α n s L 2n−1 , etc. A NLL resummation in the exponent automatically guarantees NLL Σ accuracy. However it is also possible (if not entirely straightforward), given the information at our disposal, to obtain NNLL Σ accuracy. To see how, observe that the terms that we neglect in eqs. (3.9) and (3.10) are an overall α s correction without logarithms, as well as terms α n s L n−1 in the exponent, starting at α 2 s L. The latter, if they multiply α n s L 2n when expanding the exponent, lead at most to α n+2 s L 2n+1 ∼ α n s L 2n−3 , i.e. they are NNNLL Σ . We can therefore ignore them. As for the overall α s correction, when multiplied by the double logarithms, it gives us terms α n+1 s L 2n ∼ α n s L 2n−2 , which are NNLL Σ and therefore cannot be neglected. This means that we need to determine the coefficient of the pure O (α s ) term.
To do so, let us define the NLL resummed cross section as NLL,B(v) containing only the LL and NLL resummation terms, Then we can determine the coefficient C (δ) 1,B in terms of the first order expansion of the exact and resummed distributions, and in particular their difference as v → 0, The C (δ) 1,B constant involves many contributions, including parts that cancel the µ R and µ F dependence present in the Born cross section, parts that are sensitive the observable's exact behaviour with respect to soft large-angle emission and hard collinear splitting and parts related to the exact structure of the 1-loop 2 → 2 scattering diagram. Now we can write the NNLL Σ resummed distribution as The fact that we may multiply (1 + α s C (δ) 1,B ) and dΣ (δ) in order to get NNLL Σ accuracy is a consequence of the property that soft-collinear virtual corrections, which give powers of α s L 2 , affect neither the flavour, the momentum nor the colour involved in the hard scattering or the PDFs and therefore straightforwardly multiply all the more complicated contributions that are present in C gives a distribution that is still correct to NNLL Σ accuracy, because the LL, exp(Lg In contrast, if one considers C 1 averaged additionally over subprocesses is not accurate to NNLL Σ , because

Matching of NLL to NLO
While the resummation of logarithms is necessary in the region where event-shape values are small and their logarithms large, the region of large values of V is dominated by events with three or more well separated jets. Those types of events are described more reliably by fixed order calculations. It has therefore become standard to match resummed calculations to next-to-leading order (NLO) to have a reliable prediction over a larger range of values of V . In this section we will present the formulae we use to perform the matching.
In the following we will denote with f (v) = Σ(v)/σ the integrated event-shape fraction, where Σ(v) and σ are defined in eqs. (3.2) and (3.1) respectively. After a NLL+NLO matching this quantity should satisfy the following requirements 1. it should respect the physical constraints that, when the event shape reaches its maximum value v max , we have f (v max ) = 1 exactly and df (v) dv v=vmax = 0; 2. its expansion up to relative O (α 2 s ) should reproduce the exact NLO result for the corresponding differential distribution; 3. one should obtain NNLL Σ accuracy, i.e. all logarithms O (α n s L m ) with m ≥ 2n − 2 should be correctly accounted for, which implies that the matching formula should reduce to eq. (3.19) in the limit of small v. Preferably this should be the case without having to go through the tedious procedure of manually determining C (δ) 1 separately for each event shape.
There are various matching procedures that satisfy these requirements and therefore formally have the same accuracy. We consider here the so-called log-R [32] and multiplicative [42] matching schemes, adapted to hadron-hadron collisions. In particular both need to be modified in accordance with the need, section 3.
term multiply the resummation separately for each subprocess. Actually, what matters is not so much the subprocess but the colour charges of the incoming and the outgoing Born partons. Therefore we can consider all subprocesses qq → qq, qq ′ → qq ′ , qq → qq, etc., with the same incoming and the same outgoing colour charges as belonging to a single colour channel a = qq → qq. The other colour channels are qg → qg, qq → gg, gg → qq and gg → gg. 6 We denote by Σ We have obtained fixed order cross-sections using the code nlojet++ [21]. The publicly available version computes cross-sections summed over the flavour of outgoing partons. We therefore extended it so as to have access to the flavour of both incoming and outgoing partons in the calculation of σ 0 , σ 1 and Σ 1 (v), though not forΣ 2 (v) since its colourchannel separation is not needed for NNLL Σ accuracy. To assign events with more than two outgoing partons to a definite 2 → 2 colour channel we used the exclusive flavour-k t algorithm of [54] to cluster events to a 2 → 2 topology while keeping track of flavour in an infrared safe manner. During the clustering procedure, quarks of different flavour might end up in the same jet, giving rise to multi-flavoured jets, i.e. jets whose flavour does not correspond to any QCD parton. These events, which do not correspond to any Born 2 → 2 processes and have vanishing weights for v → 0, will be labelled as having a = other.
The matching equations are defined in terms of the following resummed distribution (and its fixed order expansionsΣ where v max is the maximum kinematically allowed value of the event shape, so thatL(v = v max ) = 0. We take the values of v max from the NLO calculation, which is sensible since we want the differential distributions to reproduce the NLO result at high v. The factors x V and p modify the definition of the logarithm that one is resumming. The main effect of x V is to modify the logarithm at small values of V , and it will therefore affect subleading logarithmic terms (the change at NLL is cancelled via a suitable compensatory term in g 2 (α s L)). The main effect of p on the contrary is to modify L at large values of V , it will therefore mainly affect power suppressed terms. Our default values for x V , X V are given by (see Appendix A of ref. [23]) are fixed by setting X = 1 and In the same way as one varies renormalization and factorization scales around a central value by a factor of 2, we will probe the x V dependence through a variation of X in the range 1/2 ≤ X ≤ 2. This will provide an estimate of the error associated with unknown NNLL contributions to the resummation. In principle, one can also vary the power p around the value 1 (as a probe of terms that are suppressed by powers of v), though for simplicity in the following we just fix p = 1 and therefore do not include any uncertainty related to its variation. We now introduce the log-R matching formula whereΣ 2 (v) has been introduced after eq. (3.6). It is straightforward to verify that with this matching equation f (v) satisfies all three requirements listed at the beginning of this section.
The alternative, multiplicative matching (mod-R) scheme that we use is This matching equation has the same matching accuracy as eq. (3.26), so that using both matching procedures provides an additional way of estimating the uncertainty in the matched distributions.
In both matching formulae,Σ r,1 (v) andΣ r,2 (v) require the calculation of the order α s and α 2 s expansions of the ratios of PDFs at different scales that appear in eq. (3.12). These have been obtained using hoppet [64].
In the following we will present results for normalised differential distributions Notice that two-loop corrections to σ 2 , currently unknown, are not needed for a second order matching, as they do not contribute to the differential distribution within the target accuracy.

Coherence-violating (super-leading) logarithms
One of the assumptions that enters into the derivation of the generalised resummations of [65] is "coherence" [66], the property that real emissions and virtual corrections at large angles are independent of the structure of real emissions that have occurred at small angles (with respect to any of the incoming and outgoing legs). Physically this can be understood as arising because a large-angle emission (or virtual correction) sees only the sum of colour charges of a bunch of collinear partons and that sum of colour charges is conserved under collinear splitting. 7 The assumption of coherence is challenged by the results of ref. [39,40], which found "super-leading logarithms" (SLL), terms that go as α 4 s L 5 , when calculating the probability of there being no soft radiation (above scale Qe −L ) in a finite patch of rapidity and azimuth. Based on coherence, one would have expected only terms α n s L m with m ≤ n for such an observable. Therefore, one might also call the terms of [39,40] "coherence-violating logarithms" (CVL), a name that is suitable also in the case of observables whose leadinglogarithmic structure involves double logarithmic terms α n s L 2n (for which α 4 s L 5 is not superleading).
The interpretation of the result in ref. [39,40] is that one specific class of (soft) single logarithmic virtual correction, "Coulomb-gluon exchange," can be affected by small-angle (collinear) initial-state gluon emission, independently of how small that angle is. This is because in the calculation of [39,40] Coulomb gluons are exchanged either between two incoming partons or between two outgoing partons but not between one incoming and one outgoing parton (whereas other classes of soft contribution treat incoming and outgoing partons on an equal footing); real initial-state splittings, however small in angle, lead to a redistribution in colour between incoming and outgoing states and therefore Coulomb-gluon exchange cares about them (but not about the corresponding collinear virtual initial-state corrections). This means that the coefficient of the Coulomb single logarithms α n s L n is proportional to the probability of soft-collinear initial state emission, α m s L 2m and hence one obtains terms α n+m s L n+2m , which are super-leading with respect to the expected α n s L n . The calculation of the impact of this effect requires that one follow through the soft colour evolution of the 2 → 2 scattering [56,57,58,59,60], for which the Coulombexchange terms provide imaginary contributions. For the purposes of our discussion here it is not necessary to enter into the full detail of the soft colour evolution. Rather it suffices to be aware, following [40], that for the case of vetoing emissions into a finite patch (gap) the lowest order coherence-violating terms are contributions with structures such as where we have shown only one of the two orderings given in [40] (the other gives either the same number or fewer logarithms, depending on the observable). In the integration measures, we have labelled each momentum with (v) if it can only be virtual, and (r/v) if we are considering the difference between real and virtual cases. Gluon 2, the collinear, possibly real gluon, can have an angle corresponding to anywhere outside the gap region, down to the smallest kinematically allowed angles θ 2 ∼ k t2 /Q. Gluons 1, 3 and 4 have only transverse momentum integrations because they are either Coulomb exchange gluons or the virtual counterparts of large-angle soft-gluon emission. (In the other ordering it is gluon 1 that is collinear and possibly real). The integral for gluon 4 is limited to be above Qe −L because below that scale the observable places no constraint on real emissions and so all real and virtual effects should cancel, by virtue of unitarity. Finally, the constant C depends on the kinematics of the hard scattering and the definition of the gap region. The extension to the event-shapes case involves restricting the (r/v) integration for gluon 2 to regions of phase-space that are consistent with the the real gluon's contribution to the event shape being e −L . Eq. (3.30) therefore becomes where the coefficients a and b are those that appear in eq. (3.7) for the incoming legs, ℓ = 1, 2, for simplicity we have neglected the d l and g ℓ factors there, and the constant C may differ somewhat from that in eq. (3.30) (since there it could depend on the gap definition). Three cases arise: where for the cases b ≤ 0 the number of powers of L is that obtained when relaxing the constraint (consistent with the fact that we have ignored factors of d ℓ , g ℓ ). In words, the CVL contributions only appear with a large number of logarithms when the collinear gluon, 2, if real, is allowed to be harder than the virtual Coulomb gluons (3,4). 8 For observables with b < 0 the one logarithm arises from the integration over k t1 , while the Qe −L/a k t2 θ b/a 2 constraint forces k 2 to be at large angles with k t2 ∼ Qe −L/a , with the knock-on effect that k t3 and k t4 should also be ∼ Qe −L/a . For b = 0, we instead have just Qe −L/a k t2 for gluon 2, and an extra logarithm then arises from the integration of θ 2 in the collinear region.
The expectation for yet higher order terms is that for b ≤ 0 (all observables of this paper except the exponentially-suppressed ones) the results in eq. (3.32) could be multiplied by additional powers of α s L giving at worst a series α n s L n−2 , which is subleading both with respect to our NLL accuracy in the exponent and to our NNLL Σ accuracy in its expansion. 9 For b > 0 (the exponentially suppressed observables) one would obtain terms α n s L 2n−3 . It is not clear how they would fit into the exponential resummation, except that they would certainly destroy NLL accuracy in the exponent; in the expansion of the resummation they would represent terms NNNLL Σ and therefore be subleading with respect to our accuracy.
One caveat with regard to the above discussion is that the results of [39,40] have been obtained in a strongly ordered eikonal approximation, with the assumption that the "strong ordering" parameter is transverse momentum. With different assumptions, the results change. For example, if the correct ordering parameter were energy, then eq. (3.30) would be modified in such a way as to give an infinite result. If instead one considered emissiontime (or virtuality) ordering, 10 which leads to k t1 ≫ k t2 θ 2 ≫ k t3 , then the contribution of eq. (3.30) would be halved, 33) and the corresponding result for the event shapes case would become In this case, for all the observables being discussed in this paper, the CVL terms would be subleading relative to our accuracy.
To conclude: given today's knowledge it is not clear whether coherence-violating terms matter at our accuracy for event-shape resummations. The critical issue is that of the appropriate ordering parameter (k t , time or virtuality ordering, or some other ordering). The correct ordering needs to be derived (by going beyond the eikonal approximation), unless of course there exists some yet-to-be found contribution that cancels the CVL terms. If CVL terms do exist and k t ordering is correct, then they will invalidate our statement of NLL accuracy in the exponent for the exponentially suppressed class of observables, though not our statement of NNLL Σ accuracy in the expansion of the distribution. In practice we have reason to believe that their numerical impact will still be small: partly because the CVL terms were already not very large in [39]; and partly because the large colour factors multiplying the double logarithms of our resummation force the majority of events to be in a region where the logarithms are not actually all that large (a reflection of this will appear in section 4.5, where we will see that naive exponentiation of the NLO calculation is not too different from the full NLO+NLL result, even though it misses classes of LL terms in the exponent and LL Σ terms in its expansion). particular attention to the estimation of uncertainties on our predictions, and comparisons to separate pure NLO and NLL calculations. We will also compare our results to Monte Carlo parton-shower results, with and without tree-level matrix element matching.

Event selection cuts
The Tevatron scenarios involve pp collisions at centre-of-mass energy √ s = 1.96 TeV.
Events are clustered with the SISCone jet algorithm [67] (similar to the MidPoint algorithm [49] that is in widespread use at the Tevatron, but infrared safe), with a jet radius R = 0.7 and a split-merge overlap threshold f = 0.75. The two hardest (highest-p t ) jets in the event should have rapidities |y| < 0.7. Events are accepted for a low-p t sample if the hardest jet has p t > 50 GeV, while they are accepted for a high-p t sample if the hardest jet has p t > 200 GeV. As concerns the event shapes, the central region is defined by η C = 1.
The LHC scenarios involve pp collisions at a centre-of-mass energy √ s = 14 TeV. 11 Events are clustered with the k t jet algorithm [47,48], with a jet radius R = 0.7. The two hardest jets in the event should have rapidities |y| < 1. Events are accepted for a low-p t sample if the hardest jet has p t > 200 GeV, while they are accepted for a high-p t sample if the hardest jet has p t > 1 TeV. As concerns the event shapes, the central region is defined by η C = 1.5. The larger choice than at the Tevatron reflects the LHC detectors' larger overall rapidity coverage. The cross sections for the different selections are given in table 2. These, and all other NLO calculations presented here, have been carried out with nlojet++ 3.0 [21] (modified to provide access to parton flavour information up to O (α 3 s )), with CTEQ6M Parton Distribution Functions (PDFs) [68] and FastJet 2.3 [69] for the jet clustering. The renormalisation and factorisation scales have central values µ R , µ F = p t ≡ (p t1 + p t2 )/2, where p t1 and p t2 are the transverse momenta of the hardest and second hardest jet respectively. The quoted errors correspond to the uncertainties due to scale variation Given the cross sections in table 2 one easily concludes that with available (anticipated) luminosities at the Tevatron (LHC) there will be large event samples on which to study event shapes. One also observes that NLO corrections are larger than one is used to seeing for (say) inclusive jet cross sections. This is a consequence of our selection based on the value of p t1 . A selection based on the average p t of the two hardest jets would instead have given K-factors rather similar to those for the inclusive cross section. 12 11 The LHC will initially run at centre-of-mass energies that are below √ s = 14 TeV, though the exact energy of collisions is subject to uncertainty and will vary over the course of the initial runs. Given that the generation of the NLO results for a single combination of collider energy and event-selection cuts requires many CPU-years of computing time, we have decided to remain with √ s = 14 TeV as our default choice for the time being. The general picture as it applies to other centre-of-mass energies can be largely understood by interpolation between the Tevatron and 14 TeV LHC results. 12 The choice of a cut on p t1 was originally motivated by the observation in the context of HERA [70,71,72] that identical simultaneous cuts on p t1 and p t2 led to poor convergence of the perturbative series, for reasons discussed in [73,74]. Cutting on p t1 was intended as a way of avoiding this problem, but, as we see here, seems to introduce issues of its own. Note that it is probably not advisable to introduce a  Table 2: Cross section for events that pass the selections cuts described in the text. The uncertainty is that due to scale variation with the choice where p t is the average of the transverse momenta of the two hardest jets. Also shown is the breakdown (at LO) into the main scattering channels; q denotes both quarks and antiquarks, and channels that contribute negligibly, such as gg → qq, are not shown. Table 2 also shows the breakdown into the 3 main partonic scattering channels (as calculated at LO). At each of the colliders, for the lower p t cut, channels involving gluons are dominant, while for the higher p t cut channels involving quarks play a bigger role. This difference between low and high-p t samples will be clearly visible in the final results.

Resummed results and uncertainty studies
Here we present resummed results for the global thrust minor, T m,g together with a study of its perturbative uncertainties at the Tevatron for the high-p t sample (p t1 > 200 GeV). Fig. 1a illustrates the NLL+NLO matched distribution obtained in the log-R matching scheme with X = 1. The renormalisation and factorisation scales are chosen, eventby-event, to be µ F = µ R = p t = (p t1 + p t2 )/2 in both the resummation and the NLO calculation. The distribution has a peak at small event-shape values that is characteristic of all resummed (and physical) event-shape distributions.
The uncertainty on the prediction is almost as important as the result itself, especially as it will allow us to gauge the significance of any disagreements that we will see with other predictive methods.
The most widely used form of uncertainty estimate is the variation of renormalization and factorization scales. The solid (red) curves in fig. 1b illustrate the effect of varying these scales simultaneously, showing the ratio of results with µ F = µ R = p t /2 and µ F = µ R = 2p t to the default result. Except at very small event-shape values or at very large ones, where the distribution vanishes, one sees that the impact of symmetric scale variation is only about 5%. Asymmetric scale variations are shown by the dashed (green) curves, corresponding to µ F = {p t /2, 2p t } while keeping µ R = p t , and µ R = {p t /2, 2p t } while keeping µ F = p t . For moderate and large values of the event shape they have a significantly larger impact than symmetric scale variations, of the order of 10% for moderate T m,g . This staggered cut on p t1 and p t2 , e.g. p t1 > 50 GeV and p t2 > 40 GeV, because that introduces an extra small parameter in the problem, related to the difference between the two p t cuts. highlights the importance of considering both symmetric and asymmetric variations. Fig. 1c shows the impact of varying X in eq. (3.24), with the line thickness increasing from X = 0.5 to X = 2. As discussed in Sec. 3.2, this variation can be used to estimate the effect of higher order logarithms not included in our NLL resummation. We find that for moderate and large values of T m,g the effect is similar in size to the asymmetric renormalization/factorization scale variation. Closer to the peak of the distribution (where the bulk of events sits), the impact of the X-scale variation is mildly larger. We also note that the variation is quite asymmetric: smaller X values distort the central distribution much more than larger values.
Finally, in Fig. 1d we estimate uncertainties that arise from the details of the matching procedure. In particular we show the ratio of the mod-R matched distribution to the log-R (see eqs. (3.26), (3.27) and (3.28)). It is clear that at large values of T m,g , the difference between the two matched distributions is large, with differences of up to 45% for T m,g ∼ 0.5. These very large discrepancies occur however only in the tail of the distribution, where few events are present. Comparison to NLO at high T m,g (not shown) indicates that of the two matching schemes, log-R matching is the one with smaller higher order terms (i.e. its NLO+NLL result is closer to NLO) at large T m,g . In the following we will therefore use log-R matching as our default.
The above findings are representative of the results for the other event shapes considered here, both at the Tevatron and at the LHC and for the low-and high-p t samples (further NLO+NLL results are shown in Sec. 4.4 and on the website associated with this article [41]). In particular symmetric renormalization and factorization scale variation (as is currently done in many phenomenological studies) systematically underestimates the true size of theoretical uncertainties. While detailed error-estimate studies such as variation of Xscale and matching procedure are possible only for specific (resummed) calculations, an asymmetric µ R and µ F variation can be carried out for generic observables. This is perhaps most relevant for multi-scale observables, where the scale of α s is a priori not clear.

Comparison of resummed, NLO and matched results
In this section we compare various levels of fixed order calculations (LO, NLO), pure resummed ones (NLL) and matched ones (NLO+NLL) at the Tevatron for the high-p t sample. Because NLL+NLO resummations are rarely available (e.g. they are currently not available for non-global observables), we discuss in particular the extent to which NLO alone can be used to describe event shape distributions. As in the previous section we will use T m,g to illustrate our findings, but results are fairly independent of the specific event shape, the collider and the details of the hard cuts. Fig. 2 shows the result for the log-R matched T m,g distribution compared to pure resummation (left) and pure NLO, and LO for reference (right). For the matched resummed result, the band corresponds to the span of all the uncertainties shown in fig. 1, while in the fixed order calculations it corresponds to the variations of just the renormalisation and factorisation scales (p t /2 ≤ µ R , µ F ≤ 2 p t , with µ F /2 ≤ µ R ≤ 2 µ F ); for the pure resummed result the band corresponds to the renormalisation, factorisation and X-scale uncertainties. As expected, the matched distribution agrees with the NLO results at large values of T m,g . However for the pure NLL resummation without any coefficient function obtained from eq. (3.23), the level of agreement with NLO+NLL is quite poor even at fairly small values of T m,g . For example, the position of the peak of the distribution is not all that well predicted (at T m,g ∼ 0.09 − 0.11 rather than at T m,g ∼ 0.08). As far as the height of the peak is concerned, both NLL and NLO+NLL distributions are normalized to one, however the NLL distribution becomes negative at T m,g > 0.35, and this negative tail causes the distribution to be far too high at low T m,g . It is on the other hand reassuring that these large differences with the matched distribution are reflected in the very large uncertainty band of the NLL distribution.
As far as the fixed-order results are concerned, Fig. 2b, they are as expected divergent at small T m,g . The LO distribution essentially never agrees with the matched distribution, while the NLO does within uncertainties for T m,g 0.2. It is on the other hand evident that scale uncertainties of the NLO results at small T m,g underestimate the size of higher order corrections not included in the fixed order calculations.
Altogether, figure Fig. 2 highlights how neither NLO nor resummation alone can provide a sensible prediction, while the combination of NLO and resummation gives significantly reduced scale-dependence compared to either on its own. Furthermore the matching procedure gives the general shape that is associated with the resummation, while maintaining the large-v behaviour of the NLO prediction. We finally note that for the event shapes

NLL+NLO matched results for a range of observables
In the previous section we established that contrary to NLO or NLL alone, NLO+NLL provides robust theoretical predictions for event shapes distributions over a large range of the event-shape values. This section contains the bulk of results of the present work: we discuss NLL+NLO resummed distributions for a number of event shapes variables, both at the Tevatron and at the LHC and for both, low-and high-p t samples. We start by looking at the effect of changing the hard selection cuts and the collider (Tevatron vs LHC) for the same observable discussed previously, T m,g , Fig. 3. There is a striking similarity between the Tevatron (left) and the LHC plot (right), both for the low-p t (Tevatron, p t1 > 50 GeV and LHC, p t1 > 200 GeV) and for the high-p t (Tevatron, p t1 > 200 GeV and LHC, p t1 > 1000 GeV) samples. We also notice that low-p t curves are broader and peaked at a higher value of T m,g . This is a consequence of the higher prevalence of gluons in both the initial and final states of the hard scattering.
Because of this similarly between low-and high-p t samples at the two colliders, we examine results for a large range of observables, as defined in Sec. 2, just for the high-p t cuts at the Tevatron, Fig. 4, and for the low-p t cuts at the LHC, Fig. 5. For each observable, we give two uncertainty bands, one corresponding to a symmetric scale variation (hatched, dark blue) and one defined in terms of all theoretical uncertainties as discussed in Sec. 4.2 (solid, light blue). Comparing Fig. 4 and Fig. 5 we see that, as observed earlier for T m,g , the peaks of the distributions are further to the right and the distributions are broader for the LHC (low p t ) than for the Tevatron (high p t ). Looking at specific observables we see that, as already remarked in the case of T m,g for all observables the symmetric scale variation uncertainties are considerably smaller than the full uncertainties, and we stress that only the latter are really indicative of the size of all kinds of neglected higher order   terms. Some final remarks concerns the NLO+NLL results for T m,R , B T,R and B W,R . As discussed in [23] and at the end of Sec. 3.1.2, recoil variables are more difficult to resum than other variables, because in the caesar approach, the NLL term g 2 (α s L) has an unphysical divergence at small values of the observable. This difficulty is reflected in the substantially larger uncertainty bands for these observables than for the directly global variants and those with an exponentially-suppressed forward term. Among the recoil variables, the thrust minor and broadenings were the only ones for which an even partially acceptable result could be obtained. In order to obtain results for recoil variables of similar quality to those for the other observables requires a resummation of initial state emissions in appropriate Fourier transform variables, as done e.g. for the Drell-Yan p t resummation [75], mixed with a Sudakov type resummation, as was done for the DIS broadening [76]. This is beyond the scope of caesar. Another characteristic to be commented on is the spike for B T,R and B W,R near 0.37. We believe this could be related to a Sudakov shoulder [77] type phenomenon, and similar (though less pronounced) artifacts have also been observed in DIS event-shape distributions. Again, it is beyond the scope of caesar to resum the enhanced higher-order terms associated with these structures.

Naive exponentiation of NLO
In the previous Section we presented full NLO+NLL resummations for a range of event shapes. Both the NLO Monte Carlo calculation and the NLL resummation are highly CPU intensive and are usually both run across many CPUs. While the NLO part of the calculation is the most computer intensive, this is to some degree counterbalanced by the fact that many observables can be computed in the same NLO run. The NLL resummation on the other hand requires essentially a separate run with caesar for each observable. Altogether, a single combination of collider energy and event-selection cuts requires many CPU-years of computing time. NLO+NLL resummation also requires that the NLO total cross section and LO distributions be decomposed into flavour channels, and this information is not available in the public version of nlojet++ (nor in most other public NLO codes). Furthermore, caesar is currently not public, the range of observables that can be resummed with caesar is not as broad as one might like (see Sec. 3.1.1), and the matching procedure at hadron colliders is not as straightforward as in e + e − , as discussed in Sec. 3.2. For the above reasons, it is interesting to explore the possibility of obtaining predictions with accuracy close to NLO+NLL using publicly available NLO results only. For instance the following combination of LO and NLO integrated, flavour summed distributions, Σ 1 (v) andΣ 2 (v), and the corresponding total cross sections, σ 0 and σ 1 , as obtained directly from nlojet++, • the fixed order expansion of the corresponding differential distribution, up to relative order α 2 s reproduces the normalized NLO differential distribution; • the formal accuracy is not even LL Σ : starting from order α 3 s the terms α n s L 2n are only correct in the limit in which the event sample is dominated by a single colour channel.
Though the method does guarantee any formal resummation accuracy, it is still instructive to see how it fares in practice. We therefore show in fig. 6 a comparison of NLO+NLL matched results with full uncertainties and the naive exponentiation of NLO, as defined in eq. (4.1) with full scale uncertainties for a representative set of observables for the Tevatron with a 200-GeV jet p t cut (similar results hold also at the LHC and for other cuts as well as other observables [41]). We see that for the directly global thrust minor (top, left) the exponentiation result is well-contained in the uncertainty band of the NLO+NLL matched result, suggesting that the naive exponentiation of NLO is indeed a quite reasonable procedure (similar results hold in general for global observables). The same observation is true also for the total broadening with recoil term (top, right) but with one important difference. In this case the NLO+NLL uncertainty band at small values of the observable is divergent, signaling the breakdown of the resummation (as discussed in [23] and at the end of Sec. 3.1.2). The scale uncertainty band of the naive exponentiation does of course not account for this and is small for all values of the observable. Therefore for observables like the recoil event shapes, whose double logarithms do not fully exponentiate, the naive exponentiation of the NLO results can only be used as long as one is far from the divergence of the pure NLL resummation (whose position is one of the pieces of information provided by caesar).
Finally, we show in fig. 6 (bottom, left) how well the naive exponentiation does in the case of the three-jet resolution parameter with exponentially suppressed term. We see that in this case the naive exponentiation result is not contained in the full uncertainty band of the NLO+NLL resummation. This is true for the tail of the distribution, where there seems to be too little radiation suppression and, similarly, for the peak, whose position is slightly displaced to the left. This softer spectrum is a feature of other exponentially suppressed observables as well (though the effect is not as remarkable).
Altogether it seems that naive exponentiation is a sensible procedure to extend the range of validity of pure NLO predictions. However, since it is not guaranteed to provide any formal logarithmic accuracy, one should rely on full NLO+NLL predictions for precision studies. In any case, we stress that before carrying out this exponentiation procedure, one should understand the basic soft/collinear properties of the observable (be it with caesar or in whatever other way).

Comparison with (matched) parton showers
For most practical applications, it is far more convenient to use parton-shower Monte Carlo event generators, like Herwig [3] or Pythia [4], or event generators merged with LO matrix elements, rather than a full NLL+NLO calculation. It can however be difficult to estimate the accuracy of these tools and the reliability of the error estimates that come with them. The purpose of this section is therefore to compare the NLL+NLO results with partonshower based predictions (at parton level, in order to avoid non-perturbative corrections from hadronisation and the underlying event).
We will start with Herwig (v6.5) events showered from exact tree-level matrix elements for 2 → 2, 2 → 3 and 2 → 4 partonic scatterings, as generated with Alpgen [78]. We use the MLM prescription [79] to avoid double counting between emissions generated in the hard 2 → n scattering and those generated by the parton shower. This combination of parton-shower and tree-level matrix elements is the standard tool for many Tevatron and LHC predictions.
To gauge uncertainties in the resulting matched samples we shall simultaneously vary the renormalisation and factorisation scales in the tree-level matrix elements by a factor of two around their default settings in the MLM procedure (which are taken as in the CKKW procedure [80]). The MLM procedure also involves a separation scale between the region of phase-space to be accounted for by the tree-level matrix elements or by the parton shower. We take this separation scale to be 0.5p t,min when we look at event samples whose hardest jet has p t > p t,min and the angular distance for a jet and a parton to be matched is restricted to be ∆R < 1.05. The hard events have been generated with a p t threshold 0.4p t,min for all partons (constrained to have |y| < 5), which must separated from each other by a distance ∆R > 0.7. In principle the MLM separation scale should also be varied in order to gauge uncertainties. However, the generation threshold should also be kept lower than the separation scale and the production of the 2 → 4 tree-level event samples with the 0.4p t,min threshold already turned out to have very low efficiency and would have become prohibitive with much lower a threshold. Therefore we will only show results with a fixed separation scale. The PDF that we used was CTEQ5L [81], the default choice in Alpgen, but we also verified that the effect of switching to CTEQ6M (as used in our NLO+NLL calculations) was small. A comparison of the parton-level Alpgen+Herwig (Tree+PS) results with the NLL+NLO results is given in figs. 7 and 8, for pp collisions at Tevatron and LHC, with a cut of 200 GeV on the transverse momentum of the hardest jet. The Alpgen+Herwig results are shown as red cross-hatched bands. The NLO+NLL results are shown as two bands: a blue hatched band whose width corresponds to the uncertainty from just the symmetric variation of renormalisation and factorisation scales; and a cyan, solid band corresponding to the full set of uncertainties represented in fig. 1.
Generally, there is reasonable agreement between the Tree+PS and the NLO+NLL results. One feature of note is that the Tree+PS uncertainty band is significantly narrow than the NLO+NLL band (even that with just symmetric scale variation). It is not immediately obvious that this truly reflects smaller uncertainties in the Tree+PS case, which, based as it is on LO calculations, would be expected to show larger uncertainties than the NLO+NLL prediction. We tend instead to interpret this as indicating that symmetric scale variation does not provide a good estimate of the true uncertainty on Tree+PS predictions. 13 There does not seem to be a clear pattern to the cases where there are significant differences between the two kinds of predictions. For example, for the τ ⊥,g and the T min variables the Tree+PS predictions seem harder than the NLL+NLO results. Instead, for y 3,E , the Tree+PS results are generally softer.
We also show in figs. 9 and 10 a comparison between results obtained from different shower Monte Carlo event generators with and without matching, for the same subset of the observables at the Tevatron and the LHC with a minimum p t on the hardest jet of 200 GeV. Pythia 6.4 is shown both for the old (virtuality ordered) and new (transverse-momentum ordered) showers. All results are shown at parton-level, without multiple interactions (i.e. In general Herwig's angular-ordered shower and Pythia's virtuality ordered (old) shower give results that are quite similar (or slightly harder) to the full matched results, with deviations visible in some cases, e.g. for B W,E at the LHC.
It is perhaps surprising that unmatched, plain parton shower results, are often harder and sometimes even closer to the NLL+NLO band than the Tree+PS matched ones, this is particularly evident for y 3,g at the LHC. This is an unexpected result as the motivation for carrying out Tree+PS merging is that parton showers are unable to reproduce the structure of hard large-angle emissions.
What is also evident from figs. 9 and 10 is that there are big discrepancies between the newer, transverse-momentum ordered shower in Pythia 6.4 (in the S0A tune) and, essentially, everything else. These distributions appear to be significantly softer than those from other parton showers, with the difference most visible in the case of the y 3 variables, and inconsistent with the NLL+NLO calculation.
It is therefore useful to try to understand whether the origin of the discrepancies lies in the new shower or in the tuning of the Pythia parameters. To further probe this issue, we show in Figure 11 a comparison among plain Herwig, virtuality ordered (DW) Pythia and different tunings of the new transverse-momentum ordered shower, S0A [82], as used above, and two more recent tunes, Perugia0 [83] Pro-pt0 [84], (shown with version 6.421; version 6.412 yields identical results). Of these the two more recent p t -ordered tunes, Pro-pt0 gives results very similar to S0A, while Perugia0 is closer (though not identical) to the Herwig and virtuality-ordered Pythia results. The conclusion to be drawn from these results is that for transverse-momentum ordered showers, the shower parameters can have major implications for the reliability of the results and a consensus has yet to emerge among current tunes for the choices of these parameters.

Non-perturbative effects
So far we discussed only perturbative effects, however, before any comparison to data can be done, non-perturbative effects have to be included. As in e + e − annihilation, there are non-perturbative effects due to hadronisation, i.e. related to the transition of partons to hadrons. In hadron-hadron collisions there are also interactions between the two beamremnants, the so-called underlying event (UE). Both effects are suppressed by inverse powers of the hard scale p t of the high-p t scattering. For event shapes the dependence is linear in 1/p t [85], while for jet-resolution parameters the dependence is even more suppressed, as will become evident also from the plots presented in this section. 14 At hadron colliders, analytical predictions for non-perturbative (NP) effects are available only for a limited number of jet-observbales [85,86,87,88]. A more general way to estimate those effects is to use event generators, such as Herwig [3]or Pythia [4], which can be run at parton-or hadron-level with or without underlying event. The default Pythia underlying event model includes multi-parton interactions, while Herwig needs to be interfaced to Jimmy [89] to have a realistic modelling of the underlying event.
In Figs. 12-15, we compare parton level and hadron level results without and with UE for the set of observables discussed previously for the low-and high-p t samples at the Tevatron and at the LHC, as obtained with virtuality ordered Pythia 6.4 (DW tune).
As far as hadronisation corrections are concerned, one notices immediately that for event shapes these effects are quite large at the Tevatron for the p t1 > 50GeV sample. They systematically shift the distributions to the right and distort them (mostly squeeze them). As expected these effects decrease considerably at p t1 > 200GeV. Going from the Tevatron to the LHC, keeping the p t1 > 200GeV cut, hadronisation correction are comparable, while again they decrease when going to the p t1 > 1000GeV sample, where they are completely negligible. Since the majority of events in a sample will have jets with p t close to the p t -cut, this patten confirms the expected 1/p t scaling of hadronisation corrections. For y 3 distributions hadronisation effects follow the same pattern but are much smaller, and are already very small at the Tevatron for p t1 > 50GeV.
The effect of the underlying event on these distributions is quite different. For event shape distributions at the Tevatron for p t1 > 50 GeV, the UE broadens significantly the distributions and moves them systematically to the right. For fixed center of mass (c.o.m.) 14 A further potential non-perturbative effect is that due to in-time pileup, the additional (usually) soft pp collisions that take place during the same beam crossing as the hard pp collision of interest. Its impact can largely be eliminated by considering only events with a single primary interaction in the beam crossing. Given the huge cross-sections for the event selections outlined in table 2, this should not pose too much of a problem except, possibly, at the higher LHC p T cut.  energy, the UE decreases with p t , but contrary to the hadronisation corrections, increasing the c.o.m. energy, at fixed p t results in an increased UE activity. As for hadronisation corrections, UE effects are much smaller for y 3 distributions than for event shapes. This fact means that at sufficiently high p t one can compare perturbative predictions directly to data, without additional NP corrections. This also makes y 3 distributions (in particular the global version) suitable for direct tuning of shower parameters. We note also that different event shapes have different NP sensitivities, broadenings seem to have smaller corrections, while masses and thrust distributions tend to have larger ones. Therefore the latter seem better suited to study NP effects and to tune models of hadronisation and underlying event.
To address this last issue further, we show in Fig. 16 how different Monte Carlo showers and tunes to the same Tevatron data differ from each other for the same set of observables. Specifically we use Herwig+Jimmy 15 Pythia's virtuality ordered shower with the DW tune and Pythia with the p t -ordered shower (S0A tune). It is noticeable how the DW, S0A and Herwig+Jimmy tunes differ (sometimes substantially), despite the fact that all have been tuned to Tevatron data. For S0A in particular this is not really surprising: if perturbative predictions already differ substantially, so will full results at hadron level. However it is nevertheless instructive, because it illustrates to what extent hadron-level event-shape distributions can help constrain perturbative aspects of the shower.
Finally, in Fig. 17 we show what happens for the same p t cut at the LHC. Discrepancies between Herwig and Pythia survive (but are maybe reduced). In addition to the DW Pythia tune, we also show DWT (which was identical at Tevatron energies) and see sizable differences between them. All this suggests that event shapes have significant scope for tunes of event generators.

Multi-jet limit
One common use of event shapes is to distinguish two different classes of multi-jet events: those in which the jets cover phase space quite uniformly, as in multi-body heavy-particle decays; and those in which the jets are relatively collimated in few bunches, as induced by the collinear singularities of massless QCD multi-particle emission. For example, the ATLAS [51] and CMS [52] discussions of prospective analyses mention the use of event shapes, most notably the transverse sphericity (see below), in physics studies that range from tt analyses to searches for supersymmetric particle decays and black-hole decays (yet other applications include hidden-valley searches [91,92]).
The purpose of this section is to compare various event shapes' ability to distinguish characteristically different event topologies. Firstly we shall examine to what extent current event shapes are able to distinguish symmetric 3-jet topologies from symmetric multijet topologies. The main finding will be that they discriminate principally between 2-jet pencillike events and symmetric events, regardless of the number of jets in the latter. We shall then study the robustness of the identification of symmetric events: both with respect to parton showering and to the orientation of the multijet system. The results from these two studies will then lead us to propose event shapes that should have enhanced sensitivity to the symmetric multijet limit. 16

The transverse sphericity
One event shape that we have not considered so far is the sphericity. Since it is by far the most widely used for discriminating symmetric multi-jet topologies, let us briefly examine its properties. It is defined in terms of the eigenvalues λ 1 ≥ λ 2 of the transverse-momentum tensor: It has the property that it tends to 1 for events with circular symmetry in the transverse plane, and is 0 for pencil-like events. However the appearance of a sum of squared momentum components in M xy makes this observable collinear unsafe, as is the case [44] for the related variable in e + e − : for example, if a hard momentum along the x direction is split into two equal collinear momenta, then their combined contribution to i p 2 xi will be half that of the original momentum. Therefore collinear splittings change the sphericity by a factor of order 1. One consequence of this is that it is impossible to make perturbative predictions for the sphericity beyond leading order. Another consequence is that parton showering and hadronisation significantly alter the value of the observable, limiting its ability to discriminate different topologies (at least when the input momenta are particles; often it is jets that are used as inputs). We shall see this explicitly in section 6.3.

The circular limit
The simplest instructive study that comes to mind for event shapes intended to distinguish symmetrical multi-jet events from dijet events is to examine their values V (N) for perfectly symmetrical planar transverse events with varying numbers N of momenta, This is illustrated in fig. 18. For uniformity the results have been normalised to the value V circ obtained for perfectly circular planar events (N → ∞), as given in table 3.
The only two observables with a monotonic (and trivial) behaviour are S pheri ⊥,g and F g . The remaining ones have been grouped into the left and right hand plots according to whether they peak for n = 3 (thrust-like) or n = 4 (broadening-like) and one sees that the perfectly circular limit does not give the largest value for all observables. Perhaps the most interesting observation from these plots is the modest difference between the S pheri Table 3: V circ , the values of various observables in the transverse circularly symmetric limit.
The events have been chosen to be planar -the variables other than ρ S,C and B T,C are however insensitive to the longitudinal components of the momenta.
3-particle and fully circular events -thus they are sensitive to the absence of a unique preferred transverse direction, but not to the overall degree of symmetry of the event.

Collinear safety and showered events
Let us now ask the question of how much the collinear unsafety of S pheri ⊥,g matters in practice. To investigate this we have taken a number of 2 → 3 partonic events and showered them with Herwig, using the "inclusive" MLM prescription [78] to reject events in which the showering introduces a fourth, harder jet, or other strong modifications of the event. Figure 19 shows the distribution of S pheri ⊥,g , F g and B T,C for two kinds of input 2 → 3 partonic event, which are both planar with all particles at rapidity y = 0: Event 1 (generic) Event 2 (Mercedes) p t1 = 828 GeV, φ 1 = 0 p t1 = 666 GeV, φ 1 = 0 p t2 = 588 GeV, φ 2 = 3π/4 p t2 = 666 GeV, φ 2 = 2π/3 p t3 = 588 GeV, φ 3 = −3π/4 p t3 = 666 GeV, φ 3 = −2π/3 One clearly sees that the collinear unsafe S pheri ⊥,g has much less discriminating power between the two events than do F g or B T,C (or for that matter any of the other event shapes that were shown in fig. 18). Furthermore it is clear for the Mercedes event that the peak at  Figure 19: Distribution of event shape values after showering and hadronisation for the "generic" and Mercedes input 3-parton events (for further details of the event generation, see text). The arrows indicate the values for the 3-parton events. Small overlap between the two distributions and good correspondence with the arrows are signs of a good observable.
contrast, the distributions for the other two observables are peaked close to the expected values (indicated by the arrows). This should of course be of no surprise given the collinear unsafety of S pheri ⊥,g . However, in view of the latter's widespread current use (albeit with jets, rather than particles, as inputs), we feel that the point is worth noting.

Impact of event orientation
One of the interests of event-shape studies is in identifying massive particle decays. Most of the event shapes above have the counterproductive characteristic that they give very different results for particles that decay with just transverse components (in the particle's centre of mass) or with both longitudinal and transverse components. To illustrate this, we take the generic event given above and rotate it by π/2 around the axis of particle 1, giving p t2 = p t3 ∼ 416 GeV, φ 2 = φ 3 = π and rapidities y 2 = −y 3 ≃ 0.88. We shower it, as explained above, and the resulting distributions for three event shapes (normalised to their values in the circular limit) are shown in fig. 20.
For S pheri ⊥,g and F g there is a large difference between the distributions for the generic and rotated-generic events (and similarly for e.g. S phero ⊥,g , τ ⊥,g and T min,g ). For the broadening in contrast, which we recall involves both the y and φ dispersions of particles with respect to axes in each of the two central half-regions, the generic and rotated-generic events give rather similar distributions. A similar phenomenon occurs with the invariant masses of those regions, in that masses too are sensitive to both directions of dispersion, though their intrinsic rotational invariance is in part spoiled when one normalises to Q ⊥,C as in eq. (2.17). The rotational invariance is probably in part the origin of the usefulness of "cluster-masses" in the context of hidden-valley studies [92]. Note however that masses are significantly more sensitive to (initial-state) forward semi-hard radiation than are broadenings.  Figure 20: Distribution of event shape values after showering and hadronisation for the "generic" and rotated generic input 3-parton events, as described in the text. Observables that have similar distributions for the two sets of events are likely to be more effective for identifying massive-object decays.

Increasing sensitivity to the spherical limit
As is clear from figure 18, none of the event shapes above are particularly effective at distinguishing truly spherical events from simpler multi-jet topologies, like symmetric transverseplanar events. What one has in mind when discussing spherical events is that they have significant "volume", symmetrically distributed around the event. One way of obtaining sensitivity to this is to consider the following matrix, separately for the up and down central regions of an event: where ∆y iU = y i − y C U and ∆φ iU = φ i − φ C U , and similarly for the central down region, C D . The eigenvalues λ U 1 > λ U 2 of M U have the property that λ U 1 is non-zero if there are two non-collinear particles in the hemisphere, while λ U 2 is non-zero if there are three non-coplanar particles in the hemisphere. The observable which we name "supersphero", is therefore non-zero only if there are 3 non-coplanar particles in each of the hemispheres of the event -i.e. for events that truly bear some resemblance to spherical events. For a perfectly spherical event the two eigenvalues in each hemisphere are λ 1 = π 2 /24 ≃ 0.411 and λ 2 = π 2 /8 − 1 ≃ 0.234. 17 The S 6 observable, in terms of its use of eigenvalues of a 2 × 2 matrix, relates of course to the F g -parameter of eq. (2.4) and to the event shapes studied for boosted top-quark identification [26,28]. The latter's use of a matrix defined in the plane transverse to a jet is actually quite similar to our use of a matrix defined in a central half-region.
A detailed study of the S 6 observable would benefit from comparisons of high jetmultiplicity QCD samples and multijet samples from new-physics scenarios. Such a study is beyond the scope of this paper, but would, we believe, be of interest.

Summary of main results
Given the length of the paper, and the fact that we have addressed quite a range of issues, we find it useful, before concluding, to summarise here the main results of the paper.
There are a number of reason why event shapes provide a powerful laboratory for studying of a range of aspects of strong-interaction physics at hadron colliders. From an experimental point of view, cross sections for the QCD (dijet) events on which one carries out event-shape studies are very large both at the Tevatron and at the LHC. This means that high-statistics event samples are already available at the Tevatron and can be expected early on at the LHC. Since event shapes are defined as dimensionless ratios of combinations of hadron-momenta, and since their differential distributions are also dimensionless, many experimental uncertainties are reduced. From a theoretical point of view, one of the attractive characteristics of event-shape studies is that different variables provide complementary sensitivities to a broad variety of features of hadronic events, such as the topology of the final state, the nature of initial and final-state jet-fragmentation, hadronisation, and the underlying event.
The above points motivated us to study many hadron-collider event shapes in a single context. We exploited the automated NLL resummation procedure implemented in caesar to obtain next-to-leading logarithmic resummed distributions matched to nextto-leading order exact predictions from nlojet++ for a large number of event shapes at the Tevatron and at the LHC.
The matching procedure is conceptually simple, but technically involved, as discussed in detail in Sec. 3.2. One issue is that the best logarithmic accuracy achievable once NLO and NLL predictions are available, namely α n s L 2n−2 in the expansion of the integrated distribution, NNLL Σ , can be obtained only if the NLO code provides full information about the flavour of all incoming and outgoing partons. This decomposition into flavour channels is not present in the publicly available version of nlojet++, so we used the extended version developed in [94] in order to extract this information. We also needed to use the flavour-k t jet-algorithm of ref. [54] to map the flavour of 2 → 3 events into that of an underlying 2 → 2 Born-like event. Additionally we needed the order-by-order expansion of PDF evolution, which was provided by hoppet [64]. The computing effort should also not be neglected: our directory of resummed results contains O (10000) files, and we estimate that several tens of years of CPU time have gone into the NLO and NLL calculations used here. Part of this complexity stems from our choice to consider several different classes of uncertainties associated with uncalculated higher-order terms: those from separate variation of renormalisation and factorisation scales; redefinition of the argument of the logarithm being resummed (X-scale); and two choices of schemes for combining (matching) the NLO and NLL results.
Among the questions we asked was whether this considerable complexity is needed. It turned out that the flavour decomposition had only a modest effect (cf. Appendix A). We also found that a simple exponentiation of the NLO result, as presented in Sec. 4.5, not even correct to LL or LL Σ accuracy, comes remarkably close to reproducing most of the NLO+NLL distributions (albeit not close enough that one would forgo NLO+NLL if it is available). One interpretation is that the large amount of radiation that comes from the 4 Born legs in a 2 → 2 process causes event-shape distributions to be dominated by regions where the logarithm that is being resummed is not all that large. Note, however, that plain (unexponentiated) NLO predictions are very inadequate substitutes for the full NLL+NLO result and their uncertainty bands are misleadingly small.
We studied three generic classes of event shapes: the directly global ones, those with exponentially suppressed terms and those with a recoil term. The definition of the observables is recaleld in sections. 2.1, 2.2, and 2.3 respectively. While stable numerical results could be obtained for observables belonging to the first two classes, for the last of these it was sometimes impossible to obtain numerically sensible results. This is in part due to cancellations among contributions from multiple emissions in the recoil term, which cause the resummation provided by caesar to have a divergence at small v, as explained in Sec. 3.1.2. It is also due to structures in the middle of the physical region, akin to Sudakov shoulders [77], which would require an additional resummation. Such shoulders are visible e.g. for the broadenings with recoil term in Figs. 4 and 5. These observables are also challenging experimentally because the measurement of the recoil term is affected by cancellations between large transverse momenta of the two hard jets.
A question mark that hangs over NLL resummations is that of coherence-violating logarithms (CVL, referred to as super-leading logarithms (SLLs) in the context of interjet energy flow) [39,40], terms potentially starting at α 4 s L 5 , related to a violation of coherence, whose validity was a crucial assumption in the resummations of [31]. There is a risk that this could therefore invalidate our claim of NLL accuracy for some observables. We investigated this point in Sec. 3.3, and found that the answer depends critically on the ordering parameter used in the calculation of the SLL terms. If, as in [39,40], one makes the assumption that the ordering parameter is transverse momentum, then the claim of NLL accuracy breaks down for our "exponentially suppressed" class of observable (not for the others), while NNLL Σ remains valid for all observables. If one instead assumes virtuality ordering, then both NLL and NNLL Σ accuracies should be valid for all our observables. This highlights the importance of understanding the question of ordering for SLLs, which also affects the coefficient of the α 4 s L 5 terms in [39,40] and probably requires that one go beyond the eikonal approximation that was used there. Nevertheless, practically we tend to believe that SLLs will not seriously affect our results, one reason being that we still retain NNLL Σ accuracy.
Turning to our phenomenological results, a feature common to all observables is that the shape of the distributions is strongly influenced by the ratio of quarks to gluons among the incoming partons. This is because the double-logarithmic Sudakov exponent, responsible for the position and width of the peak of the distribution for each underlying subprocess, is determined by the total color charge of the hard emitting partons. Event samples dominated by gluon scattering (Tevatron with p t1 > 50 GeV, LHC with p t1 > 200 GeV) have broader distributions than those dominated by quark scattering (Tevatron with p t1 > 200 GeV, LHC with p t1 > 1000 GeV). This is evident e.g. in Fig. 3 in the case of our representative observable T m,g and is discussed in Sec 4.3. We remark that dijet eventshapes at hadron colliders are the first case in which a change in a kinematical cut modifies the double logarithmic behaviour of the event-shape distributions. This would not be the case for event shapes in hadron-collider processes such as Drell-Yan production, or Z+jet or W +jet.
In the absence of data on the event shapes discussed here, one of the interesting uses of our NLL+NLO results is to compare them to the results of two Monte Carlo parton shower programs, this is discussed in Sec. 4.6. We considered Pythia 6.4 and Herwig 6.5 both without and (in the case of Herwig) with matching to multi-parton tree-level matrix elements (Alpgen, MLM prescription). The quality of the agreement between plain parton showers and the resummations depends significantly on the quark/gluon admixture: in quark-dominated event samples it is often adequate, while in gluon-dominated samples it is somewhat poorer. This may be a reflection of the extensive tuning of quark parton showers carried out with LEP data, while gluon parton showers have seen fewer constraints. The importance of tuning parton showers in a context with incoming beams is highlighted particularly strongly by the results of the newer p t -ordered shower in Pythia 6.4. In two tunes, S0A and Pro-Pt0, the agreement both with NLO+NLL and with other showers is quite poor; in the Perugia0 tune it improves, as can be seen from Fig. 11.
One might expect that supplementing parton showers with matching to multi-parton tree-level events (Tree+PS) should improve the agreement with NLO+NLL results. This is the case only for some of the observables. We also examined the impact of (simultaneous) renormalisation and factorisation scale variation on the Tree+PS results and found that it leads to an uncertainty estimate that is far smaller than the actual differences between Tree+PS and NLO+NLL results, as can be seen in Figs. 7 and 8. 18 This should not be surprising: in the NLO+NLL calculations simultaneous scale variation represented only a small part of the full uncertainties. Questions that remain open therefore are whether in the Tree+PS approach uncertainties can be more faithfully estimated if one examines further "handles" (independent scale variation, matching scale, etc.), and whether we would have reached similar conclusions with other matching schemes (e.g. CKKW) and programs.
From a non-perturbative point of view, we estimated both hadronisation and UE corrections using Monte Carlo event generators, as discussed in Sec. 5. As expected, hadronisation corrections decrease when increasing the p t1 -cut on the jets. They are fairly negligible with cuts of the order of 200 GeV both at the Tevatron and at the LHC, as can be seen in Figs. 13 and 14. For lower p t cuts, they shift the distributions to the right and, for some observables they squeeze them, see e.g. Figs. 12. For jet resolution parameters (y 3,g and y 3,E ) hadronisation effects are always small, just a few percent correction for p t1 > 50 GeV at the Tevatron, much smaller in all other cases. These observations are consistent with the experience obtained from e + e − and DIS event-shape studies.
As concerns the UE, there are observables for which it has a sizable effect even at p t1 > 1 TeV, most notably for the thrusts and jet-masses, as can be seen in Fig. 15. This means that these event shapes are particularly good for tunes of the UE. Jet-resolution parameters are the only observables for which the UE effects remain consistently small (a few percent for the lower p t -cut samples, even smaller for the large p t ones). They are therefore well suited for tunes of perturbative parameters of showers and in general for perturbative studies.
Finally, in Sec. 6 we examined how well event shapes can discriminate QCD-like twojet events from BSM-like multi-jet events, and how robust this discriminating power is with respect to parton shower (radiative) corrections. In general we find that event-shapes discriminate well between events with two or more than two jets, but they do not discriminate well between three or any large number of jets: the value of event-shapes does not even increase monotonically with the number of jets for symmetric events, see Fig. 18 in Sec. 6.2. On the other hand it is possible to design new event shapes, which start with six jets in the final state, as is the case for our "supersphero" S 6 event shape defined in Sec. 6.5. We believe these might be particularly promising for extracting new-physics signals that involve relatively isotropic events with high jet multiplicity. Other considerations that we examined in Sec. 6.3 include how well event shapes retain their discriminatory power after parton showering (the collinear-unsafe, but widely used transverse sphericity, whose definition is recalled in Sec. 6.1 is particularly poor in this respect); and also their sensitivity not just to transverse event structure, but also to longitudinal event structure (the broadenings do well at treating both on an equal footing).

Conclusions
In this article we have shown the first NLO+NLL (NNLL Σ ) predictions, with full uncertainty bands, for hadronic observables at pp and pp colliders. We opted to make these predictions for event shapes in the context of dijet production, bringing together calculations with caesar and a specially adapted version of nlojet++, despite the fact that the NLO+NLL matching is technically more challenging than for event shapes in other hadron-collider processes such as Drell-Yan [24] or W/Z+jet [22] production.
Several properties of the dijet process motivated our choice: it involves both initial and final-state partons; it offers the freedom to vary the proportion of quarks and gluons involved in the Born process, through the cut on the hard jets; when that cut is placed at moderate p t , dijet production involves a substantial gg → gg scattering component, offering the most accessible example of a gluon-dominated process; and the cross sections imply large event samples.
Comparisons of our results with parton-shower Monte Carlo predictions revealed adequate agreement for historic showers (Herwig 6.5, virtuality-ordered Pythia 6.4) in quark-dominated cases, while the showers were generally too hard in gluon-dominated processes. Some common tunes of the newer, p t -ordered shower in Pythia 6.4 fared noticeably worse than the historic showers. We also examined one framework for matching to multi-parton tree-level matrix elements (MLM matching of Alpgen+Herwig 6.5). Though it led to some improvements, it was not immediately sufficient to bring about systematic agreement with the NLO+NLL results. These findings illustrate how event shapes can provide substantial input to the quest of understanding perturbative QCD at hadron colliders.
At hadron level, some event shapes are subject to significant non-perturbative corrections from hadronisation and the underlying event. We saw this to be the case, for example, for the thrusts and jet masses, while other observables, notably the y 3 variants, were largely unaffected by non-perturbative effects. Studying a broad range of event shapes, as done here, therefore provides complementary information on QCD phenomena at hadron colliders at many different physical scales.
Event shapes are of interest not just for constraining QCD dynamics, but also for discriminating BSM-like multi-jet topologies from more QCD-like events. There are many interesting questions to ask about event shapes in this context. Some that we addressed here include their robustness to parton showering (the widely used transverse sphericity fares poorly), their sensitivity to longitudinal versus transverse event structure and their behaviour in the high jet multiplicity limit, where new dedicated event shapes, like the supersphero variable introduced here, can have particular advantages.
These first steps of ours in exploring the phenomenology of event shapes at hadron colliders open a window onto a broad range of possible new studies, both theoretical and experimental. We look forward to their future development.

Acknowledgments
We thank Mrinal Dasgupta, Günther Dissertori, Pino Marchesini, Lester Pinera, Peter Skands, Mike Seymour, Matt Strassler, Jesse Thaler and Matthias Weber for fruitful discussions on this subject. We thank CERN, Milano-Bicocca University, the LPTHE (UPMC Univ. Paris 6) and ETH Zürich for hospitality while part of this work was carried out. A.B. would like to thank Craig Prescott and the High Performance Computer center at the University of Florida for the use of computing facilities in an earlier stage of this work. G.Z. is supported by the British Science and Technology Facilities Council. The work of G.P.S. is supported in part by the Agence Nationale de la Recherche under contract ANR-09-BLAN-0060.

A Cross-checking fixed order and resummation
Part of the value of having separate resummed and fixed-order calculations for event-shape distributions is that they provide cross-checks as to the validity of each of the approaches. This check is usually performed by a comparison of the exact fixed order results Σ i (v) in eq. (3.5) with the expansion of the resummed result Σ r (v) from section 3.1. At small v the two results should differ order-by-order only by terms suppressed by powers of v or by logarithmically enhanced terms that are neglected within the resummation accuracy of Σ r (v).
At order α s , the distribution Σ(v) of eq. (3.2) has the expansion at small v Since fig. 21a contains large logarithms, a better visual constraint can be obtained by plotting the difference between Σ 1 (v) and its logarithmically-enhanced part H 12 L 2 + H 11 L, which should go to a constant at small v, and indeed does. By performing this exercise separately for each colour channel one can obtain the H (a) 10 individually, and can also verify that Σ obtained in this way are not precisely the ones that multiply the resummed distribution according to either of the two matching procedures described in section 3.2, because there one resums not logarithms of v but of a rescaled quantity X V v, eq. (3.24). To get an idea of the size of the O (α s ) term as it is relevant in the matched resummations, instead of plotting Σ 1 (v)−(H 12 L 2 +H 11 L), in fig. 21b we plot the difference between Σ 1 (v) and the distribution whereX V is the constant X V of eq. (3.25) computed for the reference Born configuration used for the analysis of the event-shape properties in caesar (two hard jets in the centre-ofmass frame with an angle θ * with respect to the beam corresponding to cos θ * = 0.2). The constantsH 1m and H ′ 10 are defined in terms of the H 1m so as to give equality between the middle and right-hand sides of eq. (A.3). One observes that the difference Σ 1 (v) −Σ r,1 (v) in fig. 21b (normalised to σ 0 ) goes, as expected, to a constant. That constant should be of order α s , whereas numerically it is O (1). However, we also know from table 2 that the order α s corrections can come with large coefficients. Given the α s C (a) 1 and the corresponding NLL resummations f (a) (v), one can predict the NNLL Σ terms in the α s expansion of Σ(v), i.e. terms α n s L p with n − 2 ≤ p ≤ n. Specifically, to second order in α s , we have where again the constantsH 2m and H ′ 2m are defined so as to given agreement between the middle and right-hand sides of eq. (A.5). Fig. 21c shows the exact second-order differential distribution v[dΣ 2 (v)/dv], 19 compared to v[dΣ 2,r (v)/dv] obtained from eq. (A.5). Again one sees good agreement, which is more readily verified by examining the difference between the two distributions, fig. 21d, which is supposed to be (and is) flat (the constant results from differentiation of theH 21 L term in eq. (A.5)). We also include the result that is obtained (lower points with errorbars) if one does not carry out the colour decomposition for α s C (a) 1 , but just computes α s C 1 = H 10 /σ 0 . This gives rise to a different expansion, Σ ′ r,2 (v), whose coefficient of L 2 is different from that ofΣ r,2 (v). For τ ⊥,g one notices that the corresponding difference between the exact result v[dΣ 2 (v)/dv] and the distribution v[dΣ ′ 2,r (v)/dv] exhibits a hint of a slope at small τ ⊥,g , indicating a missing α 2 s L 2 term in Σ ′ 2,r (v). Fig. 22 shows the same comparison of fig. 21 for the global thrust-minor T m,g . In this case one is not able, within errorbars, to see any difference between a resummed prediction containing α s C (a) 1 , giving the correct H 22 , and one based on α s C 1 , as is evident from Fig. 22d. This is possibly due to the fact that the difference between the full O (α s ) results and the first order expansion of the resummation, shown in Fig. 22a, is small.

A.1 Weighted recombination in NLO calculations
NLO Monte Carlo calculations for multi-jet processes are highly CPU intensive. Consequently, one carries out multiple calculations (runs), spread across many CPUs, and averages them so as to get the final result. The correct way of determining the average is to weight each run in proportion to its number of events. In practice, however, it is common for the distribution of each run to contain one or two bins that are "outliers",  The left-hand plot corresponds to the case with (event) numberbased averaging of NLO results from separate runs. In the right-hand plot, for each bin of the NLO results, each run has been given a weight inversely proportional to the square root of the error on that bin in the run. The results are for the τ ⊥,g observable, with Tevatron energy and cuts and p t1 > 200 GeV.
obviously inconsistent with the distribution as a whole, and which are a consequence of a handful of real and subtracted NLO events with very large opposite-sign weights that end up in different bins. These outliers lead to visible anomalies also in the number-weighted average and make it almost impossible to use the final distribution without some (often questionable) prescription to deal with the outlying bins.
A common alternative to number-weighted averaging is, for each bin of a run, to choose a weight that is inversely proportional to the square of the bin's error in that run. This is an option for example in nlojet++ (and is implicit also for the total cross section in programs like MCFM [95] that use VEGAS). Since outlier bins tend to have much larger errors than normal bins, they contribute little to the average, resulting in much smoother final distributions. However, the error-weighted averaging procedure introduces a bias, because there tends to be a correlation between the value in a bin and its error: for example, in event samples with positive-definite weights, it is well known that runs with larger bin values also have larger errors, and the final error-weighted average systematically undershoots the correct result. Fig. 23 shows the analogue of fig. 21, comparing event number-based and error-based weighting. At large negative values of L there is a clear slope, i.e. the bias in the errorweighted procedure causes the result to disagree with the expectations based on resummation. Only with number-weighted averaging does one obtain results like fig. 21, which show agreement between the logarithmic structure of the NLO and resummed calculations.
So as to deal systematically with the issue of outlying bins figs. 21 and 22 use a modified version of the number-based weighting, as follows. One first determines an error-based average for a bin b (w) , and a corresponding uncertainty on its contents σ (w) -this provides an estimate for the correct value. One then carries out the number-based average with the following modification: for a given bin, one excludes runs whose result is further than Nσ (w) from the b (w) (we use N = 100 for 15 runs; N should scale as the square-root of the number of runs). This then gives us a final result that is smooth and with a substantially reduced bias relative to an error-weighted recombination.
Note that in the phenomenological plots of sections 3 and 4 we have used the error-based recombination weights. On one hand the bias that it introduces is modest compared to uncertainties from subleading effects. On the other, some of our runs used Rambo [96] phase space and others the dipole [97] phase space, and this automatically privileges whichever of the two gives best convergence in a given phase-space region.

B Comment on effect of forward rapidity cut
For both generic global event shapes and those with an exponentially suppressed forward term, in order to satisfy the globalness requirement needed for the NLL resummation, we included all particles in the event, including those in the forward/backward regions. Experimentally however, it is not possible to perform measurements up to infinite rapidity. At the Tevatron the forward detector coverage goes up to y ≃ 3.5 and, at the LHC, measurements up to y ≃ 5 are viable. Theoretical arguments suggest that as long as the event-shape's value is not too small, the effect of not including forward emissions should be negligible [22], specifically if v v min , with v min given by [23] v min ∼ e −(a+b 1,2 )ηmax , (B.1) where the a and b i parameters were discussed in section 3.1.1. Examining the pure resummed distributions in [23], we came to the conclusion that the result in eq. (B.1) for v min ensured that the cutoff would usually have an impact only well below the maximum of the distributions. Here we supplement this analysis with a numerical study that investigates the impact of the rapidity cut in practice. For this purpose we compute the NLO+NLL prediction using a rapidity cut on input particles for the NLO part of the calculation and compare this to the full NLO+NLL without forward rapidity cuts. In parallel we carry out an estimate using a Monte Carlo event generator, since it is straightforward to run it with a rapidity cut. Fig. 24 shows comparisons between NLO+NLL with (solid line) and without the cut (full uncertainty band), as well as the corresponding Monte Carlo predictions obtained with Herwig (without UE) at parton level at the Tevatron (p t1 > 200 GeV). Fig. 25 shows the corresponding results at LHC (p t1 > 200 GeV). The results with the rapidity cut are always contained in the full uncertainty band of the results without. Furthermore, there is in general very little difference between the two Monte Carlo predictions, with the exception of the directly global transverse thrust, which is the observable most sensitive to forward emissions, as the weight of emissions in the forward region is large compared to that of emissions inside the jets. We note that this is also the one observable where the difference between Monte Carlo predictions and NLO+NLL is largest.