PDF dependence of Higgs cross sections at the Tevatron and LHC: response to recent criticism

We respond to some criticism questioning the validity of the current Standard Model Higgs exclusion limits at the Tevatron, due to the significant dependence of the dominant production cross section from gluon-gluon fusion on the choice of parton distribution functions (PDFs) and the strong coupling (alpha_S). We demonstrate the ability of the Tevatron jet data to discriminate between different high-x gluon distributions, performing a detailed quantitative comparison to show that fits not explicitly including these data fail to give a good description. In this context we emphasise the importance of the consistent treatment of luminosity uncertainties. We comment on the values of alpha_S obtained from fitting deep-inelastic scattering data, particularly the fixed-target NMC data, and we show that jet data are needed for stability. We conclude that the Higgs cross-section uncertainties due to PDFs and alpha_S currently used by the Tevatron and LHC experiments are not significantly underestimated, contrary to some recent claims.


Introduction
Discovery or exclusion of the Standard Model Higgs boson (H) at the Tevatron and Large Hadron Collider (LHC) requires precise knowledge of the theoretical cross section; see, for example, refs. [1][2][3], and references therein. Cross-section predictions for the dominant production channel of gluon-gluon fusion (gg → H) are strongly dependent on both the gluon distribution in the proton and the strong coupling α S , which enters squared at leading-order (LO) with sizeable next-to-leading order (NLO) and next-to-next-to-leading order (NNLO) corrections. In particular, the Tevatron Higgs analysis [4,5], with current exclusion at 95% confidence-level (C.L.) for a Standard Model Higgs boson mass M H ∈ [158, 173] GeV [5], requires knowledge of the gluon distribution at relatively large momentum fractions x 0.1 where constraints from data on deep-inelastic scattering (DIS) or Drell-Yan production are fairly weak. In this paper, which accompanies a separate paper [6], we respond to several (related) issues which have been raised in recent months [7][8][9][10][11][12], particularly regarding the use of parton distribution functions (PDFs) determined from limited data sets in making predictions for the Tevatron (and LHC) Higgs cross sections, as alternatives to the most common choice of the MSTW 2008 PDFs [13] used in the Tevatron [4,5] and LHC [3] Higgs analyses.
First in section 2 we demonstrate explicitly how the gg → H cross sections depend on the Standard Model Higgs boson mass M H , the gluon-gluon luminosity function and the choice of α S (M 2 Z ), by comparing predictions obtained using PDFs (and α S values) from various different PDF fitting groups. In section 3 we present a detailed quantitative comparison of the quality of the description of Tevatron jet data using different PDF sets. The MSTW 2008 analysis [13] is the only current NNLO PDF fit which includes the Tevatron jet data, providing the only direct constraint on the high-x gluon distribution. In section 4 we examine the different values of the strong coupling α S used by the different PDF groups, particularly those values mainly extracted from DIS data, and we look at the constraints arising from different sources. In section 5 we respond to recent claims [11] that the theoretical treatment of the longitudinal structure function F L for the NMC data [14] can explain the bulk of the difference between predictions for Higgs cross sections calculated using either the MSTW08 [13] or ABKM09 [15] PDFs. Finally we conclude in section 6 that MSTW08 is presently the only fully reliable PDF set for calculating Higgs cross sections at NNLO, particularly if sensitive to the high-x gluon distribution, and that the recent exclusion bounds [4,5] obtained by the Tevatron experiments are robust based upon this choice.
2 Dependence of Higgs cross sections on PDFs and α S

Dependence on Higgs mass
We show the NLO and NNLO gg → H total cross sections (σ H ) versus the Standard Model Higgs boson mass M H in figure 1 at the Tevatron (centre-of-mass energy, √ s = 1.96 TeV) and the LHC ( √ s = 7 TeV) for different PDF sets and a fixed scale choice of µ R = µ F = M H , calculated with settings given in section 4.2 of ref. [6]. At NLO [16], we use the corresponding NLO PDFs (and α S values) from MSTW08 [13], CTEQ6.6 [17], CT10 [18] and NNPDF2.1 [19], all of which are fully global fits to HERA and fixed-target DIS data, fixed-target Drell-Yan production, and Tevatron data on vector boson and jet production. At NNLO [20], we use the corresponding NNLO PDFs (and α S values) from MSTW08 [13], ABKM09 [15], JR09 [21,22] and HERAPDF1.0 [23], where in the last case no uncertainty PDF sets are provided and the two curves correspond to α S (M 2 Z ) = 0.1145 and α S (M 2 Z ) = 0.1176, with the larger α S value giving the larger Higgs cross section. For the other PDF sets, we compute the "PDF+α S " uncertainty at 68% C.L. according to the recommended prescription of each group, summarised in ref. [6]. The data sets included in the MSTW08 fit at NNLO are the same as at NLO, with the omission of HERA data on jet production, while the ABKM09 and JR09 fits only include DIS and fixed-target Drell-Yan data. The HERAPDF1.0 fit only includes combined HERA I inclusive DIS data, while the other NNLO fits (MSTW08, ABKM09, JR09) instead include the older separate data from H1 and ZEUS. However, including the combined HERA I data [23] in a variant of the MSTW08 fit was found to have little effect on predictions for Higgs cross sections [24]. The NNPDF fits parameterise the starting distributions at Q 2 0 = 2 GeV 2 as neural networks, whereas other groups all use the more traditional approach of parameterising the input PDFs as some functional form in x, each with a number of free parameters, which varies significantly between groups. Contrary to the "standard" input parameterisation at Q 2 0 ≥ 1 GeV 2 , the JR09 set uses a "dynamical" parameterisation of valence-like input distributions at an optimally chosen Q 2 0 < 1 GeV 2 , which gives a slightly worse fit quality and lower α S values than the corresponding "standard" parameterisation, but is nevertheless favoured by the JR09 authors. More details on differences between PDF sets are given in section 2 of ref. [6]; see also the descriptions in refs. [25][26][27].
The size of the higher-order corrections to the gg → H total cross sections is substantial. Taking the appropriate MSTW08 PDFs and α S values consistently at each perturbative order for σ H with M H = 160 GeV, then the NLO/LO ratio is 2.1 (Tevatron) or 1.9 (LHC), the NNLO/LO ratio is 2.7 (Tevatron) or 2.4 (LHC), and so the NNLO/NLO ratio is 1.3 (Tevatron and LHC). The perturbative series is therefore slowly convergent, mandating the use of (at least) NNLO calculations together with the corresponding NNLO PDFs and α S values. The convergence can be improved by using a scale choice µ R = µ F = M H /2, which mimics the effect of soft-gluon resummation. However, the goal of this paper is to available PDFs, specifically the HERAPDF1.0 NNLO set with α S (M 2 Z ) = 0.1145, can lower the Tevatron Higgs cross section by up to 40% compared to MSTW08 for M H ≈ 160 GeV, requiring more than twice as much Tevatron data to recover the same sensitivity as the 2010 analysis by the Tevatron experiments [4], which used MSTW08 for the central prediction. This is obviously potentially very worrying. However, figure 2(c,d) shows that the lowest cross section occurs not with either of the HERAPDF sets, but with ABKM09, where the central cross section is ≈ 75% that of MSTW08 at M H ≈ 160 GeV. The cross-section ratios for ABKM09 and JR09 in figure 2(d) seem close to those in the inset of figure 1 of ref. [9], but we do not reproduce the extreme behaviour of the HERAPDF1.0 sets. Our results are supported by those in ref. [10] where it is also observed that ABKM09 gives lower Higgs cross sections at the Tevatron than the HERAPDF1.0 set with α S (M 2 Z ) = 0.1145. One obvious difference is the scale choice µ R = µ F = M H /2 used in ref. [9] rather than µ R = µ F = M H used here and in ref. [10]. However, we have checked that the ratio of cross sections with respect to MSTW08 is largely independent of the different scale choice. The detailed arguments of ref. [9] assume the "worst-case scenario" of a 40% reduction in σ H at M H ≈ 160 GeV from the central value of HERAPDF1.0 with α S (M 2 Z ) = 0.1145, and therefore the conclusions require modification if there is a mistake in their HERAPDF1.0 calculations. 1 Nevertheless, even the 25% reduction in σ H at M H ≈ 160 GeV from the central value of ABKM09 is still a problem, as it lies well outside both the MSTW08 PDF+α S uncertainty at 90% C.L. used in ref. [4] and the PDF4LHC 2 uncertainty used in ref. [5]. (These two prescriptions for uncertainties give similar results, but the former is clearly much simpler; see section 5 of ref. [6] for more discussion.) We note that in justifying the use of the HERAPDF set, BDFG [9] make the statement: "However, HERAPDF describes well not only the Tevatron jet data but also the W , Z data. Since this is a prediction beyond leading order, it has also the contributions of the gluon included. This gives an indirect test that the gluon densities are predicted in a satisfactory way." This statement is very misleading: the W charge asymmetry and the Z rapidity distribution at the Tevatron, used as a PDF constraint, are almost insensitive to the gluon distribution, and the statement makes no reference to the quantitative comparison of PDFs to jet data. In the rest of this paper we will present a number of arguments to show that, of all the currently available NNLO PDF sets, only MSTW08 provides a fully reliable estimate of the Higgs cross sections at the Tevatron and LHC.

Dependence on gg luminosity
At LO, the PDF dependence of the gg → H total cross section is simply given by the gluon-gluon luminosity evaluated at a partonic centre-of-mass energy √ŝ = M H , where f g (x, µ 2 =ŝ) is the gluon distribution and τ ≡ŝ/s. In figure 4 we show the gluongluon luminosities calculated using different PDF sets and taken as the ratio with respect to the MSTW 2008 value, at centre-of-mass energies corresponding to the (a,b) Tevatron and (c,d) LHC. The relevant values of √ŝ = M H = {120, 180, 240} GeV are indicated, along with the threshold for tt production at the LHC, √ŝ = 2m t with m t = 171.3 GeV, where this process is predominantly gg-initiated at the LHC. Indeed, tt production at the LHC is strongly correlated with gg → H production at the Tevatron, with both processes probing the gluon distribution at similar x values, as seen from figure 4. We point out in ref. [6] that the current tt cross-section measurements at the LHC [29,30] seem to distinctly favour MSTW08 over ABKM09.
The NLO luminosities in figure 4(a,c) are shown for the global fits from MSTW08 [13], CTEQ6.6 [17], CT10 [18] and NNPDF2.1 [19]. The NNLO luminosities in figure 4(b,d) are shown for MSTW08 [13], HERAPDF1.0 [23], ABKM09 [15] and JR09 [21,22]. The two HERAPDF1.0 NNLO curves shown are for both α S (M 2 Z ) = 0.1145 and 0.1176, where the latter gives the smaller gg luminosity at lowŝ values and the larger gg luminosity at highŝ values. The larger α S value means that less gluon is required at low x to fit the scaling violations of HERA data, ∂F 2 /∂ ln(Q 2 ) ∼ α S g, therefore more gluon is required at high x from the momentum sum rule. Both these effects, larger α S and more high-x gluon, raise the Tevatron Higgs cross section and improve the quality of the description of Tevatron jet data, as we will see in section 3. The NNLO trend between groups is similar to at NLO [6]. There is reasonable agreement for the global fits, but more variation for the other sets, particularly at largeŝ, where HERAPDF1.0 and ABKM09 have much softer high-x gluon distributions, and this feature has a direct impact on the gg → H cross sections, particularly at the Tevatron (see figure 2).

Dependence on strong coupling α S
The various PDF fitting groups take different approaches to the values of the strong coupling α S and, for consistency, the same value as used in the fit should be used in subsequent cross-section calculations. The values of α S (M 2 Z ), and the corresponding uncertainties, for MSTW08, ABKM09 and GJR08/JR09 are obtained from a simultaneous fit with the PDF parameters. Other groups choose a fixed value, generally close to the world average [31], and for those groups we assume a 1-σ uncertainty of ±0.0012 [26], very similar to the MSTW08 uncertainty. The central values and 1-σ uncertainties are depicted in figure 5 as the larger symbols and error bars, while the smaller symbols indicate the PDF sets with alternative values of α S (M 2 Z ) provided by each fitting group. The fitted NLO α S (M 2 Z ) value is always larger than the corresponding NNLO α S (M 2 Z ) value in an attempt by the fit to mimic the missing higher-order corrections, which are generally positive. The world average α S (M 2 Z ) [31], shown in figure 5, combines determinations made at a variety of perturbative orders, but in most cases an increase in the order corresponds to a decrease in the value of α S (M 2 Z ) obtained. The gg → H cross sections at the Tevatron and LHC start at O(α 2 S ) at LO, with anomalously large higher-order corrections, therefore they are directly sensitive to the value of α S (M 2 Z ). Moreover, there is a known correlation between the value of α S and the gluon distribution, which additionally affects the gg → H cross sections. In figures 6 and 7 we show this sensitivity by plotting the Higgs cross sections versus α S (M 2 Z ) at the Tevatron and LHC for Higgs masses M H = {120, 180, 240} GeV. We plot both NLO and NNLO predictions for a fixed scale choice µ R = µ F = M H . The format of the plots is that the markers are centred on the default α S (M 2 Z ) value and the corresponding predicted crosssection of each group. The horizontal error bars span the α S (M 2 Z ) uncertainty, the inner vertical error bars span the "PDF only" uncertainty where possible (i.e. not for ABKM09 or GJR08/JR09, where α S is mixed with the input PDF parameters in the error matrix), and the outer vertical error bars span the PDF+α S uncertainty. The effect of the additional α S uncertainty is sizeable. The dashed lines at NLO or the solid lines at NNLO interpolate the cross-section predictions calculated with the alternative PDF sets provided by each group, represented by the smaller symbols in figure 5. The NNLO plots in figure 7 also show the NLO predictions (open symbols and dashed lines) together with the corresponding NNLO predictions (closed symbols and solid lines) to explicitly demonstrate how the size of the NNLO corrections depends on both the α S (M 2 Z ) choice and the PDF choice. It is apparent from the plots that at least part of the MSTW08/ABKM09 discrepancy for Higgs cross sections is due to using quite different values of α S (M 2 Z ) at NNLO, specifically α S (M 2 Z ) = 0.1135 ± 0.0014 for ABKM09 [15]  for MSTW08 [13,32]. Comparing cross-section predictions at the same value of α S (M 2 Z ) would reduce the MSTW08/ABKM09 discrepancy at the LHC, but there would still be a significant discrepancy at the Tevatron (see also the later table 5 in section 5).

Theoretical uncertainties on α S
In ref. [32] we gave a prescription for calculating the "PDF+α S " uncertainty on an observable such as a hadronic cross section, due to only experimental errors on the data fitted. An estimate of the theoretical uncertainty on α S was given as ±0.003 at NLO and at most ±0.002 at NNLO, where these values should be interpreted as roughly 1-σ (68% C.L.). However, this additional uncertainty was not recommended to be propagated to the "PDF+α S " uncertainty on cross sections, in the same way that theoretical errors on PDFs are not generally provided and propagated to uncertainties on cross sections. It was intended simply to be an estimate of how much the value of α S (M 2 Z ) might change if extracted at even higher orders. It has subsequently been proposed (by Baglio and Djouadi) to include the theoretical uncertainty on α S in the cross-section calculation for the gg → H process at the Tevatron [7] and LHC [8], which somewhat reduces the ap- parent inconsistency between MSTW08 and ABKM09 seen in figures 2 and 3. In figure 8 we show the effect of adding in quadrature an additional theoretical uncertainty on α S to the 90% C.L. MSTW08 PDF+α S uncertainty for the gg → H cross sections at both the Tevatron and LHC, at both NLO and NNLO, plotted as a function of the Higgs mass M H . 3 If a similar theoretical uncertainty on α S was also added to the ABKM09 uncertainty band, which includes only experimental uncertainties, then the MSTW09 and ABKM09 uncertainty bands would overlap at the Tevatron, at least in the M H range shown here. However, even if the additional α S uncertainty is applied in this manner, it is misleading to claim that it leads to more of an agreement in the predictions obtained using the two PDF sets, since variations of cross sections with α S are very highly correlated between different PDF sets. We will see in the rest of this paper that differences between groups in α S values, gluon distributions and Higgs cross sections are largely due to the selection of data fitted, and it is not the case that the discrepancies should be attributed to unaccounted theoretical uncertainties. 3 We calculate the cross sections evaluated with αS(M 2 Z ) = 0.120 ± 0.003 at NLO and αS(M 2 Z ) = 0.117 ± 0.002 at NNLO, to determine the variation due to the additional theoretical uncertainty on αS at 68% C.L., then we scale this uncertainty by 1.64485 to get the 90% C.L. theoretical uncertainty.

Constraints from jet production at the Tevatron
Here we present a quantitative study of the description of the Tevatron Run II inclusive jet data [33][34][35] and dijet data [36] by different PDF sets. The goal is to compare the description of Tevatron jet data in a similar manner to the benchmark cross-section study of ref. [6], i.e. we use the same code and settings for all NLO and NNLO PDF sets (with the correct α S value for each set) to ensure that observed differences are only due to the PDF choice rather than any other factor. We do not consider the less reliable Tevatron Run I data, which prefer a much harder high-x gluon distribution [13], and are obtained using less sophisticated jet algorithms. The three data sets on inclusive jet production from the Tevatron Run II [33][34][35] were all found to be compatible [13]. The MSTW 2008 analysis [13] included the CDF Run II inclusive jet data using the k T jet algorithm [33] and the DØ Run II inclusive jet data using a cone jet algorithm [35]. Consistency was checked with the CDF Run II inclusive jet data using the cone-based Midpoint jet algorithm [34], but this data set was not included in the final MSTW08 fit, since it is essentially the same measurement (using 1.13 fb −1 ) as ref. [33] (using 1.0 fb −1 ), differing mainly by the choice of jet algorithm. The k T jet algorithm is theoretically preferred due to its property of infrared safety, and the corresponding CDF Run II data [33] was already published and implemented in the MSTW08 analysis by the time the CDF Run II Midpoint data [34] appeared. The DØ Run II inclusive jet data [35] and dijet data [36], both defined using a cone jet algorithm, are also measured from essentially the same 0.7 fb −1 of data, differing mainly by the kinematic binning, so as with the two CDF data sets it would be doublecounting to include both in the same PDF extraction. We will concentrate on the inclusive jet data (section 3.2), but we will also make a first quantitative comparison to the more recent DØ dijet data (section 3.3). However, first in section 3.1 we precisely define the goodness-of-fit measure used for the comparison of data and theory.
One obvious problem is that the complete NNLO partonic cross section (σ) for inclusive jet production is currently unknown, and needs to be approximated with the NLOσ supplemented by 2-loop threshold corrections [37], while even these 2-loop threshold corrections are unavailable for the dijet cross section. We calculate jet cross sections using fastnlo [38] (based on nlojet++ [39,40]), which includes these 2-loop threshold corrections. Following the usual way of estimating theoretical uncertainties due to unknown higher-order corrections, we take different scale choices µ R = µ F = µ = {p T /2, p T , 2p T } as some indication of the theoretical uncertainty. Smaller scale choices raise the partonic cross section, so favour softer high-x gluon distributions [13], and the central µ = p T was chosen for the final MSTW08 fit [13]. We comment on the scale dependence in section 3.4, we present distributions of pulls and systematic shifts in section 3.5, we briefly discuss other collider data on jet cross sections in section 3.6, then finally we summarise our findings in section 3.7.

Definition of goodness-of-fit, χ 2
It is important to account for correlated systematic uncertainties of the experimental data points. The full correlated error information is accounted for by using a goodness-of-fit (χ 2 ) definition given by [41,42] where T i are the theory predictions and are the data points allowed to shift by the systematic errors in order to give the best fit.
Here, i = 1, . . . , N pts. labels the individual data points and k = 1, . . . , N corr. labels the individual correlated systematic errors. The data points D i have uncorrelated (statistical and systematic) errors σ uncorr.
i and correlated systematic errors σ corr. k,i . Minimising the χ 2 in eq. (3.1) with respect to the systematic shifts r k gives the analytic result that [41,42] and δ kk ′ is the Kronecker delta. Therefore, the optimal shifts of the data points by the systematic errors, eq. (3.2), are solved for analytically. Here we use the same notation 4 as in the MSTW08 paper [13]. We treat the luminosity uncertainty as any other correlated systematic. However, we find that the relevant systematic shift r lumi. ∼ 3-5 for some PDF sets with soft high-x gluon distributions (e.g. ABKM09 and HERAPDF1.0), which is clearly completely unreasonable, as it means that the data points are normalised downwards by 3-5 times the nominal luminosity uncertainty (around 6% for both CDF and DØ). The penalty term r 2 lumi. will contribute only 9-25 units to the total χ 2 given by eq. (3.1), which can therefore still lead to reasonably low overall χ 2 values (see appendix A for details).
It is the usual situation at collider experiments that the luminosity determination is common to all cross sections measured from a given data set (see, for example, refs. [44,45]), so the requirement of a single common luminosity is mandatory when fitting multiple measurements taken during a single running period. In figure 9 we compare NNLO predictions for the W and Z total cross sections at the Tevatron Run II, calculated in the zero-width approximation with settings described in ref. [6]; see also similar comparisons in ref. [10]. The format of the plots in figure 9 is the same as for the gg → H cross sections in section 2.3, i.e. we show the cross-section predictions plotted against α S (M 2 Z ). We compare to CDF Run II data on W [46] and Z [47] total cross sections, and to DØ Run II data on the Z total cross section [48]. The thicker horizontal lines in figure 9 indicate the central value of each experimental measurement, the thinner horizontal lines indicate the statistical and systematic (excluding luminosity) uncertainties added in quadrature, while the shaded regions indicate the total uncertainty obtained by also adding the luminosity uncertainty in quadrature. The plotted CDF Z measurement with 2.1 fb −1 [47] supersedes the earlier Z measurement with 72 pb −1 [46], but both measurements are dominated by the (common) luminosity uncertainty. The DØ experiment has not published any dedicated W and Z total cross-section measurements from Run II at the Tevatron. The DØ Z total cross section shown in figure 9(b) was obtained as part of the Z+jet measurement [48]. The CDF measurement [47] is defined as the Z/γ * → ee cross section in an invariant mass range M ee ∈ [66, 116] GeV, while the DØ measurement [48] is defined as the Z/γ * → µµ cross section in an invariant mass range M µµ ∈ [65, 115] GeV. We have therefore multiplied the CDF and DØ data by factors of 1.006 and 1.004, respectively, derived using the vrap code [49] at NNLO with MSTW08 PDFs, to correct to the Z-only cross section with M ℓℓ = M Z . We note from figure 9 that the MSTW08, ABKM09 and JR09 NNLO predictions for the W and Z total cross sections at the Tevatron are in good agreement with the CDF data [46,47], and lie around 1-σ above the DØ data [48]. In the MSTW08 fit [13], the luminosity shift for the CDF jet data was correctly tied to be the same as for the more-constraining CDF Z rapidity distribution, dσ Z /dy [47], which therefore effectively acted as a luminosity monitor. The optimal CDF normalisation in the MSTW08 NNLO fit [13] was found to be very close to the nominal value, therefore it is not surprising that the CDF Z total cross section is well described in figure 9(b). The DØ experiment instead measured the Z rapidity shape distribution, (1/σ Z )dσ Z /dy [50], also included in the MSTW08 fit, which is one reason why the DØ jet data were found to be less constraining than the CDF jet data; see ref. [32]. The optimal DØ normalisation in the MSTW08 NNLO fit [13], determined only from jet data, was around 1-σ above the nominal value, consistent with the DØ Z total cross section shown in figure 9(b). If the Tevatron jet data were normalised downwards by 20-30% (i.e. 3-5 times the luminosity uncertainty), the Tevatron W and Z total cross sections would need to normalised downwards by the same amount, resulting in complete disagreement with all theory predictions shown in fig-ure 9. This example illustrates the utility of simultaneously fitting W and Z cross sections together with jet cross sections at the Tevatron (and LHC). The luminosity shifts, common to both data sets, are effectively determined by the more precise W and Z cross sections. The luminosity uncertainty is then effectively removed from the jet cross sections, thereby allowing the jet data to provide a tighter constraint on the gluon distribution (and α S ).
To avoid these completely unrealistic luminosity shifts, r lumi. ∼ 3-5, without going into the complication of simultaneously including W and Z cross sections in the χ 2 computation, we will calculate the χ 2 values for the Tevatron jet data using eq. (3.1), but with the simple restriction that the relevant systematic shift |r lumi. | ≤ 1. More practically, this means that if |r lumi. | > 1 for any particular PDF set, we fix r lumi. at ±1 and reevaluate eq. (3.1) with the luminosity removed from the list of correlated systematics. However, we note from figure 9 that the ABKM09 predictions are slightly above the central value of the CDF W and Z data, and the HERAPDF1.0 predictions are higher by around 1-σ, while both ABKM09 and HERAPDF1.0 lie above the 1-σ limit of the DØ Z data. Allowing luminosity shifts downwards by even 1-σ is therefore distinctly generous, particularly for HERAPDF1.0, and upwards luminosity shifts would bring the ABKM09 and HERAPDF1.0 predictions into better agreement with the CDF W and Z data, and especially the DØ Z data. Therefore, it should be understood that the χ 2 values quoted in the tables we will present in section 3.2 and 3.3 are rather optimistic for ABKM09 and HERAPDF1.0, and more realistic constraints in the luminosity shifts would result in even worse χ 2 values.
The form of eq. (3.1) is slightly different from the treatment of normalisation uncertainties adopted in eq. (38) of the MSTW08 paper [13], but is the form used, for example, in the CT10 analysis [18]. Rescaling only the central value of the data in eq. (3.1), but not the uncertainties, leads to so-called "d'Agostini bias" [51,52]. However, since we are only comparing and not fitting PDFs, we use the simpler form of eq. (3.1) which has the major advantage that all shifts r k can be solved for analytically. A more sophisticated approach to the treatment of normalisation uncertainties may somewhat lessen the preference of some PDF sets for large downwards luminosity shifts, but should not affect our main conclusions. The normalisation uncertainties were treated as multiplicative rather than additive in the MSTW08 fit [13], i.e. the uncertainties were correctly rescaled to reduce bias. Moreover, large normalisation shifts for any experiment were discouraged through use of a quartic penalty term rather than the usual quadratic penalty term in eq. (3.1). These small differences in χ 2 definition mean that the MSTW08 χ 2 values we quote here will be slightly different from the values quoted in ref. [13].
Even considering the constraint on the CDF and DØ luminosities from the comparison to the weak boson cross sections (see figure 9), it might be considered that imposing |r lumi. | < 1 is too restrictive if the luminosity uncertainty is assumed to be Gaussian. However, as another reason for limiting the luminosity shifts to some extent, we note that it has been claimed (see section 6.7.4 on "Normalizations", pg. 170" in [53]) that, for many experiments, quoted normalisation uncertainties represent the limits of a box-shaped distribution rather than the standard deviation of a Gaussian distribution. This was one motivation for the more severe quartic penalty term for normalisation uncertainties in the MSTW08 analysis; see discussion in section 5.2.1 of ref. [13]. Nevertheless, if we instead impose |r lumi. | < 2 rather than |r lumi. | < 1, then the change in the χ 2 /N pts. values for the most relevant ABKM09 NNLO PDF set with µ = p T is {2. 76 [34], DØ inclusive [35], DØ dijet [36]} data, respectively, so there is not a significant improvement in the χ 2 values. However, as discussed above, our main argument does not rely on the precise form of the uncertainty on the luminosity determination, but that we can use the W and Z cross sections as a luminosity monitor, where the predictions have small theoretical uncertainties, effectively providing an accurate luminosity determination independently of the CDF and DØ values. Combining these arguments, we consider allowing luminosity shifts downwards by more than 1-σ to be excessively generous.
There is a clear trade-off between the systematic shifts r k and the parameters of the gluon distribution. Deficiencies in the theory calculation can be masked to some extent by large systematic shifts, therefore it is important to check that the optimal r k values are not unreasonable. This is straightforward when using a χ 2 definition like eq. (3.1), but is more difficult using an equivalent form written in terms of the experimental covariance matrix, Then eq. (3.1) is equivalent [41] to the more traditional χ 2 form written in terms of the inverse of the experimental covariance matrix: Npts.
as used by the ABKM and NNPDF fitting groups. More precisely, NNPDF use a refinement to treat normalisation errors as multiplicative [52], while Alekhin (ABKM) treats all correlated systematic errors as multiplicative [54,55]. It can easily be seen from eqs. (3.5) and (3.6) that treating the correlated errors as uncorrelated (V ii ′ ∝ δ ii ′ ) leads to the familiar form of where the total error is simply obtained by adding all errors in quadrature,

Inclusive jet production
In tables 1, 2 and 3 we give the χ 2 per data point, calculated using eq. (3.1) with the restriction |r lumi. | < 1, for the Tevatron Run II data on inclusive jet production [33][34][35], for different PDF sets and different scale choices  Table 1. Values of χ 2 /N pts. for the CDF Run II inclusive jet data using the k T jet algorithm [33] with N pts. = 76 and N corr. = 17, for different PDF sets and different scale choices The χ 2 values are calculated accounting for all 17 sources of correlated systematic uncertainty, using eq. (3.1), including the 5.8% normalisation uncertainty due to the luminosity determination. At most a 1-σ shift in normalisation is allowed. We highlight in bold those values lying inside the 90% C.L. region, defined by eq. (3.9), which gives χ 2 /N pts. < 0.83. The values of χ 2 /N pts. computed using eq. (3.7), simply adding all experimental uncertainties in quadrature (including luminosity), are shown in brackets in the table. If the theory prediction was identically zero, the χ 2 /N pts. values would be 25.0 (37.5) with (without) accounting for correlations between systematic uncertainties. 100 replica sets. We give the χ 2 /N pts. values defined by simply adding all uncertainties in quadrature, eq. (3.7), in brackets in the tables. In this case many PDF sets and scale choices give a χ 2 /N pts. ≪ 1, so the consistent treatment of correlated uncertainties is vital for the jet data to discriminate. In the table captions we give the χ 2 values with an identically zero theory prediction, T i ≡ 0, just to illustrate how the correlated systematic shifts can partially accommodate a clearly inadequate theory prediction. We highlight in bold the χ 2 values lying inside the 90% C.L. region defined as where ξ 50 and ξ 90 are the 50th and 90th percentiles of the χ 2 -distribution with N pts. degrees of freedom. (These quantities are defined in detail in section 6.2 of ref. [13].) Here, χ 2 0 is defined as the lowest χ 2 value of all theory predictions in each table, i.e. assumed to be close to the best possible fit, so that the rescaling factor χ 2 0 /ξ 50 in eq. (3.9) empirically accounts for  any unusual fluctuations preventing the best possible fit having χ 2 ≃ ξ 50 ≃ N pts. [41]. The 90% C.L. region given in this way is used to determine the PDF uncertainties according to the "dynamical tolerance" prescription introduced in ref. [13], so PDF sets with χ 2 values far outside this region cannot be considered to give an acceptable description of the data. We consider NLO PDFs from MRST04 [56], MSTW08 [13], CTEQ6.6 [17], CT10 [18], NNPDF2.1 [19], HERAPDF1.0 [23], HERAPDF1.5 (preliminary) [57], ABKM09 [15] and GJR08 [58,59]. We consider NNLO PDFs from MRST06 [60], MSTW08 [13], HER-APDF1.0 [23], ABKM09 [15] and JR09 [21,22]. The MRST04 and MRST06 fits only included Tevatron Run I data [61,62], and were superseded by the MSTW08 fits, but we show the χ 2 values here just to demonstrate that these older fits do not give a good description of the newer Tevatron Run II data due to their harder high-x gluon distribution. The CTEQ6.6 fit includes only the Tevatron Run I data [61,62], while the CT10 fit includes Run II data [34,35] in addition to the Run I data [61,62], contrary to the MSTW08 and NNPDF2.1 fits which include only Run II data [33,36]. The GJR08 fit included some Run NLO PDF (with NLOσ)  Table 3. Values of χ 2 /N pts. for the DØ Run II inclusive jet data using a cone jet algorithm [35] with N pts. = 110 and N corr. = 23, for different PDF sets and different scale choices The χ 2 values are calculated accounting for all 23 sources of correlated systematic uncertainty, using eq. (3.1), including the 6.1% normalisation uncertainty due to the luminosity determination. At most a 1-σ shift in normalisation is allowed. We highlight in bold those values lying inside the 90% C.L. region, defined by eq. (3.9), which gives χ 2 /N pts. I [62] and Run II [63] data, while the JR09, ABKM09 and HERAPDF fits did not include any Tevatron jet data. The most constraining data set appears to be the CDF Run II inclusive jet data using the k T jet algorithm [33] (see table 1) where, other than MSTW08, only NNPDF2.1 gives an acceptable description for µ = p T , while HERAPDF1.0 and ABKM09 typically give χ 2 /N pts. ∼ 2-3, and CTEQ6.6/CT10 give better values but still much worse than MSTW08 (and NNPDF2.1). The GJR08/JR09 sets and the HERAPDF1.0 NNLO set with α S (M 2 Z ) = 0.1176 give a reasonable description, at a similar level to CT10, and give predictions for gg → H cross sections at the Tevatron which are much closer to the MSTW08 predictions than those from ABKM09 and the HERAPDF1.0 NNLO set with α S (M 2 Z ) = 0.1145. The same trend is apparent, but to a somewhat lesser extent, for the CDF Run II inclusive jet data using the cone-based Midpoint jet algorithm [34] (see table 2) and the DØ Run II inclusive jet data using a cone jet algorithm [35] (see table 3).
In figures 10, 11 and 12 we compare the description of the Tevatron inclusive jet data by  Figure 10. Data/theory ratios for the CDF Run II inclusive jet data using the k T jet algorithm [33] with N pts. = 76 and N corr.  Figure 11. Data/theory ratios for the CDF Run II inclusive jet data using the cone-based Midpoint jet algorithm [34] with N pts. = 72 and N corr. = 25, for MSTW08 and ABKM09 NNLO PDFs with NLO partonic cross sections supplemented by 2-loop threshold corrections, with scale choice µ R = µ F = p T , and (a) all experimental errors added in quadrature, then (b) accounting for correlated systematic uncertainties using eq.  Figure 13. Data/theory ratios for the DØ Run II dijet data using a cone jet algorithm [36] with N pts. = 71 and N corr. = 70, for MSTW08 and ABKM09 NNLO PDFs and NLO partonic cross sections, with scale choice µ R = µ F = p T , where p T ≡ (p T 1 +p T 2 )/2, with (a) all experimental errors added in quadrature, then (b) accounting for correlated systematic uncertainties using eq. (3.1) and showing only the uncorrelated experimental errors. the MSTW08 and ABKM09 NNLO PDFs (recall that the latter give the lowest predictions for Tevatron Higgs cross sections) by showing the ratio of data to theory defined in two different ways: (a) first we use the original data points D i /T i with uncertainties given by adding all errors in quadrature (including luminosity), σ tot.

CDF Run II inclusive jet data (Midpoint) (data after systematic shifts, show uncorrelated errors)
i /T i , with the appropriate χ 2 value in the plot legends obtained using eq. (3.7), then (b) we use the shifted data pointsD i /T i with uncertainties given by σ uncorr.
i /T i , with the χ 2 calculated according to eq. (3.1) and showing the two terms separately in the plot legends. The p T values for ABKM09 are slightly offset for clarity in the plots. The size of the second penalty term in eq. (3.1) is some measure of how much the data points are shifted compared to their systematic errors. For example, if the penalty term Ncorr. k=1 r 2 k > N corr. , then the data points are shifted by, on average, more then 1-σ for each systematic source k. In general, a poor description of data before the systematic shifts leads to a large penalty term and a poor description also after the systematic shifts, although this general statement is not universally true. We note that the shape of the data/theory ratio, both before and after the systematic shifts, looks remarkably similar as a function of both transverse momentum p T and rapidity y in figures 10 and 11. This demonstrates very clearly that the two CDF inclusive jet measurements [33,34] each contain the same data, but simply analysed in a different way, and the change in analysis method is accounted for extremely well by the change in the theory. Hence, it is not at all surprising that the two data sets can be well described by the same PDF set. Indeed, it was explicitly demonstrated in ref. [34] that the ratios of the cross sections measured with the two jet algorithms were in reasonable agreement with theoretical expectations.

Dijet production
In table 4 and figure 13 we show similar results for the DØ Run II dijet data [36], measured as a function of the dijet invariant mass, M JJ , and the largest absolute rapidity, |y| max , of the two jets with the largest transverse momentum. Again, the M JJ values for ABKM09 are slightly offset for clarity in figure 13. The fastnlo grids are provided with a scale choice proportional to the mean transverse momentum of these two jets, p T ≡ (p T 1 + p T 2 )/2, and we show results with µ R = µ F = µ = {p T /2, p T , 2p T } in table 4. Taking µ = p T /4 leads to negative cross sections at large M JJ and large |y| max . We multiply the fastnlo predictions by a factor 4 to account for a mismatch in the bin width factors of the provided grids. There are no 2-loop threshold corrections available, so we are forced to use only the pure NLO partonic cross sections with the NNLO PDFs. It can be seen that the trend in the χ 2 values for the dijet data shown in table 4 appears to be rather different from the inclusive jet data shown in tables 1, 2 and 3. In particular, in contrast to the case for inclusive jets, the ABKM09 set gives the best description for µ = p T , whereas MSTW08 and NNPDF2.1 have χ 2 /N pts. ∼ 2 and CTEQ6.6/CT10 has χ 2 /N pts. ∼ 4-5. For µ = 2p T there is a significant improvement in χ 2 for MSTW08 and NNPDF2.1, and MSTW08 NNLO for µ = 2p T gives the best description out of all PDF sets and scale choices, while the CTEQ6.6/CT10 sets still have χ 2 /N pts. ∼ 3 even for the larger scale choice. However, it is interesting to note that while figures 10(a) and 11(a) show a very similar trend for the data/theory ratios, figures 12(a) and 13(a) show quite a different trend, NLO PDF (with NLOσ)  Table 4. Values of χ 2 /N pts. for the DØ dijet data using a cone jet algorithm [36] with N pts. = 71 and N corr. = 70, for different NLO PDF sets and different scale choices implying that the change in theory in using the NLO dijet cross section at the same scale as the inclusive jet cross section does not account for the difference in the data produced by the two methods [35,36] of binning and analysis. At LO we have M JJ = 2p T cosh y * where y * = |y 1 − y 2 |/2, with y 1,2 the rapidities of the two jets. It is clear that p T is a better measure of the "hardness" of the process than M JJ and therefore µ = p T is the most common scale choice for dijet production. (Consider, for example, the extreme case of elastic pp scattering where each final-state proton is considered to be a "jet", then M JJ ≈ √ s, but p T ≈ 0.) More generally, typical scale choices in fixed-order perturbative QCD calculations are usually, for example, the mass or transverse momentum of a produced particle, or a scalar sum of such scales added either linearly or in quadrature. However, it is clear that choices of scale involving both p T and y * are perfectly feasible for dijets, whereas some multiple of p T seems more obviously the scale choice for inclusive jets. There is no reason that the choice which best mimics the full calculation at fixed order for inclusive jets need be the same as for dijets binned in M JJ , i.e. the structure of higher order corrections is not automatically the same. Indeed, a hybrid scale choice was proposed in ref. [64] to interpolate between a scale choice based on p T and one based on M JJ , namely µ = AM JJ /(2 cosh By * ), with the two adjustable parameters chosen to be A = 0.5 and B = 0.7 so that the difference between the O(α 3 S ) calculation and the Born calculation was small over the angular region of interest [64]. It would be interesting to investigate whether such a scale choice could resolve the somewhat different conclusions reached from the Tevatron Run II inclusive and dijet data. There is no requirement that the scale choice for dijets be the same as for inclusive jets. Taking µ = p T for inclusive jets and µ = 2p T = p T 1 + p T 2 for dijets, then the MSTW08 (and NNPDF2.1) PDFs would give a good description of all four Tevatron Run II jet data sets [33][34][35][36].
Another difference, possibly correlated to the issue of scale choice, is that the dijet data may probe higher x values than the inclusive jet data. If there are two jets labelled "1" and "2", and jet "1" has high p T in the forward region, then the phase space for the jet "2" is integrated over in the inclusive jet cross section, but will typically lie in the central region, creating an imbalance in the x values of the two initial partons. On the other hand, for the dijet cross section at high M JJ values, if jet "1" lies in the forward region, then jet "2" will typically lie at the same absolute rapidity in the opposite direction, giving similarly large x values of the two initial partons. Since high-x PDFs evolve very quickly, probing two high-x PDFs increases sensitivity to (factorisation) scale choices. This sensitivity will be most extreme when both PDFs are evolving quickly in the same direction (for example, both getting smaller with increasing scale), rather than one PDF getting smaller and one PDF getting larger as would be the case with one high-x parton and one low-x parton. This effect automatically means that the higher-order corrections must be slightly different in the two cases of inclusive jet and dijet production.

Scale dependence of jet cross sections
In figure 14 we compare the K-factors for the DØ inclusive and dijet data, defined as the ratio of the NLO (both with/without the 2-loop threshold corrections) jet cross sections to the LO jet cross sections, computed with the same MSTW08 NNLO PDFs (and α S ) in the numerator and denominator of the ratio. Using another PDF choice, such as ABKM09 NNLO, makes little difference to the K-factors. The choice µ = p T /2 has historically been favoured in MRST/CTEQ fits because the K-factor is close to 1 at central rapidity. However, going to forward rapidities with the choice µ = p T /2, the K-factor decreases substantially with increasing p T . The K-factor with the choice µ = p T is more uniform (with moderate size) across all rapidity bins and p T values, hence µ = p T was chosen for the MSTW08 analysis [13]. It is striking, however, that although the NLO corrections are ∼ 60% for µ = 2p T , and a further 20% or more with the 2-loop threshold corrections, the shape of the K-factor is rather more stable across all rapidity bins and p T with this choice. In figure 15 we show the ratio of the NLO (both with/without the 2-loop threshold corrections) jet cross sections with different scale choices to the NLO jet cross section with µ R = µ F = p T , again computed with the same MSTW08 NNLO PDFs (and α S ) in the numerator and denominator of the ratio. It can be seen that the use of the 2-loop threshold corrections for the inclusive jet cross sections stabilises the scale dependence (except at the very highest rapidity and p T values where the low scale choice still leads to a large variation). To some extent, different scale choices will be compensated by different systematic shifts, particularly for the luminosity (see appendix A). The predictions for µ = p T are generally in the middle of the other two choices, but this breaks down at high rapidity and p T values. Indeed, for dijets µ = p T ceases to be the central prediction at nearly all p T in the two highest rapidity bins, and is progressively less so in the middle rapidity bins than for the case of inclusive jets. This supports the idea that the optimal     i , for each of the four Tevatron data sets on jet production, with theory predictions calculated using either MSTW08 or ABKM09 NNLO PDFs, compared to the expectation of a Gaussian distribution with unit width. choice for dijets might be p T multiplied by a function f (y * ), growing with increasing y * , so that µ = p T · f (y * ) would be the central prediction over all |y| max bins. (The y * variable is closely related to the |y| max variable used by the DO dijet data [36].) In the absence of readily-available theory predictions for such a scale choice, the best description of dijet data by PDFs obtained from fitting to inclusive jet data seems to be given, as a compromise, by a scale of µ = c p T with c > 1, with our specific example being c = 2.

Distributions of pulls and systematic shifts
In figure 16 we show the distributions of pulls, (D i − T i )/σ uncorr.
i , for all four Tevatron Run II data sets on jet production, with theory predictions calculated using either MSTW08 or ABKM09 NNLO PDFs and a scale choice µ R = µ F = µ = p T . We show the expected behaviour of a Gaussian distribution with unit width, and the first χ 2 term in eq. (3.1) given simply by the sum of pulls over all data points. The histogram error bars are simply given by the square root of the number of entries. We see that the distribution of pulls is fairly close to the expected Gaussian behaviour for all four data sets, although the tails for the inclusive jet data with ABKM09 are somewhat broader than expected, leading to larger χ 2 contributions than for MSTW08, particularly for the CDF data using the k T  Figure 17. Distributions of the systematic shifts, r k given by eq. (3.3), for each of the four Tevatron data sets on jet production, with theory predictions calculated using either MSTW08 or ABKM09 NNLO PDFs, compared to the expectation of a Gaussian distribution with unit width. jet algorithm [33] shown in figure 16(a). However, it is clear that this source does not account for the complete differences in χ 2 seen previously. In figure 17 we show the similar distributions of the systematic shifts, r k , again for all four Tevatron Run II data sets on jet production. We show the expected behaviour of a Gaussian distribution with unit width and the penalty χ 2 term simply given by the sum of the r 2 k values. For the inclusive jet data, the systematic shifts for MSTW08 show the expected Gaussian behaviour, with small penalty terms Ncorr. k=1 r 2 k < N corr. . On the other hand, the systematic shifts for ABKM09 deviate substantially from Gaussian behaviour, with much larger penalty terms, in particular for the CDF inclusive jet data using the k T algorithm shown in figure 17(a). The systematic shifts for the dijet data shown in figure 17(d) have a much narrower distribution than the expected Gaussian behaviour for both MSTW08 and ABKM09, suggesting that the systematic errors are overestimated, are non-Gaussian, or are not independent (or a combination of these three explanations). Note that the number of systematic sources (N corr. = 70) for the dijet data is much greater than for any of the inclusive jet data sets. Indeed, this allows the value of χ 2 for the description of data by an identically zero theory prediction to be lower than for some of the PDF sets; see table 4.
The presentation of the results in figures 16 and 17 enables a separation between contributions to the χ 2 definition, eq. (3.1), from uncorrelated and correlated errors, respectively. This allows a more informed assessment of the fit quality compared to the more traditional χ 2 definition of eq. (3.6) in terms of the experimental covariance matrix, used by the ABKM and NNPDF fitting groups; see also appendix B.3 of ref. [42] and section 4 of ref. [18].

Other jet cross sections from collider experiments
The DØ Collaboration has recently made a measurement [65] of the three-jet differential cross section as a function of the invariant mass of the three jets with the largest transverse momentum in an event. An exercise has been carried out similar to the one presented here, where the χ 2 has been evaluated for different PDF (and α S ) choices and scale choices The trend is that MSTW08 and NNPDF2.1 are favoured, as for the inclusive jet study presented here, while ABKM09 is worse, and CT10 and HERAPDF1.0 are still poorer. We have followed a similar approach to that of ref. [65] by evaluating the χ 2 only for the central PDF fit, without accounting for PDF uncertainties. Since the Tevatron jet data provide by far the most direct constraint on the high-x gluon distribution, the agreement of the central PDF fit is more important and relevant than obtaining agreement only within possibly large PDF uncertainties. However, the potential choice of scales for the three-jet cross section is even broader than for the dijet cross section. The LHC data on jet production [66][67][68][69] are becoming more precise and show some sensitivity to the PDF choice. However, these data are still being understood and are not presented with separated correlated systematic uncertainties which would allow a quantitative χ 2 comparison. Moreover, the general sensitivity is to lower x T ∼ 2p T / √ s, and so less relevant for Higgs production at the Tevatron. Isolated photon production at the LHC may also provide a direct constraint on the gluon distribution [70]. The HERA jet data are less sensitive to the gluon distribution at high x values, being more of a constraint for x ∼ 0.001-0.1, and there is no NNLO calculation, or any approximation such as the 2-loop threshold corrections available for the Tevatron inclusive jet data.

Summary
Comparison with Tevatron jet data is subtle because of the large correlated systematic uncertainties. The systematic shifts, eq. (3.2), can compensate for inadequacies in the theory calculation. The traditional χ 2 definition in terms of the experimental covariance matrix, eq. (3.6), can hide such systematic shifts. In particular, we find that the Tevatron jet data need to be normalised downwards by typically between 3-σ and 5-σ (see appendix A) to achieve the best agreement with some PDF sets, particularly the ABKM09 predictions. Even if the luminosity shift is artificially constrained, the other systematic shifts move by large amounts for the inclusive jet data, incompatible with the Gaussian expectation. No such problems are observed for the MSTW08 predictions. It can also be seen from the plots in ref. [12] that the unshifted Tevatron jet data lie significantly above the theory predictions even after including these data in variants of the ABKM09 fit. Constraining the Tevatron luminosity shifts, for example, so that the predicted W and Z cross sections agreed with Tevatron data, would increase the constraining power of the Tevatron jet data and thereby very likely give a larger α S and high-x gluon distribution than the current studies of Alekhin, Blümlein and Moch (ABM) [12]. Even with the existing treatment, the NNLO Tevatron gg → H cross section for M H = 165 GeV goes up by {15, 12, 17, 11}% when including the {CDF k T [33], CDF Midpoint [34], DØ inclusive [35], DØ dijet [36]} data set in variants of the ABKM09 fit [12]. The dijet data has a potentially wider range of allowed scale choices than the inclusive jet data. We conclude that the data on inclusive jet production therefore provide the cleanest probe of different PDF sets.

Value of strong coupling α S from DIS
There is a common lore (see, for example, ref. [71]) that DIS-only fits prefer low α S (M 2 Z ) values, but ref. [32] showed that not all DIS data sets prefer low α S (M 2 Z ) values. In particular, this was found to be true only for BCDMS data, and for E665 and SLAC ep data, while NMC, SLAC ed and HERA data preferred high α S (M 2 Z ) values within the context of the global fit [32]. (See also the recent NNPDF study at NLO using an "unbiased" PDF parameterisation [72].) It is well known that α S is highly anticorrelated with the low-x gluon distribution through scaling violations of HERA data: ∂F 2 /∂ ln(Q 2 ) ∼ α S g. Then α S is correlated with the high-x gluon distribution through the momentum sum rule; see, for example, figure 14(b) of ref. [32]. Restrictive gluon parameterisations, without the negative small-x term allowed by MSTW [13], can therefore bias the extracted α S value. For example, the default MSTW08 NNLO fit obtained α S (M 2 Z ) = 0.1171 ± 0.0014, while imposing the restriction of a positive input gluon at Q 2 0 = 1 GeV 2 gave a best-fit α S (M 2 Z ) = 0.1157, but with a χ 2 worse by 63 units for the global fit to 2615 data points [32]. 5 What is α S from only DIS data in the MSTW08 NNLO fit? 6 Recall that the global fit gave α S (M 2 Z ) = 0.1171 ± 0.0014 [32]. To expand on the studies made in ref. [32], we performed a new NNLO DIS-only fit, which gave a best-fit α S (M 2 Z ) = 0.1104, but with an input gluon distribution which went negative for x > 0.4 due to lack of any data constraint. This implies a negative charm structure function, F charm 2 , and a terrible description (χ 2 /N pts. ∼ 10 including correlated systematic errors) of Tevatron jet data using the obtained PDFs. A DIS-only fit fixing the high-x gluon parameters to prevent such bad behaviour gave α S (M 2 Z ) = 0.1172, i.e. very similar to the global fit. However, a NNLO fit which imposed the condition of the positive low-x gluon, which stopped the gluon from going negative at high x values, and which also omitted the Tevatron jet data, gave α S (M 2 Z ) = 0.1139, rather closer to the ABKM09 value. The very low value of α S (M 2 Z ) = 0.1104 found in the DIS-only fit is due to the dominance of BCDMS data. We can show this explicitly by removing the BCDMS data from the DIS-only fit, then the best-fit α S (M 2 Z ) moves from 0.1104 to 0.1193. Repeating the global fit with BCDMS data removed gives α S (M 2 Z ) = 0.1181, i.e. a change by less than the quoted experimental uncertainty of ±0.0014. The conclusion is that the Tevatron jet data are vital to pin down the high-x gluon, giving a smaller low-x gluon and therefore a larger α S in the global fit compared to a DIS-only fit, at the expense of some deterioration in the fit quality of the BCDMS data. 7 The benefits of including the Tevatron jet data to obtain sensible results in a simultaneous fit of PDFs and α S therefore greatly outweighs any disadvantage such as lack of complete NNLO corrections.
The only input DIS value to the current world average α S (M 2 Z ) [31] is the BBG06 value [74], which is from a non-singlet analysis and therefore in principle free of assumptions made about the gluon distribution. A value of However, using the MSTW08 NNLO central fit, contributions other than valence quarks are found to make up about 10% (2%) of F p 2 at x = 0.3 (x = 0.5). As an exercise we performed the MSTW08 NNLO DISonly fit just to F p 2 and F d 2 for x > 0.3 (comprising 282 data points, 160 of these from BCDMS), which gave α S (M 2 Z ) = 0.1103 (0.1130) without (with) the singlet contribution included. This is even lower than the BBG06 value presumably due to lack of the y > 0.3 cut on BCDMS data applied in the BBG06 analysis. The low value of α S (M 2 Z ) found by BBG06 [74] is therefore due to both dominance of BCDMS data and by what we conclude is the unjustified neglect of the singlet contribution to F p 2 and F d 2 for x ≥ 0.3. Given that it was argued above that the Tevatron jet data are needed to pin down the high-x gluon, we conclude that an extraction of α S (M 2 Z ) only from inclusive DIS data is not meaningful, and the closest possible to a reliable extraction is the MSTW08 NNLO combined analysis of DIS, Drell-Yan and jet data [13,32]: This value is the only NNLO determination, from a simultaneous fit with PDFs, which is in agreement with the current world average α S (M 2 Z ) = 0.1184 ± 0.0007 [31]; see figure 5(b).

Treatment of NMC data and stability to low Q 2 data
A recent claim has been made [11] that the bulk of the MSTW08/ABKM09 difference in both the extracted α S (M 2 Z ) value and the gg → H predictions is explained by the treatment of NMC data [14]. The differential cross section for DIS of charged leptons off nucleons, ℓN → ℓX, neglecting the nucleon and lepton masses, and assuming single-photon exchange, is The low y data points from BCDMS are strongly affected by the energy scale uncertainty of the scattered muon. It has been advocated to impose a cut of y > 0.3 on the BCDMS data, which caused αS(M 2 Z ) to increase by about 0.004 in a fit to only BCDMS data and by about 0.002 in a combined fit to H1 and BCDMS data [73].  Figure 18. [14], the Q 2 -dependent SLAC R 1990 parameterisation [75], and the MSTW08 NNLO calculation including 1-σ PDF uncertainties [13]. (b) F 2 versus Q 2 (in units of GeV 2 ) at x = 0.025 comparing the two NMC extractions [14] using either R NMC or R 1990 .
is the ratio of the γ * N cross sections for longitudinally and transversely polarised photons, Q 2 is the photon virtuality, x is the Bjorken variable and y ≃ Q 2 /(x s) is the inelasticity (with √ s the ℓN centre-of-mass energy). The ABKM09 [15] analysis fitted the NMC differential cross sections directly, calculating F L to O(α 2 S ) and including empirical higher-twist corrections. The MSTW08 [13] analysis instead fitted the NMC F 2 values corrected for R, where [14] R(x, Here, R NMC (x) was a (Q 2 -independent) value extracted from NMC data, while R 1990 (x, Q 2 ) was a Q 2 -dependent empirical parameterisation of SLAC data dating from 1990 [75]. By replacing the NMC differential cross-section data by NMC F 2 data, ABM [11] find that their best-fit α S (M 2 Z ) moves from 0.1135 to 0.1170 and their gg → H cross sections at the Tevatron and LHC move closer to the MSTW08 values. ABM [11] therefore conclude that the use of NMC F 2 data in the MSTW08 fit rather than the differential cross section is the main reason for the higher α S (M 2 Z ) and Higgs cross sections obtained with MSTW08. We agree that it is more consistent to fit directly to the NMC differential cross-section data, so here we respond to this rather dramatic assertion made by ABM [11], which would obviously be very worrying if correct. However, rather than repeat the MSTW08 analysis by fitting the NMC differential cross sections, we note that the original NMC paper [14] made an alternative extraction of F 2 values using the SLAC R 1990 parameterisation [75]. In figure 18(a) we compare R NMC with R 1990 in the most affected bin of x = 0.025, i.e. a low x value where there are a reasonable number (7) of NMC data points surviving the cut on Q 2 ≥ 2 GeV 2 and where the difference between R NMC and R 1990 is at its largest. Recall that a low x value means a high y value and from eq. (5.1) the correction term from R is only important at large y. In figure 18(a) we also show the MSTW08 NNLO prediction, including PDF uncertainties at 68% C.L., with F L calculated to O(α 3 S ) and without  Table 5. Effect of NMC treatment on α S (M 2 Z ) and Higgs cross sections (M H = 165 GeV). We also show the effect of raising the cuts imposed on the DIS data compared to the default of removing data with Q 2 < 2 GeV 2 and W 2 < 15 GeV 2 . Finally, we show the effect of simply fixing α S (M 2 Z ) to be close to the ABKM09 value, or performing a fit with a positive-definite input gluon distribution and no jet data, and we compare directly to ABKM09. any higher-twist corrections. We see that it gives a good description of the SLAC R 1990 parameterisation, with any differences being very much smaller than those between R NMC and R 1990 . We note that NMC/BCDMS/SLAC F L data are included in the MSTW08 fit and are well-described at NNLO but less well at NLO (see figure 5 of ref. [32]), so the O(α 3 S ) coefficient functions are needed for a good description and the larger MSTW08 α S (M 2 Z ) perhaps explains why there is less room for higher-twist corrections, contrary to the findings of the ABM analysis. Nevertheless, figure 18(a) demonstrates that fitting the alternative NMC F 2 data extracted using the SLAC R 1990 parameterisation will give very similar results to fitting the NMC differential cross sections. In fact, given that R 1990 in figure 18(a) generally has a slightly steeper Q 2 dependence than the MSTW08 parameterisation, using this will slightly overestimate the true impact of fitting the NMC differential cross sections. In figure 18(b) we compare the two different NMC F 2 extractions, again for the most affected bin of x = 0.025, and we see that there is little difference, certainly nothing that seems likely to change α S (M 2 Z ) by 0.0035 in a fit where it is constrained with an uncertainty of about 0.0014 by over 2000 other data points.
In table 5 we show the effect of repeating the MSTW08 NNLO fit with the NMC F 2 data extracted using R 1990 on α S (M 2 Z ) and the Higgs cross sections (for M H = 165 GeV) at the Tevatron and LHC, and in figure 19 we show the change in the gluon distribution at the corresponding scale. We make other fits either cutting the NMC F 2 data for x < 0.1, above which the R correction in eq. (5.1) is very small indeed, or completely removing all NMC F 2 data. In all cases there is very little change in α S (M 2 Z ), the gluon distribution, and the Higgs cross section. We conclude that the treatment of NMC data cannot explain the difference between the MSTW08 and ABKM09 results. Similar stability has been found by the NNPDF group [76], but in a less relevant study at NLO with fixed α S .
The cuts on DIS data are not explicitly given in the ABKM09 paper [15], but the previous AMP06 paper [77] mentions that DIS data are removed with Q 2 < 2.5 GeV 2 and We also show the effect of raising the cuts imposed on the DIS data compared to the default of removing data with Q 2 < 2 GeV 2 and W 2 < 15 GeV 2 . Finally, we show the effect of simply fixing α S (M 2 Z ) to be close to the ABKM09 value, or performing a fit with a positive-definite input gluon distribution and no jet data, and we compare directly to ABKM09. W 2 < (1.8 GeV) 2 = 3.24 GeV 2 , compared to the MSTW08 fit which removes DIS data with Q 2 < 2 GeV 2 and W 2 < 15 GeV 2 . The much weaker cut on the hadronic invariant mass (squared), W 2 ≃ Q 2 (1/x − 1), clearly explains why higher-twist corrections are more important in the ABKM09 analysis. To investigate the possible effect of neglected highertwist corrections on the MSTW08 NNLO fit we raised the cuts to remove DIS data with W 2 < 20 GeV 2 and either Q 2 < 5 GeV 2 or Q 2 < 10 GeV 2 . The results are shown in table 5 and figure 19. The changes in α S , the gluon distribution and the Higgs cross sections are generally small and within uncertainties, although with the strongest Q 2 cut there is no data constraint below x = 10 −4 and little just above, so the PDFs differ but have large uncertainties at low x values. 8 In table 5 and figure 19 we show the results of the MSTW08 NNLO fit with a fixed 8 We also investigated the effect of increasing the cuts on W 2 and Q 2 in variants of the MSTW NLO fit. The changes were slightly bigger, with αS(M 2 Z ) changing from 0.1202 to 0.1192 and 0.1175 with Q 2 cuts of 5 and 10 GeV 2 , respectively. Similarly, the changes in PDFs and cross-section predictions are generally slightly greater at NLO than at NNLO, i.e. as expected there is some improved stability at higher orders. α S (M 2 Z ) = 0.113 [32] (slightly below the ABKM09 value), and even in this case the gluon distribution and Higgs cross sections move only part of the way towards the ABKM09 result, as already seen in figure 7. The MSTW08 input gluon parameterisation is [13] xg(x, compared to the much more restrictive functional forms of the other NNLO fits, namely: JR09 [21]: xg(x, Q 2 0 = 0.55 HERAPDF1.0 [23]: xg(x, The normalisation A g is determined from the momentum sum rule constraint, leaving 7 free parameters for MSTW08 compared to only 3 for ABKM09 and only 2 for JR09 and HERAPDF1.0 (although the value of Q 2 0 is optimised in the case of JR09). In the lack of any direct data constraint on the high-x gluon distribution, the other fits are therefore constrained by the form of the input parameterisation, avoiding the pathological behaviour of the negative high-x gluon distribution seen for the MSTW08 NNLO DIS-only fit described in section 4. As already mentioned in that section, in an attempt to mimic the ABKM09 fit we performed a variant of the MSTW08 NNLO fit without jet data and with the second term of eq. (5.3) set to zero. The ǫ g and γ g parameters were fixed in the fit iteration before the high-x gluon distribution went negative. The results of this fit are shown in table 5 and figure 19 and it goes some way towards reproducing the high-x gluon of the ABKM09 fit and the corresponding Tevatron gg → H prediction, certainly closer than we come with other modifications. Finally, we then investigated the effect of using NMC data corrected using R 1990 rather than R NMC in this fit. Similar to our default fit all changes were at the percent level, or less, so we do not explicitly show them, although the gluon does move marginally closer again to that of ABKM09.
Other differences between the two analyses are that ABKM09 used the NMC data for separate muon beam energies, whereas MSTW08 used the NMC data averaged over beam energies, which reduces the maximum effect of the change in R for a particular data point, i.e. at a given x and Q 2 , a data point at high y, and so very sensitive to R at a low beam energy, is at lower y for a higher beam energy. In the case of the averaged NMC data, correlated systematic uncertainties are unavailable, so the MSTW08 fit simply added errors (other than normalisation) in quadrature similar to the simple χ 2 form of eq. (3.7). As with the Tevatron jet data, deficiencies in the theory calculation may be hidden, without much trace, by large systematic shifts implicit in the χ 2 definition, eq. (3.6), similar to that used in the ABKM09 analysis. We conclude that the greater sensitivity to the treatment of NMC data found by ABM [11] is due to a variety of reasons, but perhaps most significantly, the inclusion of higher-twist corrections due to the weaker cuts on DIS data, and, as we have repeatedly emphasised, the lack of additional constraints provided by the Tevatron jet data to pin down the high-x gluon distribution.

Conclusions
The anomalously large higher-order QCD corrections to Higgs production at the Tevatron and LHC, via the dominant production channel of gluon-gluon fusion through a top-quark loop, mandate the use of (at least) NNLO calculations, together with corresponding NNLO PDFs and α S values. The Tevatron Higgs cross section, in particular, requires knowledge of the gluon distribution at large x 0.1 where constraints from DIS or Drell-Yan data are weak and the only direct constraint comes from Tevatron inclusive jet production. The MSTW08 fit [13] is currently the only public NNLO PDF set including the Tevatron jet data, and is used in the analyses of the Tevatron [4,5] and LHC [3] experiments, while other NNLO PDF fitting groups (ABKM09 [15], JR09 [21,22], HERAPDF1.0 [23]) choose to omit it, finding quite different results for the predicted Higgs cross sections. This common choice to use only the MSTW08 set, and not the other publicly available NNLO PDF sets, has faced a barrage of recent criticism [7][8][9][10][11][12], which we have responded to in detail in this paper. We summarise our main findings below: • We do not recommend that the (experimental) PDF+α S uncertainty be supplemented with an additional theoretical uncertainty on α S when calculating uncertainties on predicted cross sections, contrary to the approach taken in refs. [7,8].
• The claim [9] that the HERAPDF1.0 NNLO set with α S (M 2 Z ) = 0.1145 lowers the Higgs cross section compared to MSTW08 by ≈ 40% for M H ≈ 160 GeV at the Tevatron is due to a mistake in the calculation, and therefore the conclusions in the published version of ref. [9] are flawed. On the other hand, the observed 25% reduction with the central value of ABKM09 is still a serious problem and we give evidence in this paper that the ABKM09 set is not consistent enough with existing Tevatron data to be used for the calculation of Higgs cross sections.
• Comparison with Tevatron jet data is subtle because of the large correlated systematic uncertainties and the need to make choices in luminosity which are consistent with the predictions for W and Z cross sections. The traditional χ 2 definition in terms of the experimental covariance matrix, eq. (3.6), can hide large systematic shifts, which can compensate for inadequacies in the theory calculation. In particular, we find that the Tevatron jet data need to be normalised downwards by typically between 3-σ and 5-σ to achieve the best agreement with the ABKM09 (and some HERAPDF) predictions; see appendix A. Even if the luminosity shift is artificially constrained, the other systematic shifts move by large amounts for the inclusive jet data, incompatible with the Gaussian expectation. No such problems are observed for the MSTW08 predictions and good agreement is found with all Run II inclusive jet data, and also with the dijet data if taking a larger scale choice than for the inclusive jet data.
• We have demonstrated that the MSTW08 fit is stable to the treatment of NMC F 2 data, unlike the ABKM09 fit [11], most likely because of the averaging over muon beam energies, because the Tevatron jet data pin down the high-x gluon distribution, and also due to the stronger cuts reducing the need for large higher-twist corrections.  Table 6. Values of χ 2 /N pts. for the CDF Run II inclusive jet data using the k T jet algorithm [33] with N pts. = 76 and N corr. Moreover, the MSTW08 NNLO determination of the strong coupling α S is compatible with the world average value, unlike other NNLO determinations shown in figure 5(b).
We conclude that the current Tevatron Higgs exclusion bounds [4,5] are robust, at least with respect to the treatment of PDFs and α S in the calculation of the Higgs cross section. Similar remarks hold for the Higgs cross sections at the LHC recently calculated in ref. [3].
A Appendix: χ 2 tables with unrestricted luminosity shifts For completeness, in tables 6, 7, 8 and 9 we show χ 2 /N pts. values without the restriction in the luminosity shifts of |r lumi. | ≤ 1 imposed in the main tables given in section 3. Recall from eq. (3.2) that a positive value of r lumi. means a downwards shift in the luminosity, so we choose to give in brackets the values of "−r lumi. ", i.e. negative numbers correspond to downwards shifts in the luminosity. In the table captions we give the χ 2 values with an identically zero theory prediction (T i ≡ 0) just to illustrate an extreme case of how large downwards luminosity shifts can partially accommodate an inadequate theory prediction.  Table 7. Values of χ 2 /N pts. for the CDF Run II inclusive jet data using the cone-based Midpoint jet algorithm [34] with N pts. = 72 and N corr.   Table 8. Values of χ 2 /N pts. for the DØ Run II inclusive jet data using a cone jet algorithm [35] with N pts. = 110 and N corr.   Table 9. Values of χ 2 /N pts. for the DØ dijet data using a cone jet algorithm [36] with N pts. = 71 and