Combination of D * ± Differential Cross-Section Measurements in Deep-Inelastic ep Scattering at HERA

H1 and ZEUS have published single differential cross sections for inclusive D*±-meson production in deep-inelastic electron-proton scattering at HERA from their respective final data sets. The cross sections are here combined in common visible phase space region of photon virtuality Q2 > 5 GeV, electron inelasticity 0.02 < y < 0.7 and D*± meson's transverse momentum pT(D *±) > 1.5 GeV and pseudorapidity |η(D*±)| < 1.5. The combination procedure takes into account all correlations, yielding significantly reduced experimental uncertainties. Double differential cross-sections d2σ/dQ2dy are combined with earlier D*± data, extending the kinematic range down to Q2 > 1.5 GeV . Perturbative next-to-leading-order QCD predictions are compared with the experimental results obtained.


Introduction
Measurements of open charm production in deep-inelastic electron 1 -proton scattering (DIS) at HERA provide important input for stringent tests of the theory of strong interactions, quantum chromodynamics (QCD). Previous measurements [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] have demonstrated that charm quarks are predominantly produced by the boson-gluon-fusion process, γg → cc, whereby charm production becomes sensitive to the gluon distribution in the proton. Measurements have been obtained both from the HERA-I (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) and HERA-II (2003 data-taking periods. Different techniques have been applied at HERA to measure open-charm production in DIS. The full reconstruction of D or D * ± mesons [1-6, 10-12, 15, 16, 18, 20], the long lifetime of heavy flavoured hadrons [7-9, 12, 14, 17, 19] or their semi-leptonic decays [13] are exploited. After extrapolation from the visible to the full phase space, most of these data have already been combined [21] at the level of the reduced cross-sections and have provided a consistent determination of the charm contribution to the proton structure functions, a measurement of the charm-quark mass m c (m c ) and improved predictions for W -and Z-production cross sections at the LHC. However, the extrapolation procedure requires theoretical assumptions, which lead to theoretical uncertainties comparable in size to the experimental uncertainties [21]. Moreover, this combination was restricted to inclusive DIS variables only, such as the photon virtuality, Q 2 , and the inelasticity, y. Alternatively, the measured cross sections can be combined directly in the visible phase 1 In this paper, 'electron' is used to denote both electron and positron if not stated otherwise. space. In this case, dependences on the theoretical input are minimised and the charm production mechanism can be explored in terms of other variables. Such a combination, however, is possible only for data with the same final state, covering a common visible phase space. The recently published differential cross-section measurements by H1 [15,18] and ZEUS [20] for inclusive D * ± -meson production fulfil this requirement. The analysis of fully reconstructed D * ± mesons also offers the best signal-to-background ratio and small statistical uncertainties.
In this paper, visible D * ± -production cross sections [6,15,18,20] at the centre-of-mass energy √ s = 318 GeV are combined such that one consistent HERA data set is obtained and compared directly to differential next-to-leading-order (NLO) QCD predictions. The combination is based on the procedure described elsewhere [21][22][23][24], accounting for all correlations in the uncertainties. This yields a significant reduction of the overall uncertainties of the measurements. The possibility to describe all measurements both in shape and normalisation with a single set of theory parameter values is also investigated and interpreted in terms of future theory improvements.
The paper is organised as follows. In section 2 the theoretical framework is briefly introduced that is used for applying phase-space corrections to the input data sets prior to combination and for providing NLO QCD predictions to be compared to the data. The data samples used for the combination are detailed in section 3 and the combination procedure is described in section 4. The combined single-and double-differential cross sections are presented in section 5 together with a comparison of NLO QCD predictions to the data.

Theoretical predictions
The massive fixed-flavour-number scheme (FFNS) [25][26][27][28] is used for theoretical predictions, since it is the only scheme for which fully differential NLO calculations [29] are available. The cross-section predictions for D * ± production presented in this paper are obtained using the HVQDIS program [29] which provides NLO QCD (O(α 2 s )) calculations in the 3-flavour FFNS for charm and beauty production in DIS. These predictions are used both for small phase-space corrections of the data, due to slightly different binning schemes and kinematic cuts, and for comparison to data.
The following parameters are used in the calculations and are varied within certain limits to estimate the uncertainties associated with the predictions: • The renormalisation and factorisation scales are taken as µ r = µ f = Q 2 + 4m 2 c . The scales are varied simultaneously up or down by a factor of two for the phasespace corrections where only the shape of the differential cross sections is relevant. For absolute predictions, the scales are changed independently to 0.5 and 2 times their nominal value.
• The pole mass of the charm quark is set to m c = 1.50 ± 0.15 GeV. This variation also affects the values of the renormalisation and factorisation scales.
• For the strong coupling constant the value α • The proton parton density functions (PDFs) are described by a series of FFNS variants of the HERAPDF1.0 set [24] at NLO determined within the HERAFitter [30] framework, similar to those used in the charm combination paper [21]. Charm measurements were not included in the determination of these PDF sets. For all parameter settings used here, the corresponding PDF set is used. By default, the scales for the charm contribution to the inclusive data in the PDF determination were chosen to be consistent with the factorisation scale used in HVQDIS, while the renormalisation scale in HVQDIS was decoupled from the scale used in the PDF extraction, except in the cases where the factorisation and renormalisation scales were varied simultaneously. As a cross check, the renormalisation scales for both heavy-and light-quark contributions are varied simultaneously in HVQDIS and in the PDF determination, keeping the factorisation scales fixed. The result lies well within the quoted uncertainties. The cross sections are also evaluated with 3-flavour NLO versions of the ABM [31] and MSTW [32] PDF sets. The differences are found to be negligible compared to those from varying other parameters, such that no attempt for coverage of all possible PDFs is made.
The NLO calculation performed by the HVQDIS program yields differential cross sections for charm quarks. These predictions are converted to D * ± -meson cross sections by applying the fragmentation model described in a previous publication [21]. This model is based on the fragmentation function of Kartvelishvili et al. [33] which provides a probability density function for the fraction of the charm-quark momentum transferred to the D * ± meson. The function is controlled by a single fragmentation parameter, α K . Different values of α K [21] are used for different regions of the invariant mass,ŝ, of the photon-parton centre-of-mass system. The boundaryŝ 1 = 70 ± 40 GeV 2 between the first two regions is one of the parameter variations. The boundaryŝ 2 = 324 GeV 2 between the second and third region remains fixed. The model also implements a transverse-fragmentation component by assigning to the D * ± meson a transverse momentum, k T , with respect to the charm-quark direction [21]. The following parameters are used in the calculations together with the corresponding variations for estimating the uncertainties of the NLO predictions related to fragmentation: • The fragmentation parameter α K , the bin boundaryŝ 1 and the average k T are varied according to a prescription described elsewhere [21].
The small beauty contribution to the D * ± signal needs a detailed treatment of the B hadron decay to D * ± mesons and is therefore obtained from NLO QCD predictions for beauty hadrons convoluted with decay tables to D * ± mesons and decay kinematics obtained from EvtGen [35]. The parameters for the calculations and the uncertainties are: • The renormalisation and factorisation scales µ r = µ f = Q 2 + 4m 2 b are varied in the same way as described above for charm. The variations are applied simultaneously for the calculation of the charm and beauty cross-section uncertainties. • The fragmentation model for beauty quarks is based on the Peterson et al. [36] parametrisation using b = 0.0035 ± 0.0020 [37].
• The proton structure is described by the same PDF set (3-flavour scheme) used for the charm cross-section predictions.
The total theoretical uncertainties are obtained by adding all individual contributions in quadrature.

Data samples for cross-section combinations
The H1 [39][40][41] and ZEUS [42] detectors were general purpose instruments which consisted of tracking systems surrounded by electromagnetic and hadronic calorimeters and muon detectors. The most important detector components for the measurements combined in this paper are the central tracking detectors (CTD) operated inside solenoidal magnetic fields of 1.16 T (H1) and 1.43 T (ZEUS) and the electromagnetic sections of the calorimeters. The CTD of H1 [40] (ZEUS [43][44][45]) measured charged particle trajectories in the polar angular range 2 of 15 • < Θ < 165(164) • . In both detectors the CTDs were complemented with highresolution silicon vertex detectors: a system of three silicon detectors for H1, consisting of the Backward Silicon Tracker [46], the Central Silicon Tracker [47] and the Forward Silicon Tracker [48], and the Micro Vertex Detector [49] for ZEUS. For charged particles passing through all active layers of the silicon vertex detectors and CTDs, transverse-momentum resolutions of σ(p T )/p T 0.002p T / ⊕ 0.015 (H1) and σ(p T )/p T 0.0029p T / ⊕ 0.0081 ⊕ 0.0012/p T (ZEUS), with p T in units of GeV, have been achieved.
Each of the central tracking detectors was enclosed by a set of calorimeters comprising an inner electromagnetic and an outer hadronic section. The H1 calorimeter system consisted of the Liquid Argon calorimeter (LAr) [50] and the backward lead-scintillator calorimeter (SpaCal) [41] while the ZEUS detector was equipped with a compensating uranium-scintillator calorimeter (CAL) [51][52][53][54]. Most important for the analyses combined in this paper is the electromagnetic part of the calorimeters which is used to identify and measure the scattered electron. Electromagnetic energy resolutions σ(E)/E of 0.11/ √ E (LAr) [55], 0.07/ √ E (SpaCal) [56] and 0.18/ √ E (CAL), with E in units of GeV, were achieved.
The Bethe-Heitler process, ep → eγp, is used by both experiments to determine the luminosity. Photons originating from this reaction were detected by photon taggers at about 100 m downstream of the electron beam line. The integrated luminosities are known 2 In both experiments a right-handed coordinate system is employed with the Z axis pointing in the nominal proton-beam direction, referred to as "forward direction", and the X axis pointing towards the centre of HERA. The origin of the coordinate system is defined by the nominal interaction point in the case of H1 and by the centre of the CTD in the case of ZEUS.  Table 1. Data sets used in the combination. For each data set the respective kinematic range and the integrated luminosity, L, are given.
Combinations are made for single-and double-differential cross sections. In table 1 the datasets 3 used for these combinations are listed together with their visible phasespace regions and integrated luminosities. The datasets I-III are used to determine singledifferential combined cross sections as a function of the D * ± meson's transverse momentum, p T (D * ), pseudorapidity, η(D * ), and inelasticity, z(D * ) = (E(D * ) − p Z (D * ))/(2E e y), measured in the laboratory frame, and of Q 2 and y. The variables E(D * ), p Z (D * ) and E e denote the energy of the D * ± meson, the Z component of the momentum of the D * ± meson and the incoming electron energy, respectively. Owing to beam-line modifications related to the HERA-II high-luminosity running [58] the visible phase space of these cross sections at HERA-II is restricted to Q 2 > 5 GeV 2 , which prevents a combination with earlier D * ± cross-section measurements for which the phase space extends down to Q 2 = 1.5 GeV 2 .
In the case of the double-differential cross section, d 2 σ/dydQ 2 , the kinematic range can be extended to lower Q 2 using HERA-I measurements [4,6,10]. In order to minimise the use of correction factors derived from theoretical calculations, the binning scheme of such measurements has to be similar to that used for the HERA-II data. One of the HERA-I measurements, set IV of table 1, satisfies this requirement and is therefore included in the combination of this double-differential cross section. The visible phase spaces of the combined single-and double-differential cross sections are summarised in table 2.
The measurements to be combined for the single-and double-differential cross sections are already corrected to the Born level with a running fine-structure constant α and include both the charm and beauty contributions to D * ± production. The total expected beauty contribution is small, varying from ∼ 1% at the lowest Q 2 to ∼ 7% at the highest Q 2 . The cross sections measured previously [6,15,18] are here corrected to the PDG value [38] of the D 0 branching ratio.

Treatment of data sets for single-differential cross sections
In order to make the input data sets compatible with the phase space quoted in table 2 and with each other, the following corrections are applied before the combination: 3 Of the two sets of measurements in [18], that compatible with the above cuts is chosen. Table 2. Visible phase space of the combined cross sections.
• The H1 collaboration has published measurements of D * ± cross sections separately for 5 GeV 2 < Q 2 < 100 GeV 2 (set I) and for 100 GeV 2 < Q 2 < 1000 GeV 2 (set II). Due to the limited statistics at high Q 2 , a coarser binning in p T (D * ), η(D * ), z(D * ) and y was used in set II compared to set I. Therefore the cross section in a bin i of a given observable integrated in the range 5 GeV 2 < Q 2 < 1000 GeV 2 is calculated according to . (3.1) Here σ int denotes the integrated visible cross section and σ NLO stands for the NLO predictions obtained from HVQDIS. In this calculation both the experimental uncertainties of the visible cross section at high Q 2 and the theoretical uncertainties as described in section 2 are included. The contribution from the region 100 GeV 2 < Q 2 < 1000 GeV 2 to the full Q 2 range amounts to 4% on average and reaches up to 50% at highest p T (D * ).
• The bin boundaries used for the differential cross section as a function of Q 2 differ between sets I, II and set III. At low Q 2 this is solved by combining the crosssection measurements of the first two bins of set I into a single bin. For Q 2 > 100 GeV 2 no consistent binning scheme could be defined directly from the singledifferential cross-section measurements. However, the measurements of the doubledifferential cross sections d 2 σ/dQ 2 dy have been performed in a common binning scheme. By integrating these cross sections in y, single-differential cross sections in Q 2 are obtained also for Q 2 > 100 GeV 2 which can be used directly in the combination.
• The cross-section measurements in set III are restricted to p T (D * ) < 20 GeV while there is no such limitation in the phase space of the combination. Therefore these cross sections are corrected for the contribution from the range p T (D * ) > 20 GeV using HVQDIS. This correction is found to be less than 0.1%.

Treatment of data sets for double-differential cross sections
Since the restriction to the same phase space in Q 2 does not apply for the combination of the double-differential cross sections in Q 2 and y, the HERA-I measurement, set IV, is also included in the combination. This allows an extension of the kinematic range down to Q 2 > 1.5 GeV 2 . The p T (D * ) ranges of the measurements of sets III and IV are corrected in the same way as for the single-differential cross sections.
To make the binning scheme of the HERA-I measurement compatible with that used for the HERA-II datasets, the binning for all datasets is revised. Cross sections in the new bins are obtained from the original bins using the shape of the HVQDIS predictions as described in section 2. The new binning is given in section 5 (table 9). Bins are kept only if they satisfied both of the following criteria: • The predicted fraction of the cross section of the original bin contained in the kinematical overlap region in Q 2 and y between the original and corrected bins is greater than 50% (in most bins it is greater than 90%).
• The theoretical uncertainty from the correction procedure is obtained by evaluating all uncertainties discussed in section 2 and adding them in quadrature. The ratio of the theoretical uncertainty to the uncorrelated experimental uncertainty is required to be less than 30%.
This procedure ensures that the effect of the theoretical uncertainties on the combined data points is small. Most of the HERA-II bins are left unmodified; all of them satisfied the criteria and are kept. Out of the 31 original HERA-I bins, 26 bins satisfy the criteria and are kept. The data points removed from the combination mainly correspond to the low-y region where larger bins were used for the HERA-I data.

Combination method
The combination of the data sets uses the χ 2 minimisation method developed for the combination of inclusive DIS cross sections [22,24], as implemented in the HERAverager program [59]. For an individual dataset e the contribution to the χ 2 function is defined as (4.1) Here µ i,e is the measured value of the cross section in bin i and γ i,e j , δ i,e,stat and δ i,e,uncor are the relative correlated systematic, relative statistical and relative uncorrelated systematic uncertainties, respectively, from the original measurements. The quantities m i express the values of the expected combined cross section for each bin i and the quantities b j express the shifts of the correlated systematic-uncertainty sources j, in units of the standard deviation. Several data sets providing a number of measurements (index e) are represented by a total χ 2 function, which is built from the sum of the χ 2 exp,e functions of all data sets The combined cross sections m i are obtained by the minimisation of χ 2 tot with respect to m i and b j .
The averaging procedure also provides the covariance matrix of the m i and the uncertainties of the b j at the minimum. The b j at the minimum and their uncertainties are referred to as "shift" and "reduction", respectively. The covariances V of the m i are given in the form V = V uncor + k V k sys [23]. The matrix V uncor is diagonal. Its diagonal elements correspond to the covariances obtained in a weighted average performed in the absence of any correlated systematic uncertainties. The covariance matrix contributions V k sys correspond to correlated systematic uncertainties on the averaged cross sections, such that the elements of a matrix V k sys are obtained as (V k sys ) ij = δ sys,k i δ sys,k j , given a vector δ sys,k of systematic uncertainties. It is worth noting that, in this representation of the covariance matrix, the number of correlated systematic sources is identical to the number of correlated systematic sources in the input data sets.
In the present analysis, the correlated and uncorrelated systematic uncertainties are predominantly of multiplicative nature, i.e. they change proportionally to the central values.
In equation (4.1) the multiplicative nature of these uncertainties is taken into account by multiplying the relative errors γ i,e j and δ i,e,uncor by the cross-section expectation m i . In charm analyses the statistical uncertainty is mainly background dominated. Therefore it is treated as being independent of m i . For the minimisation of χ 2 tot an iterative procedure is used as described elsewhere [23].
The 55 systematic uncertainties obtained from the original publications were examined for their correlations. Within each data set, most of the systematic uncertainties are found to be point-to-point correlated, and are thus treated as fully correlated in the combination. In total there are 23 correlated experimental systematic sources and 5 theory-related uncertainty sources. A few are found to be uncorrelated and added in quadrature. For the combination of single-differential cross sections the uncorrelated uncertainties also include a theory-related uncertainty from the corrections discussed in section 3, which varies between 0 and 10% of the total uncertainty and is added in quadrature. Asymmetric systematic uncertainties were symmetrised to the larger deviation before performing the combination. Except for the branching-ratio uncertainty, which was treated as correlated, all experimental systematic uncertainties were treated as independent between the H1 and ZEUS data sets. Since the distributions in p T (D * ), η(D * ), z(D * ), Q 2 and y are not statistically independent, each distribution is combined separately.

Combined cross sections
The results of combining the HERA-II measurements [15,18,20] as a function of p T (D * ), η(D * ), z(D * ), Q 2 and y are given in tables 3-7, together with their uncorrelated and correlated uncertainties. 4 The total uncertainties are obtained by adding the uncorrelated and correlated uncertainties in quadrature.   Table 4. The combined differential D * ± -production cross section in the phase space given in table 2 as a function of η(D * ), with its uncorrelated (δ uncor ), correlated (δ cor ) and total (δ tot ) uncertainties.   Table 7. The combined differential D * ± -production cross section in the phase space given in table 2 as a function of y, with its uncorrelated (δ uncor ), correlated (δ cor ) and total (δ tot ) uncertainties.

JHEP09(2015)149
The individual data sets and the results of the combination are shown in figures 1-5. The consistency of the data sets as well as the reduction of the uncertainties are illustrated further by the insets at the bottom of figures 1 and 4. The combinations in the different variables have a χ 2 probability varying between 15% and 87%, i.e. the data sets are consistent. The systematic shift between the two input data sets is covered by the respective correlated uncertainties. The shifts and reductions of the correlated uncertainties are given in table 8. The improvement of the total correlated uncertainty is due to small reductions of many sources. While the effective doubling of the statistics of the combined result reduces the uncorrelated uncertainties, the correlated uncertainties of the combined cross sections are reduced through cross-calibration effects between the two experiments. Typically, both effects contribute about equally to the reduction of the total uncertainty.
The combined cross sections as a function of p T (D * ), η(D * ), z(D * ), Q 2 and y are compared to NLO predictions 5 in figures 6-10. In general, the predictions describe the data well. The data reach an overall precision of about 5% over a large fraction of the measured phase space, while the typical theoretical uncertainty ranges from 30% at low Q 2 to 10% at high Q 2 . The data points in the different distributions are statistically and systematically correlated. No attempt is made in this paper to quantify the correlations between bins taken from two different distributions. Thus quantitative comparisons of theory to data can only be made for individual distributions.
In order to study the impact of the current theoretical uncertainties in more detail, the effect of some variations on the predictions is shown separately in figure 11, compared to the same data as in figures 6, 8 and 10. Only the variations with the largest impact on the respective distribution are shown in each case.
1. The NLO prediction as a function of p T (D * ) (figure 11, top) describes the data better 5 The NLO QCD prediction for the beauty contribution to D * ± production, calculated as described in section 2, can be found on http://www.desy.de/h1zeus/dstar2015/.         Table 9. The combined double-differential D * ± -production cross section in the phase space given in table 2 as a function of Q 2 and y, with its uncorrelated (δ uncor ), correlated (δ cor ) and total (δ tot ) uncertainties.

JHEP09(2015)149
by either • setting the charm-quark pole mass to 1.35 GeV or • reducing the renormalisation scale by a factor 2 or • increasing the factorisation scale by a factor 2.
Simultaneous variation of both scales in the same direction would largely compensate and would therefore have a much smaller effect.
2. The prediction for the z(D * ) distribution ( figure 11, bottom left) describes the shape of the data better if the fragmentation parameters are adjusted such that the boundary between the two lowest fragmentation regions [21] is varied from the default of 70 GeV 2 to its lower boundary of 30 GeV 2 .
3. The preference for a reduced renormalisation scale already observed for p T (D * ) is confirmed by the z(D * ) distribution (figure 11, bottom right). However, the shape of the z(D * ) distribution rather prefers variations of the charm mass and the factorisation scale in the opposite direction to those found for the p T (D * ) distribution. The distributions of the other kinematic variables do not provide additional information to these findings [60].
As stated before, within the large uncertainties indicated by the theory bands in figures 1-5, all distributions are reasonably well described. However, the above study shows that the different contributions to these uncertainties do not only affect the normalisation, but also change the shape of different distributions in different ways. It is therefore not obvious that a variant of the prediction that gives a good description of the distribution in one variable will also give a good description of the distribution in another.
Based on the study of items 1.-3. above, a 'customised' calculation is performed to demonstrate the possibility of obtaining an improved description of the data in all variables at the same time, both in shape and normalisation, within the theoretical uncertainties quoted in section 2. For this calculation, the following choices were made: • From the three options discussed in item 1. above, the second is chosen, i.e. the renormalisation scale is reduced by a factor 2, with the factorisation scale unchanged.
• The change of the fragmentation parameterŝ 1 = 30 GeV 2 , as discussed in item 2. above, is applied.
• At this stage, the resulting distributions are still found to underestimate the data normalisation. As the renormalisation and factorisation scales are recommended to differ by at most a factor of two [61], the only significant remaining handle arising from items 1. and 3. is the charm-quark pole mass. This mass is set to 1.4 GeV, a value which is also compatible with the partially overlapping data used for a previous dedicated study [21] of the charm-quark mass.
• All other parameters, which have a much smaller impact [60] than those discussed above, are left at their central settings as described in section 2.

JHEP09(2015)149
The result of this customised calculation is indicated as a dotted line in figures 6-10. A reasonable agreement with data is achieved simultaneously in all variables. This a posteriori adjustment of theory parameters is not a prediction, but it can be taken as a hint in which direction theoretical and phenomenological developments could proceed. The strong improvement of the description of the data relative to the central prediction through the customisation of the renormalisation scale is in line with the expectation that higher-order calculations will be helpful to obtain a more stringent statement concerning the agreement of perturbative QCD predictions with the data. The improvement from the customisation of one of the fragmentation parameters and the still not fully satisfactory description of the z(D * ) distribution indicate that further dedicated experimental and theoretical studies of the fragmentation treatment might be helpful.
In general, the precise single-differential distributions resulting from this combination, in particular those as a function of p T (D * ), η(D * ) and z(D * ), are sensitive to theoretical and phenomenological parameters in a way which complements the sensitivity of more inclusive variables like Q 2 and y.
The combined double-differential cross sections with the uncorrelated, correlated and total uncertainties 6 as a function of Q 2 and y are given in table 9. The total uncertainty is obtained by adding the uncorrelated and correlated uncertainties in quadrature. The individual data sets as well as the results of the combination are shown in figure 12. Including data set IV slightly reduces the overall cross-section normalisation with respect to the combination of sets I-III only. The pull distribution of the combination is shown in figure 13. The combination has a χ 2 probability of 84%, i.e. all data sets are consistent. The shifts and reductions of the correlated uncertainties are given in table 8.
These combined cross sections are compared to NLO predictions 7 in figure 14. The customised calculation is also shown. In general the predictions describe the data well. The data have a precision of about 5-10% over a large fraction of the measured phase space, while the estimated theoretical uncertainty ranges from 30% at low Q 2 to 10% at high Q 2 . As well as the single-differential distributions, these double-differential distributions give extra input to test further theory improvements.

Conclusions
Measurements of D * ± -production cross sections in deep-inelastic ep scattering by the H1 and ZEUS experiments are combined at the level of visible cross sections, accounting for their systematic correlations. The data sets were found to be consistent and the combined data have significantly reduced uncertainties. In contrast to the earlier charm combination at the level of reduced cross sections, the present combination does not have significant theory-related uncertainties and in addition distributions of kinematic variables of the D * ± mesons are obtained. The combined data are compared to NLO QCD predictions. 6 A detailed breakdown of correlated uncertainties can be found on http://www.desy.de/h1zeus/ dstar2015/. 7 The NLO QCD prediction for the beauty contribution to D * ± production, calculated as described in section 2, can be found on http://www.desy.de/h1zeus/dstar2015/. The predictions describe the data well within their uncertainties. Higher order calculations would be helpful to reduce the theory uncertainty to a level more comparable with the data precision. Further improvements in the treatment of heavy-quark fragmentation would also be desirable.