text:
# Contents

## 1 Introduction

All calculations of high energy processes with initial hadrons, whether within the Standard Model (SM) or exploring New Physics, require parton distribution functions (PDFs) as an essential input. The reliability of these calculations, which underpins both future theoretical and experimental progress, depends on understanding the uncertainties of the PDFs. The assessment of PDF uncertainties has, therefore, become an important challenge to high energy physics in recent years.

The PDFs are derived from global analyses of experimental data from a wide range of hard processes in the framework of perturbative quantum chromodynamics (PQCD). Quantifying the uncertainties in a global QCD analysis is far from being a straightforward exercise in statistics. There are non-gaussian sources of uncertainty from perturbation theory (e.g., higher order and power law corrections), from choices of parametrization of the nonperturbative input (i.e., initial parton distributions at a low energy scale), from uncertain nuclear corrections to experiments performed on nuclear targets, and from normal experimental statistical and systematic errors. These sources of error need to be studied individually, and eventually combined in a systematic way.

We shall be concerned in this paper with uncertainties of PQCD predictions due to uncertainties of PDFs arising from experimental measurement errors. This problem is considerably more complicated than it appears on the surface. The reason is that in a *global analysis*, the large number of data points ($`1300`$ in our case) do not come from a uniform set of measurements, but consist of a collection of measurements from many experiments ($`15`$) on a variety of physical processes ($`56`$) with diverse characteristics, precision, and error determination. The difficulty is compounded by a large number of fitting parameters ($`16`$) which are not uniquely specified by the theory. Several approaches to this problem have been proposed, with rather different emphases on the rigor of the statistical method, scope of experimental input, and attention to various practical complications . Our group has initiated one of these efforts, with the emphasis on utilizing the full constraints of the global data . This work has motivated a closer examination of the standard techniques of error analysis, and necessary improvements and extensions to these techniques, as applied to a complex real world problem such as global QCD analysis of PDFs .

In this paper we present a detailed analysis of uncertainties of physical observables due to parton distribution functions, using the Lagrange Multiplier method proposed in . This method explores the entire multi-dimensional parton parameter space, using an effective $`\chi ^2`$ function that conveniently combines the global experimental, theoretical, and phenomenological inputs to give a quantitative measure of the goodness-of-fit for a given set of PDF parameters. (Cf. Sec. 2.) The method probes directly the variation of the effective $`\chi ^2`$ along a specific direction in the PDF parameter space—that of maximum variation of a specified physical variable. The result is a robust set of *optimized sample PDFs* (or “alternative hypotheses”) from which the uncertainty of the physical variable can be assessed quantitatively without the approximations inherent in the traditional error matrix approach. For concreteness, we consider the cross section $`\sigma _W`$ of $`W`$ boson production at the Tevatron as the archetypal example. (Cf. Sec. 3.)

The definition of the effective $`\chi ^2`$ function, and the inputs that go into it, do not permit a direct statistical interpretation of its numerical value. To obtain meaningful confidence levels for the optimized sample PDF sets, it is necessary to conduct a series of likelihood analyses of these sample PDFs, using all available information on errors for the individual experiments. The results from these analyses serve as the basis to assign an overall uncertainty range on the physical variable, and a corresponding tolerance measure for the effective $`\chi ^2`$ function used in the analysis, that are consistent with the experiments used in the current global QCD analysis. (Cf. Sec. 4.)

This method can be applied to any physical variable, or combination of physical variables, in precision QCD phenomenology. In Sec. 5 we present results on $`W`$ production at the LHC, and $`Z`$ production at the Tevatron and the LHC. We compare the uncertainties obtained in all cases, and comment on previous estimates in the context of these results. In Sec. 6 we present parton distribution sets that are optimized to give high/low values of the $`W`$ and $`Z`$ cross sections, while remaining consistent with current experiments according to our analysis.

The Lagrange Multiplier method provides a useful tool to test the reliability of the more traditional method of error propagation via the error matrix , which relies on the quadratic expansion of the $`\chi ^2`$ function around its minimum. In a companion paper  we perform an in-depth analysis of the uncertainties of the PDFs in the error matrix approach, using the much improved numerical method for calculating the Hessian that was developed in . There we demonstrate how the more specialized Lagrange Multiplier method can set useful benchmarks for the general purpose error matrix approach.

## 2 The Global QCD Analysis

We adopt the same experimental and theoretical input as the CTEQ5 analysis : 15 data sets from 11 experiments on neutral-current and charged-current deep inelastic scattering (DIS), lepton-pair production (Drell-Yan), lepton asymmetry in $`W`$-production, and high $`p_T`$ inclusive jet production processes are used. (Cf. Table 1 in Sec. 4.) The total number of data points is $`N=1295`$. We denote the experimental data values by $`\{D\}=\{D_I;I=1,\mathrm{},N\}`$. The theory input is next-leading-order (NLO) PQCD, and the theory value for the data point $`I`$ will be denoted by $`T_I`$. The theory depends on a set of parameters $`\{a\}\{a_i;i=1,\mathrm{},d\}`$. These parameters characterize the nonperturbative QCD input to the analysis; they determine the initial PDFs $`\{f(x,Q_0;\{a\})\}`$ defined at a low energy scale $`Q_0`$, below the energy scale of the data, which we choose to be $`Q_0=1`$ GeV. When we need to emphasize that the theoretical values depend on the PDF parameters we write $`T_I(a)`$ to indicate the dependence on $`\{a\}`$.

The parametrization of $`\{f(x,Q_0)\}`$ is somewhat arbitrary, motivated by physics, numerical considerations, and economy. Another parametrization might be employed, and differences among the possible parametrizations are in principle a source of theoretical uncertainty in their own right. For most of this study we focus on a single parametrization, but we comment on the effect of changing the parametrization at the end of Sec. 4. The number $`d`$ of the parameters $`\{a\}`$ is chosen to be commensurate with current experimental constraints. For this study we use $`d=16.`$ The detailed forms adopted for the initial functions $`\{f(x,Q_0;\{a\})\}`$ are not of particular concern in this study, since we shall be emphasizing results obtained by ranging over the full parameter space.aaaIn other words, for this paper, the PDF parameters $`\{a\}`$ play mostly the role of “internal variables”. In contrast, they occupy the center stage in the companion paper . The explicit formulas are given in Appendix C (where relevant PDFs from the results of our study are presented). The $`T_I(\{a\})`$ are calculated as convolution integrals of the relevant NLO QCD matrix elements and the universal parton distributions $`\{f(x,Q;\{a\})\}`$ for all $`Q`$. The latter are obtained from the initial functions $`\{f(x,Q_0;\{a\})\}`$ by NLO QCD evolution.

The *global analysis* consists of a systematic way to determine the best values for the $`\{a\},`$ and the associated uncertainties, by fitting $`\{T(a)\}`$ to $`\{D\}`$. Because of the wide range of experimental and theoretical sources of uncertainty mentioned in the Introduction, there are a variety of strategies to deal with the complex issues involved . In the next two sections, the primary tool we employ is conventional $`\chi ^2`$ analysis. The important task is to define an effective $`\chi ^2`$ function, called $`\chi _{\mathrm{global}}^2(a)`$, that conveniently combines the theoretical and global experimental inputs, as well as relevant physics considerations based on prior knowledge, to give an overall measure of the goodness-of-fit for a given set of PDF parameters.

Experience in global analysis of PDFs during the past two decades has demonstrated that the PDFs obtained by the minimization of such a suitably chosen $`\chi _{\mathrm{global}}^2`$ provide very useful up-to-date hadron structure functions which, although not unique, are representative of good fits between theory and experiments. Now we must quantify the uncertainties of the PDFs and their predictions; i.e., we must expand the scope of the work from merely identifying typical solutions to systematically mapping the PDF parameter space in the neighborhood around the minimum of $`\chi ^2`$.

The simplest possible choice for the $`\chi ^2`$ function would be

$$\chi ^2(a)=\underset{I=1}{\overset{N}{}}\frac{\left[D_IT_I(a)\right]^2}{\sigma _I^2}$$
(1)
where $`\sigma _I`$ is the error associated with data point $`I`$. Through $`T_I(a)`$, $`\chi ^2(a)`$ is a function of the theory parameters $`\{a\}`$. Minimization of $`\chi ^2(a)`$ would identify parameter values for which the theory fits the data. However, the simple form (1) is appropriate only for the ideal case of a uniform data set with uncorrelated errors. For data used in the global analysis, most experiments combine various systematic errors into one effective error for each data point, along with the statistical error. Then, in addition, the fully correlated normalization error of the experiment is usually specified separately. For this reason, it is natural to adopt the following definition for the effective $`\chi ^2`$ (as done in previous CTEQ analyses):

$`\chi _{\mathrm{global}}^2(a)`$ $`=`$ $`{\displaystyle \underset{n}{}}w_n\chi _n^2(a)(n\text{labels the different experiments})`$ (2)
$`\chi _n^2(a)`$ $`=`$ $`\left({\displaystyle \frac{1𝒩_n}{\sigma _n^N}}\right)^2+{\displaystyle \underset{I}{}}\left({\displaystyle \frac{𝒩_nD_{nI}T_{nI}(a)}{\sigma _{nI}^D}}\right)^2`$ (3)
For the $`n^{\mathrm{th}}`$ experiment, $`D_{nI}`$, $`\sigma _{nI}^D`$, and $`T_{nI}(a)`$ denote the data value, measurement uncertainty (statistical and systematic combined), and theoretical value (dependent on $`\{a\}`$) for the $`I^{\mathrm{th}}`$ data point; $`\sigma _n^N`$ is the experimental normalization uncertainty; $`𝒩_n`$ is an overall normalization factor (with default value $`1`$) for the data of experiment $`n`$. The factor $`w_n`$ is a possible weighting factor (with default value $`1`$) which may be necessary to take into account prior knowledge based on physics considerations or other information. The a priori choices represented by the $`w_n`$ values are present, explicitly or implicitly, in any data analysis. For instance, data inclusion or omission (choices which vary for different global analysis efforts) represent extreme cases, assigning either $`100\%`$ or $`0\%`$ weight to each available experimental data set. Similarly, choices of various elements of the analysis procedure itself represent subjective input. Subjectivity of this kind also enters into the analysis of systematic errors in experiments.

The function $`\chi _{\mathrm{global}}^2(a)`$ allows the inclusion of all experimental constraints in a uniform manner while allowing flexibility for incorporating other relevant physics input. We will make use of this function to explore the neighborhood of the best fit, and to generate sample PDFs pertinent to the uncertainty of the prediction of a specific physical variable of interest. However, the numerical value of this effective $`\chi ^2`$ function should not be given an a priori statistical interpretation, because correlations between measurement errors, and correlated theoretical errors, are not included in its definition. In particular, the likelihood of a candidate PDF set $`\{a\}`$ cannot be determined by the value of the increase $`\mathrm{\Delta }\chi _{\mathrm{global}}^2(a)`$ above the minimum.bbbThe often quoted theorem of Gaussian error analysis, that an increase of $`\chi ^2`$ by 1 unit in a constrained fit to data corresponds to 1 standard deviation of the constrained variable, is true only in the absence of correlations. When existing correlations are left out, the relevant size of $`\mathrm{\Delta }\chi ^2`$ can be much larger than $`1`$. Appendix A discusses this point in some detail. Instead, the evaluation of likelihoods and estimation of global uncertainty will be carried out in a separate step in Sec. 4, after sets of optimal sample PDFs for the physical variable of interest have been obtained.

## 3 The Lagrange Multiplier Method

The Lagrange Multiplier method is an extension of the $`\chi ^2`$ minimization procedure, that relates the range of variation of a physical observable $`X`$ dependent upon the PDFs, to the variation of the function $`\chi _{\mathrm{global}}^2(a)`$ that is used to judge the goodness of fit of the PDFs to the experimental data and PQCD.

### 3.1 The Method

The method has been introduced in . The starting point is to perform a global analysis as described in Sec. 2, by minimizing the function $`\chi _{\mathrm{global}}^2(a)`$ defined by Eq. (2), thus generating a set of PDFs that represents the best estimate consistent with current experiment and theory. We call this set the “standard set”cccThis standard set is very similar to the published CTEQ5M1 set ., denoted $`S_0`$. The parameter values that characterize this set will be denoted by $`\{a^{(0)}\}\{a_i^{(0)};i=1,\mathrm{},d\}`$; and the absolute minimum of $`\chi _{\mathrm{global}}^2`$ will be denoted by $`\chi _0^2`$. Now, let $`X`$ be a particular physical quantity of interest. It depends on the PDFs, $`X=X(a)`$, and the best estimate (or prediction) of $`X`$ is $`X_0=X(a^{(0)})`$. We will assess the uncertainty of this predicted value by a two-step analysis. First, we use the Lagrange Multiplier method to determine how the minimum of $`\chi _{\mathrm{global}}^2(a)`$ increases, i.e., how the quality of the fit to the global data set decreases, as $`X`$ deviates from the best estimate $`X_0`$. Second, in Section 4, we analyze the appropriate tolerance of $`\chi _{\mathrm{global}}^2`$.

As explained in , the first step is taken by introducing a Lagrange multiplier variable $`\lambda `$, and minimizing the function

$$\mathrm{\Psi }(\lambda ,a)=\chi _{\mathrm{global}}^2(a)+\lambda X(a)$$
(4)
with respect to the original $`d`$ parameters $`\{a\}`$ for fixed values of $`\lambda `$. In practice we minimize $`\mathrm{\Psi }(\lambda ,a)`$ for many values of the Lagrange multiplier: $`\lambda _1,\lambda _2,\mathrm{},\lambda _M`$. For each specific value $`\lambda _\alpha `$, the minimum of $`\mathrm{\Psi }(\lambda _\alpha ,a)`$ yields a set of parameters $`\{a_{\mathrm{min}}(\lambda _\alpha )\}`$, for which we evaluate the observable $`X`$ and the related $`\chi _{\mathrm{global}}^2`$. We use the shorthand $`(X_\alpha ,\chi _{\mathrm{global},\alpha }^2)`$ for this pair. $`\chi _{\mathrm{global},\alpha }^2`$ represents the lowest achievable $`\chi _{\mathrm{global}}^2`$, for the global data, for which $`X`$ has the value $`X_\alpha `$, taking into account all possible PDFs in the full $`d`$-dimensional parameter space of points $`\{a\}`$. In other words, the result $`\{a_{\mathrm{min}}(\lambda _\alpha )\}`$ is a constrained fit—with $`X`$ constrained to be $`X_\alpha `$. We can equivalently say that $`X_\alpha `$ is an extremum of $`X`$ if $`\chi _{\mathrm{global}}^2`$ is constrained to be $`\chi _{\mathrm{global},\alpha }^2`$. We denote the resulting set of PDFs by $`S_\alpha `$.

We repeat the calculation for many values of $`\lambda `$, following the chain

$$\lambda _\alpha \mathrm{min}\left[\mathrm{\Psi }(\lambda _\alpha ,a)\right]a_{\mathrm{min}}(\lambda _\alpha )X_\alpha \mathrm{and}\chi _{\mathrm{global},\alpha }^2$$

for $`\alpha =1,2,3,\mathrm{},M`$. The result is a parametric relationship between $`X`$ and $`\chi _{\mathrm{global}}^2`$, through $`\lambda `$. We call this function $`\chi _{\mathrm{global}}^2(X)`$; so $`\chi _{\mathrm{global}}^2(X_\alpha )=\chi _{\mathrm{global},\alpha }^2`$ is the minimum of $`\chi _{\mathrm{global}}^2(a)`$ when $`X`$ is constrained to be $`X_\alpha `$. The absolute minimum of $`\chi _{\mathrm{global}}^2`$, which we denote $`\chi _0^2`$, is the minimum of $`\mathrm{\Psi }(\lambda =0,a)`$, occurring at $`\{a\}=\{a^{(0)}\}`$. Thus the procedure generates a set of optimized sample PDFs along the curve of maximum variation of the physical variable $`X`$ in the $`d`$-dimensional PDF parameter space (with $`d=16`$ in our case). These PDF sets $`\{S_\alpha \}`$ are exactly what is needed to assess the range of variation of $`X`$ allowed by the data. In other words, the Lagrange Multiplier method provides optimal PDFs tailored to the physics problem at hand, in contrast to an alternative method  that generates a large sample of PDFs by the Monte Carlo method. The underlying ideas of these two complementary approaches are illustrated in the plot on the left side of Fig. 1.

$`\chi _{\mathrm{global}}^2(X)`$ is the lowest achievable value of $`\chi _{\mathrm{global}}^2(a)`$ for the value $`X`$ of the observable, where $`\chi _{\mathrm{global}}^2(a)`$ represents our measure of the goodness-of-fit to the global data. Therefore the allowed range of $`X`$, say from $`X_0\mathrm{\Delta }X`$ to $`X_0+\mathrm{\Delta }X`$, corresponding to a chosen tolerance of the goodness of fit $`\mathrm{\Delta }\chi _{\mathrm{global}}^2=\chi _{\mathrm{global}}^2\chi _0^2`$, can be determined by examining a graph of $`\chi _{\mathrm{global}}^2`$ versus $`X`$, as illustrated in the plot on the right side of Fig. 1. This method for calculating $`\mathrm{\Delta }X`$ may be more robust and reliable than the traditional error propagation because it does not approximate $`X(a)`$ and $`\chi _{\mathrm{global}}^2(a)`$ by linear and quadratic dependence on $`\{a\}`$, respectively, around the minimum.

Although the parameters $`\{a\}`$ do not appear explicitly in this analysis, the results do depend, in principle, on the choice of parameter space (including the dimension, $`d`$) in which the minimization takes place. In practice, if the degrees of freedom represented by the parametrization are chosen to match the constraining power of the global data sets used, which must be true for a sensible global analysis, the results are quite stable with respect to changes in the parametrization choices. The sensitivity to these choices is tested, as part of the continuing effort to improve the global analysis.

The discussion so far has left open this question: What is the appropriate tolerance $`\mathrm{\Delta }\chi _{\mathrm{global}}^2`$ to define the “error” of the prediction $`X_0`$? This question will be addressed in Sec. 4.

Our method can obviously be generalized to study the uncertainties of a collection of physical observables $`(X_1,X_2,\mathrm{},X_s)`$ by introducing a separate Lagrange multiplier for each observable. Although the principle stays the same, the amount of computational work increases dramatically with each additional observable.

### 3.2 A case study: the $`W`$ cross section

In this subsection we examine the cross section $`\sigma _W`$ for inclusive $`W^\pm `$ production at the Tevatron ($`p\overline{p}`$ collisions at $`\sqrt{s}=1.8`$ TeV) to illustrate the method and to lay the ground work for the quantitative study of uncertainties to be given in Sec. 4. Other examples will be described in Sec. 5. Preliminary results of this section have been reported previously .

Until recently the only method for assessing the uncertainty of $`\sigma _W`$ due to PDFs has been to compare the calculated values obtained from a number of different PDFs, as illustrated in Fig. 2, in which the plots are taken from existing literature.dddThese plots show the product of $`\sigma _W`$ times a leptonic branching ratio, which is what is measured experimentally. The branching ratio $`B`$ has some experimental error. For studying the uncertainties of $`\sigma _W`$, we will focus on $`\sigma _W`$ itself in the rest of the paper. The PDFs used in these comparisons are either the “best fits” from different global analysis groups  (hence are not pertinent to uncertainty studies) or are chosen by some simple intuitive criteria . The meaning and reliability of the resulting range of $`\sigma _W`$ are not at all clear. Furthermore, these results do not provide any leads on how the uncertainties can be improved in the future. The Lagrange Multiplier technique provides a systematic method to address and remedy both of these problems.

Let the physical quantity $`X`$ of the last subsection be the cross section $`\sigma _W`$ for $`W^\pm `$ production at the Tevatron. Applying the Lagrange method, we obtain the constrained minimum of $`\chi _{\mathrm{global}}^2`$ as a function of $`\sigma _W`$, shown as solid points in Fig. 3. The best estimate value, i.e., the prediction for the standard set $`S_0`$, is $`\sigma _{W0}=21.75`$ nb.

The curve is a polynomial fit to the points to provide a smooth representation of the continuous function $`\chi _{\mathrm{global}}^2(X)`$. We see that all the sample PDF sets obtained by this method lie on a smooth quasi-parabolic curve with the best-fit value at the minimum.

As discussed earlier (in Fig. 1) points on the curve represent our sample of optimal PDFs relevant to the determination of the uncertainty of $`\sigma _W`$. To quantify this uncertainty, we need to reach beyond the effective $`\chi _{\mathrm{global}}`$ function, and establish the confidence levels for these “alternative hypotheses” with respect to the experimental data sets used in the global analysis.

## 4 Quantifying the Uncertainty

Consider a series of sample PDF sets along the curve $`\chi _{\mathrm{global}}^2(X)`$ of Fig. 3 denoted by $`\{S_\alpha ;\alpha =0,1,\mathrm{},M\}`$ where $`S_0`$ is the standard set. These represent “alternative hypotheses” for the true PDFs, and we wish to evaluate the likelihoods associated with these alternatives. To do so, we go back to the individual experiments and, in each case, perform as detailed a statistical analysis as is permitted with available information from that experiment. After we have obtained meaningful estimates of the “errors” of these candidate PDFs with respect to the individual experiments, we shall try to combine this information into a global uncertainty measure in the form of $`\mathrm{\Delta }X`$ and $`\mathrm{\Delta }\chi _{\mathrm{global}}^2`$.

The experimental data sets included in our global analysis are listed in Table 1. For some of these experiments, information on correlated systematic errors is available (albeit usually in unpublished form). For these, statistical inference should be drawn from a more accurate $`\chi _n^2`$ function than the simple formula Eq. (3) used for the global fit. In particular, if $`\sigma _{nI}`$ is the uncorrelated error and $`\{\beta _{kI};k=1,2,\mathrm{},K\}`$ are the coefficients of $`K`$ distinct correlated errors associated with the data point $`I`$, then an appropriate formula for the $`\chi _n^2`$ function is

$$\chi _n^2=\underset{I}{}\frac{(D_{nI}T_{nI})^2}{\sigma _{nI}^2}\underset{k=1}{\overset{K}{}}\underset{k^{}=1}{\overset{K}{}}B_k\left(A^1\right)_{kk^{}}B_k^{}$$
(5)
where $`B_k`$ is a vector, and $`A_{kk^{}}`$ a matrix, in $`K`$ dimensions:

$$B_k=\underset{I}{}\beta _{kI}(D_{nI}T_{nI})/\sigma _{nI}^2;A_{kk^{}}=\delta _{kk^{}}+\underset{I}{}\beta _{kI}\beta _{k^{}I}/\sigma _{nI}^2.$$
(6)
(The sum over $`I`$ here includes only the data from experiment $`n`$.) Traditionally, $`\chi _n^2`$ is written in other ways, e.g., in terms of the inverse of the ($`N\times N`$) variance matrix. For experiments with many data points, the inversion of such large matrices may lead to numerical instabilities, in addition to being time-consuming. Our formula (5) has a significant advantage in that all the systematic errors are first combined (“analytically”) in the definitions of $`B_k`$ and $`A_{kk^{}}`$. Equation (5) requires only the inverse of the much smaller ($`K\times K`$) matrix $`A_{kk^{}}`$. ($`K`$ is the number of distinct systematic errors.) The derivation of these formulas is given in Appendix B. Equation (5) reduces to the minimum of $`\chi _n^2`$ in Eq. (3) with respect to $`𝒩_n`$ if the only correlated error is the overall normalization error for the entire data set; in that case $`\beta _I=\sigma _n^ND_{nI}`$.

By using Eq. (5), or Eq. (3) for cases where the correlations of systematic errors are unavailable, we obtain the best estimate on the range of uncertainty permitted by available information on each individual experiment. We should note that the experimental data sets are continuously evolving. Some data sets in Table 1 will soon be updated (Zeus, H1) or replaced (CCFR).eeeCf. Talks presented by these collaborations at DIS2000 *Workshop on Deep Inelastic Scattering and Related Topics*, Liverpool, England, April 2000. In addition, most information on correlated systematic errors is either unpublished or preliminary. The results presented in the following analysis should therefore be considered more as a demonstration of principle—as the first application of our proposed method—rather than the final word on the PDF uncertainty of the $`W`$ cross section.

### 4.1 Uncertainty with respect to individual experiments

As an example, we begin by comparing the $`\{S_\alpha \}`$ series for $`\sigma _W`$ at the Tevatron to the H1 data set . Results on correlated systematic errors are available for this data set,fffThese systematic errors are unpublished results, but are made available to the public on the H1 Web page. For convenience, we have approximated each of the pair of 4 non-symmetrical errors by a single symmetric error. The size of the resulting error on $`\sigma _W`$ inferred from this evaluation is not affected by that approximation. and are incorporated in the calculation using Eq. (5). The number of data points in this set is $`N_{H1}=172`$. The calculated values of $`\chi _{H1}^2/N_{H1}`$ are plotted against $`\sigma _W`$ in Fig. 4. The curve is a smooth interpolation of the points. The value of $`\chi _{H1}^2/N_{H1}`$ for the standard set $`S_0`$ (indicated by a short arrow on the plot) is $`0.975`$; and it is $`0.970`$ at the minimum of the curve. These values are quite normal for data with accurately determined measurement errors. We can therefore apply standard statistics to calculate the $`90`$% confidence level on $`\chi ^2/N`$ for $`N=172`$. The result is shown as the dashed horizontal line in Fig. 4.

We have similarly calculated $`\chi _n^2/N_n`$ including information on the correlations of systematic errors for the BCDMSp data set. The results are similar to the H1 results, except that the absolute values are all larger than 1.12, a large value for $`N=168`$ data points. This is a familiar problem in data analysis, and it is encountered in several other data sets in this global analysis (cf. below). The $`\chi _n^2/N_n`$ calculation including correlations of the errors is also done for the D0 and CDF jet cross sections.gggThe measurement errors of the jet cross sections are dominated by systematic errors, so the error correlation matrices are used for $`\chi _n^2`$ of these experiments even in $`\chi _{\mathrm{global}}^2`$. For those experiments that have only provided (effective) uncorrelated errors, we must rely on Eq. (3) for our error calculation, since that represents the best information available.

In order to obtain usable likelihood estimates from all the data sets, one must address the problem mentioned in the previous paragraph: Even in a “best fit”, the values of $`\chi ^2`$ per data point, $`\chi _n^2/N_n,`$ for individual experiments vary considerably among the established experiments (labeled by $`n`$). Specifically, $`\chi _n^2/N_n`$ ranges from $`1.51.7`$ (for ZEUS and CDFjet) on the high end to $`0.50.7`$ (for some D-Y experiments) on the low end in all good fits. Considering the fact that some of these data sets contain close to $`200`$ points, the range of variation is huge from the viewpoint of normal statistics: Experiments with $`\chi _n^2`$ $`/N_n`$ deviating from $`1.0`$ by a few times $`\sqrt{2/N_n}`$ in either direction would have to be ruled out as extremely unlikely .

The reasons for $`\chi _n^2/N_n`$ to deviate from $`1.0`$ in real experiments are complex, and vary among experiments. They are, almost by definition, not understood, since otherwise the errors would have been corrected and the resulting $`\chi ^2`$ would become consistent with the expected statistical value. Under these circumstances, a commonly adopted pragmatic approach is to focus on the relative $`\chi ^2`$ values with respect to some appropriate reference $`\chi ^2`$.hhhThe alternative is to take the *absolute* values of $`\chi _n^2`$ seriously, and hence only work with alternative hypotheses and experiments that are both self-consistent (i.e., have $`|\chi _n^2/N_n1|\sqrt{2/N_n}`$) and mutually compatible in the strict statistical sense (i.e., have overlapping likelihood functions). Since few of the precision DIS experiments are compatible in this sense, one must then abandon global analysis, and work instead with several distinct (and mutually exclusive) analyses based on different experiments. Accordingly, in the context of performing global QCD analysis, we adopt the following procedure. For each experiment (labeled by $`n`$):

* (i) Let $`\chi _{n,0}^2`$ denote the value of $`\chi _n^2`$ for the standard set $`S_0`$. We assume $`S_0`$ is a viable reference set. Because $`\chi _{n,0}^2`$ may be far from a likely value for random errors, we rescale the values of $`\chi _{n,\alpha }^2`$ (for $`\alpha =0,1,2,\mathrm{},M`$) by a factor $`C_{n0}`$, calling the result $`\overline{\chi }_{n,\alpha }^2`$

$$\chi _{n,\alpha }^2\overline{\chi }_{n,\alpha }^2C_{n0}\chi _{n,\alpha }^2.$$
(7)
The constant $`C_{n0}`$ is chosen such that, for the standard set, $`\overline{\chi }_{n,0}^2`$ assumes the most probable value for a chi-squared variable: $`\overline{\chi }_{n,0}^2=\xi _{50}`$ the 50th percentile of the chi-squared distribution $`P(\chi ^2,N_n)`$ with $`N_n`$ degrees of freedom, defined by

$$_0^{\xi _{50}}P(\chi ^2,N_n)𝑑\chi ^2=0.50.$$
(8)
(If $`N_n`$ is large then $`\xi _{50}N_n`$.) The rescaling constant $`C_{n0}`$ is thus $`\xi _{50}/\chi _{n,0}^2`$. For random errors the probability that $`\chi ^2<\xi _{50}`$ (or $`>\xi _{50}`$) is 50%. For those experiments whose $`\chi _{n,0}^2`$ deviates significantly from $`\xi _{50}`$, this rescaling procedure is meant to provide a simple (but crude) way to correct for the unknown correlations or unusual fluctuations.
* (ii) We then examine the values of $`\overline{\chi }_{n,\alpha }^2`$ for the alternative sets $`S_\alpha `$ with $`\alpha =1,2,\mathrm{},M`$, using $`\overline{\chi }_{n,\alpha }^2\overline{\chi }_{n,0}^2`$ to compute the statistical likelihood of the alternative hypothesis $`S_\alpha `$ with respect to the data set $`n`$, based on the chi-squared distribution with $`N_n`$ data points.

This procedure does not affect the results presented earlier for the H1 experiment, since $`\chi _{n,0}^2/N_n`$ is already very close to $`1`$ for that experiment.

Before presenting the results of the likelihood calculation, it is interesting to examine, in Fig. 5, the differences $`\mathrm{\Delta }\chi _{n,\alpha }^2=\chi _{n,\alpha }^2\chi _{n,0}^2`$ (before rescaling) versus $`\sigma _W`$ for the 15 data sets. (N.B. The vertical scales of the various plots are not the same, due to the large variations in the value of $`\mathrm{\Delta }\chi _{n,\alpha }^2`$ for different experiments.) The ordering of the experiments in Fig. 5 is the same as in Table 1, with experiments ordered by process (DIS, DY, $`W`$ and jet production). It is clear from these graphs that the DIS experiments place the strongest constraints on $`\sigma _W`$, because they have the largest $`\mathrm{\Delta }\chi _n^2`$ for the same $`\mathrm{\Delta }\sigma _W`$. This is to be expected since quark-antiquark annihilation makes the dominant contribution to $`\sigma _W`$. We also observe that most experiments place some constraint on $`\sigma _W`$ on both sides, but a few bound it on one side only. Globally, as shown in Fig. 3, the combined constraints give rise to a classic parabolic behavior for $`\chi _{\mathrm{global}}^2(\sigma _W)`$.

To estimate the statistical significance of the individual $`\chi _n^2`$ increases, we assume that the rescaled variable $`\overline{\chi }_n^2`$ obeys a chi-squared distribution $`P(\chi ^2,N_n)`$ for $`N_n`$ data points. Thereby, we estimate the value of $`\overline{\chi }_n^2`$ that corresponds to the 90% confidence level (CL) uncertainty for $`\sigma _W`$ (with respect to experiment $`n`$) from the formula $`\overline{\chi }_n^2=\xi _{90}`$, where $`\xi _{90}`$ is the 90th percentile defined by

$$_0^{\xi _{90}}P(\chi ^2,N_n)𝑑\chi ^2=0.90.$$
(9)
For example, Fig. 6 shows the chi-squared distribution $`P(\chi ^2,N_n)`$ for $`N_n=172`$, the number of data points in the H1 data set. The 50th and 90th percentiles are indicated. We choose a conservative 90% CL because there are other theoretical and phenomenological uncertainties not taken into account by this analysis.

To summarize our procedure, an alternative PDF set $`S_\alpha `$ lies within the 90% CL for experiment $`n`$ if it has $`\overline{\chi }_{n,\alpha }^2<\xi _{90}`$; that is, if

$$\frac{\chi _{n,\alpha }^2}{\chi _{n,0}^2}<\frac{\xi _{90}}{\xi _{50}}.$$
(10)
We judge the likelihood of $`S_\alpha `$ from the ratio of $`\chi _{n,\alpha }^2`$ to the reference value $`\chi _{n,0}^2`$, rather than from the absolute magnitude. The horizontal lines in Fig. 5 correspond to the values of $`\mathrm{\Delta }\chi _n^2`$ obtained in this way. Finally, from the intercepts of the line with the interpolating curve in each plot in Fig. 5, we obtain an estimated uncertainty range of $`\sigma _W`$ from each individual experiment. The results are presented collectively in Fig. 7, where, for each experiment, the point ($``$) is the value of $`\sigma _W`$ for which $`\chi _n^2`$ is minimum, and the error bar extends across the 90% CL based on that data set.

The uncertainty ranges shown in Fig. 7 with respect to individual experiments represent the most definitive results of our study, in the sense that the input and the assumptions can be stated clearly and the analysis is quantitative within the stated framework. It is natural to proceed further and estimate a global measure of $`\mathrm{\Delta }\sigma _W`$ and the corresponding $`\mathrm{\Delta }\chi _{\mathrm{global}}^2`$. This last step is, however, less well-defined and requires some subjective judgement.

### 4.2 The Global Uncertainty

It should be emphasized that the ranges shown by the error bars in Fig. 7 are not errors determined independently by each experiment; rather they represent the ranges allowed by the experiments for alternative global fits $`\{S_\alpha \}`$. For this reason, and others related to the rescaling of $`\chi ^2`$ mentioned earlier as well as approximations inherent in many of the original published errors,iiiFor instance, the single uncorrelated systematic error associated with each data point, which is the only systematic error given for most experimental data sets, is clearly only an “effective uncorrelated error” which qualitatively represents the effects of the many sources of systematic error, some of which are really correlated. it is not obvious how to combine these errors. We refer to the ranges in Fig. 7 by the generic term local (i.e., single-experiment) uncertainties. On a qualitative level, Fig. 7 exhibits the same features seen earlier in Fig. 5: (i) the quark dominated DIS experiments give the smallest error bars; and (ii) a few experiments only set bounds on one side, while the rest limit the range in both directions. In addition, Fig. 7 gives us an overall view which clearly shows that $`\sigma _W`$ is well constrained in the global analysis, and the experimental bounds are consistent with each other.

The important question is how to provide a sensible measure of the overall uncertainty in view of the complexity of the problem already described. The situation here is not unlike the problem of assigning an overall systematic error to an experimental measurement. Figure 7 shows a set of 90% CL ranges for $`\sigma _W`$ from different sources, but these ranges are highly correlated, because the alternative hypotheses being tested come from global fits. The final uncertainty must be a reasonable intersection of these ranges.

We will state an algorithm for obtaining the final uncertainty measure of $`\sigma _W`$ based on Fig. 7. The same algorithm can be applied in the future for predictions of other observables. It has two parts: (1) Determine the central value using all the experiments; that is the solid line in Fig. 7. (2) Then take the intersection of the error ranges as the combined uncertainty. But in calculating the intersection, experiments below the mean are used only for setting the lower bound, and experiments above the mean are used only for setting the upper bound. With this algorithm, experiments that permit a large range of $`\sigma _W`$, i.e., that depend on aspects of the PDFs that are not sensitive to the value of $`\sigma _W`$, will not affect the final uncertainty measure (as they should not). According to this algorithm, the result for the uncertainty of $`\sigma _W`$ is $`20.9\mathrm{nb}<\sigma _W<22.6\mathrm{nb}`$. These bounds are approximately $`\pm 4`$% deviations from the prediction (21.75 nb) and so we quote a $`\pm 4`$% uncertainty in $`\sigma _W`$ due to PDFs.

Now we may determine the increase in $`\chi _{\mathrm{global}}^2`$ that corresponds to our estimated uncertainty $`\mathrm{\Delta }\sigma _W`$ in the $`\sigma _W`$ prediction. Referring to Fig. 3, a deviation of $`\sigma _W`$ by $`\pm 4`$% from the minimum corresponds to an increase $`\mathrm{\Delta }\chi _{\mathrm{global}}^2180`$. That is, $`\mathrm{\Delta }\chi _{\mathrm{global}}^2`$ in Fig. 1 is 180. In other words, along the direction of maximum variation of $`\sigma _W`$ a PDF set with $`\mathrm{\Delta }\chi _{\mathrm{global}}^2180`$ is found to violate some experimental constraints by this analysis.

### 4.3 Comments

We should point out that the above uncertainty estimate, $`\mathrm{\Delta }\sigma _W/\sigma _W4`$%, represents only a lower bound on the true uncertainty, since many other sources of error have not yet been included in the analysis: theoretical ones such as QCD higher order and resummation effects, power-law corrections, and nuclear corrections. These need to be taken into consideration in a full investigation of the uncertainties, but that goes beyond the scope of this paper.jjjBecause there are these additional sources of uncertainty, we have used 90% CL’s, rather than 68% CL’s, to calculate the error. We shall add only two remarks which are more directly related to our analysis.

The first concerns a technical detail. In the results reported so far, we have fixed the normalization factors {$`𝒩_n`$} in the definition of $`\chi _{\mathrm{global}}^2`$ (Eq. (2)) at their values determined in the standard fit $`S_0`$. If we let these factors float when we perform the Lagrange Multiplier analysis, $`\mathrm{\Delta }\sigma _W`$ will increase noticeably compared to Fig. 3 for the same $`\mathrm{\Delta }\chi _{\mathrm{global}}^2`$. However, upon closer examination, this behavior can be easily understood and it does not imply a real increase in the uncertainty of $`\sigma _W`$. The key observation is that the additional increase (or decrease) in $`\sigma _W`$ is entirely due to a *uniform* increase (or decrease) of {$`𝒩_n`$} for all the DIS experiments. There is a simple reason for this: The values of the $`q`$ and $`\overline{q}`$ distributions in the relevant $`x`$ range (which determine the value of $`\sigma _W`$) are approximately proportional to {$`𝒩_n`$}DIS. Although every experiment does have a normalization uncertainty, the probability that the normalization factors of all the *independent* DIS experiments would shift in the *same* direction by the *same* amount is certainly unlikely. Hence we avoid this artificial effect by fixing {$`𝒩_n`$} at their “best values” for our study. Allowing the factors {$`𝒩_n`$} to vary randomly (within the published experimental normalization uncertainties) would not change our estimated value of $`\mathrm{\Delta }\sigma _W`$ significantly.

The second remark concerns the choice of parametrization. We have mentioned that even the robust Lagrange Multiplier method depends in principle on the choice of the parton parameter space, i.e., on the choice of the functional forms used for the nonperturbative PDFs at the low momentum scale $`Q_0`$. To check how our answers depend on the choice of parametrization in practice, we have done many similar calculations, using different numbers of free parameters within the same functional form (cf. Appendix C), and using different functional forms for the factor multiplying $`x^a\left(1x\right)^b`$. We have not seen any dependence of the uncertainty estimates on these changes. Although more radical ways of parametrizing the nonperturbative PDFs might affect the result more, there is no known example of such a parametrization, which at the same time still provides an equally good fit to the full data set.

## 5 Further Examples

### 5.1 $`W^\pm `$ production at the LHC

A study similar to the last section has been carried out for inclusive $`W^\pm `$ production at the LHC. Figure 8

shows $`\chi _{\mathrm{global}}^2`$ versus $`\sigma _W`$ for the process $`ppW^\pm X`$ at $`\sqrt{s}=14`$ TeV, summed over the two final states. The curve is a smooth interpolation of a series of PDF sets $`\{S_\alpha \}`$ generated by the Lagrange Multiplier method. The best estimate value of the LHC cross section is $`\sigma _W=189.7`$ nb.

Comparing Figs. 3 and 8, one immediately notices that the uncertainty of $`\sigma _W`$(LHC) is greater than that of $`\sigma _W`$(Tevatron) for the same $`\mathrm{\Delta }\chi _{\mathrm{global}}^2`$. This indicates that because $`W`$ production in $`pp`$ collisions at the LHC and $`\overline{p}p`$ collisions at the Tevatron involve different mixtures of parton subprocesses as well as different kinematic ranges, the constraints imposed by current experiments included in the global analysis are also different for the two cases. Referring to the map of the $`d`$-dimensional PDF parameter space on the left side of Fig. 1, we are generating sample PDFs along different directions $`L_X`$ in the two cases. Therefore it is not surprising that the rate of variation of $`\chi _{\mathrm{global}}^2`$ is also different.

To demonstrate this point, and to quantify the uncertainty on the LHC prediction, we have carried out the same error analysis as in Sec. 4, i.e., comparing the alternative PDFs to the individual experiments. Figure 9

gives the final overview of the 90% CL ranges of $`\sigma _W`$ obtained from these comparisons, analogous to Fig. 7 for the Tevatron cross section. There are some differences compared to the Tevatron case. The LHC prediction is more tightly constrained by experiments that are sensitive to PDFs at small $`x`$. This makes sense, because $`W`$ production at the LHC is not dominated by valence $`q\overline{q}`$ interactions. We note in particular that the two inclusive jet production experiments place significant constraints on $`\sigma _W`$ at the LHC.

We can combine the individual error bars in Fig. 9 according to the algorithm proposed in Sec. 4 to produce a global uncertainty measure for $`\sigma _W`$ at the LHC. The lower bound on $`\sigma _W`$(LHC) obtained by the intersections of the individual ranges is $`\sigma _W=175.3`$ nb; the upper bound is $`\sigma _W=204.6`$ nb. These bounds, shown as the dashed lines in Fig. 9, correspond to $`\pm 8`$% deviations from the prediction (189.7 nb). The global uncertainty on $`\sigma _W`$(LHC) is thus significantly larger than that on $`\sigma _W`$(Tevatron). Reinforcing this conclusion is the fact that the scatter of the points in Fig. 9 is larger than in Fig. 7.

We can again inspect the increase in $`\chi _{\mathrm{global}}^2`$ from the quoted range of alternative fits. A $`\pm 8`$% deviation from the minimum, symmetrized for simplicity, corresponds to the increase $`\mathrm{\Delta }\chi _{\mathrm{global}}^2200`$. This number is similar to the increase in $`\chi _{\mathrm{global}}^2`$ for our estimated uncertainty of $`\sigma _W`$ at the Tevatron.

In the companion paper , we make some process-independent estimates of $`\mathrm{\Delta }\chi _{\mathrm{global}}^2`$ based on completely different considerations. Those arguments also yield the same order-of-magnitude estimates of $`\mathrm{\Delta }\chi _{\mathrm{global}}^2`$ (in the range from 100 to 200) for acceptable PDFs around the global minimum. Since the effective $`\chi _{\mathrm{global}}^2`$, as a measure of goodness-of-fit, does not have a normal statistical implication, points on a constant $`\chi _{\mathrm{global}}^2`$ surface in the PDF parameter space do not necessarily correspond to a constant likelihood. Some variation with the direction in the multi-dimensional space is to be expected.

### 5.2 Uncertainties on $`Z^0`$ production

We conclude this section by presenting results from applying the Lagrange Multiplier method to $`Z^0`$ production at the Tevatron and the LHC.

Figure 10a shows the minimum $`\chi _{\mathrm{global}}^2`$ as a function of $`\sigma _Z`$ at the Tevatron. The global prediction is $`\sigma _Z=6.55`$ nb. The experimental measurement by the D0 collaboration is $`\sigma _ZB=0.221\pm 0.003\pm 0.011`$ nb; the result from CDF (all data from Run I) is $`\sigma _ZB=0.250\pm 0.004\pm 0.010`$ nb. (Here $`B`$ is the branching ratio for $`Z^0e\overline{e}`$, which is $`(3.367\pm 0.005)\times 10^2`$.) The comparison of the prediction to Tevatron data is discussed below. Analyzing the local $`\chi _n^2`$ in the manner of Sec. 4, in order to quantify the uncertainty of the prediction, we find that the uncertainty of $`\sigma _Z`$(Tevatron) due to PDFs is $`\pm 3`$% of the prediction. The corresponding increase in $`\chi _{\mathrm{global}}^2`$, symmetrized for simplicity, is approximately 130.

For the LHC process $`ppZ^0X`$, Fig. 10b shows the minimum $`\chi _{\mathrm{global}}^2`$ as a function of $`\sigma _Z`$. The dependence of $`\chi _{\mathrm{global}}^2`$ on $`\sigma _{W/Z}`$(LHC) exhibits a behavior departing from quadratic over the full range of $`\sigma _{W/Z}`$ under study. This is evidence that the Lagrange multiplier method can go beyond the traditional error matrix approach (which depends on the quadratic approximation) in exploring the neighborhood of the minimum.

The global prediction is $`\sigma _Z=58.0`$ nb. Analyzing the local $`\chi _n^2`$ as in the other cases, we find that the uncertainty of $`\sigma _Z`$(LHC) due to PDFs is approximately $`\pm 10`$%. As in the case of $`W^\pm `$ production, the PDF uncertainty for $`Z^0`$ production at the LHC is significantly larger than that at the Tevatron. Measurement of $`W^\pm `$ and $`Z^0`$ production at the LHC will therefore provide significant information on PDFs.

### 5.3 Comparison with existing data

For $`W`$ and $`Z`$ production at the Tevatron, we can compare our calculated cross sections $`\sigma _W`$ and $`\sigma _Z`$, with their ranges of uncertainty $`\pm 4`$% and $`\pm 3`$% respectively, to the measurements of CDF and D0 from Run I.  The comparison is shown in Fig. 11. The two experiments do not measure $`\sigma _W`$ and $`\sigma _Z`$ per se, but rather $`\sigma _WB_W`$ and $`\sigma _ZB_Z`$ where $`B_W`$ is the branching ratio for $`W^{}e\overline{\nu }`$ and $`B_Z`$ is the branching ratio for $`Z^0e\overline{e}`$. We have used the values $`B_W=0.106`$ and $`B_Z=0.0337`$ for the calculations . The bands in Fig. 11 show the ranges of $`\sigma _WB_W`$ and $`\sigma _ZB_Z`$ from our PDF uncertainty study (but no uncertainty included from $`B_W`$ and $`B_Z`$). The two measurements of $`\sigma _WB_W`$ are consistent with the uncertainty range. The two measurements of $`\sigma _ZB_Z`$ are not.

It should be noted that CDF and D0 use different normalizations for their luminosity determinations. The CDF collaboration bases its luminosity purely on its own measurement of the inelastic $`\overline{p}p`$ cross section , while D0 uses the world average for this cross section. Thus current luminosities quoted by CDF are $`6.2`$% lower than those quoted by D0. Consequently, all CDF cross section measurements are ab initio $`6.2`$% higher than those of D0. If the CDF/D0 measurements of $`\sigma _WB_W`$ and $`\sigma _ZB_Z`$ are rescaled by $`6.2`$% with respect to each other, they are in excellent agreement.

Because of the uncertainty in the inelastic $`\overline{p}p`$ cross section, it has been proposed to normalize future Tevatron (and LHC) physics cross sections to the measured $`W`$ cross section (or rather $`\sigma _WB_W`$). This makes the determination of the uncertainty of $`\sigma _W`$ due to PDFs even more important.

### 5.4 Comparison with previous uncertainty estimates

It is interesting to contrast our results to existing estimates of the uncertainties of $`\sigma _W`$ and $`\sigma _Z`$ at the Tevatron and LHC colliders based on the traditional method of comparing results obtained from somewhat ad hoc PDFs. Some of these previous comparisons for $`\sigma _W`$(Tevatron) between historical PDFs as well as various trial up/down sets obtained by the CTEQ and MRST groups were shown in Fig. 2. We will briefly comment on the results of  in the context of this paper.

Reference  constructs an extended set of MRST PDFs, of which the most important for $`\sigma _W`$ and $`\sigma _Z`$ are the standard set MRST99 and three pairs of up/down sets designated

$$\{\alpha _S,\alpha _S\},\{g,g\},\{q,q\}$$

in which some aspect of the parton distributions is either raised ($``$) or lowered ($``$) by an amount that represents an educated guess of a “standard deviation”. The predictions of $`\sigma _W`$ and $`\sigma _Z`$ are then compared for these alternative PDF sets to get an idea of the uncertainty due to PDFs.

In the case of the Tevatron processes, the deviations of $`\sigma _W`$ or $`\sigma _Z`$ from the value for MRST99 for sets $`\{\alpha _S,\alpha _S\},\{g,g\},\{q,q\}`$ were found to be $`\pm 2`$%, $`\pm 1`$%, $`\pm 3`$% respectively. From these results, the authors of  concluded that the uncertainties of $`\sigma _W`$ and $`\sigma _Z`$ at the Tevatron are no more than about $`\pm 4`$%, and mainly attributable to the normalization uncertainty in the input $`u`$ and $`d`$ distributions. This conclusion appears to be quite consistent with the results of the previous sections based on exploring the variation of the cross section over the entire PDF parameter space. (This range of uncertainty also happens to coincide with what one would get by comparing historical PDF sets, as shown in the right plot of Fig. 2.)

For the LHC, the MRST study found that the uncertainty of $`\sigma _W`$ and $`\sigma _Z`$ at the LHC is only slightly larger than at the Tevatron; the uncertainty was estimated to be $`\pm 5`$%. The largest observed variations came from the sets $`\alpha _S`$ and $`\alpha _S`$, differing from the standard prediction by $`\pm 4`$$`5`$%. This estimate is *considerably lower* than the $`\pm 8`$$`10`$% result obtained by the detailed analysis of the previous sections. We have verified that the PDF sets that give $`\pm 8`$$`10`$% deviations of $`\sigma _W`$(LHC) and $`\sigma _Z`$(LHC) from the standard prediction (represented by the points at the outer edges of the corresponding plots in Figs. 8 and 10) provide equally good or better fits to the global data sets compared to the fits of {$`\alpha _S`$, $`\alpha _S`$}.kkkThe values of $`\chi _{\mathrm{global}}^2`$ for {$`\alpha _S`$, $`\alpha _S`$} are {$`1531,\mathrm{\hspace{0.33em}1356}`$} compared to $`1400`$ for the outermost LHC sets shown in Figs. 8 and 10 and presented in the next section.,lllThe global fits used at first in our exploration of PDF uncertainty were conducted with fixed $`\alpha _s`$. To make sure that this restriction does not result in an underestimate of the uncertainties of $`\sigma _W`$ and $`\sigma _Z`$, we have examined the effect of freeing $`\alpha _s`$ in the analysis (but imposing the known constraints from the world average of $`\alpha _s`$). The results on the size of the uncertainties are not changed noticeably. This is because the full variations in the PDFs (particularly the gluon) allowed in the Lagrange approach can absorb the added degree of freedom. Thus, it is clear that the Lagrange Multiplier method can generate optimal PDFs, i.e., having the largest excursion of the variable $`X`$ of interest, which are difficult to discover by ad hoc trial-and-error methods used in the past.

## 6 PDF Sets for exploring W and Z Physics

The parton distribution sets used in the above calculations are useful for exploring some aspects of $`W`$ and $`Z`$ physics at the Tevatron and LHC, since they provide much more reliable estimates on the PDF uncertainties than existing ones in the literature, which are not designed to probe the full range of possibilities in the parton parameter space. With this in mind, we present in this section some representative PDF sets for applications to the rate of $`W`$ and $`Z`$ production at the Tevatron and LHC. These PDFs are relevant to the total cross sections $`\sigma _W`$ and $`\sigma _Z`$, and each corresponds to a particular direction ($`L_X`$) in the PDF parameter space (see Fig. 1). Therefore they are not suitable for estimating the PDF uncertainties of other observables that are sensitive to other aspects of the PDFs. Other PDF sets can be obtained, using the method introduced in this paper, to probe the range of other variables, such as rapidity ($`y`$) or transverse momentum ($`p_T`$) distributions (hence relevant to the measurement of $`W`$ mass). These will be investigated in subsequent work. Also, the companion paper  supplies information from the Hessian method that can be used to construct the optimal PDFs for any observable $`X`$.

The PDF set that yields the “best estimate” for all of the physical cross sections covered in this paper is our standard set $`S_0`$. The parametrization of the initial distribution is given in Appendix C. In the following, we present two sets of PDFs that bound the likely range for each of the cross sections.

To exemplify the PDFs that characterize the range of uncertainty of $`W`$ production at the Tevatron, we use two representative sets, labeled $`S_{W,\mathrm{TeV}}^\pm `$, which correspond to $`\sigma _W=`$ $`\sigma _W(S_0)\pm \mathrm{\Delta }\sigma _W`$ (with $`\mathrm{\Delta }\sigma _W/\sigma _W0.04`$) respectively. These two sets are extreme fits obtained by the Lagrange Multiplier method. The parameters $`\{a\}`$ for these sets are given in Appendix C. We now compare some of the parton distributions from the three sets ($`S_{W,\mathrm{TeV}}^{},S_{0,}S_{W,\mathrm{TeV}}^+`$), to examine the ranges of variation of the PDFs themselves.

Figure 12a shows $`u(x,Q)`$, $`d(x,Q)`$ and $`g(x,Q)`$ for $`S_{W,\mathrm{TeV}}^\pm `$, compared to the standard set $`S_0`$, at $`Q=80`$ GeV. The function $`x\mathrm{\Delta }f(x)`$ is plotted for each parton flavor, where $`\mathrm{\Delta }f`$ is $`ff_0`$. The gluon function has been divided by 10 to fit on the same graph. The solid curves ($`\mathrm{\Delta }u_+`$, $`\mathrm{\Delta }d_+`$, $`\mathrm{\Delta }g_+`$) correspond to $`S_{W,\mathrm{TeV}}^+`$, and the dashed curves to $`S_{W,\mathrm{TeV}}^{}`$. For $`S_{W,\mathrm{TeV}}^+`$, requiring $`\sigma _W`$ to be larger than $`\sigma _W^{(0)}`$ makes the $`u`$ and $`d`$ distributions larger than for the standard fit ($`u_0`$ and $`d_0`$) so $`\mathrm{\Delta }u_+`$ and $`\mathrm{\Delta }d_+`$ are positive. Then the gluon distribution must be smaller than the standard because of the momentum sum rule. In the case of $`S_{W,\mathrm{TeV}}^{}`$, the reverse is true, resulting in almost a mirror behavior. At the Tevatron, a typical $`x`$ for the parton-level process $`q_1\overline{q}_2W^\pm `$ is $`M_W/\sqrt{s}=0.04`$. The differences $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ are significant in the range $`0.01x0.04`$. The magnitude of $`\mathrm{\Delta }f(x)`$ in this range is a few percent of the standard $`f_0(x)`$, which makes sense since $`\mathrm{\Delta }\sigma _W`$ is a $`4`$% shift of $`\sigma _W`$ for these PDF sets.

We can carry out the same comparison for $`W`$ production at the LHC. The PDFs that bound the range of uncertainty are designated as $`S_{W,\mathrm{LHC}}^\pm `$, which correspond to $`\sigma _W=\sigma _W(S_0)\pm \mathrm{\Delta }\sigma _W`$ with $`\mathrm{\Delta }\sigma _W/\sigma _W0.08`$. The PDF parameters are given in Appendix C. Figure 12b shows parton distributions from $`S_{W,\mathrm{LHC}}^\pm `$. (Again, the gluon has been divided by 10.) In the LHC case, the typical $`x`$ for the process $`q_1\overline{q}_2W^\pm `$ is $`M_W/\sqrt{s}=0.006`$. The region where $`\mathrm{\Delta }u`$ and $`\mathrm{\Delta }d`$ are significant is seen accordingly lower in $`x`$ than for the Tevatron case.

In Appendix C we also present PDF sets $`S_{Z,\mathrm{TeV}}^\pm `$ and $`S_{Z,\mathrm{LHC}}^\pm `$ that characterize the range of uncertainties of $`Z`$ production at the Tevatron and LHC. These correspond to the outlying points on Figs. 10a,b. They are similar to $`S_{W,\mathrm{TeV}}^\pm `$ and $`S_{W,\mathrm{LHC}}^\pm `$, with small differences in the flavor dependence.

## 7 Summary

We have developed the Lagrange Multiplier method to calculate uncertainties of physical observables due to PDFs, and we have used the method to determine the uncertainty of the total cross sections for $`W^\pm `$ production and $`Z^0`$ production at the Tevatron and LHC. The method is more reliable than past estimates because: (i) it explores all the possibilities in the parameter space of the input PDFs, independent of other assumptions; and (ii) it produces the maximum allowed range for the specified physical variables. This is in contrast to previous attempts which relied on varying certain features of the parton distributions chosen in some ad hoc way.

From this analysis, we find that the uncertainty of the prediction for $`\sigma _W`$ or $`\sigma _Z`$ at the Tevatron with current experimental constraints is approximately $`\pm 3`$$`4`$%; and at the LHC the uncertainties are approximately $`\pm 8`$$`10`$%. These numbers do not include other uncertainties associated with theoretical approximations, nuclear corrections, and other unexpected sources. We have explored to some extent the possible effects due to the choice of parametrization of the nonperturbative input PDFs, and found them to be small. The current work should be considered exploratory in nature, as a first application of this improved approach to error estimates. A more comprehensive study, based on soon to be improved data sets, and including other sources of uncertainties, will produce better overall estimates of the physical predictions.

This study should be regarded as the precursor for many interesting applications to come, on physical processes of interest to the precision study of the Standard Model, and on predictions for New Physics at future colliders. Some examples are rapidity distributions of $`W^\pm `$ and $`Z^0`$ production, which contain a wealth of information on parton structure of the nucleon; the $`W`$ mass measurement; top and Higgs cross sections, etc.

There are other approaches to error estimates in global QCD analysis . In general, if greater emphasis is placed on the “rigor” of the statistical method, then the range of experiments that can be included in the analysis is narrower. We have chosen to emphasize the inclusion of the full range of experimental constraints, and adapt the statistical analysis to deal with the practical problems that one faces in such a system. Within our general framework, there is an alternative, complementary approach based on the conventional error matrix method . We explore this latter method, as applied to global QCD analysis of PDFs, in a companion paper . We mention briefly the contrasting features and relative merits of the two approaches here.

The Lagrange Multiplier method focuses on a given physical observable $`X`$ (or a set of observables $`\{X_k\}`$) and determines the uncertainty $`\mathrm{\Delta }X`$ allowed by the global data set within a specified tolerance for the global fit. The error matrix approach, using the Hessian matrix, focuses instead on the uncertainties of the PDFs as represented by the parameters $`\{a_i;i=1,\mathrm{},d\}`$. It is in principle universal because, once determined, these errors can be propagated to *any* physical variable $`X`$. However, the results are reliable only if the function $`\chi _{\mathrm{global}}^2(a)`$ and the observable $`X(a)`$ can be approximated by lowest order expansions in the parameters $`\{a\}`$, and if the numerical computation of the derivatives (the Hessian matrix) is under control. The latter problem is surprisingly difficult for global QCD analysis, because the eigenvalues of the error matrix vary by many orders of magnitude. This problem has been solved , and the error matrix results are consistent with the constrained fitting results . Thus, at present, both methods appear to be applicable to the study of uncertainties in global QCD analysis.

In Figures 8 and 10b there is a significant cubic term in the dependence of $`\chi _{\mathrm{global}}^2`$ on $`\sigma _W`$(LHC) and $`\sigma _Z`$(LHC), respectively. To calculate $`\chi _{\mathrm{global}}^2`$ versus $`X`$ accurately in such cases, the Lagrange Multiplier method is necessary. Traditional linear error analysis based on the Hessian matrix can only produce a quadratic approximation to the dependence.

When both methods are applicable, the Hessian method is more flexible and easier to apply computationally. But generally the Lagrange method is more robust and reliable. As we expand the investigation to other physical processes of interest, we will continue to test the efficacy of both methods and cross check the results.

#### Acknowledgement

This work was supported in part by NSF grant PHY-9802564. We thank Michiel Botje and Fabian Zomer for several interesting discussions and for valuable comments on our work. We also thank John Collins, Dave Soper, and other CTEQ colleagues for discussions.

## Appendix A The effect of correlated errors on $`𝚫\chi ^\mathrm{𝟐}`$

The global fitting function $`\chi _{\mathrm{global}}^2`$ defined in (2) resembles the standard statistical variable $`\chi ^2`$, so it is tempting to try to apply theorems of Gaussian statistics to analyze the significance of the fit between theory and experiment. However, the familiar theorems do not apply, because of correlations between measurement errors. The purpose of this Appendix is to explore this issue. The effect of correlated errors is potentially a source of confusion.

For simplicity we describe the simplest case: measurement of a single observable. The arguments can be extended to cases where multiple quantities are measured, such as the determination of parton distribution functions.

Consider an observable $`m`$ that is measured $`N`$ times. We shall refer to $`N`$ measurements of $`m`$ as one “experiment”. Let the true value of $`m`$ be $`m_0`$. The measurements are $`m_1,m_2,m_3,\mathrm{},m_N`$. The deviations from the true value are $`\alpha _1,\alpha _2,\alpha _3,\mathrm{},\alpha _N`$, where $`\alpha _i=m_im_0`$. In general the measurement errors are correlated, so in the Gaussian approximation the probability distribution of the fluctuations is

$$dP=𝒩\mathrm{exp}\left\{\frac{1}{2}\underset{i,j=1}{\overset{N}{}}\alpha _iC_{ij}\alpha _j\right\}d^N\alpha .$$
(11)
Here $`C_{ij}`$ is a real symmetric matrix, and $`𝒩=\sqrt{\mathrm{Det}C}/(2\pi )^{N/2}`$ ensures the normalization condition $`𝑑P=1`$.

We will need the variance matrix $`\alpha _i\alpha _j`$, where the notation $`Q`$ means the average of $`Q`$ in the probability distribution (11). For this Gaussian distribution,

$$\alpha _i\alpha _j=\left(C^1\right)_{ij}.$$
(12)
The mean square fluctuation $`E_i`$ of the $`i^{\mathrm{th}}`$ measurement $`m_i`$ is

$$E_i\alpha _i^2=\left(C^1\right)_{ii}.$$
(13)
To find the best estimate of the value of $`m`$ from these $`N`$ measurements, ignoring the correlations in the measurement errors, we define a chi-squared function $`\chi _u^2(m)`$ by

$$\chi _u^2(m)=\underset{i=1}{\overset{N}{}}\frac{\left(m_im\right)^2}{E_i}.$$
(14)
The value of $`m`$ that minimizes $`\chi _u^2(m)`$, call it $`\overline{m}`$, is then the best estimate of $`m_0`$ based on this information. The function $`\chi _u^2(m)`$ is analogous to the fitting function $`\chi _{\mathrm{global}}^2`$ in the CTEQ program, in the sense that it does not include information about the correlations between errors. The minimum of $`\chi _u^2(m)`$ occurs at a weighted average of the measurements,

$$\overline{m}=\frac{_{i=1}^Nm_i/E_i}{_{i=1}^N1/E_i}.$$
(15)
If all the $`E_i`$’s are equal then $`\overline{m}`$ is just the average of the measurements.

Now, what are the fluctuations of the mean $`\overline{m}`$? That is, if the “experiment” consisting of $`N`$ measurements could be replicated many times, what would be the distribution of $`\overline{m}`$’s obtained in those many trials? It turns out that $`\overline{m}`$ has a Gaussian distribution

$$\frac{dP}{d\overline{m}}=\frac{1}{\sqrt{2\pi \mathrm{\Sigma }^2}}\mathrm{exp}\left[(\overline{m}m_0)^2/(2\mathrm{\Sigma }^2)\right].$$
(16)
The standard deviation $`\mathrm{\Sigma }`$ of $`\overline{m}`$ is the RMS fluctuation; that is,

$$\mathrm{\Sigma }^2=\left(\overline{m}m\right)^2𝑑P=\frac{1}{D^2}\underset{ij}{}\frac{\left(C^1\right)_{ij}}{E_iE_j}$$
(17)
where

$$D=\underset{i}{}\frac{1}{E_i}.$$
(18)
The question we wish to answer is this: How much does $`\chi _u^2(m)`$ increase, when $`m`$ moves away from the minimum (at $`\overline{m}`$) by the amount $`\pm \mathrm{\Sigma }`$ that corresponds to one standard deviation of the mean? The answer to this question is

$$\mathrm{\Delta }\chi _u^2=\mathrm{\Sigma }^2D.$$
(19)
This result follows easily from the definition (14), because

$$\chi _u^2(\overline{m}+\mathrm{\Sigma })\chi _u^2(\overline{m})=2\mathrm{\Sigma }\underset{i}{}\frac{m_i\overline{m}}{E_i}+\mathrm{\Sigma }^2\underset{i}{}\frac{1}{E_i},$$
(20)
and the term linear in $`\mathrm{\Sigma }`$ is $`0`$ by the definition of $`\overline{m}`$. So far the discussion has been quite general. We will now examine some illustrative special cases.

Example 1: Suppose the measurement errors are uncorrelated; that is,

$$C_{ij}=\delta _{ij}/E_i.$$
(21)
Then the standard deviation of the mean $`\overline{m}`$ is $`\mathrm{\Sigma }=1/\sqrt{D}`$. Thus for the uncorrelated case, the increase of $`\chi _u^2`$ corresponding to one standard deviation of the mean is $`\mathrm{\Delta }\chi _u^2=1`$. This is the “normal” statistical result: The $`1\sigma `$ range corresponds to an increase of $`\chi ^2`$ by 1.

An even more special case is when the errors are uncorrelated and constant: $`E_i=\sigma ^2`$ independent of $`i`$, where $`\sigma `$ is the standard deviation of single measurements. The correlation matrix is $`C_{ij}=\delta _{ij}/\sigma ^2`$. In this case $`D`$ is $`N/\sigma ^2`$, and the standard deviation of the mean is $`\mathrm{\Sigma }=\sigma /\sqrt{N}`$.

The criterion $`\mathrm{\Delta }\chi ^2=1`$ for one standard deviation of a measured quantity is a standard result, often used in the analysis of precision data. But if $`\chi ^2`$ is defined ignoring the correlations between measurement errors, then the criterion $`\mathrm{\Delta }\chi ^2=1`$ is only valid for uncorrelated errors. We will next consider two examples with correlated errors, to show that $`\mathrm{\Delta }\chi _u^2`$ is not $`1`$ for such cases.

Example 2: Suppose measurements 1 and 2 are correlated, 3 and 4 are correlated, 5 and 6 are correlated, etc. Then the correlation matrix is

$$C_{ij}=\{\begin{array}{c}1/\sigma ^2\mathrm{for}i=j\hfill \\ c/\sigma ^2\mathrm{for}ij=12\mathrm{or}21,34\mathrm{or}43,\mathrm{𝑒𝑡𝑐}\hfill \\ 0\mathrm{otherwise}\hfill \end{array}$$
(22)
where $`1<c<1`$ since the determinant of $`C`$ must be positive. The inverse matrix $`C^1`$ can be constructed using the fact that $`C`$ is block diagonal, consisting of $`N/2`$ $`2\times 2`$ blocks. Then it can be shown that

$$\mathrm{\Sigma }=\frac{\sigma }{\sqrt{N}\sqrt{1+c}}\mathrm{and}\mathrm{\Delta }\chi _u^2=1c.$$
(23)
The increase of $`\chi _u^2`$ for one standard deviation of the mean ranges from $`0`$ to $`2`$, depending on $`c`$. The criterion $`\mathrm{\Delta }\chi ^2=1`$ does not apply to this example with correlated errors. A standard increase of $`\chi _u^2`$ may be smaller or larger than $`1`$.

Example 3: For an even more striking example, suppose the $`N`$ measurements that constitute a single “experiment” are, for $`i=1,2,3,\mathrm{},N`$,

$$m_i=m_0+y_i+\beta $$
(24)
where the $`y_i`$ are randomly distributed with standard deviation $`\sigma `$, and the measurements are systematically off by the amount $`\beta `$. Suppose that $`\beta `$ has a Gaussian distribution with standard deviation $`s`$ for replications of the “experiment”. In this example,

$`C_{ij}`$ $`=`$ $`{\displaystyle \frac{1}{\sigma ^2}}\left(\delta _{ij}{\displaystyle \frac{s^2}{Ns^2+\sigma ^2}}\right),`$ (25)
$`(C^1)_{ij}`$ $`=`$ $`\sigma ^2\delta _{ij}+s^2.`$ (26)
The variance of the individual measurements ($`m_i`$) is

$$m^2m^2=\sigma ^2+s^2.$$
(27)
Therefore our uncorrelated chi-squared variable $`\chi _u^2(m)`$, defined ignoring the correlations, is

$$\chi _u^2(m)=\underset{i=1}{\overset{N}{}}\frac{\left(mm_i\right)^2}{\sigma ^2+s^2}.$$
(28)
The minimum of $`\chi _u^2(m)`$ occurs at $`\overline{m}`$, which is just the average of the individual measurements. The variance of $`\overline{m}`$, averaged over many replications of the “experiment”, is

$$\mathrm{\Sigma }^2=\overline{m}^2\overline{m}^2=s^2+\frac{\sigma ^2}{N}.$$
(29)
The increase of $`\chi _u^2`$ as $`m`$ moves from $`\overline{m}`$ to $`\overline{m}\pm \mathrm{\Sigma }`$, i.e., by one standard deviation of the mean, is

$$\mathrm{\Delta }\chi _u^2\chi _u^2(\overline{m}+\mathrm{\Sigma })\chi _u^2(\overline{m})=\frac{\sigma ^2+Ns^2}{\sigma ^2+s^2}.$$
(30)
In the limit $`s/\sigma 1`$, the error correlations in this model become negligible and $`\mathrm{\Delta }\chi ^2`$ reduces to the conventional value of $`1`$. But in the limit $`s/\sigma 1`$ where the error correlations are dominant, $`\mathrm{\Delta }\chi ^2`$ approaches $`N`$.

Thus for Example 3—a systematic error with 100% correlation between measurements—the increase of $`\chi _u^2`$ for a standard deviation of $`\overline{m}`$ is much larger than 1. If $`s`$ and $`\sigma `$ are comparable, then $`\mathrm{\Delta }\chi _u^2`$ is of order $`N`$.

If the correlation matrix $`C_{ij}`$ is known accurately, then the correlation information can be incorporated into the definition of the $`\chi ^2`$ function, in the manner of Appendix B. For the full list of experiments in the global analysis of parton distribution functions, however, the correlations of systematic errors have not been published, so the fitting function $`\chi _{\mathrm{global}}^2`$ has only uncorrelated systematic errors.

We described above the measurement of a single quantity. The determination of parton distribution functions seeks to measure many quantities, i.e., the $`16`$ parameters $`\{a\}`$. The above arguments can be extended to measurements of multiple quantities. If the measurement errors are uncorrelated, then the increase of $`\chi _u^2`$ by 1 from the minimum defines a hyperellipse in parameter space—the error ellipse—corresponding to one standard deviation of linear combinations of the parameters. However, if the errors are correlated then $`\mathrm{\Delta }\chi _u^2=1`$ is not the correct criterion for a standard deviation.

The Lagrange Multiplier method finds the best fit to the data, subject to a constrained value of some quantity $`X`$. The prediction of $`X`$ is at the absolute minimum of $`\chi _u^2`$. Again, if the errors are uncorrelated then one standard deviation of the constrained quantity corresponds to an increase of $`\chi ^2`$ by 1 from the absolute minimum. But if the errors are correlated then $`\mathrm{\Delta }\chi _u^2=1`$ is not the correct criterion for one standard deviation of $`X`$.

One reason for describing this familiar, even elementary, statistics, is to avoid certain misconceptions. Our standard PDF set $`S_0`$ is a parametrized fit to 1295 data points with 16 fitting parameters. The minimum value of $`\chi _{\mathrm{global}}^2`$ is approximately 1200. Naively, it seems that an increase of $`\chi _{\mathrm{global}}^2`$ by merely 1, say from 1200 to 1201, could not possibly represent a standard deviation of the fit. Naively one might suppose that a standard deviation would have $`\mathrm{\Delta }\chi ^2\sqrt{1295}`$ rather than 1. However, this is a misconception. If the errors are uncorrelated (or if the correlations are incorporated into $`\chi ^2`$) then indeed $`\mathrm{\Delta }\chi ^2=1`$ would represent a standard deviation. But this theorem is irrelevant to our problem, because the large correlations of systematic errors are not taken into account in $`\chi _{\mathrm{global}}^2`$.

## Appendix B $`\chi ^\mathrm{𝟐}`$ function including correlated systematic errors

The purpose of this appendix is to derive the appropriate definition of $`\chi ^2`$ for data with correlated systematic errors. The defining condition is that $`\chi ^2`$ should obey a chi-squared distribution.

Let $`\left\{m_i\right\}`$ be a set of measurements, where $`i=1,2,3,\mathrm{},N`$. Let $`t_i`$ be the true, i.e., theoretical value of the $`i^{\mathrm{th}}`$ measured quantity. Several kinds of measurement errors will contribute to the difference between $`m_i`$ and $`t_i`$. The uncorrelated error of measurement $`i`$ is denoted by $`\sigma _i`$. There are also correlated errors, $`K`$ in number, denoted $`\beta _{1i},\beta _{2i},\mathrm{},\beta _{Ki}`$. Thus the $`i^{\mathrm{th}}`$ measurement can be written as

$$m_i=t_i+\mathrm{errors}=t_i+\sigma _ir_i+\underset{j=1}{\overset{K}{}}\beta _{ji}r_j^{}$$
(31)
where $`r_i`$ and $`r_j^{}`$ are independently fluctuating variables. We assume that each of these fluctuations has a Gaussian distribution with width $`1`$,

$$p(r)=\frac{e^{r^2/2}}{\sqrt{2\pi }}.$$
(32)
Note that $`r_j^{}`$ is independent of $`i`$; that is, the errors $`\beta _{j1},\beta _{j2},\mathrm{},\beta _{jN}`$ are 100% correlated for all $`N`$ data points.

The probability distribution of the measurements is

$`dP`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}p(r_i)dr_i\underset{j=1}{\overset{K}{}}p(r_j^{})dr_j^{}}`$ (33)
$`\times `$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\delta \left(m_it_i\sigma _ir_i{\displaystyle \underset{j=1}{\overset{K}{}}}\beta _{ji}r_j^{}\right)d^Nm.`$
Now we will evaluate the integrals over $`r_i`$ and $`r_j^{}`$ in two steps. First evaluate the $`r_i`$ integrals using the delta functions,

$$dP=\underset{j=1}{\overset{K}{}}dr_j^{}𝒞_1e^{\chi _1^2/2}d^Nm$$
(34)
where $`𝒞_1`$ is a normalization constant and

$$\chi _1^2=\underset{i=1}{\overset{N}{}}\left(\frac{m_it_i_j\beta _{ji}r_j^{}}{\sigma _i}\right)^2+\underset{j=1}{\overset{K}{}}r_{j}^{}{}_{}{}^{2}.$$
(35)
Note that $`\chi _1^2`$ is a function of $`r_1^{},\mathrm{},r_K^{}`$. These variables $`\{r_j^{}\}`$ could be used as fitting parameters to account for the systematic errors: Minimizing $`\chi _1^2`$ with respect to $`r_j^{}`$ would provide the best model to correct for the systematic error of type $`j`$. Because $`\chi _1^2`$ is only a quadratic polynomial in the $`r_j^{}`$ variables, the minimization can be done analytically.

To continue evaluating (33) we now do the integration over $`\{r_j^{}\}`$. Write $`\chi _1^2`$ in the form

$$\chi _1^2=\underset{i=1}{\overset{N}{}}\frac{\left(m_it_i\right)^2}{\sigma _i^2}\underset{j=1}{\overset{K}{}}2B_jr_j^{}+\underset{j,j^{}=1}{\overset{K}{}}A_{jj^{}}r_j^{}r_j^{}^{}$$
(36)
where $`B_j`$ is a vector with $`K`$ components

$$B_j=\underset{i=1}{\overset{N}{}}\beta _{ji}(m_it_i)/\sigma _i^2,$$
(37)
and $`A_{jj^{}}`$ is a $`K\times K`$ matrix

$$A_{jj^{}}=\delta _{jj^{}}+\underset{i=1}{\overset{N}{}}\beta _{ji}\beta _{j^{}i}/\sigma _i^2.$$
(38)
Then the integration over $`d^Kr^{}`$ is an exercise in Gaussian integration, with the result

$$dP=𝒞\mathrm{exp}\left[\frac{1}{2}\chi ^2\right]d^Nm$$
(39)
where $`𝒞`$ is a normalization constant and

$$\chi ^2=\underset{i=1}{\overset{N}{}}\frac{(m_it_i)^2}{\sigma _i^2}\underset{j=1}{\overset{K}{}}\underset{j^{}=1}{\overset{K}{}}B_j\left(A^1\right)_{jj^{}}B_j^{}.$$
(40)
This equation is the appropriate definition of $`\chi ^2`$ for data with correlated systematic errors. The correlated errors are defined by the coefficients $`\beta _{ji}`$ in (31), which determine the vector $`B_j`$ and matrix $`A_{jj^{}}`$. An interesting relation is that the $`\chi ^2`$ quantity in (40) is the minimum of $`\chi _1^2`$ with respect to the parameters $`r_1^{},\mathrm{},r_K^{}`$.

Another expression for $`\chi ^2`$, which may be derived from (33) by Gaussian integration, is

$$\chi ^2=\underset{i=1}{\overset{N}{}}\underset{i^{}=1}{\overset{N}{}}\left(m_it_i\right)\left(V^1\right)_{ii^{}}\left(m_i^{}t_i^{}\right)$$
(41)
where $`V_{ij}`$ is the variance matrix

$$V_{ii^{}}=\sigma _i^2\delta _{ii^{}}+\underset{j=1}{\overset{K}{}}\beta _{ji}\beta _{ji^{}}.$$
(42)
It can be shown that the inverse of the variance matrix is

$$\left(V^1\right)_{ii^{}}=\frac{\delta _{ii^{}}}{\sigma _i^2}\underset{j,j^{}=1}{\overset{K}{}}\frac{\beta _{ji}\beta _{j^{}i^{}}}{\sigma _{i}^{}{}_{}{}^{2}\sigma _i^{}^2}\left(A^1\right)_{jj^{}}.$$
(43)
Therefore (40) and (41) are equivalent. However, there is a real computational advantage in the use of (41) because it does not require the numerical inversion of the $`N\times N`$ variance matrix.

To check that (40) makes sense we can consider a special case. Suppose the number $`K`$ of systematic errors is $`N`$, and each systematic error contributes to just one measurement. Then the matrix of systematic errors has the form

$$\beta _{ji}=\delta _{ji}b_i.$$
(44)
This situation is equivalent to an additional set of uncorrelated errors $`\left\{b_i\right\}`$. The vector $`B_j`$ is then

$$B_j=\frac{b_j\left(m_jt_j\right)}{\sigma _j^2}$$
(45)
and the matrix $`A_{jj^{}}`$ is

$$A_{jj^{}}=\delta _{jj^{}}\left[1+\frac{b_j^2}{\sigma _j^2}\right].$$
(46)
Substituting these results into (40) we find

$$\chi ^2=\underset{i}{}\frac{\left(m_it_i\right)^2}{\sigma _i^2+b_i^2},$$
(47)
which makes sense: the uncorrelated errors just combine in quadrature.

The statistical quantity $`\chi ^2`$ has a chi-squared distribution with $`N`$ degrees of freedom. Thus this variable may be used to set confidence levels of the theory for the given data. But to use this variable, the measurement errors $`\sigma _i`$ and $`\beta _{ji}`$, for $`i=1,2,\mathrm{},N`$ and $`j=1,2,\mathrm{},K`$, must be known from the experiment. A chi-squared distribution with many degrees of freedom is a very narrow distribution, sharply peaked at $`\chi ^2=N`$. Therefore small inaccuracies in the values of the $`\sigma _i`$’s and $`\beta _{ji}`$’s may translate into a large error on the confidence levels computed from the chi-squared distribution.

It is equation (40) that we use in Section 4 to compare the constrained fits produced by the Lagrange multiplier method to data from the H1 and BCDMS experiments. Correlated systematic errors are also used to calculate $`\chi ^2`$ for the CDF and D0 jet experiments.

## Appendix C Parton Distribution Sets

We give here the PDFs described in Section 6. $`S_0`$ is the standard set, defined by the absolute minimum of $`\chi _{\mathrm{global}}^2`$. $`S_{W,\mathrm{Tev}}^\pm `$ are fits to the global data sets with extreme values of $`\sigma _W`$(Tevatron), i.e., the outermost points on Fig. 3, generated by the Lagrange multiplier method. $`S_{Z,\mathrm{Tev}}^\pm `$, $`S_{W,\mathrm{LHC}}^\pm `$, and $`S_{Z,\mathrm{LHC}}^\pm `$ are analogous for $`Z`$ production and $`W`$ and $`Z`$ production at the LHC.

The functional form of the initial parton distributions and the definitions of the PDF parameters at the low-energy scale $`Q_0=1`$ GeV are

$$f(x,Q_0^2)=A_0x^{A_1}(1x)^{A_2}(1+A_3x^{A_4})$$

for $`u_v,d_v,g,\overline{u}+\overline{d},s(=\overline{s})`$; and for the ratio

$$\frac{\overline{d}(x,Q_0^2)}{\overline{u}(x,Q_0^2)}=A_0x^{A_1}(1x)^{A_2}+(1+A_3x)(1x)^{A_4}.$$

The tables of coefficients follow.

$$S_0:\begin{array}{cccccc}& & & & & \\ & \hfill A_0& \hfill A_1& \hfill A_2& \hfill A_3& \hfill A_4\\ & & & & & \\ d_v& \hfill 0.5959& \hfill 0.4942& \hfill 4.2785& \hfill 8.4187& \hfill 0.7867\\ & & & & & \\ u_v& \hfill 0.9783& \hfill 0.4942& \hfill 3.3705& \hfill 10.0012& \hfill 0.8571\\ & & & & & \\ g& \hfill 3.3862& \hfill 0.2610& \hfill 3.4795& \hfill 0.9653& \hfill 1.\\ & & & & & \\ \overline{d}/\overline{u}& \hfill 3.051E4& \hfill 5.4143& \hfill 15.& \hfill 9.8535& \hfill 4.3558\\ & & & & & \\ \overline{u}+\overline{d}& \hfill 0.5089& \hfill 0.0877& \hfill 7.7482& \hfill 3.3890& \hfill 1.\\ & & & & & \\ s& \hfill 0.1018& \hfill 0.0877& \hfill 7.7482& \hfill 3.3890& \hfill 1.\end{array}$$

$$S_{W,\mathrm{TeV}}^\pm :\begin{array}{cccccc}& & & & & \\ & \hfill A_0& \hfill A_1& \hfill A_2& \hfill A_3& \hfill A_4\\ & & & & & \\ d_v& \hfill 0.2891& \hfill 0.5141& \hfill 3.8555& \hfill 10.9580& \hfill 0.4128\\ & \hfill 0.2184& \hfill 0.2958& \hfill 4.6267& \hfill 35.7229& \hfill 1.0958\\ & & & & & \\ u_v& \hfill 1.0142& \hfill 0.5141& \hfill 3.3614& \hfill 9.2995& \hfill 0.8053\\ & \hfill 0.2979& \hfill 0.2958& \hfill 3.3279& \hfill 32.8453& \hfill 0.9427\\ & & & & & \\ g& \hfill 4.6245& \hfill 0.4354& \hfill 3.4795& \hfill 0.9728& \hfill 1.\\ & \hfill 1.8080& \hfill 0.0458& \hfill 3.4795& \hfill 0.0519& \hfill 1.\\ & & & & & \\ \overline{d}/\overline{u}& \hfill 5.908E4& \hfill 5.6673& \hfill 15.& \hfill 9.8535& \hfill 4.7458\\ & \hfill 2.041E4& \hfill 5.1506& \hfill 15.& \hfill 9.8535& \hfill 4.8320\\ & & & & & \\ \overline{u}+\overline{d}& \hfill 0.4615& \hfill 0.0108& \hfill 6.6145& \hfill 0.92784& \hfill 1.\\ & \hfill 1.2515& \hfill 0.3338& \hfill 7.5216& \hfill 0.0570& \hfill 1.\\ & & & & & \\ s& \hfill 0.0923& \hfill 0.0108& \hfill 6.6145& \hfill 0.9278& \hfill 1.\\ & \hfill 0.2503& \hfill 0.3338& \hfill 7.5216& \hfill 0.0570& \hfill 1.\end{array}$$

$$S_{Z,\mathrm{TeV}}^\pm :\begin{array}{cccccc}& & & & & \\ & \hfill A_0& \hfill A_1& \hfill A_2& \hfill A_3& \hfill A_4\\ & & & & & \\ d_v& \hfill 0.6061& \hfill 0.5502& \hfill 4.0017& \hfill 5.8346& \hfill 0.5343\\ & \hfill 0.3427& \hfill 0.3728& \hfill 4.5166& \hfill 19.8510& \hfill 0.9966\\ & & & & & \\ u_v& \hfill 1.2159& \hfill 0.5502& \hfill 3.3347& \hfill 7.3386& \hfill 0.7711\\ & \hfill 0.5247& \hfill 0.3728& \hfill 3.3905& \hfill 20.1006& \hfill 0.9556\\ & & & & & \\ g& \hfill 4.4962& \hfill 0.4321& \hfill 3.4795& \hfill 0.9023& \hfill 1.\\ & \hfill 2.3113& \hfill 0.1032& \hfill 3.4795& \hfill 0.6349& \hfill 1.\\ & & & & & \\ \overline{d}/\overline{u}& \hfill 4.321E4& \hfill 5.4724& \hfill 15.& \hfill 9.8535& \hfill 4.6298\\ & \hfill 2.818E4& \hfill 5.4540& \hfill 15.& \hfill 9.8535& \hfill 4.4376\\ & & & & & \\ \overline{u}+\overline{d}& \hfill 0.4609& \hfill 0.0103& \hfill 6.6671& \hfill 0.9822& \hfill 1.\\ & \hfill 0.9900& \hfill 0.2926& \hfill 8.3205& \hfill 2.1648& \hfill 1.\\ & & & & & \\ s& \hfill 0.0921& \hfill 0.0103& \hfill 6.6671& \hfill 0.9823& \hfill 1.\\ & \hfill 0.1980& \hfill 0.2926& \hfill 8.3205& \hfill 2.1648& \hfill 1.\end{array}$$

$$S_{W,\mathrm{LHC}}^\pm :\begin{array}{cccccc}& & & & & \\ & \hfill A_0& \hfill A_1& \hfill A_2& \hfill A_3& \hfill A_4\\ & & & & & \\ d_v& \hfill 0.7326& \hfill 0.5008& \hfill 4.6393& \hfill 10.8532& \hfill 1.0595\\ & \hfill 0.5671& \hfill 0.4771& \hfill 4.2615& \hfill 8.8355& \hfill 0.8130\\ & & & & & \\ u_v& \hfill 1.0608& \hfill 0.5008& \hfill 3.4023& \hfill 9.6622& \hfill 0.8968\\ & \hfill 0.9142& \hfill 0.4771& \hfill 3.3761& \hfill 10.9138& \hfill 0.8809\\ & & & & & \\ g& \hfill 2.2379& \hfill 0.0733& \hfill 3.4795& \hfill 0.9860& \hfill 1.\\ & \hfill 2.5021& \hfill 0.3981& \hfill 3.4795& \hfill 1.6229& \hfill 1.\\ & & & & & \\ \overline{d}/\overline{u}& \hfill 2.178E4& \hfill 5.2576& \hfill 15.& \hfill 9.8535& \hfill 4.4810\\ & \hfill 4.531E4& \hfill 5.4979& \hfill 15.& \hfill 9.8535& \hfill 4.6585\\ & & & & & \\ \overline{u}+\overline{d}& \hfill 1.1980& \hfill 0.2952& \hfill 6.9475& \hfill 0.5442& \hfill 1.\\ & \hfill 0.2759& \hfill 0.0918& \hfill 8.2045& \hfill 6.3950& \hfill 1.\\ & & & & & \\ s& \hfill 0.2396& \hfill 0.2952& \hfill 6.9475& \hfill 0.5442& \hfill 1.\\ & \hfill 0.0552& \hfill 0.0918& \hfill 8.2045& \hfill 6.3950& \hfill 1.\end{array}$$

$$S_{Z,\mathrm{LHC}}^\pm :\begin{array}{cccccc}& & & & & \\ & \hfill A_0& \hfill A_1& \hfill A_2& \hfill A_3& \hfill A_4\\ & & & & & \\ d_v& \hfill 0.5659& \hfill 0.4616& \hfill 4.5297& \hfill 12.3685& \hfill 0.9836\\ & \hfill 0.4585& \hfill 0.4496& \hfill 4.2122& \hfill 10.3850& \hfill 0.7760\\ & & & & & \\ u_v& \hfill 0.8344& \hfill 0.4616& \hfill 3.3847& \hfill 12.1129& \hfill 0.8872\\ & \hfill 0.7640& \hfill 0.4496& \hfill 3.3566& \hfill 12.8253& \hfill 0.8701\\ & & & & & \\ g& \hfill 2.3282& \hfill 0.0918& \hfill 3.4795& \hfill 0.9837& \hfill 1.\\ & \hfill 2.9475& \hfill 0.4219& \hfill 3.4795& \hfill 0.9447& \hfill 1.\\ & & & & & \\ \overline{d}/\overline{u}& \hfill 2.421E4& \hfill 5.3032& \hfill 15.& \hfill 9.8535& \hfill 4.5341\\ & \hfill 4.416E4& \hfill 5.4708& \hfill 15.& \hfill 9.8535& \hfill 4.7925\\ & & & & & \\ \overline{u}+\overline{d}& \hfill 1.1130& \hfill 0.2698& \hfill 6.8490& \hfill 0.5330& \hfill 1.\\ & \hfill 0.2719& \hfill 0.0899& \hfill 8.1492& \hfill 6.5300& \hfill 1.\\ & & & & & \\ s& \hfill 0.2226& \hfill 0.2698& \hfill 6.8490& \hfill 0.5330& \hfill 1.\\ & \hfill 0.0544& \hfill 0.0899& \hfill 8.1492& \hfill 6.5300& \hfill 1.\end{array}$$

With a program to solve the PDF evolution equations, the PDFs for an arbitrary momentum scale $`Q`$ can be generated.