Amplitude analysis and branching-fraction measurement of $D^{+}_{s}\rightarrow K^{0}_{S}K^{-}\pi^{+}\pi^{+}$

Using 6.32 fb$^{-1}$ of $e^+e^-$ collision data collected by the BESIII detector at the center-of-mass energies between 4.178 and 4.226 GeV,~an amplitude analysis of the $D^{+}_{s}\rightarrow K^{0}_{S}K^{-}\pi^{+}\pi^{+}$ decays is performed for the first time to determine the intermediate-resonant contributions. The dominant component is the $D_s^+ \to K^*(892)^+\overline{K}^*(892)^0$ decay with a fraction of $(40.6\pm2.9_{\rm stat}\pm4.9_{\rm sys})$%. Our results of the amplitude analysis are used to obtain a more precise measurement of the branching fraction of the $D^{+}_{s}\rightarrow K^{0}_{S}K^{-}\pi^{+}\pi^{+}$ decay, which is determined to be $(1.46\pm0.05_{\rm stat}\pm0.05_{\rm sys}$)%.

Using 6.32 fb −1 of e + e − collision data collected by the BESIII detector at the center-of-mass energies between 4.178 and 4.226 GeV, an amplitude analysis of the D + s → K 0 S K − π + π + decays is performed for the first time to determine the intermediate-resonant contributions. The dominant component is the D + s → K * (892) + K * (892) 0 decay with a fraction of (40.6 ± 2.9stat ± 4.9sys)%. Our results of the amplitude analysis are used to obtain a more precise measurement of the branching fraction of the D + s → K 0 S K − π + π + decay, which is determined to be (1.46 ± 0.05stat ± 0.05sys)%.

I. INTRODUCTION
The decay D + s → K 0 S K − π + π + is usually used as a "tag mode" for measurements related to the D + s meson [1][2][3][4][5] due to its large branching fraction and low background contamination. The inclusion of chargeconjugate states is implied throughout the paper. In 2013 the CLEO Collaboration reported its branching fraction B(D + s → K 0 S K − π + π + ) to be (1.64 ± 0.07 ± 0.08)%, based on a data sample corresponding to an integrated luminosity of 586 pb −1 of e + e − collisions at a center-ofmass energy (E cm ) of 4.17 GeV [6]. The measurement was limited by the sample size and lack of knowledge of the intermediate processes. In addition, the branching fraction of D + s → K * (892) + K * (892) 0 was determined by the ARGUS Collaboration [7] more than twenty years ago, who claimed the contribution of D + s → K * (892) + K * (892) 0 in the D + s → K 0 S K − π + π + decays is almost 100%. The ARGUS measurement suffers from low statistics and large uncertainties in the branching frac-tion of the reference decay D + s → φ(1020)π + . An amplitude analysis of the D + s → K 0 S K − π + π + decays is necessary to investigate the resonant contributions, and thereby reduce the systematic uncertainties of its branching fraction and for providing input to measurements where amplitude information is essential.
It is well known that two-body modes dominate D + s decays [8]. The majority of the observed two-body decays have pseudoscalar-pseudoscalar or pseudoscalar-vector mesons in the final states. Among various kinds of D + s decay modes, vector-vector final states are of special interest. The ratios between different orbital angular momenta of the two vector mesons for the dominant quasi-two-body decay D + s → K * (892) + K * (892) 0 provide valuable information on CP violation with Tviolating triple-products [9]. In addition, several mesons with J P = 0 − , 1 + are reported in the mass region between 1.2 and 1.6 GeV/c 2 and decay to the (KKπ) 0 final state [10][11][12][13]. These are the η(1295), η(1405), η(1475), f 1 (1285), f 1 (1420) and f 1 (1510), although many of these states are not well established. This paper presents the first amplitude analysis and an improved branching-fraction measurement of the D + s → K 0 S K − π + π + decay with data samples corresponding to a total integrated luminosity of 6.32 fb −1 collected by the BESIII detector at E cm between 4.178 and 4.226 GeV.

II. DETECTOR AND DATA SETS
The detailed description of the BESIII detector can be found in Ref. [14]. It is a magnetic spectrometer located at the Beijing Electron Positron Collider (BEPCII) [15]. 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 charged-particle momenta resolution at 1.0 GeV/c is 0.5%, and the specific energy loss (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. The endcap TOF was upgraded in 2015 with multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [16].
The data samples used in this paper were accumulated in the years 2013, 2016 and 2017 with E cm of 4.226, 4.178 and 4.189−4.219 GeV, respectively. Generic Monte Carlo samples (GMC) that are 40 times larger than the data sets are produced with the GEANT4-based software [17]. The production of open-charm processes directly via e + e − annihilation is modeled with the generator conexc [18], which includes the effects of the beam energy spread and initial state radiation (ISR). The ISR production of vector charmonium states and the continuum processes are incorporated in kkmc [19]. The known decay modes are generated using evtgen [20], which assumes the branching fractions reported by the Particle Data Group (PDG) [8]. The remaining unknown decays from the charmonium states are generated with lundcharm [21]. The final state radiation from charged tracks are simulated by the photos package [22]. The GMC is used to estimate background and optimize selection criteria.
More than 10 million simulated events are generated with an uniform distribution in the phase space of the D + s → K 0 S K − π + π + decay to perform the normalization in the amplitude fit. Preliminary parameters of the amplitude model are obtained from an initial fit to the data. A signal Monte Carlo (SMC) sample is generated according to the preliminary parameters and is used to validate the fit performance and to estimate the detector efficiency. A final determination of the fit parameters is obtained by fitting the data using the SMC sample for the normalization.

III. EVENT SELECTIONS
The production of D ± s candidates is dominated by the process e + e − → D * + s D − s , where the D * + s meson decays to either γD + s or π 0 D + s with branching fractions of (93.5±0.7)% and (5.8±0.7)% [8], respectively. A sample of D − s mesons is reconstructed first, with nine D − s prominent hadronic decay modes, as shown in Table I, and is referred to as the "single tag (ST)" candidates. The signal decay D + s → K 0 S K − π + π + is reconstructed by selecting two π + , one K − and one K 0 S candidates from the unused tracks in each ST event, and is referred to as the sample of "double tag (DT)" candidates.
All charged tracks reconstructed in the MDC must satisfy | cos θ| < 0.93, where θ is the polar angle with respect to the direction of the positron beam. Except for K 0 S daughters, they must originate from the interaction point with a distance of closest approach less than 1 cm in the transverse plane and less than 10 cm along the beam direction. The dE/dx information in the MDC and the time-of-fight information from the TOF are combined and used for particle identification (PID) by forming confidence levels CL K(π) for kaon (pion) hypotheses. Kaon (pion) candidates are required to satisfy CL K(π) > CL π(K) .
The K 0 S candidates are selected by looping over all pairs of tracks with opposite charges, whose distances to the interaction point along the beam direction are within 20 cm. A primary vertex and a secondary vertex [23] are reconstructed and the decay length between the two is required to be greater than twice its uncertainty. Since the combinatorial background is low, this requirement is not applied for the D − s → K 0 S K − decay. The invariant mass M π + π − is required to be in the region [0.487, 0.511] GeV/c 2 . To prevent an event being retained by both the D − s → K 0 S K − and D − s → K − π + π − selections, the value of M π + π − is required to be outside of the mass range [0.487, 0.511] GeV/c 2 for the D − s → K − π + π − decay.

A. Selections for Amplitude Analysis
The tagged D − s candidates are constructed from the π + , K − , η, η ′ , K 0 S and π 0 mesons, while the signal D + s candidates are reconstructed from the K 0 S , K − and two π + mesons. The requirements on the recoiling mass of the D + s , M rec , and the mass of the tagged D − s , M tag , are summarized in Table I. The recoiling mass is calculated as follows: (1) Here, p D + s is the three momentum of the D + s candidate and m D + s is its nominal mass [8]. Mtag for individual single-tagged modes. The K 0 S , π 0 (η) and η ′ mesons decay to π + π − , γγ and π + π − η final states, respectively.

Ecm (GeV)
Mrec ( Kinematic fits are performed of the process e + e − → D * ± s D ∓ s → γD ± s D ∓ s with the photon assigned to each charm meson in turn, and the χ 2 of the fit being used to decide between the D * + s and D * − s hypotheses. The fits include constraints from four-momentum conservation in the e + e − system, and also constrain the invariant masses of K 0 S , D * ± s and tag-side D ± s candidates to their nominal masses [8]. In order to ensure that all candidates fall within the kinematic boundary of the phase space, we perform a further kinematic fit in which the signal D ± s mass is constrained to its nominal value, and the updated four-momenta are used for the amplitude analysis. To suppress the background where the π − from the signal decay and the π + from the tag modes are exchanged, which fakes the signal and the same tag mode but with opposite charges, we perform kinematic fits with D + s and D − s mass constraints for the two cases, and select the one with the smaller χ 2 . To reduce the background coming from Figure 1 shows the invariant mass distributions of the signal D + s , M sig , in data and the fit results. The signal distribution is modeled with the simulated shape convolved with a Gaussian function and the background is described by a first-order Chebychev polynomial. The fitted yields for the signal are 744±28, 415±21 and 159±13 in the invariant mass range [1.951, 1.987] GeV/c 2 , with purities of (94.7±0.5)%, (96.2±0.7)% and (93.9±1.2)% for the data samples taken at E cm = 4.178, 4.189 − 4.219 and 4.226 GeV, respectively. The candidates falling in the D + s mass region are retained for the amplitude analysis. We compare the background yield and various distributions of the events outside the signal region between data and GMC. The yield and distributions are found to be consistent within the statistical uncertainties. The background events in the signal region from GMC are used to estimate the background contributions in data.

B. Likelihood Function
An unbinned-maximum-likelihood method is applied to determine resonant contributions in the D + s → K 0 S K − π + π + decays. The likelihood function is constructed with a probability density function (PDF) of the momenta of the four daughter particles. The amplitude of the n th intermediate state (A n ) is where the indices 1, 2 and 3 correspond to the two subsequent intermediate resonances and the D + s meson. S is the spin factor constructed with the covariant tensor formalism [24], F is the Blatt-Weisskopf barrier factor and P is the propagator of the intermediate resonance. The total amplitude M is a coherent sum of the amplitudes of intermediate processes, where c n = ρ n e iφn are complex coefficients to be determined from the fit to data. The signal PDF f S (p j ) is given by where ǫ(p j ) is the detection efficiency parameterized in terms of the final four-momenta p j and j refers to the different particles in the final states. R 4 (p j ) is the standard element of the four-body phase space. The normalization is determined from the simulated events, where k sim runs from 1 to N sim , the total number of simulated events. M gen (p j ) is the PDF used to generate the simulated samples. The normalization takes into account the difference in detector efficiencies for PID and tracking between data and simulation by assigning a weight to each simulated event where i denotes the four daughter particles. The normalization is then given by The total PDF f T (p j ) is where w is the purity of the signal described by a constant parameter in the fit. We factorize out ǫ(p j ) from f T (p j ) as ǫ(p j ) is independent of the fit variables. Its contribution enters into the normalization and the background PDF. As a consequence, the combined PDF becomes where B ǫ (p j ) = B(p j )/ǫ and the background PDF B(p j ) is parameterized using RooNDKeysPdf [25]. The normalization in the denominator of the background term is calculated as Finally the log-likelihood is written as and data samples collected at different E cm are fitted simultaneously.

C. 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 [24] are defined as The pure orbital angular-momentum covariant tensors are given byt where r a = p b − p c . The spin factors S(p) used in this paper are constructed from the spin projection operators and pure orbital angular-momentum covariant tensors and are listed in Table II.

D. Blatt-Weisskopf Barrier Factors
For the process a → bc, the Blatt-Weisskopf barrier factor F L (p j ) is parameterized as a function of the angular momentum L and the momentum q of the daughter b or c in the rest system of a, where z = qR. R is the effective radius of the barrier, which is fixed to 3.0 GeV −1 × c for the intermediate resonances and 5.0 GeV −1 × c for the D + s meson [26]. The momentum-transfer squared is where s a,b,c are the invariant-mass squared of particles a, b, c, respectively. The X L (q) factors are given by

E. Propagators
The propagators for the resonances K * (892) + , K * (892) 0 , η(1295), η(1405), η(1475), f 1 (1285), f 1 (1420) and f 1 (1510) are modeled by the relativistic Breit-Wigner function, which is given by where m 0 and Γ(m) are the mass and width of the intermediate resonance, and q 0 is the value of q when s a = m 2 0 . The a 0 (980) contribution is parameterized as the Flatté formula , (18) where ρ ηπ (s a ) and ρ KK (s a ) are the Lorentz-invariant phase-space factors defined as 2q/ √ s a , and the coupling constants g 2 ηπ = 0.341 ± 0.004 GeV 2 /c 4 and g 2 KK = (0.892 ± 0.022)g 2 ηπ [27]. We use the same parameterization to describe the Kπ S-wave as Ref. [28], which is extracted from scattering data [29]. The model is built with a Breit-Wigner shape for the K * (1430) 0 and an effective range parameterization for the non-resonant component, with where a and r are the scattering length and effective interaction length, respectively. The parameters F (φ F ) and R(φ R ) are the magnitudes (phases) for the nonresonant term and the resonant contribution, 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 Collaborations [28].

F. Fit Fraction
The fit fraction (FF) for a quasi-two-body contribution is independent of any normalization and phase conventions in the amplitude formalism, and hence provides a more useful measure of amplitude strengths than the magnitudes of each contribution. The definition of the FF for the n th contribution is where the integration is approximated by the sum of the simulated events generated flatly over the phase space and without any efficiency effects included.
To estimate the statistical uncertainties on the FFs, the calculation is repeated by randomly varying the fit parameters according to the error matrix. The resulting distribution of each FF is fitted with a Gaussian function, whose width gives the corresponding statistical uncertainty.

G. Fit Results
In the fit, the magnitude (ρ) and phase (φ) of D +  Table III. The magnitude and correlation matrix are provided in the supplemental material [31]. The projections for the nine invariant-mass distributions are shown in Fig. 2.
To validate the fit performance, 300 sets of SMC samples with the same size as the data samples are generated according to the nominal fit results in this analysis. Each sample is analyzed with the same method as for data. The pull value is given by V pull = (V fit − V input )/σ fit , where V input is the input value in the generator, V fit and σ fit is the output value and the corresponding statistical uncertainty, respectively. The resulting pull distributions are fitted with Gaussian distributions. The fitted mean value of the pull distribution for the FF of D + s → η(1475)π + , η(1475) → (K 0 S π + ) S-wave K − deviates from zero by more than 3.0σ, we correct its FF according to the deviation and the uncertainty of the FF.

H. Systematic Uncertainties
The systematic uncertainties for the amplitude analysis are studied in the following categories.
• Kπ S-wave model. The fixed parameters of the model are evaluated by varying the input values within ±1σ according to Ref. [28].
• Lineshape of a 0 (980). The Flatté parameters are shifted by ±1σ based on the values given in Ref . [27].
• Effective barrier radius. The barrier radius are varied within ±1 GeV −1 × c for intermediate resonances and the D + s meson.
• Masses and widths of the resonances considered. The masses and widths are shifted by ±1σ based on their values from the PDG [8].
• Background estimation. We shift the fractions of the signal in Eq. 9 according to the uncertainty associated with the background estimation and take the largest shift as the systematic uncertainty.
• Experimental effects. To determine the systematic uncertainty due to tracking and PID efficiencies, we alter the fit by shifting the γ ǫ in Eq. 6 within its uncertainty, and the change of the nominal fit result is taken as the systematic uncertainty.
• Neglected resonances. The intermediate processes with statistical significance less than four standard deviations are added one-by-one to the nominal contributions. For each parameter, the maximum difference with respect to the nominal fit result is taken as the corresponding systematic uncertainty.
• Fit uncertainties. The fitted widths from the pull distributions described in Sec. IV G are consistent with 1.0 within 2.0σ. Therefore, the fit uncertainties are estimated properly and no systematic uncertainty is assigned from this source.
All of the systematic uncertainties of the φ and FFs are listed in Table IV. The total systematic uncertainties are obtained by adding the above systematic uncertainties in quadrature.

A. Yields and Efficiencies
The selection criteria of the tagged D − s and signal D + s candidates are the same as in Sec. III, except for the following requirements: (I) the requirement of the secondary vertex fit for K 0 S from the tag modes is removed, while that for the signal is retained; (II) a further requirement of p π ± /π 0 > 0.1 GeV/c is added to remove the soft π ± /π 0 directly from D * ± /D * 0 decays; (III) the tagged D − s candidates are reconstructed by looping over all their daughter tracks to form different combinations. If there are multiple candidates from the same event, the one with M rec closest to the D * ± s mass is accepted; (IV) at least one of the D + s /D − s candidates must satisfy III. The φ, FFs and significances for different resonant contributions, labeled as I, II, III, · · · , XIII, respectively. The first and second uncertainties are the statistical and systematic uncertainties, respectively. Here K * (892) + , K * (892) 0 and a0(980) − denote K * (892) + → K 0 S π + , K * (892) 0 → K − π + and a0(980) − → K 0 S K − , respectively, while K(892) * K indicates K * 0 K 0 S and K * (892) + K − . The FF of IV (IIX) term is the sum of I, II and III (VIII and IX) terms after considering the interference.

Label
Component φ FF(%) Significance (σ) In each fit, the signal shape is modeled using the simulated shape convolved with a Gaussian function, whose resolution and mean are free parameters, and the background is described with a second-order Chebychev polynomial. These fits give a total ST yield of N ST = 550496 ± 2411. The M tag distributions at E cm = 4.178 GeV are shown in Fig. 3 as an example. The total DT signal yield, N tot DT , is determined to be 1332 ± 42, as shown in Fig. 4. The fits to the M sig distribution for GMC are performed to estimate the corresponding ST efficiencies (ǫ ST ). The DT efficiencies (ǫ DT ) are determined by GMC, in which our amplitude analysis model is taken for the generation of the signal mode.

B. Tagging Technique and Branching Fraction
The branching fraction for the signal mode is given by where the indices i and j denote the i th tag mode and the j th center-of-mass energy point. The N ij ST and ǫ ij

ST(DT)
are the number of the ST candidates and the corresponding ST (DT) detection efficiency.

C. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties for the branching fraction measurement are studied in the following categories.
• K 0 S reconstruction efficiency. The uncertainty for the K 0 S reconstruction efficiency is assigned as 1.5% per K 0 S , obtained using control samples of J/ψ → K 0 S K ± π ∓ and φK 0 S K ± π ∓ decays.
• Fit to the DT M sig distribution. The uncertainty associated with the modeling of the DT M sig distribution is studied with alternative models for signal and background. The uncertainties are estimated by comparing with the fit results obtained using the signal and background shapes directly from the simulated samples.
for the nominal amplitude fit are shown from data samples at Ecm between 4.178 and 4.226 GeV. The black points with error bars are data, the red histograms are the results of the nominal amplitude fit, the green shaded histograms are the scaled GMC combinatorial background. For the identical pions, the one giving a lower K 0 S π + invariant mass is denoted as π + 1 , the other is denoted as π + 2 .
• Fit to the ST M tag distribution. We change the background shape from the second-order Chebychev polynomial to a third-order Chebychev polynomial, causing a 0.18% relative change of the branching fraction. The systematic uncertainty due to the modeling of the signal distribution is determined to be 0.16% by performing an alternative fit using the shape directly obtained from the simulated sample. The quadratic sum of these terms, 0.24%, is assigned as the systematic uncertainty.
• Measurement method. The possible bias due to the measurement method is estimated to be 0.3% by comparing the measured branching fraction in the SMC, using the same method as in data analysis, to the value input in the SMC generation.
• Statistics of simulated events. The uncertainty associated with the limited statistics of GMC for the detection efficiency is 0.3%.
• Amplitude analysis model. The uncertainty from the amplitude analysis model is 0.6%, estimated from the efficiency difference obtained by varying the fitted parameters c n in Eq. 3 according to the error matrix.
All the systematic uncertainties of the branching fraction measurement are listed in Table V. When added in quadrature they sum to a relative uncertainty of 3.3%, which is the same as the statistical uncertainty on the measurement.