Measurement of charm and beauty production in deep inelastic ep scattering from decays into muons at HERA

The production of charm and beauty quarks in ep interactions has been measured with the ZEUS detector at HERA for squared four-momentum exchange Q^2>20 GeV^2, using an integrated luminosity of 126 pb^{-1}. Charm and beauty quarks were identified through their decays into muons. Differential cross sections were measured for muon transverse momenta p_T^{\mu}>1.5 GeV and pseudorapidities -1.6<\eta^{\mu}<2.3, as a function of p_T^{\mu}, \eta^{\mu}, Q^2 and Bjorken x. The charm and beauty contributions to the proton structure function F_2 were also extracted. The results agree with previous measurements based on independent techniques and are well described by QCD predictions.


Introduction
The measurement of charm and beauty production in deep inelastic scattering (DIS) provides a stringent test of quantum chromodynamics (QCD) since the large quark masses provide hard scales that make perturbative calculations applicable. At leading order, heavy quarks (HQs) are produced in DIS via boson-gluon fusion (BGF) (γ * g → qq). A precise measurement of HQ production in DIS therefore provides a direct constraint on the gluon parton density function (PDF) of the proton.
Charm production in DIS at HERA has been measured previously using reconstructed charmed mesons [1,2] or inclusively by exploiting the long lifetime of charmed hadrons [3]. Beauty production in DIS has been studied in events with muons and jets [4,5] and from lifetime information [3]. The existing data are generally in good agreement with nextto-leading-order (NLO) QCD predictions. The largest differences were observed in the muon analyses [4,5] where the measured beauty cross section was about two standard deviations above the theoretical expectation.
In this paper, a simultaneous measurement of beauty and charm production using semileptonic (SL) decays into muons is presented. The fractions of muons originating from charm, beauty and light flavours (LF) were extracted by exploiting three discriminating variables: the muon impact parameter, the muon momentum component transverse to the associated jet axis and the missing transverse momentum, which is sensitive to the neutrino from SL decays.
The analysis focused on data with large squared four-momentum exchange at the electron vertex, Q 2 , where charm measurements based on muons are competitive with those based on identified charmed mesons.
The cross sections for muons from charm and beauty decays were measured for Q 2 > 20 GeV 2 , muon transverse momenta p µ T > 1.5 GeV and pseudorapidities 1 −1.6 < η µ < 2.3 as a function of p µ T , η µ , Q 2 , and of the Bjorken scaling variable x [6] and compared to QCD predictions. The muon cross sections, measured in bins of x and Q 2 , were used to extract the heavy quark contributions to the proton structure function F 2 which were compared to previous results and to QCD predictions. The data used in this analysis were collected with the ZEUS detector in the 2005 running period during which HERA collided electrons with energy E e = 27.5 GeV with protons with E p = 920 GeV corresponding to a centre-of-mass energy √ s = 318 GeV. The corresponding integrated luminosity was L = 126.0 ± 3.3 pb −1 .
1 The ZEUS coordinate system is a right-handed Cartesian system, with the Z axis pointing in the proton beam direction, referred to as the "forward direction", and the X axis pointing towards the centre of HERA. The pseudorapidity is defined as η = − ln tan θ 2 , where θ is the polar angle. 1

Theoretical predictions
Heavy quark production in DIS has been calculated at next-to-leading order (O(α 2 s )) in the so-called fixed flavour number scheme (FFNS) in which only light flavours are present in the proton and heavy quarks are produced in the interaction [7]. The results of this analysis have been compared to NLO calculations performed with the Hvqdis program [8,9]. The renormalisation and factorisation scales were set to µ 2 R = µ 2 F = Q 2 + 4m 2 q and the quark masses to m c = 1.5 GeV and m b = 4.75 GeV. The PDFs were obtained by repeating the ZEUS-S [10] PDF fit in the FFNS with quark masses set to the same values as in the Hvqdis calculation.
To calculate muon observables, the partonic results were interfaced to a model of HQ fragmentation into weakly decaying heavy hadrons and of the decay of heavy hadrons into muons. The hadron momentum was obtained by scaling the quark momentum according to the fragmentation function of Peterson et al. [11] with the parameter ǫ c = 0.055 for charm and ǫ b = 0.0035 for beauty. This choice of ǫ c corresponds to ǫ c = 0.035 for D * mesons [12] since kinematic considerations [13] and direct measurements [14] show that, on average, the momentum of the weakly decaying hadrons is ≈ 5% lower than that of D * mesons.
The semileptonic decay spectrum for charm was taken from a recent CLEO measurement [15]. The decay spectrum for beauty hadrons was taken from the Pythia [16] Monte Carlo (MC), mixing direct SL decays and cascade decays through charm according to the measured branching ratios [17]. It was checked that the MC described BELLE and BABAR data [18] well. The branching ratios were set to B(c → µ) = 0.096 ± 0.004 and B(b → µ) = 0.209 ± 0.004 [17].
The uncertainty on the theoretical predictions was evaluated by independently varying µ R and µ F by a factor two; by varying the HQ masses simultaneously to (m c , m b ) = (1.3, 4.5), (1.7, 5.0) GeV in the calculation and in the PDF fit; by varying the proton PDFs by their experimental uncertainty and by varying the fragmentation parameters within 0.04 < ǫ c < 0.12 (corresponding to 0.025 < ǫ c < 0.085 for D * mesons [19]) and 0.0015 < ǫ b < 0.0055. As a further check, the fragmentation was performed by scaling the sum of the energy and the momentum parallel to the HQ direction, E+p || , rather than the HQ momentum. The total theoretical uncertainty was obtained by adding in quadrature the effects of each variation. In the beauty case, the total uncertainty is dominated by the variation of µ R and of the mass while for charm the variation of ǫ c also gives a large contribution.
The calculations of F cc 2 and F bb 2 in the FFNS were performed using Hvqdis and cross checked with the QCD evolution code [20] used in the ZEUS PDF fit.

Monte Carlo samples
Charm and beauty MC samples were generated using Rapgap 3.00 [21] to simulate the leading order BGF process. Parton shower techniques were used to simulate higher order QCD effects. Higher order QED effects were included through Heracles 4. 6 [22]. The CTEQ5L [23] PDFs were used and the HQ masses were set to m c = 1.5 GeV and m b = 4.75 GeV.
Light flavour MC events were extracted from an inclusive DIS sample generated with DjangoH 1.3 [24] which is interfaced to Lepto 6.5 [25] to simulate the hadronic final state with the matrix element plus parton shower (meps) model and to Heracles 4.6 to include electroweak radiative corrections. The CTEQ5D [23] parton density was used.
Inelastic J/ψ production was simulated with Cascade [26] since that model generally describes the DIS data of a previous publication [27].
The above samples corresponded to at least five times the luminosity of the data. A smaller light quark sample was generated with Rapgap and mixed with the heavy quark Rapgap samples for the study of the inclusive DIS control sample (Section 6).
Fragmentation and particle decays were simulated using the Jetset/Pythia model [28,16]. The lepton energy spectrum from charm decays was reweighted to agree with CLEO data [15]. The MC events were passed through a full simulation of the ZEUS detector based on Geant 3.21 [29]. They were then subjected to the same trigger criteria and reconstructed with the same programs as used for the data.

Experimental set-up
A detailed description of the ZEUS detector can be found elsewhere [30]. A brief outline of the components that were most relevant for this analysis is given below.
Charged particles were tracked in the silicon microvertex detector (MVD) [31] and in the central tracking detector (CTD) [32], which operated in a magnetic field of 1.43 T provided by a thin superconducting solenoid. The MVD consisted of a barrel (BMVD) and a forward (FMVD) section with three cylindrical layers and four vertical planes of singlesided silicon detectors, respectively. The CTD consisted of 72 cylindrical drift chamber layers, organised in 9 superlayers covering the polar-angle region 15 • < θ < 164 • . After alignment, the single-hit resolution of the BMVD was 25 µm and the impact parameter resolution of the CTD-BMVD system for high-momentum tracks was ≈ 100 µm.
The high-resolution uranium-scintillator calorimeter (CAL) [33] consisted of three parts: the forward (FCAL), the barrel (BCAL) and the rear (RCAL) calorimeters. Each part was subdivided transversely into towers and longitudinally into one electromagnetic section and either one (in RCAL) or two (in BCAL and FCAL) hadronic sections. Under testbeam conditions, the CAL single-particle relative energy resolutions were σ(E)/E = 0.18/ √ E for leptons and σ(E)/E = 0.35/ √ E for hadrons, with E in GeV. The energy of electrons hitting the RCAL was corrected for the presence of dead material using the rear presampler detector (PRES) [34] and the small angle rear tracking detector (SRTD) [35].
The muon system consisted of rear/barrel (R/BMUON) [36] and forward (FMUON) [30] tracking detectors. The B/RMUON consisted of limited-streamer (LS) tube chambers placed behind the BCAL (RCAL), inside and outside a magnetised iron yoke surrounding the CAL. The barrel and rear muon chambers cover polar angles from 34 o to 135 o and from 135 o to 171 o , respectively. The FMUON consisted of six trigger planes of LS tubes and four planes of drift chambers covering the angular region from 5 o to 32 o . The muon system exploited the magnetic field of the iron yoke and, in the forward direction, of two iron toroids magnetised to ≈ 1.6 T to provide an independent measurement of the muon momentum.
The luminosity was measured using the Bethe-Heitler reaction ep → eγp with the luminosity detector which consisted of two independent systems, a photon calorimeter [37] and a magnetic spectrometer [38].

Event reconstruction and selection
A three-level trigger was used to select events online [30,39]. DIS events were selected by requiring a scattered electron in the CAL.
A scattered electron with energy E ′ e > 8 GeV was required offline. The primary vertex had to be within ±30 cm in Z from the nominal interaction point.
Muons were reconstructed by matching a CTD+MVD track to a track segment in the inner or outer B/RMUON chambers or to an FMUON track crossing at least four FMUON planes. This B/RMUON selection was looser than in some previous analyses, which required the muons to reach the external chambers [4,40], allowing a lower threshold for the muon transverse momentum.
The central track associated to a B/RMUON candidate was required to pass at least three CTD superlayers and to have at least four hits in the MVD to allow a good impact parameter measurement. The tracks associated to FMUON candidates were required to pass at least one CTD superlayer, corresponding to at least four degrees of freedom in the track fit.
Muons were accepted in the kinematic region defined by The hadronic system (including the muon) was reconstructed from energy flow objects (EFOs) [41] that combine the information from calorimetry and tracking, corrected for energy loss in the dead material. The EFOs were corrected using the measured momenta of identified muons [40,42].
To select a clean DIS sample, the following cuts on global variables were applied: [43], and θ e is the electron polar angle. These cuts restricted the accessible inelasticity y = Q 2 /(xs) and Q 2 to 0.01 < y < 0.7 and Q 2 > 20 GeV 2 . The DIS variables x and Q 2 were reconstructed using the Σ estimators Q 2 Σ and x Σ = Q 2 Σ /(sy Σ ) [43]. To remove background events with isolated muons (γγ → µ + µ − , J/ψ and Υ decays) and residual cosmic muons, an anti-isolation cut was applied by requiring that the hadronic energy in a cone of radius 1 in the η − φ plane around the muon candidate, excluding the muon itself, was E iso > 0.5 GeV. From MC studies this cut was 98% (90%) efficient for charm (beauty).
Jets were reconstructed from EFOs using the k T algorithm [44] in the longitudinally invariant mode [45]. About 96% of the muon candidates were associated to a jet with transverse momentum (including the muon) p jet T > 2.5 GeV and kept for further analysis. After the above selection, the final sample contained 11126 muons. A subsample of 35 events with more than one muon was found, 28 of which consisted of µ + µ − pairs. A J/ψ signal of 9 events was observed in the µ + µ − invariant mass distribution. The total contamination from J/ψ production was estimated with the Cascade MC, normalised to the observed J/ψ signal. It was found to be (0.9 ± 0.3)% and was neglected in the analysis.

Extraction of the charm and beauty fractions
The sample of selected muon candidates contained signal muons from charm and beauty decays and background from in-flight π ± and K ± decays and from the punch through of hadronic jets in the muon chambers. Candidates from in-flight decays and punch through, which are subsequently denoted as "false muons", were present both in the LF events and in events containing HQs.
The fractions of muons originating from charm, beauty or LF events were determined from a simultaneous fit of three discriminating variables sensitive to different aspects of HQ decays: • p rel T , the muon momentum component transverse to the axis of the associated jet, Due to the large b mass, muons from beauty hadron decays have a harder p rel T spectrum than those from charm or light quarks; • δ, the distance of closest approach of the muon track to the centre of the interaction region (beam spot) in the X, Y plane. A positive sign was assigned to δ if the muon track crosses the axis of the associated jet in the jet hemisphere, negative otherwise. The beam spot position was obtained by fitting the reconstructed primary vertex distribution for every 2000 ep events. The size of the interaction region was 80×20 µm 2 in X×Y . Muons from decays of long-lived heavy quarks tend to have positive δ while tracks originating from the primary interaction have a symmetric δ distribution around zero, corresponding to the experimental resolution.
• p miss||µ T , the missing transverse momentum parallel to the muon direction. The missing transverse momentum vector was calculated using the electron and the EFOs. The p miss||µ T distribution has a positive tail of events containing semileptonic HQ decays due to the presence of the neutrino.
A control sample of inclusive DIS data, selected similarly to the muon sample but without any muon requirement, was used to test the quality of the simulation of these variables. The control sample is dominated by LF events, containing, according to MC, about 18% (1%) of c (b) events. The p rel T distribution of inclusive tracks in the control sample was reasonably well reproduced by both the DjangoH and the Rapgap inclusive DIS samples. The small differences (at most 10% at p rel T > 2 GeV) were corrected for by applying a bin-by-bin correction to the p rel T distribution of the LF and charm MC samples similarly to a previous publication [40]. The quality of the MC description of p miss||µ T was also evaluated in the control sample by studying a similar p T -balance variable: the missing transverse momentum parallel to the electron p miss||e T . The best description of the p miss||e T distribution of the inclusive DIS sample was obtained by shifting the hadronic transverse momentum by (0.1 ± 0.1) GeV in the MC and by increasing the hadronic transverse momentum resolution by (5 ± 5)% in the case of Rapgap and by (0 ± 5)% in the case of DjangoH. The resolution on δ was studied using tracks in the inclusive DIS sample. Since it was underestimated in the MC by ≈ 15%, a p T -dependent smearing [46] was applied to the MC, similarly to what was done in a previous publication [47]. 6 The fractions of b, c and LF events were obtained by fitting a combination of MC distributions to the measured three-dimensional distribution of the discriminating variables [48]. The fit range was |δ| < 0.1 cm, p rel T < 2.5 GeV and |p miss||µ T | < 10 GeV. A precise measurement of δ was only possible inside the region covered by the BMVD. Hence for events with muons reconstructed in the FMUON (4% of the total) only p rel T and p miss||µ T were used in the fit. A Poisson likelihood fit was used, taking into account the limited MC statistics.
The global charm and beauty fractions resulting from the fit were with a correlation coefficient ρ cb = −0.43. Figure 1(a-c) shows the distributions of the three discriminating variables compared to the MC distributions with the normalisation corresponding to the fit. While δ and p miss||µ T provide discrimination between LF and HQs, p rel T discriminates between beauty and the other components. Figure 1 Σ and x Σ for the data and for the MC samples normalised according to the fit are shown in Fig. 2. The overall agreement is satisfactory.
The cross sections were calculated using where f q is the HQ fraction from the fit, N is the number of reconstructed muons, A q is the acceptance, C r is the QED radiative correction, and q = c, b. Differential cross sections were measured by repeating the fit in bins of the reconstructed variable V as dσ/dV = σ q i /∆V i , where σ q i is the cross section in the bin and ∆V i is the bin width. The acceptance A q was evaluated from the MC simulation as the number of reconstructed muons divided by the number of true muons from decays of the quark q. This definition takes into account the charm and beauty events in which a "false muon" is reconstructed rather than a signal muon from a HQ decay. The acceptance included the efficiency of muon reconstruction (which in turn includes the efficiency of the muon chambers and of the matching with central tracking) that was evaluated from an independent exclusive dimuon sample as explained in previous publications [40,42]. The muon reconstruction efficiency was around 50% for central muons with p µ T > 2 GeV. The acceptance for c (b) ranged from 23% (16%) at 1.5 < p µ T < 2.5 GeV to ≈ 35% (25%) at p µ T > 2.5 GeV. The difference in acceptance between c and b was mainly due to the different contribution from "false muons" which was ≈ 25% for c and ≈ 3% for b.
According to the MC simulation, the probability to find a "false muon" in a DIS event (before any muon selection) is P false ≈ 0.1%, almost independently from the event being c, b or LF. The ability of the MC to reproduce P false was studied by comparing the number of LF events in the data, as obtained from the fit, to the absolute prediction by DjangoH. The data/MC ratio, P data false /P MC false , was estimated as 0.80 ± 0.20 in the RMUON, 1.10 ± 0.20 in the BMUON, and 1.05 ± 0.40 in the FMUON, in agreement with previous studies [49].
The cross sections were corrected to the QED Born level, calculated using a running coupling costant α em , such that they can be compared directly to the QCD predictions from the Hvqdis program (Section 2). The radiative corrections were obtained as C r = σ Born /σ rad , where σ Born is the Rapgap cross section with the QED corrections turned off but keeping α em running and σ rad is the Rapgap cross section with the full QED corrections, as in the standard MC samples. The corrections were typically C r ≈ 1.05 and at maximum 1.10 in the highest Q 2 bin.

Systematic uncertainties
The following systematic uncertainties were considered (the effects on the total visible cross section for c and for b is given in parentheses): 8. resolution on δ: the smearing applied to the MC was varied by ±25% as allowed by the control sample ( −3 +2 , +11 −9 )%; 9. p rel T shape of LF and charm: it was evaluated by varying the p rel T correction by ±50% (∓1.5, +8 −5 )%; 10. hadronic energy flow near the muon: it was evaluated by varying the cut on E iso by +0.50 −0.25 GeV (0, −1 0 )%; 11. jet description: the cut on p jet T was varied by ±0.5 GeV (±2.5, −3.5 +2.5 )%; 12. charm SL decay spectrum: the reweighting to the CLEO model was varied by ±50%, ( −4 +3 , +3 −2 )%; 13. MC model dependence: the Rapgap c and b samples were reweighted to reproduce the corresponding measured differential cross sections in Q 2 or in p µ T and the largest deviation from the nominal cross section was taken (+6, +20)%; 14. higher order effects: this uncertainty was evaluated by varying the HQ distribution before parton showering in Rapgap by the difference between NLO and leading order, as evaluated with Hvqdis ( +6 −10 , +2 −3 )%; 15. MVD efficiency: the efficiency of the cut on the number of MVD hits of (90 ± 3)% was varied by its uncertainty (∓3, ∓3)%; 16. CTD simulation: tracks were required to pass ≥ 4 superlayers in the B/RMUON region and to have ≥ 7 degress of freedom in the FMUON region (+1, 0)%; 17. integrated luminosity: measurement uncertainty (∓2.6, ∓2.6)%.

Cross sections
The visible cross sections for muons from charm and beauty decays in the kinematic region of Eq. The differential cross sections as a function of p µ T , η µ , Q 2 , and x are presented in Table 1 and compared in Fig. 3 to the NLO QCD predictions based on Hvqdis. The Rapgap MC predictions are also shown, normalised according to the result of the global fit. Charm and beauty cross sections are similar for p µ T > 3.5 GeV. The charm cross sections are in good agreement with the Hvqdis calculations. The tendency of the beauty cross section to lie above the central NLO prediction is concentrated at low p µ T and low Q 2 . The statistical significance of the difference between the data and the NLO predicitons is similar to that obtained for the total visible cross section since the uncertainties are dominated by correlated systematics.
Both NLO calculations and the Rapgap MC give in general a good description of the shape of the differential cross sections. The Q 2 distributions for beauty is somewhat steeper than predicited by Rapgap and Hvqdis.

2
The heavy quark contribution to the proton structure functions, F qq 2 , F qq L and the reduced cross sectionσ qq are defined in analogy with the inclusive case from the double differential cross section in x and Q 2 for the production of the quark q: where K = Y + (2πα 2 em )/(xQ 4 ) and Y + = 1 + (1 − y) 2 . The muon cross sections, σ q , measured in bins of x and Q 2 , were used to extract F qq 2 at a reference points in the x, Q 2 plane by: where F qq,th 2 (x, Q 2 ) and σ q,th were calculated at NLO in the FFNS using the Hvqdis program. The reference points were chosen close to the average x and Q 2 of the events within each bin. Charm produced in b decays was not included in the definition of σ c . This procedure contains several corrections: the extrapolation from the restricted muon kinematic range (p µ T > 1.5 GeV, −1.6 < η µ < 2.3) to the full muon phase space; the q → µ branching ratio; the correction for the longitudinal structure function F qq L and the correction from a bin-averaged cross section to a point value (bin centring).
The largest uncertainty is related to the extrapolation to the full muon phase space. The kinematic acceptance, A, defined as the fraction of muons from HQ decays that was generated in the restricted kinematic region is, on average, A = 13%(27%), for charm (beauty). According to Hvqdis, in the charm case, A becomes sizeable (A > 0.25 A ) when one of the two charm quarks in the event has p T > 3 GeV and its rapidity is in the range (−1.5 : 2.5), which corresponds to the phase space containing 88% of the cross section. In the beauty case, A is sizeable over the full HQ phase space.
The theoretical uncertainty in the extraction of F qq 2 was evaluated by varying the Hvqdis parameters as explained in Section 2 and by using a different PDF set (CTEQ5F). It is dominated by the fragmentation uncertainty. As a further check, F qq 2 was also evaluated taking A from Rapgap and found to be consistent within the quoted uncertainties.
The muon cross sections in bins of x and Q 2 are given in Table 2. The extracted F cc 2 and F bb 2 are presented in Tables 3 and 4 Tables 3 and 4 are the factor A and the correction for the longitudinal structure function C L =σ qq /F qq 2 as obtained from the NLO theory. The effects of the individual sources of systematic and theoretical uncertainty are given in Tables 5 and 6 in the Appendix 2 . Figure 4 also contains a comparison of F cc 2 with previous results based on the measurement of D * mesons from ZEUS [2] and to results from the H1 collaboration based on inclusive lifetime tagging (VTX) [3]. The previous results were corrected to the Q 2 values of the present analysis, using the NLO theory. The agreement of the different data sets, obtained with different charm tagging techniques, is good. At high Q 2 , the precision of present data is similar or better than for the previous results. The NLO QCD calculations are also shown. Figure 5 shows the extracted F bb 2 from this analysis and also a previous H1 result [3], corrected to the reference Q 2 values used in the present analysis. The two data sets are in good agreement. The precision of the present measurement is similar to that of the H1 data at high Q 2 . The QCD calculations are also shown.

and shown in Figs. 4 and 5. Also given in
The structure functions F cc 2 and F bb 2 are also presented in Figs. 6 and 7 as functions of Q 2 for fixed values of x, compared to previous results corrected to the same reference x used in the present analysis.

Summary
The production of charm and beauty quarks was measured in DIS using their decay into muons. Total and differential cross sections for muons from c and b decays were measured in the kinematic region Q 2 > 20 GeV 2 ; 0.01 < y < 0.7; p µ T > 1.5 GeV; −1.6 < η µ < 2.3 and compared to NLO QCD calculations. The agreement is good for charm. Beauty is about a factor two above the central QCD prediction although still compatible within statistical and systematic uncertainties. The heavy quark contribution to the proton structure function F 2 was also measured and found to agree well with other measurements based on independent techniques. For Q 2 ≥ 60 GeV 2 the present results are of comparable or higher precision than those previously existing.    Table 1: Muon differential cross sections for charm and beauty as a function of η µ , p µ T , Q 2 and x. The last column shows the statistical correlation coefficient between charm and beauty.