Charged particle multiplicities in heavy and light quark initiated events above the Z 0 peak

We have measured the mean charged particle multiplicities separately for b¯b, c¯c and light quark (u¯u , d¯d , s¯s) initiated events produced in e + e − annihilations at LEP. The data were recorded with the OPAL detector at eleven diﬀerent energies above the Z 0 peak, corresponding to the full statistics collected at LEP1.5 and LEP2. The diﬀerence in mean charged particle multiplicities for b¯b and light quark events, δ bl , measured over this energy range is consistent with an energy independent behaviour, as predicted by QCD, but is inconsistent with the prediction of a more phenomenological approach which assumes that the multiplicity accompanying the decay of a heavy quark is independent of the quark mass itself. Our results, which can be combined into the single measurement δ bl = 3 . 44 ± 0 . 40 (stat) ± 0 . 89 (syst) at a luminosity weighted average centre-of-mass energy of 195 GeV, are also consistent with an energy independent behaviour as extrapolated from lower energy data.

p now at Physics Dept Southern Methodist University, Dallas, TX 75275, U.S.A. q now at IPHE Université de Lausanne, CH-1015 Lausanne, Switzerland r now at IEKP Universität Karlsruhe, Germany s now at Universitaire Instelling Antwerpen, Physics Department, B-2610 Antwerpen, Belgium t now at RWTH Aachen, Germany

Introduction
The study of heavy quark pair production in e + e − collisions with centre-of-mass energies greatly exceeding the heavy quark masses provides important tests of perturbative QCD. In particular mass effects are expected to induce substantial differences in the accompanying soft gluon radiation in heavy or light quark initiated events. A prediction of QCD concerns the difference in charged particle multiplicity, δ bl , between bb events and light quark ll (≡ uū, dd, ss) events, which is expected to be almost energy independent [1,2,3]. For a recent review see for example [4]. This prediction is in striking contrast with that from a more phenomenological approach, the so-called naïve model [5,6], which assumes that the hadron multiplicity accompanying the heavy hadrons in bb events is the same as the multiplicity in light quark events at the energy left to the system once the heavy quarks have fragmented. This naïve model predicts that δ bl decreases with increasing centre-of-mass energy.
Experimental results published at √ s = 91 GeV [7,8,9,10] and at lower energies [5,11], were not conclusive. Since the difference between the two theoretical predictions increases with increasing energy, a more powerful discrimination can be attempted at LEP2. Unfortunately the numbers of events collected at these higher energies are much smaller than at LEP1, and the results will therefore be much less precise.
In this analysis we have measured the mean charged particle multiplicity separately for b, c and uds initiated events observed with the OPAL detector at all eleven LEP energies above the Z 0 resonance, ranging from √ s = 130 GeV to √ s = 206 GeV, and derived at each energy the difference δ bl . This independent set of measurements covers fairly uniformly an energy interval of almost 80 GeV, and provides a clear discrimination between the two theoretical predictions.
The DELPHI Collaboration has published its first results [12] at three LEP2 energies, √ s = 183, 189 and 200 GeV. Their measured values of δ bl were found to be consistent with an energy independent extrapolation from lower energy data, and inconsistent with the prediction of the naïve model by more than three standard deviations.

Data sample and event simulation
The OPAL detector has been described in detail elsewhere [13]. The analysis presented here relies mainly on the reconstruction of charged particles in the central detector, which consisted of a silicon microvertex detector, a precision vertex drift chamber, a large volume jet chamber and drift chambers measuring the coordinate along the beam axis as they leave the jet chamber 1 . A solenoid providing a field of 0.435 T along the beam axis surrounded all tracking detectors. The analysis is based on data recorded with the OPAL detector between 1995 and 2000 at eleven different centre-of-mass energies, namely √ s = 130, 136, 161, 172, 183, 189, 192, 196, 200, 202 and 206 GeV. The data recorded in the year 2000 were mostly taken between 205 GeV and 207 GeV, with a weighted mean value of 206.1 GeV, and were analysed as a single energy. The data at 130 and 136 GeV were recorded in two different years, 1995 and 1997. We have checked, at both energies, that the two data sets give completely consistent results within the expected statistical uncertainties and therefore we have combined them. The integrated luminosity collected at the different centre-of-mass energies is detailed in Table 1.
High statistics samples of Monte Carlo events were generated at each energy to simulate the relevant physics process and the potential background. All generated events were passed through a detailed simulation of the OPAL detector [14] and processed using the same reconstruction and selection algorithms as the real data. To simulate signal events of the type e + e − → Z 0 /γ * → qq, we used the PYTHIA 5.7 parton shower model with fragmentation provided by the JETSET 7.4 routines [15] up to 189 GeV, and the PYTHIA 6.1 Monte Carlo program [16] interfaced with the KK2f program [17] at higher energies to obtain a more accurate description of initial state radiation. Both models have been tuned to describe the OPAL data at the Z 0 peak energy [18]. As an alternative fragmentation model we used events generated with HERWIG 6.2 [19] interfaced to KK2f, also tuned to the OPAL data 2 . We have observed that version 6.2 of HERWIG provides a substantial improvement with respect to version 5.9 in the description of the heavy quark sector, in particular the charged particle multiplicities. Both the PYTHIA and HERWIG models now provide an adequate description of multihadronic final states up to the highest LEP energies. Two-photon processes were simulated using PYTHIA, HERWIG and PHOJET [20], and τ pair production using KORALZ [21]. The 4-fermion background was studied using high statistics samples of Monte Carlo events generated with the grc4f 2.1 model [22], interfaced to JETSET 7.4 using the same parameter set for the parton shower, fragmentation and decays as mentioned above.

Event and track selection
The selection of non-radiative qq events was the same as in previous OPAL analyses [23,24]. In a first step, hadronic events were identified using criteria described in [25], optimised for running at energies above the Z 0 peak. The efficiency for selecting non-radiative hadronic events is greater than 98%, where simulated non-radiative events are defined as those with an invariant mass at the generator level, not considering photons radiated from the initial state, within 1 GeV of the nominal centre-of-mass energy.
In addition, it was required that charged tracks had transverse momentum p T > 150 MeV/c with respect to the beam axis, a minimum number of 40 hits in the jet chamber, a maximum allowed distance of the point of closest approach to the collision point in the r − φ plane, d 0 , of 2 cm and that this point should lie within 25 cm of the origin in the z direction.
To ensure a good containment in the detector and to reject background events of the type γγ → qq and e + e − → τ + τ − , we required that the polar angle of the thrust axis θ T , computed using charged tracks which passed the above mentioned criteria, satisfied the condition | cos θ T | < 0.9 and that there were at least seven accepted tracks. The residual background from these two processes was estimated to be less than 0.3% and consequently neglected. To reject events with large initial-state radiation we determined the effective centre-of-mass energy of the observed hadronic system, √ s ′ , as described in [26]. Events were rejected if √ s ′ < √ s − 10 GeV, where √ s is the nominal centre-of-mass energy. Above 160 GeV the 4-fermion background becomes significant and is dominated by hadronic decays of W pairs. This background was reduced to a manageble level by cutting on the QCD event weight variable, W QCD , as in [23]. This variable tests the compatibility of the events with QCD-like processes and details of its definition and performance for 4-fermion background rejection can be found in [24]. We accepted events with W QCD ≥ −0.5. After all cuts we found an overall efficiency for non-radiative qq events of about 78%. The residual fraction of events with a true effective centre-of-mass energy below √ s − 10 GeV is about 5%. The estimated residual 4-fermion background varies from 2.4%, at 161 GeV, to 11.6% at 206 GeV and was subtracted directly from the observed distributions, as described in the next section. The number of events selected at each energy after all cuts is presented in Table 1.

Experimental method and results
As in previous OPAL studies at LEP1 energies [8] and as proposed for LEP2 energies in [27], we used a method based on the simultaneous analysis of event samples with different quark flavour compositions to extract charged particle multiplicities separately for each flavor. At each centre-of-mass energy we selected three independent samples, one highly enriched in uū, dd, ss events (Sample 1), one slightly enriched in cc events with respect to an inclusive sample (Sample 2) and one highly enriched in bb events (Sample 3). The selection of these three samples was made using a well tested b-tagging technique [28] which uses information from b hadron lifetime, transverse momentum of leptons with respect to the jet axis, and kinematic observables. This b-tag uses artificial neural networks (ANNs) to combine optimally lifetime-sensitive tagging variables and also to combine kinematic variables in the jet-kinematics part of the tagging. The outputs of the lifetime ANN, the kinematic ANN and the lepton tag are finally combined, by using an unbinned likelihood method, into a single-valued variable L, which reflects the likelihood that a multihadronic event originates from a bb pair. The distribution of the event likelihood L at the two energies with highest statistics is shown in Figure 1. In the same figure the contributions from the different quark flavours and the residual 4-fermion background, as predicted from fully simulated events, are also shown. We show in Figure 2 the bb and uū, dd, ss event purities and efficiencies as a function of the cut on the L tagging variable, obtained from Monte Carlo events generated with PYTHIA at 189 GeV. The purity for b quarks at a given value X of L is defined as the fraction of genuine bb events with L ≥ X with respect to all events tagged with L ≥ X. Similarly, for uds quarks the purity at a given value Y of L is defined as the fraction of genuine uū, dd, ss events tagged with L ≤ Y with respect to all events tagged with L ≤ Y. The efficiencies are defined as the fractions of b or uds events tagged with that particular cut value of L with respect to the total number of produced b or uds events.
We used L ≥ 0.80 to select samples enriched in bb events (Sample 3), and L ≤ 0.05 to select samples enriched in uū, dd, ss events (Sample 1). The remaining events, namely those with 0.05 < L < 0.80, comprise Sample 2. The two vertical lines shown in Figure 1 show the cuts used to define the three samples. We could not select a sample highly enriched in cc events since the L variable is not sufficiently sensitive to the c quark fragmentation properties, and the limited statistics available at each energy does not allow direct c quark tagging using exclusive reconstruction of charm mesons, as was done at the Z 0 peak [8]. The low c quark purity of this sample translates into large uncertainties in the measurements ofn cc . Sample 2 contains, however, a slightly higher fraction of cc events compared to an inclusive sample, (30-33% against 22-26%, depending on the exact centre-of-mass energy) and has the advantage of being completely independent of the other two samples.
At a given energy, and after the subtraction of the residual 4-fermion background, the mean charged particle multiplicity measured in each sample,n corr i (i = 1, 3), corrected for detector effects, event selection cuts, residual contamination of radiative events and biases introduced by the tagging procedure, is a linear combination of the unknownsn bb ,n cc andn ll (l = u,d,s), the true mean multiplicities of the corresponding qq events. One can extractn bb ,n cc andn ll by solving the system of equationsn n corr are the fractions of bb, cc and uū, dd, ss events in the i th sample, evaluated from fully simulated e + e − → Z 0 /γ * → qq events. The experimentally corrected mean charged particle multiplicity is defined [29] as the total number of all promptly produced stable charged particles and those produced in the decay of particles with lifetimes shorter than 3 · 10 −10 sec. For each sample, the multiplicity distribution of the residual 4-fermion background after all selection cuts, normalised to the data integrated luminosity, was estimated from Monte Carlo events and directly subtracted from the experimentally measured distribution. We show in Table 1 the number of events left after all cuts and background subtraction, together with the corresponding 4-fermion background which was subtracted. The fractions of bb, cc and uū, dd, ss events as predicted by PYTHIA are also shown. The fractions predicted by HERWIG typically differ by around 1%.
To calculate the values ofn corr i , the mean valuesn obs i of the measured distributions were corrected for detector effects, event selection cuts, biases introduced by the b-tagging technique and by the residual contamination of radiative events by applying a multiplicative correction factor, C i , wheren obs i is the observed uncorrected mean charged particle multiplicity measured for the i th sample in the data andn MC−obs i is the same quantity obtained from a high statistics sample of fully simulated events.n MC−had i is the mean charged particle multiplicity obtained from a Monte Carlo sample with the quark flavour fractions as expected in sample i, without detector simulation and without simulation of the initial state radiation process. We have checked that the accuracy and the precision on the determination of the mean values obtained by using the correction method based on a simple multiplicative factor are completely equivalent to those obtained by a full matrix unfolding [30].
A reliable correction requires a good simulation of the data. As an example, we show in Figure 3 the observed charged particle multiplicity distributions, background subtracted, as measured in the three samples at 189 GeV. In the same figure we also show the corresponding distributions obtained from fully simulated Monte Carlo events generated with PYTHIA and HERWIG, normalised to the data integrated luminosity. Qualitatively both models reproduce the measured distributions reasonably well at this energy and at all the other energies considered in this paper.
The agreement between data and models is also shown in Figure 4, where the mean charged particle multiplicities as determined from the observed distributions after background subtraction are compared to those predicted by PYTHIA and HERWIG at three different energies. The errors are statistical only. Again, within the statistical uncertainties, we observed a satisfactory agreement between data and predictions from both models at all centre-of-mass energies. We therefore decided to correct our measured mean values using the coefficients C i computed, separately, with PYTHIA and HERWIG and quote the average of the two results as our reference value. For PYTHIA the value of these coefficients varies from 1.19 to 1.22 in Sample 1, from 1.05 to 1.10 in Sample 2 and from 1.08 to 1.10 in Sample 3, depending on energy. The difference between these values and those predicted by HERWIG is typically less than 1% and never exceeds 2%.
The system of equations (1) was then solved, at each energy, using the flavour fractions, f i , predicted by the corresponding model. The average between the two sets of results is presented in Table 2 and defines our reference values. The first uncertainty is statistical and the second is systematic. The statistical uncertainties on δ bl =n bb −n ll were computed taking into account the correlations between the measurements ofn bb andn ll . The correlation coefficients are positive and vary from 0.30 to 0.44, generally increasing with increasing energy. As anticipated, the large uncertainties on the measured mean multiplicities of the cc events reflect our experimental inability to efficiently select c quark enriched samples.

Systematic uncertainties
Several sources of possible systematic effects were considered.
• We have investigated a possible bias induced by the event selection. In particular an efficient suppression of the 4-fermion background was achieved by cutting on the W QCD variable which, however, also removes some qq events. We repeated the analysis using an alternative and independent cut. By rejecting events with a thrust value T < 0.83 we obtained the same level of background suppression. Half of the difference between the reference and the varied results was taken as a systematic uncertainty.
• The subtraction of the residual 4-fermion background relies on cross sections and charged particle multiplicities as predicted by the grc4f/JETSET models. By varying the predicted amount of background to be subtracted up and down by 5%, slightly more than its measured uncertainty at √ s = 189 GeV of 4% [31], we have checked that the differences with respect to the reference results are negligible. In previous OPAL analyses at 189 and 183 GeV [32], it was demonstrated that the mean charged particle multiplicity of hadronic W decays measured in the data and that predicted by the most commonly used hadronization models agreed within 1.1 times the total experimental error, corresponding to ±0.44 on the multiplicity. We varied the predicted mean charged particle multiplicity of the background by that amount, and take the difference between the reference and the varied results as a systematic error on our measurements.
• We have tested the stability of our results with respect to variations of the cuts on the variable L which determine the high b quark purity in Sample 3 and the high uds quark purity in Sample 1.
The analysis was repeated, selecting Sample 3 with cuts at 0.7 and 0.9, or selecting Sample 1 with cuts at 0.01 and 0.2, which lead to total variations of about 10% in the absolute b or uds quark purities of the samples. For each case, the magnitude of the larger variation is taken to be the systematic trend. In order to reduce statistical fluctuations in the estimate of the systematic uncertainty, we take a statistically weighted average value as a common systematic error at all energies.
• We have studied the systematic error associated with the simulation of the track resolution and its potential effect on our analysis through a change in performance of the b-tagging algorithms. A conservative estimate of this uncertainty was made by applying a smearing factor of 1.10 to the reconstructed r − φ and z projections of the track's impact parameter and the polar and azimuthal angles, in the PYTHIA samples used for the standard analysis [33]. The differences between the results with and without applying the smearing factor were taken as the systematic error.
• Uncertainties associated with the modelling of b hadron production and decay in the simulation affect the predicted efficiencies and purities of the b-tagging procedure. We have changed the ǫ b parameter in the fragmentation function of Peterson et al., so as to vary the average scaled energy < x E > of b hadrons in the range < x E >= 0.702 ± 0.008, as suggested by the LEP electroweak working group [34]. The lifetimes of b mesons and baryons were varied by ±0.02 ps and ±0.05 ps, respectively, based on the uncertainties on the measured values [35]. Finally, the average charged decay multiplicity of b hadrons was varied by ±0.062, reflecting the accuracy of the measurements by LEP experiments [34]. The three effects considered were treated independently and for each of them the largest difference between the reference result and the varied result was taken as a systematic uncertainty. At each energy the three uncertainties were added quadratically.
• To test the dependence of our results on the Monte Carlo model used to correct our data, we performed the analysis using, separately, both PYTHIA and HERWIG generators. As already mentioned in the previous section, our reference results were taken as the average between the two sets of results, and we assign half of the difference as systematic uncertainty.
The separate contributions to the systematic error on δ bl are summarized in Table 3. The total systematic uncertainty at a given energy was evaluated by adding all sources in quadrature. In the last column of Table 3 we show the systematic uncertainty of each source averaged over all energies. The errors from each source were assumed to be completely correlated at different energies and were weighted by the corresponding statistical uncertainty. The total averaged systematic error was again evaluated by adding all sources in quadrature.

Comparison with QCD predictions and models
According to perturbative QCD, soft gluon radiation from an energetic massive quark Q is suppressed inside a forward cone of half angle aperture Θ 0 = M Q /E Q , the so-called Dead Cone [36]. Here M Q is the heavy quark mass and E Q its energy, and the relation holds if E Q ≫ M Q ≫ Λ, where Λ represents the energy scale at which perturbation theory breaks down. This phenomenon produces significant differences in the structure of the soft gluon radiation emitted in light and heavy quark initiated jets. Assuming the validity of the Local Parton Hadron Duality (LPHD) concept a corresponding difference is also expected to be present at the hadronic level. For e + e − annihilations a QCD calculation within the Modified Leading Log Approximation (MLLA) [37], assuming a b quark mass of 4.8 GeV/c 2 , predicts a difference in mean charged particle multiplicity between bb and ll events of δ bl = 5.5 ± 0.8, independent of energy [1]. The quoted uncertainty is of experimental origin, while the uncertainty due to (energy-independent) missing higher order corrections is estimated to be about one unit.
More recently, the result of an improved calculation was published [2] in terms of a conservative upper bound for δ bl , which was found to lie in the range 3.7 to 4.1 depending on the b quark mass, m b , assumed to be between 5.3 and 4.7 GeV, respectively. In the same publication there was an attempt to evaluate more precisely the value of δ bl at √ s = 91 GeV, which gave δ bl = 3.68 for m b = 4.8 GeV/c 2 . However, there is no general consensus [4,38] on the theoretical consistency of the approach followed in [2], and it is still unclear whether a real improvement of the MLLA prediction has been achieved.
An independent upper limit of δ bl < 4 was obtained from phenomenological arguments and published in [3].
A more phenomenological approach, the naïve model [5,6], assumes that in heavy quark initiated events the multiplicity of light hadrons produced along with the heavy hadrons is the same as the total multiplicity of light quark initiated events produced at a centre-of-mass energy corresponding to the energy left behind after the heavy quarks have fragmented. In this framework one expects the value of δ bl to decrease with increasing energy. There are several variations of this model which lead, however, to only slightly different predictions. We have used a form where N ( √ s) = 2.554 + 0.1252 × exp(2.317 ln √ s) is a parameterization of the world mean charged particle multiplicity data, corrected to remove the effects of heavy quark production [39], x Q and x Q are the fractions of the beam energy carried by the heavy hadrons, N decay Q is the decay mean multiplicity of the heavy hadrons, and √ s is the centre-of-mass energy. We approximated the fragmentation function f (x Q ) for b quarks by a normalized Peterson function with a mean of 0.70, a conservative uncertainty of ±0.02 and assumed 2N decay Q = 11.0 ± 0.2 [1]. In Figure 5 we show our results on δ bl as a function of energy, together with all previously published results [5,7,8,9,10,11,12]. Statistical and systematic uncertainties have been added in quadrature.
In the same figure, the theoretical predictions are also shown. The striking difference between the QCD predictions (shaded area [1] and cross-hatched area [2]) and the prediction of the naïve model (single hatched area) is particularly evident at the highest LEP2 energies. One can see that the previously published results at the Z 0 peak energy and below did not allow a clear discrimination between the models. Overall they are consistent with an energy independent behaviour, but the naïve model could not be ruled out. The results published by the DELPHI Collaboration at three different LEP2 energies [12] showed a clear inconsistency with the predictions of the naïve model.
Our new results are consistent with those published by DELPHI [12], cover a much wider energy range and provide even stronger evidence of the inadequacy of this model.
A linear fit to our eleven data points, considering only their statistical uncertainties, yields a slope of 0.000 ± 0.018 (χ 2 /dof = 0.59), completely consistent with the QCD prediction of energy independence. Repeating the fit to a constant value or, equivalently, combining our results at a luminosity weighted average energy of 195 GeV, gives δ bl = 3.44±0.40 (stat) ±0.89 (syst). The overall systematic uncertainty of 0.89 was calculated assuming that each source of systematic uncertainty is fully correlated between energy points (see Table 3). This value is consistent with the published OPAL result at 91 GeV [8] of δ bl = 2.79 ± 0.30 (total error) and with the value of 2.99 ± 0.20 (χ 2 /dof = 0.79) obtained from the corresponding fit to all published results up to and including 91 GeV, assuming that the measurements are completely uncorrelated. A weighted average including results from low energy data, LEP1 and LEP2 gives δ bl = 3.05 ± 0.19, which is shown in Figure 5 as a dash-dotted line. In this average we have assumed that the systematic errors of the DELPHI measurements at LEP2 are completely correlated between energy points.
One can also see from Figure 5 that the MLLA+LPHD prediction [1] of 5.5 ± 0.8 (exp) is higher than the experimental results, even considering the additional theoretical uncertainty of about 1 unit due to missing higher order corrections. The upper bounds calculated in [2] and in [3] are consistent with the measurements.

Conclusions
We have measured the mean charged particle multiplicities for bb, cc and uū, dd, ss events at all energies collected by OPAL above the Z 0 peak, and in particular we have determined the differences between the mean multiplicity of b and uds initiated events, δ bl =n bb −n ll .
Our results are presented in Table 2 and are in agreement with previously published results [12] at three LEP2 energies. Our data alone, which fairly uniformly cover a wide energy range, strongly support the QCD prediction of the energy independence of δ bl , leading to a combined result of δ bl = 3.44 ± 0.40(stat) ± 0.89(syst) at a luminosity weighted average centre-of-mass energy of 195 GeV. The consistency of the experimental results over the entire range from √ s = 29 to √ s = 206 GeV strengthens this conclusion even further. The naïve model, which assumes that the multiplicity accompanying the decay of a heavy quark is independent of the mass of the quark itself, is strongly disfavoured.  Table 1: The integrated luminosity, L int , collected at each energy and the total number of events after all selection cuts are shown in the first three columns. The number of events in each sample after the subtraction of the residual 4-fermion background, the fraction of events due to background which was subtracted (in parentheses) and the percentage of bb, cc and uū, dd, ss events as predicted by PYTHIA are also shown.    Table 3: Contributions from different sources to the systematic uncertainty on the δ bl measurements.
In the last column we show the average of the systematic uncertainties over all energies. The errors from each source were assumed to be fully correlated at different energies, and were weighted by the corresponding statistical uncertainty. The total averaged uncertainty was obtained by adding all the contributions in quadrature.

22.
24.  Average measurement Figure 5: The difference in mean charged particle multiplicity between bb and uū, dd, ss events, δ bl , as a function of centre-of-mass energy. The data points show the experimental measurements and the total error, and those around √ s = 91, 183, 189 and 200 GeV have been separated horizontally for clarity. The original MLLA prediction [1] is shown as a shaded area to include the errors of experimental origin on this prediction, not including missing higher order corrections. The crosshatched area corresponds to the QCD upper limits as calculated in [2]. The single hatched area represents the naïve model prediction [5,6], while the dash-dotted line is the combined result from all the measurements, as discussed in section 6.