Measurement of e+e- -->D Dbar Cross Sections at the psi(3770) Resonance

We report new measurements of the cross sections for the production of D Dbar final states at the psi(3770) resonance. Our data sample consists of an integrated luminosity of 2.93/fb of e+e- annihilation data produced by the BEPCII collider and collected and analyzed with the BESIII detector. We exclusively reconstruct three D0 and six D+ hadronic decay modes and use the ratio of the yield of fully reconstructed D Dbar events ("double tags") to the yield of all reconstructed D or Dbar mesons ("single tags") to determine the number of D0 D0bar and D+D- events, benefiting from the cancellation of many systematic uncertainties. Combining these yields with an independent determination of the integrated luminosity of the data sample, we find the cross sections to be \sigma(e+e- -->D0 D0bar)=(3.615 +- 0.010 +- 0.038) nb and \sigma(e+e- -->D+D-)=(2.830 +- 0.011 +- 0.026) nb, where the uncertainties are statistical and systematic, respectively.


Introduction
The ψ(3770) resonance is the lowest-energy charmonium state above the threshold for decay to charmed meson pairs. The expectation that the ψ(3770) should decay predominantly to D 0D0 and D + D − has been validated by experiment [1], although inconsistent results for the branching fraction of ψ(3770) to non-DD final states have been reported [2,3]. The cross sections σ(e + e − →D 0D0 ) and σ(e + e − →D + D − ) at center-of-mass energy E cm =3.773 GeV, the peak of the ψ(3770) resonance, can be measured precisely and are necessary input for normalizing some measurements of charmed meson properties in ψ(3770) decays. The most precise determinations to date are from the CLEO-c Collaboration [4] using 818 pb −1 of e + e − annihilation data at E cm = 3774±1 MeV, σ(e + e − →D 0D0 )=(3.607±0.017±0.056) nb and σ(e + e − → D + D − ) = (2.882±0.018±0.042) nb. In this paper we report measurements of the DD cross sections using fully reconstructed D 0 and D + mesons in a ψ(3770) data sample that is approximately 3.6 times larger than CLEO-c's. Here and throughout this paper, charge-conjugate modes are implied unless explicitly stated.
Our procedure is an application of the D-tagging technique developed by the MARK III Collaboration [5], exploiting the kinematics of DD production just above threshold at the ψ(3770) resonance. We use ratios of fully reconstructed D mesons ("single tags") and DD events ("double tags") to determine the total numbers of DD pairs. This procedure benefits from the cancellation of systematic uncertainties associated with efficiencies and input branching fractions, giving better overall precision than measurements based on single tags. The production of D 0D0 pairs in a pure C =−1 state complicates the interpretation of measurements at ψ(3770) by introducing correlations between the D 0 andD 0 decays. We apply corrections derived by Asner and Sun [6] to remove the bias introduced by these correlations.

BESIII detector
Our measurement has been made with the BESIII detector at the BEPCII collider of the Institute for High Energy Physics in Beijing. Data were collected at the ψ(3770) peak, with E cm = 3.773 GeV. The integrated luminosity of this sample has previously been determined with large-angle Bhabha scattering events to be 2.93 fb −1 [7,8], with an uncertainty of 0.5% dominated by systematic effects. An additional data sample of 44.9 pb −1 at E cm = 3.650 GeV has been used to assess potential background from continuum production under the ψ(3770).
BESIII is a general-purpose magnetic spectrometer with a geometrical acceptance of 93% of 4π. Charged particles are reconstructed in a 43-layer helium-gas-based drift chamber (MDC), which has an average single-wire resolution of 135 µm. A uniform axial magnetic field of 1 T is provided by a superconducting solenoid, allowing the precise measurement of charged particle trajectories. The resolution varies as a function of momentum, and is 0.5% at 1.0 GeV/c. The MDC is also instrumented to measure the specific ionization (dE/dx) of charged particles for particle identification. Additional particle identification is provided by a time-of-flight system (TOF) constructed as a cylindrical ("barrel") structure with two 5-cm-thick plastic-scintillator layers and two "end caps" with one 5-cm layer. The time resolu-tion in the barrel is approximately 80 ps, and in the end caps it is 110 ps. Just beyond the TOF is an electromagnetic calorimeter (EMC) consisting of 6240 CsI(Tl) crystals, also configured as a barrel and two end caps. For 1.0-GeV photons, the energy resolution is 2.5% in the barrel and it is 5% in the end caps. This entire inner detector resides in the solenoidal magnet, which is supported by an octagonal flux-return yoke instrumented with resistive-plate counters interleaved with steel for muon identification (MUC). More detailed information on the design and performance of the BESIII detector can be found in Ref. [9].

Technique
To select a DD event, we fully reconstruct a D using tag modes that have sizable branching fractions and can be reconstructed with good efficiency and reasonable background. We use three D 0 and six D + tag modes: S π + π + π − , and D + →K − K + π + . When both the D andD in an event decay to tag modes we can fully reconstruct the entire event. These double-tag events are selected when the event has two single tags and satisfies the additional requirements that the reconstructed single tags have opposite net charge, opposite-charm D parents and no shared tracks. The yield X i for single-tag mode i is given by Eq. (1): where N DD is the total number of DD events, B(D → i) is the branching fraction for decay mode i, and i is the reconstruction efficiency for the mode, determined with Monte Carlo (MC) simulation. Extending this reasoning, the yields forD decaying to mode j and for ij doubletag events, in which the D decays to mode i and theD decays to mode j, are given as follows: and In these equations, Z ij is the yield for the double-tag mode ij, and ij is the efficiency for reconstructing both tags in the same event. Combining Eqs. (1), (2) and (3), N DD can be expressed as The cancellation of systematic uncertainties occurs through the ratio of efficiencies ij /( i · j ). The measured N DD from each combinations of i and j are then averaged, weighted by their statistical uncertainties. Finally, to determine cross sections we divide N DD by the inte-
Data and MC samples are treated identically for the selection of D tags. All particles used to reconstruct a candidate must pass requirements specific to the particle type. Charged particles are required to be within the fiducial region for reliable tracking (|cosθ| < 0.93, where θ is the polar angle relative to the beam direction) and to pass within 1 cm (10 cm) of the interaction point in the plane transverse to the beam direction (along the beam direction). Particle identification is based on TOF and dE/dx measurements, with the identity as a pion or kaon assigned based on which hypothesis has the higher probability. To be selected as a photon, an EMC shower must not be associated with any charged track [16], must have an EMC hit time between 0 and 700 ns to suppress activity that is not consistent with originating from the collision event, must have an energy of at least 25 MeV if it is in the barrel region of the detector (|cosθ|<0.8), and 50 MeV if it is in the end cap region (0.84<|cosθ|<0.92) to suppress noise in the EMC as a potential background to real photons. Showers in the transition region between the barrel and end cap are excluded.
K 0 S mesons are reconstructed from the decay into π + π − . Because of the cleanliness of the selection and the possibility of a measurably displaced decay vertex, the pions are not required to pass the usual particle identification or interaction-point requirements. A fit is performed with the pions constrained to a common vertex and the K 0 S candidate is accepted if the fit satisfies χ 2 < 100 and the candidate mass is within ∼ 3σ of the nominal K 0 S mass (487−511 MeV/c 2 ). The momentum of the K 0 S that is obtained from the constrained-vertex fit is used for the subsequent reconstruction of D-tag candidates. π 0 mesons are reconstructed through the decay into two photons. Both photons for a π 0 candidate must pass the above selection criteria, and at least one of them must be in the barrel region of the detector. To be accepted a π 0 candidate must have an invariant mass between 115 MeV/c 2 and 150 MeV/c 2 . The photons are then refitted with a π 0 mass constraint and the resulting π 0 momentum is used for the reconstruction of D-tag candidates.

Event selection
In addition to the requirements on the final-state particles, the reconstructed D-tag candidates must pass several additional requirements that ensure the measured candidate energy and momentum are close to the expected values for production via ψ(3770) → DD. The first of these requirements is ∆E =E D −E beam 0, where E D is the energy of the reconstructed D candidate and E beam is the beam energy. In calculating ∆E we use the beam energy calibrated with D 0 and D + decays, combining groups of nearby runs to obtain sufficient statistics. Selection requirements on ∆E are determined separately for each tag mode for data and MC to account for differing resolutions. As shown in Table 1, for modes decaying into all charged tracks, the requirements are set to ±3σ about the mean, while for modes with a π 0 , the requirements are asymmetric about the mean, extending on the low side to −4σ to accommodate the tail from the photon energy resolution. Figure 1 shows the data and MC overlays of the ∆E distributions by mode. Table 1. The selected range on ∆E is ±3σ about the mean, except that for modes with a π 0 an extended lower bound of −4σ is used. The resolutions and means are extracted by fitting with a double Gaussian, weighted by the two Gaussian yields, and determined separately for data and MC.
The second variable used in selecting D tags is the beam-constrained mass M BC c 2 = E 2 beam −|p tag c| 2 , where p tag is the 3-momentum of the candidate D. We use M BC rather than the invariant mass because of the excellent precision with which the beam energy is known. The requirement that M BC be close to the known D mass ensures that the D tag has the expected momentum. After application of the ∆E requirement to single-tag candidates of a given mode, we construct an M BC dis-083001-5 These plots overlay the 3.773 GeV data (blue dashed histograms) and the corresponding narrower-width MC (red solid histograms). Only requirements on the constituent particles and a very loose MBC requirement (1.83 GeV/c 2 MBC 1.89 GeV/c 2 ) have been applied. tribution in the region of the known masses of charmed mesons (1.83−1.89 GeV/c 2 ). For the MC a small upward shift of just under 1 MeV/c is applied to the measured D momentum for the calculation of M BC to compensate for input parameters that do not precisely match data. Initial inspection of the distribution in data for the twobody mode D 0 →K − π + exhibited peaking near the high end of the M BC range not seen in MC. We demonstrated this to be background from cosmic ray and QED events. To eliminate it from the distribution, additional requirements are applied in selecting D 0 → K − π + candidates with exactly two charged tracks. We veto these events if they satisfy at least one of the following conditions: TOF information consistent with a cosmic ray event, particle identification information consistent with an e + e − hypothesis, two tracks with EMC energy deposits consistent with an e + e − hypothesis, or either track with particle identification and MUC information consistent with being a muon.

Yields and efficiencies
The M BC distribution for single-tag candidates for each mode is fitted with a MC-derived signal shape and an ARGUS function background [17]. The signal shape is convolved with a double Gaussian with a common mean to allow for differences in M BC resolution between data and MC. Charge-conjugate modes are fitted simultaneously with the double-Gaussian signal-shape parameters constrained to be the same and the normalizations and background parameters allowed to vary independently in the fit. Peaking backgrounds contributed by decay modes that have similar final states to the signal mode are included in the signal shape, although the yields are corrected after the fit to count only true signal events.
An example M BC fit is shown in Fig. 2. (The full set of fits is provided in Appendix A.) In events with multiple single-tag candidates, the best candidate is chosen per mode and per charm to be the one with the smallest |∆E|. Based on the fit results tight mode-dependent requirements on ∆E are applied. To determine the tag yield, the M BC histogram is integrated within the signal region, 1.8580 GeV/c 2 M BC 1.8740 GeV/c 2 for D 0 modes and 1.8628 GeV/c 2 M BC 1.8788 GeV/c 2 for D + modes, and then the analytic integral of the ARGUS function in this region is subtracted. The efficiency for each of the 18 single-tag modes is found by using MC truth information to determine the total number generated for the denominator and using the same cut-andcount method as used for data to determine the numerator. The single-tag yields and efficiencies are summarized in Table 2, where the efficiencies include branching fractions for π 0 →γγ and K 0 S →π + π − decays. MBC fit for single-tag mode D + → K − π + π + π 0 , from data. Blue dash-dot (green dashed) line represents the total fit (the fitted background shape) and the red solid curve corresponds to the fitted signal shape. Table 2. Single-tag yields after subtracting their corresponding peaking backgrounds from data and efficiencies from MC, as described in the text. The uncertainties are statistical only.
tag mode yield efficiency(%) D 0 →K − π + 260,915 ± 520 63.125 ± 0.007 D 0 →K + π − 262,356 ± 522 64.272 ± 0.006 D 0 →K − π + π 0 537,923 ± 845 35.253 ± 0.007 D 0 →K + π − π 0 544,252 ± 852 35.761 ± 0.007 D 0 →K − π + π + π − 346,583 ± 679 38.321 ± 0.007 D 0 →K + π + π − π − 351,573 ± 687 39.082 ± 0.007 D + →K − π + π + 391,786 ± 653 50.346 ± 0.005 D − →K + π − π − 394,749 ± 656 51.316 ± 0.005 D + →K − π + π + π 0 124,619 ± 529 26.138 ± 0.014 Double tags are fully reconstructed events in which both the D and theD pass the selection criteria for one of the tag modes. In events with multiple double-tag candidates, the best candidate per mode combination per event is chosen with the [M BC (D)+M BC (D)]/2 closest to the known D mass. Following a procedure similar to the single-tag counting, we fit the two-dimensional distribution of M BC (D) vs. M BC (D) for the selected single-tag modes to define the signal region for a cut-and-count determination of the double-tag yield. A more sophisticated treatment of the background is required because of the correlations between the tags. The signal shape is again derived from MC, using truth information and including peaking backgrounds with the signal. We found that convolving the MC shape with smearing functions to account for the small data/MC resolution difference did not appreciably improve the accuracy of the tag yields, so no signal smearing is included in the double-tag fits.
The background shapes in the double-tag fits correspond to four possible ways of mis-reconstructing an event, as shown in Fig. 3. A direct product of a MCderived signal shape with an analytic ARGUS function background, with shape parameters fixed to those of the corresponding single-tag fit, is used to represent the background contributed by events with a correctly reconstructed D and incorrectly reconstructedD. The background shape for the charm-conjugate case is similarly constructed. For completely reconstructed continuum events or fully reconstructed but mispartitioned DD events (with particles assigned incorrectly to the D andD), a direct product of a double-Gaussian function and an ARGUS function rotated by 45 • is used. The kinematic limit and exponent parameters of the rotated ARGUS function are fixed, while the slope parameter is allowed to be free in the fit. Finally, the remaining background events with neither D norD correctly reconstructed are modeled with a direct product of two AR-GUS functions, with parameters taken from the corresponding single-tag fits. An example fit to data is shown in    Fig. 4. (color online) Example two-dimensional MBC double-tag fit from data as described in the text, for tag mode K + π − π − vs. K − π + π + π 0 . The top left figure is a scatter plot of the data and the top right is a scatter plot of the fit to the data. The bottom two plots are overlays of data and the fit projected onto the positive and negative charm MBC axes. The red dashed (blue solid) lines represent the total fits (the fitted signal shapes) and the solid green curves are the fitted background shapes. The magenta curve corresponds to the case when D − → K + π − π − is reconstructed correctly, while D + →K − π + π + π 0 is not.
After the two-dimensional fit is performed, the M BC histogram is integrated within the same signal region as the single-tag fits, and the integrals of the four background shapes are subtracted from this total. The resultant double-tag yields and efficiencies, which include branching fractions for π 0 →γγ and K 0 S →π + π − decays, are summarized in Tables 3 and 4.
We must correct the yields determined with the M BC fits (data and MC) for contributions from background processes that peak in the signal region. Such backgrounds come from other D decays with similar kinematics and particle compositions as the specific signal mode. We rely on MC, generated with world-average branching fractions [1], to determine the fraction of peaking background events, as well as to calculate their selection efficiencies. We apply MC-determined corrections for these in every case where more than 0.01% of the fitted yield is attributable to peaking background. The largest contribution of peaking background is for D + → K 0 S π + π + π − , approximately 2.5% of the fitted yield. D 0 → K − π + π + π − and D + → K 0 S π + π 0 both have ∼ 2.0% of their fitted yields from peaking backgrounds, and all other modes have less than 1.0%. Because the peaking backgrounds come from well understood processes, like doubly Cabibbo-suppressed modes, simultaneous misidentification of both a pion and a kaon in an event, and charged pion pairs not from K 0 S decays that pass the K 0 S invariant mass requirement, we are confident that they are well modeled by the MC.
The analysis described above results in a set of measured values of N DDij , the number of DD events determined with the single-and double-tag yields of positive tag mode i and negative tag mode j. The uncertainties are highly mode dependent because of branching fractions, efficiencies and backgrounds, so these measurements must be combined into an uncertainty-weighted mean taking into account correlations within and between the mode-specific measurements. We use an analytic procedure for this and demonstrated its reliability with a toy MC study.
For our full 2.93 fb −1 ψ(3770) data sample we find N D 0D0 =(10,621±29)×10 3 and N D + D − =(8,296±31)×10 3 . Using the integrated luminosity from Ref. [8], we obtain observed cross sections for DD production at the 083001-8 Table 3. D 0D0 double-tag yields from data and efficiencies from MC, as described in the text. The uncertainties are statistical only.

Effects of quantum correlations
As mentioned earlier in this paper, the D 0D0 yield and cross section must be corrected for correlations introduced by production through a pure C =−1 state at the ψ(3770). Asner and Sun [6] provide correction factors that can be applied directly to our measured yields with Eq. 5 for D 0 →f andD 0 →f and Eq. 6 for the case f =f .
The quantities appearing in these equations can be expressed in terms of measured parameters of D 0 decays and D 0D0 mixing, with v − jk =(z j z k −w j w k )/2, z j =2cosδ j and w j =2sinδ j . r j and δ j are defined by j|D 0 / j|D 0 = r j e iδ j , where r j = | j|D 0 / j|D 0 |, and δ j is the average strong phase difference for the Cabibbo-favored tag mode. The usual mixing parameters x and y, which are related to the differences in masses and lifetimes of the two mass eigenstates, enter throughỹ j =ycosδ j +xsinδ j . The D 0 → K − π + π 0 and K − π + π + π − tag modes require a slightly more complicated treatment because they are mixtures of modes with different phases. This requires introducing coherence factors R j to characterize the variation of δ j over phase space, with z j and w j being redefined as z j =2R j cosδ j and w j =2R j sinδ j [18]. Table 5 shows the input parameters that are used to obtain the correction factors and Fig. 5 shows the corrections to σ(e + e − →D 0D0 ) for each of the nine double-tag modes, along with the average. The overall effect is a relative change in N D 0D0 of approximately −0.2%, with final corrected values of N D 0D0 = (10,597±28)×10 3 and σ(e + e − → D 0D0 ) = (3.615±0.010) nb. The uncertainties are statistical only. The summed χ 2 value relative to the mean for all pairs of tag modes is 11.8 for D 0D0 (9 modes).

Systematic uncertainties
The sources of systematic uncertainty that have been considered for the D 0D0 and D + D − cross section measurements are listed in Table 6.
The double-tag technique used to determine the event yields and cross sections σ(e + e − →D 0D0 ) and σ(e + e − → D + D − ) has the benefit of substantial cancellation of systematic uncertainties. Detector effects including tracking, particle identification, and π 0 and K 0 S reconstruction, along with tag-mode resonant substructure and the ∆E requirement, all affect both single and double tags.

083001-10
There are, however, event-dependent effects that do not cancel in the efficiency ratio ij /( i · j ). The event environment in which D mesons are tagged affects the efficiency because higher multiplicities of charged tracks or π 0 s lower the tagging efficiency. This can arise due to three possible sources: (1) differences in multiplicitydependent efficiencies between data and MC, (2) differences between the other-side multiplicities in data and MC due to imperfect knowledge of D meson decay modes and rates, and (3) sensitivity of the best candidate selection to the number of fake-tag background events. To assess a possible uncertainty due to the first source, we study efficiencies of tracking and particle identification for charged pions and kaons, as well as π 0 reconstruction, based on doubly tagged D 0D0 and D + D − samples. We estimate uncertainties while observing how well our MC simulates these efficiencies in data with different particle multiplicities.
We evaluate the effect of the second source for both tracks and π 0 s by reweighting the MC to better match the multiplicities in data. In this we assume that data and MC are consistent in the single track and π 0 reconstruction efficiencies. We obtain corrected efficiencies separately for each tag mode, and the difference with the nominal efficiency is used as the systematic uncertainty. The effect is larger for tag modes with greater multiplicity, and so the overall effect on D + D − is greater than that on D 0D0 .
The third source arises due to the fact that we resolve multiple-candidate events when choosing single tags based on the smallest |∆E|. This selection is imperfect and sometimes the wrong candidate is chosen, lowering the efficiency for multiple-candidate events relative to single-candidate events. Although a bestcandidate selection is also applied to double tags, the number of multiple candidates in this case is small and the selection based on two beam-constrained masses is more reliable, so only the systematic uncertainty of best-candidate selection for single tags is considered. Such uncertainty only arises when both the multiple-candidate rate is different between data and MC and the singleand multiple-candidate efficiencies are different. These quantities can be measured both in data and MC, and the observed differences are propagated through to the systematic uncertainties in the cross sections.
Even though we fit both single and double tags to obtain the yields and efficiencies, the differences between one-and two-dimensional fits and the much lower background levels of the double-tag M BC distributions limit the cancellation. We consider several variations of the fitting procedures and use the changes in efficiencycorrected yields to estimate the systematic uncertainties.
The uncertainty due to the single-tag background shape is probed by substituting a MC-derived background for the ARGUS function. The uncertainty due to the signal shape is assessed by altering the smearing of the MC-derived shape (single-Gaussian-convolved instead of the double-Gaussian-convolved). To assess the uncertainty in the double-tag fitting procedure, we obtain double-tag yields and efficiencies with an alternative sideband-subtraction method, dividing the twodimensional M BC plane into sections representing the signal and various background components, as shown in Fig. 3. The signal area is the same as that used when fitting. Horizontal and vertical bands are used to represent combinations with one correctly and one incorrectly reconstructed D; a diagonal band represents the background from completely reconstructed continuum events or mispartitioned DD events; and two triangles are used to represent the remaining background, which is mostly flat. An estimate of the flat background is scaled by the ratios of the sizes of each of the other background regions and subtracted to obtain estimates of the non-flat backgrounds. These backgrounds are then scaled with area and ARGUS background parameters obtained from single-tag fits to determine the overall background subtraction and yield for the signal region for a specific tag mode. The difference in efficiency-corrected double-tag yields for each mode between this method and the standard procedure is taken as the systematic uncertainty associated with the double-tag fitting method.
The cosmic and lepton veto suppresses cosmic ray and QED background in the single-tag selection for the D 0 → K − π + mode. A cosmic ray background event is produced by a single particle that is incorrectly reconstructed as two oppositely charged tracks. The net momentum of the two tracks is therefore close to zero, and typical QED events also have small net momentum. This small momentum produces M BC values close to the beam energy, so that residual cosmic ray and QED events passing the veto distort the M BC distribution. Because the processes responsible are not included in our MC samples or well described by the ARGUS background function, the fit results may be affected. To assess this effect, we performed alternative single-tag fits for D 0 →K − π + with a cut-off in M BC at 1.88 GeV/c 2 , excluding the range where cosmic and QED events can contribute. We found the resulting difference from the standard fit procedure to be 0.18%, which we take as the systematic uncertainty due to this effect.
The line shape of the ψ(3770) affects our analysis through the modeling of initial-state radiation (ISR) at the peak of the resonance. The cross section for ψ(3770) production in radiative events depends on the cross section value at the lower effective E cm that results from ISR. While this may partially cancel in the ratio, we treat it separately for single and double tags because yields and efficiencies are affected with opposite signs, and because correlations are introduced for the doubletag fits that are not present in the single-tag fits. The MC-determined efficiencies are affected through the ∆E requirements, which select against large ISR because the ∆E calculation assumes that the energy available to the D is the full beam energy. The data yields are affected via the M BC fit shape, which acquires an asymmetric high-side tail through the contribution of ψ(3770) production via ISR. More ISR causes a larger high-side tail in both the single-and double-tag signal shapes. Additionally, because both D mesons lose energy when ISR occurs, double-tag events that include ISR will have a correlated shift in M BC , causing such events to align with the diagonal to the high-side of the signal region in the two-dimensional M BC plane. We use a preliminary BE-SIII measurement of the ψ(3770) line shape to re-weight the MC and repeat the D-counting procedure. Combining the mode-by-mode variations in N DD leads to the systematic uncertainty associated with the ψ(3770) line shape given in Table 6.
The MC modeling of final-state radiation (FSR) may lead to a systematic difference between data and MC tag-reconstruction efficiencies. FSR affects our measurement from the tag-side, so any systematic effect will also have some cancellation. To assess the uncertainty due to FSR we created signal MC samples with and without modeling of FSR and measured the changes in tag reconstruction efficiencies. The largest difference was for D 0 → K − π + , where the relative change in single-tag reconstruction efficiency was 4%. The D 0 → K − π + ,D 0 → K + π − double-tag reconstruction efficiency also changed when FSR was turned off, but the cancellation was not complete, with the ratio of efficiencies changing by 1.2%. Because the variation of turning on and off FSR modeling is judged to be too extreme (FSR definitely happens), we take 25% of this difference as our systematic uncertainty due to FSR modeling, a 0.3% relative uncertainty on the MC reconstruction efficiency ratio. To be conservative, we take the largest change, for the D 0 →K − π + mode, as the systematic uncertainty for all modes.
The correction in the D 0D0 cross section due to the treatment of quantum correlations incurs systematic uncertainty associated with the parameters x, y, δ Kπ , and r 2 Kπ , for which Ref.
[19] provides correlation coefficients. Ref. [20] provides a similar coefficient table for the rest of the variables. In evaluating our systematic uncertainty, we have doubled the reported uncertainties and treated them incoherently. Toy MC calculations were used to propagate these uncertainties to N D 0D0 , giving a systematic uncertainty in the D 0D0 cross section of 0.2%.
Finally, for the calculation of cross sections, the relative systematic uncertainty due to the integrated luminosity measurement is determined in Ref. [7,8] to be 0.5%.

Results and conclusions
The separate sources of systematic uncertainty given in Table 6 are combined, taking correlations among them into account, to give overall systematic uncertainties in the D 0D0 and D + D − cross sections of 1.05% and 0.93%, respectively. Including these systematic uncertainties, the final results of our analysis are as follows: where the uncertainties are statistical and systematic, respectively. In the determinations of σ(e + e − → DD) and σ(e + e − → D + D − )/σ(e + e − → D 0D0 ), the uncertainties of the charged and neutral cross sections are mostly uncorrelated, except the systematic uncertainties due to the assumed ψ(3770) line shape, the FSR simulation, and the measurement of the integrated luminosity.
In conclusion, we have used 2.93 fb −1 of e + e − annihilation data at the ψ(3770) resonance collected by the BESIII detector at the BEPCII collider to measure the cross sections for the production of D 0D0 and D + D − . The technique is full reconstruction of three D 0 and six D + hadronic decay modes and determination of the number of D 0D0 and D + D − events using the ratio of singletag and double-tag yields. We find the cross sections to be σ(e + e − → D 0D0 ) = (3.615±0.010±0.038) nb and σ(e + e − → D + D − ) = (2.830±0.011±0.026) nb, where the uncertainties are statistical and systematic, respectively. These results are consistent with and more precise than the previous best measurement by the CLEO-c Collaboration [4] and are necessary input for normalizing some measurements of charmed meson properties in ψ(3770) decays.