Production of excited charm and charm-strange mesons at HERA

The production of excited charm, D_1(2420)^0 and D_2^*(2460)^0, and charm-strange, D_{s1}(2536)^+-, mesons in ep collisions was measured with the ZEUS detector at HERA using an integrated luminosity of 126 pb^-1. Masses, widths and helicity parameters were determined. The measured yields were converted to the rates of c quarks hadronising as a given excited charm meson and to the ratios of the dominant D_2^*(2460)^0 and D_{s1}(2536)^+- branching fractions. A search for the radially excited charm meson, D^{*'}(2640)^+-, was also performed. The results are compared with those measured previously and with theoretical expectations.


Introduction
Heavy-quark spectroscopy has recently undergone a renaissance with the discovery of several new states [1]. The properties of these states challenge the theoretical description of heavy-quark resonances. Therefore, further measurements of excited charm and charmstrange mesons are important.
The lowest-mass states of the cq (cq) system (q = u, d, s) with spin zero (D mesons) and spin one (D * mesons) and with orbital angular momentum L = 0 are well established [1]. A singlet and a triplet of states with L = 1 are expected. These P -wave (L = 1) mesons can decay to charm mesons with L = 0 by emitting a pion or a kaon. Heavy Quark Effective Theory [2] (HQET) predicts that, in the heavy-quark limit (m Q →∞), the properties of the P -wave mesons are determined mainly by the total angular momentum of the light quark, j = L + s, where s denotes the spin of the light quark. Consequently, the four states are grouped in two doublets with j = 3/2 or 1/2. Only D-wave decays are allowed for the members of the j = 3/2 doublet; therefore they are supposed to be narrow. On the other hand, the members of the j = 1/2 doublet decay through S-wave only and therefore are expected to be broader [3]. Due to the finite charm quark mass a separation of the two doublets is only an approximation and amplitudes of two observable states with J P = 1 + can be mixtures of D-and S-wave amplitudes. Here J and P are the total angular momentum and parity of the cq system. Two pairs (neutral and charged) of narrow non-strange excited charm mesons, D 1 (2420) 0,± and D * 2 (2460) 0,± , and a pair of narrow charm-strange excited mesons, D s1 (2536) ± and D s2 (2573) ± , were observed and tentatively identified as the members of the j = 3/2 doublets with J P = 1 + and 2 + , respectively [1]. Recently, the HQET expectations were supported by the first measurements of the broad non-strange excited charm mesons: neutral and charged D * 0 (2400) 0,± with J P = 0 + [4,5], and D 1 (2430) 0 with J P = 1 + [4]. The predicted broad non-strange charged excited charm meson with J P = 1 + has not yet been observed. The recent discovery of two additional charm-strange excited mesons, D * s0 (2317) ± with J P = 0 + and D s1 (2460) ± with J P = 1 + reported initially by BABAR [6] and CLEO [7], respectively, revealed their surprisingly small masses and narrow widths [1]. The small mass values forbid their decay into D ( * ) K final states.
In addition to the orbital excitations, radially excited charm mesons D ′ (J P = 0 − ) and D * ′ (J P = 1 − ) were predicted with masses of about 2.6 GeV and dominant decay modes to Dππ and D * ππ, respectively [8,9]. An observation of a narrow resonance in the final state D * ± π + π − at 2637 MeV was reported and interpreted as the radially excited D * ′± meson by DELPHI [10]. However, OPAL found no evidence for this narrow resonance in an analogous search [11].
HERA [12,13]. The large charm production cross section at HERA also provides a means to study excited charm and charm-strange mesons produced in ep collisions. The first such study is reported in this paper. It is restricted to decays, for which significant signals were identified: The corresponding antiparticle decays were also measured 1 . A search for the radially excited charm meson, D * ′ (2640) + , in the D * + π + π − final state was also performed.

Experimental set-up
The analysis was performed using data taken with the ZEUS detector from 1995 to 2000. In this period, HERA collided electrons or positrons 2 with energy E e = 27.5 GeV and protons with energy E p = 820 GeV (1995)(1996)(1997) or E p = 920 GeV (1998)(1999)(2000). The data used in this analysis correspond to an integrated luminosity of 126.5 ± 2.4 pb −1 .
A detailed description of the ZEUS detector can be found elsewhere [14]. A brief outline of the components most relevant to this analysis is given below.
Charged particles were tracked in the central tracking detector (CTD) [15], which operated in a magnetic field of 1.43 T provided by a thin superconducting solenoid. The CTD consisted of 72 cylindrical drift chamber layers, organized in nine superlayers covering the polar-angle 3 region 15 • < θ < 164 • . The transverse-momentum resolution for full-length tracks was σ(p T )/p T = 0.0058p T ⊕ 0.0065 ⊕ 0.0014/p T , with p T in GeV. To estimate the energy loss per unit length, dE/dx, of charged particles in the CTD [16,17], the truncated mean of the anode-wire pulse heights was calculated, which removes the lowest 10% and at least the highest 30% depending on the number of saturated hits. The measured dE/dx values were corrected for a number of effects [18] and normalised such that the corrected value was one for a minimum ionising particle. The resolution of the dE/dx measurement for full-length tracks was about 9%.
The high-resolution uranium-scintillator calorimeter (CAL) [19] consisted of three parts: the forward (FCAL), the barrel (BCAL) and the rear (RCAL) calorimeters. Each part was 1 Hereafter, charge conjugation is implied. 2 From now on, the word "electron" is used as a generic term for electrons and positrons. 3 The ZEUS coordinate system is a right-handed Cartesian system, with the Z axis pointing in the proton beam direction, referred to as the "forward direction", and the X axis pointing left towards the centre of HERA. The coordinate origin is at the nominal interaction point.
subdivided transversely into towers and longitudinally into one electromagnetic section (EMC) and either one (in RCAL) or two (in BCAL and FCAL) hadronic sections (HAC). The smallest subdivision of the calorimeter was called a cell. The CAL energy resolutions, as measured under test-beam conditions, were σ(E)/E = 0.18/ √ E for electrons and σ(E)/E = 0.35/ √ E for hadrons, with E in GeV.
The luminosity was determined from the rate of the bremsstrahlung process ep → eγp, where the photon was measured with a lead-scintillator calorimeter [20] located at Z = −107 m.

Event simulation
Monte Carlo (MC) samples of charm and beauty events were produced with the Pythia 6.156 [21] and Rapgap 2.0818 [22] event generators. The Rapgap MC used Heracles 4.6.1 [23] in order to incorporate first-order electroweak corrections. The generation included direct photon processes, in which the photon couples directly to a parton in the proton, and resolved photon processes, where the photon acts as a source of partons, one of which participates in the hard scattering process. The CTEQ5L [24] and GRV LO [25] parametrisations were used for the proton and photon structure functions, respectively. The charm and bottom quark masses were set to 1.5 GeV and 4.75 GeV, respectively. Events for all processes were generated in proportion to the MC cross sections. The Lund string model [26] as implemented in Jetset [21] was used for hadronisation in Pythia and Rapgap. The Bowler modification [27] of the Lund symmetric fragmentation function [28] was used for the charm and bottom quark fragmentation. To generate D * ′+ mesons, the mass of a charged charm meson in the Jetset particle table was set to 2.637 GeV, its width was set to 15 MeV and the decay channel was set to D * + π + π − [10].
The Pythia and Rapgap generators were tuned to describe the photoproduction and the deep inelastic scattering (DIS) regimes, respectively. Consequently, the Pythia events, generated with Q 2 < 0.6 GeV 2 , were combined with the Rapgap events, generated with Q 2 > 0.6 GeV 2 , where Q 2 is the exchanged-photon virtuality. Diffractive events, characterised by a large rapidity gap between the proton at high rapidities and the centrallyproduced hadronic system, were generated using the Rapgap generator in the diffractive mode and combined with the non-diffractive MC sample. The contribution of diffractive events was estimated by fitting the η max distribution 4 of the data with a linear combination of the non-diffractive and diffractive MC samples.
To ensure a good description of the data, the transverse momenta, p T (D * + , D + , D 0 ), and pseudorapidity, η(D * + , D + , D 0 ), distributions were reweighted to the data for the combined Pythia+Rapgap MC sample. The reweighting factors, tuned using a large D * + sample (Section 4), were used for D + and D 0 mesons relying on the MC description of the differences between the D * + and D + or D 0 distributions. The effect of the reweighting on the measured values was small; the reweighting uncertainty was included when evaluating systematic uncertainties (Section 8).
The generated events were passed through a full simulation of the detector using Geant 3.13 [29] and processed with the same reconstruction program as used for the data.

Event selection and reconstruction of lowest-mass charm mesons
Events from both photoproduction [30] and DIS [13] were selected online with a threelevel trigger [14,31]. The first-and second-level trigger used CAL and CTD data to select ep collisions and to reject beam-gas events. At the third level, where the full event information was available, the nominal charm-meson trigger branches required the presence of a reconstructed D * + , D + or D 0 candidate. The efficiency of the online charmmeson reconstruction, determined relative to the efficiency of the offline reconstruction, was above 95%. Events missed by the nominal charm-meson triggers but selected with any other trigger branch, dominantly from an inclusive DIS trigger and a photoproduction dijet trigger, were also used in this analysis.
In the offline analysis, only events with |Z vtx | < 50 cm, where Z vtx is the primary vertex position determined from the CTD tracks, were used. The D * + , D + and D 0 mesons were reconstructed using tracks measured in the CTD and assigned to the reconstructed primary event vertex. To ensure both good track acceptance and good momentum resolution, each track was required to have a transverse momentum greater than 0.1 GeV and to reach at least the third superlayer of the CTD.
To suppress the combinatorial background, a cut on the ratio p T (D * + , D + , D 0 )/E θ>10 • T , motivated by the hard character of charm fragmentation, was applied. The transverse energy, E θ>10 • T , was calculated as Σ i,θ i >10 • (E i sin θ i ), where the sum runs over all energy deposits in the CAL with the polar angle θ outside a cone of θ = 10 • around the forward direction. Moreover, the measured dE/dx values of those tracks that were candidates to come from D * + , D + and D 0 were used. The parametrisations of the dE/dx expectation values and the χ 2 probabilities l K and l π of the kaon and pion hypotheses, respectively, were obtained in the same way as described in previous publications [30,32]. To maximise the ratios of the numbers of correctly assigned kaons and pions to the square roots of the numbers of background particles, the cuts l K > 0.03 and l π > 0.01 were applied.
The measurements were done in the full kinematic range of Q 2 . Events produced in the photoproduction regime with Q 2 < 1 GeV 2 contributed 70 − 80 % of the selected D * + , D + and D 0 samples.

Reconstruction of D * + mesons
The D * + mesons were identified using the two decay channels The pion from the D * + decays is referred to as the "soft" pion, π s , because it is constrained to have limited momentum by the small mass difference between the D * + and D 0 [1].
Selected tracks were combined to form D 0 candidates assuming the decay channels (1) or (2). For both cases, D 0 candidates were formed by calculating the invariant mass M(Kπ) or M(Kπππ) for combinations having a total charge of zero. The soft pion was required to have a charge opposite to that of the particle taken as a kaon and was used to form a D * + candidate having mass M(Kππ s ) or M(Kππππ s ). To reduce the combinatorial background, requirements (see Table 1) similar to those used in a previous publication [32] were applied. To determine the background under the peaks, wrong-charge combinations were used. For both channels (1) and (2), these are defined as combinations with total charge ±2 for the D 0 candidate and total charge ±1 for the D * + candidate. The histograms in Fig. 1 show the ∆M distributions for the wrong-charge combinations, normalised to the distributions of D * + candidates with the appropriate charges in the range 0.15 < ∆M < 0.1685 GeV for channel (1) and 0.15 < ∆M < 0.16 GeV for channel (2). The upper ends of the normalisation ranges correspond to the trigger selections of D * + candidates in the two decay channels. The multiple counting of a D * + candidate produced by D 0 candidates formed by the same tracks was excluded [32].
To improve the signal-to-background ratio, only D * + candidates with 0.144 < ∆M < 0.147 GeV for channel (1) and 0.1445 < ∆M < 0.1465 GeV for channel (2) were kept for the excited charm and charm-strange meson studies. After background subtraction, signals of 39500±310 D * + mesons in channel (1) and 17300±210 D * + mesons in channel (2) were found in the above ∆M ranges.
The ∆M distributions were also fitted to a sum of a modified Gaussian function describing the signal and a background function. The modified Gaussian function was defined as where x = |(∆M − M 0 )/σ|. This functional form described both data and MC signals well. The signal position, M 0 , and width, σ, as well as the numbers of D * + mesons in the signal window were free parameters of the fit. The background function had a form where m π + is the pion mass [1] and A, B and C were free parameters. The fit yielded mass difference values of 145.46 ± 0.01 MeV for channel (1) and 145.45 ± 0.01 MeV for channel (2), in agreement with the PDG value [1]. The widths of the signals were 0.59±0.01 MeV and 0.51±0.01 MeV, respectively, reflecting the detector resolution.

Reconstruction of D + mesons
The D + mesons were reconstructed from the decay D + → K − π + π + . In each event, two tracks with the same charges and p T > 0.5 GeV and a third track with opposite charge and p T > 0.7 GeV were combined to form D + candidates. The pion masses were assigned to the two tracks with the same charges and the kaon mass was assigned to the third track, after which the candidate invariant mass, M(Kππ), was calculated. To suppress the combinatorial background, a cut of cos θ * (K) > −0.75 was imposed, where θ * (K) is the angle between the kaon in the Kππ rest frame and the Kππ line of flight in the laboratory frame. To further suppress the combinatorial background, a cut p T (D + )/E θ>10 • T > 0.25 was applied. To suppress background from D * + decays, combinations with M(Kππ) − M(Kπ) < 0.15 GeV were removed. The background from D + s → φπ + with φ → K + K − was suppressed by requiring that the invariant mass of any two D + candidate tracks with opposite charges was not within ±8 MeV of the nominal φ mass when the kaon mass was assigned to both tracks. Only D + candidates in the kinematic range p T (D + ) > 2.8 GeV and −1.6 < η(D + ) < 1.6 were kept for further analysis. Figure 2a shows the M(Kππ) distribution for the D + candidates after all cuts. Reflections from D + s and Λ + c decays to three charged particles were subtracted using the simulated reflection shapes normalised to the D + s and Λ + c production rates previously measured by ZEUS [30]. A clear signal is seen at the nominal value of the D + mass. To improve the signal-to-background ratio, only D + candidates with 1.850 < M(Kππ) < 1.890 GeV were kept for the excited charm meson studies. The mass distribution was fitted to a sum of a modified Gaussian function describing the signal and a linear function describing the non-resonant background. The fit yielded a D + mass value 1867.9±0.5 MeV in agreement with the PDG value [1]. The width of the signal was 12.9±0.5 MeV, reflecting the detector resolution. The number of D + mesons yielded by the fit in the above M(Kππ) range was N(D + ) = 20430 ± 510.

Reconstruction of D 0 mesons
The D 0 mesons were reconstructed from the decay D 0 → K − π + . In each event, tracks with opposite charges and p T > 0.8 GeV were combined in pairs to form D 0 candidates. To suppress the combinatorial background, a cut of | cos θ * (K)| < 0.85 was imposed, where θ * (K) is the angle between the kaon in the Kπ rest frame and the Kπ line of flight in the laboratory frame. To further suppress the combinatorial background, a cut p T (D 0 )/E θ>10 • T > 0.25 was applied.
For selected D 0 candidates, a search was performed for a track that could be the soft pion in a D * + → D 0 π + s decay. The soft pion was required to have p T > 0.1 GeV and a charge opposite to that of the particle taken as a kaon. The corresponding D 0 candidate was rejected if the mass difference, ∆M = M(Kππ s ) − M(Kπ), was below 0.15 GeV. All remaining D 0 candidates were considered "untagged", i.e. not originating from identified D * + decays. Only D 0 candidates in the kinematic range p T (D 0 ) > 2.8 GeV and −1.6 < η(D 0 ) < 1.6 were kept for further analysis. Figure 2b shows the M(Kπ) distribution for untagged D 0 candidates after all cuts. A reflection, produced by D 0 mesons with the wrong (opposite) kaon and pion mass assignment, was subtracted using the rejected sample of the D 0 mesons originating from D * + decays [30]. A clear signal is seen at the nominal value of the D 0 mass. To improve the signal-to-background ratio, only D 0 candidates with 1.845 < M(Kπ) < 1.885 GeV were kept for the excited charm-strange meson studies. The mass distribution was fitted to a sum of a modified Gaussian function describing the signal and a background function. Monte Carlo studies showed that the background shape was compatible with being linear in the mass range above the signal. For smaller M(Kπ) values, the background shape exhibits an exponential enhancement due to contributions from other D 0 decay modes and other D mesons. Therefore the background shape in the fit was described by the form 5 Study of the excited charm mesons D 0 1 and D * 0 To reconstruct the D 0 1 , D * 0 2 → D * + π − decays, an excited charm meson candidate was formed by combining each selected D * + candidate (Section 4.1) with an additional track, assumed to be a pion (π a ), with a charge opposite to that of the D * + candidate. The additional track was required to satisfy the pion dE/dx hypothesis with l π > 0.01 (Section 4). To reduce the combinatorial background, the following requirements were applied: for the D * + decay channel (1), and for the D * + decay channel channel (2). The decay angle θ * (D * + ) is the angle between the D * + in the D * + π a rest frame and the D * + π a line of flight in the laboratory frame. A cut η(π a ) < 1.1 was applied to exclude the region of large track density in the forward (proton) direction.   (1) and (2). A clear enhancement is seen in the range 2.4 < M(D * + π a ) < 2.5 GeV, where contributions from D 1 (2420) 0 and D * 2 (2460) 0 mesons are expected. The wide D 1 (2430) 0 meson, which is also expected to contribute to the M(D * + π a ) distribution, is not distinguishable from background due to its large width (384 +107 −75 ± 74 MeV [1]). No enhancement is seen in the M(D * + π a ) distribution for wrong charge combinations (histogram) formed by combining a D * + candidate and π a with the same charges. The wrong charge distribution lies generally below the distribution for the combinations with the appropriate charges, in agreement with MC predictions; this is expected near threshold since, due to charge conservation, the invariant mass distribution for random track combinations with total charge ±2 should lie below that for track combinations with total charge zero.

Reconstruction of
To reconstruct the D * 0 2 → D + π − decays, an excited charm meson candidate was formed by combining each selected D + candidate (Section 4.2) with an additional track, assumed to be a pion (π a ), with a charge opposite to that of the D + candidate. The additional track was required to satisfy the pion dE/dx hypothesis with l π > 0.01 (Section 4). To reduce the combinatorial background, the following requirements were applied: where θ * (D + ) is the angle between the D + in the D + π a rest frame and the D + π a line of flight in the laboratory frame.  Figure 3b shows the M(D + π a ) distribution for the selected excited charm meson candidates. A small excess is seen around the nominal mass of the D * 0 2 meson. The wide D * 0 (2400) 0 meson, which is also expected to contribute to the M(D + π a ) distribution, is not distinguishable from background due to its large width (261 ± 50 MeV [1]). As expected from parity and angular momentum conservation for a 1 + state, no indication of the D 0 1 decay to D + π − is seen. Feed-downs from the D 0 1 and D * 0 2 mesons decaying to D * + π − with a consequent D * + decay to a D + and undetected neutrals, predicted by MC at M(D + π a ) ∼ 2.3 GeV, are not seen, probably due to the large combinatorial background. No signal is seen in the M(D + π a ) distribution for wrong charge combinations (histogram) formed by combining a D + candidate and a π a with the same charges.

Mass, width and helicity parameters
To distinguish the D 0 1 (1 + state from j = 3/2 doublet) and D * 0 2 (2 + state from j = 3/2 doublet) mesons from each other and from the wide D 1 (2430) 0 (1 + state from j = 1/2 doublet) meson, the helicity angular distribution was used. The helicity angle (α) is defined as the angle between the π a and π s momenta in the D * + rest frame. The helicity angular distribution can be parametrised as where h is the helicity parameter. HQET predicts h = 3 (h = 0) for the 1 + state from the j = 3/2 (j = 1/2) doublet, and h = −1 for the 2 + state from the j = 3/2 doublet. Figure 4 shows the M(D * + π a ) distribution in four helicity intervals. The D 0 1 -meson contribution is increasing with | cos(α)| and dominates the excess in the M(D * + π a ) distribution for | cos(α)| > 0.75. The dependence of the D * 0 2 -meson contribution on the helicity angle is less pronounced; it is consistent with the expected slow decrease with | cos(α)|.
To extract the D 0 1 and D * 0 2 yields and properties, a minimal χ 2 fit was performed using simultaneously the M(D + π a ) distribution ( Fig. 3b) and the M(D * + π a ) distributions in four helicity intervals (Fig. 4). Each of the D 0 1 → D * + π − , D * 0 2 → D * + π − and D * 0 2 → D + π − signals was represented in the fit by a relativistic D-wave Breit-Wigner function (see Appendix) convoluted with a Gaussian resolution function with a width fixed to the corresponding MC prediction. The dependence of the detector acceptance and resolution on the M(D * + π a ) or M(D + π a ) was obtained from MC and corrected for in the fit function. Equation (4) was used to describe the helicity distributions. The acceptance dependence on the helicity angle, found from MC to be very weak, was corrected for in the fit function. Yields of all three signals, the D 0 1 and D * 0 2 masses, and the D 0 1 width and helicity parameters were free parameters of the fit. Since the data were not able to constrain reliably the D * 0 2 width and helicity parameter, the D * 0 2 width was fixed to the recently updated world average value of 43 ± 4 MeV [1] and the HQET prediction, h(D * 0 2 ) = −1, was used for the helicity parameter.
To describe backgrounds in the M(D * + π a ) and M(D + π a ) distributions, a functional form with three shape parameters x A exp(−Bx + Cx 2 ), where x = ∆M ext − m π + , was used. It was checked that such a functional form describes the wrong charge distributions well. The yields and shape parameters of the M(D * + π a ) and M(D + π a ) background functions were independent free parameters of the fit. Since neither data nor MC demonstrated a sizeable background dependence on the helicity angle, the same background function was used for the M(D * + π a ) distributions in the four helicity intervals.
The expected feed-downs from D 0 1 , D * 0 2 → D * + π − → D + π − + neutrals (Section 5.2) were included in the M(D + π a ) fit function; the effect on the fit results was small. Contributions from the wide D 1 (2430) 0 and D * 0 (2400) 0 states were added to the M(D * + π a ) and M(D + π a ) fit, respectively. Their shapes were described with a relativistic S-wave Breit-Wigner function (see Appendix) convoluted with a Gaussian resolution function with widths fixed to the MC prediction. The masses and widths of the wide excited charm mesons were set to the world-average values [1]. The D 1 (2430) 0 yield was set to that of the narrow D 1 (2420) 0 meson since both have the same quantum numbers. The D * 0 (2400) 0 yield was set to 1.7 times the D * 0 2 → D + π − yield as observed by the FOCUS collaboration [5]. The yield measured by FOCUS covers both a direct signal from the D * 0 (2400) 0 and a feed-down from the D 1 (2430) 0 , decaying to D * + π − with a consequent D * + decay to a D + and undetected neutrals [5].  Table 2.
The differences between the D 0 1 and D * 0 2 masses and M(D * + ) PDG were and, hence, the masses of the D 0 1 and D * 0 2 were The . The observed difference can be a consequence of differing production environments. The D 0 1 width can have a sizeable contribution from the broad S-wave decay even if the S-wave admixture is small [33,34]. A larger S-wave admixture at ZEUS with respect to that in measurements with restricted phase space, which can suppress production of the broad state, could explain why the measured D 0 1 width is larger than the world average value. The D 0 1 helicity parameter was This is inconsistent with the prediction for a pure S-wave decay of the 1 + state, h = 0. It is consistent with the prediction for a pure D-wave decay, h = 3.
In the general case of D-and S-wave mixing, the helicity angular distribution form of the 1 + state is: where r = Γ S /(Γ S + Γ D ), Γ S/D is the S-/D-wave partial width and φ is the relative phase between the two amplitudes. Using Eqs. (4) and (5), cos φ can be expressed in terms of r and the measured value of the helicity parameter, h: Figure 5 compares with previous measurements the range restricted by the measured h(D 0 1 ) value and its uncertainties in a plot of cos φ versus r. The ZEUS range has a marginal overlap with that restricted by the CLEO measurement of h(D 0 1 ) = 2.74 +1.40 −0.93 [35]. BELLE performed a three-angle analysis and measured both the cos φ and r values [4]. The BELLE measurement, which suggested a very small admixture of S-wave to the D 1 (2420) 0 → D * + π − decay and almost zero phase between two amplitudes, is outside the ZEUS range; the difference between the two measurements, evaluated with Eq. (6), is ∼ 2 standard deviations.

Fragmentation and branching fractions
The numbers of reconstructed D 0 1 , D * 0 2 → D * + π − and D * 0 2 → D + π − decays were divided by the numbers of reconstructed D * + and D + mesons, yielding the rates of D * + and D + mesons originating from the D 0 1 and D * 0 2 decays. To correct the measured rates for detector effects, the relative acceptances were calculated using the MC simulation as ratios of acceptances for the D 0 1 , D * 0 2 → D * + π − and D * 0 2 → D + π − states to the inclusive D * + and D + acceptances, respectively. The acceptance of the requirement l π > 0.01 for the additional track was calculated with data using identified pions from D * + decays (Section 4.1), to be (98.9 ± 0.1)%; only pions in the kinematic range of the additional pion selection were used.
Charm production at HERA is larger than beauty production by two orders of magnitude. The small b-quark relative contributions, predicted by the MC simulation using branching fractions of b-quark decays to the charm hadrons measured at LEP [36][37][38][39] 5 , were subtracted when calculating the relative acceptances; the subtraction changed the relative acceptances by less than 1.5% of their values. The relative acceptances were 52% for the D 0 1 , D * 0 2 → D * + π − and 47% for D * 0 2 → D + π − in the kinematic ranges described in Section 4.
The fractions, F, of D * + mesons originating from D 0 1 and D * 0 2 decays were calculated in the kinematic range |η(D * + )| < 1.6 and p T (D * + ) > 1.35 GeV for the D * + decay channel (1), combined with channel (2) for p T (D * + ) > 2.8 GeV: The fraction of D + mesons originating from D * 0 2 decays, calculated in the kinematic range p T (D + ) > 2.8 GeV and |η(D + )| < 1.6 is The fractions measured in the restricted p T (D * + , D + ) and η(D * + , D + ) kinematic ranges were extrapolated to the fractions in the full kinematic phase space using the Bowler modification [27] of the Lund symmetric fragmentation function [28] as implemented in Pythia [21]. Applying the estimated extrapolation factors, ∼ 1.1 for F D 0 In the full kinematic phase space, the extrapolated fractions of D * + originating from D 0 1 and D * 0 2 and of D + originating from D * 0 2 can be expressed as where the fragmentation fractions f (c → D 0 1 ), f (c → D * 0 2 ), f (c → D * + ) and f (c → D + ) are the rates of c quarks hadronising as a given charm meson, and B D 0 1 →D * + π − , B D * 0 2 →D * + π − and B D * 0 2 →D + π − are the corresponding branching fractions. These expressions provide a means to calculate the fragmentation fractions f (c → D 0 1 ) and f (c → D * 0 2 ), and the ratio of the two branching fractions for the D * 0 2 meson: The f (c → D * + ) and f (c → D + ) values, previously measured by ZEUS [30], were recalculated with the updated PDG values of the branching fractions [1] to be where the third uncertainties are due to the branching-fraction uncertainties. This yields in agreement with the world average value of 2.3 ± 0.6 [1]. Theoretical models [34,40,41] predict the ratio to be in the range from 1.5 to 3.
Assuming isospin conservation, for which Table 3). In order to check fragmentation universality for the excited charm mesons, the measured fragmentation fractions are compared and found to be consistent with those obtained in e + e − annihilations. The measured f (c → D 0 1 ) and f (c → D * 0 2 ) values are above the predictions of the thermodynamical model [42] (Table 3). The sum of the two fragmentation fractions, agrees with the prediction of the tunnelling model of 8.5% [43]. The predictions of both models are based on fits to the production rates of light-flavoured hadrons at LEP.
The ratio is consistent with the simple spin-counting prediction of 3/5. Both thermodynamical and tunnelling models suggest the ratio should exceed the spin-counting prediction due to the difference between the D 0 1 and D * 0 2 masses.
6 Study of the excited charm-strange meson D + s1 6.1 Reconstruction of D + s1 → D * + K 0 S decays The K 0 S mesons were reconstructed in their charged-decay mode, K 0 S → π + π − , for those events containing a D * + candidate. To identify K 0 S candidates, displaced secondary vertices reconstructed from pairs of oppositely charged tracks [44] were used. The identification efficiency degraded for the displaced secondary vertices close to the primary vertex. Therefore, additional secondary vertices were formed from pairs of oppositely charged tracks that were not assigned to one of the displaced secondary vertices. This was done by calculating the intersection points of the two tracks in the XY plane and requiring |∆Z| < 3 cm between the two tracks at the intersection point. To reduce the combinatorial background originating from tracks from the primary vertex, the additional secondary vertices with distances between the primary and secondary vertices in the XY plane of less than 0.5 cm were removed.
To reduce the combinatorial background, it was required that p T > 0.15 GeV for each track from any K 0 S candidate, cos α XY > 0.97 and cos α φZ > 0.85, where α XY and α φZ are the projected angles in the XY and φZ planes, respectively, between the K 0 S -candidate momentum and the line joining the primary to the secondary vertex. Figure 6 shows the invariant-mass, M(π + π − ), distribution for all remaining K 0 S candidates. Only K 0 S candidates with 0.480 < M(π + π − ) < 0.515 GeV were kept for the reconstruction of excited charm-strange mesons. The mass distribution was fitted to a sum of a modified Gaussian function describing the signal and a linear function describing the non-resonant background. The fit yielded the K 0 S mass value 497.8 ± 0.1 MeV, in agreement with the PDG value [1]. The width of the signal was 4.1±0.1 MeV reflecting the detector resolution. The number of reconstructed K 0 S mesons in the range 0.480 < M(π + π − ) < 0.515 GeV yielded by the fit was N(K 0 S ) = 8540 ± 120. To reconstruct the D + s1 → D * + K 0 S decays, a D + s1 -meson candidate was formed by combining each selected D * + candidate (Section 4.1) with the K 0 S candidates reconstructed in the same event. For each D + s1 candidate, the extended mass difference,

Reconstruction of D
Monte Carlo studies show that a signal from the D + s1 → D * 0 K + decay, with a consequent D * 0 decay to a D 0 and undetected neutrals, should be seen in the M(D 0 K + ) distribution with an average negative shift of 142.4±0.2 MeV with respect to the nominal D + s1 mass [1], and that the shape of the signal can be reasonably well described by the modified Gaussian function (Eq. 3) with a width of 3.1 MeV.
To reconstruct the D + s1 → D * 0 K + decays, an excited charm-strange meson candidate was formed by combining each selected untagged D 0 candidate (Section 4.3) with an additional track, assumed to be a kaon (K a ), with a charge opposite to that of the particle taken as a kaon to form the D 0 candidate. The additional track was required to satisfy the kaon dE/dx hypothesis with l K > 0.03 (Section 4). To reduce the combinatorial background, the following requirements were applied: where θ * (D 0 ) is the angle between the D 0 in the D 0 K a rest frame and the D 0 K a line of flight in the laboratory frame.  Figure 7b shows the M(D 0 K a ) distribution for the selected excited charm-strange meson candidates. A signal is seen at the expected position of the feed-down from the D + s1 → D * 0 K + decay. No signal from the known decay D s2 (2573) + → D 0 K + [1] was observed, probably due to the large combinatorial background.

Mass, width and helicity parameters
The M(D * + K 0 S ) distribution in four helicity intervals is shown in Fig. 8, with the helicity angle (α) defined as the angle between the K 0 S and π s momenta in the D * + rest frame. The D + s1 signal decreases with | cos(α)|. To extract the D + s1 yields and properties, an unbinned likelihood fit was performed using simultaneously values of M(D 0 K a ), M(D * + K 0 S ), and cos(α) for D * + K 0 S combinations. The observed narrow signals in the M(D * + K 0 S ) and M(D 0 K a ) distributions were described in the fit by a Gaussian function and a modified Gaussian function, respectively. Equation (4) was used to describe the helicity distribution. The acceptance dependence on the helicity angle, found from MC to be very weak, was corrected for in the fit function. The average shift of the signal in the M(D 0 K a ) distribution with respect to the mass of D + s1 meson was fixed to the MC prediction (Section 6.2). Yields and widths of both signals, the D + s1 mass and the D + s1 helicity parameter were free parameters of the fit. To describe the background in the M(D * + K 0 S ) distribution, a function x A , where x = ∆M ext , was used. The background description for the M(D 0 K a ) distribution required a functional form with two shape parameters x A exp(−Bx), where x = ∆M ext − m K + and m K + is the kaon mass [1]. The shape parameters of the M(D * + K 0 S ) and M(D 0 K a ) background functions were independent free parameters of the fit. Since neither data nor MC demonstrated a sizeable background dependence on the helicity angle, the background function for D * + K 0 S combinations was assumed to be helicity independent. The numbers of reconstructed D + s1 mesons and values of all free background parameters yielded by the fit are summarised in Table 4. ). The measured h value is inconsistent with the prediction for a pure D-wave decay of the 1 + state, h = 3, and is barely consistent with the prediction for a pure S-wave decay, h = 0. Figure 9 shows a range, restricted by the measured h(D + s1 ) value and its uncertainties, in a plot of cos φ versus r = Γ S /(Γ S + Γ D ) (Eq. 6). The measurement suggests a significant contribution of both D-and S-wave amplitudes to the D s1 (2536) + → D * + K 0 S decay. The ZEUS range agrees with that restricted by the CLEO measurement of h(D + s1 ) = −0.23 +0.40 −0.32 [45] and with the BELLE three-angle measurement of both cos φ and r values [46].

Fragmentation and branching fractions
The numbers of reconstructed D + s1 → D * + K 0 S and D + s1 → D * 0 K + decays were divided by the numbers of reconstructed D * + and untagged D 0 mesons, respectively, yielding rates of D * + and untagged D 0 mesons originating from D + s1 decays. To correct the measured rates for detector effects, the relative acceptances were calculated using the MC simulation as ratios of acceptances for the D + s1 → D * + K 0 S and D + s1 → D * 0 K − states to the inclusive D * + and untagged-D 0 acceptances, respectively. The untagged-D 0 acceptance included subtraction of a small contamination to N(D 0 untag ) from unidentified D * + mesons. The acceptance of the requirement l K > 0.03 for the additional track was calculated with data using identified kaons from D * + decays (Section 4.1), to be (95.3 ± 0.2)%; only the kaons from the kinematic range of the additional kaon selection were used. Subtraction of the small b-quark contribution changed the relative acceptances by less than 2.2% of their values. The relative acceptances were 38% for D + s1 → D * + K 0 S and 48% for D + s1 → D * 0 K + in the kinematic ranges described in Section 4.
The fraction, F , of D * + mesons originating from D + s1 decays, corrected to the fraction of K 0 mesons decaying as K 0 S (50%) and to the branching fraction of the K 0 S decay into π + π − (69.20% [1]), was calculated in the kinematic range |η(D * + )| < 1.6 and p T (D * + ) > 1. 35 GeV for the D * + decay channel (1), combined with channel (2) for p T (D * + ) > 2.8 GeV: The fraction of untagged D 0 mesons originating from D + s1 decays, calculated in the kinematic range p T (D 0 ) > 2.8 GeV and |η(D 0 )| < 1.6 is In the full kinematic phase space, the extrapolated fractions of D * + and untagged D 0 mesons originating from D + s1 can be expressed as where the fragmentation fractions f (c → D + s1 ), f (c → D * + ) and f (c → D 0 untag ) are the rates of c quarks hadronising as a given charm meson, and B D + s1 →D * + K 0 and B D + s1 →D * 0 K + are the corresponding branching fractions.
These expressions provide a means to calculate the fragmentation fraction f (c → D + s1 ) and the ratio of the two D + s1 branching fractions: .
Using f (c → D * + ) and f (c → D 0 ) [30], recalculated with the updated values of the branching fractions [1], and calculating the fragmentation fraction into untagged D 0 where B D * + →D 0 π + is the branching fraction of the decay D * + → D 0 π + (67.7% [1]) and the third uncertainty is due to the branching-fraction uncertainties, yields in comparison with the world average value of 1.27 ± 0.21 [1]. Isospin invariance requires the matrix elements of the two measured D + s1 decay modes to be the same, while an enhancement of the D * 0 K + final state is expected due to the larger phase space [41].
Assuming that the decay width of the D + s1 is saturated by the D * K final states, i.e. Table 3). The measured fragmentation fraction value agrees with those obtained in e + e − annihilations and is above the prediction of the thermodynamical model [42].
The ratio for the two 1 + states represents the strangeness-suppression factor for P -wave charm mesons. The measured value agrees with measurements of the strangeness-suppression factor for the lowest-mass charm mesons [13,30,47] and with the value of 0.3, used by default in simulations based on the Lund string fragmentation scheme [48].

Search for the radially excited charm meson D * ′+
To search for the D * ′+ → D * + π + π − decays, a D * ′+ candidate was formed by combining each selected D * + candidate (Section 4.1) with two additional tracks with opposite charges. The additional tracks were assumed to be pions (π ± a ), and were required to satisfy the pion dE/dx hypothesis with l π > 0.01 (Section 4). To reduce the combinatorial background, the cuts η(π ± a ) < 1.1 and cos θ * (D * + ) < 0.8 were imposed, where θ * (D * + ) is the angle between the D * + in the D * + π + a π − a rest frame and the D * + π + a π − a line of flight in the laboratory frame. To further reduce the combinatorial background, the following requirements were applied: for the D * + decay channel (1) and for the D * + decay channel channel (2).
For each D * ′+ candidate, the extended mass difference, ∆M ext = M(Kππ s π + a π − a ) − M(Kππ s ) or ∆M ext = M(Kππππ s π + a π − a ) − M(Kππππ s ), was calculated. The invariant mass of the D * + π + a π − a system was calculated as M(D * + π + a π − a ) = ∆M ext + M(D * + ) PDG . The resolution in M(D * + π + a π − a ) around 2.64 GeV, where a narrow signal was reported by the DELPHI Collaboration [10], was estimated from MC simulations to be 5.6 MeV. Figure 10 shows the M(D * + π + a π − a ) distribution below 2.9 GeV. The distribution was investigated in the full accessible range; no narrow resonance was observed.
An estimate of the fraction of D * + mesons originating from the D * ′+ → D * + π + π − decays was performed in the signal window of 2.59 < M(D * + π + a π − a ) < 2.69 GeV. This window covers both theoretical predictions [9] and the DELPHI measurement [10]. The M(D * + π + a π − a ) distribution was fitted outside the signal window to the background functional form with two shape parameters, x A exp(−Bx), where x = ∆M ext − 2m π + . The number of reconstructed D * ′+ mesons was estimated to be 104 ± 83 by subtracting the background function, integrated over the signal window, from the observed number of candidates in the window.
The number of reconstructed D * ′+ → D * + π + π − decays was divided by the number of reconstructed D * + mesons, yielding a fraction of D * + mesons originating from the D * ′+ decays. To correct the measured fraction for detector effects, the relative acceptance was calculated using the MC simulation (Section 3) as a ratio of an acceptance for the D * ′+ → D * + π + π − state to the inclusive D * + acceptance. The acceptance of the requirement l π > 0.01 for the additional tracks was calculated with data (Section 5.4). Subtraction of the small b-quark contribution, performed under a conservative assumption that all D * ′+ mesons are produced in charm fragmentation, changed the relative acceptance by ∼ 1.7% of its value. The relative acceptance was found to be 34% in the kinematic range described in Section 4.1.
In the full kinematic phase space, the extrapolated ratio can be expressed as where the fragmentation fraction f (c → D * ′+ ) is the rate of c quarks hadronising as D * ′+ , and B D * ′+ →D * + π + π − is the branching fraction of the decay D * ′+ → D * + π + π − .
The upper limit is the frequentist confidence bound calculated assuming a Gaussian probability function in the unified approach [49]. It is stronger than the 0.9% limit on D * ′± production in charm fragmentation obtained by OPAL [11].
The ratio of the D * ′+ → D * + π + π − to D 0 1 , D * 0 2 → D * + π − decay yields, calculated as is compared with those obtained by DELPHI [10] and OPAL [11] in Table 5. The ZEUS measurement is more sensitive to the existence of a narrow resonance decaying to D * + π + π − . However, it is sensitive only to the resonance production in charm fragmentation while the LEP measurements are also sensitive to beauty fragmentation.

Systematic uncertainties
The systematic uncertainties of the measured values were determined by varying the analysis procedure and repeating all calculations. The sizes of the variations were chosen commensurate with the estimated uncertainties of the relevant parameters and variables.
The following groups of systematic uncertainties were considered.
• {δ 1 } The uncertainties related to the signal and helicity extraction procedures were obtained as follows: for the D * + signals: the ranges for the background normalisation were reduced by 2 MeV on either side; the fit was used instead of the subtraction procedure; for the D + signal: the range for the signal fit was reduced by 20 MeV on either side; the amounts of the subtracted D + s and Λ + c reflections were varied in the range of their uncertainties; a higher-order polynomial was included in the background parametrisation; for the untagged D 0 signal: the range for the signal fit was reduced by 20 MeV on either side; the value of M(Kπ), where the background form with the exponential enhancement turns into the linear form, was varied between 1.84 GeV and 1.88 GeV; a higher-order polynomial was included in the background parametrisation; for the D 0 1 and D * 0 2 signals: the ranges for the signal fit were reduced by 20 MeV on either side; higher-order polynomials were included in the exponential of the background parametrisations; the masses and widths of the wide excited charm mesons were varied in the range of their uncertainties [1] and their yields were varied by ±50%; for the D 0 1 helicity distribution: the acceptance dependence on the helicity angle was varied in the range of its uncertainty; the background functions in the four helicity intervals were allowed to have separate normalisations; for the D + s1 signals: the ranges for the signal fit were reduced by 12 MeV on the upper side; higher-order polynomials were included in the exponential of the background parametrisations; the average shift of the signal in the M(D 0 K a ) distribution with respect to the mass of D + s1 meson was varied in the range of its uncertainty (Section 6.2); for the D + s1 helicity distribution: the acceptance dependence on the helicity angle was varied in the range of its uncertainty; the background function was allowed to have a free helicity parameter; for the D * ′+ signal search: the range for the background fit was reduced by 12 MeV on the upper side; a higher-order polynomial was included in the exponential of the background parametrisation; • {δ 2 } The uncertainty of the tracking reconstruction and simulation was taken into account by varying all momenta by ±0.1% (magnetic field uncertainty) and by changing the track momentum and angular resolutions by ±5% of their values. • {δ 5 } The uncertainty of the CAL simulation was determined by varying the CAL energy scale by ±2%.
• {δ 6 } The uncertainties of the fragmentation fractions f (c → D * + ), f (c → D + ) and f (c → D 0 untag ) were determined by adding in quadrature their statistical and systematic uncertainties and the errors originating from the branching-fraction uncertainties. The uncertainty of the branching fraction of the K 0 S decay into π + π − [1] was also taken into account. • {δ 8 } The uncertainty of the beauty subtraction was determined by varying the bquark cross section by a factor of two in the MC sample and by varying the branching fractions of b-quarks to charm hadrons by their uncertainties [36][37][38][39].
• {δ 9 } The extrapolation uncertainties were determined by varying relevant parameters of the Pythia simulation using the Bowler modification [27] of the Lund symmetric fragmentation function [28] 6 . The following variations were performed: the mass of the c quark was taken to be 1.5 ± 0.2 GeV; the strangeness suppression factor was taken to be 0.3 ± 0.1; the fraction of the lowest-mass charm mesons produced in a vector state was taken to be 0.6 ± 0.1; production rates of the excited charm and charm-strange mesons were varied by ±50% around the central values tuned to reproduce the measured fractions of c quarks hadronising into D 0 1 , D * 0 2 or D + s1 ; -the Bowler fragmentation function parameter r c was varied from the predicted value 1 to 0.5; the a and b parameters of the Lund symmetric function were varied by ±20% around their default values [21].
Contributions from the different systematic uncertainties were calculated and added in quadrature separately for positive and negative variations. The results are given in Tables 6-7.
The relatively narrow ∆M, M(Kππ) and M(Kπ) ranges, used for the excited charm and charm-strange meson studies, selected only the central parts of the D * + , D + and D 0 signals, respectively (Section 4). It was checked that increasing the narrow ranges by 25 − 50% produced no effect on the results beyond the expected statistical fluctuations. Similarly, no systematic shifts were found when removing the η(π a , K a ) < 1.1 requirement from the excited state selections (Sections 5.1, 5.2, 6.2 and 7). It was also checked that the D 0 1 width value cannot be significantly reduced by including an interference between the signal and background. The measured D 0 1 helicity parameter is which is inconsistent with the prediction of h = 0 for a pure S-wave decay of the 1 + state, and is consistent with the prediction of h = 3 for a pure D-wave decay. In the general case of D-and S-wave mixing, the allowed region of the mixing parameters is consistent with the CLEO measurement [35] and marginally consistent with the BELLE result [4].
The measured D + s1 helicity parameter is This value is inconsistent with the prediction of h = 3 for a pure D-wave decay of the 1 + state, and is barely consistent with the prediction of h = 0 for a pure S-wave decay. The measurement suggests a significant contribution of both D-and S-wave amplitudes to the D s1 (2536) + → D * + K 0 S decay. The allowed region of the mixing parameters is consistent with the CLEO measurement [45] and with the BELLE result [46].
The ratios of the dominant D * 0 2 and D + s1 branching fractions are −0.6 (syst.), in agreement with the world average values [1].
The fractions of c quarks hadronising into D 0 1 , D * 0 2 or D + s1 mesons are consistent with those obtained in e + e − annihilations (Table 3), in agreement with charm fragmentation universality. Sizeable fractions of the D * + , D + and D 0 mesons emanate from these excited states.
total              The distribution of M(D * ± π + a π − a ) = ∆M ext + M(D * + ) PDG , where ∆M ext = M(Kππ s π + a π − a ) − M(Kππ s ) or ∆M ext = M(Kππππ s π + a π − a ) − M(Kππππ s ), for D * ′± → D * ± π + π − candidates (dots). The inset shows the D * ′± signal window covering both theoretical predictions and the DELPHI measurement. The solid curve is a fit to the background function outside the signal window. The shaded histogram shows the Monte Carlo D * ′± signal, normalised to the obtained upper limit (95% C.L.) and shown on top of the fit interpolation (dashed curve).