Inclusive Analysis of the b Quark Fragmentation Function in Z Decays at LEP

A study of b quark hadronisation is presented using inclusively reconstructed B hadrons in about four million hadronic Z decays recorded in 1992-2000 with the OPAL detector at LEP. The data are compared to different theoretical models, and fragmentation function parameters of these models are fitted. The average scaled energy of weakly decaying B hadrons is determined to be=0.7193+-0.0016(stat)+0.0036-0.0031(syst)


Introduction
Hadronisation, the transition of quarks into hadrons, is a strong interaction phenomenon which cannot yet be calculated from first principles within QCD.Monte Carlo event generators are used instead which rely on phenomenological models of this process.To some extent these models can be distinguished from each other by the shape of the predicted hadron energy distribution.Hadronisation of heavy quarks leads to a significantly harder hadron energy spectrum than for lighter quarks [1].Experimentally, heavy quark hadronisation is of special interest, because in this case the hadron containing the primary quark can easily be identified.
A precise measurement of the B hadron1 energy distribution allows the various hadronisation models available to be tested, and also helps to reduce one of the most important systematic uncertainties in many heavy flavour analyses.Earlier measurements of the B hadron energy distribution usually fell into one of three categories.1) Some analyses were based on a measurement of the energy distribution of certain exclusive B hadron decays, mostly B → D * ℓν, to constrain the B hadron energy as precisely as possible [2,3].However, this leads to small candidate samples and thus to a large statistical uncertainty.2) Other analyses attempted to increase the sample size by extracting the energy distribution of leptons from inclusive B → ℓ decays.Unfortunately the modelling of the lepton energy spectrum introduces large additional systematic uncertainties [4,5].3) The most precise results so far have been achieved through the fully inclusive reconstruction of the B hadron energy [6].The analysis presented here identifies B hadrons inclusively using secondary vertices.

Data sample and event selection
This analysis uses data taken at or near the Z resonance with the OPAL detector at LEP between 1992-2000.A detailed description of the OPAL detector can be found elsewhere [7].The most important components of the detector for this analysis are the silicon microvertex detector, the tracking chambers, and the electromagnetic calorimeter.The microvertex detector consisted of two layers of silicon strip detectors which provided high spatial resolution near the interaction region.The central jet chamber was optimised for good spatial resolution in the plane perpendicular to the beam axis 2 .The resolution along the beam direction was improved by the z information delivered by the silicon microvertex detector (except in the first version present in 1992), by a vertex drift chamber between the silicon detector and the jet chamber, and by dedicated z-chambers surrounding the other tracking chambers.The central detector provided good double track resolution and precise determination of the momenta of charged particles by measuring the curvature of their trajectories in a magnetic field of 0.435 T. The solenoid was mounted outside the tracking chambers but inside the electromagnetic calorimeter, which consisted of approximately 12 000 lead glass blocks.The electromagnetic calorimeter was surrounded by a hadronic calorimeter and muon detectors.
Hadronic events are selected as described in Ref. [8], giving a hadronic Z selection efficiency of (98.1 ± 0.5) % and a background of less than 0.1 %.Only data that were taken with the silicon microvertex detector in operation are used for this analysis.A data sample of about 3.91 million hadronic events is selected.This includes 0.41 million events taken for detector calibration purposes during the years 1996-2000, when LEP was operating at higher energies.
A total of 23.81 million Monte Carlo simulated events are used: 16.81 million events were generated with the JETSET 7.4 generator [9], 2 million events were generated with HERWIG 5.9 [10], and 5 million events were produced by HERWIG 6.2 [11].The JETSET event sample includes 4.93 million bb events and 3.19 million cc events in dedicated heavy flavour samples.All other samples are mixed five flavour Z → qq event samples.The choice of important parameters of the event generators is described in [12].All Monte Carlo simulated events are passed through a detailed detector simulation [13].The same reconstruction algorithms as for data are applied to simulated events.
The analysis is performed separately for the data of different years, where detector upgrades, in particular of the silicon microvertex detector [14], and recalibrations lead to different conditions.Separate samples of JETSET Monte Carlo are available for all years.HERWIG Monte Carlo is only available for the largest homogeneous dataset taken in 1994, and therefore HERWIG-related studies are performed exclusively for this dataset.
In the 1993 and 1995 runs, part of the data was taken at centre-of-mass energies about 1.8 GeV above and below the peak of the Z resonance.The B hadron energy distribution is sensitive to energy losses due to initial state radiation prior to the annihilation process.Initial state radiation is heavily suppressed at and just below the Z resonance, but it has significant impact in the dataset taken at an energy of m Z + 1.8 GeV.The latter samples are therefore treated separately, with Monte Carlo samples simulated for the appropriate energy, giving a total of eleven separate data and JETSET Monte Carlo samples.

Preselection of Z → bb events
The thrust axis is calculated for each event using tracks and electromagnetic clusters not associated with any tracks.To select events within the fiducial acceptance of the silicon microvertex detector and the barrel electromagnetic calorimeter, the thrust axis direction is required to satisfy | cos θ T | < 0.8, where θ T is the thrust angle with respect to the beam direction.
To achieve optimal b-tagging performance, each event is forced into a 2-jet topology using the Durham jet-finding scheme [15].In calculating the visible energies and momenta of the event and of individual jets, corrections are applied to prevent double counting of energy in the case of tracks with associated clusters [16].A b-tagging algorithm is applied to each jet using three independent methods: lifetime tag, high p t lepton tag and jet-shape tag.This algorithm was developed for and used in the OPAL Higgs boson searches.A detailed description of the algorithm can be found in [17].Its applicability to events recorded at the Z resonance peak was already shown in [18].The b-tagging discriminants calculated for each of the jets in the event are combined to yield an event b likelihood B event .Both the jet b-tagging discriminant and B event have values between zero and one and correspond approximately to the probability of a true b jet or bb event, respectively.For each event, B event > 0.2 is required.The Z → bb event purity is 83% after this requirement, and the efficiency is 54% at this stage.
The b hemisphere tag efficiency obtained from Monte Carlo simulation is compared to the actual value in data using a double tag approach as described in [19].The efficiencies obtained this way in both simulation and real data are found to agree to within 5% in all subsamples.Nevertheless a correction is applied to the Monte Carlo efficiency to further improve the agreement.

Reconstruction of B hadron energy
The primary event vertex is reconstructed using the tracks in the event constrained to the average position of the e + e − collision point.For the B hadron reconstruction, tracks and electromagnetic calorimeter clusters with no associated track are combined into jets using a cone algorithm3 [20] with a cone half-angle of 0.65 rad and a minimum jet energy of 5.0 GeV.The two most energetic jets of each event are assumed to contain the B hadrons.Only jets where the opposite hemisphere yields a b-tagging discriminant of at least 0.8, corresponding to a b probability of about 80%, are used in the analysis.The distribution of the b-tagging discriminant is shown in Figure 1.
Each remaining jet is searched for secondary vertices using a vertex reconstruction algorithm similar to that described in [21], making use of the tracking information in both the rφ and rz planes where available.If a secondary vertex is found, the primary vertex is re-fitted excluding the tracks assigned to the secondary vertex.Secondary vertex candidates are accepted and called 'good' secondary vertices if they contain at least three tracks.If there is more than one good secondary vertex attached to a jet, the vertex with the largest number of significant4 tracks is taken.If there are two or more such vertices, the secondary vertex with the larger separation significance with respect to the primary vertex is taken.Jets without an associated secondary vertex are rejected.This increases the b jet purity and improves the energy resolution of the B hadron reconstruction described in the following.
Weakly decaying B hadrons are reconstructed inclusively with a method described in an earlier publication [22].In each hemisphere defined by the positive axis of the jet found by the cone algorithm, a weight is assigned to each track and each cluster, where the weight corresponds track neural network track momentum track rapidity with respect to estimated B hadron flight direction track impact parameter with respect to primary vertex in rφ projection (d 0 ) track impact parameter with respect to primary vertex in z projection (z 0 ) d 0 impact parameter significance z 0 impact parameter significance 3d impact parameter significance with respect to the primary vertex 3d impact parameter significance with respect to the secondary vertex cluster neural network cluster energy cluster rapidity with respect to estimated B hadron flight direction Table 1: Variables used in artificial neural networks to estimate the probability that a track or calorimeter cluster originates from a B hadron decay.The impact parameter significance is defined as the impact parameter divided by its uncertainty.
to the probability that this track or cluster is a product of the B hadron decay.This weight is obtained from artificial neural networks [23] exploiting information from track impact parameters with respect to the primary and secondary vertices, and from kinematic quantities like the transverse momentum associated with a track or cluster, measured with respect to the cone jet axis.A list of all variables is shown in Table 1.The B hadron momentum is then reconstructed by summing the weighted momenta of the tracks and clusters.A beam energy constraint assuming a two-body decay of the Z and the world average B meson mass of 5.279 GeV/c 2 [24] for the B hadron is applied to improve the energy resolution.The constraints lead to a biased energy reconstruction, particularly when the true B hadron energy is very small, as can be seen in Figure 2a.However, only a small fraction of the data sample is in the low-energy region affected by a large bias.For most events, in the peak of the B hadron energy distribution, the bias is small, and all biases are taken into account by the fitting procedures used in both the model-dependent and model-independent analyses.Possible systematic uncertainties arising from the biased energy reconstruction are discussed in Section 7 of this paper.The energy of the weakly decaying B hadron is expressed in terms of the scaled energy x E = E B /E beam , where E beam = √ s/2 is the LEP beam energy for the event.The quantity x E is restricted to values above 5.279 GeV/E beam ≈ 0.1 by the B meson mass constraint, and it cannot exceed 1.0 due to the beam energy constraint.After all these requirements, the distribution of the difference between the reconstructed energy and that of generated B hadrons in simulated data has a rms width of 4.8 GeV.The energy dependence of the B hadron energy resolution is shown in Figure 2a.The complete B hadron selection applied to the full data sample results in 270 707 tagged jets with a b purity of 96%.The average B hadron selection efficiency is 16%, with an energy dependence as shown in Figure 2b.The measured B hadron energy distribution, scaled to the beam energy, is shown in Figures 3-6, and compared to the various models described in the next section.
Table 2: Fragmentation functions for the JETSET 7.4 string scheme that are fitted to data in this paper.
N is a normalisation constant, different for each fragmentation function.

Test of hadronisation models
The B hadron energy distributions predicted by the JETSET 7.4, HERWIG 5.9, and HERWIG 6.2 Monte Carlo models are compared to the OPAL data.All Monte Carlo simulated events are passed through a detailed detector simulation [13].The comparison is performed using the distribution of the reconstructed scaled energy of the weakly decaying B hadron x E .
The HERWIG Monte Carlo uses a parton shower fragmentation followed by cluster hadronisation model with few parameters.No parameters are varied in this analysis.This simplifies the model test to a mere comparison of the x E distributions obtained with data and Monte Carlo simulation.Both HERWIG versions are set up to conserve the initial b quark direction in the B hadron creation during cluster decay (cldir=1).The main difference between the two HERWIG samples used in this analysis is that Gaussian smearing of the B hadron direction around the initial b quark flight direction is applied in the HERWIG 5.9 sample (clsmr=0.35),while smearing is not used in the HERWIG 6.2 sample (clsmr(2)=0).
The JETSET Monte Carlo is based on a parton shower fragmentation followed by string hadronisation scheme.It requires a fragmentation function to describe the distribution of the fraction z of the string light cone momentum that is assigned to a hadron produced at the end of the string.The JETSET sample in this analysis is reweighted to use the fragmentation functions of Kartvelishvili et al. [25], Bowler [26], the Lund symmetric model [27], and the fragmentation functions of Peterson et al. [28], and Collins-Spiller [29].The Lund symmetric and Bowler functions are simplified by assuming the transverse mass of the b quark, m ⊥ , to be constant, which is justified by the smallness of the average transverse momentum compared to the b quark mass.A further simplification in the Bowler parametrisation is the assumption of an equality of b quark and hadron masses.The functional forms of the fragmentation functions are given in Table 2.The parameters of the respective fragmentation functions are fitted to obtain a best match of the observed x E distributions in data and Monte Carlo simulation.In the case of the Peterson et al., Collins-Spiller, and Kartvelishvili et al. models, one free parameter is available.The Lund and Bowler models each have two free fit parameters.A χ 2 fit is performed in 46 bins in the x E range of 0.5 to 0.95, where in all samples the number of candidates in each bin is large enough to justify the assumption of Gaussian errors on the bin content.The fragmentation function and its parameters are adjusted during the fit by reweighting the Monte Carlo simulated events, similar to the procedure applied in [19].
The reweighting fit is performed separately for each data sample, and the fit results are averaged with weights according to the size of the respective datasets.Consistent results are obtained for all datasets.The average parameter values are summarised in Table 3.For each parametrisation the corresponding model-dependent mean scaled energy of weakly decaying B hadrons is given.Data samples at √ s = m Z + 1.8 GeV are excluded in the calculation of the average x E .6 Model-independent measurement of x E In the previous section, information was extracted from the observed energy distribution making explicit use of a set of models to describe the data.In this section, a measurement of the mean scaled energy of B hadrons, x E , outside a specific model framework will be presented.This is accomplished by unfolding the observed energy distribution.Two complementary unfolding procedures are used to obtain an estimate of the true x E distribution from the observed distribution of the reconstructed scaled B hadron energy.In both cases the amount and energy distribution of background in the B hadron candidate sample  is estimated from the Monte Carlo simulation and subtracted from the data.
The main method starts by fitting the observed data x E distribution, and the observed and the true x E distribution in the Monte Carlo simulation with smooth functions (splines).The true and observed Monte Carlo distributions are then reweighted simultaneously until the observed x E distribution agrees in data and simulation.The reweighted true x E distribution of the Monte Carlo simulation then provides an estimate of the corresponding distribution in data.Details of how the result is stabilised are described later.This method is almost independent of the initial Monte Carlo distribution and thus reduces model-dependence in the unfolding process.Furthermore, the result is represented as an unbinned spline function, which is optimal for the calculation of the mean value of the unfolded distribution.This algorithm is coded using the software package RUN [30] and was already used in [31].
The second approach makes use of the SVD-GURU software package [32].The correspondence between the observed and true B hadron energy distributions in the Monte Carlo simulation is represented by a 20 × 20 matrix.The unfolding process comprises a matrix inversion to obtain an estimate of the true data x E distribution from the observed distribution.In this approach, the model dependence was found to be stronger than when using the RUN program.Furthermore, a coarse binning appropriately adapted to the detector's resolution and the amount of available data might lead to systematic effects when describing the energy distribution in terms of its mean value.Therefore SVD-GURU is only used to cross-check the result obtained by RUN and to provide an estimate of the systematic uncertainty due to unfolding.
Raw unfolding solutions often oscillate strongly around the correct solution.In the case of a binned representation of the data this effect can simply be understood by strong negative bin-bybin correlations introduced by the finite detector resolution.Both methods used here suppress these oscillations by limiting the number of degrees of freedom of the unfolding solution.The RUN algorithm represents the unfolding result as expansion into a set of orthogonal functions.The uncertainties on the coefficients of these functions are determined, and only those functions with coefficients significantly different from zero are taken into account.SVD-GURU rotates the unfolding matrix to estimate its effective rank.The unfolding is then performed in a rotated space with a smaller matrix including only the significant contributions.The number of degrees of freedom used for the unfolding procedure was found to agree in the RUN and SVD-GURU approaches in all subsamples described below.A further means of regularisation is available in the RUN package.Of all remaining solutions to the unfolding problem, one is chosen that minimises the integral over the squared first derivative of the unfolding solution.Monte Carlo studies show that this regularisation leads to essentially bias-free results on all samples.The performance of the unfolding algorithms is illustrated in Figure 7.
The unfolding is performed separately for data from all years of 1992-2000.The 1993 and 1995 datasets at a centre-of-mass energy below m Z show x E distributions that are compatible with those taken at the Z resonance peak in the Monte Carlo simulation.The 1993 and 1995 datasets at m Z + 1.8 GeV show a significantly lower x E , caused by a large amount of initial state radiation at this energy.In this case the quark energy prior to fragmentation is lower on average than the beam energy.As the beam energy is used as an estimator of the quark energy prior to fragmentation, the average x E value is lower than in samples without significant initial state radiation.The m Z + 1.8 GeV samples are therefore analysed separately.
Both RUN and SVD-GURU analyses are performed with a Monte Carlo simulation that is reweighted to match the best result of the model-dependent reweighting fits for the respective datasets.This procedure is also followed by SLD in their latest b hadronisation analysis [6].The goal is to reduce the dependence of the unfolding result on the Monte Carlo sample used for unfolding.The effect of not using the best parametrisation, but the second and third best instead, are studied below as a systematic effect.JETSET 7.4 Monte Carlo simulation samples are used to obtain the central result, and will be compared to unfolding results using HERWIG in the discussion of systematic effects below.
The unfolding result for all data samples at or below the Z resonance is shown in Figure 8.The mean scaled B hadron energy, obtained with the RUN unfolding algorithm and averaged over all datasets by using the subsample size for the weights, is where the uncertainty includes the statistical uncertainties due to limited data and Monte Carlo sample sizes, and the statistical uncertainty on the Monte Carlo efficiency.Consistent x E values were obtained for the individual data samples (see Figure 8).This model-independent measurement agrees with the x E values in the framework of the best models as seen in Table 3.The unfolding result spline for the full dataset is plotted in Figure 9.
The mean scaled energy observed in the m Z + 1.8 GeV samples is found to be x E = 0.7130 ± 0.0056(stat) .
The 1993 and 1995 m Z + 1.8 GeV data samples give consistent values of x E .The difference of the results for the different energies is consistent with the prediction obtained from Monte Carlo samples at similar energies.The results obtained with the unfolding program SVD-GURU (0.7195 ± 0.0015(stat) for the main dataset, 0.7152 ± 0.0053(stat) for the m Z + 1.8 GeV samples) are in very good agreement with the ones achieved with the RUN package.This is also demonstrated in Figure 9, where the results obtained from the full √ s ≤ m Z data sample by both algorithms are compared.The statistical uncertainties of both methods are also very similar.

Systematic uncertainties
Given the large data sample collected with the OPAL detector, and the inclusive character of the analysis presented here, the statistical uncertainties on the results of the previous sections are expected to be small compared to the systematic uncertainties introduced by limited knowledge of physics parameters which possibly affect the measured quantities.In this section, an overview of all systematic checks is given for the reweighting fit and for the unfolding analysis.
The distribution of the b-tagging discriminant in data agrees well with the Monte Carlo simulation.However, it is necessary for this analysis to ensure that this is separately true for different B hadron energy regions.Therefore the b-tagging discriminant was investigated in ten bins in x E .The ratio of the b-tagging discriminant distributions of data and Monte Carlo simulation is calculated and fitted by a linear function in each energy bin separately.The slope of this function is an indicator of the quality of agreement of data and simulation.All fitted slopes are compatible with zero within two standard deviations.

Systematic uncertainties for fits to various models
The systematic checks performed for the fits to models are listed below.The resulting systematic uncertainty estimates for the mean scaled energy x E in the framework of the respective models are summarised in Table 4.
• The energy resolution of the OPAL calorimeters in the Monte Carlo simulation is varied by ±10% relative to its central value.This range is motivated by jet energy resolution studies in two-jet events, where a difference between the resolution in data and the Monte Carlo simulation of 3.6% was found in some datasets.Under the assumptions that 50% of the total jet energy are contributed by neutral particles and that the observed difference can be fully accounted to the calorimeters, the above variation range covers this effect.
• Similar studies indicate a possible difference of the energy scale between data and simulation of up to 0.4% in some datasets.The energy scale is varied within this range, and the resulting difference of the fit results is taken as systematic uncertainty.
• Figure 2a shows a large energy reconstruction bias for low energy B hadrons.Both the number of low energy B hadrons and the efficiency for reconstructing them are small.This large bias therefore only affects a very small fraction of the candidates.The effect of possible mismodeling of the bias in the Monte Carlo simulation is evaluated by varying the bias around the values that describe the data best.The reconstruction bias in the high energy region cannot be changed by more than ±1% without leading to a significant degradation of the agreement between data and Monte Carlo simulation.The bias for low energy B hadrons, with a reconstructed x E below 0.6, is varied by ±10%.This range is also motivated by significant degradation of the agreement between data and Monte Carlo simulation.The largest deviations of the measured quantities are taken as systematic uncertainties.
• The relative fraction of different B hadron species in the sample of primary B hadrons influences the measurement, because different B hadron species have different energy distributions.All values obtained in this analysis are calculated with Monte Carlo samples that are reweighted to reflect the current best knowledge of the hadron fractions.The associated systematic uncertainty is estimated by varying the b baryon fraction within the range (10.3 ± 1.8)% [33] given by the average of the LEP/SLD/CDF measurements of this quantity.The fraction of B s in the sample is varied in the range of (9.8 ± 1.2)% [33].
• The amount of orbitally excited B ( * ) J mesons has been measured by all LEP collaborations [21,34].An error weighted average of the LEP measurements is (28.4 ± 3.5)%, and the fraction of orbitally excited B ( * ) J mesons is varied within this range.
• The Q-value of orbitally excited B ( * ) J mesons [24] is about 40 MeV smaller than the value in the Monte Carlo samples used in this analysis.All results are corrected for this effect, and the difference to the values obtained without correction is taken as systematic uncertainty.
• The average multiplicity of charged particles from a B hadron decay was found at LEP to be 4.955 ± 0.062 [33], and is varied within this range.
• The average lifetime of weakly decaying B hadrons affects the efficiency of secondary vertex reconstruction, and is varied in the range (1.577 ± 0.016) ps [24].
• The average lifetime of weakly decaying charm hadrons determines the amount of charm background found in the B hadron candidate sample.The D 0 , D + , D s + , and Λ c + lifetimes are varied within the uncertainties quoted in [24].
• The charm background in the Monte Carlo simulation samples is reweighted to the same hadronisation model as the B hadron distribution in the respective fits.The central values and uncertainties of the charm fragmentation function are taken from earlier OPAL measurements [35].For the evaluation of the systematic uncertainty the parameters are varied within their uncertainties.
• Jets from gluon splitting to bb quark pairs are treated as background in this analysis.
To account for the uncertainty of the average LEP result of the gluon splitting rate, 0.00254 ± 0.00050 bb pairs per hadronic event [33], the rate is varied within this range.
• The number of cc pairs from gluon splitting per hadronic event is varied within the LEP uncertainty of 0.0299 ± 0.0039 [33].
• The partial width of the Z into bb quark pairs, normalised to the total hadronic width of the Z, is measured to be R b = 0.21646 ± 0.00065 [24].Varying this fraction within the quoted uncertainty leads to varying background levels in the unfolding Monte Carlo sample.This causes negligible changes of the fit results.
• The analogous quantity for charm quark pairs, R c , is less well known, with a current best value of 0.1719 ± 0.0031 [24].The impact on the fit results of varying R c within this range is negligible.
• Limited knowledge of the LEP beam energy produces a corresponding uncertainty on x E , although the dependency of x E on the beam energy is reduced due to the fact that the beam energy also enters the calculation of the reconstructed B hadron energy via the beam energy constraint.The assumed LEP beam energy is varied within ±8 MeV, which is the largest reported uncertainty for any sample at or close to the Z resonance [36].A correlation of 100% between the resulting uncertainties for the different data taking years is assumed.
• The parameter values depend slightly on the x E range used for the fit.The lower end of the fit range is varied within x E = 0.5 ± 0.1, and the upper range is varied within x E = 0.95 ± 0.05.The largest deviation from the result obtained using the central values is taken as the systematic uncertainty.
• The bin width used in the fit is varied by ±10%.The maximum deviation from the result with standard binning is used to estimate the associated systematic uncertainty.

Systematic uncertainties of the unfolding analysis
The same systematic uncertainties as for the fragmentation function fits are evaluated also for the x E measurement, with the exception of fit range effects, which are specific to the modeldependent fit procedure.Binning effects are not present in RUN.The resulting systematic uncertainties are summarised in Table 5.
An additional uncertainty is introduced by the dependence of the x E measurement on detector and acceptance modelling in the Monte Carlo simulation.As in all unfolding problems, one needs the resolution (or spectral) function g(x 0 , x) where x 0 is the energy of the hadrons entering the detector and x their measured energy.This function is not measured, but calculated by the OPAL detector simulation [13].Since the detector simulation is generally made in the framework of some specific Monte-Carlo program generating hadron distributions, a small residual dependency of g(x 0 , x) on the particular Monte Carlo event generator remains.To estimate the associated systematic uncertainty, the unfolding procedure was repeated using not only the best reweighting fit result, but also all other models as initial estimators of the true distribution.This study was performed independently for all datasets.The unfolding was also  performed with the JETSET Monte Carlo sample replaced by HERWIG 5.9 and 6.2 samples.
The check using JETSET with the Collins-Spiller parametrisation dominates the modelling uncertainty in the negative direction, which is taken as the largest observed deviation from the central x E value.The uncertainty in the positive direction is dominated by the third best model, which in almost all datasets is the Kartvelishvili et al. parametrisation.The resulting model uncertainty is found to be +0.0024−0.0016 .The result of the unfolding procedure with the RUN algorithm is cross-checked with the SVD-GURU method, and the difference between the two results is assigned as systematic uncertainty.Furthermore, a difference of similar size is observed between the spline unfolding result of the RUN method and a binned representation of the unfolded distribution.This difference is also included in the unfolding method uncertainty.
Summing all systematic uncertainties in quadrature yields a total systematic uncertainty on x E of +0.0038 −0.0033 .As expected, the systematic uncertainty is larger than the statistical precision.Table 6 shows a representation of the RUN result in 20 bins in the x E range between 0.1 and 1.0.This table, along with the full correlation matrix of statistical (Table 7) and systematic (Tables 8 and 9) uncertainties can be used to compare the OPAL results with further hadronisation models not discussed here.It has to be pointed out again that the x E value obtained from the binned RUN result is smaller than the unbinned result by ∆ x E = −0.0002.This is caused by binning effects which are reduced by using a small bin width, but cannot be entirely avoided.

Summary and discussion
Using an unfolding technique to reduce the dependence on b quark hadronisation models, the mean scaled energy of weakly decaying B hadrons in Z decays is measured to be x E = 0.7193 ± 0.0016(stat) +0.0038 −0.0033 (syst) .This is the most precise available measurement of this quantity.Consistent results are obtained using an alternative unfolding method and from model-dependent reweighting fits.
The result obtained here is in good agreement with a recent result from the ALEPH Collaboration [3], x E = 0.716 ± 0.006(stat) ± 0.006(syst).ALEPH uses exclusive semileptonic B decays, leading to a smaller candidate sample and thus a larger statistical uncertainty.Another new measurement by SLD [6] gives a somewhat lower value: x E = 0.709 ± 0.003(stat) ± 0.003(syst) ± 0.002(model).The difference between the OPAL and the SLD measurements has a statistical significance of less than 2 standard deviations taking only the uncorrelated uncertainties into account.Another x E measurement was recently performed using inclusive B → ℓ decays [5].Modelling of the lepton energy spectrum introduces additional systematic errors in the lepton-based analysis.The result of x E = 0.709 ± 0.003(stat) ± 0.003(syst) ± 0.013(model) is compatible with the analysis presented here, especially given that the result in [5] is not model-independent, but based on a Peterson et al. parametrisation.The LEP average result for x E in the framework of the Peterson et al. model, obtained from earlier analyses [33], is 0.702 ± 0.008, again in excellent agreement with the value of 0.7023 ± 0.0006(stat) ± 0.0019(syst) found in this analysis.
The best description of the data with a fragmentation function with one free parameter is achieved with the Kartvelishvili et al. model.The Peterson et al. and Collins-Spiller models produce energy distributions which are too broad to describe the data.Similar features have been observed by SLD and ALEPH in their recent publications.The Bowler and Lund parametrisations, each having two free parameters, achieve a better χ 2 /d.o.f. in this analysis and are clearly compatible with the data.The same conclusion is reached by SLD while ALEPH did not test these models.The HERWIG cluster model is clearly disfavoured.The main difference of the two HERWIG versions tested in this analysis is the amount of smearing of the B hadron direction around the initial b quark direction.Significant smearing is used in the HERWIG 5.9 sample, softening the spectrum too much.The HERWIG 6.2 sample is used without any smearing, giving an x E distribution which is in much better agreement with the data, but which is still too broad.Similar results are obtained by SLD.
The fitted values of the parameters describing each hadronisation model agree less well between the different experiments than the measured x E values.The parameter values depend critically on details of the Monte Carlo tuning, which is not identical in all respects among the collaborations, although efforts have been made to correct most relevant Monte Carlo parameters to a common set of values.
A general conclusion of the analysis presented here is that the parton shower plus string hadronisation Monte Carlo models provide a good description of the current data.The fragmentation functions derived from intrinsic symmetries of the string model (Bowler, Lund symmetric) are favoured over the phenomenological approaches of Kartvelishvili et

Figure 1 :
Figure 1: Distribution of the jet b-tagging discriminant in data (points with error bars) and Monte Carlo (histograms).The contributions from b jets, c jets, and light quark or gluon jets are shown as open, hatched, and black area respectively.Jets with a b-tagging discriminant above 0.8 in a jet in the opposite hemisphere are retained for analysis.

Figure 2 :
Figure 2: a) Dependence of the B hadron energy resolution (black circles) and reconstruction bias (open circles) on the generated B hadron energy.b) Dependence of the B hadron reconstruction efficiency on the generated B hadron energy.

Figure 3 :Figure 4 :Figure 5 :
Figure 3: Results of the fit to the data of various hadronisation models for JETSET 7.4.The points with error bars are the uncorrected reconstructed scaled energy distribution in the 1994 data sample.Only statistical errors are shown.The histogram represents the best match as obtained from the respective fragmentation function fits.Background from charm jets is shown as hatched histogram, and light quark and gluon background is indicated by the black area.Charm jets are preferentially passing the selection if the c quark flight length is large due to a large boost.The mean energy of reconstructed charm candidates is therefore close to that of b jets.

9 OPALFigure 6 :
Figure 6: Comparison of different setups of the HERWIG Monte Carlo generator with data.The points with error bars are the uncorrected reconstructed scaled energy distribution in the 1994 data sample.Only statistical errors are shown.The histogram represents the HERWIG prediction.Background from charm jets is shown as hatched histogram, and light quark and gluon background is indicated by the black area.

Figure 7 :
Figure 7: Performance of the unfolding algorithms used in this analysis.The dashed line represents the generated scaled energy distribution of weakly decaying B hadrons.Open circles with error bars represent the observed x E distribution for Monte Carlo simulated events corresponding to the 1994 OPAL detector setup.Shape and normalisation are different from the generated distribution due to limited and energydependent efficiency, detector resolution and reconstruction bias.Full circles with error bars and the solid line with error band indicate the SVD-GURU and RUN unfolding results for this sample.

Figure 8 :
Figure 8: Unfolding results of all data samples with a centre-of-mass energy on or below the Z mass, obtained with the RUN unfolding package.The fraction of the total data sample contributed by each subsample is given on the left.The vertical line indicates the weighted average x E , and the shaded region represents the uncertainty of this average.

Figure 9 :
Figure 9: Spline representation of the RUN unfolding result (line with error band), and binned GURU unfolding result (points with error bars), for the full data sample at √ s ≤ m Z .The narrower error band around the RUN unfolding results corresponds to the statistical uncertainty, the broader error band represents the total uncertainty.

The JETSET 7 . 4
Monte Carlo samples used for the reweighting fit were generated using the Peterson et al. fragmentation function with ε b = 38 × 10 −4 .The fit result for the Peterson et al. function in data is ε b = (41.2± 0.7) × 10 −4 .The fact that the Monte Carlo tuning and the data fit result are close has the advantage that adverse effects due to weights far from 1.0 are not expected.However, additional studies were performed to verify that the closeness of the two values is not introduced by improper reweighting.An older sample of 4 million hadronic JETSET 7.4 events with Peterson et al. fragmentation function with ε b = 57 × 10 −4 is used to repeat the fit for the 1994 data sample.The fit result obtained with this sample (ε b = (40.3± 1.1) × 10 −4 ) is in agreement with the result obtained with the main ε b = 38 × 10 −4 1994 Monte Carlo sample (ε b = (40.6 ± 1.0) × 10 −4 ).

Table 3
Lund symmetric, and Kartvelishvili et al. models are preferred by the data, with respective χ 2 /d.o.f.values of 67/44, 75/44, and 99/45 in the 1994 sample.Figures 3-6 show that the Peterson et al. and Collins-Spiller parametrisations for JETSET, as well as the HERWIG 6.2 model, are too broad.The HERWIG 5.9 model is too soft.
also gives a comparison of the fit quality of all JETSET 7.4 fits on the 1994 data, and the HERWIG 5.9 and HERWIG 6.2 results.The fit results on the 1994 data are shown in Figures3-6.The ordering of the models according to the goodness of the fits in 1994 data agrees with all other large data samples; only in a few smaller samples is a slightly different ordering observed.The quoted χ 2 /d.o.f.values only take into account the statistical uncertainty of data and Monte Carlo simulation.Systematic uncertainties are discussed later.The Bowler,

Table 3 :
Results of the comparison of hadronisation models to OPAL data.The parameter fit results and corresponding x E values are weighted averages over all datasets from the years 1992-2000, where the weights are chosen according to the subsample size.The first errors on the parameters are statistical, the second systematic.The correlation of the statistical errors of a and bm 2 ⊥ is 98.5% for the Lund symmetric model, and 96.4% for the Bowler fragmentation function.The errors on x E are the propagated statistical parameter errors.The χ 2 /d.o.f.values are quoted for the 1994 dataset only, which is the largest sample.Only statistical errors are included.The errors of the two parameters of the Lund and Bowler models are almost fully correlated.The parameters given for the HERWIG Monte Carlo are not obtained from a fit, but are the values used for the generation of each sample.

Table 4 :
Overview of systematic uncertainty contributions to the model-dependent x E measurements.

Table 5 :
Summary of all contributions to the total systematic uncertainty of the x E measurement from the unfolding analysis.

Table 6 :
al., Peterson et al., and Collins-Spiller.Unfolded x E distribution obtained from the RUN program.Statistical and systematic uncertainties are given for each bin.The corresponding correlation matrices are given in Table7(statistical uncertainties), Table8(positive systematic uncertainties), and Table9(negative systematic uncertainties).A binned representation of the RUN result will naturally lead to a slightly different x E than that calculated from the spline result.

Table 7 :
Correlation matrix of statistical uncertainties of the distribution in Table6.

Table 9 :
Correlation matrix of negative systematic uncertainties of the distribution in Table6.