Measurements of the absolute branching fractions of $D^{0(+)}\to K\bar K\pi\pi$ decays

Based on 2.93~fb$^{-1}$ $e^+e^-$ collision data taken at center-of-mass energy of 3.773 GeV by the BESIII detector, we report the measurements of the absolute branching fractions of $D^0\to K^+K^-\pi^0\pi^0$, $D^0\to K^0_SK^0_S\pi^+\pi^-$, $D^0\to K^0_SK^-\pi^+\pi^0$, $D^0\to K^0_SK^+\pi^-\pi^0$, $D^+\to K^+K^-\pi^+\pi^0$, $D^+\to K^0_SK^+\pi^0\pi^0$, $D^+\to K^0_SK^-\pi^+\pi^+$, $D^+\to K^0_SK^+\pi^+\pi^-$, and $D^+\to K^0_SK^0_S\pi^+\pi^0$. The decays $D^0\to K^+K^-\pi^0\pi^0$, $D^0\to K^0_SK^-\pi^+\pi^0$, $D^0\to K^0_SK^+\pi^-\pi^0$, $D^+\to K^0_SK^0_S\pi^+\pi^0$, and $D^+\to K^0_SK^+\pi^0\pi^0$ are observed for the first time. The branching fractions of the decays $D^0\to K^0_SK^0_S\pi^+\pi^-$, $D^+\to K^+K^-\pi^+\pi^0$, $D^+\to K^0_SK^-\pi^+\pi^+$, and $D^+\to K^0_SK^+\pi^+\pi^-$ are measured with improved precision compared to the world-average values.

Current measurements of the D 0(+) → KKππ decays containing K 0 S or π 0 are limited [9]. The branching fractions (BFs) of D 0 → K 0 S K 0 S π + π − [10,11], D + → K 0 S K − π + π + [12], D + → K 0 S K + π + π − [12], and D + → K + K − π + π 0 [13] were only determined relative to some well known decays or via topological normalization, with poor precision. This paper presents the first direct measurements of the absolute BFs for the decays S π 0 π 0 decay is not included since it suffers from poor statistics and high background. Throughout this paper, charge conjugate processes are implied. An e + e − collision data sample corresponding to an integrated luminosity of 2.93 fb −1 [14] collected at a center-of-mass energy of √ s = 3.773 GeV with the BESIII detector is used to perform this analysis.

II. BESIII DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [15] located at the Beijing Electron Positron Collider (BEPCII) [16]. The cylindrical core of the BESIII detector consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-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. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identifier modules interleaved with steel.
The acceptance of charged particles and photons is 93% over 4π solid angle. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for the 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 of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps.
Simulated samples produced with the geant4based [17] Monte Carlo (MC) package including the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the backgrounds.
The simulation includes the beam-energy spread and initialstate radiation (ISR) in the e + e − annihilations modeled with the generator kkmc [18]. The inclusive MC samples consist of the production of DD pairs with consideration of quantum coherence for all neutral D modes, the non-DD decays of the ψ(3770), the ISR production of the J/ψ and ψ(3686) states, and the continuum processes. The known decay modes are modeled with evtgen [19] using the BFs taken from the Particle Data Group (PDG) [9], and the remaining unknown decays from the charmonium states are modeled with lundcharm [20]. The finalstate radiations from charged final-state particles are incorporated with the photos package [21].

III. MEASUREMENT METHOD
The D 0D0 or D + D − pair is produced without an additional hadron in e + e − annihilations at √ s = 3.773 GeV. This process offers a clean environment to measure the BFs of the hadronic D decay with the double-tag (DT) method. The single-tag (ST) candidate events are selected by reconstructing aD 0 or D − in the following hadronic final states:D 0 → K + π − , K + π − π 0 , and K + π − π − π + , and The event in which a signal candidate is selected in the presence of an STD meson, is called a DT event. The BF of the signal decay is determined by where N tot ST = i N i ST and N net DT are the total yields of the ST and DT candidates in data, respectively. N i ST is the ST yield for the tag mode i. For the signal decays involving K 0 S meson(s) in the final states, N net DT is the net DT yields after removing the peaking background from the corresponding non-K 0 S decays. For the other signal decays, the variable corresponds to the fitted DT yields as described later. Here, ǫ sig is the efficiency of detecting the signal D decay, averaged over the tag mode i, which is given by: where ǫ i ST and ǫ i DT are the efficiencies of detecting ST and DT candidates in the tag mode i, respectively.

IV. EVENT SELECTION
The selection criteria of K ± , π ± , K 0 S , and π 0 are the same as those used in the analyses presented in Refs. [22][23][24][25][26][27][28][29]. All charged tracks, except those from K 0 S decays, are required to have a polar angle θ with respect to the beam direction within the MDC acceptance |cosθ| < 0.93, and a distance of closest approach to the interaction point (IP) within 10 cm along the beam direction and within 1 cm in the plane transverse to the beam direction. Particle identification (PID) for charged pions, kaons, and protons is performed by exploiting TOF information and the specific ionization energy loss dE/dx measured by the MDC. The confidence levels for pion and kaon hypotheses (CL π and CL K ) are calculated. Kaon and pion candidates are required to satisfy CL K > CL π and CL π > CL K , respectively.
The K 0 S candidates are reconstructed from two oppositely charged tracks to which no PID criteria are applied and which masses are assumed to be that of pions. The charged tracks from the K 0 S candidate must satisfy |cosθ| < 0.93. In addition, due to the long lifetime of the K 0 S meson, there is a less stringent criterion on the distance of closest approach to the IP in the beam direction of less than 20 cm and no requirement on the distance of closest approach in the plane transverse to the beam direction. Furthermore, the π + π − pairs are constrained to originate from a common vertex and their invariant mass is required to be within (0.486, 0.510) GeV/c 2 , which corresponds to about three times the fitted resolution around the nominal K 0 S mass. The decay length of the K 0 S candidate is required to be larger than two standard deviations of the vertex resolution away from the IP.
The π 0 candidate is reconstructed via its γγ decay. The photon candidates are selected using the information from the EMC shower. It is required that each EMC shower starts within 700 ns of the event start time and its energy is greater than 25 (50) MeV in the barrel (end cap) region of the EMC [15]. The energy deposited in the nearby TOF counters is included to improve the reconstruction efficiency and energy resolution. The opening angle between the candidate shower and the nearest charged track must be greater than 10 • . The γγ pair is taken as a π 0 candidate if its invariant mass is within (0.115, 0.150) GeV/c 2 . To improve the resolution, a kinematic fit constraining the γγ invariant mass to the π 0 nominal mass [9] is imposed on the selected photon pair.

V. YIELDS OF STD MESONS
To selectD 0 → K + π − candidates, the backgrounds from cosmic rays and Bhabha events are rejected by using the same requirements described in Ref. [30]. In the selection ofD 0 → K + π − π − π + candidates, theD 0 → K 0 S K ± π ∓ decays are suppressed by requiring the mass of all π + π − pairs to be outside (0.478, 0.518) GeV/c 2 .
The taggedD mesons are identified using two variables, namely the energy difference and the beam-constrained mass Here, E b is the beam energy, p tag and E tag are the momentum and energy of theD candidate in the rest frame of e + e − system, respectively. For each tag mode, if there are multiple candidates in an event, only the one with the smallest |∆E tag | is kept. The taggedD candidates are required to satisfy ∆E tag ∈ (−55, 40) MeV for the tag modes containing π 0 in the final states and ∆E tag ∈ (−25, 25) MeV for the other tag modes, thereby taking into account the different resolutions.
To extract the yields of STD mesons for individual tag modes, binned-maximum likelihood fits are performed on the corresponding M tag BC distributions of the accepted ST candidates following Refs. [22][23][24][25][26][27][28]. In the fits, theD signal is modeled by an MC-simulated shape convolved with a double-Gaussian function describing the resolution difference between data and MC simulation. The combinatorial background shape is described by an ARGUS function [31] defined as BC , E end is an endpoint fixed at 1.8865 GeV, A f is a normalization factor, and ξ f is a free parameter. The resulting fits to the M BC distributions for each mode are shown in Fig. 1. The total yields of the STD 0 and D − mesons in data are 2327839 ± 1860 and 1558159 ± 2113, respectively, where the uncertainties are statistical only.

VI. YIELDS OF DT EVENTS
In the recoiling sides against the taggedD candidates, the signal D decays are selected by using the residual tracks that have not been used to reconstruct the tagged D candidates. To suppress the K 0 S contribution in the individual mass spectra for the D 0 → K + K − π 0 π 0 , D 0 → K 0 S K 0 S π + π − , and D + → K 0 S K + π + π − decays, the π + π − and π 0 π 0 invariant masses are required to be outside (0.468, 0.528) GeV/c 2 and (0.438, 0.538) GeV/c 2 , respectively. To suppress the background from These requirements correspond to at least five times the fitted mass resolution away from the fitted mean of the mass peak.
The signal D mesons are identified using the energy difference ∆E sig and the beam-constrained mass M sig BC , which are calculated with Eqs. (3) and (4) by substituting "tag" with "sig". For each signal mode, if there are multiple candidates in an event, only the one with the smallest |∆E sig | is kept. The signal decays are required to satisfy the mode-dependent ∆E sig requirements, as shown in the second column of Table 1. To suppress incorrectly identified DD candidates, the opening angle between the taggedD and the signal D is required to be greater than 160 • , resulting in a loss of (2-6)% of the signal and suppressing (8-55)% of the background. The events smeared along the diagonal, named BKGII, are mainly from the e + e − → qq processes. The events with uncorrelated and incorrectly reconstructed D andD, named BKGIII, disperse across the whole allowed kinematic region.
For each signal D decay mode, the yield of DT events (N fit DT ) is obtained from a two-dimensional (2D) unbinned maximum-likelihood fit [32] on the M tag BC versus M sig BC distribution of the accepted candidates. In the fit, the probability density functions (PDFs) of signal, BKGI, BKGII, and BKGIII are constructed as • signal: a(x, y), • BKGII: c z (z; √ 2E b , ξ z ) · g(k), and respectively. Here, x = M sig BC , y = M tag BC , z = (x + y)/ √ 2, and k = (x − y)/ √ 2. The PDFs of signal a(x, y), b(x), and b(y) are described by the corresponding MC-simulated shapes. c f (f ; E end , ξ f ) is an ARGUS function [31] defined above, where f denotes x, y, or z; E b is fixed at 1.8865 GeV. g(k) is a Gaussian function with mean of zero and standard deviation parametrized by Combinatorial π + π − pairs from the decays S π + π + π − π 0 [and D + → 2(π + π − )π + π 0 ] may also satisfy the K 0 S selection criteria and form peaking backgrounds around M D in the M sig BC distributions for D 0 → K 0 respectively. This kind of peaking background is estimated by selecting events in the K 0 S sideband region of (0.454, 0.478) ∪ (0.518, 0.542) GeV/c 2 . For D 0 → K 0 S K − π + π 0 , D 0 → K 0 S K + π − π 0 , D + → K 0 S K − π + π + , D + → K 0 S K + π + π − , and D + → K 0 S K + π 0 π 0 decays, one-dimensional (1D) signal and sideband regions are used. For D 0 → K 0 S K 0 S π + π − and D + → K 0 S K 0 S π + π 0 decays, 2D signal and sideband regions are used. The 2D K 0 S signal region is defined as the square region with both π + π − combinations lying in the K 0 S signal regions. The 2D K 0 S sideband 1 (2) regions are defined as the square regions with 1 (2) π + π − combination(s) located in the 1D K 0 S sideband regions and the rest in the 1D K 0 S signal region. Figure 3 shows 1D and 2D π + π − invariant-mass distributions as well as the K 0 S signal and sideband regions. (a) The π + π − invariant-mass distributions of the D + → K 0 S K − π + π + candidate events of data (points with error bars) and inclusive MC sample (histogram). Pairs of the red solid (blue dashed) arrows denote the K 0 S signal (sideband) regions. (b) Distribution of M π + π − (1) versus M π + π − (2) for the D 0 → K 0 S K 0 S π + π − candidate events in data. Red solid box denotes the 2D signal region. Pink dotdashed (blue dashed) boxes indicate the 2D sideband 1 (2) regions.
For the signal decays involving K 0 S meson(s) in the final states, the net yields of DT events are calculated by subtracting the sideband contribution from the DT fitted yield by Here, N = 1 for the decays with one K 0 S meson while N = 2 for the decays with two K 0 S mesons. The combinatorial π + π − backgrounds are assumed to be uniformly distributed and double-counting is avoided by subtracting (2) yields from (1) yields appropriately. N fit DT and N fit sidi are the fitted D yields in the 1D or 2D signal region and sideband i region, respectively. For the other signal decays, the net yields of DT events are N fit DT . Figure 4 shows the M tag BC and M sig BC projections of the 2D fits to data. From these 2D fits, we obtain the DT yields for individual signal decays as shown in Table 1.
For each signal decay mode, the statistical significance is calculated according to −2ln(L 0 /L max ), where L max and L 0 are the maximum likelihoods of the fits with and without involving the signal component, respectively. The effect of combinatorial π + π − backgrounds in the K 0 S -signal regions has been considered for the decays involving a K 0 S . The statistical significance for each signal decay is found to be greater than 8σ.

VII. RESULTS
Each of the D 0 → K 0 S K − π + π 0 , D + → K + K − π + π 0 , D + → K 0 S K − π + π + , and D + → K 0 S K + π + π − decays is modeled by the corresponding mixed signal MC samples, in which the dominant decay modes containing resonances of K * (892), ρ(770), and φ are mixed with the phase space (PHSP) signal MC samples. The mixing ratios are determined by examining the corresponding invariant mass and momentum spectra. The other decays, which are limited in statistics, are generated with the PHSP generator. The momentum and the polar angle distributions of the daughter particles and the invariant masses of each two-and three-body particle combinations of the data agree with those of the MC simulations. As an example, Fig. 5 shows the invariant mass distributions of two or three-body particle combinations of D + → K + K − π + π 0 candidate events for data and MC simulations.
The measured values of N net DT , ǫ sig , and the obtained BFs are summarized in Table 1. The current worldaverage values are also given for comparison. The signal efficiencies have been corrected by the necessary data-MC differences in the selection efficiencies of K ± and π ± tracking and PID procedures and the π 0 reconstruction. These efficiencies also include the BFs of the K 0 S and π 0 decays. The efficiency for D + → K 0 S K + π − π 0 ) due to the K 0 S (ω) rejection in the π + π − (K 0 S π 0 ) mass spectrum.

VIII. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties are estimated relative to the measured BFs and are discussed below. In BF determinations using Eq. (1), all uncertainties associated with the selection of taggedD canceled in the ratio. The systematic uncertainties in the total yields of STD mesons related to the M BC fits to the STD candidates, were previously estimated to be 0.5% for both neutral and chargedD [22][23][24].
The systematic uncertainty in the 2D fit to the M tag BC versus M sig BC distribution is examined via the repeated measurements in which the signal shape and the endpoint of the ARGUS function (±0.2 MeV/c 2 ) are varied. Quadratically summing the changes of the BFs for these two sources yields the corresponding systematic uncertainties.
The systematic uncertainty due to the ∆E sig requirement is assigned to be 0.3%, which corresponds to the largest efficiency difference with and without smearing the data-MC Gaussian resolution of ∆E sig for signal MC events. Here, the smeared Gaussian parameters are obtained by using the samples of DT events D 0 → K 0 S π 0 , D 0 → K − π + π 0 , D 0 → K − π + π 0 π 0 , and D + → K − π + π + π 0 versus the sameD tags in our nominal analysis. The systematic uncertainties due to K 0 S sideband choice and K 0 S rejection mass window are assigned by examining the changes of the BFs via varying nominal K 0 S sideband and corresponding rejection window by ±5 MeV/c 2 .
For the decays whose efficiencies are estimated with mixed signal MC events, the systematic uncertainty in the MC modeling is determined by comparing the signal efficiency when changing the percentage of MC sample components. For the decays whose efficiencies are estimated with PHSP-distributed signal MC events, the uncertainties are assigned as the change of the signal efficiency after adding the possible decays containing K * (892) or ρ(770). The imperfect simulations of the momentum and cos θ distributions of charged particles are considered as a source of systematic uncertainty. The signal efficiencies are re-weighted by those distributions in data with background subtracted. The largest change of the re-weighted to nominal efficiencies, 0.9%, is assigned as the corresponding systematic uncertainty.
The measurements of the BFs of the neutral D decays are affected by quantum correlation effect. For each neutral D decay, the CP -even component is estimated by the CP -even tag D 0 → K + K − and the CP -odd tag D 0 → K 0 S π 0 . Using the same method as described in Ref. [34] and the necessary parameters quoted from Refs. [35][36][37] S K − π + π 0 , and D 0 → K 0 S K + π − π 0 , respectively. After correcting the signal efficiencies by the individual factors, the residual uncertainties are assigned as systematic uncertainties.
The uncertainties due to the limited MC statistics for various signal decays, (0.4-0.8)%, are taken into account as a systematic uncertainty. The uncertainties of the quoted BFs of the K 0 S → π + π − and π 0 → γγ decays are 0.07% and 0.03%, respectively [9].
The efficiencies of DD opening angle requirement is studied by using the DT events of D 0 → K − π + π + π − , D 0 → K − π + π 0 π 0 , and D + → K − π + π + π 0 tagged by the same tag modes in our nominal analysis. The difference of the accepted efficiencies between data and MC simulations, 0.4% for the decays without π 0 , 0.8% for the decays involving one π 0 and 0.3% for the decays involving two π 0 s, is assigned as the associated systematic uncertainty. Table 2 summarizes the systematic uncertainties in the BF measurements. For each signal decay, the total systematic uncertainty is obtained by adding the above effects in quadrature to be (2.6-6.0)% for various signal decay modes.

IX. SUMMARY
In summary, by analyzing a data sample obtained in e + e − collisions at √ s = 3.773 GeV with the BESIII detector and corresponding to an integrated luminosity of 2.93 fb −1 , we obtained the first direct measurements of the absolute BFs of nine D 0(+) → KKππ decays containing K 0 S or π 0 mesons. The D 0 → K + K − π 0 π 0 , S K + π 0 π 0 , and D + → K 0 S K 0 S π + π 0 decays are observed for the first time. Compared to the world-average values, the BFs of the D 0 → K 0 S K 0 S π + π − , D + → K + K − π + π 0 , D + → K 0 S K − π + π + , and D + → K 0 S K + π + π − decays are measured with improved precision. Our BFs of D + → K 0 S K − π + π + and D + → K 0 S K + π + π − are in agreement with individual world averages within 1σ while our BFs of D 0 → K 0 S K 0 S π + π − and D + → K + K − π + π 0 deviate with individual world averages by 2.3σ and 2.8σ, respectively. The precision of the BF of D + → K + K − π + π 0 is improved by a factor of about seven. Future amplitude analyses of all these D 0(+) → KKππ decays with larger data samples foreseen at BESIII [39], Belle II [40], and LHCb [41] will supply rich information of the two-body decay modes containing scalar, vector, axial and tensor mesons, thereby benefiting the understanding of quark SU(3)-flavor symmetry.