Angular analysis of the $e^+ e^- \to D^{(*) \pm} D^{* \mp}$ process near the open charm threshold using initial-state radiation

We report a new measurement of the exclusive $e^+ e^- \to D^{(*) \pm} D^{*\mp}$ cross sections as a function of the center-of-mass energy from the $D^{(*) \pm} D^{* \mp}$ threshold through $\sqrt{s}=6.0$ GeV with initial-state radiation. The analysis is based on a data sample collected with the Belle detector with an integrated luminosity of $951$ fb$^{-1}$. The accuracy of the cross section measurement is increased by a factor of two over the first Belle study. We perform the first angular analysis of the $e^+ e^- \to D^{* \pm} D^{* \mp}$ process and decompose this exclusive cross section into three components corresponding to the $D^*$ helicities.


INTRODUCTION
Parameters of the vector charmonia (ψ's) lying above the open-charm threshold were obtained from the total e + e − cross section to hadronic final states in measurements by MARK-I, DELCO, DASP, MARK-II, Crystal Ball, and BES [1]. In 2008, BES remodeled their data, including in particular the interference and the relative phases but using model predictions for the ψ decays to two-body open-charm final states [2]. As a result, the obtained parameters of the ψ states cannot be regarded as universal and model-indepenent. Notably, BES did not extract the parameters of other vector charmoniumlike resonances, such as Y(4008), Y(4260), Y(4360) and Y(4660), which were observed by BaBar [3][4][5] and Belle [6][7][8][9][10] in e + e − annihilation. Y(4008) and Y(4260) decay into the J/ψ π + π − final state while Y(4360) and Y(4660) decay into ψ(2S) π + π − . The Y states do not appear as peaks either in the total hadronic cross section nor in the exclusive e + e − cross sections to opencharm final states that were measured later. In addition, there exist predictions in the literature that some of the Y states could manifest themselves as coupled-channel effects [11].
A detailed study of the exclusive e + e − cross sections to open-charm final states could help establish parameters of the vector charmonia and charmoniumlike states in a model-independent way and, therefore, to shed light on the nature of the Y-family. Such exclusive e + e − cross sections to various open-charm final states were first measured at B-factories by Belle [12][13][14][15][16][17] and BaBar [18][19][20] using the initial-state radiation (ISR) processes, and by CLEO [21] using an energy scan in the range , Λ + c Λ − c ) and threebody DD ( * ) π cross sections almost saturates the total hadronic cross section after the subtraction of the u-, d-, and s-continuum in the region √ s ≤ 5 GeV [22]. The main contribution to the inclusive open-charm cross section comes from the DD, DD * , and D * D * final states.
The first attempt to extract the parameters of the ψ states (in particular, their couplings to the open-charm channels) from a combined coupled-channel fit for all exclusive open-charm cross sections was performed in Ref. [23]. At the time, although the suggested approach provides a good overall description of the line shapes, reliable conclusions could not be made because of the limited statistical accuracy of the data and because of the absense of experimental information on each of the three helicity amplitudes of the D * D * system. For further details, see Ref. [23] and references therein.
Here, we report a new measurement of the exclusive e + e − → D ( * )+ D * − cross sections as a function of the center-of-mass energy near the D ( * )+ D * − threshold in the initial-state radiation processes 1 . Compared to Ref. [13], the larger data set, the improved track reconstruction, and the additional modes used in the D and D * reconstruction allow one to obtain more precise determination of these cross sections. We also perform the first angular analysis of the e + e − → D * ± D * ∓ processes and explicitly decompose these exclusive cross sections into the three components corresponding to the D * 's helicities. Employing a coupled-channel technique -following the lines of Ref. [23], for example-a future study could use these results to extract the parameters of the vector charmonium states.

DATA SAMPLE AND BELLE DETECTOR
The analysis reported in this work is based on the data sample with an integrated luminosity of 951 fb −1 collected with the Belle detector [24] at the KEKB asymmetric-energy e + e − collider near the energies of the Υ(4S) and Υ(5S) resonances [25].
The Belle detector is a large-solid-angle magnetic spectrometer that consists of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), a barrellike arrangement of time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter (ECL) composed of CsI(Tl) crystals located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux-return located outside of the coil is instrumented to detect K 0 L mesons and to identify muons (KLM). A detailed description of the detector can be found, for example, in Ref. [24].

METHOD
To select the e + e − → D ( * )+ D * − γ ISR signal, we use the method described in Ref. [13]. We require full reconstruction of only one of the D * + (for e + e − → D * + D * − ) or D + mesons (for e + e − → D + D * − ), the γ ISR photon, and the slow pion from the other D * − . The partial D * − reconstruction without reconstruction of theD 0 daughter of the D * − increases the overall efficiency by a factor of ∼ 20 for e + e − → D * + D * − and ∼ 10 for e + e − → D + D * − , while suppressing the backgrounds enough to be able subtract them reliably using the data. Unlike the usual method for reconstruction of ISR processes, where the hadronic final state is fully reconstructed and γ ISR is inferred from the spectrum of masses recoiling against the hadronic system, we require here that the γ ISR to be reconstructed. This requirement does not significantly de-crease the efficiency as the slow pion from D * ± decay has a low reconstruction efficiency when γ ISR is outside the detector acceptance (because of the very low transverse momentum of γ ISR in this case).
For the signal candidates, the spectrum of the mass recoiling against the D ( * +) γ ISR system, 2 peaks at the D * − mass. Here, E cm is the e + e − centerof-mass energy, and E D ( * )+ γISR and p D ( * )+ γISR are the center-of-mass energy, and momentum, respectively, of the D ( * )+ γ ISR system. According to a Monte Carlo (MC) study, this peak is wide (σ ∼ 150 MeV/c 2 ) and asymmetric due to the asymmetric photon energy resolution function and higher-order ISR corrections (i.e., more than one γ ISR in the event). Because of poor M rec (D ( * )+ γ ISR ) resolution, the signals from DD, DD * , and D * D * overlap strongly overlap; hence, one cannot distinguish these processes using this measure alone. For illustration, Fig. 1 shows the MC M rec (D * + γ ISR ) spectra for the e + e − → D * + D * − γ ISR and e + e − → D + D * − γ ISR processes. To resolve this, we use the information provided by the slow pion from the unreconstructed D * − meson. The distribution of the difference between the masses recoiling against the D ( * )+ γ ISR and D ( * )+ π − slow γ ISR , termed the recoil-mass difference, has a narrow peak for the signal process around the m D * + − m D 0 mass difference (Fig. 2, histogram). The resolution of this peak is below 2 MeV/c 2 as the uncertainty of the γ ISR momentum is mostly canceled out for 2 Hereinafter, the speed of light is set to unity for convenience. this variable. Thus, the existence of a partially reconstructed D * − in the event is identified by the presence of this peak. The method does not exclude contributions in the ∆M rec signal window from processes with extra neutrals in the final state (e.g., e + e − → D ( * )+ D * − π 0 ). However, this background is suppressed and its residual contribution can be reliably determined using the data, as discussed in the data analysis section.
To measure the exclusive cross sections as a function of √ s, one needs to obtain the D ( * )+ D * − mass spectrum despite one of D * mesons being unreconstructed. In the absence of higher-order QED processes, the D ( * )+ D * − mass corresponds to the mass recoiling against the single ISR photon: M (D ( * )+ D * − ) ≡ M rec (γ ISR ). However, the poor M rec (γ ISR ) resolution (σ ∼ 120 MeV/c 2 ) precludes the study of relatively narrow charmonium states in the D ( * )+ D * − mass spectra. To improve the M rec (γ ISR ) resolution, we refit the recoil mass against the D ( * )+ γ ISR system, constrained to the D * − mass. This procedure utilizes the well-measured momentum of the reconstructed D ( * )+ meson and the signal kinematics to better determine the momentum of the ISR photon. It works well even in case of a second ISR photon, as checked with the MC. As a result, the M rec (γ ISR ) resolution is drastically improved: near the threshold, the resolution is better than 3 MeV/c 2 , and smoothly increases to 15 MeV/c 2 at √ s ∼ 6 GeV (Fig. 3). The resolution of the recoil-mass difference after refit, ∆M fit rec , improves by a factor of ∼ 2 (Fig. 2, points with error bars); this is exploited for more effective suppression of the combinatorial background.

MONTE CARLO STUDY AND CALIBRATION OF γISR ENERGY
The simulation of the signal and background processes up to the second-order ISR corrections and vacuumpolarization corrections is performed using the PHOKHARA MC generator [26]. In the signal-MC samples, the D ( * )+ D * − characteristics (mass spectrum and angular distributions) are tuned to those measured in the data. As the measured distributions are extracted from the data using the MC simulation, this tuning is repeated until the difference between successive iterations is negligibly small.
To improve the accuracy of the MC simulation, we calibrate the photon energy. For soft and medium energy photons, this is done using π 0 → γγ [27]. For energetic (∼ 4 GeV) photons, where there are few energetic π 0 in the data, we select a clean sample of fully reconstructed e + e − → ψ(2S)γ ISR → J/ψ π + π − γ ISR events and study the spectrum of the mass recoiling against γ ISR ; a broad peak around the ψ(2S) mass is expected. A small shift between MC and data is observed, from which a correction factor for the photon energy is determined to be 0.9980 ± 0.0004. We apply this correction throughout the analysis.

DATA SAMPLE AND EVENT SELECTION
In the first Belle analysis [13], the strategy was to select a clean sample of the studied process with minimal background contribution to provide the most reliable result. It was demonstrated there that all backgrounds were well under control and were subtracted reliably using the data. The aim of the present analysis is to improve the accuracy of the cross section measurement. Therefore, we reoptimize the selection criteria and add more D ( * ) decay modes.
All charged tracks are required to be consistent with origination from the interaction point (IP): we require dr < 2 cm and |dz| < 4 cm, where dr and |dz| are the impact parameters perpendicular to and along the beam direction, respectively, with respect to the IP. Information from the TOF, the number of the photoelectrons from the ACC and the dE/dx measurement in the CDC are combined to form a likelihood L for hadron identifi-cation. Charged kaon candidates are required to have a kaon/pion likelihood ratio P K/π = L K /(L K + L π ) > 0.6. The identification efficiencies typically exceed 90%, while misidentification probabilities are less than 10% [28]. No identification requirements are applied for pion candidates, as the pion multiplicity is much higher than that of other hadrons. K 0 S candidates are reconstructed by combining π + π − pairs with an invariant mass within ±15 MeV/c 2 (≈ ± 5σ) of the K 0 S mass. The distance between the two pion tracks at the K 0 S vertex must be less than 1 cm. The transverse flight distance from the interaction point is required to be greater than 0.1 cm, and the angle between the K 0 S momentum direction and the line joining the decay vertex and the IP in the transverse plane should be smaller than 0.1 rad. Photons are reconstructed in the electromagnetic calorimeter as showers with energy greater than 50 MeV that are not associated with charged tracks. ISR photon candidates are required to have a center-of-mass energy above 3.0 GeV. Combinations of two photons are treated as a π 0 candidates if their invariant mass lies within 15 MeV/c 2 (≈ ± 3σ) of the nominal π 0 mass. Such combinations are then refitted with a π 0 mass constraint to improve the π 0 momentum resolution.
D 0 candidates are reconstructed in nine decay modes: The mass window for all modes without (with) π 0 is ±15 MeV/c 2 (± 20 MeV/c 2 ) of the nominal D 0 mass [1], corresponding to ≈ ± 3σ in each case. A mass-constrained vertex fit is applied to the D 0 and D + candidates to improve their momentum resolution.
D * + candidates are selected in the D 0 π + decay mode; the D * 0 candidates are selected in the D 0 π 0 mode and are used for background studies. In both cases, the signal window for the D * mass is chosen to be ±3 MeV/c 2 (≈ 3σ) of the nominal D * + mass [1].
To increase the efficiency, the M rec (D ( * )+ γ ISR ) signal region is defined to be ±300 MeV/c 2 of the D * − meson mass (see Fig. 1).
The process e + e − → D ( * )+ D * − π 0 miss γ ISR has the M rec (D ( * )+ γ ISR ) spectrum shifted to higher values, but still overlaps this signal window. We eschewed a further extension of the signal window to the higher masses since this would lead to an to increased contribution from the background process with an extra π 0 . The ∆M fit rec value is required to be within ±3 MeV/c 2 of the m D * + − m D 0 value (see Fig. 2).
The fraction of events with more than one candidate passing all selection criteria is small (3%). For these, we select the one with the smallest value of where χ 2 i is defined as a squared ratio of the difference between the measured and expected observable i to the corresponding resolution. For the background studies, we use M (D ( * )+ ) and ∆M fit rec (D ( * )+ γ ISR ) sidebands. In the sideband regions, a similar selection of a single candidate per event is performed based on a χ 2 calculated with respect to the center of the corresponding sideband interval.

DATA ANALYSIS
As the process e + e − → D * + D * − γ ISR contributes to the e + e − → D + D * − γ ISR when one of the D * ± mesons decays into D + π 0 or D + γ, we study the process e + e − → D * + D * − first. Then we use the obtained result to correct the e + e − → D * + D * − γ ISR Monte Carlo and use the latter to extract the information that is needed to subtract the D * + D * − contribution from the e + e − → D + D * − γ ISR process.
The D * + D * − mass spectrum, after applying all the selection criteria, is shown in Fig. 4. Taking advantage of the fine M (D * + D * − ) resolution at the threshold and higher statistics in comparison with the first Belle paper [13], we study the mass structure near the threshold with finer bins. The D * + D * − mass spectrum is presented in the inset of Fig. 4. We consider the following background contributions: 1. combinatorial background under the reconstructed D * + peak, when the slow pion is truly a D * + daughter; 2. real D * + mesons combined with a combinatorial (random) slow pion; 3. both the D * + meson and slow pion are combinatorial; 4. reflections from the processes e + e − → D * + D * − π 0 γ ISR where the π 0 is lost; fast where the hard π 0 fast is misidentified as γ ISR . This can happen if the two photons from the π 0 fast decay merge into one ECL cluster or, in the case of an asymmetric π 0 fast decay, if one of the photons carries a large fraction of the π 0 fast energy. The contribution from the combinatorial backgrounds, i.e., the first three background sources, is extracted using the M (D * + ) and ∆M fit rec sidebands. The M (D * + ) vs ∆M fit rec scatter plot is shown in Fig. 5 for the data. The signal region (indicated as box A) is defined by the requirement that M (D * + ) and ∆M fit rec lie within ±3 MeV/c 2 of the corresponding nominal values. The contributions of combinatorial D * + candidates or random slow pions -backgrounds (1) and (2) -are estimated from the regions shown by boxes B and C, respectively. The sidebands of double the width are shifted from the signal region to avoid signal over-subtraction. Background (3) is present in both M (D * + ) and ∆M fit rec sidebands and therefore is subtracted twice. To correct for this over-subtraction, the double sideband region D is used.  (4) where M B, C, D are the D * + D * − mass spectra from the B, C, and D sidebands, respectively, and the scaling factors are calculated to provide normalization of the corresponding background contributions within the signal window. To obtain these scaling factors, the distributions of M (D * + ) and ∆M fit rec are fit using signal shapes fixed to those from the MC simulation and the following background parameterization: where x is the M (D * + ) or ∆M fit rec , x thr is a corresponding threshold values, and a, b, and c are free parameters. The fit results are shown in Fig. 6 To estimate the contribution from background (4), the process e + e − → D * + D * − π 0 miss γ ISR , we study the isospin-conjugated process e + e − → D * 0 D * − π + miss γ ISR using a similar partial reconstruction method. Because of charge imbalance of the reconstructed combination (D * 0 π − slow γ ISR ), only events with a missing π + miss can contribute to the recoil mass difference peak (indicating the presence of a D * − meson). The D * 0 D * − mass spectrum is obtained by applying all requirements and subtracting combinatorial backgrounds. To recalculate the background (4) contribution from the isospin-conjugated M (D * 0 D * − ) spectrum, one needs to correct the latter for the ratio of the D * + and D * 0 reconstruction efficiencies (equal to 1.77, according to the MC) and the isospin factor 1/2. As in the previous Belle analysis [13], this contribution is found to be consistent with zero (Fig. 7). However, to avoid additional systematic errors, the spectrum of background (4) obtained from the data is subtracted bin-by-bin, and the statistical errors of the subtraction of this contribution are included in the final result.
The D * + D * − mass spectrum after bin-by-bin combinatorial background subtraction is shown in Fig. 8. Background (5) is also estimated from the data us- ing the same method of partial reconstruction, substituting γ ISR for a fully reconstructed energetic π 0 . From the fit to the π 0 mass distribution, we find 56 ± 12 events corresponding to the process e + e − → D * + D * − π 0 with M (D * + D * − ) ≤ 6 GeV/c 2 . From the MC simulation, we estimate the ratio of efficiencies for background process (5) to be reconstructed as D * + D * − γ ISR and D * + D * − π 0 to be equal to 0.39. Thus, we conclude that only 22 ± 5 events contribute to our reconstructed ISR spectrum in the region M (D * + D * − ) ≤ 6 GeV/c 2 , i.e., under 0.3% of the signal. As the M (D * + D * − ) spectrum from e + e − → D * + D * − π 0 is distributed uniformly, far less than one event per bin from this background is expected. Therefore, we incorporate this background contribution into the systematic uncertainty.
Measurement of e + e − → D + D * − γISR The D + D * − mass spectrum, after applying all the selection criteria, is shown in Fig. 9.
The contribution from the combinatorial backgrounds (1)-(3) is estimated using two-dimensional sideband regions of the D + candidate mass versus the recoil mass difference (Fig. 10).  Fig. 11. In the case of the M (D + ) sidebands, the scaling factor turns out to be exactly related to the ratio of widths for the sideband and signal windows, as the background is well described by a linear function. We note that background (4) contributes to the ∆M fit rec spectrum peak. However, this peak is wider than the signal's as the refitting procedure of the recoil mass against D ( * )+ γ ISR into the D * − mass for this background process works improperly and does not improve the ∆M fit rec resolution. We fix this contribution (shown in Fig. 11  To estimate the contribution of the background process (4a), the isospin-conjugated process e + e − → D 0 D * − π + γ ISR is studied using the same method of partial reconstruction by replacing D + with D 0 . Using the MC simulation, we verify that this method gives an accurate subtraction of the background without bias. The measured M (D 0 D * − ) spectrum is shown in Fig. 12a as the points with error bars. We repeat the procedure of the subtraction of combinatorial background (shown by the hatched histogram in Fig. 12a similarly to the studied process and obtain the net M (D 0 D * − ) spectrum as shown in Fig. 12b. The background process e + e − → D + D * − π 0 γ ISR includes contributions from e + e − → D * + D * − γ ISR with D * + → D + π 0 as well as higher resonance D * * + → D + π 0 miss or non-resonant D + π 0 miss contributions. For the two latter sources, the ratio D + π 0 to D 0 π + is related by the isospin factor 1/2. A similar factor is valid for the ratio of B(D * + → D + π 0 )/B(D * + → D 0 π + ) = 0.453 ± 0.011 [1]. We do not study all of these contributions separately, but rather scale the measured D 0 D * − spectrum for the ratio of efficiencies for D + and D 0 (0.28 according to MC) and the factor 0.453, assuming the main contribution has π 0 miss coming from D * + . The difference between B(D * + → D + π 0 )/B(D * + → D 0 π + ) and the isospin factor is taken into account in the systematic error of the result. The measured contribution from the process e + e − → D + D * − π 0 γ ISR to the D + D * − mass spectrum is shown in Fig. 13 as the open circles.
The small contribution of background (4b), from e + e − → D * + D * − γ ISR with D * + → D + γ, is estimated using the MC simulation. In the e + e − → D * + D * − γ ISR MC sample, the D * + D * − spectrum is set according to our measurement. Background (4b) is tiny in comparison with (4a) because of the small B(D * + → D + γ) branching fraction and a greatly smeared ∆M fit rec distribution in the case of a lost γ.
Background (5) is estimated similarly to the study of e + e − → D * + D * − and found to be negligibly small. Its contribution is incorporated into the systematic uncertainty.

CROSS-CHECKS
We performed several checks to determine that the background subtraction procedure does not bias the measured spectra.
To verify that the combinatorial background subtraction is performed correctly, each of B and C sideband regions is divided into two equal intervals: B 1 -B 2 and C 1 -C 2 , respectively, as shown in Fig. 10. Then, the consistency of both the M (D * + D * − ) shapes and normalizations in the pairs of subintervals is checked. The differences between these two pairs of control spectra are consistent with zero. The same procedure, repeated for the processes e + e − → D + D * − , also demonstrates good agreement between the M (D + D * − ) spectra from differ- ent sidebands regions. We thus conclude that the combinatorial background shapes and normalizations are well understood.
We check that, after the energy correction of fast photons, the spectra of masses recoiling against D ( * )+ γ ISR combinations in the data are consistent with the MC expectations. Figures 14a and 14b show the M rec (D ( * )+ γ ISR ) spectra in the data (points with errors) and signal MC events (hatched histogram) for the processes e + e − → D + D * − γ ISR and e + e − → D * + D * − γ ISR , respectively. The good agreement between MC and data spectra demonstrates not only the correctness of the photon energy calibration but also the proper background subtraction.
where M is the D ( * )+ D * − mass, equivalent to √ s, dN/dM is the measured mass spectrum, η tot is a M -dependent total efficiency, and dL/dM is the differential luminosity.
The dependence of the efficiency on M is calculated using the MC simulation and is defined as the ratio of reconstructed mass spectrum for true MC candidates after applying all requirements to the generated spectrum (Fig. 15). Alternatively, we calculate the efficiency without the requirement for the selected MC candidates to be the true combination, but repeating in this case the procedure of the background subtraction. The latter method checks for possible oversubtraction of the signal by the applied procedure of the combinatorial background subtraction due to tails in the M (D * + ) and recoil mass difference resolution functions. It is found that the oversubtraction is negligibly small, and both methods result in the same η tot (M ) dependence. The differential luminosity is calculated as a sum over all energy points -Υ(4S), Υ(5S), and continuumusing the known luminosities of each data sub-sample. We use the dL/dM formula that includes the secondorder QED corrections [29]. The latter varies from 2.5 to 3.5% of the leading contribution in the studied √ s interval. In the previous Belle paper, this was treated as a systematic uncertainty.

Efficiency
Finally, the obtained exclusive e + e − → D + D * − and e + e − → D * + D * − cross sections are shown in Fig. 16.

ANGULAR ANALYSIS
We study the D * ± helicities from both processes. The D * ± helicity angle, θ, is defined as the angle between the π ± slow from D * ± decay and the D ( * )+ D * − system, seen in the D * ± rest frame. The angular distributions of the D * − decays can be studied even in the case of partial D * − reconstruction. The refit procedure of M rec (D * + γ ISR ) to m D * − provides sufficient accuracy in the unreconstructed D * − momentum determination. The helicity angle resolution of the partially reconstructed D * ± mesons (σ θ ∼ 0.06 rad) is slightly worse than for fully reconstructed D * ± mesons (σ θ ∼ 0.05 rad).
For the e + e − → D + D * − process, the helicity of the D * − meson is uniquely defined by the angular momentum and parity conservation: the D * − meson polarization should be transverse. Thus, we perform the D * − angular analysis for this process to verify the method only.
The D * − helicity angle distribution is analyzed in each bin of M (D + D * − ). Figure 17  combinatorial background contribution, taken from the M (D + ) and ∆M fit rec sidebands, is shown as the hatched histogram, and the sum of the combinatorial background and D * + D * − feed-down is shown as the open histogram; the latter contribution is obtained from the D 0 D − data, corrected for the D 0 and D − efficiency ratio and the ratio of B(D * + → D + π 0 )/B(D * + → D 0 π + ). In each bin of M (D + D * − ), the cos θ distribution after subtraction of all background contributions is fitted with the function F (cos θ) = η(cos θ) · dM/dL · (f L + f T ), (8) where η(cos θ) is the efficiency depending on the helicity of the D * − , dL/dM is calculated in each M bin, and f L = σ L · cos 2 θ and f T = σ T · (1 − cos 2 θ) are the contributions from distinct D * − helicity states. The subscripts L and T refer to longitudinally and transversely polarized D * ± mesons, correspondingly. The fitting procedure scans over M bins and the fits return the cross sections σ T,L ( √ s) for T and L components. Figure 18    We analyze the two-dimensional distribution of c 1 ≡ cos θ f vs. c 2 ≡ cos θ p , where the first helicity angle, θ f , corresponds to the fully reconstructed D * + , and the second, θ p , is calculated for the partially reconstructed D * − . As an illustration, the two-dimensional distributions of c 1 versus c 2 for four mass ranges are shown in Fig. 19.
We perform binned maximum likelihood fits to these distributions in bins of M (D * + D * − ). The fitting function is an incoherent sum of the LL and T L and T T contributions and the background component: Here, η(c 1 , c 2 ) is the efficiency, depending on the two D * helicity angles, and f LL , f T L , f T T , and f bg are the contributions from the three mutually orthogonal signal components and from the background, respectively. In this study, we take into account the combinatorial background only, and ignore the contributions of backgrounds (4) and (5) because they are very small. We define: where σ LL , σ T L , and σ T T are free parameters in the fit that represent the cross sections of each component in M bins. The f B, C, D (c 1 , c 2 ) are elements of the matrices that are constructed from the binned c 1 vs. c 2 histograms for the corresponding sideband regions; the background normalization, N bg , is fixed in each M (D * + D * − ) bin within the error.
The efficiency map is shown in Fig. 20 a. The efficiency is almost symmetric with respect to c 1 vs. c 2 , as it depends mainly on the momenta of the two slow pions. The combinatorial background f bg is presented in Fig. 20 b.
The angular-fit results are plotted in Fig. 21. We observe that the cross sections corresponding to the different D * + D * − helicities have distinct dependencies on √ s. Near the threshold, the T T and T L components have a similar sharp rise, while the LL component rises slowly. This can be explained by the high centrifugal barrier for the LL component, which originates from the F wave (chapter 48 of Ref. [1]). All three components reach the same value of ∼ 1 nb at √ s ∼ 4.15 GeV and fall into a common dip at √ s ∼ 4.25 GeV. The LL and T T components are attenuated and only the T L component survives in the region of high √ s 4.5 GeV, which is in good agreement with theoretical predictions [31].

SYSTEMATIC ERRORS
The systematic errors in the cross section calculation for the studied processes are summarized in Table I.
The systematic errors arising from the background subtractions include the uncertainty in the calculation of the scaling factors for the sideband distributions in Eqs. 4 and 6, systematic errors in the determination of background (4) with missing neutrals and unsubtracted background (5). To estimate the uncertainty in the calculation of the scaling factors, we perform fits to Figs. 6 and 11 with different parameterizations and in different M (D * + D * − ) intervals. As a result, the scaling factor extracted from the integral under the signal and sideband regions varies within ±15%. In spite of large uncertainties in this scaling factor, the final systematic error due to the subtraction of combinatorial backgrounds (1)-(3) is estimated to be only 2% as these backgrounds are small (only 15% of the signal). The systematic error associated with subtraction of background (4) includes the uncertainty of the isospin-conjugated spectra and efficiency ratio for the isospin-conjugated final states. The upper limit on background (5) is considered as an extra contribution to the systematic error. The main contribution to the total systematic error comes from the uncertainties in track and photon reconstruction, estimated as 0.35% per track and 1.5% per photon. An extra uncertainty of 1% is ascribed to slow pion(s), and 2% to K 0 S reconstruction. Uncertainty of the kaon identification efficiency is 1%, evaluated by the study of inclusive D * mesons.
The systematic error of the selection efficiency comes from the possible difference of resolution and calibration in the MC and data. Particularly, the calibration of the fast photon can influence the efficiency of the requirement |M rec (D ( * )+ γ ISR ) − m D * − | < 300 MeV/c 2 . From the error in the photon energy correction factor obtained in our study, we estimate this uncertainty as ±1%. Other contributions are estimated by comparison of M (D ( * ) ) and ∆M fit rec distributions in the MC and data.
The uncertainty due to the D * meson helicity distributions are reduced in comparison to the first Belle analysis [13], because the angular distributions are analyzed. We set the measured D * helicity in the MC simulation to reduce the uncertainty of the slow pion reconstruction due to angular distribution in D * decays.
The systematic error ascribed to the cross section calculation is estimated from a study of the cos θ 0 dependence (θ 0 defines the polar angle range for γ ISR in the e + e − c.m. frame: θ 0 < γ ISR < 180 • − θ 0 ) of the final result and includes a 1.4% error on the total luminosity. We add uncertainties of the same origin linearly, while independent uncertainties for different types of particles are summed quadratically. Different D modes have slightly different reconstruction uncertainties; the average error is calculated according to the weight of each mode in the final sample.
Other contributions come from the uncertainty in absolute the D ( * ) branching fractions [1] and MC statistics.
We divide the total systematic errors in two parts: correlated and uncorrelated.
The √ s-independent correlated errors come from track reconstruction, selection efficiency, cross section calculation and uncertanty from the D ( * ) branching fractions. These errors influence the normalization of the measured cross section as a whole. Therefore, to take into account the correlated errors, the measured cross section should be multiplied by a factor of 1± 0.05 for e + e − → D + D * − and 1 ± 0.06 for e + e − → D * + D * − .
On the contrary, the uncorrelated errors vary in each √ s bin. These errors come from the background subtraction (for e + e − → D + D * − ) and from calculating efficiency dependence on the helicity of the D * + D * − final state (for e + e − → D * + D * − ). Uncorrelated systematic errors are given in Tables II and III.

SUMMARY
In summary, we report measurements of the exclusive e + e − → D + D * − and e + e − → D * + D * − cross sections at √ s near D + D * − and D * + D * − thresholds. These results supersede those of Ref. [13]. Due to the increased size of the data sample, an improved track reconstruction efficiency, and additional modes for the charmed meson reconstruction, the accuracy of the cross section measurements is increased by a factor of two compared to Ref. [13]. The systematic uncertainties are also significantly improved. We also extend the studied region up to √ s = 6 GeV and, taking advantage of the improved resolution and high statistics, we halve the size of the √ s steps close to threshold.
The complex shape of the e + e − → D * + D * − cross sections can be explained by the fact that its components can interfere constructively or destructively. The fit of this cross section is not trivial, because it must take into account the threshold and coupled-channels effects [23].
The first angular analysis of the e + e − → D * + D * − process allows us to decompose the corresponding exclusive cross section into three possible components for the longitudinally, and transversely-polarized D * ± mesons. The obtained components have distinct behavior near the D * + D * − threshold. The only non-vanishing component at higher energy is the T L helicity of the D * + D * − final state. The measured decomposition allows the future measurement of the couplings of vector charmonium states into different helicity components, useful in identifying their nature and in testing the heavy-quark symmetry [32].
We thank the KEKB group for the excellent operation of the accelerator; the KEK cryogenics group for the efficient operation of the solenoid; and the KEK computer group, the National Institute of Informatics, and the PNNL/EMSL computing group for valuable computing and SINET5 network support.