Exclusive Production of $D^+_sD^-_s$, $D^{*+}_sD^-_s$, and $D^{*+}_sD^{*-}_s$ via $e^+ e^-$ Annihilation with Initial-State-Radiation

We perform a study of exclusive production of$D^+_sD^-_s$, $D^{*+}_sD^-_s$, and $D^{*+}_sD^{*-}_s$ final states in initial-state-radiation events from $e^+ e^-$ annihilations at a center-of-mass energy near 10.58 GeV, to search for charmonium $1^{--}$ states. The data sample corresponds to an integrated luminosity of 525 $fb^{-1}$ and was recorded by the BaBar experiment at the PEP-II storage ring. The $D^+_sD^-_s$, $D^{*+}_sD^-_s$, and $D^{*+}_sD^{*-}_s$ mass spectra show evidence of the known $\psi$ resonances. Limits are extracted for the branching ratios of the decays $X(4260)\to D^{(*)+}_sD^{(*)-}_s$.


I. INTRODUCTION
The surprising discovery of new states decaying to J/ψπ + π − [1,2] has renewed interest in the field of charmonium spectroscopy, since not all the new resonances are easy to accommodate in the quark model. Specifically, the BABAR experiment discovered a broad state, X(4260), decaying to J/ψπ + π − , in the initial- * Now at Temple University, Philadelphia, Pennsylvania 19122, USA † Also with Università di Perugia, Dipartimento di Fisica, Perugia, Italy ‡ Also with Università di Roma La Sapienza, I-00185 Roma, Italy § Now at University of South Alabama, Mobile, Alabama 36688, USA ¶ Also with Università di Sassari, Sassari, Italy state-radiation (ISR) reaction e + e − → γ ISR X(4260). Its quantum numbers J P C = 1 −− are inferred from the single virtual-photon production. Enhancements in the ψ(2S)π + π − mass distribution at 4.36 GeV/c 2 [3,4] and 4.66 GeV/c 2 [4] have been observed for the reaction e + e − → γ ISR ψ(2S)π + π − . Charmonium states at these masses would be expected [5,6] to decay predominantly to DD, D * D, or D * D * . It is peculiar that the decay rate to the hidden charm final state J/ψπ + π − is much larger for the X(4260) than for the higher-mass charmonium states (radial excitations) [7]. Many theoretical interpretations for the X(4260) have been proposed, including unconventional scenarios: quark-antiquark gluon hybrids [9], baryonium [10], tetraquarks [11], and hadronic molecules [12]. If the X(4260) were a diquark-antidiquark state [cs][cs], as proposed by L. Maiani et al. [11], this state would predominantly decay to D + s D − s . For a discussion and a list of references see, for example, Ref. [13].
In this paper, we present a study of the ISR production of D + s D − s , D * + s D − s , and D * + s D * − s [14] final states, and search for evidence of charmonium states and resonant structures. This follows earlier BABAR measurements of the cross section of DD [15] and of D * D and D * D * production [16] and studies of these final states [17,18] by the Belle Collaboration. Recently the CLEO Collaboration [19] studied e + e − annihilation to D + s D − s , D * + s D − s , and D * + s D * − s final states at center of mass energies from threshold up to 4.3 GeV GeV/c 2 . In the present analysis we extend these measurements up to 6.2 GeV/c 2 .
This paper is organized as follows. A short description of the BABAR experiment is given in Section II, and data selection is described in Section III. In Sections IV, V, and VI we present studies of the D + s D − s , D * + s D − s , and D * + s D * − s final states, respectively. Fits to the three final states are described in Section VII, and in Section VIII we present limits on the decay of the X(4260) to D ( * )+ s D ( * )− s . A summary and conclusions are found in Section IX.

II. THE BABAR EXPERIMENT
This analysis is based on a data sample of 525 fb −1 recorded mostly at the Υ(4S) resonance and 40 MeV below the resonance by the BABAR detector at the PEP-II asymmetric-energy e + e − storage rings. The sample includes also 15.9 fb −1 and 31.2 fb −1 collected at the Υ(2S) and Υ(3S) respectively, and 4.4 fb −1 above the Υ(4S) resonances. The BABAR detector is described in detail elsewhere [20]. We mention here only the components of the detector that are used in the present analysis. Charged particles are detected and their momenta measured with a combination of a cylindrical drift chamber (DCH) and a silicon vertex tracker (SVT), both operating within a 1.5 T magnetic field of a superconducting solenoid. Information from a ring-imaging Cherenkov detector is combined with specific ionization measurements from the SVT and DCH to identify charged kaon and pion candidates. The efficiency for kaon identification is 90% while the rate for a kaon being misidentified as a pion is 2%. Photon energies are measured with a CsI(Tl) electromagnetic calorimeter (EMC).

III. DATA SELECTION
For each candidate event, we first reconstruct a D + s D − s pair. While one of the D + s is required to decay to K + K − π + , we include three different decay channels for the second D − s (see Table I). D * + s decays are reconstructed via their decay D * + s → D + s γ. For all final states, events are retained if the number of well-measured charged tracks having a transverse momentum greater than 0.1 GeV/c is exactly equal to the total number of charged daughter particles. EMC clusters with a minimum energy of 30 MeV that are not as- sociated with a charged track are identified as photons.
Candidates for the decay π 0 → γγ are kinematically constrained to the π 0 mass. For K 0 S → π + π − candidates we apply vertex and mass constraints. The tracks corresponding to the charged daughters of each D + s candidate are constrained to come from a common vertex. Reconstructed D + s candidates with a fit probability greater than 0.1% are retained. Each D + s D − s pair is refit to a common vertex with the constraint that the pair originates from the e + e − interaction region. Only candidates with a χ 2 fit probability greater than 0.1% are retained. For each event we consider all combinations.
ISR Monte Carlo (MC) [21] events for each final state are fully simulated using the GEANT4 detector simulation package [22], and they are processed through the same reconstruction and analysis chain as the data.
We select D + s and D * + s candidates using the reconstructed D + s mass and the mass difference, which for D + s parameters are obtained by fitting the relevant mass spectra using a polynomial for the background and a single Gaussian for the signal. For D * + s , we use the PDG [8] mass and a Gaussian width σ = 6 MeV/c 2 obtained by MC simulations. Events are selected within ±2.0σ from the fitted central values. The D ( * )+ s candidate three-momentum is determined from the summed three-momenta of its decay particles. The nominal D ( * )+ s mass [8] is used to compute the energy component of its four-momentum.
The ISR photon is preferentially emitted at small angles with respect to the beam axis, and it escapes detection in the majority of ISR events. Consequently, the ISR photon is treated as a missing particle. We define the squared mass M 2 rec recoiling against the D + s D − s , D * + s D − s , and D * + s D * − s systems using the four-momenta of the beam particles p e ± and of the reconstructed D This quantity should peak near zero for both ISR events and for exclusive production of e + e − → D  We allow additional π 0 and photon candidates due to radiative or background photons. This introduces multiple candidates in the reconstruction of the different channels. For channel (1)-(2) ambiguities, each [K + K − π + (π 0 )] combination is considered as a candidate for both channels. For all channels, each To discriminate among the different D + s channels and D ( * )+ s D ( * )− s final states, and to separate signal from background, we make use of a likelihood ratio test: where N is the number of discriminating variables, while P DF S i and P DF B i are normalized distributions describing signal and background, respectively. Signal P DF S i are obtained from MC simulations. Background P DF B i are obtained from the data. Since the ISR signal is very small compared to the entire data set of candidates (< 0.1 %), we use the data as the background model by relaxing all the selection criteria, except m(D The discriminating variables used in the likelihood ratio are the following. • The number of additional π 0 candidates in the event. For decay channel (2) this number is computed after removing the π 0 from D + s decay. This distribution is expected to peak at zero for signal events.
• The residual energy in the calorimeter, which is computed after removing any ISR photon candidate, identified by a center-of-mass energy greater s decays is excluded from the residual energy calculation. This distribution is expected to peak at zero for signal events.
• The distribution of cos θ * , where θ * is the polar angle of the D ( * )+ s D ( * )− s system in the center-ofmass frame which peaks at ±1 for ISR events.
• The momentum distribution of the π 0 from the D + s for decay channel (2).   polynomial. The signal lineshape is taken from D + s D − s MC simulations. The resulting yield and the fitted purity P , defined as P = N signal /(N signal + N background ) are summarized in Table II.
The D + s D − s mass spectrum, presented in Fig. 2(a), shows a threshold enhancement at the position of the ψ(4040) and a small enhancement around 4.26 GeV/c 2 .
We make use of the Gaussian functions to describe the presence of the peaking backgrounds. The D + s D − s background, taken from M 2 rec sideband events (1.5 < |M 2 rec | < 3.5 GeV 2 /c 4 ), is fitted to a sum of a Gaussian function and a 3 rd order polynomial. The fitted D + s D − s mass spectrum for these events, normalized to the background estimated from the fit to the M 2 rec distribution, is presented as the shaded distribution in Fig. 2(a).
The D + s D − s reconstruction efficiency and the mass resolution for each channel have been studied in the mass region between 4.25 and 6.25 GeV/c 2 . The D + s D − s mass resolution is similar for decay channels (1) Table I, The ǫ B function for D + s D − s is shown in Fig. 3(a). The three D + s D − s decay channels, after correcting for efficiency and branching fractions, have yields that are consistent within statistical errors. The D + s D − s cross section is computed using s is the background-subtracted yield. The differential luminosity is computed as [23] dL dm D +  Fig. 4(a). This result can be compared with QCD calcu- lations ref. [24], which predict a vanishing cross section near 5 GeV/c 2 . The list of systematic uncertainties for the D + s D − s cross section is summarized in Table III and it is evaluated to be 23%. It includes contributions from particle identification, tracking, photon and π 0 reconstruction efficiencies, background estimates, branching fractions, and the criteria to select the final state. All contributions are added in quadrature. The D + s D − s systematic error is dominated by the uncertainty in the veto of the D * + s D − s events. The M 2 rec distribution for D * + s D − s candidates is shown in Fig. 1(b) where a clear signal of ISR production is observed. The number of ISR candidates and sample purity are summarized in Table II.
The D * + s D − s mass spectrum and background are shown in Fig. 2(b) and is dominated by the ψ(4160) resonance. The D * + s D − s mass resolution is similar for the three decay channels and increases with D * + s D − s mass from 7 to 8 MeV/c 2 in the mass region of the ψ resonances. The weighted efficiency ǫ B is shown in Fig. 3(b). The D * + s D − s cross section is calculated using the method described in Sec. IV for D + s D − s . The result is shown in Fig. 4(b). The overall systematic error for the cross section is 13% and is dominated by the uncertainties in the branching fractions [8] (see Table III).

FINAL STATE
For the selection of D * + s D * − s candidates, we do not make use of the likelihood test described in the previous sections because no improvement for the signal to background ratio is obtained. Instead, we require the two photon invariant mass m(γγ) to lie outside the π 0 window.
The ambiguity in the γ's assignment to D * + s or D * − s is resolved by choosing the D + s γ-D − s γ combinations with both D ( * )± s masses closest the expected value. Fig. 1(c) shows the resulting M 2 rec distribution which shows clear evidence for the signal final state produced in interactions with ISR. For the selected the ISR signal candidates, we show the ∆m(D + s γ) distribution (two combinations per event) in Fig. 5(b). The resulting event yield and purity are summarized in Table II. The D * + s D * − s mass spectrum and background are shown in Fig. 2(c). Due to the presence of structures, the background in this case is fitted using a 3 rd order polynomial and two Gaussians. Monte Carlo studies indicate that an important part of this background is due to the D * + s D − s final state plus a random background γ. The D * + s D * − s mass resolution is similar for the three decay channels. It increases with D * + s D * − s mass from 9 to 11 MeV/c 2 in the mass region of the ψ resonances. The weighted efficiency ǫ B is shown in Fig. 3(c). The D * + s D * − s cross section shown in Fig. 4(c) is calculated using the same method used to compute the D + s D − s cross section. The overall uncertainty for the cross section is 13% and is dominated by the uncertainties on the D ( * )± s branching bractions (see Table III). The where m is the D ( * )+ s D ( * )− s mass, c i and φ i are free parameters, W i (m) are P-wave relativistic Breit-Wigner distributions [8], P (m) represents the nonresonant contribution, B(m) is the background described in Sect. IV, ǫ B (m) is the weighted efficiency, and f is the signal fraction fixed to the values obtained fitting the M 2 rec distributions. In this way we allow interference between the resonances and the nonresonant contribution P (m). The shape of the nonresonant contribution P (m) is unknown; we therefore parametrize it in a simple way as where C(m) is the phase space function for D ( * )+ s D ( * )− s , and a and b are free parameters. The size of the nonresonant production is determined by the fit.
The mass and width of the ψ(4040), ψ(4160), ψ(4415) and X(4260) are fixed to the values reported in [8]. Resolution effects can be ignored since the widths of the resonances are much larger than the experimental resolution.
The three D + s D − s , D * + s D − s , and D * + s D * − s likelihood functions are computed with different thresholds, efficiencies, purities, backgrounds, and numbers of contributing resonances appropriate for each final state. The results of this fits are compared to the data in Fig. 2, both the total fitted yield as well as the coherent nonresonant contribution, |P (m)| 2 , ignoring any interference effects. The fraction for each resonant contribution i is defined by the following expression: The fractions f i do not necessarily add up to 1 because of interference between amplitudes. The error for each fraction has been evaluated by propagating the full covariance matrix obtained by the fit. The resulting fit fractions are given in Table IV mass spectra include the X(4260) resonance, which is allowed to interfere with all the other terms. In all cases, the X(4260) fraction is consistent with zero. We note that the weak enhancement around 4.26 GeV/c 2 in the D + s D − s mass spectrum is described by the fit in terms of interference between the ψ(4040) and ψ(4160) resonances.

VIII. LIMITS ON X(4260)
The X(4260) yields are used to compute the cross section times branching fraction, to be compared with a BABAR measurement of the J/ψπ + π − final state [2]. The fractions from the fits reported in Table IV are converted to yields which are divided by the mass dependent ǫ B efficiency and the integrated luminosity. Systematic errors due to the mass and the width of the ψ(4040), ψ(4160), ψ(4415), and X(4260) resonances are evaluated by varying the masses and widths by their uncertainty in the fit. The size of the background contributions is varied within the statistical error, and the meson radii in the Breit-Wigner terms [25] are varied between 0 and 2.5 GeV −1 . Statistical and systematic errors are added in quadrature. We obtain B(X(4260) → D + s D − s ) B(X(4260) → J/ψπ + π − ) < 0.7, B(X(4260) → D * + s D − s ) B(X(4260) → J/ψπ + π − ) < 44, and B(X(4260) → D * + s D * − s ) B(X(4260) → J/ψπ + π − ) < 30, (12) at the 95% confidence level.

IX. TOTAL CROSS SECTION AND CONCLUSION
The sum of the e + e − → D + s D − s , e + e − → D * + s D − s , and e + e − → D * + s D * − s cross sections is shown in Fig. 6; the arrows indicate the position of the different ψ resonances and the X(4260). At the X(4260) mass, there is a local minimum, similar to the measured cross section for hadron production in e + e − annihilation [8]. In conclusion, we have studied the exclusive ISR production of the D + s D − s , D * + s D − s , and D * + s D * − s final states. The mass spectra show production of the J P C = 1 −− states, ψ(4040), ψ(4160) and a weak indication for a smaller enhancement near 4.3 GeV. From fits to the mass spectra for the three different final states we have determined contributions by different cc resonances.
Upper limits on X(4260) decays to these final states relative to J/ψπ + π − are computed. If the X(4260) is a 1 −− charmonium state, it should decay predominantly to open charm. Within the present limited data sample