Amplitude analysis and branching fraction measurement of the decay \boldmath $D_{s}^{+} \to K^+\pi^{+}\pi^{-}\pi^{0}$

The singly Cabibbo-suppressed decay $D_{s}^{+} \to K^+\pi^{+}\pi^{-}\pi^{0}$ is observed, using a data set corresponding to an integrated luminosity of 6.32 $\rm fb^{-1}$ recorded by the BESIII detector at the centre-of-mass energies between 4.178 and 4.226 GeV. The first amplitude analysis of $D_{s}^{+} \to K^+\pi^{+}\pi^{-}\pi^{0}$ reveals the sub-structures in this decay and determines the fractions and relative phases of different intermediate processes. The dominant intermediate process is $D_s^+ \to K^{*0}\rho^+$, with a fit fraction of $(40.5\pm2.8_{\rm{stat.}}\pm1.5_{\rm{syst.}})\%$. With the detection efficiency based on our amplitude analysis, the absolute branching fraction for $D_{s}^{+} \to K^+\pi^{+}\pi^{-}\pi^{0}$ is measured to be $(9.75\pm0.54_{\rm{stat.}}\pm0.17_{\rm{syst.}})\times 10^{-3}$.


Introduction
The hadronic decays of charmed mesons have been studied extensively in both experiment and theory since the discovery of charmed mesons in 1976 by Mark I [1,2]. However, a precise theoretical description for exclusive hadronic charmed meson decays is still challenging because the mass of charm quark is too light to adopt a sensible heavy quark expansion and too heavy to apply chiral perturbation theory [3]. Amplitude analyses and measurements of the branching fractions (BFs) for hadronic decays of charmed mesons provide valuable information about the underlying mechanism of the charmed meson decays.
Four-body hadronic decays of D + s mesons can be dominated by two-body intermediate processes [4], such as D + s → V V and D + s → AP decays, where V, A, and P denote vector, axial-vector and pseudoscalar mesons, respectively. The investigations of the D + s → V V decays have attracted a great deal of attention [5][6][7][8][9], but the experimental information about the D + s → V V decays is sparse. And the improved knowledge of BFs of D + s → AP decays, such as D + s → K 1 (1270) 0 π + and D + s → K 1 (1400) 0 π + , is important to improve the understanding of the mixing of the K 1 (1270) 0 and K 1 (1400) 0 mesons [10]. The singly Cabibbo-suppressed hadronic decay of D + s → K + π + π − π 0 is expected to be dominated by the intermediate decays D + s → K * 0 ρ + and K 0 1 π + (ρ denotes ρ(770), K * denotes K * (892) and K 1 denotes K 1 (1270)/K 1 (1400)), since the decay width calculated by external Wemission process with final states of neutral kaonioc states (i.e. K * 0 , K 0 1 ) is greater than internal W-emission process with charged kaonic states (i.e. K * + , K + 1 ) and the difference between the annihilation amplitudes could be ignored [11]. Take D + s → K * 0 ρ + and K * + ρ 0 states as an example, the tree T -diagrams and annihilation A-diagrams of these two decay modes are shown in Fig. 1 and Fig. 2, respectively. More experimental information from the amplitude analysis of this decay will offer important experimental input to improve the theory predictions and explore charge-parity (CP ) violation in the charm meson decays [9,12]. The amplitude analysis of D + s → K + π + π − π 0 also provides access to D + s → V P decays, such as D + s → ωK + . Evidence for D + s → ωK + was first reported by BESIII experiment, and the BF was measured to be (0.87 ± 0.25 stat. ± 0.07 syst. ) × 10 −3 [13], which was based on 3.19 fb −1 data samples taken at the center-of-mass energy (E cm or √ s) 4.178 GeV. The predicted value of BF (2.12 × 10 −3 ) [11] was too large compared to the experimental value of (0.87 × 10 −3 ), but after taking into account SU(3) F breaking in internal W-emission, the predicted BF now is reduced to (0.99 × 10 −3 ) [14]. Therefore, the amplitude of D + s → ωK + decay is important to investigate the W-annihilation contribution in D + s → V P decays and improve the understanding of SU(3) F flavor symmetry breaking effects in hadronic decays of charmed mesons [11,14,15]. This paper reports the first amplitude analysis and BF measurement of the decay D + s → K + π + π − π 0 , using e + e − collision data samples corresponding to an integrated luminosity of 6.32 fb −1 collected at the √ s between 4.178 and 4.226 GeV with the BESIII detector. Charged-conjugate modes are always implied throughout this paper except when discussing CP violation.

Detector and data sets
The BESIII detector is a magnetic spectrometer [16,17] located at the Beijing Electron Positron Collider (BEPCII) [18]. 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 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 identification modules interleaved with steel. The resolution of charged-particle momentum at 1 GeV/c is 0.5%, and the resolution of specific 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 in the TOF barrel region is 68 ps, while that in the end cap region is 110 ps. The end cap TOF system was upgraded in 2015 using multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [19][20][21]. About 83% of the data in this analysis benefits from the upgrade.
The integrated luminosities of different centre-of-mass energies of the data samples used in this analysis are listed in Table 1 [22][23][24]. For some aspects of the analysis, these samples are organised into three sample groups, 4.178 GeV, 4.189-4.219 GeV, and 4.226 GeV, and each of them is acquired during the same year under consistent running conditions. Since the cross section of D * ± s D ∓ s production in e + e − annihilation is about a factor of twenty larger than that of D + s D − s [25], and the D * ± s meson decays to γD ± s have a dominant BF of (93.5 ± 0.7)% [4], the signal events discussed in this paper are selected from the process e + e − → D * ± s D ∓ s → γD + s D − s . Simulated samples produced with a geant4-based [26] Monte Carlo (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 in the e + e − annihilations with the generator kkmc [27,28]. The inclusive MC sample includes the production of open charm processes, the initial state radiation production of vector charmonium(-like) states, and the continuum processes incorporated in kkmc [27,28]. The known decay modes are modeled with evtgen [29,30] using BFs taken from the Particle Data Group (PDG) [4], and the remaining unknown charmonium decays are modeled with lundcharm [31,32]. Final state radiation from charged final state particles is incorporated using photos [33].

Event selection
To obtain the signal samples with high purity, we adopt the double tag method [34] in this analysis. In this method, a single-tag (ST) candidate requires only one of the D ± s mesons to be reconstructed via a hadronic decay; a double-tag (DT) candidate has both D + s D − s mesons reconstructed via hadronic decays, where one D s meson is reconstructed via the signal mode and the other via any of the tag modes. The D ± s candidates are constructed from individual π ± , π 0 , K ± , K 0 S , η and η particles, with the following selection criteria. All charged tracks reconstructed in the MDC must satisfy |cos θ| < 0.93, where θ is the polar angle of a charged track with respect to the positive direction of the MDC axis. For charged tracks not originating from K 0 S decays, the distance of closest approach to the interaction point is required to be less than 10 cm along the beam direction and less than 1 cm in the plane perpendicular to the beam. Particle identification (PID) for charged tracks is performed by using the dE/dx measured by the MDC and the flight time in the TOF. The confidence level for pion and kaon hypotheses (CL K and CL π ) 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. The distances of the charged tracks to the interaction point along the beam direction are required to be less than 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 in the interval of |M π + π − − m K 0 [4]. The decay lengths of the K 0 S candidates are required to be twice greater than its uncertainty. Photon candidates are identified by their 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). The minimum opening angle between the position of each shower in the EMC and the closest extrapolated charged track is required to be greater than 10 • to exclude the showers originating from tracks. The difference between the EMC time and the event start time is required to be within [0, 700] ns to suppress electronic noises and showers unrelated to the event.

Further selection criteria
To obtain data samples with high purities for the amplitude analysis, the following dedicated selection criteria are imposed on the signal candidates.
The seven-constraint kinematic fit to the process e + e − → D * ± s D ∓ s → γD + s D − s , where the D − s decays to one of the tag modes and the D + s decays to the signal mode, is required to converge. In addition to the constraints of four-momentum conservation in the e + e − centre-of-mass system, the invariant masses of (γγ) π 0 , tag D − s , and D * ± s candidates are constrained to the corresponding known masses [4]. The combination with the minimum χ 2 is chosen, assuming that D * + s decays to D + s γ or D * − s decays to D − s γ. In order to ensure that all candidates fall within the phase-space boundary, the constraint of the signal D + s mass is added to the kinematic fit and the updated four-momenta from this kinematic fit are used for the amplitude analysis.
The energy of the transition photon from D * + s → γD + s is required to be less than 0.18 GeV. The recoiling mass against this photon and the signal D + s candidate is required to lie within the range of [1.955, 1.995] GeV/c 2 .
There is a wrong-combination background from D 0 → K − π + π 0 versusD 0 → K + π + π − π − , which fakes D + s → K + π + π − π 0 versus D − s → K − K 0 S , K 0 S → π + π − by exchanging a π 0 and π − . It also fakes D + s → K + π + π − π 0 versus D − s → K + K − π − by identifying a π + from the D 0 as a K + and exchanging a π − from theD 0 with the π 0 from D 0 . This background is excluded by rejecting the events which simultaneously satisfy |M K − π + π 0 − M D 0 | < 75 MeV/c 2 and |M K + π + π − π − − MD0| < 50 MeV/c 2 , where M D 0 is the known D 0 mass [4]. There is also a wrong-combination background from by exchanging π + and π − , then adding a π 0 . This background is excluded by rejecting the events which simultaneously satisfy where M D + is the known D + mass [4]. Figure 3 shows the fits to the invariant-mass distributions of the accepted signal D + s candidates, M sig , for various data samples. In the fits, the signal is described by a MC-simulated shape convolved with a Gaussian function and the background is described by a simulated shape derived from the inclusive MC sample. Then, a mass window,

Fit method
An unbinned maximum likelihood fit is used in the amplitude analysis of D + s → K + π + π − π 0 , the likelihood function L is constructed with a signal-background combined probability density function (PDF), which depends on the momenta of the four final state particles. The likelihood is written as where i indicates the data sample groups. The p k denotes the four-momenta of the final state particles, where k denotes the k th event in the data sample i. The N D,i is the number of candidates in the data sample i, f S (f B ) is the signal (background) PDF and w i is the purity of the signal discussed in Sec. 4.1.
The signal PDF is given by where (p) is the detection efficiency in bins of a five-dimensional space of two-and threebody invariant masses, and R 4 is the four-body phase space. The total amplitude M is modeled with the isobar model, which is the coherent sum of the individual amplitudes of intermediate processes, given by M = ρ n e iφn A n , where the magnitude ρ n and phase φ n are the free parameters to be determined by the fit. The amplitude of the n th intermediate process (A n ) is given by where k is the index of the k th event of the signal MC sample, and N MC is the number of the selected MC events. The M g (p) is the signal PDF used to generate the signal MC sample in the MC integration. The normalization integral for background is also realised by a MC integration method like Eq. 4.4, To account for the bias caused by differences in tracking, PID efficiencies and π 0 reconstruction between data and MC simulation, each signal MC event is weighted with a ratio, γ (p), and it is calculated as where j denotes the final four daughter particles, j,data (p) and j,MC (p) are the tracking, PID and π 0 reconstruction efficiencies as a function of the momenta of the daughter particles for data and MC simulation, respectively. By weighting each signal MC event with γ , the MC integration is modified to be (4.7) The background PDF is given by where B (p) = B(p)/ (p) is the efficiency-corrected background shape. The background shape B(p) is derived by using a multi-dimensional kernel density estimator [36] named RooNDKeysPdf implemented in RooFit [37], which models the distribution of an input dataset as a superposition of Gaussian kernels using background events in the M sig signal region from the inclusive MC sample. The M K + π − , M K + π 0 , M π + π − , M π + π 0 and M K + π − π 0 distributions of the inclusive MC events outside the M sig signal region are compared to these distributions from the data to check their validity. The distributions of background events from the inclusive MC sample within and outside the M sig signal region are also examined. They are compatible with each other within statistical uncertainties.

Blatt-Weisskopf barrier factors
For the process a → bc, the Blatt-Weisskopf barrier factors [38], X L (p), are parameterised as a function of the angular momenta L and the momenta q of the final state particle b or c in the rest system of a. They are taken as where z = qR, z 0 = q 0 R and the effective radius of the barrier R is fixed to 3.0 GeV −1 for the intermediate resonances and 5.0 GeV −1 for the D + s meson. The momentum q is given by the value of q 0 is that of q when s a = m 2 a and the s a (s b , s c ) denotes the invariant-mass squared of the particle a(b, c).
Considering the obvious mass deviation reported in the PDG [4], the mass and width of K 1 (1270) 0 are fixed to 1289 MeV/c 2 and 116 MeV, respectively, from results obtained by the LHCb experiment [39]. The ρ resonances are parameterised by the Gounaris-Sakurai (GS) lineshape [40], which is given by . (4.12) The function f (m) is given by and The normalisation condition at P GS (0) fixes the parameter The K-Matrix parametrisation is used to describe the π + π − S-wave. Detailed descriptions of the K-matrix formalism can be found in various references [41][42][43][44]; parameters used are summarised in Tables 4 and 5. We use the "K-matrix amplitude" to describe the amplitude of channel u (u = 1 − 5 denote the channels ππ, KK, 4π, ηη, ηη ) in the form Here the vectorP describes the production of bare states and the non-resonant production of meson pairs, while the term (I − iKρ) −1 describes their re-scattering.
The scattering matrixK can be parameterised as a combination of the sum of N poles with real bare masses m α , together with slowly-varying non-resonant parts (SVPs): Table 5. ππ S-wave P-vector parameters obtained from the D 0 → K 0 S π + π − Dalitz plot distribution from D * + → D 0 π + . P-vector parameters f prod 1v , for v = 1, are defined as f prod 1v /f prod 11 .
Component a r φ r (deg) β 1 9.3 ± 0.4 −78.7 ± 1. describe a smooth part for the K-matrix elements and m 2 0 , s A , and s A 0 are real constants of order unity. All these parameters are taken from Ref. [42]. Here s denotes the invariant mass squared of π + π − .
The production vectorP vector is parameterised in a form analogous to theK matrix and it is given byP where β α and f prod v are complex production constants for the poles and non-resonant SVPs, respectively, both of them depend on the final state channel.
The Kπ S-wave is modeled by a parameterisation from scattering data [45], which is described by a K * 0 (1430) Breit-Wigner together with an effective range non-resonant component with a phase shift. It is given by with where the parameters F (φ F ) and R(φ R ) are the magnitude (phase) for non-resonant state and resonance terms, respectively. The parameters a and r are the scattering length and effective interaction length, respectively. We fix these parameters (M, Γ, F, φ F , R, φ R , a, r) to the results obtained from the amplitude analysis to a sample of D 0 → K 0 S π + π − by the BABAR and Belle experiments [45]; these parameters are summarised in Table 6.

Spin factors
For the process a → bc, the four-momenta of the particles a, b, and c are denoted as p a , p b , and p c , respectively. The spin projection operators [35] are defined as µνµ ν (a) = µ ν (a) . (4.21) The pure orbital angular-momentum covariant tensors are given bỹ where r a = p b − p c . The spin factors S(p) used are listed in Table 7. The tensor describing the D + s decays with orbital angular-momentum quantum number l is denoted byT (l)µ and that of intermediate a → bc decay is denoted byt (l)µ , and theT (l)µ has the same definition ast (l)µ in Ref. [35].

Fit results
Using the method described in Sec. 4.2, we perform an unbinned maximum likelihood fit to the D + s → K + π + π − π 0 decay channel. The fit is performed in steps, by adding resonances one by one. The corresponding statistical significance for the newly added amplitude is calculated with the change of the log likelihood value, taking the change of the number of the degrees of freedom into account.
For the amplitude fits, the magnitude and phase of the D + s → K * 0 ρ + reference amplitude are fixed to 1 and 0, respectively, while those of the other amplitudes are Decay chain floated in the fit. The amplitudes for D + s → K * + ρ 0 and D + s → K + ω are also included, as they are clearly observed in the corresponding invariant mass spectra. After testing each, , D + s → (K + π 0 ) V ρ 0 , and D + s → (K + π 0 ) S−wave (π + π − ) S−wave are added since each has a statistical significance greater than 3σ. Considering the isospin relationship in hadron decays, some Clebsch-Gordan relations are fixed, with details in Appendix A. A full list of other allowed contributions (based on known states) with statistical significances less than 3σ are listed in Appendix B.
The fit fraction (FF) for the n th amplitude is computed numerically with generator-level MC events with the definition as where N gen is the number of phase space MC signal events at generator level. The sum of these FFs is generally not unity due to net constructive or destructive interference. Interference IN between the n th and n th amplitudes is defined as (4.24) In order to determine the statistical uncertainties of FFs, the amplitude coefficients are randomly selected by a Gaussian-distributed set by the fit results according to their uncertainties and the covariance matrix. The distribution of each FF is fitted with a Gaussian function whose width is then taken as the uncertainty of this FF. The phases, FFs, and statistical significances(Stat.Signi) for different amplitudes are listed in Table 8. The mass projections of the nominal fit are shown in Fig. 4.

Systematic uncertainties for the amplitude analysis
The systematic uncertainties for the amplitude analysis, summarised in Table 9, are now detailed.
i Amplitude model: The masses and widths of resonances are varied by their corresponding uncertainties [4,39]. The GS lineshape of ρ is replaced with the RBW formula. The coupling constants of the ππ S-wave model are varied within their uncertainties given in Ref. [44]. The changes of the phases and FFs are assigned as the associated systematic uncertainties. Since replacing the lineshape of the Kπ S-wave model from BABAR with the K-matrix formula [46] results in different normalisation factors, the effect on the phase of the amplitude related to Kπ S-wave is not considered for this source.
ii R values: We assume the distribution of values for barrier effective radius (R), as defined in Sec. 4.2.1, as a uniform distribution. The systematic uncertainties associated with R are estimated by repeating the fit procedure by varying the R of both the intermediate state and D + s mesons by R/ iii Background: The uncertainty from background size is studied by varying the fractions of signal (equivalent to the fractions of background), i.e. w i in Eq. 4.1, within their corresponding statistical uncertainties. Another source is the simulation of background shapes. First, alternative MC shapes where the relative fractions of the dominant backgrounds from e + e − → qq (q = u, d, s) and non-D * ± s D ∓ s open charm processes are varied by the statistical uncertainties of their cross sections are used. Second, the background PDF is extracted using other five variable input combinations with varied smoothing parameters of RooNDKeysPDF [37]. iv Simulation effects: To estimate the uncertainties caused by γ , as defined in Eq. 4.6, an amplitude fit is performed by varying efficiencies of PID, tracking and π 0 reconstruction according to their uncertainties.
v Fit bias: The uncertainty from the fit process is evaluated by studying 600 signal MC samples with the size equal to the data sample size that are generated to check the pull. The pull variables, (V input −V fit )/σ fit , are defined to evaluate the corresponding uncertainty, where V input is the input value in the generator, V fit and σ fit are the output value and the corresponding statistical uncertainty, respectively. Expected to be the standard normal distribution for an unbiased fit, the distributions of pull values for the 600 sets of sample are fitted with a Gaussian function. The fitted mean values for the pulls of FFs of D + s [P ] → K * 0 ρ + and D + s → (K + π 0 ) S−wave (π + π − ) S−wave deviate from zero by larger than 3 times of the standard deviation. We correct all resonances' FFs and phases by the fitted mean values, and assign the uncertainty of the fitted mean values as the corresponding systematic uncertainties.

Branching fraction measurement
On top of the selection criteria described in Sec. 3, the momenta of all pions are further required to be greater than 100 MeV/c to exclude soft pions from D * decays. The best tag candidate is chosen with M rec closest to m Ds if there are multiple ST candidates. The yields for various tag modes are obtained from the fits to the corresponding M tag distributions and the results are summarised in Table 10. As an example, the fits to the data sample at √ s = 4.178 GeV are shown in Fig. 5. In the fits, the signal is modeled by an MC-simulated shape convolved with a Gaussian function to take into account the data-MC resolution difference. The background is described by a second-order Chebychev polynomial. Inclusive MC studies show that there is no significant peaking background in any tag mode, except for D − → K 0 S π − and D − s → ηπ + π − π − faking the D − s → K 0 S K − and D − s → π − η tags, respectively. Therefore, the MC-simulated shapes of these two peaking background sources, with the yields included as free parameters, are added to the fits, respectively.
Once a tag mode is identified, we select the signal decay D + s → K + π + π − π 0 . In the case of multiple candidates, the DT candidate with the average mass, (M sig + M tag )/2, closest to m Ds is retained.
For a given single tag mode, the ST and DT yields can be written as:   tag modes and sample groups gives the total DT yield:

Tag mode
where α represents tag modes in the i th sample group. Therefore, the BF of the signal decay can be determined by 4) where N ST α,i and ST α,i are obtained from the data and inclusive MC samples, respectively, while DT α,sig,i is determined with signal MC samples generated based on our amplitude analysis. The branching ratio B π 0 →γγ has been introduced as it is not included in the MC generation. The DT yield N DT total is found to be 776 ± 43 from the fit to the M sig distribution of the selected D + s → K + π + π − π 0 candidates. The fit result is shown in Fig. 6, where the signal shape is modeled by an MC-simulated shape convolved with a Gaussian function to take into account the data-MC resolution difference. The background shape is derived from the inclusive MC sample. After correcting for the differences in K + and π ± tracking, PID and π 0 reconstruction efficiencies between data and MC simulation, we determine the BF of D + s → K + π + π − π 0 to be (9.75 ± 0.54 stat. ± 0.17 syst. ) × 10 −3 according to Eq. 5.4.
The BFs for the charge-conjugated modes D + s → K + π + π − π 0 and D − s → K − π − π + π 0 , which are labeled as B(D + s ) and B(D − s ), are measured to be (9.10±0.71 stat. ±0.16 syst. )×10 −3 and (10.39 ± 0.79 stat. ± 0.18 syst. ) × 10 −3 , respectively. The asymmetry of the BFs, , is determined to be (6.5 ± 5.4 stat. ± 0.7 syst. )%. No significant CP violation is observed with the current sample size. Note that the systematic uncertainties due to pion tracking and PID, π 0 reconstruction are canceled in the A CP calculation. The systematic uncertainties in the BF measurement are discussed as follow.
• Background shape: Since the e + e − → qq and non-D * ± s D ∓ s open charm processes are the major background sources, an alternative MC-simulated background shape is obtained by varying the relative fractions of the background components from these processes by the statistical uncertainties of their cross sections. The largest change, 0.5%, is taken as the related uncertainty.
• π 0 reconstruction: The π 0 reconstruction efficiency is investigated by using a control sample of the pro-cess e + e − → K + K − π + π − π 0 . The data-MC efficiency ratio for π 0 reconstruction is estimated to be 0.995 ± 0.008. After correcting the efficiency by this factor, we assign 0.8% as the systematic uncertainty.
• MC sample size: The uncertainty due to the limited MC sample size is obtained by α (f α δ α α ) 2 , where f α is the tag yield fraction, and α and δ α are the signal efficiency and the corresponding uncertainty of tag mode α, respectively. The uncertainty corresponding to MC statistics is 0.2%.
• Amplitude analysis model: The uncertainty arising from the amplitude analysis model is estimated by varying the model parameters based on their error matrix. The distribution of 600 efficiencies resulting from this variation is fitted by a Gaussian function. The fitted width divided by the mean value, 0.4%, is taken as an uncertainty.
All of the systematic uncertainties are summarised in Table 11. Adding them in quadrature results in a total systematic uncertainty of 1.7%.
The dominant intermediate process is determined to be D + s → K * 0 ρ + , with a fraction of (40.5 ± 2.8 stat. ± 1.5 syst. )%. The decay D + s → K + ω is observed with a significance greater than 10σ and its BF is measured to be (0.95±0.12 stat. ±0.06 syst. )×10 −3 , which is consistent with the BESIII result (0.87 ± 0.24 stat. ± 0.08 syst. ) × 10 −3 [13] within 1σ, but the precision is improved by a factor of 2.1. Information about the two K 1 states in this decay provides inputs to further investigations of the mixing between these two axial-vector kaon states [15]. The asymmetry for the BFs of the decays D + s → K + π + π − π 0 and D − s → K − π − π + π 0 is determined to be (6.6 ± 5.4 stat. ± 0.7 syst. )%. No evidence for CP violation is found under the current sample size.

A Clebsch-Gordan relation
Considering the isospin relationship in hadron decays, some amplitudes are fixed by Clebsch-Gordan relations, as listed in Table 13. The amplitudes with fixed relations share the same magnitude (ρ) and phase (φ).

B Other intermediate processes tested
Some other tested amplitudes with significance less than 3σ are listed below, the value in each of brackets corresponds to the significance.