Amplitude analysis and branching fraction measurement of $D_{s}^{+} \rightarrow K^{+}K^{-}\pi^{+}$

We report an amplitude analysis and branching fraction measurement of $D_{s}^{+} \rightarrow K^{+}K^{-}\pi^{+}$ decay using a data sample of 3.19 $\rm fb^{-1}$ recorded with BESIII detector at a center-of-mass energy of 4.178 GeV. We perform a model-independent partial wave analysis in the low $K^{+}K^{-}$ mass region to determine the $K^{+}K^{-}$ S-wave lineshape, followed by an amplitude analysis of our very pure high-statistics sample. The amplitude analysis provides an accurate determination of the detection efficiency allowing us to measure the branching fraction ${\cal B}(D_{s}^{+} \rightarrow K^{+}K^{-}\pi^{+}) = (5.47\pm0.08_{{\rm stat}}\pm0.13_{{\rm sys}})\%$.

a Also at Bogazici University, 34342 Istanbul, Turkey b Also at the Moscow Institute of Physics and Technology, Moscow 141700, Russia c Also at the Novosibirsk State University, Novosibirsk, 630090, Russia d Also at the NRC "Kurchatov Institute", PNPI, 188300, Gatchina, Russia e Also at Istanbul Arel University, 34295 Istanbul, Turkey f Also at Goethe University Frankfurt, 60323 Frankfurt am Main, Germany The decay D + s → K + K − π + is widely used as a reference mode in D ± s analyses because of its large branching fraction (BF) and low background contamination.An amplitude analysis can reveal the intermediate states involved in this decay and thereby reduce the detection efficiency systematic uncertainties.The improved precision of the BF is important for D ± s analysis using this decay as a reference channel.Furthermore, theoretical studies [1] predict the BFs of D + s → K * (892) 0 K + and D + s → φ(1020)π + to be in the range of (3.9 − 4.2)% and (3.4 − 4.51)%, respectively.Combining the results of the amplitude analysis and the BF measurement, one can obtain the BFs of such intermediate processes, which can help to improve the theoretical model [1].
Dalitz plot analyses of the D + s → K + K − π + decay have been performed by the E687 [2], CLEO [3] and BaBar [4] collaborations.The E687 collaboration used about 700 pure signal events and did not take the f 0 (1370)π + intermediate state into account.In the CLEO analysis about 14400 events with a purity of 84.9% were selected in an untagged analysis of 0.586 fb −1 of data similar to the present analysis.The analysis of BaBar collaboration used about 100 000 events with a purity of about 95%.Table I shows the comparison of the fit fractions (FFs) from these previous Dalitz plot analyses.There are obvious differences between FFs of BaBar collaboration and CLEO collaboration.
The decay D + s → a 0 (980) 0 π + has been observed through D + s → π + π 0 η [5], and should also be present in D + s → K + K − π + , which was not taken into account before.Due to the large overlap of a 0 (980) → K + K − and f 0 (980) → K + K − and their common J P C , we do not distinguish between them in this paper and denote the combined state as S(980).A model-independent partial wave analysis (MIPWA) is performed to study this low-mass resonance.
In this paper we report an amplitude analysis and BF measurement of D + s → K + K − π + (the inclusion of charge conjugates is implied) using a 3.19 fb −1 data sample collected with the BESIII detector at a center-ofmass energy (E CMS ) of 4.178 GeV.At this energy, the cross section for the D * ± s D ∓ s final state in e + e − annihilations is one order magnitude larger than that for D + s D − s [6].Moreover, the D * ± s decays are dominated by the process D * ± s → γD ± s [7].Thus, the process e + e − → D * ± s D ∓ s → D + s γD − s is the main signal process.
Using a tagging technique [8] (described in Sec.V A), we get a nearly background-free data sample to use for an amplitude analysis and BF measurement.The process e + e − → D + s D − s also contributes to the BF measurement.For the MIPWA (Sec.IV), only the signal decay is reconstructed, while for the amplitude analysis (Sec.V) and BF measurement (Sec.VI) both the signal D s and the other D s are reconstructed.

II. BESIII DETECTOR AND DATA SETS
The BESIII detector is a magnetic spectrometer [9] located at the Beijing Electron Positron Collider (BEPCII) [10].The inner subdetectors are surrounded by a superconducting solenoidal magnet which provide a 1.0 T magnetic field.Starting from the interaction point these consist of a main drift chamber (MDC), a plastic scintillator time-of-flight (TOF) system, a CsI(Tl) electromagnetic calorimeter (EMC).Charged particle identification is performed by combining the ionization energy loss (dE/dx) measured by the MDC and the time-of-flight measured by the TOF.The EMC provides shower information to reconstruct photons.Outside the solenoidal magnet is a multi-gap resistive-plate chamber system, which provides muon identification.
Monte Carlo (MC) samples are produced with geant4-based [11] software.To assess background processes and determine detection efficiencies, we produce and analyzed an inclusive MC sample, with size that is 40 times the integrated luminosity of data.The sample includes all known open charm production processes, the continuum processes ( e + e − → q q, q = u, d and s), Bhabha scattering, µ + µ − , τ + τ − , diphoton process and production of the cc resonances J/ψ, ψ(3686) and ψ(3770) via initial state radiation (ISR).The generator conexc [12] is used to model the open charm processes directly produced via e + e − annihilation.The simulation of ISR production of ψ(3770), ψ(3686) and J/ψ is performed with kkmc [13].The known decays with BFs taken from the Particle Data Group (PDG) [7] are simulated with evtgen [14] and the unknown decays are generated with the lundcharm model [15].Finalstate radiation from charged tracks is produced by photos [16].Additionally, we generate two MC samples with e + e − → D   larger than the expected number of signal events in data.
The one with a uniform distribution of D + s → K + K − π + decays over the phase space (PHSP) is called "PHSP MC".In the second sample, called "signal MC", the D + s → K + K − π + decay is generated according to the model obtained from the amplitude analysis presented in this paper.PHSP MC is used to calculate the MC integrations, while signal MC is used to check the fit performance, calculate the goodness of fit and estimate the detection efficiency.

III. EVENT SELECTION
The polar angles (θ) of charged tracks with respect to the beam axis must satisfy |cos θ| < 0.93.Except for tracks from K 0 S decays, the distances of closest approach to the beamspot for charged tracks in the transverse plane and along the beam direction must be less than 1 cm and 10 cm, respectively.
Photons are reconstructed from showers in the EMC.The deposited energies of the photons from the endcap (0.86 < |cos θ| < 0.92) should be larger than 50 MeV and those of the photons from the barrel (|cos θ| < 0.80) should be larger than 25 MeV.Furthermore, the shower should be detected within 700 ns after a beam crossing.
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.
Any π ± and π 0 candidates with momentum less than 0.1 GeV/c are vetoed to remove soft π ± and π 0 from D * decays.
Pairs of π + π − are used to reconstruct K 0 S mesons.The polar angles θ of the two pions should satisfy |cos θ| < 0.93.The distances of closest approach to the beamspot along the beam direction should be less than 20 cm.The invariant mass m(π + π − ) of π + π − pairs should satisfy 0.487 < m(π + π − ) < 0.511 GeV/c 2 .A secondary vertex fit, constraining the pion candidate pair to a common vertex is performed to determine the decay length L of the K S .We require L/σ L > 2, where σ L is the uncertainty on L.
Tagged D s candidates are reconstructed from various combinations of K ± , π ± , η, η , K 0 S and π 0 , while the signal D + s candidates are reconstructed from K + K − π + combinations.Candidates with an invariant mass in the mass window [1.87, 2.06] GeV/c 2 and a recoiling mass M rec in the mass window [2.051, 2.180] GeV/c 2 are retained.The recoiling mass M rec is defined as: where p Ds is the momentum of D s candidate in e + e − center-of-mass system, m Ds is D s mass quoted from PDG [7].The requirement on M rec is chosen to retain both the monochromatic D s that are produced directly from the e + e − collision as well as the broader distribution that arises from D * ± s → D ± s γ decays.
A MIPWA is performed to determine the S-wave lineshape near the threshold of K + K − mass spectrum.As the background contamination is rather low in this region of the Dalitz plot, and higher statistics are needed in this MIPWA, the event selection in this section is different from those in the amplitude analysis (Sec.V) and BF measurement (Sec.VI).
In the data sample used in the MIPWA, D + s → K + K − π + candidates are reconstructed according to the selections in Sec.III.The daughter tracks are further subjected to a 1C kinematic fit constraining them to the nominal D + s mass from PDG [7]; selection of the best D + s → K + K − π + candidate is based on the smallest χ 2 in cases of multiple candidates.The best photon candidate for the decay D * ± s → D ± s γ, is obtained via the recoiling mass against the signal D s and the photon: , (2) where E γ and p γ refer to the energy and momentum of a certain photon candidate in e + e − center-of-mass system, respectively.The photon candidate resulting in the M oth closest to the nominal D s mass is chosen as the best one.
A multi-variate analysis (MVA) method is used to suppress background from the q q continuum and other open charm processes.With the gradient boosted decision tree classifier (BDTG) provided by TMVA [17], we train the MVA separately with two sets of variables for the two categories depending on the D + s origin.Two categories of events are selected in an M rec versus ∆M 2D plane, where ∆M ≡ M (D + s γ) − m sig , m sig is the invariant mass of signal D s and M (D + s γ) refers to the invariant mass of D + s and the photon from D * + s → D + s γ, as shown in Fig. 1.The events that satisfy |M rec − 2.112| < 0.02 GeV/c 2 (the region within the red solid lines in Fig. 1) are denoted as category 1, while the events that satisfy |M rec − 2.112| > 0.02 GeV/c 2 and 0.112 < ∆M < 0.167 GeV/c 2 (the region within the green dashed lines in Fig. 1) are denoted as category 2.
For category 1, the BDTG takes three discriminating variables as input: the recoiling mass M rec , the total momentum of the un-reconstructed objects in the event (not part of the D + s → K + K − π + candidate) and the energy of the photon from D * s .For category 2, the BDTG takes three additional variables as input: ∆M , M rec and the total number of charged tracks and neutrals in an event N tracks .Here, M rec is defined as , where p Dsγ is the momentum of the D s γ combination in e + e − center-of-mass system and m D * s is the nominal D * s mass.According to studies with the inclusive MC sample, the BDTG requirement gives a relatively pure sample (background less than 4%) and the background ratios of category 1 and category 2 are similar.After applying the BDTG requirement, we fit to the candidate signal D s invariant mass for both category 1 and category 2 events.The signal shape is modeled with the MC-simulated shape convolved with a double Gaussian function to account for the difference between data and MC simulation, while the background is described with a second-order Chebychev polynomial.This fit gives a background yield in signal region (1.950 < m sig < 1.986 GeV/c 2 ) of 766 ± 30 and a corresponding signal yield of 18600 ± 141, as shown in Fig. 2.
Assuming N is the number of events for a given mass interval of m(K + K − ), the angular distribution dN d cos Θ can be expanded in terms of spherical harmonic functions: where L max = 2 max , and max is the maximum orbital angular momentum quantum number required to describe the K + K − system at m(K + K − ) (e.g.max =1 when only S-, P-wave are considered), Θ is the angle between the K + direction in the K + K − rest frame and the prior direction of the The background contribution is subtracted from the Here, msig is the mass without 1C kinematic fit correction.The signal shape is the MC-simulated shape convolved with a double Gaussian function and the background shape (red dotted line) is second-order Chebychev polynomial.
selected sample using the shape of the m(K − π + ) versus m(K + K − ) distribution from the inclusive MC sample, while the background normalization is fixed according to the fit results (see Fig. 2).After that the distribution dN/d cos Θ of data are corrected for efficiency and phase space.The distribution m(K − π + ) versus m(K + K − ) of PHSP MC is used to calculate the efficiency.The correction for the Lorentz invariant phase space factor is calculated as , where m K is the nominal mass of K ± [7].The harmonic functions Y 0 k (cos Θ) are normalized as follows: Considering the orthogonality condition, we can obtain the expansion coefficients according to Eqs. ( 3) and ( 4): In this section, the formalism GeV/c 2 ), as shown in Fig. 3.
Assuming that only S-and P-wave amplitudes are necessary at the low end of K + K − mass spectrum, the distribution dN/d cos Θ can also be written in terms of the partial wave amplitudes: where S and P refer to the amplitudes of S-wave and Pwave, respectively.Comparing Eqs. ( 3) and ( 6) [18], we can obtain where φ SP = φ S − φ P is the phase difference between Swave and P-wave, φ S and φ P are the phases of S-wave and P-wave, respectively.Calculating |S| 2 , φ SP and |P| 2 in every mass interval of m(K + K − ) in the threshold region, the distribution of |S| 2 , |P| 2 , φ SP and φ S can be obtained, as shown in Fig. 4.There are two curves in Fig. 4(c) because of the sign ambiguity of φ SP extracted from cos φ SP .
We found that the Flatté parameterization [19] is insensitive to the ππ or πη coupling or the coupling induced between them while fitting the distribution of |S| 2 .Therefore, the lineshape of S(980) is empirically parameterized with the following formula: Fitting the distribution of |S| 2 in Fig. 4(a) with |A S(980) | 2 , we can obtain the values of m 0 and Γ 0 : Figure 4(a) shows the fit result.The χ 2 /NDF of the fit is 44.46/38 = 1.17.
The S(980) mass central value obtained from the fit is much lower than the threshold of m(K + K − ) (about 0.988 GeV/c 2 ).Therefore the distribution of φ S is expected to be roughly constant.The phase φ P of the φ(1020) is given by Eq. ( 21) in Sec.V B; it increases rapidly near the φ(1020) peak because of its narrow width.Then the sign ambiguity of φ SP is solved by choosing the black curve in Fig. 4(c), which decreases rapidly near the mass of the φ(1020), ensuring that φ S = φ P +φ SP is roughly constant.The resulting phase of the S(980), φ S , is shown in Fig. 4(d).The solid line in Fig. 4(d) shows the phase of S(980) amplitude obtained from the amplitude analysis (described in Sec.V).We can see that the shapes of φ S from model-independent analysis     2 (arbitrary units) and φ S in every mass interval of the threshold region are listed in Table II.
Systematic uncertainties considered for the MIPWA include: • Data-MC agreement for the BDTG output.A control sample is obtained with same event selection as that in Sec.V A due to its high purity, but without the kinematic fit criteria.The efficiency of data and MC samples from the BDTG requirement is then considered, where the efficiency of data (MC) is defined as , where N d0 (N M 0 ) and N d1 (N M 1 ) are the number of events before and after applying the BDTG requirement.We can now correct the data sample with e data eMC .We fit the corrected shape of the S(980) and take the shift of m 0 and Γ 0 as the systematic uncertainty.
• Background subtraction.We change the bin size, fit range and replace the background shape with a third-order Chebychev polynomial in the fit shown in Fig. 2. New fits are performed and we take the quadrature sum of the shifts as the uncertainty of the background fraction.Then we vary the background fraction, (3.9 ± 0.3)%, within its uncertainty and take the shift of the S(980) fit results as the systematic uncertainty related to the background fraction.The background shape of inclusive MC sample is also replaced with that of sideband (1.90 < m sig < 1.95 GeV/c 2 and 1.986 < m sig < 2.03 GeV/c 2 ) for data to perform a fit and the shift is taken as the systematic uncertainty related to the background shape.The quadrature sum of the shifts of m 0 and Γ 0 are 0.002 GeV/c 2 and 0.001 GeV, respectively.
• Particle identification (PID) and tracking efficiency difference between data and MC simulation.The PID efficiencies are studied using control samples of e + e − → K is used for the study of tracking efficiencies.In these control samples, all final particles are reconstructed with the selection criteria mentioned in Sec.III except the target kaon (pion).The total number of the target kaon (pion) are inferred by fitting the missing mass distributions while the number of the reconstructed target kaon (pion) is determined by applying corresponding selection criteria.The efficiency difference of data and MC samples is assigned as the associated systematic uncertainty.These efficiencies are also used in the amplitude analysis (Sec.V) and BF measurement (Sec.VI C).We weight each event with the data/MC efficiency differences and fit the shape of the S(980).The shift of m 0 and Γ 0 are 0.001 GeV/c 2 and 0.013 GeV, respectively.
• The f 0 (1370) contribution.The f 0 (1370) contribution in the S(980) region was subtracted according to the measured FF.The shape of f 0 (1370) at the low end of m K + K − mass spectrum was obtained from the MC simulation.The interference effect was ignored.The resulting shifts of m 0 and Γ 0 are 0.001 GeV/c 2 and 0.003 GeV, respectively.
All of the systematic uncertainties mentioned above are summarized in Table III.The quadrature sum of the uncertainties is taken as the total uncertainty.We obtain the result for m 0 and Γ 0 with statistical and systematic errors to be: m 0 = (0.919 ± 0.006 stat ± 0.030 sys ) GeV/c 2 , Γ 0 = (0.272 ± 0.040 stat ± 0.024 sys ) GeV, (10) which are consistent with the BaBar analysis [4].Note that m 0 and Γ 0 in Eq. ( 10) are only used for the parameterization of the S(980) in Sec.V.

V. AMPLITUDE ANALYSIS
An unbinned maximum likelihood method is used to determine the intermediate resonance composition in the decay D + s → K + K − π + .The likelihood function is constructed with a probability density function (PDF) which depends on the momenta of the three daughter particles.

A. Tag Technique in Amplitude Analysis
As D s mesons are produced in pairs, D s mesons can be reconstructed with a tag technique which provides both single tag (ST) and double tag (DT) samples.In the ST samples, only one D − s meson is reconstructed through selected hadronic D s decays, the so-called tag modes.The eight tag modes used in the amplitude analysis and BF measurement presented in Sec.VI are  in the corresponding mass intervals are negative and a physical solution for φSP cannot be found according to Eq. (7).four-momentum of the electron-positron system and the invariant mass of the D * ± s to the corresponding PDG value [7].This gives a total of five constraints (5C).For each D + s D − s γ candidate, the extra γ is paired with both the tag and signal D s to form the D * ± s , and the combination with the lower fit χ 2 5C is retained as the presumably correct pairing.If there are multiple candidate D * ± s D ∓ s pairs in an event, the candidate with minimum χ 2 5C is selected as the best one.The invariant mass of signal D s (m sig ) and tag D s (m tag ) candidates are required to be within the mass regions shown in Table IV.

Mode
Mass window (GeV/c 2 ) To ensure that all events fall within the physical region on the Dalitz plot, we perform a 7C fit where constraints on both signal and tag D s masses to the PDG values are added to the previous 5C constraints.The four-momenta of the tracks after 7C fit are used to perform the amplitude analysis.
The background of the DT sample in the amplitude analysis is estimated using the inclusive MC sample.The fit to the signal D s invariant mass without 7C kinematic fit gives the signal yield and purity, as shown in Fig. 5.In the fit, the signal shape is modeled with the MC-simulated shape convolved with a Gaussian function while the background is described with a second-order Chebychev polynomial.There is no obvious peaking background in the signal region (1.950 < m sig <1.986 GeV/c 2 ) and we obtain 4399 signal candidates with a purity of 99.6%.Figure 6 shows the Dalitz plot of the signal D + s → K + K − π + candidates.

B. Likelihood Function Construction
For a three-body process the amplitude A n (p) for the n th mode may be written as where p refers to the set of the three daughter particles' four-momenta, P n (p) is the propagator, S n (p) is the spin factor constructed with the covariant tensor formalism [20] decays, respectively.According to the isobar formulation, the total amplitude M (p, α) is obtained by the coherent sum of all intermediate modes: where α is a set of fit parameters, which includes the n th mode of complex coefficient c n = ρ n e iφn (ρ n and φ n are the magnitude and phase, respectively).Then the signal PDF f S (p, α) is given by: where (p) is the detection efficiency and R 3 (p) is the three-particle phase space density, which is defined as: where β = 1, 2, 3 is the index of the three daughter particles.The likelihood function is given by where N data is the number of candidate events in data.Note that the last term of Eq. ( 15) is independent of the fit parameters and dropped during the log-likelihood fit.
Background contribution is neglected in the amplitude analysis and the possible bias is included to the systematic uncertainties, see Sec.V D below.The normalization integral in Eq. ( 15) is first determined by the following equation using PHSP MC events with a uniform distribution: where p kMC is the k th MC set of four-momenta.Here, N MC,gen and N MC,sel are the numbers of generated phasespace events and selected phase-space events, respectively.A set of estimated fit parameters, denoted as α , is obtained from a preliminary fit using the phase-space MC to evaluate the normalization integral.The normalization integral is evaluated with signal MC samples: where M gen (p kMC , α ) is the PDF modeled with α to generate signal MC and N MC is the number of events in the MC sample.The computational efficiency of the MC integration is significantly improved by evaluating the normalization integral with signal MC samples, which intrinsically take into account the event selection acceptance and the detection resolution.Correction factors γ are introduced to correct for the bias caused by PID and tracking efficiency inconsistencies between data and MC simulation: where j refers to PID or tracking, j,data and j,M C refer to the PID or tracking efficiencies for data and MC, respectively.Taking the correction factors γ into account, the normalization integral can be obtained by: 1. Propagator and Blatt-Weisskopf Barrier For a given two-body decay (a → bc), p a , p b and p c are the momenta of particles a, b and c.The variables s a , s b and s c refer to the squared invariant masses of particles a, b and c.The momentum q is defined as the magnitude of the momentum of b or c in the rest system of a: The resonances K * (892), f 0 (1710), φ(1020) and f 0 (1370) are parameterized with a relativistic Breit-Wigner formula, where m 0 and Γ 0 are the mass and the width of the intermediate resonance, fixed to the PDG values [7], with the exception of f 0 (1370).The mass and width of f 0 (1370) are fixed to 1350 MeV/c 2 and 265 MeV [21], respectively.The value of q 0 in Eq. ( 21) is that of q when s a = m 2 0 , L denotes the angular momenta and Blatt-Weisskopf Barrier F L (q) is defined as: where z = qR and z 0 = q 0 R. R is the effective radius of the intermediate resonance or D s meson.The values of R are fixed to 3.0 GeV −1 for intermediate states and 5.0 GeV −1 for D s meson, respectively.The uncertainty of R values is taken into account in evaluation of the systematic uncertainties.
The K * 0 (1430) 0 is parameterized with the Flatté formula: where s is the squared K − π + invariant mass, ρ Kπ (s) and ρ η K (s) are Lorentz invariant PHSP factor, and g 1,2 are coupling constants to the corresponding final state.The parameters of the K * 0 (1430) 0 are fixed to values measured by CLEO [22].
For the resonance S(980) (representing the f 0 (980) and a 0 (980)), we use Eq. ( 8) to describe the propagator and the values of parameters are fixed to those in Eq. ( 9) obtained from the MIPWA section (Sec.IV).

Spin Factors
The spin projection operators [20] for a two-body decay are defined as µνµ ν (a) = ( The corresponding covariant tensors are expressed as follows where r a = p b −p c is the momentum difference between b and c.The spin factor for the process D s → aX (where a is a resonance and X is a direct daughter of the D s meson) with a → bc is, where T (L) µ (D s ) and t(L) µ (a) are the covariant tensors with angular momenta L for D s → aX and a → bc, respectively.

C. Fit Result
We start the fit of data by considering the amplitudes containing K * (892) 0 , φ(1020) and S(980) resonances, as these resonances are clearly seen in Fig. 6.We choose K * (892) 0 as the reference amplitude and fix the magnitude ρ and phase φ for D + s → K * (892) 0 K + to 1.0 and 0.0, respectively.The magnitudes and phases of other processes are free parameters in the fit.We then add amplitudes with resonances listed in the PDG [7] and non-resonant components until no additional amplitude has significance larger than 5σ.The statistical significance for a certain intermediate process is calculated using the change of likelihood and number of degrees of freedom between with and without this process.The six intermediate processes retained in the nominal fit are D + s → K * (892) 0 K + , φ(1020)π + , S(980)π + , K * 0 (1430) 0 K + , f 0 (1370)π + and f 0 (1710)π + .The magnitudes, phases and corresponding significances of these amplitudes are listed in Table V.Other tested amplitudes when determining the nominal fit model, but finally not used, are listed in Table VI.The systematic uncertainties are discussed in Sec.V D.
With the coefficients c n obtained from the fit, the FFs are calculated with generator-level phase-space MC as where the summation is performed over the generated PHSP MC events.
To properly treat correlations, we randomly vary the coefficients c n according to the corresponding error matrix to produce many sets of c n and then obtain a series of FFs for each intermediate process.A Gaussian function is used to fit the distribution of FF for each intermediate process and the width of the Gaussian function is taken as the corresponding statistical uncertainty of the FF.The resultant FFs are listed in Table V.
Signal MC samples modeled according to the fit result are generated to compare the projections of the Dalitz plots with data and to calculate the fit bias, which will be discussed in Sec.V D. The Dalitz plot projections are shown in Fig. 7. To evaluate the goodness of fit with a of the fit using an adaptive binning of the Dalitz plot of m 2 (K + K − ) versus m 2 (K − π + ), in which each bin has at least 10 events.Here N data , σ data and N exp refer to the number of events from data, the error of N data and the expected number obtained from signal MC in each bin, respectively.We find a χ 2 /NDF = 290.0/280.

D. Systematic Uncertainty
The following categories of systematic uncertainties are studied for the amplitude analysis: I Resonance parameters.The masses and widths of resonances are shifted by their corresponding uncertainties.For the S(980), m 0 and Γ 0 are shifted according to the errors from Eq. (10).The mass and width of f 0 (1370) are shifted according to the uncertainties from Ref. [21].The parameters of K * 0 (1430) 0 are shifted according to the errors from  Ref. [22].For other states, uncertainties are taken from the PDG [7].
II The effective radius in the Blatt-Weisskopf barrier factor is varied within the range [1.0, 5.0] GeV −1 for intermediate resonances and [3.0, 7.0] GeV −1 for D s mesons.
III Fit bias.Pull distribution checks using 300 signal MC samples are performed to obtain the fit bias.
Here the pull value for a certain parameter x is defined as (x true − x MC )/σ xMC , where x MC and σ xMC are the value and the statistical error of x obtained from the fit to a certain signal MC sample and x true refers to the true value of x used in the MC generation.The signal MC samples each have the same size as the data.Fits to the pull distributions with Gaussian functions show no obvious biases and under-or over-estimations on statistical uncertainties.We add quadrature sum of the mean value and the error of mean to get the corresponding systematic uncertainty in units of the corresponding statistical uncertainty.
IV Detector effects.These effects are related to the efficiency difference between MC simulation and data caused by PID and tracking, reflected in the γ in Eq. ( 18).The uncertainties associated with γ are obtained by performing alternative amplitude analyses varying PID and tracking efficiencies according to their uncertainties.
V Model assumptions.We replace the Flatté expression in Eq. ( 23) with the LASS model [23].For the S(980), Eq. ( 8) is replaced with the Flatté parameterization [19] to describe the lineshape of the S(980) and the parameters in the Flatté parameterization are obtained from the fit to |S| 2 in Fig. 4(a).The quadrature sum of the shifts in the results are taken as the corresponding systematic uncertainties.
VI Background estimation.The background is ignored in the nominal fit.We subtract the contribution of the background by assigning a negative weight to the background events in the likelihood calculation [24].Individual changes of the results with respect to the nominal one are taken as the corresponding systematic uncertainties.
VII Contributions with statistical significances less than 5σ.The intermediate processes with statistical significances less than 5σ are added in the nominal fit one by one.The quadrature sum of each parameter variations is taken as the corresponding systematic uncertainty.
Systematic uncertainties on the magnitudes, phases and FFs are summarized in Table VII and the total uncertainties are obtained as the sum of all the contributions in quadrature.

A. Efficiency and Data Yields
After the selection described in Sec.III, the tag technique is also used to perform the BF measurement.We use the same eight tag modes as in Sec.V.For each tag mode, if there are multiple tag D s candidates in an event, the candidate with M rec closest to the nominal mass of D * s [7] is retained.The ST yields are obtained by the fits to the D s invariant mass distributions, as shown in Fig. 8, along with the mass windows listed in Table IV.The signal shape is modeled as the MC-simulated shape convolved with a Gaussian function, while background is parameterized as a second-order Chebychev polynomial.Fits to m tag for inclusive MC are performed to estimate the corresponding ST efficiencies.The ST yields (Y ST ) and ST efficiencies ( ST ) are listed in Table VIII.
After the best candidates of ST D − s mesons are identified, we search for the D + s → K + K − π + .Only the best D + s candidate with the average mass of tag D − s and signal D + s closest to the nominal mass of D s is retained for each tag mode in an event.For all tag modes, we found a mean of 0.16% (0.17%) of DT events in data (inclusive MC) contains multiple D + s candidates.The effect due to the multiple candidate selection is, therefore, negligible.The DT efficiencies, listed in Table IX, are obtained based on the signal MC samples.As D − s → K + K − π − is not only the signal mode but also one of the tag modes, we divide the events into two categories: -Cat.A: Tag D − s decays to one of the tag modes except D − s → K + K − π − .The inclusive MC sample with the signal removed shows no peaking background around the fit range of 1.90 < m sig < 2.03 GeV/c 2 .Thus, the DT yield is determined by the fit to m sig , shown in Fig. 9(a).The background is described with a second-order Chebychev polynomial.The DT yield is 3497 ± 64.
-Cat.B: Tag D − s decays to K + K − π − .As both of the two D s mesons decay to the signal modes, we fit dM (the mass of the signal D + s minus that of the tag D − s ), which is shown in Fig. 9(b).Here, the background is described by a second-order Chebychev polynomial.The DT yield is 1651 ± 42.

B. Tagging Technique and Branching Fraction
For the DT samples with a certain tag mode α, we have where N D + s D − s is the total number of D * ± s D ∓ s produced from e + e − collision; the yields N obsA,α sig and N obsB,α sig refer to the yields with tag mode α for Cat.A and Cat.B, respectively; B tag and B sig are the BFs of a specific tag mode and the signal mode, respectively; tag is the efficiency to reconstruct the tag mode; tag,sig is the efficiency to reconstruct both the tag and signal decay modes.
Using the above equations, one can obtain: where the yields N obsA sig , N obsB sig and Y α ST are obtained from data, while tag and tag,sig can be obtained from the updated inclusive MC samples.The process D + s → K + K − π + in the updated inclusive MC is generated with the Dalitz model obtained in Sec.V.

C. Systematic Uncertainty
Most systematic uncertainties related to the efficiency for reconstructing the tag side cancel for BF measurement due to the DT technique.The following sources are taken in account to calculate systematic uncertainty.
• Uncertainty in the number of ST D − s candidates.We perform alternative fits with different background shapes and signal shapes to obtain these uncertainties.We change the background shape from a second-order Chebychev polynomial to a third-order Chebychev polynomial and the relative change of BF is 0.18%.The systematic uncertainty in signal shape is determined to be 0.16% by performing an alternative fit without convolution with the Gaussian smearing function.The quadrature sum of these terms, that is the uncertainty in the number of ST D − s candidates, is 0.23%.• DT signal shape.The systematic uncertainty due to the signal shape is studied with the fit without the Gaussian function convolved, the DT yield shift is taken as the related uncertainty.
• DT background shape.For background shape in the fit, a third-order Chebychev polynomial is used to replace the nominal one.The quadrature sum of the BF shifts is taken as the related uncertainty.
• Fit bias.The updated inclusive 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.1%, which is negligible.
• K ± and π ± tracking/PID efficiency.The ratios between data and MC efficiencies are weighted by the corresponding momentum spectra of signal MC events.We obtain the systematic uncertainties related to tracking efficiency to be 0.5% for each kaon track and 0.2% for each pion track based on the study of the tracking efficiency.The systematic uncertainties related to PID efficiencies are estimated to be 0.5% for each K ± and 0.4% for each π ± .Tracking efficiency systematics are added linearly for the three tracks, as are the PID efficiency systematics.
• MC statistics.The uncertainty due to the MC statistics is obtained as , where f α is the DT yield fraction, α is the DT signal efficiency of the tag mode α and δ α is the error on α due to the limited MC statistics.
• Dalitz model.The uncertainty from the Dalitz model is estimated as the change of efficiency when the Dalitz model parameters (c n ) are varied according to the error matrix.
All of the systematic uncertainties mentioned above are summarized in Table X.We take the quadrature sum of the systematic uncertainties above as the total systematic uncertainty in the BF of D + s → K + K − π + .This paper presents the amplitude analysis of the decay D + s → K + K − π + .The results on FFs for D + s → f 0 (1370)π + , D + s → f 0 (1710)π + and D + s → f 0 (980)π + /a 0 (980)π + are consistent with those of BaBar and E687.In addition, our results on FFs also agree with those of CLEO, except for D + s → f 0 (980)π + /a 0 (980)π + and D + s → f 0 (1370)π + where 2.4σ and 3.4σ differences, respectively, with CLEO are observed.
In this analysis, as a 0 (980) and f 0 (980) overlap and parameters of a 0 (980) and f 0 (980) are not well measured, we have extracted the S-wave lineshape in the low end of K + K − mass spectrum with a model-independent method.
We have also measured the BF B(D + s → K + K − π + ) = (5.47 ± 0.08 stat ± 0.13 sys )% which is currently the most precise measurement.Comparisons with other results are presented in Tables XI and Tables XII.

g
Also at Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education; Shanghai Key Laboratory for Particle Physics and Cosmology; Institute of Nuclear and Particle Physics, Shanghai 200240, People's Republic of China h Also at Key Laboratory of Nuclear Physics and Ion-beam Application (MOE) and Institute of Modern Physics, Fudan University, Shanghai 200443, People's Republic of China i Also at Harvard University, Department of Physics, Cambridge, MA, 02138, USA j Currently at: Institute of Physics and Technology, Peace Ave.54B, Ulaanbaatar 13330, Mongolia k Also at State Key Laboratory of Nuclear Physics and Technology, Peking University, Beijing 100871, People's Republic of China l School of Physics and Electronics, Hunan University, Changsha 410082, China (Dated: August 6, 2021) I. INTRODUCTION ( * ) s D s , in which the D + s meson decays into K + K − π + while the D − s meson decays to one of the tag modes listed in

FIG. 1 .
FIG. 1. Two dimensional plane of Mrec versus ∆M ≡ M (D + s γ) − msig from the simulated D + s → K + K − π + decays.The red solid (green dashed) lines mark the mass window for the D + s category 1 (category 2) around the Mrec (∆M ) peak.

FIG. 2 .
FIG.2.The fit to the signal Ds invariant mass msig (the dots with error bars) after BDTG requirement.The area between the pink lines is the signal area of the sample for MIPWA.Here, msig is the mass without 1C kinematic fit correction.The signal shape is the MC-simulated shape convolved with a double Gaussian function and the background shape (red dotted line) is second-order Chebychev polynomial.

FIG. 4 . 2 (
FIG. 4. The distribution of (a) |S| 2 , (b) |P| 2 , (c) φSP and (d) φS in the threshold region of m(K + K − ).The description of the fit in (a) can be found in the text.The solid line in (d) shows the phase of S(980) amplitude obtained from the amplitude analysis (Sec.V).
In the DT samples, photons from the decay D * ± s → D ± s γ, tag mode D − s and signal D + s (i.e., decays to K + K − π + ) are all fully reconstructed.A 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.We constrain the four-momentum of the D * ± s D ∓ s system to the initial

FIG. 5 .
FIG.5.The fit to the signal Ds invariant mass msig before the 7C kinematic fit (dots with error bars).The area between the pink lines is the signal area of the sample for the amplitude analysis.

FIG. 8 .
FIG. 8. Fits to the mtag distributions of data.The points with error bars indicate data and the solid lines indicate the fit.Red short-dashed lines are signal, violet long-dashed lines are background.The region within the purple lines denotes the signal region.

FIG. 9 .
FIG. 9. Fit of msig for (a) Cat.A and dM for (b) Cat.B. The signal shapes are the corresponding simulated shapes convolved with a Gaussian function and the background shapes are described with second-order Chebychev polynomials.
Table IV, with size 600 times

TABLE I .
Comparison of FFs for different decay modes.The first and second uncertainties are statistical and systematic, respectively.

TABLE II .
The values of |S| 2 (arbitrary units), |P| 2 (arbitrary units) and φS; the units chosen preserve the relative |S| 2 and |P| 2 sizes.Uncertainties in the table are statistical only.Some values of φS are not listed in the table because the values of Y 0 2

TABLE III .
|S| 2 (arbitrary units) |P| 2 (arbitrary units) φS (degrees) Systematic uncertainties of the partial wave analysis in the low K + K − mass region.The quadrature sum of all contributions is taken as the total uncertainty.

TABLE IV .
The mass windows for the signal mode and various tag modes.

TABLE V .
The results on the magnitudes, phases, FFs and significances for the six amplitudes.The first and second uncertainties are the statistical and systematic, respectively.

TABLE VI .
The significances for other tested amplitudes.

TABLE VII .
Systematic uncertainties on the φ, ρ and FFs for different amplitudes in units of the corresponding statistical uncertainties.Here I, II, III, IV, V, VI and VII denote the propagator parameterizations of the resonances, the effective radius of Blatt-Weisskopf barrier factor, fit bias, detector effects, model assumptions, background estimation and contributions with statistical significances less than 5σ, respectively.The quadrature sums of these terms are taken as the total systematic uncertainties.

TABLE VIII .
The ST yields (YST) and ST efficiencies ( ST).The BFs of the sub-particle (K 0 S , π 0 , η and η ) decays are not included in the efficiencies.

TABLE IX .
The DT efficiencies ( DT).The BFs of the subparticle (K 0 S , π 0 , η and η ) decays are not included in the efficiencies.

TABLE X .
The relative systematic uncertainties on the BF.The quadrature sum of all contributions is taken as the total uncertainty.