1 Introduction

The production of D mesons in inelastic scattering of 160 GeV/c muons on nucleons μN→μ′DX is assumed to be dominated by a process where the exchanged virtual photon γ ∗ fuses with a gluon g into a charm anti-charm quark pair, \(\gamma^{*} g \rightarrow c \bar{c}\). The cross-section \(\sigma^{\gamma^{*} g \rightarrow c \bar{c}}\) of this photon–gluon fusion (PGF) process and its dependence on the relative polarization of photon and gluon can be calculated in perturbative QCD [1–8]. Thus, using polarized muons and polarized nucleons, a measurement of the photon nucleon cross-section asymmetry \(\varDelta\sigma^{\gamma^{*} N \rightarrow c \bar{c}X} / \sigma^{\gamma^{*} N \rightarrow c \bar{c}X}\) allows the determination of the gluon polarization Δg/g in the nucleon. With this objective, open charm production has been studied in the COMPASS experiment at CERN for longitudinally polarized muons interacting with longitudinally polarized deuterons. The incoming muon energy of 160 GeV was chosen, since the cross-section difference \(\varDelta\sigma^{\gamma^{*} N \rightarrow c \bar{c}X} \) for parallel and anti-parallel spins of photon and nucleon reaches a maximum for virtual photon energies around 80 GeV according to most models for the gluon helicity distribution function Δg(x Bj ,Q 2), and the polarization transfer from muon to virtual photon is large in the relevant photon energy range.

Final states, where the decays D 0→K − π + or D ∗+→D 0 π +→K − π + π + or the charge conjugate decays are detected, are chosen in order to achieve the best possible combination of mass resolution, signal-over-background ratio and signal statistics. Based on samples of events with these final states, extracted from data taken during the years 2002–2006, COMPASS has published results for 〈Δg/g〉 [9].

The photon–gluon cross-section asymmetry \(a_{LL} = \varDelta\sigma^{\gamma^{*} g \rightarrow c \bar{c}} / \sigma^{\gamma^{*} g \rightarrow c \bar{c}}\) needed for extracting 〈Δg/g〉 is estimated making two assumptions: only PGF contributes as calculated in leading order QCD and charm and anti-charm quarks hadronize independently of the target polarization. The parton kinematics are estimated event-by-event on the basis of the observed 3-momentum of the D 0 meson and the momentum difference of the incoming and the scattered muon using a parametrisation based on the Monte Carlo event generator AROMA with default charm quark fragmentation [10].

However, production mechanisms other than PGF with standard charm quark fragmentation may contribute to the observed events with charmed mesons. The interaction of the virtual photon with intrinsic charm [11–13] is one possible competing mechanism. The associated production of \(\varLambda_{c} \bar{D}^{0}\) [14] or, more generally, asymmetric hadronization of c and \(\bar{c}\) like in the Dual Parton model with a meson and a baryon string [15, 16] may play an important role in some regions of phase space. The unpolarized phase space distribution and semi-inclusive differential cross-sections of D 0 and D ∗ mesons contains information about the contributions of different production mechanisms. Therefore, a detailed study has been performed within the acceptance of the COMPASS spectrometer, averaging over the two polarizations of the target [17].

At HERA, i.e. at much larger center-of-mass energy, charm electro-production has been studied in detail by the H1 and ZEUS Collaborations, see [18–20] and references therein. In addition to the PGF other production mechanisms also contribute in this case, like gluon-gluon fusion to \(c \bar{c}\) from a resolved photon. The hadronizations of c and \(\bar{c}\) can more safely be assumed to be independent. COMPASS covers a complementary kinematic region with virtual photon energies in the range from threshold at about 30 GeV up to 140 GeV in the laboratory frame. Prior to COMPASS, this energy range was covered only by the EMC experiment [21], which collected about 90 D 0 meson events produced by deep inelastic scattering of 240 or 280 GeV muons on hydrogen and deuterium targets for a study of the charm production mechanism. Only one charm photo-production experiment explored the region close to threshold [22], while two concentrated on high energy photons [16, 23].

The present article shows details of the phase space distributions of D meson production as a function of various kinematic variables: the energy \(\nu= E_{\mu} - E_{\mu}^{\prime}\) of the exchanged virtual photon γ ∗ (assuming single photon exchange) with four momentum q=p−p′, the inelasticity y=ν/E μ of the event, the negative invariant γ ∗ mass squared Q 2=−q 2=−(p−p′)2 and the Bjorken scaling variable x Bj =Q 2/2P⋅q=Q 2/2M p ν. Here E μ and \(E_{\mu}^{\prime}\) are the laboratory energies, p and p′ the 4-momentum vectors of the incoming and scattered muon respectively, P is the 4-momentum of the target nucleon and M p is the proton mass.

In order to describe both D ∗ and D 0 meson production, the following kinematic variables are used: the transverse momentum p T of the Kπ pair (from the D 0 decay) with respect to the γ ∗ direction, the D 0 energy E in the laboratory system and the energy fraction z=E/ν.

The outline of the paper is as follows: after a brief overview of the experimental set-up (Sect. 2) and a detailed description of the data selection procedure (Sect. 3), the methods of signal extraction are described in Sect. 4. The kinematic distributions of events from signal and background regions are shown in Sect. 5. They are based on the entire available data sample collected during the years 2002–2006, and are not corrected for acceptance. The purpose of this section is to compare the distributions of open charm to those of background events. Section 6 describes the acceptance correction and luminosity calculation needed to extract the total and differential semi-inclusive cross sections for charm production. They are performed for the 2004 data only. In Sect. 7, the differential cross-sections as a function of the various kinematic variables and the total cross-section obtained by integration are shown and compared with the AROMA predictions for the production of D mesons by PGF. Significant differences between D ∗+ and D ∗− meson production are observed for the acceptance-corrected data from 2004. A statistically more precise comparison of D ∗+ and D ∗− production is based on the entire data sample (2002–2006). Particle-antiparticle asymmetries are determined under the assumption, which is verified for the 2004 data, that the D ∗+ and D ∗− acceptances are equal to a good approximation.

2 Experimental setup

The data were taken using the COMPASS spectrometer situated at the M2 beam line at the CERN Super Proton Synchrotron. A detailed description of the COMPASS spectrometer can be found in Ref. [24].

The momentum of the positive muon beam is about 160 GeV/c with a spread of 5 %. The momentum of each incoming muon is measured with a precision of Δp/p<1 % in the beam momentum spectrometer located upstream of the experimental hall, and its direction is measured with a precision of 30 μrad with a detector telescope in front of the target.

The (polarized) 6LiD target consisted of two 60 cm long cells during the years 2002–2004 and of 3 cells with a total length of 120 cm in 2006. The polarization is reversed regularly such that the products of integrated luminosities times acceptance are equal for both polarizations. The sum of both corresponds to the unpolarized case. Hence unpolarized distributions, which are the subject of the present analysis, are obtained from the sum of all data. The target is housed in a superconducting solenoid magnet, which determines the angular acceptance of the spectrometer. The acceptance in the polar angle, measured at the upstream edge of the target, was 70 mrad in 2002–2004, while with the new target magnet in 2006 it was increased to 180 mrad.

The 2-stage spectrometer is designed to reconstruct and identify the scattered muon and produced hadrons over a wide momentum range. It contains a large angle (LAS) and a small angle (SAS) part, each part equipped with a dipole magnet. Tracking detectors are located in front and behind each magnet, and electromagnetic and hadron calorimeters behind. The LAS covers polar laboratory angles from about 15 mrad up to 70 mrad in 2002–2004 and, with the new target magnet, up to 180 mrad in 2006. The SAS covers the polar laboratory angles below 20 mrad.

Particle tracking is done with a large variety of tracking detectors: several stations of silicon microstrip detectors, scintillating fiber detectors, high resolution micromesh gaseous chambers, gas electron-multiplier chambers, drift chambers, large area straw drift chambers, multiwire proportional chambers and muon drift tubes. The scattered muons are identified downstream of additional hadron absorbers placed behind the hadron calorimeters. Charged hadrons are identified by a Ring Imaging Cherenkov detector (RICH) in the LAS.

The trigger system [25] uses hodoscope and calorimeter information to select inelastic muon interactions with minimum bias. The overall trigger and muon track reconstruction efficiency is in the range 60 % to 80 % for most of the kinematic region covered by COMPASS.

3 Data selection

The total number of events with an incoming muon (140 GeV/c<p μ <180 GeV/c) and a scattered muon from a common vertex is 5.2×109, which corresponds to an integrated luminosity of about 2.8 fb−1. This sample is used to search for \(D^{0}(\bar{D}^{0})\) and D ∗± mesons. A fiducial volume cut makes sure that the extrapolated incoming muon trajectory traverses all target cells and that the primary vertex is located within the volume of one of the target cells.

Since the COMPASS experiment uses a large solid target, the detection of a secondary decay vertex, which is a standard method in open charm detection, is excluded and the selection of \(D^{0}(\bar{D}^{0})\) and D ∗± mesons relies on requirements on event kinematics and particle identification. The event selection is the same as used in the previous COMPASS open charm publication [9] except for stricter requirements on the selection of the incoming muon.

Cuts used to select D 0 originating from the decay of a D ∗ (referred to as D ∗ or ‘tagged’ D 0 sample) slightly differ from those used to select directly produced D 0 mesons (referred to as ‘untagged’ D 0 sample). An event from the untagged D 0 sample contains at least one candidate for the 2-body decay D 0→K − π + or its charge conjugate (c.c.), while in the tagged D 0 sample a slow pion from the decay chain \(D^{*+} \rightarrow D^{0} \pi^{+}_{s} \rightarrow K^{-}\pi^{+}\pi^{+}_{s}\) (or c.c.) has to be present in addition.

Particles are identified by using the RICH. All tracks with momentum measured in one or both spectrometer stages and falling within the geometrical acceptance of the RICH are used to calculate the likelihoods (LKs) that the Cherenkov photons detected by the RICH are due to electron, muon, pion, kaon, proton, or background. The LK for a specific particle is calculated only if the particle velocity is above the threshold for the emission of Cherenkov photons in the radiator gas. This threshold depends on the refractive index that is extracted from the data on a run-by-run basis. For pions, kaons, and protons, this gives an average momentum threshold of 2.5, 8.9 and 16.9 GeV/c respectively. At large momenta pions and kaons cannot be efficiently separated, thus it is required that the momentum of the particle is below 50 GeV/c. In the tagged D 0 sample, due to the small mass difference between the D ∗(2010) and the D 0(1865), only a limited energy is available for the pion produced in the D ∗→D 0 π s decay. In this case, the π s candidate must not have been identified as an electron by the RICH. Details on the LK requirements and the use of the RICH information can be found in Ref. [17].

For untagged D 0, the following cuts are applied to the K − π + and K + π − combinations: p π >7 GeV/c, z>0.2, |M(Kπ)−M(D 0)|<700 MeV and |cosθ K |<0.65, where θ K is the decay kaon angle in the D 0 center-of-mass system with respect to the D 0 direction of flight.

For the tagged D 0, the Kππ s invariant mass is calculated only if the Kπ system has an invariant mass in the range |M(Kπ)−M(D 0)|<700 MeV. The distribution of ΔM=M(Kππ s )−M(Kπ)−M(π) as a function of M(Kπ) is shown in Fig. 1. Here a clear spot for the D ∗ is visible at ΔM∼6 MeV in the region of the D 0 mass. The cut 3.2 MeV<ΔM<8.9 MeV improves the D 0 signal with respect to the combinatorial background by more than an order of magnitude. The Kπ system is also required to have z>0.2 and |cosθ K |<0.9.

Fig. 1
figure 1

Scatter plot for D ∗ candidates before applying the ΔM cut. Vertical axis: ΔM; horizontal axis: M(Kπ). The accumulation of events around the D 0 nominal mass of 1.864 GeV and ΔM=6.1 MeV corresponds to the decay sequence D ∗→π s D 0→π s (Kπ)

Full size image

These sets of cuts define the untagged and tagged D 0 samples, i.e. the D 0 and D ∗ candidates.

The cosθ K distribution is the only distribution where a safe theoretical prediction can be made. The uncorrected cosθ K distribution of Kπ events before any mass cuts, i.e. mostly background, shown in Fig. 2 is strongly peaked towards cosθ K =−1. For signal events, the cosθ K distribution should be flat after acceptance correction since the D 0 has spin 0. This expectation for the D 0 is confirmed in Fig. 2 where, for the tagged D 0 sample, the distribution for D 0 is shown before and after acceptance correction (the method of signal extraction and the correction for the acceptance will be described in Sects. 4 and 6).

Fig. 2
figure 2

Distribution of cosθ K in the Kπ rest frame for (mostly) background combinations (scaled by 0.001, solid line), for the D 0 signal region (scaled by 5, full circles) and for the acceptance corrected D 0 signal (open squares). The dashed lines correspond to the |cosθ K |<0.9 cut. The D 0 signal is from the 2004 tagged D 0 sample

Full size image

The so-called ambiguity cut applied in Ref. [9] is not applied in the present analysis. This cut discards an event if two D 0 or \(\bar{D}^{0}\) meson candidates are found within the mass window of ±700 MeV and removes a significant number of good events. However, the probability to find two D 0 (or two \(\bar{D}^{0}\)) mesons in the signal peak is practically zero. Hence the present analysis, which extracts the number of signal events from fitting separately the D 0 and \(\bar{D}^{0}\) peaks, does not suffer from this ambiguity.

In the mass window of ±700 MeV around the nominal D 0 mass, the tagged D 0 sample consists of 160×103 neutral Kπ combinations. In order to avoid overlapping samples, at this stage, the D ∗ candidates are taken out of the untagged D 0 sample. In the same mass window, the untagged D 0 sample comprises 17×106 neutral Kπ combinations.

The invariant Kπ mass spectra are shown in Fig. 3a for the untagged D 0 sample, for all neutral Kπ combinations and also separately for the K + π − and K − π + combinations. These spectra exhibit the D 0 peak at 1865 MeV. The prominent peak to the left is due to the decay of the narrow \(K_{2}^{*}(1430)\). In Fig. 3c invariant mass spectra are shown for the tagged D 0 sample. In this case, only some feed-through of the \(K_{2}^{*}(1430)\) resonance is seen and a pronounced, rather narrow peak about 250 MeV below the nominal D 0 mass. As shown by Monte Carlo simulations, this peak at 1620 MeV results from 3-body decays of the D 0→Kππ 0, where the π 0 escaped detection, with some contributions from D ∗ decays with more than 3 particles in the final state. The signal-over-background (S/B) ratio is about 1:1 for the events of the tagged D 0 sample. For the untagged D 0 sample, S/B is only 1:10, but the number of signal events is four times higher.Footnote 1

Fig. 3
figure 3

Invariant M(Kπ) mass spectra within a window of ±700 MeV around the nominal D 0 mass. (a) D 0 sample before and (b) after background subtraction, (c) D ∗ sample before and (d) after background subtraction. Both neutral charge combinations are shown separately, together with their sum in (a) and (c). See text for the background subtraction by fits

Full size image

Mass spectra of all the Kπ combinations are shown in Fig. 4a separately, using only data from 2004. The spectra for the two neutral charge combinations show three narrow peaks corresponding to K ∗(890), \(K_{2}^{*}(1430)\) and D 0(1865). Also, other short lived kaonic (strange) resonances are present but they superimpose together with combinatorial background to a structureless distribution that can almost perfectly be described by a single exponential function, see Fig. 4b.

Fig. 4
figure 4

Kπ invariant mass distribution for the full mass range using data from 2004 only. (a) All four charge combinations (from top to bottom + −, − +, − −, + +) are shown separately, with a linear vertical scale; (b) only neutral, non-exotic charge combinations (− +, + −) summed up, logarithmic vertical scale

Full size image

4 Method of signal extraction

The invariant Kπ mass spectra shown in Figs. 3a, c are fitted with a function given by the sum of the following elements: a Gaussian for the D 0→Kπ signal, an exponential for the background, a shape determined by Monte Carlo simulations for the peak at 1620 MeV from 3-body decays of the D 0 and by relativistic Breit-Wigner intensities for the \(K_{2}^{*} (1430)\) and \(K_{3}^{*} (1780)\). The latter K resonance is barely visible in this spectrum but shows up clearly in certain kinematic regions, see below.

The Kπ spectra are remarkably well described fitting them with 12 parameters, as described above. Figures 3b, d show the spectra after subtraction of the exponential background. From the fits one obtains (3610±90) \(D^{*+}\rightarrow(K^{-} \pi^{+}) \pi^{+}_{s}\) and (4530±100) \(D^{*-}\rightarrow(K^{+} \pi^{-}) \pi^{-}_{s}\) for the tagged sample as well as (15200±800) D 0→K − π + and (18400±900) \(\bar {D}^{0}\rightarrow K^{+} \pi^{-}\) for the untagged D 0 sample.

The dependence on kinematic variables of the production rate of D 0 and D ∗, together with those of the neighbouring \(K_{2}^{*} (1430)\) resonance and the background, is extracted by fitting the mass spectra for each kinematic bin. Alternatively, the signal distributions of the D 0 and the D ∗ are obtained by side-band subtraction.

Using the first method, the fit yields in every bin of a given kinematic variable the number of D, \(K_{2}^{*} (1430)\) and \(K_{3}^{*} (1780)\) together with the background. In Figs. 5a–f examples of the Kπ invariant mass spectra for different intervals in z are shown before and after the subtraction of the fitted exponential background. These fits did not include the \(K_{3}^{*} (1780)\). The fitting method allows monitoring of all details of the fit, as illustrated in the inserts of Fig. 5. The broad structure showing up for z>0.75 in Fig. 5 is attributed to the \(K_{3}^{*} (1780)\) resonance. This resonance follows the same behaviour as the \(K_{2}^{*} (1430)\) resonance, i.e. it is produced at larger values of z than the D 0 (see Sect. 5). The introduction of the \(K_{3}^{*} (1780)\) resonance in the final fit also removed a small but statistically significant and unexplained discrepancy between fit and data on the left side of the D 0→Kπ peak in the z-integrated spectrum, where the fit before the inclusion of the \(K_{3}^{*} (1780)\) was systematically below the data (see Fig. 5g to be compared with Fig. 3b, where the \(K_{3}^{*} (1780)\) has been included).

Fig. 5
figure 5

(a–f) Invariant Kπ mass spectra in bins of the energy fraction z for the untagged D 0 sample. The vertical (red) lines indicate the nominal positions of \(K_{2}^{*} (1430)\) and D 0. The inserts demonstrate the signal behaviour after the removal of the fitted exponential background. The fit contains D 0→Kπ at 1865 MeV, D 0→Kππ 0 at 1620 MeV, the \(K_{2}^{*} (1430)\) and an exponential background. (g) shows the signal behaviour after removal of the fitted exponential background for the entire z range (no \(K_{3}^{*} (1780)\) assumed). This figure has to be compared with Fig. 3b where the \(K_{3}^{*} (1780)\) was included in the fit (Color figure online)

Full size image

The second method for signal extraction is the standard side-band subtraction. This method can only be applied to the D 0 and the D ∗ signals, due to the limited mass range (±700 MeV around the nominal D 0 mass). Three Kπ mass windows are chosen. The central one, which is 100 MeV wide and centered at the nominal D 0 mass, contains the D 0→Kπ signal plus background. The two side-bands contain only background events. They are 50 MeV wide and centered at ±100 MeV above or below the nominal D 0 mass. Thus three independent distributions are obtained as a function of each kinematic variable. The sum of the side-band distributions is subtracted from the central distribution, assuming that the side-band distributions correctly represent the distribution of background under the signal. This assumption is supported by the observed similar behaviour of the distributions in the two side-bands.

Usually, the background below the signal is obtained by linear interpolation between the side-bands. Such a linear interpolation overestimates the background under the signal. Therefore it cannot be applied for the untagged D 0 sample, where S/B∼1/10. Instead, an estimate of the background under the signal is obtained from the fit. The total number of background events in the two side-bands is correspondingly rescaled.

For the chosen width of the central window, about 5 % of the signal is found outside. Hence the number of signal events obtained by side-band subtraction is expected to be lower by 5 % than that obtained with the signal fitting method.

5 Comparison of kinematic distributions

In this section, event distributions are shown as a function of the relevant kinematic parameters, for both the tagged and untagged D 0 samples as well as for the \(K_{2}^{*} (1430)\) and background. The data collected in 2002–2006 are used, and the distributions are not corrected for acceptance. However, the geometric acceptances for the various compared Kπ systems are similar.

The distributions of the \(K_{2}^{*} (1430)\) signal are obtained from the untagged D 0 sample using the signal fitting method. The distributions of the Kπ-background combinations are extracted from the two Kπ side bands of the tagged D 0 sample, at invariant masses of 1765±50 and 1965±50 MeV. The kinematic distributions of D 0 and D ∗ are obtained by applying both signal extraction methods described above, allowing to cross-check the stability of the result. While for the tagged D 0 sample perfect agreement is found between the two methods, for the untagged sample some disagreement beyond the statistical error is observed, for instance at low values of z or low Kπ energy E. This is the result of strongly varying background shapes with additional broad resonances emerging below the \(K_{2}^{*} (1430)\). The corresponding data points for D 0 and \(K_{2}^{*} (1430)\) are omitted, since a more complex background description would be needed.

In Fig. 6, the distributions of the D 0, the \(K_{2}^{*} (1430)\) and the background under the D 0 are compared, showing their different behaviour. The distributions as a function of the inclusive variables Q 2 and x Bj are displayed in Figs. 6a, b. For the tagged sample, the average values of Q 2 and x Bj extracted from these distributions are about 0.5 (GeV/c)2 and 0.005, respectively. Some differences between signal and background events are observed at large values of Q 2 and x Bj . As a function of ν, the distributions for the various Kπ systems are significantly different (see Fig. 6c). The \(K_{2}^{*} (1430)\) distribution peaks at lower values than that of D 0, and the rise at low ν that is caused by the increase of both acceptance and cross-section starts at lower ν. The background peaks at a somewhat higher values, but has a similar rise with ν as D 0 and D ∗.

Fig. 6
figure 6

Measured kinematic distributions of various (Kπ) systems before acceptance correction as a function of (a) Q 2, (b) x Bj , (c) ν, (d) E, (e) \(p_{T}^{2}\), and (f) z. The symbols D 0 and \(K^{*}_{2}\) denote D 0 and \(K_{2}^{*} (1430)\) from the untagged sample. The symbols D ∗ and bkg denote D ∗ and background from the two side-band windows for the D ∗ sample. The data are from the years 2002 to 2006

Full size image

No clear differences are observed between the distributions as a function of the Kπ energy E (see Fig. 6d). Given the reason described above, the \(K_{2}^{*} (1430)\) and the D 0 data points at lower values of E are omitted from the untagged sample. The distribution for the D ∗ signal as a function of \(p_{T}^{2}\) (Fig. 6e) shows an almost single-exponential decrease, while the distribution for the D 0 flattens above 3 (GeV/c)2. The difference between D 0 and D ∗ may be related to the fact that for the D ∗ only the \(p_{T}^{2}\) of the 2-body subsystem is shown. Both distributions are significantly different from those of background and \(K_{2}^{*} (1430)\). From a fit of an exponential function up to \(p_{T}^{2} = 2~(\mbox{GeV/c})^{2}\), the following slopes are obtained in units of (GeV/c)−2: −0.84±0.03 for D ∗, −0.96±0.06 for D 0, −1.94±0.01 for \(K_{2}^{*} (1430)\) and −1.69±0.01 for background. The distributions in z show significant differences, too. The background is concentrated at smaller values of z than the D 0 signal. Moreover, the distribution of the \(K_{2}^{*} (1430)\) is peaked at significantly higher values of z than that of the D 0.

In conclusion of the comparison: remarkable differences are observed between the distributions of the D meson signals, the \(K_{2}^{*} (1430)\), and the background as a function of the kinematic variables ν, p T and z. This clearly points to different production mechanisms for D mesons and the \(K_{2}^{*} (1430)\). The observed differences between D mesons and \(K_{2}^{*} (1430)\) agree qualitatively with the differences expected if the D mesons result from the fragmentation of a pair of charm quarks and the \(K_{2}^{*} (1430)\) from the fragmentation of a quark knocked out in a leading order process.

The interpretation of the kinematic distributions of the background is more complex, since this background is dominated by combinatorial entries. No attempt is made to interpret it. However, one should mention that other background events of non-combinatorial origin (e.g. in the untagged sample the background taken from side bands has also large contributions from resonances or from πK correlated production in the fragmentation) have been observed to behave very similar to the background shown in Fig. 6.

6 Acceptance and integrated luminosity

Acceptances and integrated luminosity, which are needed to extract semi-inclusive total and differential cross-sections, are calculated only for the tagged D 0 sample of the year 2004. Since this is the first detailed acceptance calculation for this particular final state at COMPASS, the present section also aims at illustrating the acceptances of the COMPASS spectrometer for the detection of the scattered muon and the D ∗±. For this reason, 2-dimensional acceptances will be shown as a function of selected variables.

Acceptance calculations are done using a complete Monte Carlo simulation of the detector configuration, including the triggers and the track reconstruction code for the 2004 data. Events are generated using AROMA 2.2.4 [10], which assumes photon–gluon fusion into \(c \bar{c}\) to be the dominant underlying mechanism for D ∗ production. Default fragmentation functions are used and parton showers are generated. The charm quark mass is set to 1.35 GeV. Produced D ∗s are forced to decay to D 0 π +→K − π + π + for D ∗+ or to \(\bar{D}^{0}\pi^{-} \rightarrow K^{+} \pi^{-} \pi^{-}\) for D ∗−. Trigger conditions and data selection criteria applied to the Monte Carlo events are the same as for real data. In total, 107 events were generated for both decays. The acceptances are calculated as a function of the reconstructed values of the kinematic variables, thus accounting for experimental resolution and bin-to-bin smearing.

Figure 7 shows the number of generated events (a) as a function of x Bj , y and (b) as a function of p ⊥, E of the D meson. The transverse momentum p ⊥ is measured with respect to the direction of the incoming muon beam. In both pictures the generated events are mainly concentrated in the lower left corner.

Fig. 7
figure 7

Number of generated (AROMA) events as a function of (a) x Bj and y and (b) p ⊥ and E

Full size image

For illustration, the acceptance for D ∗ production is shown at two stages, i.e. after requiring the reconstruction of the scattered muon and after the additional reconstruction of the three hadrons from the D ∗ decay. The ‘inclusive’ acceptance A μ (x Bj ,y) is shown in Fig. 8a, and the overall acceptance \(A_{D^{*}}(x_{Bj},y)\) in Fig. 8b. In the kinematic region relevant for charm production, the inclusive acceptance A μ (x Bj ,y) is fairly homogeneous and ranges between 50 % and 80 %. The overall acceptance \(A_{D^{*} }(x_{Bj},y)\) is also homogeneous for y>0.2 and ranges from 1 % to 5 %. The cut-off at y=0.2 is due to the momentum selection for the RICH identification.

Fig. 8
figure 8

(a) ‘Inclusive’ acceptance A μ (x Bj ,y) and (b) overall acceptance \(A_{D^{*}}(x_{Bj},y)\) before applying the E cut

Full size image

The overall acceptance \(A_{D^{*}}\) as a function of E and p ⊥ (i.e. transverse momentum with respect the incoming muon) is shown in Fig. 9. The upper limit of about 100 mrad for the spectrometer acceptance in the year 2004 can be seen at low energy and large p ⊥. For 20 GeV<E<80 GeV the acceptance ranges between 5 % and 13 %. Outside this energy region the acceptance drops to zero due to the lack of particle identification and therefore 20 GeV<E<80 GeV is required in the further analysis. The one-dimensional acceptances used below to determine the differential inclusive cross-sections are limited to this range of D 0 energies.

Fig. 9
figure 9

Overall acceptance \(A_{D^{*}}( p_{\perp},E)\)

Full size image

The one-dimensional acceptance functions \(A_{D^{*+}}\) and \(A_{D^{*-}}\) are shown in Fig. 10 as a function of ν, E, z and \(p_{T}^{2}\). In addition, the ratio of the D ∗+ acceptance over that of D ∗− is shown in each case. Note that within the statistical uncertainty of the Monte Carlo simulation the acceptances are almost equal for D ∗+ and D ∗−, with the acceptance for D ∗+ being slightly higher than that for D ∗−.

Fig. 10
figure 10

One-dimensional acceptances for D ∗± production as a function of (a) ν, (b) E, (c) z and (d) \(p_{T}^{2}\). Boxes (red) correspond to D ∗+, triangles (blue) to D ∗− events. The circles (black) show the ratio of acceptances \(A_{D^{*+}} / A_{D^{*-}}\), the ordinates for the ratios are drawn on the right-hand side of the figures (Color figure online)

Full size image

The integrated luminosity \(\mathcal{L}\) is determined by a comparison of the measured number of inclusive inelastic muon scattering events with the best available measurement of the corresponding cross-section. The differential number of events is the product of integrated luminosity, inclusive muon acceptance and inclusive differential cross-section:

$$ \frac{d^2 N(x_{Bj},y)}{dx_{Bj} \, dy} = { \mathcal{L}} \cdot A_{\mu}(x_{Bj},y) \cdot\frac{d^2\sigma^{\mu N \rightarrow \mu'X}(x_{Bj},y)}{dx_{Bj} \, dy}. $$
(1)

The inclusive inelastic muon-deuteron cross-section was measured by the NMC Collaboration for various muon energies between 90 and 280 GeV and published as a parameterization of the structure function F 2 [26]. Thus the cross-section has to be reconstructed based on this F 2 parameterization. The measured cross-section is connected with the one-photon exchange cross-section via a radiative correction factor η(x Bj ,y):

$$ \frac{d^2\sigma^{\mu N \rightarrow\mu'X}(x_{Bj},y)}{dx_{Bj} \, dy} =\frac{1}{\eta(x_{Bj},y)}\frac{d^2\sigma_{1 \gamma}(x_{Bj},y)}{dx_{Bj} \, dy}. $$
(2)

The one-photon exchange cross-section is connected with F 2 by:

(3)

where m is the muon mass. The factor R(x Bj ,Q 2) is the cross-section ratio for longitudinal over transverse photons:

$$ R\bigl(x_{Bj},Q^2\bigr) = \frac{\sigma_L}{\sigma_T}. $$
(4)

The radiative correction factor η(x Bj ,y) is calculated with codes based on [27]. The ratios R(x Bj ,Q 2) are determined as in Ref. [28]. Given the light material composing the target (Li, D and He), nuclear effects have been neglected.

The integrated luminosity is determined in bins of (x Bj ,y) as:

$$ {\mathcal{L}} = \frac{1}{A_{\mu}(x_{Bj},y)}\cdot \frac{d^2 N(x_{Bj},y)/(dx_{Bj} \, dy)}{ d^2\sigma^{\mu N \rightarrow\mu'X}(x_{Bj},y)/(dx_{Bj} \, dy)} . $$
(5)

The integrated luminosity on the left-hand side of Eq. (5) has to be constant, while all terms on the right-hand side depend on x Bj and y. As a side product of extracting the integrated luminosity, this equation can be used to evaluate the uncertainty of the muon acceptance calculation for Q 2 values larger than about 0.6 (GeV/c)2, where the NMC parameterization is valid. The values of \({\mathcal{L}}\) obtained for different (x Bj ,y) bins vary indeed by up to 20 % over the relevant (x Bj ,y) range, so that an overall systematic uncertainty of 20 % is attributed to the product of integrated luminosity and inclusive muon acceptance. The average value of the integrated luminosity is calculated as a weighted mean of the luminosities determined in (x Bj ,y) bins, using the data at Q 2>0.6 (GeV/c)2. For a given bin the weight is the number of events in that bin. The result for the integrated luminosity of the 2004 data is 0.71±0.14 fb−1. Since the statistical uncertainty is negligible, only the systematic one is quoted.

7 D ∗± production cross-sections

The acceptance uncorrected distributions presented in Sect. 5 were given for all data taken in 2002–2006. The signals for D 0 and \(\bar{D}^{0}\) were summed up, and so were those for D ∗+ and D ∗−. In the following, the semi-inclusive differential cross-sections for D ∗± production, determined for data from the year 2004 only, will be obtained separately for D ∗+ and D ∗−. The acceptances, the integrated luminosity and the known branching ratio (2.6 %) of D ∗ to Kππ are taken into account. At the end of this section, D ∗+ and D ∗− asymmetries will be shown for all 2002 to 2006 data, since integrated luminosity and also the acceptances cancel in these asymmetries to a good approximation.

Figure 11 displays the semi-inclusive differential cross-sections of D ∗+ and D ∗− events as a function of ν,E,z and \(p_{T}^{2}\). The numerical values of the measured differential cross-sections are compiled in Table 1. These cross-sections are compared with the theoretical predictions obtained from the AROMA generator, which assumes \(c\bar{c}\) production via photon gluon fusion, uses CTEQ2L as default LEPTO [29] PDFs and includes parton showers. The AROMA total cross-section is rescaled to the value of 1.9 nb measured by COMPASS, see below. They were calculated using the same program package parameters as those for the determination of acceptances. The \(p_{T}^{2}\) and z distributions are also compared with results published by the EMC Collaboration 20 years ago [21], based on 92 events, obtained with higher muon beam energy and a cut on Q 2>3 (GeV/c)2. EMC combined \(\bar{D}^{0}\) and D 0 as within the statistical precision no differences were observed. In order to compare with the present data, their measured values and uncertainties are divided by a factor of 2.

Fig. 11
figure 11

Semi-inclusive differential cross-sections for D ∗+ and D ∗− production as a function of (a) virtual photon energy ν, (b) D 0 energy E, (c) fractional energy z and (d) squared transverse momentum \(p_{T}^{2}\). For all distributions, the squares (red) correspond to D ∗+ and triangles (blue) to D ∗− events (2004 data, D ∗ sample). The circles (green) are semi-inclusive differential cross-sections for D 0 from the EMC experiment, see text. The curves represent AROMA predictions, dashed for D ∗− and dotted for D ∗+ (Color figure online)

Full size image
Table 1 Semi-inclusive differential cross-sections for D ∗+ and D ∗− production as a function of (a) γ ∗ energy ν, (b) D 0 energy E, (c) fractional energy z and (d) squared transverse momentum \(p_{T}^{2}\) of the D 0. The central values and bin sizes of ν and E are given in units of GeV, those of \(p_{T}^{2} \) in (GeV/c)2. The last two lines show the integrated cross-sections. Statistical uncertainties are given
Full size table

Good agreement is observed between the shapes of the measured distributions and the corresponding AROMA predictions. The distributions of D ∗+ and D ∗− as a function of ν show that the points for D ∗− are systematically higher than those for D ∗+. The effective threshold of D ∗+ appears to be about 10 GeV higher than that of D ∗−. The AROMA generator produces also somewhat more D ∗− than D ∗+ but the differences at threshold are far less pronounced. A similar feature can be observed for the z distribution. In the large-z region, that has a large contribution from low-ν events, the cross-section of D ∗− becomes significantly larger than that of D ∗+. The AROMA calculations predicts more D ∗− than D ∗+ as well, but the size of the effect is smaller. For the semi-inclusive differential cross-sections as a function of E and \(p_{T}^{2}\), no remarkable differences are observed between the shapes of the distributions of D ∗+ and D ∗−.

The total cross-sections for D ∗+, D ∗− and D ∗± production are extracted by integration of the differential ones. The differences between the results from the integration over ν, E, z and \(p_{T}^{2}\) (see Table 1) are used to evaluate the systematic uncertainty of acceptance corrections. Using the RMS of the four results (from ν, E, z and \(p_{T}^{2}\)) one obtains a systematic contribution of 0.05 for both D ∗+ and D ∗− and 0.10 for the sum D ∗±, i.e. at the level of the statistical uncertainty. In the ratio of D ∗+ over D ∗− the acceptances almost cancels. The values of the ratio vary between 0.77 and 0.81, with an average of 0.80 and a RMS of <0.02, i.e. two to three times smaller than the statistical uncertainty of ∼0.05.

The final result for the D ∗ meson production cross-section is then σ(μN→μ′D ∗± X)=1.86±0.06 (stat)±0.10 (sys)±0.37 (luminosity) nb. The only cut applied is the energy window for the D 0 meson between 20 GeV<E<80 GeV in the laboratory frame, corresponding to 22 GeV<E<86 GeV for the D ∗ energy.

For charm-anticharm production, AROMA gives a cross-section of 7.2 nb with 1.35 GeV chosen as the default charm quark mass. Using the common assumption of 0.6 D ∗ mesons per charm event and accounting for the energy cut 20 GeV<E<80 GeV, which reduces the number of charm Monte Carlo events by another factor of 0.6, the corresponding AROMA cross-section predicted for COMPASS is 2.6 nb. Given the number of assumptions which underlie the AROMA default options (charm quark mass, fragmentation, no radiative corrections, leading order QCD apart from parton showers) the agreement with the above experimental result is considered to be good.

However, deviations from the AROMA predictions are observed in the data with respect to D ∗+ and D ∗− production. These may provide valuable insight into their production mechanisms. In a simple LO approach, assuming photon–gluon fusion with independent fragmentation of the charm and anti-charm quarks to be the relevant production mechanism, no differences should be observed between D ∗+ and D ∗−. Differences may occur for all processes where the quark content of the target nucleon matters. The quark content of D ∗ mesons indicate that only D ∗− may contain a valence quark from the target nucleon. Furthermore, instead of fragmenting into D ∗+ the c quark, together with a diquark of the target nucleon, can hadronize into a charmed baryon, leading to associated production of e.g. D ∗− Λ c . Thus the D ∗− may result from a valence quark and/or associated production. If parton showers are included in AROMA the flavour dependent quark distribution functions of the nucleon come into play. Processes like associated production of D ∗− Λ c lead to differences between kinematic distributions of D ∗+ and D ∗−. The same happens for processes where an initial quark in the nucleon absorbs the virtual photon and radiates a heavy gluon which then decays to \(c \bar{c}\), or where in the course of fragmentation the \(\bar{c}\) quark picks up quarks from the nucleon.

In order to provide statistically more precise information on the potentially interesting differences between D ∗+ and D ∗− production, Fig. 12 shows particle-antiparticle asymmetries of the semi-inclusive cross-sections,

$$ A(X) =\frac{d\sigma^{D^{*+}}(X) -d\sigma^{D^{*-}}(X)}{d\sigma^{D^{*+} }(X) + d \sigma^{D^{*-}}(X)} $$
(6)

as a function of X=ν, E, z and \(p_{T}^{2} \) for both the D ∗ sample and Monte Carlo events generated by AROMA. Here the full statistics of the years 2002–2006 is used. It is assumed that the acceptances for the two charge combinations are equal. In the previous section it was shown that for the year 2004 this is indeed approximately true. The numerical values of the measured asymmetries are given in Table 2, where only statistical uncertainties are shown, based on the assumption that acceptance cancels. A small cross-section assymmetry between D + and D − production has been observed recently in a different energy range by the LHCb experiment [30].

Fig. 12
figure 12

Measured D ∗+ and D ∗− asymmetries for data (stars) and AROMA generator (crosses) events as a function of X=ν, E, z and \(p_{T}^{2}\). All 2002–2006 data are used (Color figure online)

Full size image
Table 2 Measured asymmetry A(X) as a function of X=ν, E, z and \(p_{T}^{2} \). The central values and bin sizes of ν and E in (a) and (b) are given in units of GeV, those of \(p_{T}^{2} \) in (d) in units of (GeV/c)2
Full size table

As one can see from the figure, the measured asymmetry decreases significantly stronger than that predicted by AROMA when ν decreases below 40 GeV and/or when z increases above 0.6. The distributions shown as a function of ν clearly exhibit different thresholds for D ∗+ and D ∗− production, which supports a stronger presence of mechanisms other than PGF with independent fragmentation. As a function of z, the most pronounced differences between D ∗+ and D ∗− are seen at large values of z, whereas at z values lower than 0.5 the production rates are nearly equal. Values of z larger than 0.5 indicate an asymmetric sharing of the energies between a D meson and its associated partner with opposite charm content. Since the cross-section of D ∗−, which contains a down and an anti-charm quark, increases with increasing z stronger than that of D ∗+, this observation suggests processes where the anti-charm quark is fast and the charm-quark is slow. Here, a candidate process is again associated production of a D ∗− along with a charmed baryon, i.e. D ∗− Λ c . Alternatively, since the D ∗− may also contain a valence quark of the nucleon whereas the D ∗+ does not, one may think of processes other than associated production, which involve valence quarks of the nucleon.

Asymmetries between the production of D 0 and \(\bar{D}^{0}\) or D ∗+ and D ∗− were already observed in numerous earlier experiments (see e.g. [22] for charm photoproduction and [31–37] for charm production by hadrons), although not as pronounced as in the present experiment that covers the region of virtual-photon energies from threshold up to 140 GeV.

8 Summary and conclusions

The observed total cross section of (1.9±0.4) nb for the production of D ∗+ and D ∗− mesons in inelastic muon nucleon interactions at 160 GeV incident muon energy within the COMPASS acceptance (20 GeV<E<80 GeV and 22 GeV<E<86 GeV, for D 0 and D ∗ respectively) lies within the range of values expected if the dominant process is photon–gluon fusion to open charm production. The total error is dominated by the uncertainty on the luminosity.

The detailed comparison of the measured differential cross sections of D ∗ production as a function of the variables ν, E, z and \(p_{T}^{2}\) shows good agreement with those expected from the model underlying the AROMA generator used to produce the theoretical distributions. This is remarkable as most of the kinematic distributions of D mesons are quite different in shape compared to those of the background and the neighbouring \(K_{2}^{*} (1430)\) resonance.

The observed large asymmetries between D ∗+ and D ∗− production for ν<40 GeV and z>0.6 can only partially be described by the model used in AROMA, which predicts differences of the same sign but of smaller magnitude. This indicates that the hadronization processes of charm and anti-charm quarks differ more significantly than expected or/and processes other than PGF contribute by a larger amount than assumed.

The observed dependences of these differences on the kinematic variables, in particular on the photon energy ν and the fractional energy z, suggest that associated production (e.g. DΛ c ) plays a dominant role at low photon energy. Also, D ∗− production involving valence quarks of the nucleon may contribute to the observed asymmetries between D ∗+ and D ∗− production.