Study of Charm Fragmentation into D^{*\pm} Mesons in Deep-Inelastic Scattering at HERA

The process of charm quark fragmentation is studied using $D^{*\pm}$ meson production in deep-inelastic scattering as measured by the H1 detector at HERA. Two different regions of phase space are investigated defined by the presence or absence of a jet containing the $D^{*\pm}$ meson in the event. The parameters of fragmentation functions are extracted for QCD models based on leading order matrix elements and DGLAP or CCFM evolution of partons together with string fragmentation and particle decays. Additionally, they are determined for a next-to-leading order QCD calculation in the fixed flavour number scheme using the independent fragmentation of charm quarks to $D^{*\pm}$ mesons.


Introduction
The production of charm quarks is expected to be well described by perturbative Quantum Chromodynamics (pQCD) due to the presence of a hard scale provided by the charm mass. The evolution of an "off-shell" charm quark via gluon radiation until it is "on-shell" can be calculated in pQCD in fixed order of the strong coupling or by summing all orders in the leading-log approximation. The transition of an on-shell charm quark into a charmed hadron is, however, not calculable within the framework of pQCD and is thus usually described by phenomenological models. One of the major characteristics of this transition process is the longitudinal momentum fraction transferred from the quark to the hadron, the distribution of which is parametrised by a fragmentation function.
Several phenomenological models are available to describe the transition of a quark into hadrons, for example the independent fragmentation [1], the string [2], and the cluster model [3]. The fragmentation function is unambiguously defined in a given context of a phenomenological model together with a pQCD calculation. Universality is then only expected to hold within this context.
The fragmentation function is not a directly measurable quantity as the momentum of the heavy quark is experimentally not accessible. Also the momentum distribution of the heavy hadron can only be measured within a restricted phase space. The momentum spectrum is further affected by the fact that some heavy hadrons are not produced directly, but are the result of decays of still heavier excited states, whose contribution is not well known.
The production of charmed hadrons has been measured and parameters of fragmentation functions have been determined in e + e − annihilation experiments [4]. The H1 and ZEUS collaborations have published total cross sections for the production of various charmed hadrons in deep-inelastic ep scattering (DIS) [5] and in photoproduction [6]. These data show that the probabilities of charm quarks to fragment into various final state hadrons are consistent, within experimental uncertainties, for e + e − and ep collisions.
In this paper the transition of a charm quark into a D * ± meson in DIS is further investigated. The normalised differential cross sections are measured as a function of two observables with different sensitivity to gluon emissions. The momentum of the charm quark is approximated either by the momentum of the jet, which includes the D * ± meson, or by the sum of the momenta of particles belonging to a suitably defined hemisphere containing the D * ± meson. The measurements are performed for two different event samples. The DIS phase space and the kinematic requirements on the D * ± meson are the same for both samples. In the first sample, referred to as the "D * ± jet sample", the presence of a jet containing the D * ± meson and exceeding a minimal transverse momentum is required as a hard scale. In the second sample no such jet is allowed to be present. This sample, referred to as the "no D * ± jet sample", allows the investigation of charm fragmentation in a region close to the kinematic threshold of charm production.
The normalised differential cross sections are used to fit parameters of different fragmentation functions: for QCD models as implemented in the Monte Carlo (MC) programs RAP-GAP [7] and CASCADE [8], which use the Lund string model for fragmentation as implemented in PYTHIA [9], and for a next-to-leading order (NLO) QCD calculation as implemented in HVQDIS [10] with the addition of independent fragmentation of charm quarks to D * ± mesons. The paper is organised as follows. Section 2 gives a brief description of the H1 detector. It is followed by the details of the event selection, the D * ± meson signal extraction and the jet selection in section 3. The experimental fragmentation observables are defined in section 4. The QCD models and calculations used for data corrections and for the extraction of fragmentation functions are described in section 5. The data correction procedure and the determination of systematic uncertainties is explained in section 6. In section 7 the results of the measurements and of the fits of the fragmentation parameters are given.

H1 Detector
The data were collected with the H1 detector at HERA in the years 1999 and 2000. During this period HERA collided positrons of energy E e = 27.5 GeV with protons of energy E p = 920 GeV, corresponding to a centre-of-mass energy of √ s = 319 GeV. The data sample used for this analysis corresponds to an integrated luminosity of 47 pb −1 .
A right handed Cartesian coordinate system is used with the origin at the nominal primary ep interaction vertex. The direction of the proton beam defines the positive z-axis (forward direction). Transverse momenta are measured in the x-y plane. Polar (θ) and azimuthal (φ) angles are measured with respect to this reference system. The pseudorapidity is defined as η = − ln(tan θ 2 ). A detailed description of the H1 detector can be found in [11]. Here only the components relevant for this analysis are described. The scattered positron is identified and measured in the SpaCal [12], a lead-scintillating fibre calorimeter situated in the backward region of the H1 detector, covering the pseudorapidity range −4.0 < η < −1.4. The SpaCal also provides information to trigger on the scattered positron in the kinematic region of this analysis. Hits in the backward drift chamber (BDC) are used to improve the identification of the scattered positron and the measurement of its angle [13]. Charged particles emerging from the interaction region are measured by the Central Silicon Track detector (CST) [14] and the Central Tracking Detector (CTD), which covers a range −1.74 < η < 1.74. The CTD comprises two large cylindrical Central Jet drift Chambers (CJCs) and two z-chambers situated concentrically around the beam-line, operated within a solenoidal magnetic field of 1.16 T. The CTD also provides triggering information based on track segments measured in the r-φ-plane of the CJCs and on the z-position of the event vertex obtained from the double layers of two Multi-Wire Proportional Chambers (MWPCs). The tracking detectors are surrounded by a finely segmented Liquid Argon calorimeter (LAr) [15]. It consists of an electromagnetic section with lead absorbers and a hadronic section with steel absorbers and covers the range −1.5 < η < 3.4.
The luminosity determination is based on the measurement of the Bethe-Heitler process ep → epγ, where the photon is detected in a calorimeter close to the beam pipe at z = −103 m.

Data Selection and Analysis
The events selected in this analysis are required to contain a scattered positron in the SpaCal and at least one D * ± meson candidate. The scattered positron must have an energy above 8 GeV.
The virtuality of the photon Q 2 and the inelasticity y are determined from the measured energy E ′ e and the polar angle θ ′ e of the scattered positron via the relations: In addition, the energy W of the γ * p rest-frame is determined using: where s = 4E e E p is the centre-of-mass energy squared of the ep system. The photon virtuality is required to be in the range 2 < Q 2 < 100 GeV 2 . This kinematic range is determined by the geometric acceptance of the SpaCal. The inelasticity is required to lie in the region 0.05 < y < 0.7. The difference between the total energy E and the longitudinal component P z of the total momentum, as calculated from the scattered positron and the hadronic final state, is restricted to 40 < E −P z < 75 GeV. This requirement suppresses photoproduction background, where a hadron is misidentified as the scattered positron. It also reduces the contribution of DIS events with hard initial state photon radiation, where the positron or photon escapes in the negative z-direction. This leads to values of E − P z significantly lower than the expectation 2E e = 55 GeV.
The D * ± mesons are reconstructed from tracks using the decay channel D * + → D 0 π + s → (K − π + )π + s and its charge conjugate, where π s denotes the low momentum pion from the D * ± meson decay. Requirements on the transverse momentum and pseudorapidity of the D * ± meson candidate and its decay products, as well as on particle identification using dE/dx, are similar to those used in previous H1 analyses [16]. A summary of the most important requirements is given in table 1. To select D * ± meson candidates the invariant mass difference method [17] is used. The distribution of ∆M D * ± = M(Kππ s ) − M(Kπ) is shown in figure 1 for the full data sample, together with the wrong charge K ± π ± π ∓ s combinations, using K ± π ± pairs in the accepted D 0 mass range. Detailed studies show that the wrong charge ∆M D * ± distribution provides a good description of the right charge K ∓ π ± π ± s combinatorial background. The signal is extracted using a simultaneous fit to the ∆M D * ± distribution of the right and wrong charge combinations. The signal is fitted using a modified Gaussian function [18] where x = |∆M D * ± − M 0 |/σ. The signal position M 0 and width σ as well as the number of D * ± mesons N D * ± are free parameters of the fit. The background is parametrised as a power function of the form N(a+1)(∆M D * ± −m π ) a /(M max −m π ) a+1 , with the fit boundaries given by the charged pion mass m π and M max = 0.17 GeV. The two free parameters a and N determine the shape and normalisation of the background, respectively. The total event sample is fitted using six free parameters: three for the modified Gaussian, two for the normalisation of the right and wrong charge ∆M D * ± background distributions and one for the background shape, common for the right and wrong charge combinatorial background. In total 2865 ± 89 (stat.) D * ± mesons are obtained. For the differential distributions, the number of D * ± mesons in each measurement bin is extracted using the same procedure, except that the position of the signal peak and its width are fixed to the values determined from the fit to the total sample.
The hadronic final state is reconstructed in each event using an energy flow algorithm. The algorithm combines charged particle tracks and calorimetric energy clusters, taking into account their respective resolution and geometric overlap, into so called hadronic objects while avoiding double counting of energy [19]. The hadronic objects corresponding to the three decay tracks forming the D * ± meson are removed from the event and replaced by one hadronic object having the four-momentum vector of the reconstructed D * ± meson candidate. The energy of the D * ± meson is calculated using M(D * ± ) = 2.010 GeV [20].
Jets are found in the γ * p rest-frame using the inclusive k T cluster algorithm [21] with the distance parameter R = 1 in the η-φ plane. In order to combine hadronic objects into jets, the E-recombination scheme is applied using the four-momenta of the objects. The jet containing the D * ± meson candidate is referred to as the D * ± jet and is required to have a jet transverse energy E * T > 3 GeV in the γ * p rest-frame 1 . According to MC simulations, the D * ± jet is found to be well correlated with the original direction of the charm or anti-charm quark. The distance in azimuth-pseudorapidity, ∆r = ∆η 2 + ∆φ 2 , between the charm quark jet, found using final state partons ("parton level"), and the D * ± jet, found using final state hadrons ("hadron level"), is below 0.3 for 90% of all events. The correlation between the D * ± jet at hadron level and the D * ± jet found using charged particle tracks and calorimetric clusters ("detector level") is even better, since most of the energy of these jets is reconstructed from tracks, which are well measured in the tracking system. The number of D * ± mesons is 1508 ± 68 (stat.) in the D * ± jet sample and 1363 ± 54 (stat.) in the no D * ± jet sample.

Definition of Experimental Observables
A standard method to study fragmentation is to measure the differential production cross section of a heavy hadron (H) as a function of a scaled momentum or energy. In e + e − experiments a customary experimental definition of the scaled energy is z e + e − = E H /E beam , where E beam is the energy of the beams in the centre-of-mass system. In leading order (LO), i.e. without gluon emissions, the beam energy is equal to the energy of the charm or anti-charm quark, which are produced in a colour singlet state. The differential cross section of heavy hadron production as a function of z e + e − is directly related to the fragmentation function.
In the case of ep interactions the situation is more complex. In DIS the dominant process for D * ± meson production at HERA is photon-gluon fusion γ * g → cc [16]. In this case the cc pair is produced in a colour octet state. The energy of the charm quark pair depends on the energy of the incoming photon and gluon. Hadrons produced by initial state gluon emissions and by fragmentation of the proton remnant are also present in the final state.
In this analysis charm fragmentation is studied by measuring the differential cross sections of D * ± meson production as a function of two different observables z hem and z jet defined below, which are sensitive to the fraction of momentum inherited by the D * ± meson from the initial charm quark [22].

The hemisphere method: z = z hem
In the LO photon-gluon fusion process, which dominates charm production at HERA, the charm and anti-charm quarks are moving in the direction of the virtual photon in the γ * p rest-frame of reference. This is due to the fact that the photon is more energetic than the gluon, which typically carries only a small fraction of the proton's momentum. Assuming no further gluon radiation in the initial and final state, the charm and anti-charm quarks are balanced in transverse momentum (figure 2, left). This observation suggests to divide the event into hemispheres, one containing the fragmentation products of the charm quark, the other one those of the anti-charm quark. In order to suppress contributions from initial state radiation and the proton remnant, particles pointing in the proton direction of the γ * p rest-frame (η * < 0) are discarded. The projections of the momenta of the remaining particles onto a plane perpendicular to the γ * paxis are determined. Using the projected momenta, the thrust-axis in this plane, i.e. the axis maximising the sum of the momenta projections onto it, is found. A plane perpendicular to the thrust-axis allows the division of the projected event into two hemispheres, one of them containing the D * ± meson and usually other particles (figure 2, right). The particles belonging to the same hemisphere as the D * ± meson are attributed to the fragmentation of the charm or anti-charm quark. The fragmentation observable is defined as: where in the denominator the energy E * and the momentum P * of all particles of the D * ± meson hemisphere are summed. The longitudinal momentum P * L D * ± is defined with respect to the direction of the three-momentum of the hemisphere, defined as the vectorial sum of the three-momenta of all particles belonging to the hemisphere. The variable z hem is invariant with respect to boosts along the direction of the sum of the momenta of all particles in the hemisphere. Neglecting the mass of the D * ± meson and of the hemisphere, this definition of z hem simplifies to the ratio of their energies.

The jet method: z = z jet
In the case of the jet method the energy and direction of the charm quark are approximated by the energy and direction of the reconstructed jet, which contains the D * ± meson. The fragmentation observable is defined in analogy to z hem as: where the longitudinal momentum P * L D * ± is defined with respect to the direction of the threemomentum of the jet. The jet finding and the determination of z jet are performed in the γ * p rest-frame.
Both fragmentation observables are defined in such a way that they would lead to similar distributions, assuming independent fragmentation and no gluon radiation. The measured distributions, however, are expected to differ, as they have different sensitivities to gluon radiation and charm quarks, which are colour connected to the partons of the proton remnant. The hemisphere method typically includes more energy around the charm quark direction than the jet method. The parameters of fragmentation functions should however be the same, if extracted for a QCD model, which provides a very good description of the underlying physics over the full phase space of this analysis. A comparison of both methods thus may provide a consistency check and a test of the perturbative and non-perturbative physics as encoded in the models.
The measurement is restricted to the regions 0.2 < z hem ≤ 1.0 and 0.3 < z jet ≤ 1.0, as at lower z it is not possible to separate the D * ± meson signal from the large combinatorial background. In order to minimise the sensitivity of the analysis to the total D * ± meson cross section, and to reduce systematic errors, normalised differential cross sections are measured as a function of the fragmentation observables z hem and z jet . The normalisations are chosen such that their integrals over the respective z-regions yield unity.

QCD Models and Calculations
The MC programs RAPGAP and CASCADE are used to generate events containing charm and beauty quarks, which are passed through a detailed simulation of the detector response, based on the GEANT simulation program [23]. They are reconstructed using the same software as used for the data. These event samples are used to determine the acceptance and efficiency of the detector and to estimate the systematic errors associated with the measurements. In addition, these models are fitted to the data in order to determine the parameters of the fragmentation functions.
The Monte Carlo program RAPGAP [7], based on collinear factorisation and DGLAP [24] evolution, is used to generate the direct process of photon-gluon fusion to a heavy cc pair, where the photon acts as a point-like object. In addition, RAPGAP allows the simulation of charm production via resolved processes, where the photon fluctuates into partons, one of which interacts with a parton in the proton, and the remaining partons produce the photon remnant. The program uses LO matrix elements with massive (massless) charm quarks for the direct (resolved) processes. Parton showers based on DGLAP evolution are used to model higher order QCD effects.
The CASCADE program [8] is based on the k T -factorisation approach. Here, the calculation of the photon-gluon fusion matrix element takes into account the charm quark mass and the virtuality and the transverse momentum of the incoming gluon. Gluon radiation off the incoming gluon as well as parton showers off the charm or anti-charm quark are implemented including angular ordering constraints. The gluon density of the proton is evolved according to the CCFM equations [25]. The k T -unintegrated gluon density function A0 [26], extracted from inclusive DIS data, is used.
In both RAPGAP and CASCADE the hadronisation of partons is performed using the Lund string model as implemented in PYTHIA [9]. In the Lund model, the heavy hadron is produced in the process of string breaking. The fraction of the string longitudinal momentum z carried by the hadron is generated according to different choices of adjustable fragmentation functions D H Q (z). Within this analysis three widely used parametrisations are employed, of which two depend on a single free parameter, and one depends on two free parameters. The parametrisation suggested by Peterson et al. [27] has the functional form: and the one by Kartvelishvili et al. [28] is given by: The free parameters ε and α determine the "hardness" of the fragmentation function and are specific to the flavour of the heavy quark, i.e. charm in the case of D * ± meson production. The parametrisation inspired by Bowler and Morris [29] (referred to as the Bowler parametrisation) has the functional form: The shape of the fragmentation function is determined by two free parameters a and b, m Q is the mass of the heavy quark, M T = M 2 H + P 2 T the transverse mass of the heavy hadron, and r Q = 1 as default in PYTHIA.
For data corrections the parameter setting tuned by the ALEPH collaboration [30] together with the Peterson fragmentation function is used for the fragmentation of partons in PYTHIA. It includes higher excited charm states, of which some also decay to D * ± mesons and contribute significantly to the D * ± meson yield. When extracting parameters of the fragmentation functions also the default parameter setting of PYTHIA is used as an alternative. In this case no higher excited charm states are produced. Both parameter settings are indicated in table 2.
The parameters for the Kartvelishvili and Peterson fragmentation functions are also extracted for the HVQDIS program [10]. HVQDIS is based on the NLO, i.e. O(α 2 s ), calculation in the fixed flavour number scheme, with three light active flavours as well as gluons in the proton. The proton parton density functions (PDFs) of the light quarks and the gluon are evolved according to the DGLAP equations. Massive charm quarks are assumed to be produced only perturbatively via photon-gluon fusion and higher order processes. The final state charm quarks are fragmented independently into D * ± mesons in the γ * p rest-frame. Kartvelishvili and Peterson parametrisations are used to generate the charm quark's momentum fraction transferred to the D * ± meson. The energy of the charm quark is calculated using the on-mass-shell condition. In addition, the D * ± meson can be given a transverse momentum P T with respect to the charm quark, according to the function P T exp(−βP T ). The value used for the parameter β corresponds to an average P T (D * ± ) of 350 MeV.
The Monte Carlo programs RAPGAP and HERWIG [31] are used to estimate the size of the hadronisation corrections to the data for comparison with HVQDIS predictions. While the perturbative QCD model of HERWIG is similar to the one of RAPGAP, the HERWIG program employs the cluster hadronisation model, which is quite different from the Lund string model used by PYTHIA.
The basic parameter choices for the QCD models and the NLO calculation are summarised in table 3.

Data Corrections and Systematic Errors
In this analysis, the differential cross section for the production of D * ± mesons, which result from the fragmentation of charm quarks either directly or via decays from higher excited charm states, is measured. The small contribution of D * ± mesons originating from B-hadron decays is estimated with RAPGAP and is subtracted from the data. It is on the level of 1 to 2%. The data are corrected for detector and QED radiative effects. The transverse momentum and pseudorapidity distributions of the D * ± mesons of the Monte Carlo event samples, which are used to correct the data samples for detector effects, are reweighted to the corresponding distributions of the data in order to achieve an improved description. The η and P T dependent reweighting factors differ from unity by typically 10 − 30%. After this reweighting, both RAPGAP and CASCADE provide a good description of the data as shown in figure 3. The description of the no D * ± event sample by the reweighted MC models, as shown in figure 4, is worse.
The measurement bins are defined in such a way that the purity in each bin, defined as the fraction of events reconstructed in a z hem or z jet bin that originate from that bin on hadron level, is between 40 and 70%.
The correction for the detector effects is done using regularised deconvolution, taking into account migrations between measurement bins [32]. The detector response matrix is generated using RAPGAP, and the value of the regularisation parameter is determined through decomposition of the data into eigenvectors of the detector response matrix. As a check, the detector response matrix was also generated using CASCADE and found to be consistent with the one from RAPGAP. Statistical errors are calculated by error propagation using the covariance matrix. The data are then corrected for migrations into the visible phase space using RAPGAP and CASCADE. The effects of QED radiation are corrected for using the HERACLES [33] program, which is interfaced with RAPGAP. Correction factors are calculated from the ratio between cross sections obtained from the model including and not including QED radiation.
The corrections are applied bin-by-bin in z hem and z jet . In the case of z jet the corrections are about 2%. In the case of z hem , where photons radiated into the D * ± meson hemisphere can be mistaken as fragmentation products of the charm quark, the corrections reach 10% for the lowest value of z hem .
In contrast to the QCD models discussed so far, the HVQDIS program provides only a partonic final state with the exception of the additional D * ± meson from charm fragmentation. In order to compare the NLO predictions to the measurements, the data are corrected to the parton level by means of hadronisation corrections, which are estimated using the MC generators RAP-GAP and HERWIG. While the quantity (E * + P * L ) D * ± in equation 4 and 5 is calculated using the momentum of the D * ± meson, the jet finding and the calculation of the jet and hemisphere quantities, the denominators in equations 4 and 5, are performed using the partonic final state. All partons after parton showering are considered, and the same jet and hemisphere finding algorithms are applied at parton and hadron level. For each z-bin the hadronisation correction factor is calculated as the ratio of parton to hadron level cross section. The arithmetic mean of the hadronisation correction factors of both models is used to multiply the data cross section. In case of z hem the hadronisation corrections differ from unity by typically ±40%. For z jet they differ from unity by typically ±20%, except for the highest z-bin, where they are about 50%. The hadronisation corrections as determined by HERWIG and RAPGAP are similar for most of the measurement bins, with the exception of the lowest bin in z jet , where they differ by about 60%.
The following systematic uncertainties on the normalised differential cross sections are considered: • The energy uncertainty of the scattered positron varies linearly from ±3% for an energy of 8 GeV to ±1% for 27 GeV.
• The polar angle of the scattered positron has an estimated uncertainty of ±1 mrad.
• The uncertainty of the energy scale of the hadronic objects is made up of ±0.5% due to tracks and ±4% (±7%) due to LAr (Spacal) clusters.
• The effect of the uncertainty of the tracking efficiency on reconstructing the D * ± meson is determined by changing the nominal efficiency in the simulation as a function of track η and P T . In the central region of the accepted η-P T phase space the estimated uncertainty of the nominal efficiency is ±2%, in the regions of large |η| but not small P T it is ±3%, and for large |η| and small P T it is ±4%.
• The value of dE/dx of the D * ± meson decay products has an estimated uncertainty of ±8%, which is of similar size as the experimental resolution in dE/dx.
• The uncertainty of the D * ± meson signal extraction is estimated using different D * ± counting techniques and by using different fit functions for the background parametrisation. The largest uncertainty comes from the background description, which determines the systematic error on the signal extraction.
• The uncertainty of beauty production by the RAPGAP MC is assumed to be ±100%. The resulting small uncertainty of the normalised D * ± meson cross sections is taken to be symmetrical.
• The effect of using different MC models for the small correction for migrations into the visible phase space is studied using RAPGAP and CASCADE. The factors used to correct the data are determined as the average of the correction factors obtained from the two models. Half the difference is taken as systematic uncertainty.
• For parton level corrected distributions half the difference between the hadronisation correction factors of RAPGAP and HERWIG is taken as the uncertainty due to the different fragmentation models.
Other systematic effects, which are investigated and found to be negligible, are: the effect of reflections, i.e. wrongly or incompletely reconstructed D * ± meson decays, on the shape of the fragmentation observables, the effect on acceptance and reconstruction efficiency from including diffractive events, the effect of using different MC models for the deconvolution of the data and the uncertainty of the QED radiative effects.
Each source of systematic error is varied in the Monte Carlo simulation within its uncertainty. In each measurement bin, the corresponding deviation of the normalised cross sections from the central value is taken as the systematic error. Among the systematic errors the uncertainties due to the scattered positron energy scale, the hadronic energy scale, and the beauty fraction are correlated amongst the bins in z. In the extraction of the parameters of the fragmentation functions, the statistical and systematic errors with their correlations are taken into account. The average effect of various systematic errors on the z hem and z jet distributions is summarised in table 4. Since the distributions of z hem and z jet are normalised, the effect of many systematic uncertainties is reduced and the statistical error dominates the uncertainty of the measurement.

Normalised differential cross sections and comparison with predictions
The differential cross sections of D * ± meson production as a function of the fragmentation observables z hem and z jet are shown in figure 5 for the D * ± jet sample. They refer to the visible phase space given by 2 < Q 2 < 100 GeV 2 , 0.05 < y < 0.7, 1.5 < P T (D * ± ) < 15 GeV and |η(D * ± )| < 1.5. In addition, a D * ± jet with E * T > 3 GeV in the γ * p rest-frame is required in order to have the same hard scale in the event for both distributions z hem and z jet . The measurements and the corresponding predictions are normalised such that their integrals over the respective z-regions yield unity. The normalised cross sections and their errors are given in table 5 for the hemisphere observable and in table 6 for the jet observable. Figure 5 also includes predictions of RAPGAP with three commonly used fragmentation parameter settings for PYTHIA (described in table 2), obtained from e + e − annihilation. The settings and the corresponding values of χ 2 /n.d.f., as calculated from the data and the model predictions, are summarised in table 7. In general, there is reasonable agreement between data and the QCD model with all settings for both the jet and the hemisphere observables. The large difference between the two distributions observed in the highest z jet bin is mainly due to a significant fraction of D * ± jets consisting of a D * ± meson only, for which z jet equals unity. CASCADE provides a similar description of the data as RAPGAP.

Extraction of parameters for the Kartvelishvili and Peterson fragmentation functions
The normalised D * ± meson differential cross sections as a function of z hem and z jet are used to extract optimal parameters for the Peterson and Kartvelishvili fragmentation functions described in section 5. Both parametrisations have a single free parameter.
The parameter extraction is done by comparing different model configurations with the data. A configuration is defined by one of the QCD calculations (RAPGAP, CASCADE or HVQDIS), by one of the fragmentation functions (Peterson or Kartvelishvili) and by a possible value for the corresponding fragmentation parameter, ε or α. For RAPGAP and CASCADE the configuration also depends on the PYTHIA parameter settings used (ALEPH and default, see table 2). In order to be able to compare all configurations to the data, a reweighting procedure is applied. For each of the QCD calculations large event samples with D * ± mesons are generated using the Peterson fragmentation function. For these events the z-value of the fragmentation function, used by the model to generate the fraction of charm quark or string momentum transferred to the D * ± meson, is stored such that each event can be reweighted to another fragmentation function or any other parameter value. For each configuration the predicted and measured distributions of the fragmentation observables are used to determine a χ 2 as a function of the fragmentation parameter. In the calculation of the χ 2 the full covariance matrix is used, taking into account correlated and uncorrelated statistical and systematic errors. The fragmentation parameter is determined at the minimum of the χ 2 . The shape of the χ 2 distribution is used to determine the ±1σ error (using χ 2 min + 1) of the extracted parameter. As an example, in figure 6 the data are compared to the prediction of RAPGAP with the ALEPH setting for PYTHIA as given in table 2 but using the Kartvelishvili parametrisation. The two lines indicate the ±1σ total uncertainty around the best fit value of α. The description of the data by CASCADE is similar.
The parameters α and ε, which are extracted using RAPGAP and CASCADE, with and without higher excited charmed hadrons, are summarised in table 8 together with their corresponding values of χ 2 /n.d.f.. With the fitted parameters the model predictions using either the Peterson or the Kartvelishvili parametrisations describe the data reasonably well, with the Kartvelishvili parametrisation being in all cases slightly preferable, as indicated by the values of χ 2 /n.d.f.. When using the same PYTHIA parameter setting, the fragmentation parameters α and ε, extracted from the z hem and z jet observables, are in good agreement. Both RAPGAP and CASCADE lead to statistically compatible parameters ε and α. A priori, agreement in the fragmentation function parameters for RAPGAP and CASCADE is not required, since the models differ in terms of simulated processes (direct and resolved in case of RAPGAP compared to direct only for CASCADE) and in their implementation of perturbative QCD.
The fragmentation parameters α and ε depend significantly on the PYTHIA parameter settings used, i.e. whether D * ± mesons are assumed to be produced only via direct fragmentation of charm quarks or additionally originate from decays of higher excited charm states. In the latter case the D * ± mesons carry a smaller fraction of the original charm or anti-charm quark momentum in comparison with the directly produced ones. Both the default PYTHIA setting and the setting containing higher excited charm states describe the data equally well. The values of the Peterson parameter ε extracted for the PYTHIA setting containing higher charm states, see table 8, are in agreement with the value ε = 0.04 tuned by ALEPH [30]. This result is consistent with the hypothesis of fragmentation universality in ep and e + e − processes.
The NLO calculation as implemented in HVQDIS with the Kartvelishvili fragmentation function leads to a good fit of the data, corrected for hadronisation effects, as shown in figure 7. On the other hand HVQDIS provides a rather poor description of z hem and z jet , if the Peterson fragmentation function is used (the χ 2 values of the fit are shown in table 8). Simulating a P T of the D * ± meson with respect to the charm quark direction, as explained in section 5, has only a little effect on the extracted value of α.
The distributions of z hem and z jet are also measured in two bins of Q 2 and W . In the case of Q 2 , the bins are defined as 2 < Q 2 < 10 GeV 2 and 10 < Q 2 < 100 GeV 2 . The accessible range in W , determined by the cuts on Q 2 and y, corresponds to 70 < W < 270 GeV, and the bins are defined as 70 < W < 170 GeV and 170 ≤ W < 270 GeV. Correction factors and systematic uncertainties for these samples are determined in the same way as for the full D * ± jet sample. The data are compared to the QCD models with the PYTHIA parameter setting including higher excited charm states and using the Kartvelishvili fragmentation function. For the low and high Q 2 bins the measured distributions are found to be almost the same and well described by the QCD models. The distributions of z jet for the low and high W regions are also similar. A difference is observed for the z hem distribution, which is softer at high W as shown in figure 8. RAPGAP and CASCADE show the same behaviour as a function of W as observed in data. This behaviour can be understood as being partly due to enhanced gluon radiation at high W and partly due to the kinematic effect of the requirement P T (D * ± ) > 1.5 GeV. For events at low W , where the charm quark tends to have smaller energy than at high W , a D * ± meson needs to carry a large fraction of the original quark momentum in order to pass the P T requirement.
The hemisphere observable allows an investigation of charm fragmentation close to the kinematic threshold, at the limit of applicability of the concept of fragmentation functions, by selecting events without a D * ± jet with E * T > 3 GeV. As estimated by MC, the mean centreof-mass energy squared of the γ * g systemŝ for this sample is about 36 GeV 2 , to be compared with about 100 GeV 2 for the D * ± jet sample. This event sample (no D * ± jet sample) has no overlap with the D * ± jet sample investigated so far. The normalised cross section as a function of z hem for the no D * ± jet sample is shown in figure 9 and listed in table 9. Predictions of RAPGAP with the three commonly used fragmentation parameter settings for PYTHIA (see table 2), which provide a reasonable description of the D * ± jet sample (see table 7), fail for this sample. Also the prediction using the fragmentation parameters obtained from the D * ± jet sample is not able to describe these data.
The fragmentation parameters for RAPGAP, CASCADE and the NLO calculation are extracted from the no D * ± jet sample using the same procedure as for the D * ± jet sample. The fit results are summarised in table 10. The predictions of RAPGAP, showing the ±1σ total uncertainty around the fitted value of α are also presented in figure 9. The fragmentation parameters obtained for RAPGAP and CASCADE are statistically compatible. The fragmentation parameters fitted to the no D * ± jet sample are found to be significantly different from those for the D * ± jet sample. They indicate that the fragmentation function for an optimal description of the sample without a D * ± jet needs to be significantly harder than for the D * ± jet sample. The NLO calculation as implemented in HVQDIS fails to describe the no D * ± jet sample as shown in figure 10.
Several parameters of the QCD models, for example those influencing parton showers, have been varied, in trying to describe both samples using the same value for the fragmentation function parameter. However, it was not possible to find MC parameters leading to a consistent fragmentation function for the two samples. Furthermore, the effect of diffractive production of D * ± mesons was not able to explain the difference between the fragmentation parameters observed for the two samples. These investigations indicate that QCD models, together with simple parametrisations of the fragmentation functions, are not able to describe charm fragmentation consistently in the full phase space down to the kinematic threshold.

Conclusions
The fragmentation of charm quarks into D * ± mesons in DIS is studied using the H1 detector at the HERA collider. The normalised D * ± meson differential cross sections are measured as a function of two observables sensitive to fragmentation, the hemisphere observable z hem and the jet observable z jet , in the visible DIS phase space defined by 2 < Q 2 < 100 GeV 2 and 0.05 < y < 0.7, and the D * ± meson phase space defined by 1.5 < P T (D * ± ) < 15 GeV and |η(D * ± )| < 1.5. An additional jet with E * T > 3 GeV, containing the D * ± meson, is required in the γ * p rest-frame in order to provide a hard scale for the events.
The data are compared with predictions of RAPGAP with three widely used PYTHIA parameter settings and the Peterson and the Bowler parametrisations for the fragmentation of heavy flavours obtained from e + e − annihilation. They provide a reasonable description of the ep data presented.
The normalised differential cross sections are used to fit the parameters of the Kartvelishvili and Peterson fragmentation functions within the framework of the QCD models RAPGAP and CASCADE. The fragmentation parameters extracted using the z hem and z jet observables are in good agreement with each other. Both QCD models lead to statistically compatible parameters. The value of the Peterson parameter ε extracted for the PYTHIA parameter setting, which includes not only D * ± mesons from direct fragmentation of charm quarks but also from the decays of higher excited charm states, is in agreement with the value of ε = 0.04 tuned by ALEPH. This result is consistent with the hypothesis of fragmentation universality between ep and e + e − collisions.
The QCD models, with the fragmentation parameters fitted to the data, also provide a good description of the Q 2 and W dependence of the fragmentation observables.
The data, corrected to the parton level, are also compared to the NLO calculation as implemented in HVQDIS, with the addition of independent fragmentation of charm quarks to D * ± mesons. A good fit to the data is obtained when using the fragmentation function by Kartvelishvili et al., while using the one of Peterson et al. results in a poor fit.
Finally, the hemisphere method is used to study the fragmentation of charm produced close to the kinematic threshold, by selecting a data sample fulfilling the nominal requirements on the DIS and D * ± meson phase space, but without a D * ± jet having E * T > 3 GeV in the event. The fragmentation parameters extracted for the QCD models, using this sample of events, are significantly different from those fitted to the D * ± jet sample. Furthermore, the fit for the NLO calculation using the no D * ± jet sample fails. Both observations can be interpreted as an inadequacy of the QCD models and the NLO calculation to provide a consistent description of the full phase space down to the kinematic threshold.
PYTHIA ALEPH Default Description parameter setting setting MSTJ(12) 2 2 baryon model option MSTJ(46) 0 3 parton shower azimut. corr. MSTJ(51) 0 0 Bose-Einstein correlations off PARJ (1) 0.108 0.100 P(qq)/P(q) PARJ (2) 0.286 0.300 P(s)/P(u) PARJ (3) 0.690 0.400 P(us)/P(ud)/P(s)/P(d) PARJ (4) 0.050 0.050 (1/3)P(ud 1)/P(ud 0) PARJ (11) 0.553 0.500 P(S=1)d,u PARJ (12) 0.470 0.600 P(S=1)s PARJ (13) 0.650 0.750 P(S=1)c,b PARJ (14) 0.120 0.000 P(S=0,L=1,J=1) AXIAL PARJ (15) 0.040 0.000 P(S=1,L=1,J=0) SCALAR PARJ (16) 0.120 0.000 P(S=1,L=1,J=1) AXIAL PARJ (17) 0.200 0.000 P(S=1,L=1,J=2) TENSOR PARJ (19) 0.550 1.000 extra baryon suppression PARJ (21) 0.366 0.360 σ q PARJ(25) Fragmentation model Lund string Lund string cluster independent Table 3: Parton density functions (PDFs), fragmentation models and basic parameters used in the QCD models and the NLO calculation. The mass m c of the charm quark is 1.5 GeV in all cases. The transverse momentum of the charm quark in the γ * p rest-frame is given by P * T . The invariant mass squared and the transverse momentum squared of the cc pair are denoted byŝ and Q * 2 T , respectively.   Table 5: Normalised D * ± meson differential cross sections as a function of z hem for the D * ± jet sample, in the visible phase space described in section 7.1. The measurements are normalised such that their integral over the z hem range yields unity. All errors are considered to be symmetric in each bin. For correlated systematic errors a relative sign is indicated.  Table 6: Normalised D * ± meson differential cross sections as a function of z jet for the D * ± jet sample, in the visible phase space described in section 7.1. The measurements are normalised such that their integral over the z jet range yields unity. All errors are considered to be symmetric in each bin. For correlated systematic errors a relative sign is indicated.    Table 9: Normalised D * ± meson differential cross sections as a function of z hem for the no D * ± jet sample, in the visible phase space described in section 7.1. The measurements are normalised such that their integral over the z hem range yields unity. All errors are considered to be symmetric in each bin. For correlated systematic errors a relative sign is indicated.    Figure 4: Comparison on detector level between the no D * ± jet sample and reweighted Monte Carlo models (see section 6) used to correct the data for detector effects for the no D * ± jet sample. Shown are E * T and η * of the D * ± meson hemisphere, calculated from the sum of momenta of all particles in the hemisphere. All quantities are calculated in the γ * p rest-frame.  Figure 5: Normalised D * ± meson cross sections as a function of z jet and z hem for the D * ± jet sample. The measurements are normalised to unity in the displayed range of z jet and z hem , respectively. The data are compared with MC predictions of RAPGAP, using PYTHIA default settings with Peterson or Bowler parametrisations and the ALEPH setting, which includes the production of higher excited charm states (see table 2). The ratio R = MC/data is shown as well as the relative statistical uncertainties (inner error bars) and the relative statistical and systematic uncertainties added in quadrature (outer error bars) for the data points put to R = 1.  Figure 6: Normalised D * ± meson cross sections as a function of z jet and z hem for the D * ± jet sample. The measurements are normalised to unity in the displayed range of z jet and z hem , respectively. The same data as in figure 5 are compared to the predictions of the MC program RAPGAP with the ALEPH setting for PYTHIA and Kartvelishvili parametrisation using the fragmentation function parameter α fitted according to the procedure described in section 7. The full and dashed lines indicate a variation of the fragmentation parameter by ±1σ around the best fit value of α. The ratio R = MC/data is described in the caption of figure 5.  Figure 10: Normalised D * ± meson cross sections as a function of z hem for the no D * ± jet sample. The data are corrected for hadronisation effects (see section 5). The measurements are normalised to unity in the displayed range of z hem . The data are compared to NLO predictions of HVQDIS with the Kartvelishvili fragmentation function using the fragmentation parameter α fitted according to the procedure described in section 7. The full and dashed lines indicate a variation of the fragmentation parameter by ±1σ around the best fit value of α. The ratio R = MC/data is described in the caption of figure 5.