A Simultaneous Measurement of the QCD Colour Factors and the Strong Coupling

Using data from e+e- annihilation into hadrons, taken with the OPAL detector at LEP at the Z pole between 1991 and 1995, we performed a simultaneous measurement of the colour factors of the underlying gauge group of the strong interaction, CF and CA, and the strong coupling, alpha(s). The measurement was carried out by fitting next-to-leading order perturbative predictions to measured angular correlations of 4-jet events together with multi-jet related variables. Our results, CA = 3.02 +/- 0.25 (stat.) +/- 0.49 (syst.), CF = 1.34 +/- 0.13 (stat.) +/- 0.22 (syst.), alpha(s)(M_Z) = 0.120 +/- 0.011 (stat.) +/- 0.020 (syst.), provide a test of perturbative QCD in which the only assumptions are non-abelian gauge symmetry and standard hadronization models. The measurements are in agreement with SU(3) expectations for CF and CA and the world average of alpha(s)(M_Z).


Introduction
Electron-positron annihilation into hadrons at high energies provides a precise means to test Quantum Chromodynamics (QCD), the theory of strong interactions. Due to the purely leptonic initial state, there are many experimental quantities for which the long-distance (non-perturbative) effects are modest, especially when compared to similar quantities in hadron-hadron collisions or deep inelastic scattering. These quantities, for instance the total cross section and jet-related correlations, can be calculated in perturbative QCD as a function of a single parameter, the strong coupling strength, α s . Therefore, many QCD tests based on measurements of α s have been carried out at LEP, the Large Electron-Positron collider at CERN (see, for example, [1][2][3][4][5][6]).
A key ingredient of QCD is the underlying gauge group. The simultaneous measurement of the strong coupling and the eigenvalues of the quadratic Casimir operator of the underlying gauge theory, the C F and C A colour factors, provides a more general and comprehensive test of QCD than measurements of α s alone. The possible existence of a light gluino 1 influences both α s and the measured value of the colour factors (or, assuming SU(3) dynamics, the number of light fermionic degrees of freedom N f ). A simultaneous fit of these parameters to data provides a means to investigate whether the data favour or exclude the additional degrees of freedom.
Previous analyses of LEP data have been performed to probe the underlying gauge structure of the theory [9][10][11][12][13]. The leading order perturbative prediction was fitted to distributions of four-jet angular correlations (normalised to the total number of the selected four-jet events), leaving the ratios of the colour factors as free parameters. Due to the normalization, the leading order prediction of these distributions is independent of the strong coupling. However, the independence of the leading order prediction from the strong coupling does not imply that the radiative corrections are negligible. In fact, in Ref. [14] the value of T R /C F (see the Appendix) was estimated to increase about 25% if next-to-leading order (NLO) theoretical predictions are used instead of the leading order (LO) ones. It is therefore desirable to explicitly verify the effect of the next-to-leading order corrections on these normalised angular distributions.
Beyond leading order in the theoretical predictions, the strong coupling and the colour factors are interdependent. On the one hand, even the normalised next-to-leading order predictions for the angular distributions depend upon α s [14], while on the other hand, α s itself depends on the colour factors through its running. In this sense, beyond leading order, the only consistent way to determine the three free parameters is a simultaneous measurement of the strong coupling and the colour factors.
Predictions beyond leading order are made possible by theoretical developments achieved in the last few years. For multi-jet rates as well as numerous event shape distributions with perturbative expansions starting at O(α s ), matched next-to-leading order and nextto-leading logarithmic approximations provide very precise descriptions of the data over a wide range of the available kinematic region [15][16][17][18]. Also, recently, the next-to-leading order predictions for the distributions of four-jet angular correlations have been calculated [14,19]. To make use of these developments, we perform a simultaneous fit of the strong coupling and ratios of colour factors using next-to-leading order predictions of four-jet rates, differential two-jet rates and four-jet angular correlations. The data were collected by the OPAL Collaboration at LEP.
The outline of the paper is as follows: In Section 2, we present the observables used in the analysis and describe their best available perturbative predictions. The analysis procedure is explained in detail and the fit results are presented in Section 3. Section 4 contains the discussion of the systematic checks which were performed and the corresponding systematic errors. We collect our results and compare them to the results of previous studies in Section 5. Section 6 contains our conclusions.

Observables
To perform a simultaneous measurement of the strong coupling and the colour factors we use multi-jet related variables to gain sensitivity to α s and four-jet angular correlations to gain sensitivity to the colour factors.
There are many ways of defining jets. For the present analysis we use the Durham scheme [16,20]. Starting by defining each particle to be an individual jet, a resolution variable y ij is calculated for each pair of jets i and j: where E i and E j are the energies of jets i and j, θ ij is the angle between them and E vis is the sum of the energies of the visible particles in the event or of the partons in a theoretical calculation. If the smallest value of y ij is smaller then a predefined value y cut , the pair is replaced by a pseudo-jet with four momentum p µ ij = p µ i + p µ j and the clustering starts again with momenta p µ i and p µ j dropped and p µ ij added to the final state. Clustering ends when the smallest value of y ij is larger than y cut .
In our analysis we use the differential two-jet rate, where y 23 is the y cut value for which the two-and three-jet configurations are separated in a given event, and the four-jet rate, to gain sensitivity to α s . For theory, σ tot is the total hadronic cross section at O(α s ) accuracy. For experiment, σ tot is the total visible hadronic cross section. The D 2 distribution has been extensively used at LEP to determine α s [4,21]. Use of R 4 allows us to obtain sensitivity to α s using a four-jet based variable. For these observables both next-to-leading order (NLO) [17,18] and next-to-leading logarithmic (NLL) [16] perturbative predictions are known. The former provide a fixed order approximation valid for large values of the resolution variable, while the latter provide an all order summed but approximate cross section valid for small values of the resolution variable. The explicit NLO and NLL formulae and their matched expressions used in our fits are given in the Appendix. Using those formulae, we obtain various theoretical predictions for the four-jet rates (R 4 ) and the differential two-jet rates (D 2 ) which are shown in Fig. 1. In these plots, the theoretical predictions with the exception of the fitted curves were obtained using the world average of the strong coupling α s (M Z ) = 0.118 [22] and standard QCD values for the colour factors, C F = 4/3 and C A = 3. The theoretical predictions in Fig. 1 are shown in comparison to OPAL data corrected to the parton level. The manner in which the data are corrected to the parton level is explained in Section 3.3. It is seen that the theoretical predictions provide a good description of the data. For the theoretical plots in Fig. 1, we used the R-matched expression for R 4 and ln R-matched expression for D 2 (see the Appendix).
The other class of observables we employ in our measurements is the four-jet angular correlations. The angular variables are • the Bengtsson-Zerwas angle [23]: • the modified Nachtmann-Reiter angle [24]: • the Körner-Schierholtz-Willrodt angle [25]: • the angle between the two lowest energy jets [10]: where the p i denote the three-momenta of the energy-ordered jets (E 1 > E 2 > E 3 > E 4 ). We defined the jets using the Durham clustering and selected four-jet events at y cut = 0.008. In order to have a high statistics four-jet sample, y cut should be chosen in the range from 10 −2 to 10 −3 . However, we cannot choose y cut to be much smaller than 0.01 for technical reasons (see Section 3.3.1). For the four-jet events, used in obtaining the angular correlations, an additional cut was placed on the energy of each jet: E > E min = 3 GeV.
To obtain as much information as possible, we used all four angular variables in our analysis. The shape of the distributions of these angular correlations at LO and NLO are hardly distinguishable. Fig. 2 shows the normalised NLO theoretical predictions fitted to the data.
Although the various angular correlations are not entirely independent, the correlations (see Section 3.4) are not large as can be seen from Table 1, where the maximum values of the magnitudes of the bin-wise correlations between the observables are given. We also performed fits using a single angular variable together with one of the jet-related variables and we observed that these fits give very scattered central values. However, in these cases we found very strong correlations between the fitted parameters (≃ 1), indicating that the distributions do not contain sufficient information to constrain all of the QCD parameters simultaneously.

3
Analysis Procedure and Results  Table 1: Maximum values of the bin-wise correlations between the angular correlation variables.

The OPAL detector
A detailed description of the OPAL detector can be found in Ref. [26]. This analysis relies mainly on the reconstruction of charged particle trajectories and on the measurement of energy deposited in the electromagnetic and hadronic calorimeters. The central detector contains a silicon micro-vertex detector and three drift chamber devices: an inner vertex chamber, a large jet chamber and a surrounding z-chamber. The central detector is located inside a solenoidal magnet which provides a uniform axial 2 magnetic field of approximately 0.435 T. Tracking of the charged particles is performed with this detector. Most of the tracking information is obtained from the jet chamber which provides up to 159 measured space points per track, and close to 100 % tracking efficiency in the region | cos θ| < 0.98. The average angular resolution is approximately 0.1 mrad in φ and 1 mrad in θ.
Electromagnetic energy is measured by lead glass calorimeters surrounding the magnet coil, separated into a barrel (| cos θ| < 0.82) and two end-cap (0.81 < | cos θ| < 0.98) sections. The electromagnetic calorimeter consists of 11 704 lead glass blocks with a depth of 24.6 radiation lengths in the barrel and more than 22 radiation lengths in the end-caps. This calorimeter is surrounded by the hadronic calorimeter of the sampling type measuring the energy of hadrons emerging from the electromagnetic calorimeter.

Data samples
The analysis presented in this paper is based on 4.1 million hadronic events recorded within 3 GeV of the Z peak by the OPAL detector between 1991 and 1995. The descriptions of the OPAL trigger system and the offline multihadronic event selection are given in Refs. [27] and [28]. In this analysis, additional criteria were applied to eliminate poorly measured tracks and obtain well contained events. Each track was required to have (i) at least 40 measured points in the jet chamber, (ii) transverse momentum in the r-φ plane greater than 0.15 GeV/c, (iii) measured momentum less than 60 GeV/c, (iv) | cos θ| < 0.97, (v) distance of closest approach to the origin in the r-φ plane of no more than 5 cm, (vi) and less than 25 cm in the z direction. Hadronic events were required to contain at least five tracks to reduce contamination from e + e − → τ + τ − events.
Clusters in the electromagnetic calorimeter were accepted if they had more than 0.1 GeV energy in the barrel section or more than 0.25 GeV in the end-cap section. The corresponding clusters were required to span at least two lead glass blocks. Clusters in the hadronic calorimeter were required to have energies larger than 1 GeV.
For the calculation of the visible energy, E vis , the tracks are assumed to have the pion mass, and clusters are treated as photons. To determine the energy of each cluster measured in the electromagnetic and hadronic calorimeters, we used a matching algorithm [29]. This algorithm reduces double counting of energy and gives a better resolution both in energy and angle than the calorimeters alone. We did not impose energy-momentum conservation, but checked that such a constraint does not affect our results.
The event thrust axis was determined using the accepted tracks and clusters. In order that the event be well contained, we required | cos θ thrust | < 0.9.
The final event sample contained about 3.6 million hadronic Z decays from which about 250 000 were identified as four-jet events.

Data Correction
We choose as the theoretical reference the distributions given in Eqs. (22), (28) of the Appendix which are obtained from a parton level calculation. In order to compare our (detector level) data to these (parton level) theoretical predictions the measured distributions were corrected for the effects of the detector and hadronization.
First, we combined the distributions of the six variables cos χ BZ , cos Θ NR , cos Φ KSW , cos α 34 , D 2 and R 4 into a one-dimensional distribution, with 128 bins, 20 for each angular correlation and 24 each for the four-jet rates (covering the range of 0.001 < y cut < 1.0) and the differential two-jet rates (in the range of 1.0 < − ln y 23 < 5.8). Next, we applied bin-by-bin corrections for detector distortions and hadronization, as described below. The correction factors were obtained from the same distributions computed from Monte Carlo events.

Monte Carlo simulation
The correction factors were estimated using two different modes of the JETSET Monte Carlo program [30,31]: the parton shower mode and the matrix element mode. In the parton shower approach, the evolution of the parton system is treated as a branching process based on the leading logarithmic approximation starting from the original qq pair. In the matrix element mode, up to four partons are generated according to the second-order ERT calculation [32].
It is well known that the parton shower versions of the available Monte Carlo programs provide a good description of jet rates, but cannot describe the structure of the four-jet events properly [33]. Therefore, to determine correction factors for the four-jet angular correlations, we used Monte Carlo events with only four partons (qqgg and qqqq) in the final state. These events were generated using the second-order matrix element (ME) calculation, subjected to hadronization, as implemented in JETSET version 7.4 assuming standard QCD colour factors and the set of parameters ERT-MC-1 of Ref. [34]. 560 000 Monte Carlo events were generated. To avoid singularities in the four-parton generation of the LO matrix elements, an intrinsic cutoff y 0 has to be applied in the Monte Carlo simulation. In JETSET, y 0 = 0.01 is used, calculated according to the invariant mass resolution definition, which for four-jet events corresponds to a Durham resolution of about 0.004 (yielding the same four-jet rates). Our resolution criterion of y cut = 0.008 is safely larger than the intrinsic cut-off. For normalised angular correlations, the distributions at leading and next-to-leading order are almost the same [14,19]. Thus for the purpose of determining the corrections the ERT Monte Carlo events are sufficient.
For the correction of the differential two-jet rates and the four-jet rate, we used 2 million events generated with the parton shower (PS) version of JETSET, using parameters tuned to OPAL data at √ s =91 GeV described in Ref. [35]. For both models, hadronization of the parton system is based on the Lund string-fragmentation scheme [36]. The events of both types of samples were passed through the OPAL detector simulation [37], processed by the same reconstruction programs, and subjected to the same event selection criteria as the data.

Correction procedure
To correct the data for the finite acceptance and resolution of the detector, we determined bin-by-bin correction factors for every bin of the distribution from Monte Carlo samples with hadronization included, both before (hadron level, H i MC ) and after (detector level, D i MC ) the detector simulation and selection cuts. The hadron level treats all charged and neutral particles with lifetimes greater than 3 × 10 −10 s as stable and has no initial state photon radiation. The correction factor for detector effects for the i th bin of the distribution is given by: We examined the effects of initial state photon radiation by switching the radiation on in the hadron level Monte Carlo simulation. The correction factors were found to be identical within statistical uncertainties with or without photon radiation over the entire phase space.
We estimated the hadronization effects by computing the distributions of Monte Carlo events at the parton (P i MC ) and hadron levels (H i MC ) (before and after hadronization) and defined the correction factor for hadronization by: Using the two correction factors C i det and C i had , we correct the measured distribution to the parton level according to We show the distributions of the total correction factors,

Fit procedure
Having prepared the corrected distribution of six variables, we performed a χ 2 minimization to determine the most probable values of the variables η = α s C F /(2π), x = C A /C F and y = T R /C F (see Section A) with the program MINUIT [38]. The function χ 2 was defined by: where δ i is the bin-wise difference between the theoretical QCD prediction and the corrected data distribution, while σ ij is the covariance matrix taking into account the statistical error of the data, the finite Monte Carlo statistics, and the correlations between the bins of the distributions. The data and Monte Carlo event samples were each divided into 90 subsamples 3 to compute σ ij , using the standard formula: where D (n) i is the value of the corrected distribution in bin i of the nth subsample, D i is the value in that bin for the whole sample and N is the number of subsamples. We included bins in the fit only if the total correction factors were in the range 0.9 < C tot i < 1.1, except for the four-jet rate, where 0.85 < C i tot < 1.15 was used. 4 This bin choice is motivated by our expectation that the Monte Carlo estimate is less reliable if hadronization and/or detector distortions are large. The total number of bins used in the fit was 82, leaving 79 degrees of freedom since there are three fitted parameters (η, x and y, see Section 3.5).

Fit results
The central values of the fit results for the parameter η and colour factor ratios x and y are listed in the upper section of Table 2. The renormalization scale was fixed at x µ = 1. The errors and correlations (ρ) in the table are purely statistical.
We observe a strong correlation between x and y. We tested other variables (for instance, the three-jet rate) in an attempt to decrease these correlations, but did not observe a significant improvement. The correlations were taken into account when the fit results were converted to the standard QCD parameters α s (M Z ) and colour factors using the definitions in Eqs. (10) and (12) and T R = 1/2. These latter results are given in the lower section of Table 2.

Systematic uncertainties
The systematic uncertainties were evaluated by considering the following effects: • the measurement process and accuracy of the Monte Carlo detector simulation, • dependence on the model of hadronization, • theoretical uncertainties arising from the choice of the renormalization scale, • variation of the matching scheme in the case of the differential two-jet distribution, • variation of the fit range, • the background from five parton events, • variation of the parameter y cut .
In the systematic checks we always used the same covariance matrix that was optimized for the standard analysis.

Detector effects
The uncertainty related to the detector simulation and the measurement process was estimated by repeating the analysis using either tracks only or clusters only for the data and the detector level Monte Carlo samples. For the tracks and clusters we applied the same selection procedure as described in Section 3.2, except that all cluster related criteria were left out in the former case and all track related criteria were left out in the latter. Furthermore, the entire analysis was repeated using | cos θ thrust | < 0.7 (rather than | cos θ thrust | < 0.9), or N ch ≥ 8 (rather than N ch ≥ 5). The results obtained for these analyses are given in Table 3.

Hadronization models
Several combinations of different Monte Carlo models were considered to estimate the uncertainty associated with the hadronization correction. Since we used two different Monte Carlo samples to determine the correction factors (ME and PS), we examined the effects induced by variations of these models separately.
The angular correlations were corrected using the ME version of JETSET with the Lund string fragmentation scheme, as stated in Section 3.3.1. In the systematic studies we used the ME version of the HERWIG generator [39], which produces four-parton final states according to the leading order QCD matrix elements, and then a parton shower is started from each of the four partons followed by hadronization based on the cluster fragmentation. 5 We also varied the principal hadronization parameters of the JETSET ME Monte Carlo program, σ q and a. As the hadronization parameter b is strongly correlated with a we kept it fixed. We increased a or decreased σ q , alternately, by 10 %. These changes correspond to about one standard deviation in the parameter values allowed for the JETSET ME generator at Z peak, as found by the LEP experiments [40]. The fits were redone using data corrected by the factors calculated from these Monte Carlo samples at the generator level (reevaluating C had i only) and using the standard C det i correction factors for detector effects.
The effect of the hadronization process on the four-jet rates and differential two-jet rates was studied by using the PS version of HERWIG instead of the PS version of JET-SET and by varying the Lund hadronization parameters σ q and a as for the angular correlations. The parameter Q 0 defining the end of the parton shower cascade was also varied from its standard value of 1.9 GeV to 2.4 GeV. The theoretical predictions assume zero quark masses. To assess the effect of this assumption, we determined the hadronization correction, C had , using JETSET Monte Carlo events with only light primary quarks (u, d, s, c) at the parton level but with all available flavours (u, d, s, c, b) at the hadron level. The results provided by these variations are presented in Table 4.

Renormalization scale ambiguity
The arbitrariness in the choice of the renormalization scale µ leads to an ambiguity in the theoretical prediction. We perform our default analysis at the scale factor x µ = µ/E cm = 1.
To estimate the systematic error due to the scale ambiguity, we varied the renormalization scale within the range 0.5 ≤ x µ ≤ 2. Varying x µ probes the effect of missing higher orders, and these could be very different for the NLO and NLLA calculations. We performed two variations separately: varying x µ for the jet-related variables while keeping it fixed at the default value for the four-jet angular correlations, or vice versa. Table 5 presents the results of this systematic check. As expected, the normalised angular correlations were not sensitive to the scale choice, but the jet related variables showed significant sensitivity. In Fig. 3 we show the dependence of the χ 2 of the fit on x µ , for the case of jet-related variables with x µ = 1 kept fixed for the angular correlations. The minimum of the χ 2 is seen to lie close to x µ = 1.

Matching ambiguity
For the differential two-jet rates, we performed the fit using the R-matching scheme to combine the NLL and NLO approximations rather than the ln R-matching. The result of this variation is listed in Table 5.

Fit range
The uncertainty related to the number of bins selected for the fit was evaluated by repeating the fit using those bins for which the values of the correction factors were in the range of 0.8 < C i tot < 1.2 instead of the default range 0.9 < C i tot < 1.1 or 0.85 < C i tot < 1.15. The results are presented in Table 5. The variation of the fit range for the jet-related variables does not influence the value of the χ 2 significantly. The fit-range variation for the angular variables cos χ BZ and cos α 34 is the main source of the much increased value of the χ 2 .

Five-parton background
Five-parton events which are reconstructed as four-jet events at the parton level are taken into account by the NLO theory. We examined the effect of those five-parton events which yield four jets at the detector level but five at the parton level. To reduce this five-parton background, we applied a cut on y 45 , the value of the resolution variable which separates the four-and five-jet configurations of any given event. We optimized the value of this cut by studying five-parton events generated with the PENTAJET Monte Carlo program [41]. This program generates a five-parton configuration using the same value of the intrinsic y 0 as the JETSET-ME program, then hadronizes the partons using JETSET string fragmentation with the ERT-MC-3 parameters of Ref. [34]. The events were passed through the OPAL detector simulation program. Of 20 000 five-parton events, 2 517 were found to be four-jet events at detector level, but five-jet events at parton level, for y cut = 0.008. The number of these events reduces to 376 when we keep events with y 45 < 0.003, eliminating 85 % of the five-parton background, while leaving 60 % of the data.
The results of the fit with this additional requirement are presented in Table 5.

Dependence on the parameter y cut
The value of the jet resolution parameter y cut , used in selecting four-jet events to define the angular correlations, can be chosen arbitrarily, as long as it is safely larger than the intrinsic y 0 of the JETSET-ME program. In order to have a high statistics sample, our default choice was the smallest possible value, which allows us to examine the stability of the results when y cut is increased. The four-jet event sample decreases very rapidly with increasing y cut , which constrains the region in y cut where a reliable systematic check can be performed. We repeated the entire analysis using y cut = 0.01. The results of this fit are presented in Table 5.

Additional systematic checks
We tested the influence of the two-and three-parton events that give rise to four-jet events at the detector level by repeating the analysis requiring E min = 5 GeV on the minimum energy of the jets instead of the default value E min = 3 GeV. A higher value of the jet energy reduces the number of two-and three-parton events which satisfy the four-jet selection criteria in the detector. We generated 1 million events with two, three and four partons in the initial state using the ME version of JETSET with renormalization scale set to x 2 µ = 0.002, which provides the correct proportion of two, three and four jet events for y cut = 0.008. We then counted the proportion of the two-and three-parton events which passed the four-jet selection at detector level. With E > E min = 5 GeV, there are no background events from the two-and three-parton events. With E > E min = 3 GeV, a background of 1 % appears (519 two-or three-parton events in the sample of 53 302 four-jet events selected at the detector level). We found negligible effects on the final results of the fit, therefore, we did not include the effect of this check in the systematic errors.
As stated in Section 2, the fits to individual angular correlations give rather ambiguous results. Therefore, we choose to use all four angular correlation variables in our standard analysis. As a consistency check, we performed the fit using only five variables instead of six, removing one of the six variables at a time. The results of these fits are presented in Table 6. The errors are statistical only. We found that leaving out any one of the observables yields consistent results for the measured parameters within the statistical errors of the measurements. Therefore, we did not include the effect of this check in the systematic errors.
The number of light quark flavours in a fixed order calculation is usually set to N f = 5 at LEP energies because the hard scattering scale is assumed to be the e + e − center of mass energy. However, one may argue that in the four-quark channel the g → bb splitting is suppressed compared to the g → qq (q=u, d, s, c) splittings, because the scale of that process is much lower than the scale at the primary vertex. To study the effect of suppressed g → bb splitting on our results, we also used N f = 4 light flavours at the secondary vertex in the fixed order prediction for the angular correlations. The results of this check are also shown in Table 6. The colour factor ratios increase slightly, remaining well within the statistical errors. The quality of the fit is almost unchanged. We did not include the effect of this check in the systematic errors, because it is not customary to do so.

Total systematic error
To combine the effects of the systematic variations we employed a Bayesian method that was developed in a similar analysis by the ALEPH collaboration [9], where the explicit formulae we used can be found. The basic idea is that the information obtained from a fit is accepted or rejected based on the quality of the fit, i.e. the magnitude of the χ 2 . From a Bayesian point of view we presume that all models are equally well suited to the analysis. However, a large χ 2 /d.o.f. indicates that the probability of a model is low and, therefore, it should contribute to the systematic errors with a small weight. Accordingly, the systematic error corresponds to an increase of the χ 2 by 1.
We could not use the Bayesian method in two cases. In the variation of the fit range the number of bins changed, while in checking the renormalization scale uncertainty for the angular correlations we did not find a pronounced minimum of the χ 2 curve. Furthermore we did not use the Bayesian method to estimate the systematic uncertainty from the fiveparton background because the event statistics for this check were modified significantly, by about 40%, relative to the standard analysis. In these cases we estimated the systematic error by the deviation of the results from the standard values.
The errors obtained from the different systematic checks are summarized in Table 7. The theoretical uncertainty refers to the larger of the errors from the scale uncertainty of the jet-related variables, the matching ambiguity (both of those estimate the effect of unknown higher order perturbative contributions), and the scale uncertainty of the fourjet angular variables. The different sources were added in quadrature to define the total systematic error. 6 Taking into account the systematic correlations presented in Table 7, we converted these values to the errors of the standard QCD parameters α s (M Z ) and colour factors using the definitions in Eq. (10) and (12) and T R = 1/2. The result of the conversion is presented in Table 7. The errors are dominated by the theoretical uncertainty.  Table 7: Contributions to the total systematic error. 6 If the Bayesian approach were not adopted and the systematic errors were estimated by the largest deviation from the standard values, then we obtain the following systematic errors: ∆η = 0.0027, ∆x = 0.27, ∆y = 0.08. For the basic QCD parameters these errors are: ∆α s = 0.030, ∆C A = 0.72, ∆C F = 0.31.
To compare the results to previous measurements, we show in Fig. 4 the two dimensional 68% and 95% confidence level contour plots of the colour factor ratios based on total uncertainties obtained by our analysis.
The DELPHI Collaboration [12] has performed a least-squares fit of leading order predictions to the two-dimensional distribution in the variables cos Θ NR and cos α 34 in order to measure of the colour factor ratios. The DELPHI Collaboration [42] also performed this analysis with b-tagging to separate quark jets from gluon jets which increases sensitivity. The OPAL Collaboration fitted the LO prediction to a three-dimensional distribution of the cos Θ NR , cos α 34 and cos χ BZ variables [13]. In a second measurement the OPAL Collaboration used event shape distributions of LEP1 data for fitting α s and one of the colour factors [43]. Furthermore, another recent measurement of the QCD colour factors and α s using event shape distributions at several energy points and power corrections was carried out [44]. The ALEPH Collaboration [9] was the first to present a simultaneous measurement of the strong coupling and the colour factors. The distribution D 2 and all four angular distributions defined in Section 2 have been employed. The results of the ALEPH experiment were determined using leading order predictions for the four-jet observables. The systematic uncertainties include contributions from variation of the renormalization scale, the matching ambiguity, from hadronization and detector effects, and finally from the estimation of mass effects. Our results are in agreement with these previous results. The main new feature of our analysis is that we obtained our results using NLO predictions for four-jet angular correlations. Furthermore, our systematic checks contain additional items arising from the estimate of the effects of the five-parton background and the variation of the y cut parameter. We also checked that assuming four rather than five light flavours at the secondary vertex does not influence our results. This effect was not examined in any of the previous studies.
With the inclusion of the higher order corrections, the colour factor ratios became better constrained, especially T R /C F . We found an increase of about 15% on T R /C F using the NLO theoretical predictions instead of the LO one, while C A /C F was only slightly affected (about 3%).
The diamond symbol in Fig. 4 indicates the result one would expect if a leading order perturbative theory with light gluinos was the correct theory. This point falls well outside our 95% C.L. contour, indicating that the data do not favour the existence of light gluinos. A complete analysis to test for the existence of gluinos should use hadronization models and NLO predictions that include the effect of the gluinos [7]. The latter is possible using the DEBRECEN partonic Monte Carlo program we used for generating the NLO results (see the Appendix), but presently none of the existing hadronization models include gluinos.
Finally, we converted our measured parameters to the standard QCD parameters (using T R = 1/2), which leads us to our main results: α s = 0.120 ± 0.011(stat.) ± 0.020(syst.) C A = 3.02 ± 0.25(stat.) ± 0.49(syst.) C F = 1.34 ± 0.13(stat.) ± 0.22(syst.) Other techniques [22] provide precise determinations of α s . In this study the measurement of α s serves as a crosscheck and a consistent value was found. The larger uncertainty we obtain for α s compared to other techniques can be traced to the functional dependence on the other two parameters.

Summary
In this paper, a test of perturbative QCD at LEP has been presented. With only the assumption of non-abelian gauge symmetry and standard hadronization models, we measured the eigenvalues of quadratic Casimir operators -the colour factors -of the underlying gauge group, together with the coupling of QCD. The measurement was based on the comparison of next-to-leading order perturbative predictions for four-jet angular variables, and the resummation improved next-to-leading order perturbative predictions for the multi-jet rates, to OPAL data taken between 1991 and 1995 at the Z peak. The results for the colour factors are as follows:

A Appendix: Theoretical Predictions
The NLO perturbative expansion of the differential two-jet rates has the following expression: where y 23 is the y cut value which separates the two-and three-jet configurations of an event. The corresponding formula in the case of the four-jet rates takes a similar form: In Eqs. (8) and (9) σ tot is the total hadronic cross section at O(η) accuracy and In the second equality in Eq. (10), we used the two-loop expression for the running coupling with and µ denotes the renormalization scale (x µ = µ/ √ s, where √ s is the total c.m. energy).
We define the colour factor ratios as where N f is the number of active flavours. Then the coefficients in the expansion of the β function can be written as The functions A ′ D 2 , B R 4 and B ′ D 2 , C R 4 in Eqs. (8) and (9) are the perturbatively calculable coefficient functions in the Born approximation and the radiative correction, respectively, which are independent of the renormalization scale. The A ′ D 2 function is independent of x and y. The B functions are linear forms and the C 4 function is a quadratic form of these colour charge ratios [45]: where with being a cubic Casimir invariant of the gauge group. 7 We calculated these coefficient functions of the differential two-jet rates in the −6 ≤ ln(y 23 ) ≤ 0 y 23 range and those of the four-jet rates in the −3 ≤ log 10 (y cut ) ≤ 0 y cut range using the partonic Monte Carlo program DEBRECEN [46]. All theoretical predictions were obtained for N f = 5 light-quark flavours and with the normalization T R = 1 2 . The latter choice is arbitrary and influences the values of the coefficient functions in such a way that the physical prediction remains unchanged. 7 The results of the Monte Carlo integration for C (4) z are compatible with zero [18], therefore, z is not included as a free parameter in our fits.
The resummation formulas for various event shapes are given in Ref. [47] at the NLL accuracy for the cumulative cross section, in the following form: where L = − ln(y 23 ). The functions Lg 1 (ηL) and g 2 (ηL) represent the sums of the leading and next-to-leading logarithms, respectively, to all orders in the strong coupling. The colour decomposition of the g 1 and g 2 functions is given in Ref. [48]. The coefficient C 1 = −5/2 + π 2 /6 − 6 ln 2 is known exactly [49] and C 2 can be determined from numerical integration of the O(α 2 s ) matrix elements (see below). At NLO (i.e. O(α 2 s ) accuracy), the cumulative cross section for the differential two-jet rates has the form: In order to have the best possible theoretical description, we fit the experimental data to the matched NLL and NLO results, which is expected to describe the data in a wider kinematical range than any of the two separately. In the case of D 2 , our default matching procedure is the ln R-matching, which combines the NLL and NLO approximations according to the following formula: where Σ(L, η) = Lg 1 (ηL) + g 2 (ηL) and the following two terms represent the two lowest order terms in the η-expansion of ln R. In the actual analysis, we used the differential distribution dσ/dL that is obtained from Eq. (21) as with A ′ (L) and B ′ (L) given in Eq. (8).
Eq. (22) contains the cumulative coefficients A(L) and B i (L), the latter ones implicitly in R, that cannot be directly calculated using a partonic Monte Carlo program. We know, however, that and where the colour decomposition of the R (i) 2m logarithmic coefficients can be found in Ref. [5], except for R (i) 21 . This latter coefficient as well as the C (i) 2 constants and the D 1 (L), D (i) 2 (L) functions can be determined by differentiating Eqs. (23) and (24) with respect to L and fitting to the coefficient functions as obtained from the partonic Monte Carlo program DEBRECEN. Following the fitting procedure of Refs. [47,48] we obtained R , C x 2 = −1.53 ± 1.12(stat.), C y 2 = −3.565 ± 1.17(stat.). In order to account for the dependence on the renormalization scale x µ , one has to make the following substitutions in Eq. (22): The other widely used matching scheme is R-matching, which will serve as a systematic check of our matching procedure. In this case the NLL and NLO predictions are combined according to the following formula [1,16]: where B NLL and C NLL are the coefficients in the η-expansion of R NLL . In the case of multijet rates, the NLL approximation does not exponentiate, therefore, the viable matching scheme for the four-jet rates is R-matching. The resummation formulas for the four-jet rates can be found in Refs. [16,18]. The general form of the NLO differential cross section for the four-jet angular correlations O 4 is given by [14]: where σ 0 denotes the Born cross section for the process e + e − →qq. In Eq. (29) the Born and correction functions B O 4 and C O 4 have an analogous decomposition to Eqs. (14) and (15). To obtain the distributions normalised to unity in perturbation theory, we calculated the coefficient functions B    Figure 2: Distributions of the four-jet angular correlations. OPAL data corrected to the parton level are denoted by points. The statistical uncertainty in the bins is smaller than the size of the circles. The small sections below each plot indicate the correction factors C tot i for experimental and hadronization effects with their statistical error (shaded area). The horizontal lines with arrows indicate the selected fit region. For the normalised angular correlations, the distributions at leading and next-to-leading order are almost the same [14,19]