Amplitude analysis and branching fraction measurement of $D_s^+ \to K^-K^+\pi^+\pi^0$

The first amplitude analysis of the decay $D_s^+\to K^-K^+\pi^+\pi^0$ is presented using the data samples, corresponding to an integrated luminosity of 6.32 fb$^{-1}$, collected with the BESIII detector at $e^+e^-$ center-of-mass energies between 4.178 and 4.226 GeV. More than 3000 events selected with a purity of 97.5\% are used to perform the amplitude analysis, and nine components are found necessary to describe the data. Relative fractions and phases of the intermediate decays are determined. With the detection efficiency estimated by the results of the amplitude analysis, the branching fraction of $D_s^+\to K^-K^+\pi^+\pi^0$ decay is measured to be $(5.42\pm0.10_{\rm stat.}\pm0.17_{\rm syst.})\%$.

The first amplitude analysis of the decay D + s → K − K + π + π 0 is presented using the data samples, corresponding to an integrated luminosity of 6.32 fb −1 , collected with the BESIII detector at e + e − center-of-mass energies between 4.178 and 4.226 GeV. More than 3000 events selected with a purity of 97.5% are used to perform the amplitude analysis, and nine components are found necessary to describe the data. Relative fractions and phases of the intermediate decays are determined. With the detection efficiency estimated by the results of the amplitude analysis, the branching fraction of D + s → K − K + π + π 0 decay is measured to be (5.42 ± 0.10stat. ± 0.17syst.)%.

I. INTRODUCTION
Accurate measurement of D + s decays are important for the studies of other decay processes that are dominated by final states involving D + s mesons, particularly for those of B 0 s decays [1]. The decay D + s → K − K + π + π 0 is a Cabibbo-favored decay (the inclusion of charge conjugated reactions is implied throughout this paper). Due to its large branching fraction (BF), it is usually selected as a "tag mode" for the measurement of other decays of the D + s meson [2][3][4][5][6][7]. However, the BF of the D + s → K − K + π + π 0 decay has a large systematic uncertainty due to the poor knowledge of intermediate state processes [8,9]. An amplitude analysis of this decay is expected to provide a detailed understanding of its intermediate structures and significantly improve the experimental precision of its decay BF.
The four-body hadronic decays of D + s mesons are dominated by two-body intermediate processes, for example D + s → AP and D + s → V V decays, where V , A, and P denote vector, axial-vector, and pseudoscalar mesons, respectively. Measurements of the BFs of these two-body decays are important to test theoretical calculations [10][11][12][13] and to better understand the decay mechanisms of the D + s meson. In recent years, many measurements of D + s → P P and D + s → V P decays have been reported [14]. However, there are few studies focusing on D + s → AP and D + s → V V decays. The amplitude analysis of D + s → AP decay allows the study of substructures involving K 1 (1270), K 1 (1400), and f 1 (1420) mesons. The measurements of the intermediate resonances K 1 (1270) and K 1 (1400) are also useful for understanding the mixing of these two axial-vector kaons [15]. For D + s → V V , two processes, namely D + s → φρ + and D + s →K * 0 K * + , which are represented by the decay diagrams in Fig. 1, can be studied in the D + s → K − K + π + π 0 decay. The BF of the decay D + s → φρ + was measured to be (8.4 +1. 9 −2.3 )% [16] by the CLEO experiment based on a data sample corresponding to an integrated luminosity of 689 pb −1 at the Υ(3S) and Υ(4S) resonances and at e + e − center-of-mass energies (E cm ) just below and above the Υ(4S) resonance. The previous most precise determination of the BF of D + s →K * 0 K * + decay, (7.2 ± 2.6)% [17], was performed by the ARGUS experiment using a data sample of 432 pb −1 collected at E cm = 10.4 GeV. The goal of the present analysis is to improve the precision of these measurements. Moreover, a recent study [18] points out that the measured values of the ratio of K 1 (1270) decay (R K1(1270) = B(K1(1270)→K * π) B(K1(1270)→Kρ) ), which are listed in Table I, are inconsistent between different experiments [19][20][21][22][23][24][25]. They are expected to be identical under the narrow width approximation for the K 1 (1270) meson and assuming CP conservation in strong decays [18]. The decays related to this ratio may be observed in the D + s → K − K + π + π 0 decay. Table I. Values of R K 1 (1270) determined by different experiments. Fit 1 and fit 2 refer to amplitude analyses performed with the mass and width of the K1(1270) + meson fixed or left free in the fit, respectively.

II. DETECTOR AND DATA SETS
The BESIII detector is a magnetic spectrometer [27,28] located at the Beijing Electron Position Collider (BEPCII) [29]. 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 the 4π solid angle. The resolution of charged-particle momentum at 1 GeV/c is 0.5% while that of the specific ionization energy loss (dE/dx) 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 of the TOF barrel part is 68 ps. The end-cap TOF system was upgraded in 2015 with multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [30,31].
The data samples used in this analysis contain a total integrated luminosity of 6.32 fb −1 collected at center-ofmass energies between E cm = 4.178 and 4.226 GeV with the BESIII detector. The integrated luminosity of each data sample is shown in Table II. In this energy region, pairs of D ± s D * ∓ s mesons are produced. The D * ± s meson predominantly decays to D ± s γ (93.5%), and only a small fraction decays to D s π 0 (5.8%) [14]. A doubletag (DT) technique is employed to measure the absolute BF of the D + s decays [32]. First, the D − s meson is fully reconstructed in one of the following decay modes: These are referred to as single-tag (ST) events. Second, the D + s → K − K + π + π 0 decay events are selected. Generic Monte Carlo (MC) simulated event samples are produced with the GEANT4-based [33,34] software at E cm = 4.178 − 4.226 GeV. The samples include all known open charm decays; the continuum process (e + e − → qq, q = u, d, and s); Bhabha scattering; the µ + µ − , τ + τ − , and diphoton processes; and the cc resonances (J/ψ, ψ(3686), and ψ(3770)) via initial-state radiation (ISR). The open charm processes are generated using conexc [35], and their subsequent decays are modeled by evtgen [36] with known BFs from the Particle Data Group (PDG) [14]. The simulation of ISR production of J/ψ, ψ(3686), and ψ(3770) mesons is performed with kkmc [37]. The effects of final-state radiation (FSR) from charged tracks are simulated by photos [38]. The remaining unknown decays are generated with the lundcharm model [39]. The generic MC sample is used to study backgrounds and determine the efficiencies of tag modes and the signal mode for the BF measurement, in which our amplitude analysis model is used to generate the signal mode events.
A phase-space (PHSP) MC sample is produced with the D + s meson decaying to K − K + π + π 0 generated with a uniform distribution and the D − s meson decaying to the tag modes. Initially, the PHSP MC sample is used to calculate the normalization integral used in the determination of the amplitude model parameters in the fit to data. Then, the signal MC sample is re-generated with the D + s meson decaying to K − K + π + π 0 using the amplitude model and the D − s meson decaying to the tag modes. The normalization integral performed with signal MC samples results in more accurate fit parameters of magnitudes and phases and improves the computational efficiency of the MC integration. The signal MC sample is also used to calculate the goodness of the fit in this analysis. The PHSP MC sample is used to determine the efficiency mentioned in Section IV A.

III. EVENT SELECTION
Charged tracks except for those from K 0 S decays are required to have a distance of closest approach to the interaction point (IP) within 1 cm in the transverse plane and within 10 cm along the MDC axis (z axis). The polar angle of the charged track with respect to the z axis θ is required to satisfy | cos θ| < 0.93. Kaons and pions are identified by combining the dE/dx information in the MDC and the time-of-flight from the TOF. If the probability of the kaon hypothesis is larger than that of the pion hypothesis, the track is identified as a kaon. Otherwise, the track is identified as a pion. Particle identification (PID) is not performed for the π + or π − from K 0 S decays.
The π 0 and η candidates are reconstructed via diphoton decays. The timing of the electromagnetic showers in the EMC is required to be within [0,700] ns of the trigger time, and the deposited energy must be greater than 25 (50) MeV in the barrel (endcap) region of the EMC. Good showers must satisfy | cos θ| < 0.80 (0.86 < | cos θ < 0.92) in the barrel (endcap) and be more than 20 • from the nearest charged track. The unconstrained invariant masses of π 0 , η and η ′ (η ′ → π + π − η γγ ) are required to be within [115,150] MeV/c 2 , [500, 570] MeV/c 2 , and [946, 970] MeV/c 2 , respectively. A kinematic fit is performed to constrain M γγ to the known π 0 (η) mass, and the χ 2 of the corresponding fit is required to be less than 30 (20) for π 0 (η) candidates.
The K 0 S candidates are reconstructed in the decay K 0 S → π + π − . Two oppositely charged tracks with distances of closest approach to the IP less than 20 cm along the z axis are assigned as π + π − without further PID requirements. A constrained vertex fit of each pair of tracks is performed. Candidate K 0 S particles are required to have the χ 2 of the vertex fit less than 100 and an invariant mass of the π + π − pair (M π + π − ) in the range [487, 511] MeV/c 2 . In the case of the decay modes D − the decay length of the K 0 S candidates obtained with the secondary vertex fit [40] must be at least two times its fit uncertainty. For the D − s → K − π + π − process, M π + π − is required to be outside of the range [487,511] MeV/c 2 , to remove possible misidentified events of D − s → K 0 S K − . To identify the process e + e − → D * − s D + s , the recoil mass M rec of D − s candidates is defined as [14] and p D − s is the momentum of the D − s candidate. The values of M rec are required to be in the regions depending on the center-ofmass energy as listed in Table II. The D − s mass windows for the eight tag modes are shown in Table III.
The D + s meson decays with invariant masses M D + s in the region [1.87, 2.06] GeV/c 2 are selected. Good vertex fits of all charged tracks on both the signal and the tag side are required. A multi-constraint kinematic fit of e + e − → D * ± s D ∓ s → γD ± s D ∓ s with D − s decaying to one of the tag modes and D + s decaying to the signal mode is performed. The set of constraints including four-momentum conservation in the e + e − system and the mass constraints of the π 0 meson, the D + s meson, the D − s meson and the D * ± s meson is labeled C 1 . Based on the requirements of C 1 , a set of constraints C 2 is defined by excluding the signal M D + s constraint, and C 3 is defined by excluding the mass constraints of the D ± s meson on both the signal and tag sides.
If there are multiple candidate combinations in an event, the candidate with the minimum χ 2 of the C 2 kinematic fit (χ 2 C2 ) is chosen. A good C 1 kinematic fit is required. To reduce the background while avoiding peaking background which is caused by constraining the mass of the D ± s meson (M D ± s ), the χ 2 of the C 3 kinematic fit (χ 2 C3 ) is required to be less than 250. The classes of background events, which are listed in Table IV, are rejected. For backgrounds categorized as (a), (b) and (c), a π 0 from the D − s decay is wrongly associated to the D + s meson on the opposite side. These are vetoed if the χ 2 of the C 1 kinematic fit (χ 2 C1 ) of the reconfigured combination is better than that of the original. For backgrounds categorized as (d), the events with D + → K − π + π + decay versus D − → K + K − π − π 0 decay are wrongly reconstructed as D − s → K − π + π − decay versus D + s → K + K − π + π 0 decay, when a π − meson from D − decay is exchanged with a π + meson from D + decay. If the reconstructed D ± masses of the signal and the tag modes fall in the region within 0.055 GeV/c 2 of the nominal M D ± , the events are rejected. For background categories (e) and (f), events with K 0 is the nominal mass of D 0 [14]. For background category (g), the wrong signal combination survives due to exchanging the π 0 meson from D 0 decay and the π − meson fromD 0 decay and misidentifying the π − meson as a K − meson. They are suppressed by rejecting events satisfying |M K − π + π 0 − M PDG D 0 | < 0.055 GeV/c 2 and |M K + π − π + π − − M PDG D 0 | < 0.055 GeV/c 2 . Background type (h) events are suppressed by applying a veto on events with |M K − π + π 0 − M PDG D 0 | < 0.045 GeV/c 2 .
Events containing a possible mis-formed π 0 meson on the signal side are also rejected. Events in which the invariant mass of the higher-energy photon from the signal side combined with a photon from the D * s → D s γ decay is within [0.12, 0.15] GeV/c 2 and with |dM recombined | < |dM | are rejected, where dM is the mass difference between the signal D + s meson and the tagged D − s meson, and dM recombined is the corresponding mass difference with the signal π 0 replaced by the recombined π 0 . A veto is also applied to reject events with recombined mass of the higher-energy photon from the sig-nal side and the photon from the tag side falling within [0.12, 0.15] GeV/c 2 .
After the full selection, the invariant mass spectra of the signal D + s candidates for data samples collected at center-of-mass energies 4.178-4.226 GeV are shown in Fig. 2, together with fits to the mass spectra. There are 1708, 1024, and 356 events retained in the signal region [1.935, 1.99] GeV/c 2 for the amplitude analysis with purities, w sig , of 97.7 ± 0.4%, 97.3 ± 0.5% and 97.5 ± 0.8% at E cm = 4.178 GeV, 4.189-4.219 GeV, and 4.226 GeV, respectively. Studies of the generic MC samples show that peaking background is negligible. The background description by the generic MC has been verified by comparisons of data with the generic MC samples in the sideband regions [1.88, 1.92] GeV/c 2 and [2.00, 2.04] GeV/c 2 . A good agreement is found, and the generic MC samples are used to model the residual background contamination in the signal region. The four-momenta of the final state particles after a two-constraint kinematic fit to the signal candidate, constraining the D + s mass and π 0 mass to their known values [14], are used to perform the amplitude analysis.

IV. AMPLITUDE ANALYSIS
The amplitude analysis of D + s → K − K + π + π 0 decay is performed by using an unbinned maximum likelihood fit. The likelihood function is constructed with the probability density function (PDF) described in the following, in which the momenta of the four daughter particles are used as inputs.

A. Likelihood Function Construction
The PDF used to construct the likelihood of the amplitude is given by where p j is the set of the final state particles' four momenta, and ǫ(p j ) is the detection efficiency parameterized in terms of the final state particles' four momenta. The PDF f S (p j ) is normalized by the integration. The standard element of the four-body PHSP [26] is defined as where j runs over the four daughter particles and E j is the energy of particle j. This analysis uses an isobar model formulation, where the signal decay amplitude, M (p j ), is represented as a coherent sum of many two-body amplitude modes where c n is written in the polar form as ρ n e iφn (ρ n and φ n are the magnitude and phase for the n th amplitude, respectively). A n (p j ) is the n th amplitude function modeled as where the indices 1 and 2 correspond to the two intermediate resonances, respectively. X D + s n (p j ) is the Blatt-Weisskopf barrier factor [26,[41][42][43] for the D + s meson, while P 1,2 n (m 1 , m 2 ) and X 1,2 n (p j ) are the propagators and Blatt-Weisskopf barrier factors of the intermediate resonances 1 and 2, respectively. For non-resonant states, the propagator is set to unity, which can be regarded as a very broad resonance. S n (p j ) is the spin factor which is constructed with the covariant tensor formalism [26].
The 2.5% background contribution is described by the background PDF: The background events in the signal region from the generic MC sample are used to model the corresponding background in data. The background shape B(p j ) is derived using a multi-dimensional kernel density estimator [44] with five independent Lorentz invariant variables [45], which models the distribution of an input dataset as a superposition of Gaussian kernels. This background PDF is then added to the signal PDF incoherently, and the combined PDF is written as where the factor ǫ(p j ) in the numerator can be taken out as in Eq. 8. In this way, the ǫ(p j ) term, which is independent of the fitted variables, is a constant and can be dropped in the likelihood fit. For the determination of ǫ(p j ), totally 300 million PHSP MC events at E cm = 4.178 GeV, 4.189-4.219 GeV and 4.226 GeV are generated, and near 15 million events are selected with the event selection. The background shape is determined from the selected generic MC events; hence one has to divide the background function by the efficiency, The combined PDF becomes The corresponding likelihood function is defined as where i denotes the data sample, k i runs over each event and N i data is the number of events in data sample i. The log-likelihood is used to perform the max-likelihood calculation.
The PDFs and the efficiencies are considered separately for three data samples with E cm = 4.178 GeV, 4.189-4.219 GeV, and 4.226 GeV, corresponding to how the data were collected. Therefore, the log-likelihood functions for the three samples are summed up The normalization integrals in the denominator of Eq. 8 are obtained by summing over an MC sample where N MC is the number of the selected MC events and M gen (p j ) is the amplitude that is set with the parameters used to generate the signal MC sample, which are initially obtained from the results using the PHSP MC integration. M gen (p j ) is a constant over the whole PHSP. Then with the results obtained from the fit to data, the signal MC sample is generated and used in MC integration. Totally 12 million PHSP MC events and 10 million signal MC events are selected at E cm = 4.178 GeV, 4.189-4.219 GeV and 4.226 GeV with satisfying all selection criteria as those of the data sample.
The effect from the PID, tracking and reconstruction efficiency differences between data and simulation is considered by multiplying the weight of the MC event by a factor γ ǫ , which is calculated as where j = K ∓ , K ± , π ± and π 0 . The signal MC integration becomes

Spin Factors
For a decay process of the form a → bc, p a , p b , p c are used to denote the momenta of the particles a, b, c, respectively. The spin projection operator P (S) µ1...µSν1...νS (a), for a resonance a with spin S = 0, 1, 2 and four-momentum p a is given by where g µµ ′ is the Minkowski metric.
The covariant tensorst µ1...µ l (a) [26] for the final states of pure orbital angular momentum L are constructed from the relevant momenta p a , p b , p c : where r a = p b − p c . The corresponding covariant tensors with L = 0, 1, 2 are given as The eleven types of decay modes used in this analysis are listed in Table V.

Blatt-Weisskopf Barrier Factors
The Blatt-Weisskopf barrier X(p j ) [26,[41][42][43] is a barrier function for a two-body decay process, a → bc. The Blatt-Weisskopf barrier depends on the angular momenta and the momenta of the final state particles in the rest frame of the parent particle. The definition is given by where L denotes the orbital angular momentum; R is the effective radius of the barrier; the values of R used in this analysis are taken to be 3.0 GeV −1 and 5.0 GeV −1 for intermediate resonances and the D + s meson, respectively [46]; and q is the magnitude of the momenta of the final-state particles in the rest system of the parent particle.
For a process a → bc, s i = E 2 i − p 2 i is defined, where i denotes a, b, c and E i , p i are the particle's energy and momentum, such that
The RBW function is given by Decay chain where m = E 2 − p 2 and m 0 is the nominal mass of the resonance, and Γ(m) is given by where q 0 indicates the value of q when s a = m 2 0 . Considering the obvious mass deviation, the mass and width ofK 0 1 (1270) are set to the average values (shown in Table VI) without including the results from Belle [47].
The function f (m) is given by where and m π is the charged pion mass. The normalization condition at P GS (0) fixes the parameter d = f (0)/(Γ 0 m 0 ). It is found to be The a 0 (980) meson lineshape is parameterized by the Flatté formula [49], where m 0 is the mass of a 0 (980) and g 2 ηπ(KK) is the coupling constant. These parameters are fixed at the values given in Ref. [50], in which m 0 = (0.990 ± 0.001)GeV/c 2 , g 2 ηπ = (0.341 ± 0.004)GeV 2 /c 4 and g 2 KK = (0.892 ± 0.022)g 2 ηπ . The ρ ηπ(KK) is the PHSP factor and is given by ρ ηπ(KK) = 2q/ √ s a . The Kπ S-wave is modeled by a parameterization from scattering data [51], which is built from a Breit-Wigner shape for the K * 0 (1430) resonance combined with an effective range parameterization for the non-resonant component with a phase shift given by with where a and r are the scattering length and effective interaction length, respectively. The parameters F (φ F ) and R(φ R ) are the magnitude (phase) for the nonresonant state and resonance terms, respectively. The parameters M , F , φ F , R, φ R , a and r are fixed to the results of the D 0 → K 0 S π + π − analysis by the BABAR and Belle experiments [51] and are given in Table VII. Table VII. The Kπ S-wave parameters, obtained from a fit to the D 0 → K 0 S π + π − Dalitz plot from the BABAR and Belle experiments [51]. The first uncertainties are statistical and the second systematic.

B. Fit Fractions and Statistical Uncertainty
The fit fractions of the individual components (amplitudes) are calculated according to the fit results. In the calculation, a large PHSP MC sample (twelve million events) with neither detector acceptance nor resolution included is used. The fit fraction (FF) for a component or an amplitude is defined as where the integration is approximated by the PHSP MC summation at the generator level,Ã n (p k j ) is either the n th amplitude (Ã n (p k j ) = c n A n (p k j )) or the n th component of a coherent sum of amplitudes (Ã n (p k j ) = c ni A ni (p k j )), and N g,ph is the number of PHSP MC events.
For the statistical uncertainty of FF, it is impractical to analytically propagate the uncertainties of the FFs from those of the magnitudes and phases. Instead, the variables are randomly perturbed within their uncertainties obtained from the fit, and the FFs are calculated to determine the statistical uncertainties. The distribution of each FF is fitted with a Gaussian function, and the width is the statistical uncertainty of the FF.

C. Results of the Amplitude Analysis
The amplitude of the D + s [S] → φρ + decay is expected to have the largest FF. Thus, this amplitude is chosen as the reference (its phase is fixed to 0, and the magnitude is fixed to 1). The notation [S] denotes a relative orbital angular momentum L = 0 between daughters in a decay, and similarly for [P ] (L = 1), [D] (L = 2). In addition, some necessary physical relations are fixed, as shown in Appendix A.
The fit to the data is initially performed with a baseline model including the amplitudes of D + s → φρ + , D + s → K * 0 K * + , D + s →K 0 1 (1270)K + (K 0 1 (1270) → K − ρ + and K * π) and D + s →K 0 1 (1400)K + (K 0 1 (1400) → K * π) decays, as the φ, ρ + ,K * 0 , K * + , K * − ,K 0 1 (1270), and K 0 1 (1400) resonances are clearly observed in the corresponding invariant mass spectra. The statistical significances (SSs) of the above amplitudes, which are determined from the changes in log-likelihood and the numbers of degrees of freedom (NDOF) when the fits are performed with and without the amplitude included, are all much larger than 4σ.
Starting from the baseline model above, the amplitudes involving f 1 (1420), f 1 (1510), η(1405), and η(1475) resonances are added to improve the fit quality of the K − K + π 0 invariant mass spectrum. The amplitudes with significances larger than 4σ are retained for the next iteration. The amplitudes of D + s → f 1 (1420)π + (f 1 (1420) → K * K) and D + s → η(1475)π + (η(1475) → a 0 0 (980)π 0 ) decays have significances larger than 5σ, and the amplitude of D + s → f 1 (1420)π + (f 1 (1420) → a 0 0 (980)π 0 ) decay improves the fit of the K − K + π 0 mass spectrum. Then, other amplitudes are tested, but only the D + s → a 0 0 (980)ρ + decay is significant (9σ). Finally, eighteen amplitudes are retained in the nominal fit, which are categorized into nine processes, as shown in Table VIII. The amplitudes of the nominal fit are listed in Table IX. Other possible processes with significance less than 3σ are listed in Appendix B.
The fit results with phases, FFs and SSs for each amplitude are shown in Table IX . Invariant mass distributions of (a) K − K + , (c) K + π 0 , (d) K − π 0 , (e) π + π 0 , (f) K − π + , (g) K − π + π 0 , (h) K − K + π + , (i) K + π + π 0 , and (j) K − K + π 0 , where the black points with error bars are data, and the red histograms are for the fit projections. The blue histograms are the backgrounds. Plot (b) shows the φ mass region of plot (a) with an expanded scale.

D. Goodness of Fit
The goodness of the fit is checked with five invariant masses, M K − π + , M π + π 0 , M K − K + , M K − π + π 0 and M K − K + π 0 , which are divided into cells of equal size. When cells contain fewer than 10 events, adjacent cells are combined until the number of events in each cell is given by (n − 1) − n par , where n par is the number of the free parameters in the fit. Overall, the value of χ 2 /NDOF is determined to be 316.8/273.

E. Systematic Uncertainties
Systematic uncertainties from the amplitude model, the background level, and the fit bias are considered. The systematic uncertainties of phases (φ) and FFs for different amplitudes are shown in Tables X and XI, respectively.

• Amplitude model
The systematic uncertainties related to the amplitude model involve the masses and widths of the intermediate resonances, the line shape of the ρ + meson, the parameters of a 0 0 (980) resonance, and the barrier effective radii of the D + s meson and other intermediate states.
-(3) The coupling constants and mass of a 0 0 (980) resonance are varied within the uncertainties given by Ref. [50]. • Experimental effects -(5) These effects are related to the PID, tracking, and reconstruction efficiency differences between data and MC, and are reflected in the factor γ ǫ in Eq. 13. The PID efficiencies are studied using clean samples of e + e − → K + K − K + K − , K + K − π + π − , K + K − π + π − π 0 , π + π − π + π − and π + π − π + π − π 0 decays, while two clean samples of the continuum process K + K − π + π − and the e + e − → K + K − π + π − π 0 decay are used for the studies of the tracking efficiencies and the π 0 reconstruction efficiency, respectively. These efficiencies are also used in the BF measurement (Section V C). The PID, tracking, and reconstruction systematic uncertainties are taken as the efficiency differences between data and MC simulation. The uncertainties associated with γ ǫ are obtained by performing alternative amplitude analyses varying PID, tracking, and reconstruction efficiencies according to their uncertainties.
• Background -(6) The MC background yields are varied within their uncertainties and the largest difference from the fits is taken as the uncertainty from the background level. The background shape B(p j ) mentioned in Eq. 6, derived from another combination of five variables (M K − K + , M π + π 0 , M K + π 0 , M K − π + , and M K − K + π 0 ) is considered and applied. The square root of the quadratic sum of these two uncertainties is taken as the background uncertainty.
• Fit bias -(7) The uncertainty due to the fit procedure is evaluated by studying signal MC samples. An ensemble of 300 signal MC samples are generated according to the nominal result in this analysis. After applying the selection criteria, each of these samples has the same size as the data sample and is used to perform the same amplitude analysis. The pull of each parameter is defined as Out(i)−In(i) σstat. (i) , where i denotes different parameters, In(i) denotes the input value as taken from the nominal fit to data, Out(i) is the value obtained from the fit to a signal MC sample and σ stat. (i) is the corresponding statistical uncertainty. For each parameter, 300 pull values are obtained and the deviation of their average from zero is taken as the systematic uncertainty.

V. BRANCHING FRACTION MEASUREMENT
To determine the absolute BF of the decay D + s → K − K + π + π 0 , the ST candidates with eight tag modes, as shown in Table III, are reconstructed and studied. Then the DT candidates are obtained by fully reconstructing the tag channels and the signal channel. Table IX. Phase, FF, and SS for the different amplitudes, labeled as I, II..., XIV. Groups of related amplitudes are separated by horizontal lines. The last row of each group gives the total fit fraction of the above components with interferences considered. The amplitudes VIII, IX, X, and XII are constructed by two sub-amplitudes with fixed relations (see Appendix A). The ρ + resonance decays to π + π 0 . The φ and a 0 0 (980) resonances decay to K − K + . TheK * 0 resonance decays to K − π + , and the K * ± resonance decays to K ± π 0 . K * π indicatesK * 0 π 0 and K * − π + . The first and second uncertainties are statistical and systematic, respectively.
The BF of the signal side is determined by where the N DT and N i ST yields are obtained from the data sample, while ε i DT and ε i ST are obtained from the generic MC samples, where i indicates the tag mode. In particular, ε i DT is determined by the amplitude analysis model used in the generic MC samples.
The signal BF B sig is determined by where i denotes the tag mode and n indicates the data sample at 4.178 GeV, 4.189-4.219 GeV or 4.226 GeV. For the numerator, n N nDT , the combined data sample is fitted to obtain the total DT data yield.

A. Event Selection
For the BF measurement, it is necessary to guarantee that the DT sample is a strict subset of the ST sample. Therefore, the ST candidates are selected ahead of the selection of DT candidates. For this measurement, the event selection criteria are relaxed or modified in order to increase the signal yield. Here, the background level does not play as crucial a role as in the amplitude analysis.  In order to reject the soft pions from D * decays, all the π mesons are required to satisfy P π > 100 MeV/c, and the χ 2 of the kinematic fit for the π 0 → γγ decay must be less than 20. The new criteria for selecting K 0 S are 487 < M π + π − < 511 (MeV/c 2 ) and that the vertex fit χ 2 must be less than 100.
For the ST selection, if there are multiple candidates for a tag mode, the one with M rec closest to the nominal M D * ± s [14] is retained. The M rec windows are given in Table II

B. Data Yields, Efficiencies and BFs
The ST yields are determined from fits to the M D − s distributions of data, as shown in Fig. 4. In the fits, an MC-simulated shape convolved with a Gaussian function is used to describe the signal shape of M D − s and a 2 ndorder polynomial function to describe the combinatorial background. For the tag mode D − s → K 0 S K − , there is some peaking background coming from D − → K 0 S π − . The shape of this background is taken from the generic MC samples and added to the fit, leaving its yield floating. For the tag mode D − s → π − η ′ , there is peaking background coming from D − s → ηπ + π − π − . The shape and the yield of this background are taken from the generic MC samples and added to the fit. The DT yields are obtained from an unbinned fit to the signal D + s mass spectrum of the combined data sample, which is shown in Fig. 5. The number of D + s → K − K + π + π 0 decays is determined to be n N nDT = 4365 ± 83. Tables XII-XIV summarize the ST efficiencies, DT efficiencies, and ST yields in data samples at 4.178-4.226 GeV.
Inserting the values of the ST and DT data yields and the ST and DT efficiencies into Eq. 33, the BF of the D + s → K − K + π + π 0 decay is measured to be B sig = (5.42 ± 0.10 stat. )% .

C. Systematic Uncertainties in the BF
The sources of the systematic uncertainties in the BF measurement are considered as follows.
• K ± meson and π ± meson tracking/PID efficiencies The ratios between data and MC efficiencies are weighted by the corresponding momentum spectra of signal MC events. The systematic uncertainties associated with tracking efficiency and PID for each charged particle are both estimated to be 0.5%. The samples used to estimate the uncertainties are mentioned in Section IV E.
• π 0 meson reconstruction efficiency According to the studies in Ref. [52], this systematic uncertainty about the π 0 reconstruction is assigned to be 2.0%.
• The numbers of ST D − s candidates The BF measurement is not sensitive to systematic    178 GeV data sample, where the dots with error bars are data, the solid blue curve shows the best fit, the red dotted curve shows the signal shape, the green dashed line shows the shape of the combinatorial backgrounds, the brown area shows the background estimated by the generic MC samples, and the pairs of pink arrows are the mass windows. In the plots for D − s → K 0 S K − and D − s → π − η ′ decays, the green dashed lines include contributions from D − → K 0 S π − and D − s → ηπ + π − π − backgrounds, respectively.
ǫ i is the average DT efficiency of tag mode i. The related uncertainty is determined to be 0.34%.
• The shape of the signal D + s mass The systematic uncertainty due to the shape of the signal is studied by fitting without the convolved Gaussian function. The difference of the DT yield is taken as the systematic uncertainty and is determined to be 0.5%.
• Background shape of the signal D + s meson For the background shape of the signal D + s , the MC-simulated shape is used to replace the nominal one, and an uncertainty of 0.75% is obtained.
• Bias of the measurement method Ten updated inclusive generic MC samples are used as fake data to estimate the possible fit bias. The BF for each sample is determined, and the relative difference between the average of BFs and the MC truth value is 0.16%, which is considered negligible.
• Amplitude model The parameters (magnitudes and phases) of the amplitude model are randomly perturbed 400 times within their statistical uncertainties according to the covariant matrix of the nominal fit to obtain the DT efficiency distribution. Then, the DT efficiency distribution is fitted with a Gaussian function. The fitted width divided by the fitted mean is 0.4% and assigned as the systematic uncertainty arising from the amplitude model.
The systematic uncertainties in the BF are summarized in Table XV. The total systematic uncertainty is ob-  Figure 5. Invariant mass distribution of the DT D + s → K − K + π + π 0 events. The black dots with error bars are data. The red dashed line represents the MC-simulated shape convolved with a Gaussian function. The green dashed line represents the MC background shape, which is fitted by a 1 st -order Chebychev polynomial. The blue solid line represents the total fitted shape.