Determination of the $C\!P$-even fraction of $D^0\rightarrow K_S^0\pi^+\pi^-\pi^0$

Quantum-correlated $D\bar{D}$ pairs collected by the BESIII experiment at the $\psi(3770)$ resonance, corresponding to an integrated luminosity of 2.93 fb$^{-1}$, are used to study the $D^0 \rightarrow K^{0}_S\pi^{+} \pi^{-} \pi^{0}$ decay mode. The $C\!P$-even fraction of $D^0 \rightarrow K^{0}_S\pi^{+} \pi^{-} \pi^{0}$ decays is determined to be $0.235\pm 0.010\pm 0.002$, where the first uncertainty is statistical and the second is systematic.

a Also at the Moscow Institute of Physics and Technology, Moscow 141700, Russia b Also at the Novosibirsk State University, Novosibirsk, 630090, Russia c Also at the NRC "Kurchatov Institute", PNPI, 188300, Gatchina, Russia Quantum-correlated D D pairs collected by the BESIII experiment at the ψ(3770) resonance, corresponding to an integrated luminosity of 2.93 fb −1 , are used to study the D 0 → K 0 S π + π − π 0 decay mode.The CP -even fraction of D 0 → K 0 S π + π − π 0 decays is determined to be 0.235 ± 0.010 ± 0.002, where the first uncertainty is statistical and the second is systematic.

I. INTRODUCTION
The complex phase of the Cabibbo-Kobayashi-Maskawa (CKM) matrix is the only source of the CP violation within the quark sector of the Standard Model (SM) [1].However, this source of CP violation is too weak to explain the huge matter-antimatter asym-metry in the universe.This fact motivates precise studies of flavor transitions to test the CKM paradigm and to search for evidence of CP violation beyond the SM.The unitarity triangle is a graphical representation of the CKM matrix in the complex plane, with angles denoted α, β, and γ [2].An improved measurement of the angle γ is of particular importance in studies of CP violation.Comparisons between the values of γ obtained from direct measurements [2] and global CKM fits [3] provide an important test of CKM unitarity and allow for searches for new physics beyond the SM.
The angle γ of the unitarity triangle is measured in decays which are sensitive to interference between favored b → c and suppressed b → u quark transition amplitudes [4].Typically, the interference is measured in B ± → Dh ± decays, where D is an admixture of D 0 and D0 flavor states and h ± is either a charged pion or kaon.The theoretical uncertainty of the measurement of γ through this approach is negligible [5].The current uncertainty of the γ measurement is statistically dominated [4].Therefore, including more D-decay modes is desirable to improve the precision of the γ measurement.One of the important classes of B ± → Dh ± measurement is the so-called 'GLW' strategy, in which the D decay is reconstructed in a CP eigenstate [6].This strategy can be extended to encompass quasi-CP eigenstates that are predominantly CP odd or even [7].This approach requires that the CP -even fraction, quantified by the parameter F + , in known.A promising quasi-CP eigenstate is the decay D → K 0 S π + π − π 0 , which has a higher branching fraction than decays to individual CP eigenstates and is known to be predominantly CP odd from a measurement performed with data collected by CLEO-c corresponding to an integrated luminosity of 0.82 fb −1 that yielded the result F + = 0.238 ± 0.012 ± 0.012 [8].In this paper we present an improved determination of this parameter made with quantum-correlated D D pairs, based on a data sample of 2.93 fb −1 taken at the ψ(3770) resonance by the BESIII detector.

II. FORMALISM
The wave function for the pairs of neutral charm mesons produced at the ψ(3770) resonance is antisymmetric because of their odd charge conjugation.This quantum correlation allows the CP -even fraction F + of D 0 → K 0 S π + π − π 0 decays to be determined with doubletag (DT) [9] yields, in which the D meson is reconstructed in the signal decay and the D meson is reconstructed in one of several tag modes [10].The tag modes used in this analysis are summarized in Table I.
The expected DT yield is given by [7] where N D D is the number of D 0 D0 pairs [11] produced at the ψ(3770) resonance, S (T ) denotes the signal mode ] is the branching fraction of the D 0 decaying into the S (T ), ε(S|T ) is the DT efficiency, F + (F T + ) is the CP -even fractions for the D → S (D → T ) decays.In addition, Eq. ( 1) omits terms of order O(x 2 D , y 2 D ) ∼ 10 −5 associated with D 0 D0 mixing effects, where x D and y D are the mixing parameters [2].The dependency on the N D D , B(T ), and single-tag (ST) efficiency of the tag mode is removed by measuring the ST yield of the tag mode T .The expected ST yield for the tag mode T is given by [7] where ε(T ) is the ST efficiency.A further correction of 1/ 1 − (2F T + − 1)y D is applied to N ST to account for D 0 D0 mixing effects, where y D is one of the D 0 D0 mixing parameters [2].With assumption of ε(S|T ) = ε(S) • ε(T ), the ratio of N DT to N ST can be written as For CP -even (CP -odd) tags, F T + is equal to 1 (0).Then the corresponding ratio R − (R + ) for the CP -tag mode T is given by where η + T (η − T ) is the CP eigenvalue of the CP -even (CPodd) tag mode T .With such CP -eigenstate tags, the CP -even fraction is given by which has no dependence on B(S) and ε(S).For a quasi-CP -tag mode T of known CP -even fraction F T + , the CPeven fraction of the signal mode is where R T is the ratio of the measured DT yield to the measured ST yield for the tag mode T .The signal decay D → K 0 S π + π − π 0 can also be used as a self-tag mode.Denoting the ratio of the DT yield of the self-tag mode to the corresponding ST yield by R S , the CP -even fraction is given by

Type
Tag modes CP -even There is no direct measurement of F + for the mixed CP -tag modes D → K 0 S,L π + π − .However, quantumcorrelated studies indicate that F K 0 S,L π + π − + is close to 0.5 for both decay modes [12].It follows from Eq. (1) that using these decays integrated over all phase space as tags gives very low sensitivity to F + of the signal mode.Therefore a localized measurement is pursued, in which the DT yields are measured in the pairs of eight bins in the Dalitz plot of D → K S,L π + π − .The binning scheme adopted is the "equal ∆δ binning" of Ref. [13], as shown as Fig. 1.The expected DT yield with the K 0 S π + π − final state in the i-th bin can be written as [12] and the corresponding expected yield for the K 0 L π + π − case as where for the decay ) with the final state in the i-th bin, and c i (c i ) is the amplitude-weighted average of cosine of the strong-phase difference [12].An event can migrate between the bins because of the finite detector resolution.To improve the momentum resolution of the final-state particles and suppress the migration between phase-space bins, the invariant mass of the K 0 S π + π − final state is constrained to the known D 0 mass [2], denoted as M D PDG .

III. DETECTOR AND DATA SAMPLES
The BESIII detector [14] records symmetric e + e − collisions provided by the BEPCII storage ring [15] in the  center-of-mass energy range from 2.0 to 4.95 GeV, with a peak luminosity of 1 × 10 33 cm −2 s −1 achieved at centerof-mass energy √ s = 3.77 GeV.BESIII has collected large data samples in this energy region [16].The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator timeof-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field [17].The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identification modules interleaved with steel.The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for electrons from Bhabha scattering.The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end-cap) region.The time resolution in the TOF barrel region is 68 ps, while that in the end-cap region is 110 ps.
Simulated data samples produced with a geant4based [18] MC package, which includes the geometric description of the BESIII detector and the detector response, are used to determine detection efficiencies and to estimate backgrounds.The simulation models the beam-energy spread and initial-state radiation (ISR) in the e + e − annihilations with the generator kkmc [19].The inclusive MC sample includes the production of D D pairs, the non-D D decays of the ψ(3770), the ISR production of the J/ψ and ψ(3686) states, and the continuum processes incorporated in kkmc [19].All particle decays are modeled with evtgen [20] using branching fractions taken from the PDG [2].Final-state radiation from charged final state particles is incorporated using the photos package [21].For each tag mode, signal MC samples for ST and DT are produced.The multibody decay is generated with available amplitude modes or the expected intermediate resonance involved.

IV. EVENT SELECTION
Charged tracks detected in the MDC are required to be within a polar angle (θ) range of |cos θ| < 0.93, where θ is defined with respect to the z-axis, which is the symmetry axis of the MDC.For charged tracks not originating from K 0 S decays, the distance of the closest approach to the interaction point must be less than 10 cm along the zaxis, |V z |, and less than 1 cm in the transverse plane, |V xy |.Particle identification (PID) for charged tracks combines measurements of the specific ionization energy loss in the MDC (dE/dx) and the flight time in the TOF to form likelihoods L(h) (h = K and π) for each hadron h hypothesis.Charged kaons and pions are identified by comparing the likelihoods for kaon and pion hypotheses, L(K) > L(π) and L(π) > L(K), respectively.
Photon candidates are identified using showers in the EMC.The deposited energy of each shower must be more than 25 MeV in the barrel region (|cos θ| < 0.80) and more than 50 MeV in the end-cap region (0.86 < |cos θ| < 0.92).To suppress electronic noise and showers unrelated to the event, the difference between the EMC time and the event start time is required to be within [0, 700] ns.Each K 0 S candidate is reconstructed from two oppositely charged tracks satisfying |V z | < 20 cm.The two charged tracks are assigned as π + π − without imposing further PID criteria.They are constrained to originate from a common vertex and are required to have an invariant mass within [0.485, 0.510] GeV/c 2 .The decay length of the K 0 S candidate is required to be greater than twice the vertex resolution away from the interaction point.Pairs of selected photon candidates are used to reconstruct π 0 (η) candidates.The invariant masses of the photon pairs are required to be within [0.110, 0.155] ([0.480, 0.580]) GeV/c 2 .To improve the momentum resolution, a kinematic fit is applied to constrain the γγ invariant mass to the known π 0 (η) mass [2], and the χ 2 of the kinematic fit is required to be less than 20.The fitted momenta of the π 0 (η) are used in the subsequent analysis.The ω candidates are selected by requiring the invariant mass of the π + π − π 0 combination to be within [0.750, 0.820] GeV/c 2 .The η → π + π − η and η → γπ + π − decays are used to reconstruct η mesons, with the invariant masses of the π + π − η and π + π − γ required to lie within [0.938, 0.978] GeV/c 2 .
For the D → K + K − and D → π + π − tag modes, the background from cosmic rays and Bhabha events are suppressed with the following requirements [22].First, the two charged tracks used as the CP tag must have a TOF time difference less than 5 ns and they must not be consistent with being a muon pair or an electron-positron pair.To be a muon candidate, the track must have |χ dE/dx | < 5, 0.15 < E EMC < 0.30 GeV, and d µ > 40 or d µ > 80p − 60, where χ dE/dx [23] is the χ for the hypothesis of muon in the PID with dE/dx, E EMC is deposited energy in the EMC by the track, d µ is the penetration depth in the muon identification modules with unit of center meter, and p is the momentum of the track.To be an electron (positron) candidate, the track must have greater probability (combing the dE/dx and flight time) being an electron (positron) than being a kaon or pion.If the corresponding EMC shower is reconstructed, the track must satisfy: (1) E/p > 0.8 if |cos θ| < 0.70 , or (2) E/p > −7.5|cos θ| + 6.05 if 0.7 < |cos θ| < 0.8, or (3) E/p > 0.6 if |cos θ| > 0.85, where E/p is the ratio of the deposited energy in the EMC to the momentum of the track and θ is the polar angle of the shower.Second, there must be at least one EMC shower (other than those from the CP -tag tracks) with an energy larger than 50 MeV or at least one additional charged track detected in the MDC.For the D → π + π − π 0 and D → π + π − π + π − tag modes, events are rejected if any combination of m π + π − lies in the mass window of [0.481, 0.514] GeV/c 2 in order to suppress the background from K 0 S decays.

V. DETERMINATION OF THE ST YIELD
The tag modes that do not involve a K 0 L are fully reconstructed.The signal candidate of the fully reconstructed tag mode is identified by the beam-energy constrained mass where E beam and p D are the beam energy and reconstructed momentum of the D candidate in the e + e − center-of-mass frame, respectively.If multiple combinations are reconstructed for an event, the combination with the lowest value of |∆E| is retained for further analysis, where ∆E is the difference between the E beam and the reconstructed energy.To suppress combinatorial background, the value of ∆E for each event is required to be within the ±3σ ∆E around the peak of the ∆E distribution, where σ ∆E is the resolution of the distribution.
Figure 2 shows the M BC distributions of the ST D candidates for the fully reconstructed tag modes.The ST yields for the tag modes are obtained by an unbinned maximum likelihood fit to the M BC distributions.In the fit, the signal component is described by the MC-simulated shape.To account for the difference in resolutions between the MC simulation and data, the MC-simulated shape for each mode is convolved with a Gaussian function with free parameters.The background component is modeled with an AR-GUS function [24], where the slope parameter is a free parameter and the endpoint is fixed to the beam energy [25] in the center-of-mass frame.The fit results are shown in Fig. 2. In addition to the combinatorial background, there are also peaking backgrounds that have similar final states as the signal and are included in the signal yield.The contamination rate of the peaking background is estimated by analyzing the inclusive MC sample.For the tag mode D → K 0 S π 0 (K 0 S π 0 π 0 ), the background level is 4.7% (0.4%) dominated by D → π + π − π 0 (π + π − π 0 π 0 ) decays.
The D → K 0 L π 0 and D → K 0 L ω tag modes cannot be fully reconstructed.Therefore, the effective ST yields are calculated from Eq. ( 2).The effective ST efficiencies for the D → K 0 L π 0 and D → K 0 L ω tags, which are required inputs for Eq. ( 2), are defined to be the ratios of the corresponding DT efficiencies to the efficiency for the ST of the signal mode.The effective ST yields are also summarized in Table II

VI. DETERMINATION OF THE DT YIELD
The DT candidates for the fully reconstructed tag modes are isolated with the beam-constrained masses for the signal and tag modes, denoted as M S BC and M T BC , respectively.In the case of multiple combinations, the combination with the least |M S BC + M T BC − 2M D PDG | is retained for further analysis.The ∆E variables for the signal and tag modes are required to lie within the regions given in Sec.V. To determine the DT yield, a two-dimensional (2D) unbinned maximum likelihood fit is performed to the M S BC versus M T BC distribution.An example of the M T BC versus M S BC distribution for the D → π + π − π 0 tag mode is shown in Fig. 3.The signal component in the fit is described by a 2D MCsimulated shape obtained from the signal MC sample smeared with a Gaussian function in each dimension.The Gaussian function is introduced to account for the difference in resolutions between the MC simulation and data.The mass and width of the Gaussian function for each tag mode (the signal mode) are obtained from the one-dimensional fit in the ST-yield determinations for the corresponding tag mode (the signal mode).The back-ground component with the correctly reconstructed signal mode (tag mode) and incorrectly reconstructed tag mode (signal mode) is modeled by the product of the signal and background shapes from the fits for the ST yield determinations of signal mode (tag mode) and tag mode (signal mode), respectively.The parameters of the shapes are fixed at the values obtained from the corresponding one-dimensional fit.The background component where signal and tag mode are both reconstructed incorrectly is modeled by the product of the background shapes from corresponding fits for the ST yield determinations.The backgrounds involving swapped final-state particles from the two charm mesons and continuum processes, corresponding to the diagonal band in Fig. 3, are modeled by the product of a Gaussian function and the ARGUS function rotated by 45° [26].The endpoint of the ARGUS function is fixed at the beam energy in the e + e − center-of-mass frame.Figure 4 shows the projections of 2D fits on the M S BC distributions for the fully reconstructed tag modes.According to studies performed with the inclusive MC sample, there are small contributions from peaking backgrounds for each mode.The dominant components for the signal and tag modes are the same as those in the determination of ST yields.The  peaking background yields are determined by analyzing the inclusive MC sample and are corrected for the quantum correlation with Eq. ( 1).For the D → K 0 S π + π − tag mode, the DT yields in the eight pairs of bins are determined with the same method as described above.The fitted results are shown in Fig. 5.
The DT candidates tagged by the D → K 0 L π 0 and D → K 0 L ω tag modes cannot be fully reconstructed.They are selected by the variable M 2 miss defined as where i E i is the sum of the reconstructed energies of the tag mode and i p i is the sum of the reconstructed momenta of the signal mode and tag mode.The distribution from correctly reconstructed DT candidates peaks around the squared mass of the K 0 L meson.To suppress background, events with excess reconstructed charged tracks or π 0 candidates are rejected.The DT yields are determined by fitting to the M 2 miss distribution.In this fit, the signal is described by an MC-simulated shape, which is convolved with a Gaussian with free parameters.The combinatorial background is modeled with a second-order Chebyshev function.The yield of peaking background for the D → K 0 S π 0 (D → K 0 S ω) decay is estimated with the inclusive MC sample and is corrected for quantum correlations.The peaking background in the D → K 0 L π 0 (D → K 0 S ω) tag mode originates from D → ηπ 0 (D → ηK 0 S π 0 ) decays and is modeled by an MC-simulated shape with its yield fixed according to the results from the study of the inclusive MC sample.For the D → K 0 L ω tag mode, the non-resonant background from D → K 0 S,L π + π − π 0 is estimated with the ω sideband events in the M π + π − π 0 distribution.Figure 6 shows the fit results for the two tag modes.The DT yields for the D → K 0 L π + π − tag mode in the eight pairs of bins are determined with the same method.The corresponding fit results in the eight pairs of bins are shown in Fig. 7.

VII. SYSTEMATIC UNCERTAINTIES ON THE YIELD DETERMINATIONS
A. The ST yields The ST yields for the fully reconstructed tag modes are determined by fitting the M BC distributions after subtracting the peaking background yields estimated from the inclusive MC sample.The uncertainty associated with the fit is estimated by floating the endpoint, which is fixed in the baseline fit, of the ARGUS function.The difference in the yield to that of the baseline fit is taken as the systematic uncertainty.The uncertainties for the different tag modes lie in the range of [0.1, 0.3]%.The uncertainties on the peaking background yields are estimated by varying the quoted branching fractions [2] by ±1σ and range from 0.1% to 0.5%.
For the D → K 0 L π 0 and D → K 0 L ω tag modes, the uncertainties for the effective ST yields calculated as Eq. ( 2) are associated with N D D and the product of the branching fractions and efficiencies of the two tag modes.The uncertainty of N D D which is 1.0% has been estimated in the measurement of D D cross section [11].The uncertainties of the products of the branching fractions and efficiencies for the D → K 0 L π 0 and D → K 0 L ω tag modes have been estimated in the branching fraction measurements of the two tag modes [27,28].They are 3.1% and 2.6% for the D → K 0 L π 0 and D → K 0 L ω tag modes, respectively.These uncertainties are propagated to the effective ST yields.

B. The DT yields
The DT yields for the fully reconstructed tag modes are determined by the 2D fits to the M S BC versus M T BC distributions after subtracting the estimated peaking background after correcting for the effects of quantum correlations.Since the fit strategies are the same for all tag modes, the largest DT sample, which is that involving D → π + π − π 0 tags, is adopted to estimate the uncertainties introduced by the fits in order to minimize the effects of statistical fluctuations.For the uncertainty arising from the signal models, the 2D MC-simulated shape without a smearing Gaussian resolution function is taken as the alternative model.The change of the signal yield, 0.34%, is taken as the uncertainty.For the background shapes, the fixed parameters of the ARGUS functions are changed to free parameters in the fit.The change in the signal yield, which is 0.13%, is assigned as the corresponding uncertainty.The uncertainties due to the peaking background subtraction with the inclusive MC sample are estimated with the uncertainties of the corresponding branching fractions [2] and F + for quantum correlations [7,8,12,29].The estimated uncertainties are in the range of [0.1, 1.0]%.The propagated systematic uncertainties for the tag modes are listed in Table II.
For the D → K 0 L π 0 and D → K 0 L ω tag modes, the yields are determined from the fits to the M 2 miss distribution after subtracting the peaking background yield estimated from the inclusive MC sample and data studies.The systematic uncertainties from the fits to the M 2 miss distributions have several components.To assess that uncertainty coming from the background shape, the secondorder Chebyshev function is replaced by a third-order Chebyshev function.The resulting differences in the yields, 1.1% for the D → K 0 L π 0 tag mode and 0.17% for the D → K 0 L ω tag mode, are assigned as the uncertainties from the fit procedure.The uncertainties from the peaking background subtractions are estimated by varying the assumed branching fractions [2] and F + for the quantum correlation corrections [7,8,12,29] by ±1σ.For the D → K 0 L ω tag mode, the contribution of the peaking background from non-resonant D → K 0 S,L π + π − π 0 final state is estimated by fitting the M 2 miss distribution from the ω sidebands.The uncertainty on the fitted peaking background yield is taken as the uncertainty from the sample size of sideband events.The sidebands are altered to estimate the uncertainty due to the choice of sidebands.The estimated uncertainties from the peaking background subtractions are 0.32% and 2.1% for the D → K 0 L π 0 and D → K 0 L ω tag modes, respectively.The combined systematic uncertainty of the DT yield for each tag mode is summarized in Table II.

VIII. THE F+ MEASUREMENT A. Measurement with the CP -tag modes
The expected ratio of the DT yield to the corresponding ST yield is calculated from Eq. ( 4).Implicity in this expression is the assumption of ε(S | T ) = ε(S) • ε(T ).However, studies of the signal MC samples indicate that this assumption is not always true.Therefore, a correction factor of ε(S) • ε(T )/ε(S|T ) is applied to the measurement of R ± for each tag mode.Figure 8 shows the measured R ± values for each tag mode after applying this correction.The mean values of R + and R − are determined by least χ 2 fits.The χ 2 ± for R ± in the fit is constructed as follows: where R ± is the mean value of R ± , R ± i (R ± j ) is the ratio of the DT yield to the corresponding ST yield of the i-th (j-th) tag mode, and V ± ij is the covariance between modes i and j.The measured values of R ± for the different tag modes are independent, except for the D → K 0 L π 0 and D → K 0 L ω tags, where there is a correlation coefficient of 0.02 introduced by the common use of N D D .The fitted result for R ± is shown as the yellow bands in Fig. 8.These results from the CP -tag modes lead to a value of F + = 0.229 ± 0.013, where the uncertainty includes both statistical and systematic contributions.Re-performing the fit with only the statistical uncertainties included on the inputs allows the statistical and systematic contributions on the fit uncertainty to be isolated, with the statistical uncertainty found to be 0.013 and the systematic uncertainty to be 0.001.
It is necessary to apply a correction to this result for F + to account for the fact that the signal efficiency is in principle different for DTs involving CP -even and CPodd tags.This is because the distribution over phase space of final-state particles will be different for decays of the signal mode when it is tagged as CP -odd or CPeven.For example, the intermediate process D → K 0 S η exists in signal decays tagged by the CP -odd eigenstates, but not for CP -even tags.Comparison of DT events containing CP -even and CP -odd tags show the momentum and cos θ distributions of the signal decay are very similar, apart from that of the K 0 S momentum.Studies of these K 0 S momenta distributions, together with the known variation in reconstruction efficiency with K 0 S momentum [30], is used to determine that the ratio of the efficiency of the signal mode tagged by CP -even eigenstates to that of the efficiency of the signal mode tagged by CP -odd eigenstates is 1.008±0.06.Applying this ratio as a correction leads to the result from CP -tagged events of F + = 0.229 ± 0.013 ± 0.002.An additional systematic component for the potential difference in the efficiencies has been included, which is estimated by the difference between this central value and the one obtained from the corrected fit.This result is consistent with that obtained from CLEO-c data [8] with CP -tag modes and is a factor 1.6 more precise.

B. Measurements with the quasi-CP -tag and self-tag modes
Using Eq. ( 6), F + is determined with the quasi-CP -tag modes D → π + π − π 0 and D → π + π − π + π − , where the R T is also corrected with the factor ε(S) • ε(T )/ε(S|T ) for the same reason mentioned in Sec.VIII A. The CPeven fractions of these modes, denoted as F π + π − π 0 + and F π + π − π + π − + , are taken from Refs.[7,29].The value of R + is taken from the measurement with the CP -tag modes described in Sec.VIII A. The ratios of the DT yields to the corresponding ST yields of the D → π + π − π 0 and D → π + π − π + π − tag modes are calculated with the corresponding ST and DT yields listed in Table II.After propagating the uncertainties from the input parameters, F + is determined to be 0.227 ± 0.014 ± 0.003 with D → π + π − π 0 tags and 0.227 ± 0.016 ± 0.003 with D → π + π − π + π − tags.Here the systematic uncertainties are assigned in the same way as for the measurement with the CP -tag modes.
The self-tag yield is also sensitive to F + , as shown in Eq. (7).The ratio of the DT yield for the self-tag mode to the corresponding ST yield is determined with the corresponding yields listed in Table II and is corrected with the factor ε(S) • ε(T )/ε(S|T ).The value of R − is taken from the measurement with the CP-tag modes.The CP -even fraction measured from the self-tag mode is F + = 0.244 ± 0.019 ± 0.002.Here the systematic uncertainty is assigned using the same method as for the measurement with the CP -tag modes.
The measurement of F + with the D → K 0 S π + π − and D → K 0 L π + π − tag modes is performed with the measurements of the populations in the eight bin-pairs for the two tag modes [7].The measured DT yields for the D → K 0 S π + π − and D → K 0 L π + π − tag modes after subtracting peaking background in the eight bin pairs are shown in Fig. 9, which are determined with the same strategies described in Sec.VI.Eq. 8 and 9 are modified to account for migration effects and variations in bin-tobin efficiencies, such that the expected DT yield in the i-th bin pair as a function of F + is given by for D → K 0 L π + π − tags, where ε ij (ε ij ) is the migration matrix, determined from MC simulation, describing the efficiency for an event produced in the j-th bin and reconstructed in the i-th bin.To determine the value of F + , a log-likelihood fit is performed.The likelihood is given by where N exp i is the expected yield in the i-th bin pair as a function of F + , N i obs is the observed yield with peaking background subtracted in the i-th bin, σ N i obs is the uncertainty of the observed yield in the i-th bin, K i and σ Ki are the flavor-tagged fraction in bin i and its uncertainty, respectively and c i is the strong-phase parameter of the tag mode in bin i with covariance matrix V .The K i and c i parameters are fit parameters, but constrained through Gaussian functions.The means ( K inp i and c inp i ) and covariances (σ K inp i and V ij ) of the Gaussian functions for K i and c i are taken from the combined results from the BESIII and CLEO Collaborations [12].The elements of the migration matrix are also fit parameters, but included with χ 2 constraints, with mean ε inp ij and width σ ε inp ij , where the uncertainties arise from from the finite size of MC sample.The fit is performed twice, once with the full uncertainties included, and then with only the statistical contributions.From these fits, it is found that F + is 0.244 ± 0.021 ± 0.006, where the first uncertainty is statistical, and the second is systematic.Figure 9 shows the DT yield in each bin, together with the expected yields with the fitted value of F + , and the expected yields with other values of F + .Individual fits performed with each tag mode separately return compatible results, with 0.211 ± 0.029 for D → K 0 S π + π − tags and 0.290 ± 0.037 for D → K 0 L π + π − tags, where the combined statistical and systematic uncertainties are given.

D. Combination of results
Table III summarizes the results of F + determined with different tag modes, which are seen to be consistent with each other.A least-squared fit is performed to combine the results of F + , taking the uncertainties and correlations between the results into account.The correlations are introduced because of the common use of R + and R − in Eqs. ( 5) to (7).The correlation coefficients between the correlated tag modes are summarized in Table IV.The combined result from all tags F + is 0.235 ± 0.010 ± 0.002, which is consistent with the result 0.238 ± 0.012 ± 0.012 [8] based on the CLEO-c data, but is a factor 1.7 times more precise.The result is also compatible with the value F + = 0.226 ± 0.020 deduced from the strong-phase parameters c i , also determined with CLEO-c data [8].

IX. SUMMARY
The CP -even fraction F + of the D 0 → K 0 S π + π − π 0 decay has been measured by analyzing 2.93 fb −1 of data collected at √ s = 3.773 GeV with the BESIII detector.The measurement is performed with five categories of tag mode listed in Table III, which give a consistent set of results.The combined result is F + = 0.235 ± 0.010 ± 0.002, where the first uncertainty is statistical and the second is systematic.This result is consistent with that obtained from CLEO-c data [8], but is a factor 1.7 times more precise.The measured F + is an important input for the measurement of the unitarity triangle angle γ in B → DK, D → K 0 S π + π − π 0 decays.Currently, the measurement is dominated by statistical uncertainty.A future larger data sample [31] allows us to improve the precision significantly.

Acknowledgement
The BESIII Collaboration thanks the staff of BEPCII and the IHEP computing center for their strong sup-

8 FIG. 1 .
FIG. 1.The binning of the D → K 0 S,L π + π − Dalitz plot.The color scale represents the absolute value of the bin number |i|, where i is negative for the bin with m K 0 S,L π + < m K 0 S,L π − . .

FIG. 2 .
FIG. 2. Fits to the MBC distributions from the ST D candidates.The corresponding decay modes are denoted by the labels on the plots.The black points with error bars represent data.The blue solid curves are the fit results.The red and green dashed curves represent the signal and background contributions of the fits, respectively.

FIG. 3 .
FIG. 3. The distribution of M TBC versus M S BC for the DT candidates tagged by the D → π + π − π 0 tag mode.

FIG. 4 .
FIG. 4. The projections of the 2D fits on the M S BC distribution.The black points represent the data.Overlaid is the fit projection in the continuous red line.The blue dashed line indicates the combinatorial component.

8 FIG. 5 .
FIG. 5.The projections of the 2D fits on the M S BC distribution in the eight pairs of bins for the D → K 0 S π + π − tag mode.The black points represent the data.Overlaid is the fit projection in the continuous red line.The blue dashed line indicates the combinatorial component.
FIG. 6.Fits to the M 2 miss distributions of the D → K 0 L π 0 (top) and D → K 0 L ω (bottom) tag modes.The points with error bars represent data, the blue dashed curves are the fitted combinatorial backgrounds, the dashed red and magenta lines show the MC-simulated signal and peaking background shapes, respectively, and the blue solid curves show the total fits.

8 FIG. 7 .
FIG. 7. Fits to the M 2 miss distributions in the eight pairs of bins of D → K 0 L π + π − tag mode.The points with error bars represent data, the blue dashed curves are the fitted combinatorial backgrounds, the dashed red and magenta lines show the MC-simulated signal and peaking background shapes, respectively, and the blue solid curves show the total fits.

FIG. 8 .
FIG. 8.The R + values (top) for the CP -odd tag modes and the R − values (bottom) for the CP -even tag modes.The horizontal error bars show the total uncertainty for each measurement.The yellow bands show the fitted values with uncertainties.

FIG. 9 .
FIG. 9.Predicted and measured yields for the D → K 0 S π + π − (top) and D → K 0 L π + π − (bottom) tag modes in each pair of bins.The black points with error bars show the measured values.The red lines show the predicted values from the fit, the dashed blue lines correspond to F+ = 0, the dashed-dotted cyan lines are for F+ = 0.5 corresponding to no quantum correlation, and the dashed-dotted green lines present the expected yields with F+ = 1.

d
Also at Goethe University Frankfurt, 60323 Frankfurt am Main, Germany e Also at Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education; Shanghai Key Laboratory for Particle Physics and Cosmology; Institute of Nuclear and Particle Physics, Shanghai 200240, People's Republic of China f Also at Key Laboratory of Nuclear Physics and Ion-beam Application (MOE) and Institute of Modern Physics, Fudan University, Shanghai 200443, People's Republic of China g Also at State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, People's Republic of China h Also at School of Physics and Electronics, Hunan University, Changsha 410082, China i Also at Guangdong Provincial Key Laboratory of Nuclear Science, Institute of Quantum Matter, South China Normal University, Guangzhou 510006, China j Also at Frontiers Science Center for Rare Isotopes, Lanzhou University, Lanzhou 730000, People's Republic of China k Also at Lanzhou Center for Theoretical Physics, Lanzhou University, Lanzhou 730000, People's Republic of China l Also at the Department of Mathematical Sciences, IBA, Karachi 75270, Pakistan

TABLE I .
The tag modes used in this analysis.

TABLE II .
The DT and ST yields (NDT and NST) for each tag mode.For the fully reconstructed tag modes, the first uncertainties are statistical and the second are systematic.

TABLE IV .
Correlation coefficients between the results measured by different types of tag modes.