Production rate of bb(bar) quark pairs from gluons and bb(bar)bb(bar) events in hadronic Z Decays

The rates are measured per hadronic Z decay for gluon splitting to bb(bar) quark pairs, g_bb, and of events containing two bb(bar) quark pairs, g_4b, using a sample of four-jet events selected from data collected with the OPAL detector. Events with an enhanced signal of gluon splitting to bb(bar) quarks are selected if two of the jets are close in phase-space and contain detached secondary vertices. For the event sample containing two bb(bar) quark pairs, three of the four jets are required to have a significantly detached secondary vertex. Information from the event topology is combined in a likelihood fit to extract the values of g_bb and g_4b, namely g_bb = (3.07 +- 0.53(stat) +- 0.97(syst))x10^-3 g_4b = (0.36 +- 0.17(stat) +- 0.27(syst))x10^-3


Introduction
Bottom quark pairs in Z 0 decays can be produced either directly via Z 0 → bb or indirectly, when a gluon is radiated from a quark and then splits into a bb quark pair. A special case consists of events with both direct and indirect b quark production, Z 0 → bbg → bbbb. The rates g bb and g 4b per hadronic Z 0 decay, for the reactions Z 0 → qqg with g → bb and Z 0 → bbbb are sensitive to both the b quark mass and the strong coupling constant α s . Hence measurements of these rates are tests of the theory of Quantum Chromodynamics (QCD). Events with gluon splitting into bb quark pairs are an important background for the measurement of R b , the fraction of hadronic Z 0 decays into bb quark pairs. Therefore, a more precise determination of g bb might lead to a reduction of the systematic uncertainty in R b .
The rate g bb has been calculated in [1,2], including the re-summation of leading logarithmic terms. Those authors point out that the parton-shower approach as implemented e.g. in Jetset 7.4 [3] is a good approximation. Numerical calculations of g bb are given in [1,2,4], predicting a rate in the range g theor bb = (1.8 − 2.9) × 10 −3 , depending on the b quark mass and the strong coupling constant.
An estimate of g 4b can be obtained from the rate for direct bb production, R b , multiplied by the rate of indirect bb production, g bb . This simple picture is modified because the phase-space for two bb quark pairs is smaller than that for two light and two b quarks. The interference between secondary bb production and primary bb production cancels to zero at leading order in α s , except for the bbbb final state [1,5]. In [1], the contribution of this interference term is shown to be less than 0.2% of g bb .
In this analysis, decays of the Z 0 into four-jet final states are investigated. Jets from b or b quarks are identified by reconstructing secondary decay vertices. The invariant mass of bb quark pairs originating from gluons tends to peak just above threshold. This leads to a small relative momentum of the two b hadrons produced in the fragmentation process. By contrast, directly produced bb quark pairs have high invariant masses, since they carry a large part of the Z 0 energy. This and other characteristics are used to select event samples enriched in the process g → bb. An angular correlation defined similarly to the Bengtsson-Zerwas angle [10] is used to further differentiate between qqqq and qqgg final states. In addition, events with three reconstructed secondary vertices are selected. They are used to measure the g 4b rate.
2 Event selection and reconstruction 2

.1 The OPAL detector
The OPAL detector is described in detail elsewhere [11]. Only a brief description of the detector elements relevant to this analysis is given here. Charged tracks are reconstructed in the central tracking system. It consists of a silicon microvertex detector, a vertex drift chamber equipped with axial and stereo wires, a large jet chamber and z chambers 1 . A solenoid providing a uniform magnetic field of 0.435 T parallel to the z-axis surrounds the central tracking system. The silicon microvertex detector [12] has two layers which measure tracks in (r, φ). This detector was upgraded in 1993 to provide a precise measurement of the z coordinate [13]. Before this detector was upgraded again in 1995 for the high energy operation in the LEP 2 programme [14] the inner layer covered the range | cos(θ)| < 0.83 and the outer layer the range | cos(θ)| < 0.77. The vertex chamber extends over the range | cos θ| < 0.95. The coil is surrounded by scintillators for time-of-flight measurements and a barrel lead-glass electromagnetic calorimeter. Including the endcap electromagnetic calorimeter, the lead-glass blocks cover the range | cos θ| < 0.98. The magnetic return yoke is instrumented with streamer tubes and serves as a hadron calorimeter. The return yoke in turn is surrounded by muon chambers.

Event selection and reconstruction
The analysis uses data taken with the OPAL detector in the years 1992-1995 on or near the Z 0 resonance. Hadronic Z 0 decays are selected with an efficiency of 98.4%, as described in [15]. Only events with the tracking system and the electromagnetic calorimeter fully operational are used in this analysis. A total number of 3.35 million hadronic events are selected. In these events well-measured charged tracks are used with a momentum p t > 0.15 GeV/c in the (x, y) plane, and clusters in the electromagnetic calorimeter with energies above 0.1 GeV (0.25 GeV) in the barrel (end-cap) region. The energies of clusters pointed to by charged tracks are corrected for double counting by subtracting the energy deposition expected from the track momentum [16].

Simulated events
A total number of 150 000 events generated with Pythia 6.130 [17] and 14 million events generated with Jetset 7.4 [3] are used to evaluate the efficiencies for signal and background, respectively. Within the Jetset 7.4 sample, 4.3 million (1.9 million) events were generated in special runs with a primary bb (cc) quark pair to increase the statistical significance of the description of background processes with secondary decay vertices. To study the g → bb signal process 150 000 events with gluon splitting to bb in the parton shower were generated with Pythia 6.130. This generator is used with a special option for gluon splitting to massive quarks, as explained in Section 4.1. All signal and background events were passed through a complete simulation of the OPAL detector [18].
The simulated events are weighted to correspond to the measured values of cc and bb production, R c = 0.1724 and R b = 0.21664, given in [19]. The rate of gluon splitting to cc pairs is set to g cc = 3.20 × 10 −2 as measured by OPAL [20]. The ratio of the number of events with primary produced b quarks and g → cc to the total number of events with g → cc is fixed to the value predicted by Jetset 7.4.
In the following, simulated events with gluon splitting to bb signal are referred to as g → bb.
They are further differentiated into events with a primary charm or light quark, referred to as qqbb (q = udsc), and events with two bb quark pairs, referred to as bbbb.
Background to this analysis consists of various sources of four jet events. The dominant fourjet process is the production of two gluons in addition to the primary quarks, either from double Bremsstrahlung or from the triple gluon vertex. About 7% of the four-jet events are expected to be from gluon-splitting to a quark antiquark pair [21].
As this analysis is based on the identification of secondary decay vertices, there are two sources of background which have to be studied more carefully: events with gluon splitting to cc accompanied by any flavor of primary quark, and four-jet events with a primary bb quark pair but without gluon splitting to heavy quarks. They are referred to as bbxx (x = guds).
For comparisons to the data, the rate of gluon splitting to bb pairs is set to g bb = 2.5 × 10 −3 and the rate of bbbb events is set to g 4b = 0.4 × 10 −3 . Note that the rate g bb includes bbbb events. The fit procedure, described in Section 3.3, is completely independent of these two numbers.

The four-jet selection
The four-momenta of the selected tracks and clusters are combined to form four jets, using the k ⊥ (Durham) algorithm [22]. The value y cut at which an event makes a transition between a three-jet and a four-jet assignment is denoted y 34 . Figure 1 shows the normalized distribution of the quantity y 34 for events that are selected as hadronic Z 0 decays, thereby comparing the data to the Monte Carlo prediction. In addition, the distribution predicted for g → bb events is shown, scaled by a factor of 400. These events populate mainly the region of high y 34 . A cut y 34 > 0.006 is made to define the four-jet sample. This cut rejects nearly 90% of the background events, while retaining about 60% of the signal events. The estimated signal fraction assuming a signal rate g bb = 2.5 × 10 −3 is 1.2%. In the data 443 334 events are selected, corresponding to 13.0% of all hadronic Z 0 decays. In the Monte Carlo simulation only 11.5% of the events are selected, because the prediction is slightly shifted to lower values of y 34 with respect to the data. This corresponds to a deficit of (11.7 ± 0.2 ± 0.8)% four-jet events in the simulation compared to the data, where the first uncertainty is due to the finite number of events in the data and the simulation, the second is from the uncertainty in the rate of gluon splitting to cc, g cc = (3.20 ± 0.43) × 10 −2 [20]. This deficit is attributed to missing higher orders in the parton shower simulation. To deal with this normalization problem, the number of simulated events is normalized to the number of four-jet events throughout this analysis, instead of normalizing to the number of hadronic Z 0 decays. The stability of this analysis with respect to variations of the cut in y 34 , affecting this normalization, is discussed in Section 5.

Vertex reconstruction
A primary vertex is reconstructed for each event as described in [23]. Then for each jet a secondary vertex is reconstructed, using the full three-dimensional tracking information [24]. All tracks in the jet satisfying additional quality cuts on the momentum p > 0.5 GeV/c, the distance to the fitted primary vertex in the r/φ plane d 0 < 0.3 cm and the uncertainty on d 0 of σ d 0 < 0.1 cm are fitted to a common vertex. The relatively high momentum cut results in an increased contribution of tracks from the decay of b hadrons. Tracks with a large χ 2 contribution are removed from the vertex and the fit is repeated until the χ 2 contribution of each track is smaller than 4. If three or more tracks remain, the secondary vertex is accepted.
The decay length l and its error σ are calculated as the distance of the primary to the secondary vertex. The direction is constrained to be parallel to the jet axis. The decay length l is positive if the angle between the jet axis and the vector pointing from the primary to the secondary vertex is less than 90 degrees, negative otherwise. Vertices with l > 0 are used to identify b hadrons. In this analysis, the following variables are used to identify secondary vertices originating from the decay of b hadrons: • The decay length significance l/σ, calculated from the decay length l and the experimental uncertainty σ on l. Secondary vertices significantly separated from the primary vertex are selected with a cut l/σ > 3.
• For the vertices surviving the l/σ > 3 cut the output NN of a neural network is calculated. This neural network has been trained to separate vertices originating from b hadron decays from those in charm or light quark events.
The neural network was developed for the OPAL R b analysis [24]. It has five inputs: the decay length significance l/σ, the decay length l, the number of tracks in the secondary vertex n s , the reduced decay length l r /σ r , where one well-defined track [24] has been removed from the vertex fit, and a variable x D , sensitive to the invariant mass of the tracks in the jet that have a high probability of originating from a b hadron decay. The neural network output NN lies between zero and one. Values close to one indicate a high probability that the vertex is associated with a b hadron decay.

Analysis
In Section 3.1 the selection of candidate events for gluon splitting is described. By changing the y-cut and exploiting the transition from four to three jets, two jets are selected in each event as candidates to have originated from gluon splitting. The candidate jets are checked for secondary vertices and events with two significant secondary vertices are selected. The event sample is subdivided into two distinct classes, depending on the event topology. Optimized cuts on the neural network outputs are applied for each class, to define the candidate events. In Section 3.2 a dedicated selection of candidates for the process Z 0 → bbbb is discussed, where all four jets are checked for secondary vertices. Finally, in Section 3.3 the rate of gluon splitting to bb is calculated. For each of the selected event samples, angular distributions sensitive to four-quark final states are studied. The rates g bb and g 4b are calculated using a binned maximum likelihood fit. The signal and background selection efficiencies and these angular distributions are used as input to the fit.

The qqbb event selection
In each four-jet event the y-cut is increased, until the event changes to a three-jet event (y cut > y 34 ). The two jets that are combined in this step are considered as candidates for g → bb. The y-cut is then increased further, until the event changes to a two-jet event (y cut > y 23 ). There are two distinct possibilities for this, as shown in Figure 2.
Class "2+2": Events belong to this class if none of the original four jets is identical to any of the two jets obtained after increasing the y-cut to force the event into a two-jet event.
Class "3+1": Events belong to this class if one of the original four jets is identical to one of the two jets obtained after increasing the y-cut to force the event into a two-jet event.  Table 1: Events selected at the different steps of the analysis (data), and the number of background (bgnd) and signal events expected from the simulation with g bb = 2.5 × 10 −3 and g 4b = 0.4 × 10 −3 . Also shown are the efficiencies to select qqbb (q = udsc) and bbbb signal events, predicted from the Monte Carlo simulation.

Number of events
Monte Carlo studies show that events from the process g → bb are preferentially selected in class "3+1". The corresponding selection efficiencies are given in Table 1. For a signal rate g bb = 2.5×10 −3 , in class "3+1" 1.47% of the events are from g → bb. In class "2+2" only 0.99% of the events are from g → bb. The event classification procedure selects 235 201 data events in class "2+2" and 208 133 in class "3+1". The Monte Carlo simulation predicts 240 774 and 202 560 events, respectively. Implications on the final result for g bb and g 4b from this statistically significant difference in the number of observed to expected events in those two event classes are considered in the discussion of the experimental uncertainties.
The two g → bb candidate jets selected in the first step are checked for secondary vertices. If both of these jets have a reconstructed secondary vertex with l/σ > 3, the event is selected. This cut is studied in Fig. 3a, where the decay length significance (l/σ) 2 of the second g → bb candidate jet is shown. The two jets are ordered such that (l/σ) 1 > (l/σ) 2 . If the variable (l/σ) 2 is larger than 3, the event is selected, otherwise it is rejected. Note that (l/σ) 2 > 3 implies (l/σ) 1 > 3. Because only the smaller decay length significance of the two vertices enters this distribution, there are many entries at negative (l/σ) 2 values. The cut (l/σ) 2 > 3 reduces the background from events with light flavors and it also reduces the fraction of events from the g → cc process.
In the next step, the output of a neural network trained to recognize vertices from b hadron decays is calculated for the reconstructed vertices in the selected events. The two g → bb candidates are ordered by their corresponding neural network output NN such that NN 1 > NN 2 . The sum NN 1 + NN 2 shows good sensitivity to the signal process g → bb for both event classes "2+2" and "3+1". It is depicted in figures 3b and 3c. To enrich the signal, cuts NN 1 + NN 2 > 1.7 in class "2+2" and NN 1 + NN 2 > 1.1 in class "3+1" are made. The cut values are chosen such as to obtain purities of about 40% for each class of events. Some discrepancies between data and the prediction are observed in Fig. 3 and in table 1, mainly at low values of l/σ and low values of NN 1 + NN 2 . These can be explained by uncertainties in the knowledge of the detector resolution, studied in Section 4.3.
Finally, the other two jets that were not considered to be gluon splitting candidates are checked for secondary vertices to gain some separation power for the bbbb signal events. These jets are referred to as "primary quark jet candidates". If both of these jets have a vertex, only the jet with the larger decay length significance is considered. Fig. 4 shows the distribution of the decay length significance of this reconstructed vertex, (l/σ) p , for the events selected from class "2+2" and "3+1". The variable (l/σ) p shows some separation between the bbbb events and the qqbb (q = udsc) events. A cut at (l/σ) p = 2 is made to select event samples enriched or depleted with bbbb events. The following four event samples are defined Event sample A: Events from class "2+2" with NN 1 + NN 2 > 1.7 and (l/σ) p > 2, enriched in bbbb.
Event sample B: Events from class "2+2" with NN 1 +NN 2 > 1.7 not selected in sample A, depleted in bbbb.
Event sample D: Events from class "3+1" with NN 1 +NN 2 > 0.7 not selected in sample C, depleted in bbbb.

The dedicated bbbb event selection
Starting again with the entire sample of four-jet events, selected with y 34 > 0.006, a dedicated selection of bbbb events is set up, independent of the cuts presented in the previous section. To find events with four b hadrons, events are selected where at least three decay vertices with a decay length significance l/σ > 3 are found. Fig. 5a illustrates this selection of three significant vertices, showing the third largest decay length significance, denoted (l/σ) 3 . The bbbb signal is enhanced for large values of this variable. The cut (l/σ) 3 > 3 suppresses light flavors, most of the g → cc events, and also the qqbb (q = udsc) events. Fig. 5b shows the third largest neural network output NN 3 for all selected vertices, where the events already selected in sample A-D are excluded to avoid double counting. The background dominates the region of low NN 3 , while the bbbb signal extends to high values of NN 3 . A cut NN 3 > 0.7 is chosen to select the final bbbb candidates from distribution, denoted as sample E.

Calculation of g bb and g 4b
After applying all cuts, 221 events remain in the event samples A-E, where the simulation with g bb = 2.5 × 10 −3 and g 4b = 0.4 × 10 −3 predicts 217 ± 7 events. The efficiencies for selecting signal events, qqbb (q = udsc) or bbbb events, are denoted ǫ qqbb and ǫ bbbb . The total efficiencies as obtained from the simulation are ǫ qqbb = (0.64 ± 0.02)% and ǫ bbbb = (1.64 ± 0.07)%. The uncertainties quoted on the number of expected events and the selection efficiencies are due to the limited number of Monte Carlo events. Table 1 summarizes the number of selected events for the cuts applied, the number of signal and background events expected from the simulation, and the efficiencies to select signal or background reactions.
For the 221 signal events the angle α 12−34 between two jet-jet planes is studied to distinguish signal events from the background, dominated by events with two quarks and two gluons. The first plane is spanned by the two jets that are joined into one jet by the jet-algorithm at the transition from four to three jets. The other plane is formed by the other two jets. The definition of this angle α 12−34 is similar to the angular correlation proposed in [10] to measure the QCD color factors. In Fig. 6, the α 12−34 distribution is shown to be consistent with the theoretical prediction. The signal appears preferentially at high values of α 12−34 while the background has a flatter distribution.
A maximum likelihood fit of the 26 bins in α 12−34 shown in Fig. 6 is performed to extract g bb and g 4b , assuming Poisson distributions calculated from the signal and background efficiencies for each bin. The likelihood function is given by where d i is the number of observed events in bin i and µ i is the corresponding number of predicted events (assumed to be the mean of the Poisson distributions). The latter is given by where N 4−jet is the total number of selected four-jet events used for this analysis. The rate of gluon splitting to cc per hadronic Z 0 decay is denoted g cc . The two rates g bb and g 4b are taken as the fit parameters. The probabilities to select the signal process for primary light quarks and the bbbb events in bin i are denoted ǫ i qqbb and ǫ i bbbb , while ǫ i qqcc and ǫ i qq are the probabilities to select the g → cc events or other background events in bin i. Finally ǫ 4−jet qqbb , ǫ 4−jet bbbb , ǫ 4−jet qqcc , ǫ 4−jet qq are the efficiencies for events in the signal and the background channels to pass the four-jet selection. The likelihood analysis leads to the result g bb = (3.07 ± 0.53) × 10 −3 , with a correlation coefficient between g bb and g 4b of +0.007. The errors are statistical only.

Evaluation of systematic uncertainties
Various sources of systematic uncertainties are considered, as summarized in Table 2 and discussed in the following. Additional cross-checks are presented in Section 5.

Model dependence
Though the parton shower approach as implemented in Jetset 7.4 is expected to describe the gluon splitting process quite well [1], it is desirable to look at various alternative models and at exact calculations. In Pythia 6.130 [17], using the option MSTJ(42)=3, the calculation of the opening angle in the gluon splitting process has been modified to take mass effects into account, as compared to Jetset 7.4. The Pythia 6.130 prediction for the g → bb process is used to evaluate the main  Table 2: Summary of the systematic uncertainties on the rates g bb and g 4b and the corresponding correlations results of this analysis. The Herwig 5.91 [25] Monte Carlo generator 2 provides an alternative model for the parton shower. Another approach is given by the Color Dipole Model, using the event generator Ariadne 4.08 [26].
The program Wphact [9] implements e + e − annihilation to four fermions, where the masses of b quarks are taken into account for the calculation of the matrix elements. It can be used to study bbbb and bbqq final states. Version 1.3 of Wphact 3 [27] allows for the possibility to switch off Feynman graphs where the b quarks couple to the Z 0 or γ, enabling studies of the gluon splitting process at tree level. For the case of primary light quarks, calculations are available including next-to-leading order logarithmic terms [1,2]. The shapes of the Wphact 1.3 predictions are in agreement with these calculations in the regions where such a comparison is possible. Fig. 7 shows differential distributions of kinematic variables of the gluon that splits to bb, calculated with Pythia 6.130,Jetset 7.4, Herwig 5.91, Ariadne 4.08 and Wphact 1.3. The "gluon" variables have been calculated from the bb quark pair at the end of the shower. For the case of bbbb events the bb quark pair with the lowest invariant mass is chosen. The hardest energy spectrum is predicted by Jetset 7.4, while Ariadne 4.08 leads to the softest energy spectrum, when comparing the five models. For the gluon virtuality the Ariadne 4.08 prediction leads to the hardest distribution, while Herwig 5.91 predicts the softest spectrum. For the decay angle in the gluon rest frame the models differ most significantly at high | cos θ ⋆ |. The extreme cases are covered by Jetset 7.4 with a low number of events in this region and by Ariadne 4.08 showing an enhancement of events with increasing | cos θ ⋆ |. Note that the Pythia 6.130 prediction is always well between those extreme cases, this is why we decided to use it for the main analysis results.
To study hadronisation and detector effects, we use two sets of events generated either with Jetset 7.4 or with Pythia 6.130, and having a detailed hadronisation and detector simulation. The Jetset 7.4 sample is only used for systematic checks. The differential efficiencies to select the signal events, qqbb (q = udsc) and bbbb in this analysis, are shown in figure 8. They are evaluated with the Pythia 6.130 event sample. The efficiencies are low for small gluon energies E and gluon virtualities m, as it is difficult to detect the two b hadrons in two separate jets in this region. The same is true for decay angles | cos θ ⋆ | close to one, because in this case the two b hadrons have a small transverse momentum relative to the gluon flight direction.
To evaluate model-dependent effects, the complete analysis chain is performed using either the Pythia 6.130 or the Jetset 7.4 event sets. In addition, the Pythia 6.130 events are weighted such as to reproduce the gluon distributions from Fig. 7 predicted by Herwig 5.91, Ariadne 4.08 and Wphact 1.3, and the analysis is repeated. Table 3 summarizes the individual results of these studies. To study possible deficiencies in the reweighting procedure, it is repeated based on the Jetset 7.4 event sample rather than on the Pythia 6.130 event sample. This results in slightly lower values of g bb and larger values of g 4b , as compared to the Herwig 5.91, Ariadne 4.08 and Wphact 1.3 values shown in table 3. The most significant differences to the standard reweighting procedure are observed for Herwig 5.91, namely g bb = 2.88 × 10 −3 and g 4b = 0.53 × 10 −3 .   Finally, the largest differences on g bb and g 4b to the central value are taken as systematic uncertainty due to limited knowledge of the gluon splitting mechanism. For g bb this is the Pythia 6.130-Jetset 7.4 difference, while for g 4b this is the Pythia 6.130-Herwig 5.91 difference, where the Herwig 5.91 result is calculated with the help of reweighted Jetset 7.4 events.

Dependence on the b quark mass
The influence of the b quark mass assumed for the Monte Carlo simulation has been studied by changing the b quark mass in the Pythia 6.130 program 4 and investigating the properties of the b hadrons produced in the fragmentation process. The central value for the b quark mass in Pythia 6.130 used throughout this analysis is 5 GeV/c 2 . The shapes of the distributions of various kinematic variables are used to reweight the g → bb Monte Carlo events corresponding to a b quark mass of 4.5 GeV/c 2 and 5.25 GeV/c 2 . This variation covers the uncertainty of the b quark pole mass [28] up to the B meson mass. Note that the way the quark mass parameter is used in Pythia 6.130 corresponds rather to a constituent quark mass than to a pole mass definition of the quark mass [29], which justifies the choice of 5 GeV/c 2 for the central value.
The quantities used for the reweighting process are the momenta of the two b hadrons produced in the fragmentation process and their invariant mass. These variables are chosen because the efficiency to identify b hadrons strongly depends on their momentum. The efficiency to resolve the two b hadrons in different jets depends on their invariant mass. The larger deviations from the standard results are found for a b quark mass of 4.5 GeV/c 2 . They are taken into account as systematic uncertainties.

Experimental sources of systematic uncertainties
The modeling of the OPAL detector is important because the analysis depends on an accurate understanding of the decay vertex reconstruction. The b tagging efficiencies are mainly sensitive to the parameters modeling the production and decay of b hadrons. This will be discussed in Section 4.4. The light quark tagging efficiencies in the Monte Carlo simulation are mainly sensitive to details in the modeling of the tracking system. Studies are done by increasing the difference of the reconstructed track parameters with respect to the true track parameters in the Monte Carlo simulation by 10% to cover uncertainties in the knowledge of the flavor tagging variables [24]. This smearing is applied separately for parameters defined in the (r, φ) and the (r, z) plane. In (r, φ) a simultaneous smearing of the distance of closest approach, d 0 , and the azimuthal angle φ 0 is performed. In (r, z) smearing is done simultaneously for the z coordinate of the point of closest approach in (r, φ) and the polar angle θ. Finally the gluon splitting analysis is repeated using the modified Monte Carlo sets. The uncertainties from the smearing in (r, φ) and (r, z) are added quadratically.
The normalization to the number of four-jet events revealed differences in the population of the event classes "2+2" and "3+1"when comparing the data with the Monte Carlo prediction, as discussed in Section 3.1. This is addressed by repeating the likelihood fit, using the numbers of events and selection efficiencies of the "2+2" and "3+1" classes, rather than the number of four-jet events with their corresponding selection efficiencies. The variation found is assigned as systematic uncertainty due to the event classification.

Uncertainties from heavy flavor physics
The rate of primary bb production in the four-jet sample R 4−jet b can be measured using a double-tag technique. Agreement with the Monte Carlo simulation within 2.4% is found, where the statistical uncertainties are in the order of 1.5%. The rate R 4−jet b in the Monte Carlo sample has been changed by 2.4% and the differences in g bb and g 4b to our standard results are taken as systematic uncertainties.
Varying the rate of primary cc production R c within its uncertainties only has a negligible effect on the results.
The rate of secondary cc quark pairs from gluons, g cc , has been measured by OPAL to be g cc = (3.20 ± 0.43) × 10 −2 [20]. This measurement is based on lepton identification and the reconstruction of D ⋆ mesons in three-jet events. It can be considered as statistically independent of the g bb measurement presented here, which is based on lifetime-tags in four-jet events. The rate g cc thus is varied within its uncertainty to obtain systematic uncertainties on g bb and g 4b . The explicit dependence of the results of this analysis on g cc is given by The fragmentation functions for charm and bottom quarks are varied to reflect the uncertainties in the knowledge of the average scaled energies for D and B mesons, x E c = 0.484 ± 0.008 and x E b = 0.702 ± 0.008 [31]. This is done by varying the parameter ǫ in the parameterization of the fragmentation function suggested by Peterson et. al. [32]. In addition the sensitivity to the shape of the fragmentation function is checked, using the models suggested by Collins and Spiller [33] and Kartvelishvili [34]. The parameters of these models are chosen to reproduce the measured values of the mean scaled energies of the D and B mesons in the Monte Carlo simulation, and the larger deviations in g bb and g 4b are used as systematic uncertainties from this source. The uncertainties from the knowledge of the mean scaled energies and the shape of the distributions are added in quadrature.
The charged decay multiplicities of the D mesons are varied within the errors given in [35] around the Jetset 7.4 prediction. Particles from decay-chains of b hadrons are excluded from this variation. The mean charged decay multiplicity of weakly decaying b hadrons is varied within n B = 4.995± 0.062 [31]. This multiplicity includes secondary decays of charm hadrons produced in the decay chains. Variations of the neutral decay multiplicities have not been studied, but are expected to be small.
The production rates of charm and bottom hadrons in the Monte Carlo simulation are varied within the ranges given by the LEP Electroweak Working group [31] and the Particle Data Group [28], respectively. This variation includes particles from primary bb and cc production as well as those from gluon splitting to heavy quarks.
The lifetimes of charm and bottom hadrons are varied around their central value according to the numbers given by the Particle Data Group [28].
The systematic errors on g bb and g 4b have a correlation +0.78. This large positive correlation can be understood from the fact that most of the systematic variations influence the selection efficiencies for both the g → bb and the bbbb signal in the same direction.

Additional cross-checks
The four-jet selection is defined using a cut y 34 > y min 34 . Fig. 9a shows the number of events in the data divided by the predicted number of events as a function of this cut. The disagreement in the rate of four-jet events, as discussed in Section 2.4, is clearly visible. The fit results for g bb and g 4b as a function of this cut are studied in Fig. 9b and 9c. The fit results are stable under the variation of y min 34 within the independent statistical error. Note that the variation of y min 34 from 0 to 0.014 corresponds to a variation in the number of candidate events from 1003 to 93 and a variation of the estimated g → bb signal purity from 10% to 50%. These checks show that the treatment of the normalization to the four-jet rate in the likelihood fit is correct.
To check the fitting procedure a likelihood fit is applied to the normalized signal and background shapes in α 12−34 , calculating directly the fractions of qqbb (q = udsc) and bbbb events in the final event selection from these shapes. Using the number of selected events, the total number of hadronic Z 0 decays and the signal efficiencies, g bb and g 4b can be calculated from these event fractions. This procedure avoids the need to know the absolute selection efficiency for the background. Alternatively, the likelihood fit is repeated using only the individual signal and background selection efficiencies of the event samples A-E. This analysis is thus independent of the variable α 12−34 . Both results are compatible with the default method.
The k ⊥ algorithm is used in this analysis not only to define the jets, but also to define the event topology. Therefore it is important to test the analysis, using a different algorithm. The JADE E0 [30] jet clustering scheme is used as an alternative method. Four-jet events are selected with a cut y JADE 34 > 0.015, resulting in a similar number of four-jet events as compared to the default analysis which is based on the k ⊥ algorithm. To classify the events and define the gluon splitting candidates, as described in section 3.1, and to define the angle α 12−34 , as described in Section 3.3, the new definition of the jet resolution parameter in the JADE E0 scheme is used. All other analysis cuts are unchanged with respect to the standard analysis. This results in 207 selected events. The result of the likelihood fit using the JADE E0 algorithm is g JADE bb = (3.11 ± 0.59) × 10 −3 and g JADE 4b = (0.45 ± 0.17) × 10 −3 (statistical uncertainties only). These results are compatible with the results obtained with the k ⊥ algorithm within the statistical errors, considering the small fraction of events common to both selections.

Summary
A measurement of the inclusive rate of gluon splitting to bb per hadronic Z 0 decay has been performed using data taken by OPAL. The result is g bb = (3.07 ± 0.53 ± 0.97) × 10 −3 .
The correlation of g bb and g 4b is +0.007 from statistical uncertainties and +0.78 from systematic uncertainties. The ratio g 4b /g bb is thus measured to be g 4b g bb = 0.116 ± 0.060 ± 0.065. This is in agreement with the simplified expectation g 4b /g bb ≈ R b and with a more detailed calculation performed by DELPHI [8], predicting g bb /g 4b = 0.1833 ± 0.0003. The results for g bb and g 4b are compatible with previous measurements [6,7,8] as well as theoretical predictions [1,4].