IPPP/02/49
DCPT/02/98
CavendishHEP2002/10
5 November 2002
Uncertainties of predictions from parton distributions
I: Experimental errors
A.D. Martin, R.G. Roberts, W.J. Stirling and R.S. Thorne^{1}^{1}1Royal Society University Research Fellow.
Institute for Particle Physics Phenomenology, University of Durham, DH1 3LE, UK
Rutherford Appleton Laboratory, Chilton, Didcot, Oxon, OX11 0QX, UK
Cavendish Laboratory, University of Cambridge,
Madingley Road, Cambridge, CB3 0HE, UK
We determine the uncertainties on observables arising from the errors on the experimental data that are fitted in the global MRST2001 parton analysis. By diagonalizing the error matrix we produce sets of partons suitable for use within the framework of linear propagation of errors, which is the most convenient method for calculating the uncertainties. Despite the potential limitations of this approach we find that it can be made to work well in practice. This is confirmed by our alternative approach of using the more rigorous Lagrange multiplier method to determine the errors on physical quantities directly. As particular examples we determine the uncertainties on the predictions of the chargedcurrent deepinelastic structure functions, on the crosssections for production and for Higgs boson production via gluon–gluon fusion at the Tevatron and the LHC, on the ratio of to production at the LHC and on the moments of the nonsinglet quark distributions. We discuss the corresponding uncertainties on the parton distributions in the relevant domains. Finally, we briefly look at uncertainties related to the fit procedure, stressing their importance and using , and extractions of as examples. As a byproduct of this last point we present a slightly updated set of parton distributions, MRST2002.
1 Introduction
Recently, much attention has been focused on uncertainties associated with the parton distributions that are determined in the nexttoleading order (NLO) global analyses of a wide range of deep inelastic and related scattering data. There are many sources of uncertainty, but they can be divided into two classes: those which are associated with the experimental errors on the data that are fitted in the global analysis and those which are due to what can loosely be called theory errors. In this latter category we have uncertainties due to (i) NNLO and higherorder DGLAP contributions, (ii) effects beyond the standard DGLAP expansion, such as extra and terms, higher twist and saturation contributions, (iii) the particular choice of the parametric form of the starting distributions, (iv) heavy target corrections, (v) model assumptions, such as . In order to estimate some of these ‘theory’ errors, we can also look at the uncertainties arising from different choices of the data cuts (, , ), defined such that data with values of , or below the cut are excluded from the global fit. This approach indicates where the current theory is struggling to fit the data compared to other regions.
Here we study the uncertainties due to the errors on the data, and leave the discussion of the ‘theory’ uncertainties to a second paper. Other groups [1]–[7] have also concentrated on the experimental errors and have obtained estimates of the uncertainties on parton distributions within a NLO QCD framework, using a variety of competing procedures. Of course, the parton distributions are not, themselves, physical quantities. However, using the standard approach of the linear propagation of errors, these uncertainties of the parton distributions can be translated into uncertainties on observables. Therefore, we first follow the general approach in [4] and [5], the Hessian method, and diagonalize the error matrix, parameterizing an increase in of the fit in terms of a quadratic function of the variation of the parameters away from their best fit values. This gives us a number of sets of partons with variations from the minimum in orthogonal directions which can be used in a simple manner to calculate the uncertainty on any physical quantity. However, this approach depends for its reliability on the assumption that the quadratic dependence on the variation of the parton parameters is very good. We find that this approximation, with some modifications of the precise framework, i.e. the elimination of some parameters and rescaling of others, can be made to work well. We make available 30 sets of partons – 2 for each of the 15 eigenvector directions in parton parameter space – which can be used to calculate the uncertainties on any physical quantity.
Despite its convenience, the Hessian approach does suffer from some problems if one looks at it in detail, and if one tries to extrapolate results, in particular if we consider large increases in . It is also not, in principle, the most suitable method when allowing to vary as one of the free parameters in the fit. Hence, in this paper we also investigate the uncertainties on observables directly. In order to do this we apply the Lagrange multiplier method [8] to the observables themselves, therefore avoiding some of the approximations involved in the linear propagation of errors from partons to the observables, and confirming that these approximations do not usually cause serious problems. When using this Lagrange multiplier approach, the resulting sets of parton distributions, which correspond to the extreme values of each observable, can to a certain extent be thought of as the maximum allowed variation of the dominant contributing partons in the relevant kinematic () domain. We select observables which are particularly relevant for experiments at present and future colliders, and which illustrate the uncertainties on specific partons in a variety of kinematic () domains. In order to determine the true uncertainty on quantities we also let vary along with the parameters describing the parton distributions directly, which is easy to implement in this approach. Some quantities are then far more sensitive to than others. Fortunately our global fit [9] produces a value of which is consistent with the world average, with the same type of error, i.e., . Hence, it is completely natural to simply let vary as a free parameter in the fit in the same way as all the other parameters when determining uncertainties. However, we also perform an investigation of the uncertainties with fixed at 0.119 in order to study more directly which variations in the parton distributions are responsible for extreme variations in given physical quantities, and to compare with the results of the Hessian approach.
The physical observables that we select as examples in this introductory study are, first, the chargedcurrent structure functions for deep inelastic scattering at high and at HERA. These observables almost directly represent the , valence quarks at high and , where deep inelastic data do exist [10]–[12], but have errors of or more at present. The precision on these data is expected to increase dramatically in the near future. Second, we determine the uncertainties on the crosssections and , for boson production and for the production of a Higgs boson of mass^{2}^{2}2There is nothing special about the choice of 115 GeV. We may choose different values in order to probe the gluon in different domains. GeV by gluon fusion respectively, at Tevatron and LHC energies. The cross section is sensitive to the sea quarks (and also, at the Tevatron energy, weakly sensitive to the valence quarks) in a range of rapidity centered about , and for . Similarly, is sensitive to the gluon distribution in the domain and .
As a third example we determine the uncertainty on the ratio of to production at the LHC energy. This ratio is expected to be extremely accurately measured in the LHC experiments. Other relevant examples, which we study, are the uncertainties of the moments of the nonsinglet (–) quark distributions. These are quantities for which lattice QCD predictions are becoming available, see, for example, Refs. [13, 14].
The same techniques can be easily and quickly applied to a wide variety of other physical processes sensitive to different partons and different domains. Besides giving a direct evaluation of the uncertainties on the observables, we can, in principle, unfold this information to map out the uncertainties on NLO partons over the whole kinematic domain where perturbative QCD is applicable.
The plan of the paper is as follows. In Section 2 we discuss the Hessian method, and outline our extraction of different parton distribution sets using this approach. In particular we highlight the problems encountered, and how they are dealt with in order to obtain reliable results. We make the sets of partons obtained publicly available. In Section 3 we briefly recall the elements of the Lagrange multiplier method. In the following four sections we determine the uncertainties of the observables that we have mentioned above. This will involve a series of global fits in which the observables are constrained at different values in the neighbourhood of their values obtained in the optimum global fit. In each case we explore, and discuss, the allowed variation of the dominantly contributing partons. Using this more rigorous method we also confirm the general appropriateness of the Hessian approach, but discuss where it can start to break down.
Finally, in Section 8, we summarize and briefly investigate the uncertainties associated with the initial assumptions made in performing the global fit. In order to do this we compare the and boson predictions with those obtained using both a slightly updated set of our own partons, MRST2002, and using the CTEQ6 partons [5]. (All the results in Sections 2–7 are based on MRST2001 partons [9].) We find that for the comparison with CTEQ some of the variations in predictions are surprisingly large. We also illustrate the same result for the extractions of by various different groups. This implies that uncertainties involved with initial assumptions and also with theoretical corrections can be more important than those due to errors on the data.
2 The Hessian method
The basic procedure involved in this approach is discussed in detail in [4], but we briefly introduce the important points here. In this method one assumes that the deviation in for the global fit^{3}^{3}3The data that are fitted can be found in Refs. [6, 10, 11] and [15]–[29]. We treat the errors as in [9]. from the minimum value is quadratic in the deviation of the parameters specifying the input parton distributions, , from their values at the minimum, . In this case we can write
(1) 
where is an element of the Hessian matrix, and is the number of free input parameters. In this case the standard linear propagation of errors allows one to calculate the error on any quantity using the formula
(2) 
where is the covariance, or error matrix of the parameters, and is the allowed variation in . Hence, in principle, once one has either the Hessian or covariance matrix (and a suitable choice of ) one can calculate the error on any quantity.
However, as demonstrated in [4], it more convenient and more numerically stable to diagonalize either the Hessian or covariance matrix, and work in terms of the eigenvectors. Since the Hessian and covariance matrices are symmetric they have a set of orthogonal eigenvectors defined by
(3) 
Moreover, because variations in some directions in parameter space lead to deterioration in the quality of the fit far more quickly than others, the eigenvalues span several orders of magnitude. Hence it is helpful to work in terms of the rescaled eigenvectors . Then the parameter displacements from the minimum may be expressed as
(4) 
or using the orthogonality of the eigenvectors
(5) 
i.e., the are the appropriately normalized combinations of the which define the orthogonal directions in the space of deviation of parton parameters. In practice a is often dominated by a single .^{4}^{4}4CTEQ have even implemented the diagonalization procedure in the fitting procedure itself in order to improve numerical stability [30]. We do not think this will have effects significant enough to outweigh the inherent errors in the Hessian approach described below.
The error determination becomes much simpler in terms of the . The increase in is
(6) 
i.e., the surface of constant is a hypersphere of given radius in space. Similarly the error on the quantity is now
(7) 
Thus it is convenient to introduce parton sets for each eigenvector direction, i.e., from Eq. (4) we define
(8) 
where the tolerance is defined by and is the allowed deterioration in fit quality for the error determination. Then, assuming the quadratic behaviour of about the minimum, (7) becomes the simple expression
(9) 
If everything were ideal this framework would provide us with a simple and efficient method for calculating the uncertainty due to experimental errors on any quantities, where we would use the standard choice of . However, the real situation is not so simple, and there are two major complications we must overcome in order to obtain reliable results.
Although, in principle, the uncertainty in any crosssection should be given by , the complicated nature of the global fitting procedure, where a large number of independent data sets are used, results in this being an unrealistically small uncertainty [31]. This is undoubtedly due to some failure of the theoretical approximation to work absolutely properly over the full range of data, which introduces the type of theoretical errors outlined in the Introduction, and also due to some sources of experimental error not being precisely quantified. Both problems are in practice extremely difficult to surmount. We shall implicitly ignore the potential theoretical error in this paper, but account for the lack of ideal behaviour in the framework by determining the uncertainties using a larger . We estimate to be a conservative uncertainty (perhaps of the order of a confidence level or a little less than ) due to the observation that an increase of 50 in the global , which has a value for data points, usually signifies that the fit to one or more data sets is becoming unacceptably poor. We find that an increase of 100 normally means that some data sets are very badly described by the theory. Though this estimation does not rely on any real mathematical foundation we do not think it is any less valid than the approaches used in e.g. [5] or [1, 7], both of which ultimately appeal to some value judgment rather than using all available information in a statistically rigorous manner, and ultimately give similar results. The approaches [2, 3, 6] do use but either rely on much smaller and more internally compatible data sets, or in some cases have rather small errors.
The second complication is the breakdown of the simple quadratic behaviour in terms of variations of the parameters, i.e., the fact that Eq. (1) may receive significant corrections and the simple linear propagation of errors is therefore not accurate. Of course, we expect some deviations from this simple form for very large , but unfortunately very significant deviations can occur for relatively small , as outlined below. Due to the very large amount of data in our global fit, we have a lot of parameters in order to allow sufficient flexibility in the form of the parton distributions. Each of the valence quarks and the total sea quark contribution are parameterized in the form
(10) 
where for the valence quarks the normalization is set by the number of valence quarks of each type. Because we find it necessary to have a negative input gluon at low the gluon parameterization has been extended to
(11) 
where is determined by the momentum sum rule, and can be set to some fixed large value, e.g. 10 or 20, so that the second term only influences large . The combination has a slightly different parameterization, i.e.,
(12) 
Overall, this gives 24 free parameters. In our standard fits we allow all these parameters to vary. However, when investigating in detail the small departures from the global minimum we notice that a certain amount of redundancy in parameters leads to potentially disastrous departures from the behaviour in Eq. (1). For example, in the negative term in the gluon parameterization very small changes in the value of can be compensated almost exactly by a change in and (to a lesser extent) in the other gluon parameters over the range of probed, and therefore changes in lead to very small changes in . However, at some point the compensation starts to fail significantly and the increases dramatically. Hence, this certain degree of redundancy between and leads to a severe breaking of the quadratic behaviour in . Essentially the redundancy between the parameters leads to a very flat direction in the eigenvalue space (a very large/small eigenvalue of the covariance/Hessian matrix) which means that cubic, quartic etc. terms dominate. During the process of diagonalization this bad behaviour feeds through into the whole set of eigenvectors to a certain extent.
Therefore, in order that the Hessian method work at all well we have to eliminate the largest eigenvalues of the covariance matrix, i.e., remove the redundancy from the input parameters. In order to do this we simply fix some of the parameters at their best fit values so that the Hessian matrix only depends on a subset of parameters that are sufficiently independent that the quadratic approximation is reasonable. In fact we finish up with 15 free parameters in total – 3 for each of the 5 different types of input parton. In particular, fixing the other parameters at the best fit values we find that , and are sufficient for the gluon – one for high , one for medium and one for low . However, we emphasize that we cannot simply set the other parameters to zero. For example must be of a size as to allow a sufficiently negative input gluon at low with a sensible value of , but we cannot allow it to vary simultaneously with . We could possibly allow one or two more parameters to be free, but judge that the deterioration in the quality of the quadratic approximation does not outweigh the improvements due to increased flexibility in the parton variations. We note that this problem seems to be a feature of the full global fits obtained by CTEQ and MRST, and that the other fitting groups have not yet needed to introduce enough parameters to notice such redundancy. It has clearly been noticed by CTEQ though, since in [4] they only have 16 free parameters out of a possible 22, and in [5], where they use a significantly altered type of parameterization, they have only 20 free parameters out of a possible 26.
Hence, we produce 30 sets of parton distributions labeled by to go along with the central best fit; that is 15 ‘‘+’’ sets corresponding to each eigenvector direction, and 15 ‘‘’’ sets^{5}^{5}5In order to produce the errors on the parton distributions a higher numerical accuracy was required than that used when we previously found just the “best fit”. This results in the partons from the central fit being very slightly different to the standard MRST2001 partons, and we label them by MRST2001C. In fact some of the input parameters are quite different to those in the MRST2001 default, but the partons themselves differ by fractions of a percent. This is an example of the redundancy in some input parameters noted above. The 31 parton sets (, MRST2001C) are available at http://durpdg.dur.ac.uk/hepdata/mrs.. Even though we have limited the number of free parameters in the calculation of the Hessian matrix, we note that we still have significant departure from the ideal quadratic behaviour. For the 10 or so lowest eigenvalues of the covariance matrix the quadratic approximation is very good – the distance needed to go along one of the to produce being the expected to good accuracy in both “+” and “–” directions. However, for 4 or 5 of the largest eigenvalues of the covariance matrix, corresponding mainly to the large quark, large gluon and distributions, the absolute scaling and symmetry break down somewhat. In the very worst case of the largest eigenvalue, the scale factors to produce are 9.5 and 4.5 in the two opposite directions. In order to produce the sets corresponding to we have to multiply the parton deviations required for by these scale factors rather than the expected 7.07. (In fact we do this for all 30 sets, but in most cases the scale factor is in the range 6.5–7.5.) Hence, as in [4, 5], this necessitates the supply of both “+” and “–” sets, whereas in the quadratic approximation one could easily be obtained from the other. Indeed from Fig. 9 of [4] it is clear that CTEQ encounter a breakdown of the quadratic behaviour of much the same type that we do.
Using the 30 parton sets corresponding to the 15 eigenvector directions for variations of the partons about the minimum , one can use Eq. (7) to calculate the error for any quantity, assuming an allowed . In fact it has been proposed [32] that one may also account for some of the asymmetry due to departures from quadratic behaviour by replacing Eq. (9) by the slightly more sophisticated form
(13) 
where represents the best fit set of partons. In [32] and [33] examples are discussed where the use of Eq. (13) instead of Eq. (9) leads to not only an asymmetric error, but also a larger uncertainty overall. We find only fairly minor effects, with no real evidence that Eq. (13) leads to markedly more reliable results than Eq. (9), so we use the simpler Eq. (9) in this paper.
As an example of the use of the Hessian method we show in Figs. 1–4 the uncertainty on some of the parton distributions at various values of , namely the distribution, the distribution and the gluon distribution respectively. As one sees, the distribution is very well determined, and the uncertainty shrinks slightly with increasing . The lowest uncertainty is in the region of where there are very accurate data which mainly constrain the valence quarks. At lower the direct constraint is on the sum of valence and sea quarks. The distribution is also well determined in general, but is rather more uncertain as we go to the highest values. The gluon distribution is known less well, but at the highest has an uncertainty of as little as for where it is constrained by both of the HERA data and the lowest Tevatron jet data. It becomes very uncertain for where only the relatively imprecise highest jet data provide any information. The fractional uncertainty at very small decreases very rapidly as increases because much of the small gluon at higher is generated from that at higher via evolution. We also show the gluon at explicitly in Fig. 4. At this low scale the central gluon is negative at , but we see that the gluon may be positive within the uncertainty. This just about persists if we go to as low as at this , but at our input scale the gluon would be negative for less than 0.0005, outside the level of uncertainty chosen. Also shown on the plots are the CTEQ6M partons. For the distribution the agreement is excellent. For the distribution the agreement at is very good, but there is a discrepancy below this value. However, in this range, the valence quarks become very small indeed and the data only really constrain the total distribution which is completely dominated by the sea. This apparent discrepancy is probably due to parameterization effects, and is irrelevant in practice. However, in Fig. 3 we see that the MRST and CTEQ gluons show a genuine and significant level of incompatibility. We will comment on this more in Section 8.
One might worry that the fixing of some of the parameters, that determine the input parton distributions, will cast some doubt on the error obtained. However, we stressed that these are largely redundant parameters, and we have checked that the errors obtained (when using ) are indeed compatible with the errors obtained using more rigorous means, i.e., the Lagrange multiplier method, in the following sections.^{6}^{6}6We have checked the effects of using Eq. (13) instead of Eq. (9) in these comparisons. In all cases the former only introduced a relatively small asymmetry in the uncertainty, with the average being very close indeed to the result using the latter. Also, the asymmetry was of the same sign as that found using the Lagrangian approach only half the time, i.e. the use of Eq. (13) did not reliably predict the direction of steeper increase of , even when the asymmetry was quite large. We find this surprising, and have no good explanation. However, it illustrates the semiqualitative nature of the Hessian approach compared to the more rigorous Lagrange Multiplier method. Nevertheless, it is a sign of the breakdown of the quadratic approximation. Of more practical concern is the fact that this breakdown is also exhibited in a nontrivial manner in some of the eigenvectors used – particularly those eigenvectors associated with the least known partons, e.g. the high down quark and gluon. The scaling has been designed to give correct results if is used. However, one cannot simply extrapolate to different choices of . For example if were deemed a more suitable choice, in principle the error would just be that using Eq. (7) divided by , but the breakdown of quadratic behaviour does not guarantee this, especially for some directions in parameter space. Also, if one wished to be very conservative in the estimation of an uncertainty, simple extrapolation cannot reveal when might start to increase rapidly. We will see examples of this later.
Also we note that we have performed this analysis for a fixed value of the coupling constant: . One can in principle include this as another free parameter. Indeed we then find that the behaviour obeys the quadratic approximation quite well and that gives an error of about , corresponding well to our error of obtained in [9] using . We will discuss extractions of again in Section 8. However, for the Hessian approach there is a slight difference between varying and varying the parton parameters. When is fixed the maximum error on any quantity is obtained from some linear combination of our different parton sets, and in principle one could reproduce the particular parton set which corresponds to this linear combination, which would be a perfectly welldefined set itself. However, a linear combination of coming from contributions with different does not actually correspond to one particular choice of (each contribution has a branch point at a different value of , so a linear combination will have multiple branch points), so one cannot precisely define a particular set of partons corresponding to a particular for the extreme.
Hence, although the 30 parton sets obtained using the Hessian approach provide the most convenient framework for calculating the uncertainties on a physical observable, for the reasons described above we would also like to study an alternative approach, partially just to check how well our adapted Hessian approach really works. A more robust method, which also allows us to directly investigate the partons, and , corresponding to the extreme variations of a given physical quantity is the Lagrange multiplier method. We study this in detail below.
3 Lagrange multiplier method
It is much more rigorous to investigate the allowed variation of a specific observable by using the Lagrange multiplier method. This was also one of the approaches used by the CTEQ collaboration [8]. In this, one performs a series of global fits while constraining the values of one, or more, physical quantities in the neighbourhood of their values obtained in the unconstrained global fit. To be precise, we minimize the function
(14) 
with respect to the usual set of parameters, which specify the parton distributions and the coupling . This global minimization is repeated for many fixed values of the Lagrange multipliers . At the minima, with the lowest , the observables have the values and the value of is the minimum for these particular values of . These optimum parameter sets depend on the fixed values of . Clearly, when , we have and . In this way we are able to explore how the global description of the data deteriorates as the move away from the unconstrained best fit values . Thus by spanning a range of we obtain the profile for a range of values of about the best fit values, . In this study we take the best fit values corresponding to the MRST2001 partons [9].
This procedure involves none of the approximations involved in the Hessian approach. We can use the full set of parameters in the fit, obtaining maximum flexibility in the partons without having to worry about the large correlations or anticorrelations between some parameters. We never make any assumption about quadratic dependence on the parameters, and indeed, by using different values of the Lagrange multipliers, we can map out precisely how the quadratic approximation breaks down in the uncertainty for any physical quantity. Also, one produces a particular set of partons with a particular value of at every point in the space of crosssections for the physical quantities mapped, so the interpretation of the extremes is more obvious and natural. Hence, in principle, this is a far superior method of obtaining uncertainties to the Hessian approach. However, it suffers from the large practical disadvantage that a series of global fits must be done every time one considers a new quantity. As examples we investigate a number of interesting physical cases below.
4 The chargedcurrent structure functions )
The contour plot for the variation of and about their predicted values from the unconstrained global fit is shown in Fig. 5, where we allow to be a free parameter (unstarred labels) or fix it at the best fit value of (starred labels). We show the contours for , 100, etc. Overall, the ellipses one would expect from the quadratic approximation for in Section 2 are more or less what one sees, but there is a certain asymmetry in that increases rather more rapidly for an increase in both and than for a corresponding decrease in both.
Thus, from Fig. 5, we see that the uncertainties of the and structure functions at and (due to the experimental errors on the data in the global fit) are about % and respectively. In comparison, the values using the Hessian approach are and respectively, in good agreement, although slightly smaller for . At this value of the uncertainties in and have a particularly simple interpretation since is almost exactly proportional to the valence down quark distribution, , and is almost exactly proportional to the distribution. This can clearly be seen in Fig. 6, which shows the and distributions for the extreme sets (T*,U*,V* and W*) corresponding to maximum and minimum and (for the case of fixed ). Rather obviously the distribution maximises at large for the case of maximum and minimises for minimum , with similar behaviour for the distribution and . Note however that in each case the extreme in the parton distribution is not precisely at , but at slightly higher , where the data are less constraining. There are also sum rules on the partons which must be satisfied. It is also clear that there is a strong inverse correlation between the and distributions. This is because the data which constrain the relevant partons are the structure function measurements , and which are essentially proportional to , and respectively, where at . This constrains far more than as we have seen, but means that for maximum variation in the partons a change in must be compensated by a much larger opposite change in . The result that the major axis of the ellipse for given change in is approximately aligned along (i.e., ) is therefore not at all surprising. The rate of quickest increase in is then along , where the changes in the partons add in such a way as to maximise changes in the measured structure functions.
We see that allowing to also vary allows the error ellipses to grow slightly, mainly in width. Now the maximum and minimum allowed values of (or ) correspond to and and to parton sets T and V respectively. Most of the constraining data are for , and must be well fit, but smaller means slower evolution of the quarks and thus greater values of and at . Opposite considerations lead to the maximum . Since the extrema of and are more involved, due to the negative correlation with the distribution, they are less altered by allowing to vary; see sets U and W. We see that the axes of the ellipse are essentially unchanged when is left free. Thus Fig. 6 is much the same except that the variations for parton sets T and V are a little greater than for T and V.
It is, of course, the fixed target data which constrain these crosssections and the high quarks. It is very largely the BCDMS measurements which are responsible for the upper extremum in . The best fit tends to overshoot these data in the region of , and a large increase in makes the fit to these measurements very poor. For the extrema in and , the deterioration is more evenly spread over pretty much all the fixed target data at (with the exception that the description of the BCDMS measurements improves slightly), but the cumulative result is a very poor fit. One of the worst instances of deterioration is for the NMC ratio.
5 and production at the LHC and Tevatron
The contour plot for the variation of and about their predicted values at the LHC energy from the unconstrained global fit is shown in Fig. 7, where again we allow to be a free parameter or fix it at . Again we show the contours for , 100, etc. This time the Hessian approach should work well, although the ellipses start becoming a little rectangular. Allowing to vary, we see that the uncertainties of the and crosssections at the LHC (due to the experimental errors on the data in the global fit) are about % and respectively, and are positively correlated.
Again this analysis also gives information on the uncertainties of particular parton distributions. To be specific, the parton sets which correspond to the points A,B,C,D, on the contour in Fig. 7, give the uncertainties in the parton distributions that dominantly determine and in the kinematic domain , GeV relevant to and production at the LHC. The extrema in , represented by A and C, correspond to variations in the sea quark distributions, while the extrema in , represented by B and D, correspond to variations in the gluon distribution and . The values of for sets A and C are 0.119 and 0.118 respectively, both very close to the default MRST2001 value, showing that , which begins at zeroth order, is insensitive to . However, the values of for fits B and D are 0.120 and 0.117 respectively, reflecting the fact that . This is well illustrated by repeating the entire analysis with fixed at the default value (0.119) obtained in the unconstrained global fit [9]. The and 100 contours for this additional analysis are shown by the smaller shapes in Fig. 7. We can see that the uncertainty on is almost unchanged, while that for is reduced to about . The corresponding values using the Hessian approach are and , in good agreement but slightly smaller in each case.
The up quark distribution for each ‘extrema’ set with fixed is shown in Fig. 8(a) and the gluon distribution in Fig. 8(b). We see that indeed the parton distributions do reflect the extrema in the crosssections in a fairly simple manner. The quark densities at high show almost no variation between fits since they are well constrained at high and because the and production crosssections are sensitive to the partons at an range centered at a few . Indeed the maximum and minimum crosssections correspond to the maximum and minimum sea quarks for at . The maximum and minimum Higgs crosssections correspond to the maximum and minimum gluon distributions in the same sort of range, although the large gluon must now decrease for increases in the small partons, and vice versa, in order to maintain the momentum sum rule. The strong correlation between the two crosssections is due to the fact that at high the size of the quark distribution at small is mainly determined by evolution, and the larger the small gluon the stronger the quark evolution (and vice versa). When is left free the resulting partons at the extrema are similar to the fixed results. However, in this case, their variation is a little larger at smaller , since the slight changes in lead to different rates of evolution.
For the case of fixed the main contributions to come from the HERA small structure function data and, because of the changes in the high gluon, also from the Tevatron jet data. For the upper extrema in and the slope tends to be too large for , while for the lower extrema the slope is too small. In both cases the fit to jet data deteriorates due to the shape of the high gluon becoming wrong. When is allowed to vary the data which are particularly sensitive to this also play a role, for example the BCDMS data are fitted less well when in fit B, and the NMC data are described less well when in fit D.
The corresponding contour plot for the Tevatron is shown in Fig. 9, where again we either allow to be a free parameter or fix it at . We see that the uncertainty of the crosssection at the Tevatron (due to the experimental errors on the data in the global fit) has decreased to about while that for the Higgs has increased to about % for varying , and that the correlation has disappeared. For fixed at 0.119 is again largely unaffected, but the uncertainty of now more than halves to about % , reflecting the fact that this time the maximum and minimum Higgs crosssections for variable correspond to and respectively. With fixed there is now even a very slight anticorrelation between the crosssections.
The extrema in , represented by P and R, correspond roughly to variations in the quark distributions at , while the extrema in , represented by Q and S, correspond to variations in the gluon distribution at and . The values of sampled at the Tevatron are an order of magnitude greater than at the LHC. This, coupled with the fact that it is a proton–antiproton collider, rather than a proton–proton collider, complicates the interpretation of the extremes of the crosssections in terms of partons.
The up quark distribution for each extrema set with fixed is shown in Fig. 10(a) and the gluon distribution in Fig. 10(b). The corresponding distributions obtained when is allowed to vary are shown in Fig. 11. We first consider the cases of the maximum and minimum crosssections, which are insensitive to whether is left to vary or not. For discussion purposes, let us consider only the and quark flavour contributions. Then the crosssection is roughly proportional to
(15) 
where 1 refers to the proton and 2 to the antiproton and . Hence the average value of . This is sufficiently large that there is a distinct difference between the quark and antiquark distributions, and the contribution to the crosssection from the quark contribution is the greater. Hence, one can decrease the crosssection by replacing a quark by its antiquark at , or vice versa. Of course, there is a fundamental constraint in doing this due to the sum rule for each valence quark. However, the only real experimental constraint is from the CCFR data, all other structure function data being insensitive to the distinction between quark and antiquark. In the optimum global fit most data would like there to be more quarks at high , while the CCFR data would prefer more valence quarks at . This leads to a compromise where for the best fit the CCFR data at low are undershot. The minimum is therefore achieved mainly by this exchange of quark for antiquark, which most data are happy with, and hence the deterioration in at P (and P) is almost entirely from the description of the CCFR data. Hence, both the gluon and quark distribution for P (and P) are hardly changed, as seen in Figs. 10 and 11, but and decrease for . Going in the other direction, an increase in and the consequent decrease in the valence quarks at higher causes a large penalty in and the maximum is achieved in a different manner. At the quark evolves much more slowly than at and the density at is determined largely by the input value, and modified by the rate of evolution. Hence the maximum is achieved by having a large quark distribution at at low and also by having an enhanced gluon at to increase evolution. These are displayed in Figs. 10 and 11. The deterioration in then comes mostly from the low quarks causing an overshooting of NMC structure function data, but there is also a contribution due to the enhanced gluon at causing it to be smaller for and hence fitting the Tevatron jet data less well.
The extrema of the Higgs crosssection are also slightly complicated. It is not possible to simply increase or decrease the gluon in a range centered on because this is precisely the region where the majority of the gluon’s momentum is carried, and this total is very well constrained by the momentum sum rule and the accurate high quark determination. Therefore, for fixed the change in is largely reliant on the fact that this total crosssection actually probes quarks within about an order of magnitude either side of the central production value of . Hence, as we see from Fig. 10 the maximum crosssection is obtained from the gluon in set Q which is slightly reduced for and more enhanced for and the minimum crosssection is obtained from the gluon in set S which is slightly increased for and more reduced for . In both cases those data sets sensitive to the small and large gluon, i.e., HERA structure function data and Tevatron jet data respectively, are those for which the description deteriorates. When is allowed to go free it varies by about and there is a large increase in the variation of . This is not only because but also because the HERA data anticorrelate and the small gluon. Therefore, in set Q, for example, the increased value of allows the small gluon to get much smaller, and the high gluon much larger, compared to set Q. This compensation between and the small gluon also means that HERA data remains well fit, and it is the jet data (particularly CDF), sensitive to large , and the large phobic BCDMS data, for which the description deteriorates. Similar considerations apply to set S as compared to S. Here it is the D0 jet data and the small phobic SLAC and NMC data that are badly fit.
For the Hessian approach gives an uncertainty of for and for , at the Tevatron energy. In simplistic terms this is in good agreement, but a little smaller for the gluonsensitive Higgs crosssection. However, in this case we see from Fig. 9 a very marked asymmetry on the contour plot. For fixed the ellipses are certainly not centered on the best fit values, and for varying we see that is clearly increasing far more rapidly for increases in the predicted crosssection than for corresponding decreases. Thus, it is clear that within the framework of this fit, increases of the crosssection of much more than are completely ruled out, whereas decreases of the same amount are much more acceptable. This information would be largely lost in the Hessian approach, and for these quantities the Lagrange multiplier method does supply some important additional information.
6 The ratio of to production at the LHC
The ratio of the to the production crosssections at hadron colliders is a particularly interesting observable. The measurement is expected to be quite precise (better than at the LHC, see e.g. [34]), since many of the experimental uncertainties cancel in the ratio. The uncertainty in the prediction of the ratio at the LHC can be deduced from the profile shown in Fig. 12. Taking, as before, the measure, we obtain , and the Hessian approach is in very good agreement with this. Since the ratio is sensitive to the ratio of the and quark distributions, it is not surprising that the increase in is almost entirely due to the NMC data [25].
A detailed discussion of the ratio may be found in Ref. [35]. Consider, for instance, the ratio as a function of the rapidity
(16) 
where at the LHC. In Eq. (16) we have ignored, for simplicity, the contributions involving strange and heavier quarks. Thus a measurement of the ratio at large would provide a direct determination of at large . For example, for , we measure at at the LHC. Of course, it is the decay lepton rapidity that is measured, rather than the parent rapidity, and so the ratio in a given rapidity bin will have a greater uncertainty than that for .
7 The moments of the (–) distribution
The parton distribution functions of the nucleon are fundamental quantities that should, in principle, be calculable from first principles in QCD. In particular, the moments of parton distributions at a given scale are related, by the operator product expansion, to a product of perturbatively calculable Wilson coefficients and nonperturbative matrix elements of quark and gluon operators. The latter can be computed using lattice QCD and, indeed, in recent years the precision of the lattice calculations has improved significantly. Although in principle the lattice results can be related to moments of physical structure functions, in practice it is more efficient to use parton distributions determined in a global fit to represent the physical ‘data’. Comparisons of recent lattice moment calculations with the predictions of earlier MRS parton distributions are encouraging, see for example [13, 14].
In order to quantify the agreement between the lattice calculations and the parton distribution predictions it is obviously important to know the uncertainties in the latter. It is straightforward to apply the Lagrange multiplier method used in previous sections to determine the uncertainties in observable crosssections to the moments of parton distributions.
To avoid contamination from gluon contributions, the lattice calculations focus on the moments of nonsinglet quark operators. For example, lattice results are available for the first three moments of the combination , i.e.,
(17) 
with . The predictions of the MRST2001 set (at ) for these moments are given in Table 1.
The contour plot for the (percentage) variation of the second and third moments about their predicted values is shown in Fig. 13. We again show the and contours corresponding to the fixed analysis, but there is evidently little difference between the fixed and variable coupling results in this case.
As expected, there is a strong positive correlation between the two moments. Using the , varying criterion for defining a conservative error, we obtain errors of , and for the second, third and fourth moments respectively. The corresponding predictions for the errors on the moments are also given in Table 1. The increasing relative error with increasing moment is to be expected – higher moments probe the region where there are fewer DIS structure function data. Again we notice that there is a small asymmetry in the contours – the increase in when both moments increase being less severe than when both moments decrease.
% error  

2  0.1655(70)  4.2 
3  0.0544(26)  4.8 
4  0.0232(12)  5.0 
The uncertainties on the moments using the Lagrange multiplier method with a fixed are slightly smaller: , and for the second, third and fourth moments respectively. These results are in excellent agreement with the (fixed ) Hessian approach, where the corresponding errors are , and .
Since, as we have already seen in Section 4, the quark at high is far more constrained than the quark, the allowed variation in these moments is mainly due to variations in the distribution. The minimum extremum (H in Fig. 13) of the moments is therefore due to the largest allowed distribution at high and arises from a similar set of partons to those for the maximum . Thus, as in this previous case, it is mainly the comparison to the BCDMS measurements which causes the deterioration in the quality of the fit. The maximum of the moments (G in Fig. 13) corresponds roughly to the minimum distribution at high and it is largely the fit to the ratio that breaks down.
For a number of years, lattice QCD has been used to calculate the moments of nucleon structure functions from first principles. The most recent comprehensive results are from the LHPCSESAM [13] and QCDSF [14] collaborations. Although the comparisons with experiment (via parton distributions obtained from global fits) are encouraging, there are still many problems to be overcome, for example finite lattice spacing and volume effects, renormalization and mixing of operators, unquenching and chiral extrapolation to physical quark masses. A comparison with the recent lattice results [13, 14] and the above MRST2001 moment predictions reveals that (a) the errors in the latter are at present significantly smaller than in the former, especially for the higher moments, and (b) the lattice results for the moments are systematically higher. The explanation appears to be that the linear chiral extrapolation used in the lattice determinations is not valid – nonperturbative longdistance effects in the nucleon gives rise to nonlinear, nonanalytic dependence on [36]–[40] which is particularly important at small . In the most recent analyses (see for example the comprehensive study in [41]), the experimental (i.e., pdf) values for the moments are used to constrain a priori unknown nonperturbative parameters which enter in the nonanalytic terms in the chiral extrapolation formula. It will be interesting to investigate the effect of using the new MRST2001 moment predictions and errors in such studies.
8 Comparison between different central parton sets
So far in this paper we have investigated the uncertainty on physical quantities coming from the experimental error of the measurements used to determine the parton distributions. We have discussed both the Hessian and Lagrange Multiplier approaches, concluding that the latter is in principle preferable, but recognizing the practical advantages of the former. We have compared the results each provide for the uncertainties using the criterion, noting that they are generally in good agreement. The Hessian approach does tend to give slightly smaller uncertainties for the quantities sensitive to the least welldetermined partons, i.e., which is sensitive to the gluon distribution and which is sensitive to the high down quark distribution. This is probably partly due to the neglected effect of the not entirely redundant parameters, and partly due to errors associated with those eigenvectors which do not respect the quadratic approximation for too well, which indeed are mainly concerned with the gluon and high down quark. However, the discrepancy is quite small, and we judge that we can trust the Hessian approach, at least for in the region of 50 or less, to give quantitative results. Hence, for fixed , we have made available 30 parton sets corresponding to the 15 different eigenvector directions in the space of variation of parton parameters away from their values at the minimum of the global fit, each set corresponding to an increase in of 50. These can easily be used to obtain the error on any physical quantity, as outlined in Section 2. We have also made available various parton sets with fixed and varying corresponding to extreme variations in the predictions for various important crosssections and other relevant observables.
We note that the uncertainties obtained due to the errors on the experimental data are generally very small, of the order of , except for quantities sensitive to the high down quark and gluon, where they can approach . However, in all of this we have implicitly assumed that the theoretical procedure is precisely compatible with the data used, we have not considered the uncertainties due to (i) the data sets chosen, (ii) the choice of starting parameterizations, (iii) the heavy target corrections, etc. In practice this is far from true, as discussed in the Introduction. In this final section we acknowledge this to some extent and investigate qualitatively the impact of the initial assumptions going into the fit on the uncertainty on some quantities. In order to do this we first perform a slightly updated fit of our own (which includes minor modifications in terms of parameterization and the treatment of errors and data sets) so as to produce the best set of uptodate partons. This was partially inspired by the question of why CTEQ6 [5] gives a much better fit to the Tevatron jet data than MRST2001, but also by the availability of new ZEUS data. We call the new set MRST2002 partons.^{7}^{7}7The MRST2002 parton set can be found at http://durpdg.dur.ac.uk/hepdata/mrs.
8.1 CTEQ6, MRST2001 and a new parton set (MRST2002)
We found that we can improve the fits to jets within the global fit by a couple of modifications. In order to obtain the best global fit with partons input at we had previously found that we needed a parameterization which allows the gluon to go negative at small . Hence we used
(18) 
where , and was fixed at , so as not to affect the high distribution. Unexpectedly, allowing to vary to resulted in a slight improvement in the fit to Tevatron jets. We also modified our treatment of the errors for the Drell–Yan data [28]. The fit to these data actually competes with that to the jets, and using only statistical errors, as in our previous studies (the systematic errors being defined a little vaguely), overemphasizes the effect of the Drell–Yan measurements. Adding systematic errors in quadrature to the statistical errors (which is probably the best approach [28]) also improves the fit to the jet data. Both these modifications appear appropriate and are implemented in our updated set. Also included in the new analysis is the new ZEUS high data [42], which has little effect on the partons. The only significant change in the MRST2002 partons, compared to MRST2001 partons [9], is an increase in the gluon at high , which we show in Fig. 14. The fit to the Tevatron jet data now has compared to for MRST2001, and the fit to The Drell–Yan data with systematic errors has . The quality of fit for all other data sets is almost identical to that for the MRST2001 partons.
The CTEQ6 partons are very similar to the MRST2001 (and MRST2002) partons in most aspects. However, in this CTEQ analysis [5] a number of different choices are made about the way in which the fit is implemented, which leads mainly to a significantly different gluon distribution. These differences are: the development of a different type of parameterization for the partons, which allows for a different shape at very high ; CTEQ omit data below , compared to our choice of ; they do not fit to some data sets used in [9], i.e., they omit SLAC and one H1 high set of data; they use systematic errors (in quadrature) for Drell–Yan data; moreover, CTEQ have a positivedefinite small gluon at their starting scale of . They also use a massless charm prescription and there are various other minor differences as compared with MRST.^{8}^{8}8The way in which these different assumptions lead to an improved fit to the Tevatron jet data is outlined in [43].
The CTEQ6 gluon is also shown in Fig. 14. Clearly MRST2002 has a similar high gluon to CTEQ6, both being larger than MRST2001. However, the MRST gluons are different from the CTEQ6 gluon at smaller due to their freedom to have a negative input distribution, and due to slight differences in the choice of data sets fitted. The different assumptions made in obtaining the CTEQ partons, although they improve the quality of the jet fit, do not lead to the best fit when including the data sets omitted by CTEQ and the fit is not good at all for data with . Hence, within the context of trying to obtain as inclusive a global fit as possible using NLO QCD, we take MRST2002 to be the best set of parton distributions.
8.2 Comparison of predictions for and for
The predictions for and Higgs crosssections using the different partons are shown in Fig. 15. Since MRST2002 only differs from MRST2001 in the high gluon, to which these crosssections are insensitive, the predictions for MRST2002 are very similar to those of MRST2001. (Hence our decision to keep MRST2001 partons as the base set for this paper). However, the corresponding predictions obtained using the CTEQ6 partons are quite different. At the LHC the prediction for is similar, but is towards the top of our (qualitative) confidence level. From Fig. 14 this is clearly due to the larger gluon in the region, which is due to the positive definite input for the CTEQ6 gluon. At the Tevatron the discrepancy between CTEQ6 and MRST is even larger. The CTEQ6 predictions for both and are effectively completely outside our expectations. The reason for the small prediction of is evident from Fig. 14—the CTEQ6 gluon is considerably smaller in the region of . This, in turn, is then responsible for a slower evolution of the quarks, making them smaller at high and hence making smaller. Presumably the difference comes about because CTEQ6 use a more restricted form of the gluon and omit one H1 data set and data which prefer larger . Whatever the precise reasons for the discrepancies, it is clear that different choices for the overall framework of the global fit can completely outweigh the uncertainties due to errors on the data actually chosen to go into the fit. It would be easy to illustrate similar types of discrepancy comparing to other alternative sets of partons—in particular, due to the absence of the Tevatron jets in the fits, many of the parton sets in [1]–[7] have rather smaller gluons at large , and would have different predictions for various quantities sensitive to the high gluon.
Group  variation  

of  
CTEQ6  
ZEUS  
MRST02  
MRST01  
H1 

Alekhin  
GKK 
8.3 Comparison of predictions for
We also find a large variation in the values of extracted from the fits of the different collaborations: CTEQ6 [5], ZEUS [7], MRST2001 [9], H1 [6], Alekhin [3] and Giele et al. (GKK) [2]. The resulting values of are listed in Table 2, together with the determination of this work (MRST2002), in order of decreasing tolerance (), which is reflected in the size of the corresponding experimental error. Not all are presented as determinations of , but all are extracted using the same criteria as for the uncertainty on partons in the respective fit, and hence should be as reliable. Clearly there is a very large variation, with some very low values. The uncertainties due to experimental errors are determined in different fashions in each case, and a summary can be found in [44]. We use for the ZEUS determination [7], because they use the offset method for determining uncertainties which for gives a larger uncertainty than the more common Hessian method. ZEUS estimate that this is equivalent to if they were to use the same treatment of errors as CTEQ. We also use for the GKK value [2], because the uncertainties are obtained using confidence limits, but the error quoted corresponds to the one sigma usually associated with .
The model errors incorporate such effects as the heavy quark prescription and masses, parameterizations, changes in the starting scale of evolution etc. The theory error is often determined by variation of renormalization and factorization scales, though MRST use an estimate appealing to current knowledge of NNLO and resummations, which we feel is more reliable. Since each fit is centered on NLO QCD with scales equal to , the “theory errors” are very strongly correlated, and cannot therefore be responsible for the differences. These discrepancies are undoubtedly due to the assumptions going into the fits, mainly on which data sets are included and which cuts on and are used.
MRST, who obtain the largest value of , use the widest range of data sets and also the least conservative cuts.^{9}^{9}9The slightly different treatment in this work (MRST2002) leads to a marginal raising of as compared to MRST2001 [9], as seen in Table 2. We still use for our onesigma uncertainty, since if corresponds to 90% confidence level, or 1.65 sigma, simple scaling implies that one sigma corresponds to , i.e. to a good approximation. CTEQ use only a slightly smaller number of data sets but also cut data below , as described previously. They also use a definition of the NLO coupling which truncates the solution of the renormalization group equation, whereas most other groups use the full solution of the NLO equation. Both approaches are equally correct, but the truncation of the solution leads to a slightly higher value of at scales below , for the same value of , than the other method, and thus tends to yield a lower . CTEQ also have a very conservative estimate of the error, though it is meant to be somewhat more than a onesigma error. ZEUS and Alekhin use a similar selection of data sets, i.e., HERA DIS data (only ZEUS data in the former case) and a number of fixed target DIS data sets. Hence, it is unsurprising that they obtain similar central values of , with respective errors which are easily explained by their choices of . H1 and GKK both use a small number of sets of data: the former collaboration uses H1 DIS data [6, 10] and BCDMS fixedtarget proton DIS data [22], while GKK use older H1 DIS data [45] together with BCDMS and E665 [26] fixedtarget proton DIS data. Both determinations are heavily influenced by the BCDMS proton data set which prefers rather small^{10}^{10}10Recall the determination from BCDMS data alone [46]. , and this feeds into the final values. Also, both are strict in their statistical interpretation, obtaining small uncertainties, even with relatively small data samples. Finally we note that only CTEQ and MRST include the Tevatron jet data in their analyses. This is relevant because of the –gluon correlation.
8.4 Final comment
From the discussion of the previous two subsections, it is clear that different ideas about the best way to perform a NLO fit can lead to a wide variation in both the central values and the errors of as well as in predictions for physical quantities such as and . The fact that the various ‘NLO’ fits can yield such different outputs is disturbing, and is indicative of the uncertainty arising from theoretical assumptions. Indeed, we have always believed that ‘theory’, rather than experiment, will provide the dominant source of error [44]. We have already produced approximate NNLO parton distributions and predictions [47] (based on the approximate splitting functions [48] obtained from the known NNLO moments [49]), and find, for example, that the NNLO crosssection at the Tevatron is higher than at NLO, and believe that this result is reliable. This change is somewhat larger than the uncertainty due to experimental errors shown in Fig. 9. Moreover, production is likely to be subject to smaller theoretical uncertainty than many other observables—particularly those directly related to the gluon. Our estimates for the uncertainty in at small are or more even at , and significantly larger at lower , for example. Hence, an understanding of theoretical uncertainties is clearly a priority at present, and a preliminary attempt at this will be the subject of a future publication [50].
Acknowledgements
We would like to thank Mandy CooperSarkar, Vladimir Chekelian, Dan Stump and Wuki Tung for useful discussions. RST would like to thank the Royal Society for the award of a University Research Fellowship. RGR would like to thank the Leverhulme Trust for the award of an Emeritus Fellowship. The IPPP gratefully acknowledges financial support from the UK Particle Physics and Astronomy Research Council.
References
 [1] M. Botje, Eur. Phys. J. C14 (2000) 285.

[2]
W.T. Giele and S. Keller, Phys. Rev. D58 (1998) 094023;
W.T. Giele, S. Keller and D.A. Kosower, hepph/0104052.  [3] S.I. Alekhin, Phys. Rev. D63 (2001) 094022.
 [4] CTEQ Collaboration: J. Pumplin et al., Phys. Rev. D65 (2002) 014013.
 [5] CTEQ Collaboration: J. Pumplin et al., JHEP 0207:012 (2002).
 [6] H1 Collaboration: C. Adloff et al., Eur. Phys. J. C21 (2001) 33.

[7]
A.M. CooperSarkar, hepph/0205153, J. Phys.
G28 (2002) 2669;
ZEUS Collaboration: S. Chekanov et al., Phys. Rev. D67 (2003) 012007 .  [8] CTEQ Collaboration: D. Stump et al., Phys. Rev. D65 (2002) 014012.
 [9] A.D. Martin, R.G. Roberts, W. J. Stirling and R.S. Thorne, Eur. Phys. J. C23 (2002) 73.
 [10] H1 Collaboration: C. Adloff et al., Eur. Phys. J. C13 (2000) 609.
 [11] H1 Collaboration: C. Adloff et al., Eur. Phys. J. C19 (2001) 269.

[12]
ZEUS Collaboration: J. Breitweg et al., Eur.
Phys. J. C12 (2000) 411;
ZEUS Collaboration: S. Chekanov et al., Phys. Lett. B539 (2002) 197.  [13] D. Dolgov et al., Phys. Rev. D66 (2002) 034506.
 [14] S. Capitani et al., heplat/0111012, Nucl. Phys. Proc. Suppl. 106 (2002) 299.
 [15] ZEUS Collaboration: S. Chekanov et al., Eur. Phys. J. C21 (2001) 443.
 [16] D0 Collaboration: B. Abbott et al., Phys. Rev. Lett. 86 (2001) 1707.
 [17] CDF Collaboration: T. Affolder et al., Phys. Rev. D64 (2001) 032001.
 [18] ZEUS Collaboration: J. Breitweg et al., Eur. Phys. J. C12 (2000) 35.
 [19] CCFR Collaboration: U.K. Yang et al., Phys. Rev. Lett. 86 (2001) 2742.
 [20] NuTeV Collaboration: M. Goncharov et al., Phys. Rev. D64 (2001) 112006.
 [21] E866 Collaboration: R.S. Towell et al., Phys. Rev. D64 (2001) 052002.
 [22] BCDMS Collaboration: A.C. Benvenuti et al., Phys. Lett. B223 (1989) 485.
 [23] BCDMS Collaboration: A.C. Benvenuti et al., Phys. Lett. B236 (1989) 592.
 [24] L.W. Whitlow et al., Phys. Lett. B282 (1992) 475, L.W. Whitlow, preprint SLAC357 (1990).
 [25] NMC Collaboration: M. Arneodo et al., Nucl. Phys. B483 (1997) 3; Nucl. Phys. B487 (1997) 3.
 [26] M.R. Adams et al., Phys. Rev. D54 (1996) 3006.
 [27] CCFR Collaboration: W.G. Seligman et al., Phys. Rev. Lett. 79 (1997) 1213.
 [28] E605 Collaboration: G. Moreno et al., Phys. Rev. D43 (1991) 2815.
 [29] CDF Collaboration: F. Abe et al., Phys. Rev. Lett. 81 (1998) 5744.
 [30] CTEQ Collaboration: J. Pumplin et al., Phys. Rev. D65 (2002) 014011.
 [31] R.S. Thorne et al., hepph/0205233, J. Phys. G28 (2002) 2717.

[32]
P. M. Nadolsky and Z. Sullivan,
“PDF uncertainties in production at Tevatron,”
in Proceedings of Snowmass 2001: the Future of Particle Physics,
edited by N. Graf (Stanford Linear Accelerator Center, Stanford, 2002),
eConf C010630, P510, hepph/0110378;
Z. Sullivan and P. M. Nadolsky, in Proceedings of Snowmass 2001: the Future of Particle Physics, edited by N. Graf (Stanford Linear Accelerator Center, Stanford, 2002), eConf C010630, P511, hepph/0111358.  [33] Z. Sullivan, Phys. Rev. D66 (2002) 075011.
 [34] M. Dittmar, F. Pauss and D. Zurcher, Phys. Rev. D56 (1997) 7284.
 [35] A.D. Martin, R.G. Roberts, W.J. Stirling and R.S. Thorne, Eur. Phys. J. C14 (2000) 133.
 [36] W. Detmold, W. Melnitchouk, J.W. Negele, D.B. Renner and A.W. Thomas, Phys. Rev. Lett. 87 (2001) 172001.
 [37] W. Detmold, W. Melnitchouk and A.W. Thomas, Eur. Phys. J. C13 (2001) 1.
 [38] A.W. Thomas, W. Melnitchouk and F.M. Steffens, Phys. Rev. Lett. 85 (2000) 2892.
 [39] D. Arndt and M.J. Savage, nuclth/0105045, Nucl. Phys. A697 (2002) 429.
 [40] J.W. Chen and X. Ji, Phys. Lett. B523 (2001) 107.
 [41] W. Detmold, W. Melnitchouk and A.W. Thomas, Phys. Rev. D66 (2002) 054501.
 [42] ZEUS Collaboration: S. Chekanov et al., hepex/0208040.
 [43] R.S. Thorne, A.D. Martin, R.G. Roberts, and W. J. Stirling, Acta. Physica.Polon. B33 (2002) 2927.
 [44] R.S. Thorne, hepph/0205235, J. Phys. G28 (2002) 2705.
 [45] H1 collaboration: S. Aid et al., Nucl. Phys. B470 (1996) 3.
 [46] M. Virchaux and A. Milsztajn, Phys. Lett. B274 (1992) 221.
 [47] A.D. Martin, R.G. Roberts, W.J. Stirling and R.S. Thorne, Phys. Lett. B531 (2002) 216.
 [48] W.L. van Neerven and A. Vogt, Phys. Lett. B490 (2000) 111.

[49]
S.A. Larin, P. Nogueira, T. van Ritbergen and J.A.M.
Vermaseren, Nucl. Phys. B492 (1997) 338;
A. Rétey and J.A.M. Vermaseren, Nucl. Phys. B604 (2001) 281.  [50] A.D. Martin, R.G. Roberts, W.J. Stirling and R.S. Thorne, in preparation.