{'fraction_non_alphanumeric': 0.03171834705801742, 'fraction_numerical': 0.009517995739326515, 'mean_word_length': 4.309213572326212, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Resummed event-shape variables in DIS

We complete our study of resummed event-shape distributions in DIS by presenting results for the class of observables that includes the current jet mass, the C-parameter and the thrust with respect to the current-hemisphere thrust axis. We then compare our results to data for all observables for which data exist, fitting for alpha_s and testing the universality of non-perturbative 1/Q effects. A number of technical issues arise, including the extension of the concept of non-globalness to the case of discontinuous globalness; singularities and non-convergence of distributions other than in the Born limit; methods to speed up fixed-order Monte Carlo programs by up to an order of magnitude, relevant when dealing with many x and Q points; and the estimation of uncertainties on the predictions.


Introduction
Quantum chromodynamics is special among the theories of the standard model in that its strongly and moderately coupled domains are within reach of experiment, and indeed unavoidable at current energies. This leads to a wealth of non-trivial phenomena involving multiple particle production, divergent perturbative series and large non-perturbative effects even at formally perturbative scales. A full understanding of such phenomena remains one of the main aims of QCD, both because of their fundamental interest and because they are inextricably present in any measurement involving hadrons, be it a determination of the properties of QCD itself or even searches for new physics.
A feature of QCD is that the more exclusive an observable, the more difficult it is to predict. This is because more exclusive observables often involve several disparate scales, and higher-order terms in the perturbative series have coefficients which are enhanced by (large) logarithms of ratios of scales. Furthermore whenever the observable involves integrals over phase space which extend significantly into the infrared, or when one of the scales is close to the strong-coupling region, the observable is likely to be subject to large non-perturbative contributions.
A class of observables involving an interesting compromise between infrared sensitivity and perturbative calculability is that of event shapes and their distributions. These infrared-collinear (IRC) safe observables measure properties of the flow of energy-momentum, for example the normalised squared invariant mass of a hemisphere, ρ. When pushing such observables to the exclusive limit (in the case of ρ, small invariant masses) the perturbative distribution develops double logarithms at each order of perturbation theory (α s ln 2 1/ρ) n , which must be resummed to all orders. Furthermore because such observables are usually linear in soft momenta, there are significant non-perturbative effects, formally of order 1/Q, which cannot be ignored.
This class of observables has already been widely examined for 2-jet events in e + e − collisions and resummations of large-logarithmic terms [1][2][3][4][5][6] have been in existence for some time now. The literature is at first sight slightly confusing when it comes to defining the accuracy of these resummations. Let us use R(ρ) to denote the probability that the jet mass is smaller than some value ρ. Then we can write R as One can define leading logarithms (LL) as consisting of all terms with m = 2n, nextto-leading logarithms (NLL) as terms with m = 2n − 1, and so on (e.g. [3,7]). In this terminology, the state of the art for resummations is NNLL order, and in terms of relative accuracies on R it means that for α s L 2 ∼ 1, R is determined to within relative corrections of order α 3/2 s . For α s L 2 ≫ 1 however, any fixed order p of logarithmic resummation in this hierarchy, N p LL, loses its predictive power, because there will always be neglected terms α n s L m , with 2n − p > m > n, which are formally larger than 1. Fortunately many observables can be shown to exponentiate, which means that we are entitled to write , (1.2) (in some cases R needs to be broken up into a sum of terms, each of which exponentiates, but with different exponents). With this form of expansion, leading logarithms now refer to all terms m = n + 1, while next-to-leading logs (or single logs -SL) are those with m = n (e.g. [1]). The state of the art in this form of expansion is NLL accuracy. It means that for α s L ∼ 1 the neglected terms correspond to modifications of R of relative order α s . More generally, for observables that do not exponentiate, one can use this condition on the relative size of higher-order corrections to define NLL accuracy.
There is partial overlap between the terms needed for NNLL accuracy in R and NLL accuracy in ln R. Our aim is to guarantee both, though unless we explicitly specify otherwise, our statements of accuracy will refer to the expansion of ln R.
Of course, event-shape observables are not restricted to e + e − collisions and there are good reasons to extend their study to other processes. Many of the hadronisation models embody a concept of universality, whereby all 1/Q hadronisation effects should be governed by a restricted set of non-perturbative parameters which are both observable and process independent. In e + e − event-shapes there is good evidence in favour of universality [18], but in order to be fully established it needs to be tested in processes involving incoming hadrons as well. It is also important that the perturbative predictions for exclusive final-state properties be tested in a wider context than just e + e − , both so as to establish confidence in their general applicability and in view of the fact that they may well come to play a role in searches for new physics through studies of the radiation pattern associated with the different colour flows that are relevant in 'new physics' processes as opposed to the QCD backgrounds [19].
It was to investigate these issues in more detail that a study of NLL resummed DIS Breit-frame current-hemisphere (H C ) event-shape observables was initiated in [20][21][22]. From the point of view of studying hadronisation effects, DIS observables are particularly interesting in that there is a considerable amount of data [23][24][25] over a wide range of Q values (within the same experiment), making it possible to explicitly test the predicted Q-dependence of hadronisation effects. This was part of the motivation for the original theoretical studies [26,27] of hadronisation effects in mean values of DIS event shapes.
The DIS resummations carried out in [20,22] dealt with observables defined with respect to the boson (z) axis -two variants of the thrust (τ zQ and τ zE ) and one of the jet broadening (B zE ). In the current paper we complete our study of resummed DIS eventshapes by presenting results for three observables defined independently of the boson axis (except through the restriction that only particles in H C are measured) and by comparing all of our resummed results to data. 1 The three observables resummed in this paper are the thrust with respect to the thrust axis and normalised to the current hemisphere energy, T tE , the invariant squared jet mass, ρ E and the C-parameter, C E . Their exact definitions are as follows: Additionally, we refer to τ tE = 1 − T tE . Events where the energy in the current hemisphere is less than some value E lim are discarded. This is necessary to ensure the all-order infrared and collinear (IRC) safety of the observable, in particular to protect against divergences associated with events containing hard emissions in the remnant hemisphere (H R ) alone and arbitrarily soft emissions in H C . It is important that E lim not be too small a fraction of the photon virtuality Q, in order to avoid developing large ln E lim /Q terms. The fact that the above observables do not involve the boson axis in their definitions means that they are essentially insensitive to radiation in the remnant hemispherethis is in contrast to the boson axis observables considered previously, which were all sensitive to radiation in H R through recoil effects on the current quark. This has several consequences for the resummation: the double-logs are less strong because they only involve the suppression of soft and collinear radiation from one hard parton; a fraction α s of the time, the hard parton in H C is not a quark but a gluon, implying different colour factors for the double-logs; at the single log-level the appropriate factorisation scale for the parton distributions is Q 2 , independently of the value of the observable; and the resummation has non-global single-logs, discovered recently in [28,29] which involve multiple energy-ordered large-angle gluons in H R emitting coherently into H C .
The different aspects of the resummation are discussed in sections 2, 3 and 4 with the latter section addressing also the issue of previously neglected non-global logarithms for the boson-axis observables. This leads to our introducing the concept of discontinuously global variables such as τ zE and dynamically discontinuously-global variables such as B zE . In practice the phenomenological impact of the non-global logarithms on these observables is small.
In addition to there being a need for a resummation in the single current-jet limit (V → 0), it turns out that there are poorly convergent, log-enhanced structures in other regions of phase space, associated with discontinuities or non-smoothness in low-order distributions. These are analogous to the problems seen for the e + e − C-parameter at C = 3/4, and a resummation generally leads to a smoothing-out of these structures [30]. In most cases their impact is restricted to a relatively small region, so in section 5 we limit ourselves to illustrating the origin of the problems for the the DIS C E parameter, and giving a list of the non-smooth structures that arise for the other observables. We do not carry out any resummation of these observables around the singular structures, restricting ourselves to resummation near the exclusive boundary, V → 0, where V is the event shape in question.
Another, technical, issue that needs to be dealt with is that of determining the fixedorder distributions for all the observables. This turns out to be a very computer-intensive process if one wishes to obtain a reasonable precision and resolution. So in section 6 we present a method for speeding up fixed-order Monte Carlo programs such as disent [31] or disaster++ [32] through the reuse of the same event to calculate many x and Q 2 points at the same time.
Once one has the resummed and fixed-order predictions (matched with the procedures of [1,22]) one can then address (section 7) the issue of including non-perturbative corrections. We follow the approach of [11,12,33] based on the concept of a universal infrared-finite coupling, where the power correction is essentially a shift of the distribution by an amount which depends on a moment α 0 of the coupling in the infrared and on a calculable, observable-dependent coefficient (as determined in [22,26,27]).
With all the ingredients in place we then show, in section 8 the impact of the resummation, and discuss choices and uncertainties associated with unknown higher-orders and with certain prescriptions for the non-perturbative effects. This is followed in section 9 by comparisons to data, and fits for α s and α 0 . Finally, in section 10 we give our conclusions, summarising some of the many developments made in the course of our DIS resummation project.

Kinematics
Here we establish the kinematic properties of the different variables in which we are interested. The expressions derived will be approximate but sufficient for single-log/next-toleading log (NLL) accuracy.
In order to achieve this precision we need to consider configurations of soft and/or collinear gluons in addition to the hard initiating parton (quark or gluon). To be precise we want expressions for the observable in several kinematic configurations: • A hard particle in H C , close to the boson axis, accompanied by any number of soft and collinear gluons. Here the expression for the observable has to be correct to within terms of relative order of a power of the softness so as to account for both double and single logs.
• A hard particle (close to the boson axis) in H C along with an energetic almost collinear gluon. This region is only associated with single-logs and the expression for the observable need simply be correct to within an overall factor (because replacing α n s ln n 1/ρ by α n s ln n A/ρ, with A some fixed number, is equivalent to a NNLL change).
• A hard particle (close to the boson axis) in H C accompanied by a soft wide-angle gluon in H C . This situation is also responsible only for single logs and so it is again sufficient to have an expression for the observable correct to within an overall factor.
• A hard particle in H C away from the boson axis, accompanied by any number of particles soft and collinear to the hard parton. This kind of hard configuration is suppressed by a power of α s and so only gives NNLL terms (both in R and in ln R).
There it is sufficient to have an expression for the observable valid to within an overall factor.
In addition, for our resummation approach to be valid, the conditions for exponentiation (discussed in [6,22]) must also hold. This is the case for all the observables discussed here. What we shall actually do is derive expressions for the jet mass, which has already been considered in various related cases in [28,34,35], and then show that to the accuracy that we need, other observables are related to it by a constant factor.
We first obtain an expression for ρ relevant for many soft (and optionally collinear) emissions. We take a hard parton p, close to the photon axis and with energy p 0 ≃ Q/2, and any number of soft and collinear particles k i , with energies ω i and angles θ ip to the hard parton. We denote by E = p 0 + i∈H C ω i the total energy of all particles in H C .
With this kinematics, the jet mass is given by To extend this result to the case of a single hard collinear emission, k 1 , we need to take into account the recoil of the hard parton and the result for ρ is This is the information used for example in the resummations of [28,34,35]. An expression for the thrust in terms of the jet mass, in the limit where all emissions are soft (and at a polar angle of less than π/2 from the quark direction) and/or collinear, can be obtained by noting that in this limit the thrust axis coincides with the direction given by the sum of all three-momenta in the hemisphere, as in the e + e − case [1].
Denoting by P the total vector 3-momentum sum, we can rewrite ρ as while the thrust is given by where we have made use of the condition that none of the particles is moving backwards with respect to P . In the limit ρ ≪ 1, implying E −| P | ≪ E, we immediately see that τ tE = 2ρ E +O ρ 2 , allowing us to reuse the ρ resummation for τ tE , simply replacing every occurrence of ρ E with τ tE /2. 2 In the case of the C-parameter we need to explicitly derive expressions for it before being able to compare to the jet mass. Starting from eq. (1.5) we see that for any number of soft (and optionally collinear) emissions: while for a single hard collinear emission, k 1 , we have to take into account the recoil of the hard parton to get Comparing to the corresponding expressions for the jet-mass we see that in the presence of collinear emissions (many soft or one hard) C E = 12ρ E . For large-angle soft emissions we do not have an exact proportionality, but both ρ E and C E receive contributions proportional to ω i and so should have the same large-angle single logs. Therefore the resummed result for C E is simply that for ρ E with the replacement ρ E → C E /12. One final point to note is that because ρ E and τ tE are proportional in all situations involving a single soft or collinear emission in H C , the corresponding pure O (α s ) contributions to the resummed results (the C 1 terms discussed below) are identical, and the coefficients of the power corrections are also related by the same proportionality factor. In contrast, because of the different dependence on large-angle soft emissions, for C E one of the constant terms differs (by a δ-function in x), as discussed in appendix B and the power correction is not 12 times larger than that for ρ, but rather only 3π times larger.

Resummation
To write down a resummed result to NLL accuracy for the current hemisphere jet mass ρ we simply refer to Ref. [28]. There we established a formula for the hemisphere jet mass in e + e − annihilation and an essentially identical form applies here. We define the integrated current hemisphere jet mass distribution at a given x, Q 2 , where dσ 0 /dxdQ 2 is the Born cross section. 3 This integrated jet-mass distribution is given to NLL accuracy by the expression (3. 2) The precise analytical forms for the different elements in this equation are given in the appendices. Here we limit ourselves to an explanation of the physical meaning of the various factors. The first term in eq. (3.2) accounts for situations in which the hard particle in H C is a quark; Σ q (α s , L) is the form factor for the hemisphere mass to be less than e −L assuming just independent eikonal emission off a qq line (including hard collinear corrections). It can be written as in terms of the jet-mass distribution, J q , computed several years ago [34,35]. It contains both double and single logarithms. The factor S(α s L) embodies the single-logarithmic corrections associated with the fact that ρ E is non-global, namely that it measures emissions in H C only. This means that a constraint on ρ E does not imply any explicit constraint on emissions in H R , and large-angle gluons in H C are actually radiated off an ensemble of energy-ordered large-angle gluons in H R . Thus the pattern of large-angle gluon radiation in H C differs from independent emission off an eikonal qq line at the single logarithmic level, and this needs to be accounted for with the factor S(α s L), as discussed in more detail in section 4 (together with related, but rather subtle contributions that were neglected in our previous work [20,22]). It is to be noted that beyond order α 2 s L 2 , S(α s L) is known only numerically and in the large-N c limit.
The ρ-independentᾱ s C q 1 (x, Q 2 ) contribution is the only piece in the first term of (3.2) that is sensitive to the parton distribution functions (PDFs). It has a number of separate, though not always separately well-defined, origins. One physically identifiable contribution to C q 1 is from the configuration where one has a lone hard quark in H C and a gluon in H R . This on its own is divergent and cancels against corresponding virtual corrections and against the term from the factorisation of the parton distribution functions (PDFs). This implies that C q 1 has both factorisation scheme and scale dependence, similar to what is discussed in appendix A of [20] (Eqs. (A.9) and (A.10)). It should be kept in mind that when there is a hard large-angle gluon in H R , then the hard quark in H C will not be aligned along the boson axis nor will it have energy Q/2, and a number of the approximations made at single log level in Σ q and S will be incorrect -however this is beyond our accuracy, because it is relevant only a fraction α s of the time.
The remainder of α s C q 1 is associated with approximations made in Σ q for the pattern of real and virtual corrections for a hard gluons and (in some cases) with approximations for the dependence of the observable on soft large-angle and hard collinear emissions. We note that configurations with an empty H C (strictly E C < E lim ) are considered as having an undefined value of the jet mass, and so are not included in the integral for R and correspondingly in the definition of C q 1 . A final point to note about the α s C q 1 piece is that when considering α s L ∼ 1, as we usually do, and defining our hierarchy in terms of the series expansion of ln R, the C 1 terms are not formally needed to NLL accuracy. In contrast if we are in the region α s L 2 ∼ 1 and looking at the expansion of R to NNLL accuracy, then the C 1 terms are needed.
Let us now come to the second term of eq. (3.2). This accounts for the situation in which the hard parton in H C is a gluon. This happens in only a fraction α s C g 1 (x, Q 2 ) of events, and the form factor Σ g (α s , L) (also given in [34,35]) is related at leading-log level to Σ q (α s , L) by replacing C F with C A in the exponent. At single log level in Σ g the situation would be more complicated because we cannot make the approximation of the gluon being aligned along the boson axis or of its having energy Q/2. However as in the case of C q 1 , given that this whole term is multiplied by α s , to our overall level of accuracy we are free to mistreat single-logs in Σ g . In practice we systematically neglect all single-log contributions in this term, including those from non-global effects (which could have been represented by an extra factor S g in (3.2)).
Finally we recall that in order to derive a result for the C E and τ tE distributions one just replaces ρ by C/12 and τ tE /2 respectively throughout eq. (3.2) -the Σ q , Σ g and S functions have identical forms for all the above variables. It should however be noted that while the coefficient function C q 1 is the same for the τ tE and ρ variables, it is different for the C parameter (see appendix B).

Non-global logarithms and discontinuously global observables
In a previous paper we introduced the concept of a non-global observable, defined as an observable sensitive to emissions only in a restricted angular portion of phase space. We pointed out that at single-logarithmic accuracy, such observables receive a set of contributions associated with energy-ordered large-angle emissions, termed non-global logarithms [28]. In this section, we wish to refine that definition.

Non-global logs
The difference between a global observable and a non-global one is illustrated in fig. 1, which shows in grey the phase-space regions where emissions are forbidden if the observable's value is to be below a given value τ . In the left hand figure we consider the global observable τ zQ , whose value for a given ensemble of emissions is where the sum i runs over all soft (and collinear) particles and the transverse momenta are implicitly normalised to Q for brevity. For observables like τ zQ it is safe to neglect  Figure 1: Schematic depiction of the region in rapidity and log transverse momentum, in which emissions must be forbidden (grey shaded area) in order for an observable not be be larger than a given value τ . Left: shown for τ zQ , the photon-axis thrust normalised to Q (a global observable). Right: shown for τ tE , the thrust-axis thrust normalised to the energy in H C (a non-global observable). The yellow (light grey) diagonal bands indicate the hard-collinear limits. energy-ordered large-angle contributions as can be seen by considering the case with two emissions A and B, with k t,A ≫ k t,B and η A ≃ η B : the modification of the observable's value by emission B is negligible, k t,A ≃ k t,A + k t,B , and there will be full cancellation (at the single-log level) between real and virtual contributions for B.
A non-global observable is one such as τ tE , illustrated in fig. 1b. The observable's value is given by (cf. eq. (2.4)) i.e. there is only sensitivity to emissions in the current hemisphere, and placing a constraint on the value of the observable implies no direct constraint on emissions in H R . 4 In fig. 1b, this translates to only the current-hemisphere region being shaded. Now let us consider the situation with emissions A and B as shown in the figure, k t,A ≫ k t,B and both emissions at large angles, with A in H R and B in H C . Emission A does not contribute to the observable. Emission B would contribute, since it is in the shaded region -so it must be forbidden, corresponding to a resummation of virtual corrections. However the form factor Σ q in eqs. (3.2) and (3.3) is derived with the assumption that the pattern of soft large-angle emission is that off an eikonal qq current, whereas the presence of A at large angles in H R modifies the pattern of large-angle soft emissions in H C (and the corresponding virtual corrections) to be that from a qqg antenna.
So it is necessary to correct Σ q for the fact that the emission pattern is not that from a qq pair, but in general that from an ensemble that can contain any number of energyordered large-angle gluons in H R . At the two gluon level shown in fig. 1b the correction to Σ q is of the form where the double integration over transverse momentum is associated with the relevant k t ranges of A and B. At all orders, the contribution from non-global effects can be written as a function S(t) where t is given by This function S(t) has been calculated numerically in the large N c limit in [28] and progress in understanding the dynamics associated with its all-orders form has been made in [29,36].

Discontinuously global observables
While the definition given above for Figure 2: Analogue of fig. 1, shown for a discontinuously global observable, the photon-axis thrust, normalised to the energy in H C , τ zE . a non-global observable is sufficient to ensure the presence of non-global logarithms, the converse is not true: it turns out that there exist observables that are global, but which are nonetheless subject to non-global logarithms.
An example is τ zE whose value is (4.5) The region of phase space to be excluded for the observable to have smaller than some given small threshold is shown in fig. 2. The whole angular extent of phase space is relevant, so the observable is global. However the dependence of the observable on emissions in H R and H C is parametrically different, going as k 2 t and k t respectively. To denote this kind of observable, we introduce a new term: a discontinuously global observable. This is defined as being an observable, which, while sensitive to emissions in the whole angular reach of phase space, has a parametric dependence on k t (or energy) which changes discontinuously across one or more boundaries in angle. 5 In the case of τ zE , non-global logarithms arise through emission patterns of the kind represented by gluons A and B in fig. 2. Large-angle emission of A into H R will not contribute significantly to the observable as long as k t,A ≪ √ τ . However it will modify the pattern of large-angle emissions B in H C when k t,B ≪ k t,A . On the other hand large-angle emission B in H C must be forbidden if k t,B τ . Therefore there is an energy-ordered region τ in which we should account for the modification of the emission pattern of B by the presence of A. This gives a non-global log, which at second order in α s contributes as i.e. one quarter of the contribution in the fully non-global τ tE case, eq. (4.3). Analogously, at all orders we will have the same function S(t) as for the non-global τ tE , but with a different integration range in the definition of t: This non-global contribution was erroneously omitted in the original resummation for τ zE [20]. This ought, in principle, to have been revealed by the comparison with fixed-order calculations, however in practice the numerical accuracy of the fixed-order calculations is such that it is not possible to distinguish between comparisons with and without the non-global logs (this is specific to τ zE -other observables such as τ tE have effects which are four times larger and which therefore can be more straightforwardly distinguished in comparisons to fixed-order calculations).

Dynamically discontinuously global observables
An even more subtle situation involving non-global logarithms is that of the rather clumsily named 'dynamically discontinuously global observables.' These are observables which in a two-emission analysis (such as figs. 1, 2) appear to be global without discontinuities, but which develop a discontinuity dynamically at all orders. To illustrate what we mean, we consider the broadening with respect to the photon axis, B zE . This is given by Note that the parametric dependence on k t is the same in both hemispheres (there is a change in normalisation dependence in going from one hemisphere to another, but this is not sufficient to induce single-logarithms).
There is a subtlety however in that it is possible for | i k t,i | to be parametrically smaller than some of the individual k t,i 's in H R , if there is a cancellation in the vector sum. So for emissions in H R it is in principle possible to have k t,i ≫ B zE and one is then faced with the issue that such emissions, when they are at large angles, modify the pattern of emissions into H C . At first sight this may seem quite unlikely because of the small probability of a sufficient cancellation in the vector sum. However as discussed in [22], at sufficiently small values of the broadening the favoured mechanism for obtaining a small value of B zE is the cancellation in the vector sum, rather than the suppression of emissions. In such a situation non-global logarithms will arise.
In [22] we addressed the issue of this vector cancellation in some detail, showing that it was associated with a divergence in the pure NLL resummation and discussing how to improve the resummation to maintain the accuracy of there being only relative O (α s ) corrections in the R ′ ∼ 1 region. However that treatment is incomplete since it did not take into account the non-global logarithms that arise in this context due to the dynamically discontinuous globalness.
From a phenomenological point of view, it is perhaps fortunate that these issues become relevant only considerably to the left of the Sudakov peak, a region that we ignore in the comparisons to data shown later. Indeed in the fits to data we use the pure NLL resummed formula for the broadening, for which there are no non-global logarithms.

Phenomenological impact of non-global logs
In figure 3 we illustrate the effect of  non-global logarithms in the resummation of ρ E (other non-global observables are similar). The figure compares (unmatched) resummations with and without account of non-global logs. One sees that neglecting them leads to an overestimate of the peak-height by almost 50%, though elsewhere in the distribution, the effect is smaller. In practice, matching to fixed-order calculations reduces the impact of neglecting the non-global logarithms, but the effects are still of the order of 20 to 30%. The effect of non-global logs in the discontinuously global observable τ zE (not shown), is found to be smaller, at about the 10% level. This is consistent with the naive expectation (cf. section 4.2) of a reduction by a factor of four compared to the non-global case. In section 9 we examine the effect of the non-global logarithms in fits for α s and the non-perturbative α 0 parameter.

Singularities other than at V = 0
The main aim of this paper is to discuss the resummation of the singular behaviour of distributions around V = 0. However, as was pointed out in [30] there can also be singular behaviour inside the phase space of certain observables, the notable example in e + e − scattering being the C-parameter. In DIS the problem turns out to be much more common than in e + e − , even involving (in some cases) two divergent structures inside the phase space of a single observable. This is in part because one only counts particles in H C and hence there is a sharp angular cut-off (the interface of H C and H R ) which can introduce problems of non-smoothness in the various event shapes.
Divergences inside the phase space are associated with discontinuous or non-smooth structures at some lower order of perturbation theory. Higher order corrections act to smooth out these features, but order by order are divergent. One example of this is that in e + e − the LO distribution for the C-parameter has a step function at C = 3/4; the inclusion of all higher orders smooths out the step to give a 'Sudakov shoulder,' ∼ α s e −αs ln 2 (C−3/4) [30], though order by order the shoulder is associated with (a divergent series of singular) α n+1 s ln 2n (C − 3/4) terms. Such resummations, while interesting, are beyond the scope of this article. One should nevertheless be aware of the kinds of non-smoothness that arise in the different observables. We shall discuss in particular the DIS C E -parameter, and then give a summary of the properties of other observables.

The C-parameter
It is convenient here to use the variables z and ξ defined in [26], such that the incoming quark, outgoing quark and outgoing gluon have respectively (Breit-frame) momenta p, r and k, where The left-hand part of fig. 5.1 shows the different regions in the ξ, z plane: The C E -parameter is undefined in region D, zero in regions B and C, and in region A has the value [26]: For example C E takes its maximum value of 3/4 at In the neighbourhood of this point with the second derivatives c ξξ and c zz both negative. The resulting distribution of C E in the neighbourhood of C E = C max is a step function, in close analogy with the e + e − case [30]. We note that if the maximum is not smooth, but instead sharply peaked (e.g. like the summit of a tetrahedron), then the distribution goes smoothly to zero at the maximum, because the measure associated with a given value of V goes to zero linearly with V max − V . A second kind of feature that gives non-smooth distributions is a smooth extremum along an edge, for example at where δx ⊥ and δx represent distances away from the maximum in directions perpendicular and parallel to the edge respectively. 6 In the case of the C E -parameter, c x ⊥ > 0, while c x x < 0 and where we have already performed the integration over x . The easily observed fact that the x ⊥ integral is not dominated by the small x ⊥ region suggests that it is not sufficient to be considering only the region close to the edge-maximum in (5.8). However for C E < C c , the qualitative behaviour of dσ/dC E can nevertheless be understood by ignoring the In other words below C c there is a rapid fall of the distribution, behaving as a constant minus a term proportional to √ C c − C E just below the critical point. When approaching the critical point from above rather than from below the distribution remains smooth.
Both the step function at C max and the square-root behaviour at C c can be seen in the leading-order distribution for the C E -parameter in fig. 5. The figure shows also the LO+NLO result. The pathological nature of the structures which arise at NLO is only partially evident because of the limited resolution of the NLO fixed order terms.

Other observables
The kinds of patterns that we have discussed can be relevant for other observables -i.e. a smooth (quadratic) maximum V (ξ, z) in 2-dimensions leads to a step function in the distribution of V , while a smooth maximum along an edge in ξ, z leads to a square-root behaviour for the distribution around that edge-maximum. The sign of the coefficient of the square-root is the same as the sign of c x ⊥ .
It should be noted that features such as step functions (especially at the upper boundary of phase space) can also arise from other kinds of maxima, in particular an edge of constant height (i.e. c x = c x x ≡ 0 and c x ⊥ < 0). One such example is the τ zE observable whose distribution has a step function at τ zE = 1, which is the upper phase-space limit at all orders. Note however that, at least at NLO, there seems to be no higher-order 'pathological' behaviour associated with this particular discontinuity of the distribution, perhaps because of the different origin of the discontinuity relative to those discussed above. Given that these questions are not the main topic of this article, we shall not enter into any more detail. We do however give a summary of the non-smooth features of all the different observables discussed at various stages of this article, table 1.

Speeding-up fixed-order calculations
A necessary element in the comparison of our predictions to data is the matching of the re-summed results with fixed-order calculations. However the determination of high-resolution fixed order predictions with a reasonable precision is a significant technical obstacle. High resolution is needed, even when data are coarsely binned, because the shift (or convolution) of the distribution coming from non-perturbative corrections will displace the effective position of the bin boundaries in the perturbative distribution by an a priori unknown amount. Four programs exist for NLO final-state calculations in DIS: mepjet [37], disent [31], disaster++ [32] and nlojet [38] (not available when this project was started). mepjet has been found [39] to give substantially different results from the other programs (this is unfortunate, because it is the only one with the option of calculating Z-exchange contributions) and additionally is apparently not suited to the determination of event-shape distributions.
Accordingly at the start of this project, disent and disaster++ were chosen for generating fixed-order distributions. For certain observables however, non-negligible differences had been observed between their predictions [39]. Comparisons to resummations of τ zQ established [20] that some logarithmically enhanced terms (α 2 s L 3 in the integrated shape cross-section) were incorrect in disent and the problem has yet to be fixed. This formally leaves us with only disaster++ as an option.
But disaster++ is slow. On a computer with a 2 GHz Pentium IV processor, 7 two billion events 8 represents 50 days' computing time. Within the standard approach for such programs, one does a run at each x, Q pair of experimental interest. The resulting requirement of 400 days' computing time would have been prohibitive, especially given that at the time when the calculations were carried out, the computers available to us were a factor of two or three slower than those currently on the market.
An alternative approach is based on the observation that for fixed Bjorken-y, all that changes when going to different x and Q values is the parton distributions. The matrix element and event-shape calculations remain the same, and there is no need to repeat them. In other words one can exploit factorisation to separate out the structure function (x, Q) dependence from the rest of the problem. An algorithm which does this (accounting also for Bjorken-y dependence) is presented below. It has been implemented in a wrapper program, dispatch, which is available from http://www.lpthe.jussieu.fr/~salam/dispatch/. In practice, with disaster++ it allowed us to calculate distributions at 8 x, Q points using only about 30% more computing time than would have been needed for a single x, Q point. With disent the overhead of processing several x and Q points is larger in relative terms (because the matrix-element calculations are quicker in disent ), but we still need only about 2.5 times more computing time for 8 x, Q points than for a single one.

The dispatch algorithm
The algorithm, as currently implemented, is suitable for any observable whose definition, when expressed in the Breit frame, contains no dependence on x or on any dimensionful quantity other than through its ratio with Q (so for example there must be no dependence on the incoming proton momentum, or on the proton remnant, nor any cuts on absolute momenta, nor cuts in the lab frame). Breit-frame event shapes are a good example of such observables.
We shall present our discussion in terms of weighted cross sections σ (w) , where the weighting function will for example be the value observable itself if one wants to calculate its mean value, or a product of Θ functions so as to obtain the cross section for the observable's value to be in a certain 'bin'. For such weighted cross sections we can write where we have introduced a shorthand D (w) (x, y, Q 2 ) for the weighted differential cross section. As is standard for the total cross section we have decomposed the result into transverse+longitudinal (2) and pure longitudinal (L) weighted structure functions, F A Monte Carlo program such as disent or disaster++ returns event contributions which, after taking into account the weighting function and averaging over many events with the same value of x/ξ, correspond to the following function (again a vector in flavour space) where the n refers to the order of perturbation theory. There are two points of interest. Firstly, order by order, D n is a function only of ξ and y, but not of x or Q 2 . Therefore, for a fixed y value, the event contributions provided by the Monte Carlo at a single x and Q 2 value are sufficient to allow us to calculate D n (x, y, Q 2 ) for any x and Q 2 value, through the following convolution For comparison to data, a fixed value of y is usually not suitable. However the second point of interest in eq. (6.3) is that all the y dependence comes from the linear weighting of the two y-independent functions, C 2n and C Ln . So if we have calculated D (w) n for two values of y, say y 1 and y 2 , then we can obtain it for any third value, y 3 , as follows: So one simply carries out two dispatch runs at different y values and then uses a utility program to combine them with the correct weights so as to correspond to a fixed value of ep centre of mass energy, √ s.
This means that in principle we have to do two full runs of dispatch. However in many practical situations, for example when comparing to the H1 event-shape distributions in [24] the relevant y values all satisfy y 0.4 -choosing a low value for y 1 and a high value for y 2 it turns out that the weight for the high-y run is usually somewhat smaller than that for the low-y run, allowing one to generate fewer events at the high-y point. As an example, for the fixed-order distributions needed in this paper, we used y 1 = 10 −3 and y 2 = 0.9, and generated only one tenth of the total number of Monte Carlo events at the higher y value. So there is relatively little overhead from generating the second high-y point.

Inclusion of non-perturbative effects
All of the observables that we consider in this article receive significant non-perturbative contributions. So before comparing our perturbative resummed results to data, it is necessary to correct them for non-perturbative effects. Until a few years ago, this could only be done using the difference between parton and hadron levels in Monte Carlo event generators such as Herwig [41] or Pythia [42]. However such procedures tend to be ill-defined, because of ambiguities in the separation between perturbative and non-perturbative effects. Furthermore, though very successful in describing the data, they involve a number of parameters that must be tuned to the data, as well as many model-dependent assumptions.
When studying event shapes, it can be argued that many of the details of hadronisation (hadron multiplicities, relative abundances of different hadron species, etc.) factor out from the problem, due to the semi-inclusive nature of the observable. This makes it possible to consider the more inclusive aspects of hadronisation analytically, within much simpler and more transparent theoretical approaches.
A general conclusion from such studies is that the dominant hadronisation correction to event shapes is of order Λ/Q [9][10][11][13][14][15]. Furthermore there exists a universality hypothesis whereby the size of non-perturbative corrections should involve the product of a calculable, observable-dependent parameter (which we call c V ) and an a priori unknown parameter (which we will later call α 0 ) which should be universal, i.e. independent of the observable and of the process. This hypothesis can be tested by fitting for α 0 and α s in a range of observables and processes. 9 Substantial evidence in favour of universality has been obtained with e + e − data, using both mean values and distributions [18] and also from mean values of DIS observables [24,25]. In general distributions tend to give a much more robust test of universality. This is because with mean values one has one data-point for each Q value and the fit distinguishes α 0 and α s through their different Q dependences.
This procedure can lead to strong correlations between results for α s and α 0 which can mask inconsistencies in the universality hypothesis (or errors in the assumptions behind the calculation of the coefficient c V ). Distributions of event shapes on the other hand, involve many data points at each Q value and so have far more discriminatory power. This is one of the main motivations for extending the DIS studies to distributions. It should be kept in mind however that the extension to distributions is not completely without problems, as we illustrate below.

Means
For mean values non-perturbative effects change the observable by an amount δV , which within the approach of [11,43] is assumed to be of the form with µ I an arbitrary infrared matching scale (conventionally taken to be 2 GeV) and M ≃ 1.49 the Milan factor [27,[43][44][45]. In this approach α 0 is interpreted as being the average value of an infrared finite α s between scales 0 and µ I . The perturbative subtraction terms eliminate the double counting between perturbative and non-perturbative contributions (or equivalently they cancel the perturbative renormalon contribution) and should ensure that the resulting correction is independent of the infrared scale µ I (though α 0 (µ I ) itself is not). The renormalisation scale (µ R ) dependence of the perturbative subtraction terms has the interesting property that it tends to cancel part of the renormalisation scale dependence of the perturbative answer.

Distributions
For distributions, the situation is more complicated because non-perturbative effects can modify the whole shape of the distribution -this implies a dependence on more than just a single non-perturbative parameter. A very general relation between a perturbative distribution D PT,V (v), and the full result including non-perturbative effects, D full,V , can be written as For suitable observables (strictly speaking those that are linear in the momenta of several soft particles), it is argued that in the 2-jet region non-perturbative effects are independent of v and of α s and can be encoded through a simpler shape function f V (x) [15,46], One can rewrite this in terms of integer moments f V,n of f V (x), where D (n) P T,V is the n th derivative of D P T,V . Consistency with the mean value implies The higher moments each involve further non-perturbative unknowns. As discussed in [15,16,46] the fact that they multiply D (n) P T,V (v) ∼ 1/v n implies an expansion in powers of 1/(vQ), which in the limit of small v should be resummed. In other words formally one needs complete knowledge of the shape function. This requirement is the price one has to pay to 'gain access' to all the extra information in distributions as compared to mean values.
It is of interest to consider to what extent it is truly necessary to have complete knowledge of the shape function. Let us take a generic resummation where the normalised integrated distribution is (for simplicity we neglect all NLL terms) and define R ′ (v) = dR/d ln(1/v). The resummation is valid down to α s ln 1/v = O (1), which translates to R ′ (v) = O (1). The peak of the distribution is at R ′ = 1 + O (α s ), so let us take this as the point up to which we wish to calculate the distribution. To understand the hierarchy of the expansion (7.5) one needs to know the value of v at R ′ = 1. This depends strongly on the observable being considered, but can be written fairly generally as where for B zE and the e + e − total and wide-jet broadenings, For τ zE there is no analytical solution, but numerically, So for the broadenings (for which a shape function is in any case not sufficient to describe the non-perturbative effects) the neglect of higher moments might be justifiable down to the peak region, whereas for the observables derived in this paper this is less likely to be true. Physically, an interesting approximation is to replace the shape function by a δfunction, which just shifts the perturbative distribution [15,33], (7.10) Remaining terms associated with the full shape function would then correct D shift,V (v) as follows While this does not improve the formal convergence of the expansion compared to (7.5) (we still have an expansion in powers of 1/vQ) one expects partial cancellation between f V,n (f V,1 ) n to improve the actual convergence. Indeed for some observables in e + e − , perturbative arguments [16] explicitly suggest f V,2 = (f V,1 ) 2 and this seems also to be supported by the data. The approximation of the shift is very powerful, because one has only a single nonperturbative parameter α 0 , which is supposed to be universal between different observables, different processes and both mean values and distributions. 10 So in what remains of this article we will restrict our analyses to a shift.

Broadenings
It was stated above that the use of a v and α s -independent shape function held for linear observables. One particular case where this is known not to be true is that of the the broadenings [22,49] and analogous multi-jet observables [50]. The non-linear dependence of the broadening on multiple emissions has the consequence that the typical non-perturbative correction is roughly proportional to ln 1/θ, where θ is the angle of the current-quark with respect to the photon (or thrust) axis. Because of a correlation between θ and the broadening, one finds that a non-perturbative correction proportional to ln 1/B. This is taken into account in the fits that follow.

Hadronic mass effects
When testing the universality of 1/Q corrections, care is needed in the treatment of hadronic mass effects. Most observables are defined in terms of 3-momenta. Some however, in particular jet masses, depend in their conventional definitions on differences between particle energies and 3-momenta. This leads to universality-breaking Λ/Q corrections (enhanced by powers of ln Q) [48]. Such contributions can be eliminated by suitable redefinitions of the observables. Here we shall consider observables defined in the p-scheme, meaning that in the definition of the observable all occurrences of a particle's energy are replaced by the modulus of its 3-momentum. Of the DIS observables considered in this article, the only one to be modified is the jet mass. Phenomenologically the use of a consistent scheme has been found to be of vital importance in analyses [16,48] of e + e − measurements.
Strictly speaking in the p-scheme there remain small universality-breaking mass effects, however they are in practice expected to be negligible, and given that we do not dispose of data in properly universal schemes such as the E-scheme or a decay scheme (where measurements are 'carried out' after the decay of all massive particles) we make do with p-scheme measurements.
Of more concern perhaps is that fact that there can be universal mass effects which could give logarithmic corrections to the Q-dependence of power-suppressed effects. Because we are not able to estimate the absolute normalisation of these universal mass effects (or even their sign) it is not a priori clear how to eliminate them. However it should be noted that even 'traditional' 1/Q corrections may be subject to anomalous dimension-like corrections, coming from single-logs (α s ln Q/Λ) n associated with the fact that soft particles are emitted not off a bare qq current, but rather off a 'dressed' qq pair: i.e. the qq emits large angle soft gluons with k t 's ranging between between Q and Λ (giving powers of α s ln Q/Λ) and soft gluons with transverse momenta of order Λ are emitted coherently off that whole ensemble. The approximate rapidity-independence of the 'dressing' ought to ensure however that universality is not broken.
A final point relevant to all these anomalous-dimension-like corrections (whether associated with mass effects or coherent emission off dressed qq pairs) is that in distributions, the single logs (α s ln Q/Λ) n may acquire a dependence on the value of the observable. For example, for a continuously-global observable such as τ zQ one schematically expects occurrences of (α s (Q) ln Q/Λ) n in corrections to the mean, to be replaced by (α s (τ Q) ln τ Q/Λ) n in distributions, because the hardest (perturbative) large-angle emission must have a scale which is smaller than O (τ Q). For non-global and discontinuously global observables, the modification is likely to be more complex.
In any case it should be emphasised that even the first order correction term, α s ln Q/Λ, has yet to be calculated for any observable. So in the phenomenological analyses that follow later we shall ignore these complications. In figure 6 we show pure NLO predictions, pure NLL predictions (with modified logs) and matched NLO-NLL results (ln R scheme, as described in [22]) for the distributions of two observables: ρ E and the C E -parameter. Aside from leading to the appearance of a Sudakov peak in the exclusive limit, the resummation has a non-negligible (20%-40%) effect even for intermediate values of the observable, where the logarithms being resummed, while not formally large, can still be numerically significant. In many cases this will be the region used for fits.
It is interesting also to examine the region where the value of the observable is large. For the jet mass we see that the matched-resummed curve coincides almost exactly with the NLO curve close to the kinematical limit ρ E 1/4, where the leading contribution is O (α s ). The situation for C E is more complicated because beyond C E = 3/4 the leading contribution is O α 2 s . There one sees that the NLL+NLO curve behaves very differently from the pure NLO curve. This is almost certainly an artefact of the matching procedure, which is not intended to be used in situations where the leading term is not of order α s .

Choice of default X-scale for the logs
In a resummation which is truncated to some fixed logarithmic order (e.g. NLL) there is an ambiguity associated with the choice of logarithm to be resummed. For example for the jet mass we resummed ln 1/ρ E , while for the thrust we took ln 2/τ tE . The different factors were motivated by the kinematical relation between the jet mass and the thrust, derived in section 2. But we could just as easily have chosen to resum ln 2/ρ E and ln 1/τ zE -the single-logarithmic contribution in terms of these redefined logarithms would have a different functional form, but the overall answer would be the same to NLL accuracy.
In [22] we formalised this ambiguity by introducing an "X-scale" parameter, allowing us to redefine the logarithm to be resummed, where V 0 is the default coefficient in the numerator of the logarithm (ρ E,0 = 1, τ tE,0 = 2, C E,0 = 12).
Having introduced the parameter X one should establish a consistent procedure for setting it for a range of variables. We adopt the convention of choosing the default X = X 0 such that the resulting G 11 coefficient only contains terms associated with hard collinear contributions (the yellow bands in fig. 1). For τ tE , τ zQ , ρ E and C E (as for all the usual e + e − resummations) this corresponds to X 0 = 1, while for τ zE , X 0 = 1/2 and for B zE , X 0 = 1/ √ 2.

Summary of uncertainties
There are a number of potential sources of uncertainty on matched resummed event-shape distributions, essentially associated with unknown higher-order contributions.
We can obtain an indication of the size of higher-order corrections in a number of ways, and the effects are illustrated in figs. 7 and 8. These figures show the matched resummed distributions (ρ E and τ zE respectively) for default choices of parameters, matching procedures and non-perturbative corrections, together with the ratio to the default distributions when varying different scale and matching parameters.
The first ratio plot shows the effect of varying µ R by a factor of two, with the effect being at about the 10% level, and slightly larger for ρ E . Also shown is the effect of taking µ F = 2Q (with µ R = Q), which is in comparison negligible -it should be kept in mind of course that at this x value the scale dependence of the structure function is quite weak, however it is generally true that the resummation is much more sensitive to µ R than to µ F . The next plot shows the effect of varying X. The range to choose for X is somewhat arbitrary (as is the range for the renormalisation scale). The question has been investigated in some detail in [51] for e + e − event shapes, where the effect of variations in X has been compared analytically to known higher-order effects, for example to the size of the subleadingᾱ 2 s L (G 21 ) term. As would be expected a suitable range is somewhere between 1/ √ 2 to √ 2 and 1/2 to 2, with the actual range in [51] chosen to be between 2/3 and 3/2. This particular range has the feature in e + e − event shapes that the overall impact of the X variation is similar to that of the (standard) µ R variation, the purpose of the X variation then in part being to have access to an alternative functional form for the uncertainty. This is also the range used in figs. 7 and 8 and we note that here too, the effect is of the same order as that from µ R variation, though the functional dependence on the observable is somewhat different.
Next we investigate the effect of different kinds of matching, described in [22], related to analogous procedures in e + e − [1]. In the default ln R matching (R is given by eq. (3.2)) ln R is adjusted to have the correct O α 2 s expansion and then exponentiated to give the integrated distribution. In M matching, one calculates the difference between the expansion of R and the exact distribution at first and second order in α s , and this difference (multiplied by Σ q , so as to ensure correct behaviour in the V → 0 limit, essentially the condition that dR/dV → 0 for V → 0) is then added to R. Actually the difference between first order resummed and exact distributions will automatically go to zero for V → 0it is only at second order that this does not happen (because we do not know the values of the G 21 and C 2 terms of the resummation: in e + e − these are usually fitted to the fixed order, but this is not feasible in DIS). Accordingly the multiplication by Σ q is only really necessary at second order, which then gives M 2 matching. The plot shows the effect of changing to M and M 2 matchings, and we see that it is small.
Though not strictly an uncertainty, for ρ E we show additionally the effect of replacing the fixed order calculation from disaster++ with that from disent . The effect turns out to be negligible. The same is found to be true for most of the other observables as well, despite the fact that an examination of the logarithmic structure does in general reveal differences at O α 2 s ln 3 V . There are two main reasons why the the effect of changing to disent is relatively small. Firstly the largest differences between disent and disaster++ are in the gluon induced channel, which is relatively suppressed at this Q value. Secondly because of the structure of ln R, M and M 2 matching, any differences between the disent distribution and the correct logarithmic structure are multiplied by a Sudakov form factor, and so are suppressed by the matching procedure. 11 Accordingly in some instances in this paper, where for technical reasons we did not have suitable disaster++ fixed-order results, we use instead disent results (one such case is in the estimates of the µ F dependence).
In the case of τ zE , which involves the evaluation of the PDFs at a scale √ τ zE Q we show the effect of changing from LO to NLO PDF evolution (as explained in [22], at single-log accuracy both are equally valid, though for the matching and C 1 terms to be correct, the starting distribution at scale µ F must be an NLO parton distribution). At this value of Q, x = 0.11 and scaling violations are small, and as a result there is very little effect in switching to NLO PDF evolution. At other Q values (with correspondingly different x values), the effects are visible, but still relatively unimportant. Another aspect of matching is the use of modified logs, where, in the original formulation [1] p P T was 1. This modification of the logarithms is necessary to ensure that R(V ) goes to the correct value (total cross section minus cross section for E C < E lim ) at V max and in some cases additionally that the distribution goes to zero at V max . This modification introduces terms involving products of V p P T and ln V . Such terms are known to exist (see e.g. [2]). However only terms V n ln m V with n ≥ 1 are expected to be present, so we restrict p P T ≥ 1. Two possible values p P T > 1 are shown in figs. 7 and 8, and in the region of small variable, the p P T dependence is fairly weak, though it increases, especially for ρ E towards the upper kinematical limit. Analogously, one can modify the power correction so that it vanishes at the upper limit of the observable. This is especially relevant for observables whose distribution remains finite at the upper limit. The modification that we choose [22] is By default we have taken p N P = 2 (on the grounds that we wish to modify the simple picture of a shift as little as possible). We see that the effect of varying it between 1 and 3 is small.    For reference, the last ratio plot shows the effect of varying α s by ±0.005. It enables one to make an order of magnitude estimate of the relation between uncertainties on the distributions and a corresponding uncertainty in a fit for α s . We point out however that when carrying out a simultaneous fit for both α s and α 0 , the situation is more complicated because of the correlation between the two quantities.   [24]. Q values are given in GeV.

Comparison to data
event-shape distributions to H1 data [24]. It would be beyond the scope of our expertise to carry out a 'proper' fit, taking into account for example full correlations of the systematic errors. Rather the results of this section should be taken as indicative, and in particular the errors on our fit results are likely to be underestimated. 12 We shall examine all the measured event shapes, with the jet mass taken in the p scheme, as discussed in section 7.3.2. 13 The data is binned in Q as shown in table 2.
Formally when calculating the theoretical distribution, one ought to divide the Q and x ranges into small sub-bins, calculate the full matched and non-perturbatively shifted distribution in each sub-bin and then recombine to obtain results for the whole bin. The main practical limitation on this approach is that of calculating the fixed-order perturbative distributions for all the sub-bins, which even with the methods of section 6 would probably have required computing resources beyond those available to us. 14 Instead we simply calculate the distributions for the mean x and Q values as given in table 2. Both the theoretical and experimental distributions are normalised to the total cross section for the energy in the current hemisphere to be larger than E lim = 0.1Q.
The fits will involve two free parameters: α s (M Z ) and the non-perturbative parameter α 0 (µ I ) with µ I = 2 GeV). We aim to learn several things from them: we wish to verify the consistency of α 0 (and α s ) with other measurements (from e + e − means and distributions and DIS means). And we would also like to check that the fit results are consistent across a range of Q values. To accomplish this we will carry out two sets of fits, one for all data with Q > 30 GeV and a second one with additionally the bin 20 > Q > 30 GeV. The 20 > Q > 30 GeV region on its own has a statistical weight slightly larger than the combined Q > 30 GeV region.
The upper limit on the variable in the fit is chosen to be below the O (α s ) maximum value for the observable (which in some cases coincides with the all-orders maximum), because beyond this point the fixed-order distribution is known only at leading order and because the current matching procedures make little sense in steeply falling regions beyond the O (α s ) maximum. The lower limit on the variable (applied to the bin centre) is chosen to 12 For those wishing to carry out a more sophisticated analysis, the programs used to calculate the matched resummed distributions are available from http://www.lpthe.jussieu.fr/~salam/disresum/. 13 We thank Uli Martyn for providing us with the numerical values for the distribution of this observable. 14 The problem could probably be rendered tractable using an adaptation of methods described in [52].
scale as 1/Q. The coefficient of the scaling is taken larger for observables with larger power corrections. In some cases at the highest Q point even the lowest variable-bin is included. While at first sight it may look like one is therefore probing the distribution down to zero (which would be unsafe), what one is actually doing in such a bin is probing the integrated distribution at the upper boundary of the bin, which is potentially acceptable as long as that upper boundary is in the immediate neighbourhood, or to the right of the maximum of the distribution. The fitted distributions are com-  pared to data in figures 9 (ρ E ) and 10 (C E , τ tE , τ zE and B zE ). 15 Only the red (round) points have been used in the fit (the green squares correspond to those points that will additionally be used in the fit for Q > 20 GeV). The 1-sigma contours of the fit results are shown in figure 11a and the numerical results together with the χ 2 values are given in table 3. One sees remarkably good consistency between the observables, with the possible exception of the broadening, whose α s value exceeds the world average by just over a standard deviation (but one should keep in mind that our neglect of correlations between errors may lead to an underestimate of the errors). The results are also in very good agreement with those from e + e − mean values (see for example the p-scheme results of [48]). Furthermore they are reasonably compatible with results from fits to means in DIS [24,25] and from distributions in e + e − [18], though in these two cases the internal consistency between fits to different observables is not perfect.
Returning to the DIS distributions, for all variables except ρ E the χ 2 results are reasonable. Looking at the comparison between the ρ E data and theoretical predictions, figure 9, there seem to be problems in the central region of the distribution for for the Q = 57.9 GeV bin and also somewhat for the Q = 36.7 GeV bin. There is no evidence for any such discrepancy in the Q bins immediately above and below, so it is a problem which needs to be better understood.
It is interesting also to examine what happens when we extend our fits to lower Q. Figure 11b shows the fit results when one additionally includes the Q = 23.6 GeV bin (the green points of figures 9 and 10). Because of the larger number of events at lower Q values, the statistical weight of the Q = 23.6 GeV data points is comparable to and sometimes even larger than that from all the higher-Q bins combined. Its inclusion leads to a significant reduction in the consistency between different observables, though α 0 is still variable-independent to within ±10%.
The decrease in consistency as one one goes to lower Q, both between observables and compared to the world average for α s , could well be due to the relatively more important role of higher order contributions at lower Q, whether in the form of α s (Q)/Q or 1/Q 2   non-perturbative effects or NNL perturbative contributions. It would be of considerable interest to have higher-precision data at the larger Q values so as to be able to trace such effects in a continuous fashion and perhaps even identify their origin.
In this context we note that the apparent more limited consistency between observables in fits to DIS mean values [24,25] might be due to the inclusion of very low Q data. To establish whether this is truly the explanation, there is a need for higher statistics in order to be able to carry out a meaningful fit using just the higher-Q mean values.
A point of some interest is that it is the observables measured with respect to the photon axis that are less consistent with the world average for α s (giving larger values). The fact that these observables are directly sensitive to the transverse momentum of the incoming parton leads one to wonder whether we are seeing the early onset of some small-x effect or of an 'intrinsic' transverse momentum of the proton. An extension of measurements and theoretical predictions to a wider range of observables might therefore be of interest, in particular for observables which exist in 'thrust-axis' and 'photon-axis' variants (in analogy to the two thrusts), such as the broadenings (i.e. B tE as well as B zE ) and the thrust majors. Theoretical predictions for the latter would probably need to use the semi-numerical methods developed in [6].
As an aside, for illustrative pur-  In the right-hand part of table 4 we show the changes in α s and α 0 that follow from neglecting non-global logarithms in the resummed component of matched, power-corrected fits to the data. One sees effects of between 5% and 7% on α s for non-global variables and much less of an effect for the discontinuously global τ zE . The non perturbative α 0 parameter is in contrast relatively insensitive to non-global effects. It should be kept in mind that in these fits the matching in any case fixes up the α 2 s L 2 term, so it is only the effect of terms α n s L n with n ≥ 3 (about half the total effect) that is reflected. Finally, returning to the normal matched resummed fits, we should emphasise that we have not considered the impact of the uncertainties discussed in section 8.3. These should be examined in future studies, together with an analysis of the stability of the results with respect to variations in the fit ranges (with the current statistics this is difficult unless one includes lower-Q data), and of the impact of integrating the predictions over the whole bin rather than just taking them at the mean x, Q of the bin. Finally if and when the tools for fixed-order calculations including Z exchange and interference become available, then it would be useful to evaluate their impact on the results (the resummations themselves are not affected, but the constant terms and matching are).

Conclusions
The study of DIS event shapes has been largely motivated by the vast amount of HERA data that has become available over the past few years and the notable success of a similar study for e + e − variables. While mean values were the first properties of DIS event shapes to be examined theoretically, interest is shifting to the study of differential distributions for which there are over 500 data points available from H1. Here we shall summarise the developments that have taken place during the entire course of our study involving DIS event shape distributions.
The project of studying the DIS distributions implies the calculation of resummed predictions, matched to fixed order ones to enable one to study the distribution over a wide enough range of values. While it was initially believed that the DIS resummed predictions would offer relatively few extra complications, given the techniques in place for e + e − variables (the development of which spanned nearly a decade) and the similarities in definitions of the event shapes, this turned out not to be the case.
From a purely theoretical viewpoint the first complication encountered was that of the different treatment, of collinear logarithms associated with the parton distribution functions and the presence of an incoming gluon channel which features were naturally absent in the e + e − cases. In [20,22] we demonstrated that for boson axis variables such collinear radiation (including the incoming gluon channel) is associated with a change in scale of the structure function q(x, Q 2 ) → q(x, V n Q 2 ) where V denotes the variable. This variable dependence of the scale of the structure function had been encountered prior to work on DIS event shapes in for example the study of the transverse momentum distribution of a vector boson produced in hadron-hadron collisions (see e.g [53]).More recently it was also established to be present in certain three jet event shapes, both in DIS and hadronhadron collisions, for which we refer the reader to the last two items of [50]. In contrast for two (1+1) jet DIS event-shape variables not defined using the boson axis, such as those discussed in the present article, one just needs to factorise the parton densities as for a cross section and the relevant scale of the PDFs remains Q 2 .
While straightforward to understand theoretically, the presence of parton distributions, different incoming partonic channels, transverse and longitudinal components of the result with their different y Bj dependences and the presence of anomalous dimension matrices rather than simple numbers in the coefficients of certain single logarithmic terms, creates considerable problems of notation and particularly implementation. We dealt with and explained most of these issues in Refs. [20,22] and therefore in this article presented the results in only a skeletal form.
One other problem that we encountered during the course of our study was the lack of a suitable fixed order Monte Carlo program with which to combine (match) our results. While both disent and disaster++ were available to us, we found that our analytical resummations, expanded to NLO, indicated disagreement, in certain terms, with disent and compatibility with disaster++. We would therefore have preferred to use disaster++ for our purposes but for the fact that it was slow. We therefore developed a package (dispatch), described in Sec. 6, that considerably speeds up (by an order of magnitude) disaster++, making it possible to obtain results for several x, Q points with a precision which, from the point of view of computing time, would have been unfeasible previously.
While on the subject of fixed order computations we also developed several matching schemes to combine our results with those from NLO Monte Carlos [22]. While some of these (e.g ln R) involved extensions and generalisations of the corresponding e + e − procedures [1] (albeit more subtle in the DIS case) we also proposed, for the first time, multiplicative matching schemes (M and M 2 matching) [22].
Another development that we made in the course our study of DIS event shapes was the writing of our own PDF evolution code. This was necessitated by for example the need for flexibility in the evolution, specifically to fit α s using the same value for all components of the result rather than have part of the result 'contaminated' by the α s used in the standard (e.g. [54,55]) PDF sets. Also it was found that for our particular variables, there were problems arising from the non-smoothness of the interpolation procedure used in [54]. Our PDF evolution code enabled the discovery of bugs in the standard evolution codes (MRST and CTEQ), which have as a result been fixed [56,57].
Technical spin-offs apart, one important theoretical discovery made during our investigation of DIS event shapes turned out to be that of non-global logarithms. These have been discussed in great detail elsewhere [28,29], while in the current article we have considered their extension to discontinuously global variables, correcting (phenomenologically relatively unimportant) omissions in [20,22].
Lastly, we come to the original aim of this work: the comparisons with data. These give a strong boost to the idea of a universal non-perturbative extension to the standard coupling α s . In the case of most variables we find that the agreement with data at intermediate to larger Q values is remarkable down to very small values of the variables in question. The α s and α 0 values obtained by fitting our results to H1 data are in general consistent with the world average and e + e − values respectively. From the viewpoint of theory, the noticeable worsening of the description that occurs at low Q values is an area that may require further investigation. On the experimental front, data from ZEUS is eagerly awaited as is high luminosity H1 data. This apart, a proper error analysis in the fits to existing data should also be carried out. Lastly the impact of the various subleading uncertainties involved in our resummation on actual fits to the data should also be examined.
with R ′ being the derivative where we have used λ = β 0 α s L, L = ln 1 ρ and in the derivative R ′ have dropped subleading (NNLL) pieces involving the derivative of f 2 .
The functions f 1 and f 2 are listed below and In the above results the β function coefficients β 0 and β 1 are defined as and the constant K is given by For the gluonic contribution Σ g one only requires the leading log piece of the answer, since this contribution is already suppressed by a factor of α s which comes from the probability of having a hard gluon in the current hemisphere: where f 1,g is identical to the corresponding quark result f 1 (A.4), provided one replaces the quark colour charge C F by C A . Finally we come to the non-global resummation factor S. This was parameterised as below in Ref. [28] S(α s L) ≃ exp −C F C A The above parameterisation takes into account the exact leading order O α 2 s result computed in [28] and all subsequent orders in the large N c limit (i.e. neglecting contributions O 1/N 2 c ). The above parameterisation ought to be be accurate to within a few percent for t < 0.7 which corresponds to 1 − 2β 0 α s L ≥ 0.005 (and to within a fraction of a percent up to t = 0.4 which remains beyond the largest values of t probed phenomenologically).

B. Constant pieces
The various coefficient functions that we shall require can be written in the notation C b δa where δ = 2, L shall denote the transverse or longitudinal piece, the index a = q, g shall signify the incoming particle and the upper index b = q, g shall denote the hard particle in H R off which we perform the resummation to obtain the form factor Σ b .
One has (multiplied in (3.2) byᾱ s = αs 2π ) We find it convenient to introduce the function C 2q which is the whole constant piece (sum of contributions from current quark and gluon pieces) that appears in the one-loop calculation of the probability of having jet mass above ρ. From C 2q and C g 2q we can extract C q 2q via the relation C q 2q = C 2q − C g 2q The function C 2q is Analogously one has and 2ξ)) , (B.5) The thrust variable τ tE has identical coefficient functions to those listed above for the the jet mass. For the C parameter the only difference is in the piece C 2q (which also affects C q 2q ) where one needs to add the term C F π 2 6 δ(1 − ξ) to the C 2q term listed above. Note that in order to go from the coefficient functions above to the Bjorken x dependent constant pieces relevant to (3.2) we have to take 2, L pieces with the appropriate weights, convolute the functions above with the parton density functions and normalise to the Born value for F 2 . Incoming quark C b δq pieces need to be convoluted with f (x/ξ) = q,q e 2 q q (x/ξ) where the sum runs over quarks and antiquarks, and convolutions are defined as Similarly incoming gluon pieces C δg need a convolution with q,q e 2 q g (x/ξ) where g is the gluon density. To then go to the coefficient functions used in eq. (3.2) one uses relations such as where y is Bjorken-y. Finally we recall that as in [20], the coefficient functions are given in the DIS scheme. To go to another scheme one should add to C 2q and C 2g the corresponding F 2 coefficient functions in that scheme.