A Study of the DsJ(2317) and DsJ(2460) Mesons in Inclusive ccbar Production Near sqrt(s) = 10.6 GeV

A study of the $D_{sJ}^*(2317)^+$ and $D_{sJ}(2460)^+$ mesons in inclusive $c\bar{c}$ production is presented using 232 ${\rm fb}^{-1}$ of data collected by the BaBar experiment near $\sqrt{s} = 10.6$ \gev. Final states consisting of a $D_s^+$ meson along with one or more $\pi^0$, $\pi^\pm$, or $\gamma$ particles are considered. Estimates of the mass and limits on the width are provided for both mesons and for the $D_{s1}(2536)^+$ meson. A search is also performed for neutral and doubly-charged partners of the $D_{sJ}^*(2317)^+$ meson.

and limits on the width are provided for both mesons and for the Ds1(2536) + meson. A search is also performed for neutral and doubly-charged partners of the D * sJ (2317) + meson.

I. INTRODUCTION
The D * sJ (2317) + meson, discovered by this collaboration [1] and confirmed by others [2,3], does not conform to conventional models of cs meson spectroscopy. Included with this discovery were suggestions of a second state, the D sJ (2460) + meson. This meson, observed by CLEO [2] and confirmed by this collaboration [4] and BELLE [3], has a mass that is also lower than expectations [5][6][7][8]. Because the masses of these two states are so unusual there has been speculation [9,10] that both possess an exotic, four-quark component. The possibility that the D * sJ (2317) + and D sJ (2460) + are exotic has attracted considerable experimental and theoretical interest and has focused renewed attention on the subject of charmed-meson spectroscopy in general.
Presented in this paper is an updated analysis of these two states using 232 fb −1 of e + e − → cc data collected by the BABAR experiment. From this analysis new estimates of the D * sJ (2317) + and D sJ (2460) + masses, limits on their intrinsic widths, calculations on their production cross sections, and the branching ratios of D sJ (2460) + decays to D + s γ and D + s π + π − with respect to its decay to D + s π 0 γ are presented. These measurements are performed by fitting the invariant mass spectra of combinations [11] of D + s π 0 , D + s γ, D + s π 0 γ, D + s π 0 π 0 , D + s γγ, and D + s π + π − particles. Combinations of D + s π + and D + s π − are also studied to search for new states. The mass spectrum of each final-state combination is studied in detail. In particular, features in the spectra that arise from reflections of other cs meson decays are individually identified and modeled. The analysis of the D + s π 0 γ final state includes an explicit search for the two most likely sub-resonant decay channels for D sJ (2460) + meson decay. This paper is organized as follows. First, the current status of the D * sJ (2317) + and D sJ (2460) + mesons is reviewed. The reconstruction of D + s , π 0 , γ, and π ± candidates is then described, including an estimate of D + s yield in terms of the D + s → φπ + branching fraction. Each combination of final-state particle species is then discussed individually. The paper finishes with a summary of results and conclusions. * Also at Laboratoire de Physique Corpusculaire, Clermont-Ferrand, France † Also with Università di Perugia, Dipartimento di Fisica, Perugia, Italy ‡ Also with Università della Basilicata, Potenza, Italy FIG. 1: The cs meson spectrum, as predicted by Godfrey and Isgur [5] (solid lines) and Di Pierro and Eichten [8] (dashed lines) and as observed by experiment (points). The DK and D * K mass thresholds are indicated by the horizontal lines spanning the width of the plot.
II. REVIEW OF THE D * sJ (2317) + AND DsJ (2460) + Much of the theoretical work on the cs system has been performed in the limit of heavy c quark mass using potential models [5][6][7][8] that treat the cs system much like a hydrogen atom. Prior to the discovery of the D * sJ (2317) + meson, such models were successful at explaining the masses of all known D and D s states and even predicting, to good accuracy, the masses of many D mesons (including the D s1 (2536) + and D s2 (2573) + ) before they were observed (see Fig. 1). Several of the predicted D s states were not confirmed experimentally, notably the lowest mass J P = 0 + state (at around 2.48 GeV/c 2 ) and the second lowest mass J P = 1 + state (at around 2.58 GeV/c 2 ). Since the predicted widths of these two states were large, they would be hard to observe, and thus the lack of experimental evidence was not a concern.
The D * sJ (2317) + meson has been observed in the decay D * sJ (2317) + → D + s π 0 [1-3, 12, 13]. The mass is measured to be around 2.32 GeV/c 2 , which is below the DK threshold. Thus, this particle is forced to decay either electromagnetically, of which there is no experimental evidence, or through the observed isospin-violating D + s π 0 strong decay. The intrinsic width is small enough that only upper limits have been measured (the best limit previous to this paper being Γ < 4.6 MeV at 95% CL as established by BELLE [3]). If the D * sJ (2317) + is the missing 0 + cs meson state, the narrow width could be explained by the lack of an isospin-conserving strong decay channel. The low mass (160 MeV/c 2 below expectations) is more surprising and has led to the speculation that the D * sJ (2317) + does not belong to the D + s meson family at all but is instead some type of exotic particle, such as a four-quark state [9].
The D sJ (2460) + meson has been observed decaying to D + s π 0 γ [2-4, 12, 13], D + s π + π − [3], and D + s γ [3,12,13]. The intrinsic width is small enough that only upper limits have been measured (the best limit previous to this paper being Γ < 5.5 MeV at 95% CL as established by BELLE [3]). The D + s γ decay implies a spin of at least one, and so it is natural to assume that the D sJ (2460) + is the missing 1 + cs meson state. Like the D * sJ (2317) + , the D sJ (2460) + is substantially lower in mass than predicted for the normal cs meson. This suggests that a similar mechanism is deflating the masses of both mesons, or that both the states belong to the same family of exotic particles.
The spin-parity of the D * sJ (2317) + and D sJ (2460) + mesons has not been firmly established. The decay mode of the D * sJ (2317) + alone implies a spin-parity assignment from the natural J P series {0 + , 1 − , 2 + , . . . }, assuming parity conservation. Because of the low mass, the assignment J P = 0 + seems most reasonable, although experimental data have not ruled out higher spin. It is not clear whether electromagnetic decays such as D * s (2112) + γ can compete with the strong decay to D + s π 0 , even with isospin violation. Thus, the absence of experimental evidence for radiative decays such as D * sJ (2317) + → D * s (2112) + γ is not conclusive. Experimental evidence for the spin-parity of the D sJ (2460) + meson is somewhat stronger. The observation of the decay to D + s γ alone rules out J = 0. Decay distribution studies in B → D sJ (2460) + D ( * )− s [12,13] favor the assignment J = 1. Decays to either D + s π 0 , D 0 K + , or D + K 0 would be favored if they were allowed. Since these decay channels are not observed, this suggests, when combined with the other observations, the assignment J P = 1 + . In this case, the decay to D * sJ (2317) + γ is allowed, but it may be small in comparison to the D + s γ decay mode. Table I lists various possible decay channels for the D * sJ (2317) + and D sJ (2460) + mesons. Several of these decays are forbidden assuming the spin-parity assignments discussed above.

III. THE BABAR DETECTOR AND DATASET
The data used in this analysis were recorded with the BABAR detector at the PEP-II asymmetric-energy storage rings and correspond to an integrated luminosity of 232 fb −1 collected on or just below the Υ(4S) resonance.
A detailed description of the BABAR detector is presented elsewhere [14]. Charged particles are detected with a five-layer, double-sided silicon vertex I: A list of various decay channels and whether they have been seen, are allowed, or are forbidden in the decay of the D * sJ (2317) + and DsJ (2460) + mesons. The predictions assume a spin-parity assignment of J P = 0 + and 1 + , respectively.
Decay Channel D * sJ (2317) + DsJ (2460) + D + s π 0 Seen Forbidden D + s γ Forbidden Seen D + s π 0 γ (a) Allowed Allowed D * s (2112) + π 0 Forbidden Seen D * sJ (2317) + γ -Allowed D + s π 0 π 0 Forbidden Allowed D + s γγ (a) Allowed Allowed D * s (2112) + γ Allowed Allowed D + s π + π − Forbidden Seen (a) Non-resonant only tracker (SVT) and a 40-layer drift chamber (DCH) using a helium-isobutane gas mixture, placed in a 1.5-T solenoidal field produced by a superconducting magnet. The charged-particle momentum resolution is approximately (δp T /p T ) 2 = (0.0013 p T ) 2 + (0.0045) 2 , where p T is the transverse momentum in GeV/c. The SVT, with a typical single-hit resolution of 10 µm, measures the impact parameters of charged-particle tracks in both the plane transverse to the beam direction and along the beam. Charged-particle types are identified from the ionization energy loss (dE/dx) measured in the DCH and SVT, and from the Cherenkov radiation detected in a ring-imaging Cherenkov device (DIRC). Photons are detected by a CsI(Tl) electromagnetic calorimeter (EMC) with an energy resolution σ(E)/E = 0.023 · (E/ GeV) −1/4 ⊕ 0.019. The return yoke of the superconducting coil is instrumented with resistive plate chambers (IFR) for the identification of muons and the detection of neutral hadrons.

IV. CANDIDATE RECONSTRUCTION
The goal of this analysis is to study the possible decay of the D * sJ (2317) + and D sJ (2460) + mesons into the final states listed in Table I. These decay channels consist of one D + s meson combined with up to two additional particles selected from π 0 , π ± , and γ. The first step in this analysis is to identify D + s mesons in the BABAR data. For each resulting D + s candidate a search is performed for associated π 0 , γ, and π ± particles. Signals from D * sJ (2317) + and D sJ (2460) + decay are isolated using the invariant mass of the desired combination of particle species.
The D + s → K + K − π + decay mode is used to select a high-statistics sample of D + s meson candidates. Each K ± and π ± candidate is separated from other charged particle species by a likelihood-based particle identification algorithm based on the Cherenkov-photon information from the DIRC together with dE/dx measurements from s candidates before (top histogram) and after (bottom points) applying subresonant φπ + and K * K + selection. The light (dark) areas indicate the signal (sideband) regions. (b) The invariant γγ mass of π 0 candidates before (top histogram) and after (bottom points) applying the π 0 veto described in the text. The light histogram indicates those candidates that pass the χ 2 requirement. The curve in (a) is the χ 2 fit described in the text.
the SVT and DCH. A geometrical fit to a common vertex is applied to each K + K − π + combination. An acceptable K + K − π + candidate must have a fit probability greater than 0.1% and a trajectory consistent with originating from the e + e − luminous region. To reduce combinatorial background, each K + K − π + candidate must have a momentum p * in the e + e − center-of-mass frame greater than 2.2 GeV/c 2 , a requirement that also removes nearly all contributions from B-meson decay. Background from D 0 → K + K − , which is evident from the corresponding K + K − mass distribution, is removed by requiring that the K + K − mass be less than 1.84 GeV/c 2 .
The upper histogram in Fig. 2(a) shows the K + K − π + mass distribution for all candidates. A clear D + s signal is seen. To reduce the background further, only those candidates with K + K − mass within 10 MeV/c 2 of the φ(1020) mass or with K − π + mass within 50 MeV/c 2 of the K * (892) mass are retained; these densely populated regions in the D + s Dalitz plot do not overlap (see Fig. 3). The decay products of the vector particles φ(1020) and K * (892) exhibit the expected cos 2 θ h behavior required by conservation of angular momentum, where θ h is the helicity angle. The signal-to-background ratio is further improved by requiring | cos θ h | > 0.5. The lower histogram of Fig. 2(a) shows the net effect of these additional selection criteria. The D + s signal (1.954 < m(K + K − π + ) < 1.981 GeV/c 2 ) and sideband (1.912 < m(K + K − π + ) < 1.934 GeV/c 2 and 1.998 < m(K + K − π + ) < 2.020 GeV/c 2 ) regions are shaded. This distribution can be reasonably modeled in a χ 2 fit by the sum of two Gaussian distributions with a common mean (hereafter referred to as a double Gaussian) on top of a quadratic background. The result of this fit is a D + s signal peak consisting of approximately 410 000 decays and a mass of 1967.8 MeV/c 2 with negligible statistical error. s sub-resonant selection requirements. The light shaded area is the kinematic range for D + s → K + K − π + decay. The vertical line with an arrow represents the selection requirement used to remove D 0 → K + K − decay. The four dark regions indicate those portions of the D + s phase space used for final candidate selection, corresponding to D + s → φπ + and D + s → K * K + decay.
FIG. 4: The sideband subtracted p * spectrum for D + s candidates after all selection requirements are met.
The approximate p * distribution for selected D + s mesons can be obtained by simple sideband subtraction, assuming linear background behavior under the D + s . The result is shown in Fig. 4.
The final list of D + s candidates are those that lie within the signal window. For each such candidate, the momentum vector is calculated from the simple addition of K + , K − , and π + momentum vectors. The energy is chosen to reproduce the PDG value for the D + s mass (1968.5 ± 0.5) MeV/c 2 [15].
It is assumed that the D * sJ (2317) + and D sJ (2460) + mesons have lifetimes that are too small to be resolved by the detector. Thus, the most likely point of D * sJ (2317) +

FIG. 5:
The D + s -sideband subtracted p * spectrum for associated (a) π ± , (b) π 0 , and (c) γ, after the selection requirements described in the text are fulfilled. and D sJ (2460) + decay for each D + s candidate is chosen to be the interaction point (IP), calculated from the intersection of the trajectory of the candidate and the e + e − luminous region. To produce a list of π ± particles that could arise from D * sJ (2317) + or D sJ (2460) + decay, the trajectories of all π ± candidates that are not daughters of the D + s candidate are constrained to the IP using a geometric vertex fit. The approximate p * spectrum of these candidates associated with real D + s mesons can be obtained by using simple D + s sideband subtraction. The result is shown in Fig. 5(a).
The selection of γ and π 0 candidates is a two-step process. The first step is the selection of a fiducial list of γ and π 0 candidates. The fiducial list of γ candidates is constructed from energy clusters in the EMC with energies above 100 MeV and not associated with a charged track. The energy centroid in the EMC combined with the IP position is used to calculate the γ momentum direction. Each fiducial π 0 candidate is constructed from a pair of γ particles in the γ fiducial list. This γ pair is combined using a kinematic fit assuming a π 0 mass. The resulting π 0 momentum is required to be greater than 150 MeV/c. The γγ invariant mass spectrum from this π 0 selection is shown in the top histogram of Fig. 2(b). To produce the final fiducial list of π 0 candidates, the χ 2 probability of the kinematic fit is required to be greater than 2%.
The final list of γ candidates consists of any γ in the fiducial list that is not used in the construction of any π 0 in the π 0 fiducial list. The final list of π 0 candidates consists of any π 0 in the fiducial list that does not share a γ with any other π 0 in the fiducial list. The γγ invari-ant mass distribution of the final list of π 0 candidates is shown in the bottom histograms of Fig. 2(b), both before and after applying the χ 2 probability requirement. The approximate p * spectrum of the final list of γ and π 0 candidates as determined using D + s sideband subtraction is shown in Figs. 5(b) and 5(c).

V. MONTE CARLO SIMULATION
Monte Carlo (MC) simulation is used for the following purposes in this paper: • To calculate signal efficiencies.
• To provide independent estimates of background levels.
• To characterize the reconstructed mass distribution of the signal.
• To predict the behavior of various specific types of backgrounds (commonly referred to as reflections) produced when the mass distribution from an established decay mode is distorted by the loss of one or more final-state particles or the addition of one or more unassociated particles.
Various sets of MC events were generated. For the purposes of understanding signal efficiencies, signal shapes, and reflections, individual MC sets of 500 000 decays were generated for each known decay mode of the D * sJ (2317) + and D sJ (2460) + mesons. In addition, MC sets of 250 000 decays were produced for each hypothetical D * sJ (2317) + and D sJ (2460) + decay as needed. Finally, a set of e + e − → cc events, corresponding to an integrated luminosity of approximately 80 fb −1 , was generated to study sources of combinatorial background.
Each MC set was processed by the same reconstruction and selection algorithms used for the data. Independent tests of the detector simulation have demonstrated an accurate reproduction of charged particle detection efficiency. The systematic uncertainty from these tests is estimated to be 1.3% for each charged track. Since the decay D + s → K + K − π + involves three charged tracks, the systematic uncertainty in D + s efficiency from the simulation alone is estimated to be 3.9%.
The simulation of D + s → K + K − π + decay was designed to match approximately the known Dalitz structure. MC events were reweighted to match more precisely the relative φπ + and K * K + yields observed in the data.
The simulation assumes a D * sJ (2317) + and D sJ (2460) + intrinsic width of Γ = 0.1 MeV. All D * sJ (2317) + and D sJ (2460) + final states are generated using phase space. The generated p * distribution of D * sJ (2317) + and D sJ (2460) + mesons in the MC simulation was adjusted to roughly reproduce observations. VI. ABSOLUTE D + s YIELD In order to calculate D * sJ (2317) + and D sJ (2460) + production cross sections, it is necessary to provide an estimate of absolute D + s selection efficiency. Since this analysis uses D + s → K + K − π + decay, the approach is to normalize D + s yield with respect to the D + s → φπ + , φ → K + K − branching fraction, the world average of which is 1.8 ± 0.4% [15]. To perform this normalization correctly, the following must be accounted for: • The K * K + portion of the D + s sample. • Non-resonant K + K − π + background under the φ peak.
• The fraction of the φ signal that falls outside of the K + K − mass selection and θ h requirements.
The K * K + selection represents approximately 48% of the total D + s sample. An inspection of the p * distribution of the φπ + and K * K + subsamples indicates that this fraction is, to a good approximation, independent of p * . Therefore, a constant factor is sufficient to account for the contribution from the K * K + portions of the D + s sample.
Shown in Fig. 6 is the D + s -sideband subtracted K + K − invariant mass spectrum for all D + s candidates before applying the φπ + and K * K + selection requirements. A prominent φ peak is observed. A binned χ 2 fit to this spectrum is used to extract both the fraction of φπ + decays that fall outside the φ selection window and the number of non-φ decays that leak inside. This fit is described below.
To model the signal portion of the K + K − mass spectrum, the φ → K + K − line shape σ(m) can be reasonably well described (ignoring potential interference effects) by a relativistic Breit-Wigner function: where m 0 = (1019.456 ± 0.020) MeV/c 2 is the intrinsic φ mass [15]. The mass-dependent width Γ(m) can be approximated by: where Γ 0 = 4.26 MeV is the intrinsic width, 0.493 is the branching fraction for this decay mode, R = 3 GeV −1 is an effective φ radius (to control the tails), and q (q 0 ) is the total three-momentum of the K ± decay products in the φ center-of-mass frame assuming an effective φ mass of m (m 0 ): A reasonable approximation of the total width includes the three dominant decay modes (K + K − , K S K L , FIG. 6: The D + s -sideband subtracted K + K − invariant mass spectrum near the φ mass for the D + s sample obtained before applying the φπ + and K * K + selection requirements. The curve is the fit described in the text. The dashed line is the portion of the fit attributed to contributions from other than φ decay. The vertical lines indicate the φ mass selection window. where Γ ′ (m) is calculated in the same manner as Γ(m) but using the K S /K L mass and branching fraction (0.337) and the width for the three-pion decay mode (which is well above threshold) is treated as a constant. The fit function P (m) used to describe the K + K − mass spectrum of Fig. 6 is: where S φ (m) is the φ signal shape and C(m) represents non-resonant contributions. The empirical form used for C(m) is a four-parameter threshold function: The form used for S φ (m) is the Breit-Wigner function of Eq. 1 smeared by a Gaussian: where N is an overall normalization. In the fit to this function, the value of Γ 0 is kept fixed to the PDG value but m 0 and δ are allowed to vary.
Despite its simplicity, the function of Eq. 5 describes the K + K − spectrum quite well, as shown in Fig. 6. The value of m 0 produced by the fit is slightly lower [(−56 ± 6) keV/c 2 , statistical error only] than the PDG average [15]. To determine a correction factor for the D + s yield, the total φ yield (calculated from the integral of the signal line shape determined by the fit up to a K + K − mass of 1.1 GeV/c 2 ) can be compared to the total number of candidates which fall inside the φ mass window. The result is a correction factor of 1.09, with negligible statistical uncertainty.
To test the above calculation, the fit is repeated on a D + s sample that includes the | cos θ h | > 0.5 requirement discussed earlier for the φ. The change in φ integrated yield is consistent with a cos 2 θ h distribution. The measured correction factor increases to 1.10. The difference between this value and 1.09 is taken as a systematic uncertainty.
Other systematic checks performed include increasing the range in K + K − mass of the φ line shape integration and changing the value of R from 1 to 5 GeV/c 2 . The total uncertainty in the 1.09 correction factor is found to be a 0.043 (3.9% relative), as calculated from a quadrature sum.

VII. CANDIDATE SELECTION OPTIMIZATION
This paper explores the eight final-state combinations shown in Table II, each involving a D + s meson and up to two total of π ± , π 0 , and/or γ particles. For each combination it is necessary to distinguish possible signals from D * sJ (2317) + and D sJ (2460) + decay from combinatorial background. The separation of signal and background is made more distinct if additional candidate selection requirements are imposed. This section discusses those additional requirements.
In four of the final-state combinations (D + s π 0 , D + s π 0 γ, D + s γ, and D + s π + π − ) a signal is expected. An estimate of signal significance is calculated in these cases based on expected signal and background rates, the former calculated from previously published branching ratio measurements combined with the appropriate MC sample. For the remainder of the final states, an estimate of signal sensitivity is calculated for the hypothetical D * sJ (2317) + and D sJ (2460) + meson decay. This sensitivity calculation is based on signal efficiency, determined using MC samples, and expected background levels.
To avoid potential biases, both the signal significance and sensitivity estimates are calculated solely using MC samples.
The following selection requirements are adjusted in order to produce optimal values of signal significance and sensitivity: • A minimum total center-of-mass momentum p * .
• A minimum energy (for γ) and/or momentum calculated in the laboratory (for π 0 ) or center-of-mass (for π ± ) frame of reference.
For the minimum p * , it is important to choose the same value for all final-state combinations in order to minimize systematic uncertainties in the branching ratios. A minimum value of p * > 3.2 GeV/c is chosen as a reasonable compromise. Values for the remaining selection requirements are chosen separately for each final-state combination. The results are listed in Table II. The MC samples can be used to estimate the approximate efficiency for detecting a signal with a p * of at least 3.2 GeV/c after applying the above selection requirements. The resulting efficiencies vary between 0.4% and 13% (see Table II).

VIII. CROSS SECTION NOTATION
To report production yields of a particular cs meson D Y to a particular final state D + s X, the following quan- where the cross section σ is defined for a center-of-mass momentum p * above 3.2 GeV/c. The quantityσ is calculated by taking the number of D Y → D + s X decays observed in the data, correcting for efficiency using the appropriate MC sample (restricted to p * > 3.2 GeV/c), correcting for the relative D + s → φπ + yield as calculated in Section VI, and dividing by the luminosity (232 fb −1 ). A relative systematic uncertainty of 1.2% is introduced to account for the uncertainty in the absolute luminosity.
There is no attempt to correct the cross section for radiative effects (such as initial-state radiation). Since a reasonably accurate representation of such radiative effects is included in our MC samples, the calculation of se-lection efficiencies from these samples is accurate enough for the purposes of this paper. Fig. 7 is the invariant mass distribution of the D + s π 0 combinations after all selection requirements are fulfilled. Signals from D * s (2112) + and D * sJ (2317) + decay are evident. An unbinned likelihood fit is applied to this mass distribution in order to extract the parameters and yield of the D * sJ (2317) + signal and upper limits on D sJ (2460) + decay. The likelihood fit includes six distinct sources of D + s π 0 combinations: • A reflection from D * s (2112) + → D + s γ decay in which an unassociated γ particle is added to form a false π 0 candidate.
• Combinatorial background from unassociated D + s and π 0 mesons. The probability density function (PDF) used to describe the mass distribution of each of these sources is described below.
As shown in Fig. 8(a), the reconstructed mass distribution of the D * sJ (2317) + → D + s π 0 decay, as predicted by MC, has non-Gaussian tails and is slightly asymmetric. To describe this shape, the MC sample is fit to a modified Lorentzian function F L (m): where a 1 and a 2 correspond roughly to a mean and width, respectively. This function is simply a convenient parameterization of detector resolution. The fit results are shown in Fig. 8(a). A similar procedure is used for the hypothetical D sJ (2460) + → D + s π 0 decay ( Fig. 8(b)). The 350 MeV/c π 0 momentum requirement removes the majority of D * s (2112) + → D + s π 0 decays. The remaining signal is modeled using a distribution J(m) of Gaussians: The parameters a 3 and a 4 of this function are determined using a fit to a suitable MC sample. The mean mass a 2 is set equal to 2112.9 MeV/c 2 (0.5 MeV/c 2 higher than the PDG value [15]) to match the data. A reflection in D + s π 0 produced by D * s (2112) + → D + s γ decay appears as a broad distribution peaking at a mass of approximately 2.17 GeV/c 2 . This reflection is produced by fake π 0 candidates consisting of the γ particle from D * s (2112) + decay combined with unassociated γ candidates. Kinematics limit this reflection to D + s π 0 masses above the quadrature sum of D + s and π 0 meson masses (approximately 2.1167 GeV/c 2 ). This distribution falls gradually as mass is increased due to the rapidly falling inclusive γ energy spectrum. Detector resolution tends to smear the lower kinematic mass limit. To model this reflection, a quadratic function with a sharp lower mass cut off is convoluted with a Gaussian distribution. The parameters of this function are determined directly from the D + s π 0 data sample. The D sJ (2460) + reflection requires careful attention because it appears directly under the D * sJ (2317) + signal. This reflection is produced by the D + s π 0 projection of D sJ (2460) + → D * s (2112) + π 0 decay (in which the γ from D * s (2112) + decay is ignored). A kinematic calculation ( Fig. 8(c)) of the D sJ (2460) + Dalitz distribution predicts that this reflection, at the limit of perfect resolution and efficiency, is a flat distribution in mass squared centered at a D + s π 0 mass of 2313.4 MeV/c 2 with a full width of 41.3 MeV/c 2 (assuming a D sJ (2460) + mass of 2458.0 MeV/c 2 ). The D sJ (2460) + reflection is flat in mass squared only if the π 0 efficiency is constant. In practice this is not the case, as illustrated by MC simulation (Fig. 8(d)). To accommodate the non-constant efficiency, the D + s π 0 mass distribution from the MC sample is fit to a function consisting of a bounded quadratic function smeared by a double Gaussian. The result of the fit is shown in Fig. 8d.
The threshold function C(m) of Eq. 6 is used to represent the mass spectrum from combinatorial background where the threshold value a 1 is fixed to 2103.5 MeV/c 2 , the sum of the assumed D + s and π 0 masses. The remaining parameters of C(m) are determined directly from the data.
The results of the likelihood fit to the D + s π 0 mass spectrum is shown in Fig. 7. In this fit, the size, shape, and mean mass of the D sJ (2460) + reflection are fixed to values consistent with the yield and mass results determined in Section XI of this paper. The yield of D * s (2112) + and D * sJ (2317) + decay and the D * sJ (2317) + mass is allowed to vary to best match the data. A D * sJ (2317) + mass of (2319.6 ± 0.2) MeV/c 2 is obtained (statistical error only). A total of 3180 ± 80 D * sJ (2317) + decays are found. The fit includes a hypothetical contribution from D sJ (2460) + → D + s π 0 in the form of a line shape of fixed shape and mass. The result is a yield of −40 ± 50 (statistical errors only). The size of this yield is small enough that the curve cannot be distinguished in Fig. 7.
The D sJ (2460) + reflection arises from contamination from D * s (2112) + decay. It is an interesting exercise to identify and separate some of this background. This can be accomplished by searching for any γ candidates that, when combined with the D + s in the same event, produce a D + s γ mass within 15 MeV/c 2 of the D * s (2112) + mass. Those D + s π 0 combinations in which such a match is not found will contain a smaller proportion of D sJ (2460) + reflection, whereas the remaining D + s π 0 combinations will contain fewer D * sJ (2317) + decays. This is indeed the case as illustrated in Fig. 9. The same likelihood fit procedure used for the entire sample is repeated for these D + s π 0 subsamples, including the MC prediction of the yield and shape of the D sJ (2460) + reflection. The fit results are consistent with the data.
As can be seen in Figs. 7 and 9, the D * sJ (2317) + line shape derived from MC simulation and used unchanged in the likelihood describes the data well. Since the MC simulation is configured with an intrinsic width (0.1 MeV) nearly indistinguishable from zero, it follows that the data are consistent with a zero width D * sJ (2317) + meson.
To extract a 95% CL upper limit on the intrinsic D * sJ (2317) + width, the D * sJ (2317) + line shape is convolved with a relativistic Breit-Wigner function σ J (m) with constant width Γ: The fit is then repeated in its entirety at incremental steps in Γ to produce a likelihood curve. Integrating this curve as a function of Γ produces a 95% CL upper limit of Γ < 1.9 MeV (statistical error only).
In order to produce an estimate of D * sJ (2317) + yield, the fit results must be corrected for selection efficiency. This efficiency is calculated using a D * sJ (2317) + MC sample and is p * dependent. Since the p * distribution observed in data does not exactly match the MC simulation, it is important to take into account this p * dependence. Two methods are used to do this. The first is to weight each D + s π 0 combination by the inverse of the selection efficiency before applying a likelihood fit. After correcting for absolute D + s → φπ + yield (Section VI), the result is: The second method is to divide the D + s π 0 sample into bins of p * . A likelihood fit is applied to each bin and the yield corrected for the average selection efficiency in that bin. The result is the p * distribution shown in Fig. 10. The total yield from this method is 26470 ± 660 (statistical error only).
The systematic uncertainties for the D * sJ (2317) + mass and yield are summarized in Table III. The uncertainties in D sJ (2460) + yield are calculated in the same fashion. The assumed D + s mass value and the 1% relative uncertainty in the EMC energy scale are the two largest contributors to the error on the D * sJ (2317) + mass. Uncertainties in the signal shape produce the largest uncertainties in D * sJ (2317) + yield and width. For example, the amount of D sJ (2460) + reflection is proportional to the D sJ (2460) + yield, which, as will be discussed in Section XI, has an 11% uncertainty. Adjusting the contribution of this reflection in the fit with no other changes produces a relative uncertainty of 2.0% in the D * sJ (2317) + yield with little change in the D * sJ (2317) + mass and width limit. If the likelihood fit is allowed to choose a D sJ (2460) + reflection yield that best matches the data, little change in either the D * sJ (2317) + mass or yield is observed. The limit on the intrinsic width, however, increases to 3.4 MeV, since reducing the D sJ (2460) + reflection allows the observed mass line shape to accommodate a larger intrinsic D * sJ (2317) + width.
Another uncertainty that has a similar effect is the assumed D + s π 0 mass resolution. Small variations in resolution, consistent with comparisons of data and MC simulation of other known particles, can change the D * sJ (2317) + line shape sufficiently to lower or raise yields by 3.2%. Allowing better reconstructed resolution provides more room for a large intrinsic width, raising the 95% CL for Γ to 3.0 MeV.
Other uncertainties in the D * sJ (2317) + yield include the accuracy (±3%) of the MC prediction of π 0 efficiency, the difference of the two methods for correcting for p *dependent efficiency, and the D + s → φπ + branching fraction. The total systematic uncertainty is calculated from the quadrature sum of all sources. The result is the following D * sJ (2317) + mass: m = (2319.6 ± 0.2 (stat.) ± 1.4 (syst.)) MeV/c 2 , and the following yields: where the first error is statistical and the second is systematic.
As can be seen in Table III, the determination of the D * sJ (2317) + mass is limited by the understanding of the EMC energy scale. A more primitive calculation of this energy scale was used in a previous estimate of the D * sJ (2317) + mass from this collaboration [4], resulting in an mass estimate that is 2.3 MeV/c 2 lighter then the estimate presented here. The associated systematic uncertainty in this previous work was also incorrectly calculated. The central value and systematic uncertainty in mass reported here reflects the current best understanding of these calibration issues.
For the sake of simplicity, in order to incorporate systematic effects into a limit on Γ, the least strict limit obtained from the various systematic checks is quoted. This produces a 95% CL of Γ < 3.8 MeV. The lineshape produced by this limit is illustrated in Fig. 11.
The D + s γ mass distribution using γ candidates with loose (150 MeV) and final (500 MeV) minimum energy requirements is shown in Fig. 12. Some structure in the vicinity of the D * sJ (2317) + and D sJ (2460) + masses becomes apparent once the tighter energy requirements are applied. The looser requirement is useful for studying the D * s (2112) + peak. A fit to that peak consisting of two Gaussians on top of a polynomial background function results in a peak D * s (2112) + mass of 2113.8 MeV/c 2 and a yield of 75 000 decays, both with negligible statistical uncertainties.
An unbinned likelihood fit is applied to the final D + s γ mass distribution in order to extract the parameters and yield of the D sJ (2460) + meson and upper limits on D * sJ (2317) + decay. For simplicity, the fit is performed only for masses between 2.15 and 2.85 GeV/c 2 . The likelihood fit includes five distinct sources of D + s γ combinations: • A reflection from D * sJ (2317) + → D + s π 0 decay in which only one of the γ particles from π 0 decay is included.
• A reflection from D sJ (2460) + → D * s (2112) + π 0 decay in which only one of the γ particles from π 0 decay is included.
• Background from both unassociated D + s and γ mesons and the high-mass tail from D * s (2112) + → D + s γ decay. The PDF used to describe the mass distribution of each of these sources is described below.
As shown in Fig. 13(b), the reconstructed mass distribution of the D sJ (2460) + → D + s γ decay, as predicted by MC, has a long, low mass tail. To describe this shape, the MC sample is fit to a modified Lorentzian function F L2 (m): The fit results are shown in Fig. 13(b). A similar procedure is used for the hypothetical D * sJ (2317) + → D + s γ decay ( Fig. 13(a)).
Ignoring resolution effects, the D * sJ (2317) + reflection produces an invariant D + s γ mass distribution up to a maximum of approximately 140 MeV/c 2 below the D * sJ (2317) + mass. Candidate selection requirements produce a distribution that peaks at this limit. Resolution effects smear this sharp peak producing the shape shown in Fig. 13(c), as predicted by MC. This distribution can be reasonably described (in the mass range of interest) by a bounded quadratic function convoluted with a double Gaussian. The D sJ (2460) + reflection has a similar behavior near the D * sJ (2317) + mass ( Fig. 13(d)). Both distributions overlap the direct D * sJ (2317) + → D + s γ decay.
The following function D(m) is used to represent the This includes combinatorial background along with any tail from D * s (2112) + → D + s γ decay. The MC simulation fails to reproduce the shape of this background, either due to unknown cs contributions (for example, higher resonant states) or unexpected behavior of the tail of the distribution from D * s (2112) + decay. This issue, combined with the complex shapes associated with the D * sJ (2317) + and D sJ (2460) + reflections, leads to considerable systematic uncertainty in the fit. Likelihood fits under several different conditions are attempted in order to understand the uncertainty.
One fit that produces a good representation of the data is shown in Fig. 14. In this fit, all parameters except the upper mass limit of the D * sJ (2317) + reflection are allowed to vary. The estimated raw D sJ (2460) + (D * sJ (2317) + ) yield from this fit is 920 ± 60 (−130 ± 130). The fitted D sJ (2460) + mass is (2459.5 ± 1.2) MeV/c 2 (statistical errors only).
The D sJ (2460) + signal is far enough removed from the D * sJ (2317) + and D sJ (2460) + reflections that accurate mass and yield results are obtained. The same two methods described in the previous section are used to estimate the D sJ (2460) + yield. The first method, using p * dependent weights proportional to the inverse of efficiency, produces a corrected yield of: for p * > 3.2 GeV/c (statistical error only). The second method produces the p * spectrum shown in Fig. 15. The total yield from the spectrum is 3080 ± 240 (statistical error only), approximately 6.2% lower than the first estimate.
Systematic uncertainties associated with the assumed PDFs for the background and reflection shapes are explored using different variations of the likelihood fit. Among the variations applied is an alternate description In addition, the MC predictions for the size and shape of the D * sJ (2317) + reflection are used unaltered (despite producing a fit of inferior quality). Large variations of raw D * sJ (2317) + yield of up to 490 events are observed. The shape of the D sJ (2460) + → D + s γ signal is sensitive to several factors that are difficult to simulate exactly, including EMC energy resolution. Variations in the assumed resolution are used to study the associated systematic uncertainty in yield and mass. The result is an uncertainty of 3.5% in yield and no significant change in mass.
All systematic uncertainties for the D sJ (2460) + mass and yield are listed in Table IV. The result is the following D sJ (2460) + mass: m = (2459.5 ± 1.2 (stat.) ± 3.7 (syst.)) MeV/c 2 , and the following yields: where the first error is statistical and the second is systematic.
The invariant mass spectrum for all selected D + s π 0 γ candidates is shown in Fig. 16(a). A D sJ (2460) + signal is apparent. The shape of this signal is characterized by applying the following modified Lorentzian fit function F L3 : δ ≡ (m − a 1 ) /a 2 . to a D sJ (2460) + → D * s (2112) + π 0 MC sample ( Fig. 16(b)). A binned χ 2 fit to the D + s π 0 γ spectrum that includes this shape along with a polynomial description of the background produces a D sJ (2460) + mass (2459.5 ± 2.0) MeV/c 2 and a yield of 560 ± 80 events (statistical errors only).
In the following, it is assumed that the D sJ (2460) + meson decays to D + s π 0 γ entirely through either of the two kinematically allowed sub-resonant decay modes: Due to a kinematic accident, the phase space of these two sub-resonant modes overlap, as illustrated in Fig. 17. It is therefore possible to remove background while retaining both sub-resonant decay modes by selecting D + s π 0 γ candidates in either a restricted range of D + s π 0 mass or a restricted range of D + s γ mass. This analysis uses a requirement that the D + s γ mass must be within 20 MeV/c 2 of the D * s (2112) + mass. As shown in Fig. 16(a), the resulting D + s π 0 γ mass distribution in this D + s γ signal window is considerably cleaner.
The D + s γ requirement introduces a source of background that peaks underneath the D sJ (2460) + signal. This background is a reflection from D * sJ (2317) + → D + s π 0 decays that are not associated with any D sJ (2460) + decay. The reflection arises because, as illustrated in Fig. 17, any D * sJ (2317) + signal that is combined with a γ candidate that produces a D + s γ mass near the D * s (2112) + meson results in a D + s π 0 γ mass near the D sJ (2460) + meson. If the D + s γ mass requirement is shifted upwards or downwards, this D * sJ (2317) + reflection shifts up and down in D + s π 0 γ mass by a predictable amount.
Another type of background that behaves similarly to the D * sJ (2317) + reflection is D sJ (2460) + → D + s π 0 γ decay in which the wrong γ candidate is chosen. This type of background is slightly wider and smaller than the D * sJ (2317) + reflection, but otherwise has a similar mass distribution. To describe this background contribution, it is assumed that the D sJ (2460) + decays entirely through D * s (2112) + π 0 and is produced at a rate comparable to previous measurements. Both assumptions need not be entirely accurate since this background has a relatively small contribution.
To characterize the D * sJ (2317) + and D sJ (2460) + reflections, upper and lower D + s γ mass selection windows are chosen centered at ±60 MeV/c 2 away from the D * s (2112) + mass. The mass distribution from MC samples of D * sJ (2317) + → D + s π 0 and D sJ (2460) + → D * s (2112) + γ decay are shown in Fig. 18 for the signal and two sideband D + s γ mass windows. The shape of the two combined reflections in all three cases can be successfully described by a fit to a Gaussian.
To determine the mass, width, and yield of the D sJ (2460) + meson, an unbinned likelihood fit is applied to the D + s π 0 γ mass distribution of candidates selected in the D + s γ signal window. This fit includes the following contributions: • D sJ (2460) + → D + s π 0 γ decay.
• A reflection from D * s (2112) + → D + s γ decay in which an unassociated γ candidate is added to form a fake π 0 candidate.
• Smooth background sources that do not have any peaking behavior.
The D * s (2112) + → D + s γ reflection is similar to that observed in D + s π 0 combinations (see, for example, Fig. 7). The smooth background is represented by the C(m) function described in Eq. 6.
Two similar fits excluding the D sJ (2460) + signal are applied to the upper and lower D + s γ mass samples. These fits suggest that the MC prediction of the absolute rate of the D * s (2112) + → D + s γ reflection is approximately 21% too low. The fit models are adjusted accordingly. The fit results for all three D + s γ mass ranges after this adjustment are shown in Fig. 19. The result is a D sJ (2460) + mass of (2458.6 ± 1.3) MeV/c 2 and a raw yield of 560 ± 40 events (statistical errors only).
Although the overall size of the D * sJ (2317) + → D + s π 0 reflection is allowed to vary, the MC prediction for the relative contributions in the three D + s γ mass windows is preserved. Since the size of this reflection is adequately modeled in the two D + s γ sidebands, there is some confidence that the size is well established in the D + s γ signal window.
Note that the size of the smooth background is relatively larger in the signal D + s γ window due to contributions from D * s (2112) + → D + s γ decay. In addition, if there were significant non-resonant contributions to D sJ (2460) + decay, peaks at the D sJ (2460) + mass would be visible in the two D + s γ sidebands. No such evidence is visible.
The two methods described in the two previous sections are used to estimate the D sJ (2460) + yield. The first method, using p * dependent weights proportional to the inverse of efficiency, produces a corrected yield of: method produces the p * spectrum shown in Fig. 20. The total yield from the spectrum is 9890 ± 810 (statistical error only). The systematic uncertainties in the mass and yield of the D sJ (2460) + meson are shown in Table V. As described previously, the size of the D * sJ (2317) + reflection was adjusted in the fit to match the D + s γ sideband samples. If the size of the reflection is taken unchanged from MC predictions, the D sJ (2460) + yield increases by 7.3%. A second likelihood fit described later in this section used to distinguish between the two D sJ (2460) + sub-resonant decay modes also produces an estimate of D sJ (2460) + yield. The difference between the two fits is treated as a systematic uncertainty. The total D sJ (2460) + yield, without distinguishing between the two possible subresonant decay modes, is measured to bē The D sJ (2460) + signal PDF used in the likelihood fit of Fig. 19 includes a D sJ (2460) + intrinsic width of   Table V, the result is a 95% CL limit of Γ < 6.3 MeV.
Having established a D sJ (2460) + → D + s π 0 γ signal, it is now necessary to distinguish between the two possible sub-resonant decay modes shown in Eq. 16. These two decay modes can be distinguished by their D + s π 0 and D + s γ invariant mass distributions, as shown in Fig. 21. The distributions for the D sJ (2460) + → D * s (2112) + π 0 subresonant mode are determined using a MC sample. The reconstructed D + s γ mass distribution, which is relatively narrow (as it arises from D * s (2112) + decay), is represented by a χ 2 fit to the F L3 function (Eq. 15). The wider D + s π 0 mass distribution is accurately modeled by In contrast, for the D sJ (2460) + → D * sJ (2317) + γ subresonant decay mode, the D + s π 0 mass distribution is narrow and the D + s γ mass distribution is wide. The D + s π 0 distribution is determined using a D * sJ (2317) + → D + s π 0 MC sample. The D + s γ distribution is calculated using the parameters determined from the D + s π 0 distribution from D sJ (2460) + → D * s (2112) + π 0 described above converted to the appropriate kinematic range. The shapes assumed for both D sJ (2460) + → D * sJ (2317) + γ mass projections are shown in gray in Fig. 21.
The D + s π 0 and D + s γ mass distributions of the signal cannot be explored without correctly subtracting backgrounds from unassociated D * s (2112) + → D + s γ and D * sJ (2317) + → D + s π 0 decay. This subtraction is performed by a two-dimensional unbinned likelihood fit applied to the D + s π 0 and D + s γ mass distributions of the data. The likelihood fit is restricted to the data sample contained inside the grid shown in Fig. 22. This fit includes five sources of D + s π 0 γ candidates: • Combinatorial background represented by a twodimensional quadratic function.
• Background from D * s (2112) + → D + s γ decay combined with unassociated π 0 candidates represented by a D * s (2112) + line shape in the D + s γ mass and as a linear function in D + s π 0 mass. • Background from D * sJ (2317) + → D + s π 0 decay combined with unassociated γ candidates represented by a D * sJ (2317) + line shape in the D + s π 0 mass and as a linear function in D + s γ mass. • A signal from D sJ (2460) + → D * s (2112) + π 0 with D + s π 0 and D + s γ mass distributions represented by the curves in Fig. 21.
• A signal from D sJ (2460) + → D * sJ (2317) + γ with D + s π 0 and D + s γ mass distributions represented by the gray regions in Fig. 21.  Fig. 23, divided into the regions delineated by the grid shown in Fig. 22. The fit produces an adequate model of the data in all regions. The result (statistical errors only) is a total yield of 520 ± 50 D sJ (2460) + → D + s π 0 γ decays with a fraction of (2.5 ± 8.8)% proceeding through the D * sJ (2317) + γ channel, the former number being somewhat smaller than the yield determined by the D + s π 0 γ mass fit (Fig. 19), though consistent within systematic uncertainties. Based on these results, it appears that the decay D sJ (2460) + → D + s π 0 γ can be successfully described as proceeding entirely through the channel D * s (2112) + π 0 . The systematic uncertainties listed in Table V can be applied to the above fit results. In combination with the results of the fit to the D + s π 0 γ mass distribution of Fig. 19, and treating correlated systematic uncertainties in the appropriate fashion, the following yields for the subresonant specific decays are obtained:

The result of this likelihood fit is shown in
where the first error is statistical and the second is systematic.
A simple helicity analysis is used to test the J P assignment of the D sJ (2460) + meson under the assumption that the decay D sJ (2460) + → D + s π 0 γ proceeds entirely through the subresonant mode D * s (2112) + π 0 . This analysis is performed in terms of the helicity angle ϑ h , defined as the angle of the γ in the D * s (2112) + centerof-mass frame with respect to the D * s (2112) + direction. Since the D * s (2112) + is a vector particle, the helicity distribution must consist of some combination of the zero helicity distribution H 0 : As listed in Table VI, the expected combination of H 0 and H 1 depends on the assumed D sJ (2460) + spin and parity.
To measure the helicity distribution, the D + s π 0 γ candidates are divided into five bins of cos ϑ h . The D + s π 0 γ mass fit of Fig. 19 is repeated for each of these subsamples using p * -dependent weights inversely proportional to the selection efficiency in order to correct for acceptance. The result is shown in Fig. 24. The integral of the following function is calculated in each cos ϑ h bin: A χ 2 fit is used to determine the most likely value of a 2 . The result is a 2 = 0.76 ± 0.14 (statistical errors only).
The same procedure can be repeated after each relevant systematic check listed in Table V is performed. The differences in a 2 values so obtained are added in quadrature to estimate the total systematic uncertainty. The final result is: where a 2 = 0 (a 2 = 1) corresponds to a D sJ (2460) + of helicity zero (±1). This value of a 2 deviates from zero by 5.1 standard deviations, which strongly disfavors the J P = 0 − interpretation of the D sJ (2460) + while remaining consistent with a J = 1 or higher interpretation of either parity XII. THE D + s π 0 π 0 FINAL STATE The D + s π 0 π 0 final state contains potential contributions from D sJ (2460) + → D * s (2112) + π 0 subresonant decay through the channel D * s (2112) + → D + s π 0 , which has a branching ratio of 5.8% [15]. Since this sub-resonant mode is more efficiently investigated using the D + s π 0 γ final state (as discussed in the previous section), it is removed from the D + s π 0 π 0 sample by requiring the D + s π 0 invariant mass to be greater than 2117.4 MeV/c 2 for both π 0 candidates. This requirement excludes the edges of the D * sJ (2317) + and D sJ (2460) + phase spaces, as illustrated in Fig. 25. This figure also demonstrates how little phase space is available to the D * sJ (2317) + meson in this decay in comparison to the D sJ (2460) + meson.
The invariant mass distribution of the selected D + s π 0 π 0 candidates is shown in Fig. 26. There is no evidence of D * sJ (2317) + or D sJ (2460) + meson decay, nor is there evidence of any structure in the background. The mass distribution is fit using a likelihood function that consists of the smooth background function C(m) (Eq. 6) and D * sJ (2317) + and D sJ (2460) + meson contributions, the latter two having yields that are allowed to fluctuate to negative values. The shape of the signals is modeled by double Gaussians, the parameters of which are determined by fits to MC samples. The masses of the D * sJ (2317) + and D sJ (2460) + mesons are fixed to values measured in the previous sections.
Various systematic uncertainties are considered. According to MC simulation, the selection efficiency varies by as much as 20% across the Dalitz plot. Since no specific Dalitz distribution is assumed, this variation is translated directly into a multiplicative systematic uncertainty. An alternate shape C ′ 2 (m) for the background is considered: A fit with this background increases the raw D * sJ (2317) + yield by 0.5 events. Detector resolution and D * sJ (2317) + and D sJ (2460) + meson mass variations have similar effects. The final results arē where the first error is statistical and the second is systematic.

XIII. THE D + s γγ FINAL STATE
It is assumed that there are three possible contributions to the decay of the D * sJ (2317) + and D sJ (2460) + mesons to the D + s γγ final state: • Two body decay to D + s π 0 followed by π 0 → γγ. • Subresonant decay to D * s (2112) + γ followed by D * s (2112) + → D + s γ. • Non-resonant decay directly to D + s γγ. The D + s π 0 final state, already studied in Section IX, is removed from the D + s γγ sample by the γ selection requirements. Potential background from π 0 mass tails is removed by further requiring the γγ invariant mass to be less than 100 MeV/c 2 or greater than 170 MeV/c 2 . The remaining two D + s γγ contributions are treated by dividing the sample into two portions, one rich in D * s (2112) + decay (referred to as the D * s (2112) + sample), and the remainder (referred to as the non-resonant sample). A D + s γγ candidate is placed in the D * s (2112) + sample if either γ produces a D + s γ invariant mass within 15 MeV/c 2 of the PDG value for the D * s (2112) + mass (2112.4 MeV/c 2 ) [15]. Figure 27 illustrates how the the D + s γγ candidates are divided in terms of the phase space of D * sJ (2317) + and D sJ (2460) + meson decay. The D + s γγ invariant mass distribution of the two samples is shown in Fig. 28. No clear D * sJ (2317) + or D sJ (2460) + signal is observed. The mass distributions contain structure associated both with the π 0 veto requirements and with the D * s (2112) + sample selection requirements. Much of this structure can be avoided by applying likelihood fits in a restricted mass range. Two types of background that produce structure that cannot be so avoided are described below. Either of the two γ candidates in each D + s γγ combination provides two opportunities for the D + s γγ combination to be placed into the D * s (2112) + sample. At a specific value of D + s γγ invariant mass, however, if one γ candidate falls inside the D * s (2112) + mass window, the other γ is likely to do the same, due to kinematics. This produces a deficit of candidates at this mass in the D * s (2112) + sample along with a corresponding excess of candidates in the non-resonant sample. The D + s γγ invariant mass distribution of this excess can be approximated by the triangular shape centered at 2255.9 MeV/c 2 shown in Fig. 29(a). The second background source, which only appears in the D * s (2112) + sample, is a reflection from D sJ (2460) + → D * s (2112) + π 0 decay in which one γ particle from the π 0 decay is ignored. This reflection is suppressed but not eliminated by the requirement that each γ candidate not belong to a fiducial π 0 candidate. The MC prediction of the shape of this reflection is shown in Fig. 29(b). The shape is accurately modeled by a quadratic polynomial bounded on two sides and smeared by a double Gaussian.
The unbinned maximum-likelihood fit to the D + s γγ mass distribution of the D * s (2112) + sample is performed between masses of 2.22 and 2.70 GeV/c 2 . This mass distribution has a cusp near a mass of 2.4 GeV/c 2 that is difficult to describe using a simple polynomial. Instead, the combinatorial background is parameterized using an empirical function K(m) composed of a line and a parabola that intersect at one point. To make the function smooth, a cubic spline K S (m) is used near the intersection point (m = a 1 ) to extrapolate between these two polynomials: where the value δ = 20 MeV/c 2 is chosen to approximately match the resolution.
The reconstructed mass distributions for the hypothetical decay of D * sJ (2317) + and D sJ (2460) + mesons to D + s γγ and D * s (2112) + γ are modeled by the functional form of Eq 12, the parameters of which are determined by fits to the corresponding MC samples.

A contribution from the D sJ (2460) +
→ D * s (2112) + π 0 reflection, of fixed shape based on the fit of Fig. 29(b) and with a yield consistent with the results of Section XI.

A background described by the function of Eq. 22.
The result of applying the fit function on the D * s (2112) + sample is shown in Fig. 30.
The proximity of the D sJ (2460) + → D * s (2112) + π 0 reflection to potential D * sJ (2317) + and D sJ (2460) + signals coupled with the unknown shape of the combinatorial background leads to large uncertainties in the fit results. For example, if a third-order polynomial is used in place of the function K(m) to describe the smooth background, the D * sJ (2317) + and D sJ (2460) + raw yields increase by 260 and 150 candidates, respectively. Although the result of this fit is a less faithful representation of the mass distribution of the data, as a conservative estimate the entire difference is quoted as a systematic uncertainty. Other systematic checks performed include a variation of the range in D + s γγ mass over which the fit is applied, consideration of uncertainties in detector resolution, and variations of the D * sJ (2317) + and D sJ (2460) + masses and the relative size of the D sJ (2460) + → D * s (2112) + π 0 reflection.
The raw D * sJ (2317) + and D sJ (2460) + yields are corrected for selection efficiency by weighting each D + s γγ combination by the inverse of the p * -dependent selection where the first error is statistical and the second is systematic. The D + s γγ mass distribution of the phase-space subsample has a cusp that is more pronounced than in the D * s (2112) + γ sample but at a mass well below the D * sJ (2317) + mass. This cusp is avoided by restricting the fit to masses above 2.2 GeV/c 2 . In this mass range a simple polynomial is sufficient to describe combinatorial background. The unbinned likelihood fit includes four components: 2. D sJ (2460) + → D + s γγ decay (hypothetical).
3. An excess of events of the shape described in Fig. 29(a) but of variable size.

Combinatorial background represented by a thirdorder polynomial.
The result of applying the fit function on the nonresonant sample is shown in Fig. 31. A raw D * sJ (2317) + (D sJ (2460) + ) yield of 190 ± 120 (80 ± 160) is obtained (statistical errors only).
As in the D * s (2112) + sample, raw D * sJ (2317) + and D sJ (2460) + yields are corrected for selection efficiency by weighting each D + s γγ combination by the inverse of the p * -dependent selection efficiency. This efficiency was determined using a MC sample that simulated the non-resonant decays of the D * sJ (2317) + and D sJ (2460) + mesons to D + s γγ using phase space (such that the Dalitz plot was evenly populated). To simulate an alternate decay model, the MC samples were weighted to form a cos 2 ϑ Dγ distribution, where ϑ Dγ is the angle between the D + s and each γ in the D * sJ (2317) + and D sJ (2460) + center-of-mass frame. This weighting had the effect of reducing the selection efficiency by 25% and 15% for the D * sJ (2317) + and D sJ (2460) + mesons, respectively. Other systematic checks include variations in detector resolution, the D * sJ (2317) + and D sJ (2460) + meson masses, and the fit range. The final results for the nonresonant D + s γγ final state are: where the first error is statistical and the second is systematic.
XIV. THE D + s π ± FINAL STATES The invariant mass distributions of the D + s π − and D + s π + candidates are shown in Fig. 32 for masses between 2.25 and 2.61 GeV/c 2 . No structure is apparent and both distributions can be adequately described by a second-order polynomial. Since no signal is apparent, the remaining task is to place limits on the yield of hypothetical doubly-charged (D * sJ (2317) ++ ) or neutral (D * sJ (2317) 0 ) partners of the D * sJ (2317) + meson. To do so requires a few assumptions: • The intrinsic width Γ of either the D * sJ (2317) 0 or D * sJ (2317) ++ particle is too small to be resolved by the detector.
• The mass of the D * sJ (2317) 0 or D * sJ (2317) ++ particle is within ±10 MeV/c 2 of the D * sJ (2317) + mass. A series of unbinned likelihood fits based on these assumptions are applied to the D + s π − and D + s π + mass distributions. Included in these fits is a D * sJ (2317) 0 or D * sJ (2317) ++ signal modeled using a line shape extracted from a fit to a D * sJ (2317) 0 → D + s π − MC sample. According to this fit, the mass resolution is approximately 1.3 MeV/c 2 . To avoid potential statistical biases, the mass of each hypothetical particle is fixed at a specific value for each fit. Several such fits are applied with the assumed mass placed between 2307.3 and 2327.3 MeV/c 2 at intervals of 1 MeV/c 2 . The fits that produce the largest yields are shown in Fig. 32.
The two final states discussed in this section involve only charged particles and are thus subject to relatively small systematic uncertainties. The largest uncertainty arises from the assumed shape of the background and is estimated by substituting a third-order polynomial for the second-order one. There is also a 1.3% relative uncertainty in the reconstruction efficiency of each π ± candidate. The results from the fits that produce the largest yields are: where the first error is statistical and the second is systematic. These results will be used to calculate upper limits on the cross section.
The invariant mass distribution of the D + s π + π − candidates is shown in Fig. 33. Clear peaks from D sJ (2460) + and D s1 (2536) + decay are apparent. The distribution has no additional structure of any statistical significance. To determine the mass and yield of the D sJ (2460) + meson, and to place limits on D * sJ (2317) + decay, an unbinned likelihood fit is applied to the D + s π + π − mass distribution. This fit includes the following contributions: • D * sJ (2317) + → D + s π + π − decay (hypothetical). • D sJ (2460) + → D + s π + π − decay. • D s1 (2536) + → D + s π + π − decay. • A third-order polynomial to describe the background.
Each signal decay mode is described by a PDF consisting of three Gaussians with a common mean, the parameters of which are determined using fits to MC samples (generated with Γ = 0.1 MeV). The fits for the D sJ (2460) + and D s1 (2536) + mesons are shown in Fig. 34 and correspond to mass resolutions of approximately 1.3 and 1.9 MeV/c 2 , respectively. The D sJ (2460) + and D s1 (2536) + masses are allowed to vary in the fit and the D * sJ (2317) + mass is fixed to the value determined in Section IX.
The result of the likelihood fit is shown in Fig. 33. A D sJ (2460) + mass of (2459.7 ± 0.2) MeV/c 2 and The uncertainties in the D sJ (2460) + and D s1 (2536) + masses are summarized in Table VII. The largest uncertainty (0.6 MeV/c 2 ) is in the assumed D + s mass. In comparison, uncertainties associated with the likelihood fit, such as the background shape or mass resolution, are relatively small. The remaining uncertainties are attributed to potential momentum biases in the tracking detectors.
Uncertainties in tracking have several sources, discussed in detail as part of a recent measurement of the Λ c mass from this collaboration [16]. A similar treatment of these uncertainties is reproduced for this analysis. The magnetic field of the BABAR detector is known to high precision, and variations of both the overall strength of the solenoidal field and of magnetization effects of PEP-II beam elements produce no significant variation in measured D + s π + π − mass. Residual φ-dependent momentum tracking biases are also insignificant. The largest uncertainty arises from material inside the tracking volume and the effect of this material on energy loss corrections. Studies of large samples of Λ and K S decays suggest that either the amount of material and/or atomic weight composition of the SVT is underestimated by approximately 20%. MC studies are used to estimate the bias introduced by this underestimation. A correction of 0.46 MeV/c 2 and 0.29 MeV/c 2 for the D sJ (2460) + and D s1 (2536) + masses is calculated based on these studies. The final results are: where the first error is statistical and the second is systematic.
The D sJ (2460) + and D s1 (2536) + signal PDF used in the likelihood fit of Fig. 19 was based on intrinsic widths of Γ = 0.1 MeV. After applying the same likelihoodintegration technique described in Section IX for the D + s π 0 final state, and after considering systematic uncertainties in the background shape and mass reconstruction resolution, the result is a 95% CL limit of Γ < 3.5 MeV and Γ < 2.5 MeV for the D sJ (2460) + and D s1 (2536) + mesons.

XVI. COMBINED MASS AND WIDTH RESULTS
The mass and width results for the D * sJ (2317) + , D sJ (2460) + , and D s1 (2536) + mesons are summarized in Table IX. The D * sJ (2317) + and D s1 (2536) + mesons are observed in only one decay mode covered by this analysis; those portions of this table are copied unchanged from the respective sections of this paper. The D sJ (2460) + mass is the average of that obtained from the D + s γ, D + s π 0 γ, and D + s π + π − final states, although the latter measurement dominates in the average due to superior systematic uncertainties. This average is calculated using a χ 2 method that properly accounts for correlations among the systematic uncertainties in the three measurements. The limit on the intrinsic D sJ (2460) + width Γ is taken as the best limit obtained from these three decay modes.

XVII. YIELDS AND BRANCHING RATIOS
The eighteen decay yieldsσ measured in this paper are collected in Table X. As described by Eq. 8, these numbers correspond to total yields for the decay of mesons having a center-of-mass momentum p * of at least 3.2 GeV/c to a final state that includes a D + s meson that decayed to φπ + .
A 95% CL upper limit is calculated for those yields which are not statistically significant. These limits are calculated using a frequentist approach by determining in each case the hypothetical value ofσ that is 1.96 standard deviations above the measured values. In order to calculate the systematic uncertainty associated with a hypothetical value ofσ, those uncertainties that are pro- [a] Denominator includes both D * s (2112) + π 0 and D * sJ (2317) + γ channels.
portional to the signal yield (such as uncertainties related to selection efficiency) are scaled as appropriate. The results are shown in Table X.
The yields listed in Table X are used to calculate  the branching ratios shown in Table XI. For the D * sJ (2317) + meson, only one decay mode has been observed; this is used as the denominator when calculating the D * sJ (2317) + branching ratios. For the D sJ (2460) + meson, the D + s π 0 γ decay mode (consisting of possible decay through either D * s (2112) + π 0 or D * sJ (2317) + γ) is chosen for this role. For completeness, the yield from hypothetical D * sJ (2317) ++ → D + s π + and D * sJ (2317) 0 → D + s π − decays is compared to D * sJ (2317) + → D + s π 0 decay to produce a quantity that is proportional to both the respective branching ratios and production rates.
In order to calculate the systematic uncertainty in the branching ratio results, systematic uncertainties common to the nominator and denominator are first discarded. Such common uncertainties include those associated with the D + s → φπ + branching ratio correction and, in some cases, those associated with π 0 and γ reconstruction efficiency.
A 95% CL upper limit is calculated for those branching ratios associated with decay modes that lack statistical signficance. The method used is similar to the frequentist recipe used for the upper limits on the cross section. The results are included in Table XI.

XVIII. DISCUSSION
The results in this paper confirm previous measurements, with generally higher precision. No new decay modes have been uncovered. Lacking any additional evidence to the contrary, the J P = 0 + and J P = 1 + assignments for the D * sJ (2317) + and D sJ (2460) + mesons remain a viable hypothesis. Except, perhaps, for the mass, there is currently no conflict with the interpretation of the D * sJ (2317) + and D sJ (2460) + mesons as the J P = 0 + and J P = 1 + p-wave cs states, as shown in Fig. 1.
The D sJ (2460) + meson mass measured in this analysis is an improvement over previous measurements [2][3][4]. High precision is obtained by using a decay mode that includes only charged particles coupled with a detailed understanding of the performance of the BABAR tracking detectors. The determination of the D * sJ (2317) + mass remains limited by uncertainties in EMC energy scale.
No intrinsic width of any statistical signficance is observed for either the D * sJ (2317) + or D sJ (2460) + meson. There is sufficient detector resolution in the D + s π + π − decay mode that better limits for the D sJ (2460) + meson should be attainable in the future with higher statistics. In contrast, limits on the D * sJ (2317) + width are currently limited by systematic uncertainties and will be difficult to improve.
The remaining branching ratio results shown in Table XI are generally competitive with the corresponding BELLE limits [3]. An exception is the limit on the decay D * sJ (2317) + → D + s γ. The BELLE publication, however, does not address the difficulties associated with modeling the D * sJ (2317) + → D + s π 0 reflection (see Fig. 14) which is the source of much of the systematic uncertainty quoted in this analysis.
The searches for the D + s π 0 π 0 and D + s γγ decay modes of the D * sJ (2317) + and D sJ (2460) + have been reported in this analysis for the first time. Both final states suffer from large backgrounds. More precise information on these modes will require more advanced analysis techniques, perhaps with the use of B meson decay or some other method of background suppression.
The D * sJ (2317) + and D sJ (2460) + masses are considerably lower than predictions from potential models developed before their discovery.
Since the D * sJ (2317) + /D sJ (2460) + mass splitting is so much larger than that for the D s1 (2536) + and D s2 (2573) + mesons, it would appear that this conflict is intrinsic. Studies from Cahn and Jackson [17] suggest, however, that it is possible with some adjustment to fit the D * sJ (2317) + and D sJ (2460) + into a perturbative model. In any case, given the approximate nature of these models, there is a danger that even large discrepancies can be overstated.
Perturbative calculations from Godfrey [19] predict branching ratios of 0.19 and 0.55 for the radiative transition decay of the D * sJ (2317) + and D sJ (2460) + to D * s (2112) + γ. These predictions are not consistent with the results reported here. This discrepancy could be considered evidence that both the D * sJ (2317) + and D sJ (2460) + are not cs mesons. These theoretical predictions, however, are plagued with large uncertainties in the partial widths of the isospin-violating D * sJ (2317) + → D + s π 0 and D sJ (2460) + → D * s (2112) + π 0 decays. In addition, leading-order corrections to the radiative transition partial rates might explain current experimental limits [20]. Nevertheless, the possibility that the D * sJ (2317) + and D sJ (2460) + are some type of new bound state, such as a four-quark DK molecule, should be taken seriously.
The possibility of four-quark states has long been proposed [21,22]. Perhaps the most unambiguous signature of the molecular interpretation of the D * sJ (2317) + meson would be the production of neutral and charged partners decaying to D + s π ± . For an isospin 1 DK molecule, production of these isospin partners is expected at the same rate as the D * sJ (2317) + [9]. This is clearly ruled out by our data. Other molecular interpretations, however, do not have this limitation [9,10], and it may be challenging to rule them all out with the data from the B-factories alone.
Finally, now that more precise data on the D * sJ (2317) + and D sJ (2460) + is available, measurements of the other cs states have become more important. The measurement of the D s1 (2536) + mass in this analysis will be helpful in further constraining models. With high statistics samples available from both B-factory experiments, improvements in the parameters of the D + s , D * s (2112) + , and D s2 (2573) + mesons should not be neglected and will hopefully follow soon.

XIX. CONCLUSION
An updated analysis of the D * sJ (2317) + and D sJ (2460) + mesons using 232 fb −1 of e + e − → cc data is performed. Established signals from the decay D * sJ (2317) + → D + s π 0 and D sJ (2460) + → D + s π 0 γ, D + s γ, and D + s π + π − are confirmed. A detailed analysis of invariant mass distributions of these final states including consideration of the background introduced by reflections of other cs decays produces the following mass values: where the first error is statistical and the second systematic. Upper 95% CL limits of Γ < 3.8 MeV and Γ < 3.5 MeV are calculated for the intrinsic D * sJ (2317) + and D sJ (2460) + widths. All results are consistent with previous measurements.
Since the results presented here are consistent with J P = 0 + and J P = 1 + spin-parity assignments for the D * sJ (2317) + and D sJ (2460) + mesons, these two states remain viable candidates for the lowest lying p-wave cs mesons. The lack of evidence for some radiative decays, in particular D * sJ (2317) + → D * s (2112) + γ and D sJ (2460) + → D * s (2112) + γ, are in contradiction with this hypothesis according to some calculations, but large theoretical uncertainties remain. No state near the D * sJ (2317) + mass is observed decaying to D + s π ± . If charged or neutral partners to the D * sJ (2317) + exist (as would be expected if the D * sJ (2317) + is a four-quark state), some mechanism is required to suppress their production in e + e − collisions.

Acknowledgments
We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also relies criti-cally on the expertise and dedication of the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them. This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), Institute of High Energy Physics (China), the Commissariatà l'Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), the Bundesministerium für Bildung und Forschung and Deutsche Forschungsgemeinschaft (Germany), the Istituto Nazionale di Fisica Nucleare (Italy), the Foundation for Fundamental Research on Matter (The Netherlands), the Research Council of Norway, the Ministry of Science and Technology of the Russian Federation, and the Particle Physics and Astronomy Research Council (United Kingdom). Individuals have received support from CONACyT (Mexico), the Marie-Curie Intra European Fellowship program (European Union), the A. P. Sloan Foundation, the Research Corporation, and the Alexander von Humboldt Foundation.