Study of the fragmentation of b quarks into B mesons at the Z peak

The fragmentation of b quarks into B mesons is studied with four million hadronic Z decays collected by the ALEPH experiment during the years 1991-1995. A semi-exclusive reconstruction of B->l nu D(*) decays is performed, by combining lepton candidates with fully reconstructed D(*) mesons while the neutrino energy is estimated from the missing energy of the event. The mean value of xewd, the energy of the weakly-decaying B meson normalised to the beam energy, is found to be mxewd = 0.716 +- 0.006 (stat) +- 0.006 (syst) using a model-independent method; the corresponding value for the energy of the leading B meson is mxel = 0.736 +- 0.006 (stat) +- 0.006 (syst). The reconstructed spectra are compared with different fragmentation models.


Introduction
The process of hadron production at e + e − colliders is usually modelled as the convolution of a perturbative part (hard gluon radiation for energies above approximately 1 GeV) and a non-perturbative part, called hadronisation or fragmentation, in which the quarks are confined in colourless hadrons. While the first step is in principle calculable, the fragmentation needs a phenomenological approach and is usually parametrised in terms of the variable where p is the hadron's momentum along the direction of the quark, and (E + p) quark is the sum of the quark energy and momentum just before fragmentation, i.e. taking into account initial and final state photon radiation, and hard gluon emission. With this definition, the fragmentation process can be described in terms of the probability of a hadron H to be generated with a given z, called D H q (z), where q is the flavour of the generating quark. In this paper the fragmentation of b quarks is studied.
The fraction z is not accessible experimentally, and hence a direct reconstruction of D H b (z) is not possible. The energy spectrum of b hadrons is therefore described in terms of the scaled energy x b , defined as the ratio of the heavy hadron energy to the beam energy In contrast to the z variable, the effects of initial and final state radiation and hard gluon emission are not unfolded.
In the analysis presented, the energy of B mesons is reconstructed using a partially exclusive method: semileptonic decays B → ℓνD (⋆) are identified by pairing lepton candidates with fully reconstructed D (⋆) mesons; the scaled energy of the weakly-decaying B meson is then computed adding an estimate of the neutrino energy. Five channels are chosen because of their good signal purity and statistical significance; they are shown in Table 1. In the following x obs B indicates the reconstructed energy of B meson candidates, x wd B the energy of weakly decaying B mesons, corrected for detector acceptance and resolution; x L B stands for the corrected scaled energy of the leading B meson, that is the first meson produced in the fragmentation process, which can also be a heavier resonance (B ⋆ , B ⋆⋆ ).

Channel
Decay chain The analysis uses the full LEP I statistics collected by ALEPH between 1991 and 1995, amounting to almost four million hadronic Z decays. Recently this data set has been reprocessed using improved reconstruction algorithms. The main benefits for this analysis are related to the enhanced secondary vertex reconstruction efficiency and the improved particle identification. A discussion of the reprocessing can be found in [1].
After a description of the ALEPH detector, the selection of B → ℓνD (⋆) decays is detailed in Section 3. In Section 4 the reconstruction of the B meson energy is described, followed by the extraction of the spectrum and comparison with the predictions of different models in Section 5. Systematic errors are discussed in Section 6, and checks on the selfconsistency and robustness of the analysis are presented in Section 7.

The ALEPH detector
The ALEPH detector and its performance are described in detail elsewhere [2,3]. A high resolution vertex detector (VDET) consisting of two layers of silicon with double-sided readout measures rφ and z coordinates at average radii of 6.5 cm and 11.3 cm, with 12 µm resolution at normal incidence. The VDET provides full azimuthal coverage, and polar angle coverage to | cos θ| < 0.85 for the inner layer and | cos θ| < 0.69 for both layers. Outside VDET, particles traverse the inner tracking chamber (ITC) and the time projection chamber (TPC). The ITC is a cylindrical drift chamber with eight axial wire layers with radii between 16 and 26 cm. The TPC measures up to 21 space points per track at radii between 40 and 171 cm, and also provides a measurement of the specific ionization energy loss (dE/dx) of each charged track. These three detectors form the tracking system, which is immersed in a 1.5 T axial magnetic field provided by a super-conducting solenoid. The combined tracking system yields a transverse-momentum resolution of σ(p T )/p T = 6 × 10 −4 p T (GeV/c) ⊕ 0.005. The resolution of the three-dimensional impact parameter for tracks having two VDET hits can be parametrised as σ = 25 µm + 95µm/p, (p in GeV/c).
The electromagnetic calorimeter (ECAL) is a lead/wire chamber sandwich operated in proportional mode. The calorimeter is read out in projective towers that subtend typically 0.9 • × 0.9 • , segmented in three longitudinal sections. The hadron calorimeter (HCAL) uses the iron return yoke as absorber. Hadronic showers are sampled by 23 planes of streamer tubes, with analogue projective tower and digital hit pattern readouts. The HCAL is used in combination with two double layers of muon chambers outside the magnet for muon identification.

Selection of B → ℓνD (⋆) decays
A Monte Carlo simulation based on JETSET 7.4 [4] and tuned to ALEPH data [5,6] has been used in order to extract resolution functions, acceptance corrections and background compositions. About five million bb events were simulated, and more than twice the data statistics of qq events. The present analysis uses bb events to determine the x wd B and x L B spectrum starting from observed x obs B spectra, and qq events to evaluate the non-bb component of the selected sample.
The decays B → ℓνD (⋆) are searched for in hadronic events, containing at least one lepton (electron or muon) identified using standard criteria [7]. The momentum cut used to define lepton candidates is p > 2 GeV/c for electrons and p > 2.5 GeV/c for muons. The transverse momentum p T of the lepton with respect to the nearest jet, with the lepton excluded from the jet, is required to be larger than 1 GeV/c, which helps rejecting fake candidates and leptons not coming from direct decays of b hadrons. Both electron and muon candidates are required to have a measured dE/dx compatible with the expected value.
Events are divided into two hemispheres by a plane perpendicular to the thrust axis; in each hemisphere containing a lepton a D meson reconstruction is attempted in the decay modes described in Table 1. At least two charged tracks from the D meson decay are required to have VDET hits, in order to ensure a good reconstruction of the D vertex position and to reject combinatorial background. Loose cuts are applied to track momenta, in order to minimise the bias in the B momentum distribution. Tracks are not considered as kaon candidates if their measured ionization is incompatible with the kaon hypothesis by more than three standard deviations. The charge of the kaon candidate is required to be the same as that of the lepton, as expected for semileptonic B meson decays.
Tracks assigned to a D meson decay are fitted to a common vertex, and the track combination is rejected if the χ 2 of the fit is larger than 20. If more than one combination fulfils this requirement for channels 3 and 4, the one with the smallest χ 2 is chosen. In channel 5, the π 0 closest in angle to the charged pion is selected and added to form the D 0 .
For channels 1, 2 and 5, a soft pion π s is added to the D candidate to form a D ⋆ meson; the π s momentum is required to be larger than 250 MeV/c and smaller than 3 GeV/c. The difference between the reconstructed D ⋆ and D masses is required to be within 5 MeV/c 2 of its nominal value. In the case of multiple candidates in a given hemisphere, the track combination is chosen for which the reconstructed (D ⋆ − D) mass difference is closest to the nominal value.
A vertex fit is performed using the D candidate and the lepton track, and again the combination of tracks is rejected if the χ 2 of the fit is larger than 20. The B vertex is required to lie between the interaction point, reconstructed event-by-event, and the D vertex.
Channels 3 and 4 are further enriched in signal events using harder cuts on the kaon; in addition a π s veto is applied: if a track is found which is compatible with the reconstructed B vertex and the combination track-D candidate has a mass close to the mass of the D ⋆ , the candidate is discarded. This procedure reduces the overlap between channels at the permil level. Finally, tighter cuts on the reconstructed D mass and the χ 2 of the vertex fit are imposed.
The D mass spectra are shown in Fig. 1. The reconstructed D mass peaks are fitted in a region between 1.7 and 2.0 GeV/c 2 with a Gaussian and a linear component. Table 2 shows the chosen D mass windows, the number of reconstructed candidates and the fitted Gaussian fractions.
The fractions of the Gaussian components measured in the Monte Carlo are compatible  with those in the data within statistical errors, while the widths are about 5−10% smaller; this is taken into account in the evaluation of the systematic uncertainties.

B energy reconstruction
The scaled energy of the weakly-decaying B → ℓνD (⋆) hadron is estimated as The terms E D (⋆) and E ℓ are provided by the direct reconstruction, while the neutrino energy E ν is estimated from the missing energy in the hemisphere: where E hemi tot is estimated taking into account the measured mass in both hemispheres [8]: Both charged and neutral particles are used in Eqns. (4) and (5). In the lepton hemisphere neutral hadronic energy is expected to come only from fragmentation. Therefore, in order to avoid spurious calorimetric fluctuations, its contribution is taken into account only outside a cone of 10 degrees of half opening angle around each of the B meson decay products. Table 3 shows the resolution on x B estimated on simulated bb events; the distributions are well described by two Gaussians, accounting for core and tails.
Channel Table 3: For the five channels the x B resolution on simulated events. The resolution can be parametrised with two Gaussians, describing the core and the tails.

Unfolding methods
The scaled energy of the weakly-decaying B mesons, x obs B , is reconstructed in 20 bins between 0 and 1 with a variable width. In each bin, the non-bb background is estimated using the simulation, and subtracted from the spectrum. This amounts to about 2% of the events, concentrated mostly at low x obs B . The measured spectra after subtraction are shown in Fig. 2.
With these events two different kinds of analyses can be performed: • a model-dependent analysis, in which different fragmentation models available in the literature are tuned to fit the observed spectra; • a model-independent analysis, in which the shapes of x wd B and x L B are reconstructed by correcting the observed spectra for detector acceptance, resolution and missing particles.

Model-dependent analysis
Various fragmentation functions D H b (z) are implemented in the Monte Carlo generator JETSET 7.4, which also simulates initial and final state photon radiation and hard gluon emission. The reconstructed spectra obtained from the simulated bb events are tuned to best reproduce the x obs B distribution observed in the data, by minimising the global χ 2 . The following parametrisations for D H b (z) are used: Kartvelishvili et al. [10] : Collins et al. [11] : The minimisation is performed with respect to the free parameter of each model, and the χ 2 is written as: where c runs over the channels used, i runs over the x obs B bins defined as in Table 6, n DT and n MC are the number of candidates per channel and per bin observed in the data and expected from the Monte Carlo, normalised to the same number of entries. The quantities σ are defined as statistical uncertainties. Table 4 shows the fitted values for the different model parameters, together with statistical and systematic uncertainties from the sources discussed in Section 6. Also shown are the values for the mean scaled energy.

Model
Fit results Mean energies Peterson Table 4: Fit results with different fragmentation models. The systematic errors account for the sources of uncertainties discussed in Section 6. The χ 2 /N DOF is calculated using statistical errors only.
The Kartvelishvili model describes the data slightly better than the Peterson model. The Collins model is clearly disfavoured.

Model-independent analysis
The x wd B and x L B spectra are obtained by correcting the observed x obs B spectra for acceptance, detector resolution and missing particles.
The normalised binned spectrum f i x wd B can be obtained using the relation where ǫ wd i (c) is the acceptance correction in bin i for channel c; n DT j (c) is the number of reconstructed B mesons in the data for channel c, with a measured energy falling in bin j; G wd ij (c) is the resolution matrix that links mesons with x obs B in bin j and x wd B in bin i, for channel c; T is the normalisation factor defined by the condition i f i = 1. A similar equation holds for the extraction of f i x L B , where the effect of the missing particles from B ⋆ and B ⋆⋆ decays is folded in ǫ L and G L .
The acceptance corrections ǫ i and the resolution matrix G ij are taken from the simulation. The acceptance corrections show a different behaviour among the different channels; the two extreme situations are shown in Fig. 3. A dependence on the fragmentation function present in Monte Carlo is induced in the measured spectrum through G ij ; hence, the Monte Carlo used to calculate G ij must be reweighted to the best estimate of f i x wd The whole procedure is then repeated until the change in f i in consecutive iterations is a small fraction of the statistical errors.
The statistical error matrix E ij is calculated by repeating 20 × 5 analyses, varying in each of them the quantities n DT j (c) by one standard deviation: where f (ck) i is the result of the convergence for f i when n DT k (c) is varied by its statistical error, and f STD i the nominal result. The results, together with the statistical and systematic errors, are given in Table 6. The full error matrices are shown in the Appendix.
From the binned spectra f i x wd B the mean value is calculated as where x i is the central value of bin i. The bin size chosen is such that the deviation from linearity of the distribution within a bin is negligible. The statistical error on x B is calculated using the same procedure used as for  Table 5: Exclusive branching ratios for the B → ℓνX process [12,13]. The sum is consistent with the measurement of the inclusive B → ℓνX rate.

Systematic errors
Possible systematic effects due to uncertainties on the physics parameters used in the Monte Carlo, limited accuracy in the simulation of the detector response, or effects intrinsic to the analysis method have been investigated.
The physics parameters used in the Monte Carlo simulation that are relevant for the analysis are adjusted to the most recent experimental measurements and varied within their estimated uncertainty by reweighting simulated events. The effect of the reweighting propagates to the results through the resolution matrix G ij (c) and the acceptance corrections ǫ i (c), which are taken from the simulation. The differences from the standard results are taken as systematic errors.
The sources of uncertainty considered are: • Semileptonic decays of B mesons.
The current experimental knowledge [12,13] of the semileptonic branching ratios of B mesons is summarised in Table 5. The sum of the exclusive (or semi-exclusive) rates is consistent within errors with the inclusive measurement of BR(B → ℓνX).
The analysis is not sensitive to the total BR(B → ℓνX), but is affected by a change in the relative rates of the different components, since these contribute in a different way to the average acceptance corrections and resolution matrix.
Six sources of systematic error are calculated using the values in Table 5: When deriving the energy spectrum of the leading B meson, the correction due the energy carried away by the pion produced in the B ⋆⋆ decay enters in the resolution matrix and the acceptance corrections. The rate of b → B ⋆⋆ is varied within its experimental error: f B ⋆⋆ = 0.299 ± 0.058 [14,15,16], and the resulting systematic error is ∆ x L B = 0.0025. The weakly-decaying B meson spectrum is not affected by this source of uncertainty.
• Modelling of the B ⋆⋆ production.
• Production of B ⋆ from b quarks.
Due to the small mass difference between B ⋆ and B mesons, the effect of B ⋆ production in b quark fragmentation is found to be much smaller than for B ⋆⋆ , and it is completely negligible for the present analysis.
The relevant sources of systematic uncertainties due to the detector simulation are identified to be: • Neutrino energy reconstruction.
The accuracy of the neutrino energy reconstruction is checked in hadronic events enriched in light primary quarks and in hadronic decays of τ τ events. In the first sample, a "fake" neutrino is simulated by removing a charged particle from the reconstructed event; its energy is then reconstructed using Eqns. (4) and (5). The method can be applied both to data and Monte Carlo events, determining the bias between the reconstructed energy and the momentum measured with the tracking system. Such a bias is found to be reproduced by the Monte Carlo with a precision better than 50 MeV, for all momenta of the deleted track. In the second sample, where the event topology is much simpler, the energy of the "reconstructed" ν τ is compared between data and Monte Carlo events. Also in this case the worst discrepancy observed is smaller than 50 MeV. This value is used as a conservative estimate of the systematic uncertainty on the neutrino energy, resulting in ∆ x wd B = 0.0023 and ∆ x L B = 0.0023.
If the purity and the kinematic properties of the selected candidates are not well described by the simulation, the acceptance corrections and resolution matrices can be inadequate. In order to check for these effects, the distributions of the χ 2 probability for the reconstructed D vertices are compared, channel by channel, with the simulation. Small differences are observed, and the Monte Carlo distribution is reweighted in order to reproduce the data. The shift in the corrected average energy is taken as systematic uncertainty. The resulting error estimates are ∆ x wd B = 0.0001 and ∆ x L B = 0.0001. Furthermore, the reconstructed D mass distributions in data and Monte Carlo are compared. In simulated events the widths of the mass spectra are found to be 5 − 10% smaller, while the fractions of the Gaussian components, estimated from a fit to the sidebands, are reproduced within their statistical error of about 5%. The mass cuts reported in Table 2 are adjusted in order to take into account both effects, taking the total shift in the extracted energy spectrum as systematic uncertainty. The resulting estimates are Possible systematic effects related to the analysis procedure are: • Background subtraction.
As previously explained, a bin-by-bin subtraction of candidates not coming from bb events is performed before deriving the B meson energy spectra. The efficiency for this kind of background has been extracted directly from data events, which have been enriched in background events by selecting wrong sign candidates. It is found to be compatible with The bin-by-bin results for the measured spectra, with the total systematic uncertainties, are shown in Table 6, while the statistical and total error matrices are reported in the Appendix. The spectra are also shown in Fig. 4, where they are compared with the Monte Carlo predictions from different fragmentation models, with the free parameters fitted to the data.
The models of Peterson and Kartvelishvili give the best agreement with the data, and are compared with the x wd B measurement in Fig. 5.  Bin

Systematic checks
Possible systematic effects intrinsic to the analysis method are checked by measuring the energy spectra in a sample of 8 million simulated qq events. The average values for the scaled energies of the weakly-decaying and leading B mesons are measured to be: , which compare well with the true values: Electron and muon identification are affected by different sources of background, and the selection efficiencies and purities have a different dependence upon the track kinematics and isolation. It is therefore interesting to perform the analysis using separately events with electron candidates or muon candidates. Consistent results within uncorrelated errors are found: The acceptance corrections for the five channels are significantly different (Fig. 3) and the same is true for the resolution matrices. An inaccurate description of these inputs would easily lead to incompatible results among the different channels. This is checked by performing the analysis separately in the five sub-samples. The results, reported in Table 7, are compatible within uncorrelated uncertainties.

Channel
x wd B  It has been checked that the results are independent of the choice of fragmentation functions in the Monte Carlo sample used to estimate the resolution matrix G ij and the acceptances ǫ i .
As explained in Section 5.2, the weights applied to reweight G ij to a given fragmentation function are smoothed with a polynomial function to reduce the bin-to-bin fluctuations. However, the values for the mean scaled energies move by a small fraction of the statistical errors when such smoothing is not applied, and the total statistical errors remain nearly constant: Heavy flavoured hadrons originating from gluon splitting g → bb have an energy much lower than hadrons coming from primary b quarks. A check on Monte Carlo events shows that the contribution of such events is negligible.
The analysis uses a binned representation of the fragmentation functions to compensate the relatively small statistical sample in the data. The binning chosen must not introduce biases in the measured values nor should it affect the statistical errors. This is checked by performing a number of analyses in which the binning is varied randomly around the standard one. Both the central values and the statistical uncertainties are stable.

Conclusions
Using the data collected by the ALEPH experiment at and around the Z resonance in the years 1991-1995, about 3400 semileptonic B 0 and B ± decays have been selected. The scaled energy spectra of weakly-decaying and leading B mesons have been reconstructed, and their mean values were found to be: This measurement supersedes a previous analysis from ALEPH [18], which used a different method and smaller statistics.
The present result is compatible with the published results using b hadrons from L3 [19], OPAL [20] and SLD [21], and using B mesons from OPAL [22].

Acknowledgements
We wish to thank our colleagues in the CERN accelerator divisions for the successful operation of LEP. We are indebted to the engineers and technicians in all our institutions for their contribution to the excellent performance of ALEPH. Those of us from nonmember countries thank CERN for its hospitality. Tables 8 and 9 show the statistical and total error matrices for x wd B ; Tables 10 and 11 give the same information for x L B .