Amplitude analysis of D + → K

M. Ablikim, M. N. Achasov, S. Ahmed, M. Albrecht, M. Alekseev, A. Amoroso, F. F. An, Q. An, Y. Bai, O. Bakina, R. Baldini Ferroli, Y. Ban, K. Begzsuren, J. V. Bennett, N. Berger, M. Bertani, D. Bettoni, F. Bianchi, J. Bloms, I. Boyko, R. A. Briere, H. Cai, X. Cai, A. Calcaterra, G. F. Cao, N. Cao, S. A. Cetin, J. Chai, J. F. Chang, W. L. Chang, G. Chelkov, Chen, G. Chen, H. S. Chen, J. C. Chen, M. L. Chen, S. J. Chen, Y. B. Chen, W. Cheng, G. Cibinetto, F. Cossio, X. F. Cui, H. L. Dai, J. P. Dai, X. C. Dai, A. Dbeyssi, D. Dedovich, Z. Y. Deng, A. Denig, I. Denysenko, M. Destefanis, F. De Mori, Y. Ding, C. Dong, J. Dong, L. Y. Dong, M. Y. Dong, Z. L. Dou, S. X. Du, J. Z. Fan, J. Fang, S. S. Fang, Y. Fang, R. Farinelli, L. Fava, F. Feldbauer, G. Felici, C. Q. Feng, M. Fritsch, C. D. Fu, Y. Fu, Q. Gao, X. L. Gao, Y. Gao, Y. Gao, Y. G. Gao, Z. Gao, B. Garillon, I. Garzia, A. Gilman, K. Goetzen, L. Gong, W. X. Gong, W. Gradl, M. Greco, L. M. Gu, M. H. Gu, Y. T. Gu, A. Q. Guo, L. B. Guo, R. P. Guo, Y. P. Guo, A. Guskov, S. Han, X. Q. Hao, F. A. Harris, K. L. He, F. H. Heinsius, T. Held, Y. K. Heng, Y. R. Hou, Z. L. Hou, H. M. Hu, J. F. Hu, T. Hu, Y. Hu, G. S. Huang, J. S. Huang, X. T. Huang, X. Z. Huang, Z. L. Huang, T. Hussain, N. Hsken, W. Ikegami Andersson, W. Imoehl, M. Irshad, Q. Ji, Q. P. Ji, X. B. Ji, X. L. Ji, H. L. Jiang, X. S. Jiang, X. Y. Jiang, J. B. Jiao, Z. Jiao, D. P. Jin, S. Jin, Y. Jin, T. Johansson, N. Kalantar-Nayestanaki, X. S. Kang, R. Kappert, M. Kavatsyuk, B. C. Ke, I. K. Keshk, T. Khan, A. Khoukaz, P. Kiese, R. Kiuchi, R. Kliemt, L. Koch, O. B. Kolcu, B. Kopf, M. Kuemmel, M. Kuessner, A. Kupsc, M. Kurth, M. G. Kurth, W. Kühn, J. S. Lange, P. Larin, L. Lavezzi, H. Leithoff, T. Lenz, C. Li, Cheng Li, D. M. Li, F. Li, F. Y. Li, G. Li, H. B. Li, H. J. Li, J. C. Li, J. W. Li, Ke Li, L. K. Li, Lei Li, P. L. Li, P. R. Li, Q. Y. Li, W. D. Li, W. G. Li, X. L. Li, X. N. Li, X. Q. Li, X. H. Li, Z. B. Li, H. Liang, H. Liang, Y. F. Liang, Y. T. Liang, G. R. Liao, L. Z. Liao, J. Libby, C. X. Lin, D. X. Lin, Y. J. Lin, B. Liu, B. J. Liu, C. X. Liu, D. Liu, D. Y. Liu, F. H. Liu, Fang Liu, Feng Liu, H. B. Liu, H. M. Liu, Huanhuan Liu, Huihui Liu, J. B. Liu, J. Y. Liu, K. Y. Liu, Ke Liu, Q. Liu, S. B. Liu, T. Liu, X. Liu, X. Y. Liu, Y. B. Liu, Z. A. Liu, Zhiqing Liu, Y. F. Long, X. C. Lou, H. J. Lu, J. D. Lu, J. G. Lu, Y. Lu , Y. P. Lu, C. L. Luo, M. X. Luo, P.W. Luo, T. Luo, X. L. Luo, S. Lusso, X. R. Lyu, F. C. Ma, H. L. Ma, L. L. Ma, M.M. Ma, Q. M. Ma, X. N. Ma, X. X. Ma, X. Y. Ma, Y. M. Ma, F. E. Maas, M. Maggiora, S. Maldaner, Q. A. Malik, A. Mangoni, Y. J. Mao, Z. P. Mao, S. Marcello, Z. X. Meng, J. G. Messchendorp, G. Mezzadri, J. Min, T. J. Min, R. E. Mitchell, X. H. Mo, Y. J. Mo, C. Morales Morales, N. Yu. Muchnoi, H. Muramatsu, A. Mustafa, S. Nakhoul, Y. Nefedov, F. Nerling, I. B. Nikolaev, Z. Ning, S. Nisar, S. L. Niu, S. L. Olsen, Q. Ouyang, S. Pacetti, Y. Pan, M. Papenbrock, P. Patteri, M. Pelizaeus, H. P. Peng, K. Peters, J. Pettersson, J. L. Ping, R. G. Ping, A. Pitka, R. Poling, V. Prasad, M. Qi, T. Y. Qi, S. Qian, C. F. Qiao, N. Qin, X. P. Qin, X. S. Qin, Z. H. Qin, J. F. Qiu, S. Q. Qu, K. H. Rashid, C. F. Redmer, M. Richter, M. Ripka, A. Rivetti, M. Rolo, G. Rong, Ch. Rosner, M. Rump, A. Sarantsev, M. Savrié, K. Schoenning, W. Shan, X. Y. Shan, M. Shao, C. P. Shen, P. X. Shen, X. Y. Shen, H. Y. Sheng, X. Shi, X. D. Shi, J. J. Song, Q. Q. Song, X. Y. Song, S. Sosio, C. Sowa, S. Spataro, F. F. Sui, G. X. Sun, J. F. Sun, L. Sun, S. S. Sun, X. H. Sun, Y. J. Sun, Y. K. Sun, Y. Z. Sun, Z. J. Sun, Z. T. Sun, Y. T. Tan, C. J. Tang, G. Y. Tang, X. Tang, V. Thoren, B. Tsednee, I. Uman, B. Wang, B. L. Wang, C.W. Wang, D. Y. Wang, H. H. Wang, K. Wang, L. L. Wang, L. S. Wang, M. Wang, M. Z. Wang, Meng Wang, P. L. Wang, R. M. Wang, W. P. Wang, X. Wang, X. F. Wang, Y. Wang, Y. F. Wang, Z. Wang, Z. G. Wang, Z. Y. Wang, Zongyuan Wang, T. Weber, D. H. Wei, P. Weidenkaff, H.W. Wen, S. P. Wen, U. Wiedner, M. Wolke, L. H. Wu, L. J. Wu, Z. Wu, L. Xia, Y. Xia, S. Y. Xiao, Y. J. Xiao, Z. J. Xiao, Y. G. Xie, Y. H. Xie, T. Y. Xing, X. A. Xiong, Q. L. Xiu, G. F. Xu, L. Xu, Q. J. Xu, W. Xu, X. P. Xu, F. Yan, L. Yan, W. B. Yan, W. C. Yan, Y. H. Yan, H. J. Yang, H. X. Yang, L. Yang, R. X. Yang, S. L. Yang, Y. H. Yang, Y. X. Yang, Yifan Yang, Z. Q. Yang, M. Ye, M. H. Ye, J. H. Yin, Z. Y. You, B. X. Yu, C. X. Yu, J. S. Yu, C. Z. Yuan, X. Q. Yuan, Y. Yuan, A. Yuncu, A. A. Zafar, Y. Zeng, B. X. Zhang, B. Y. Zhang, C. C. Zhang, D. H. Zhang, H. H. Zhang, H. Y. Zhang, J. Zhang, J. L. Zhang, J. Q. Zhang, J. W. Zhang, J. Y. Zhang, J. Z. Zhang, K. Zhang, L. Zhang, S. F. Zhang, T. J. Zhang, X. Y. Zhang, Y. Zhang, Y. H. Zhang, Y. T. Zhang, Yang Zhang, Yao Zhang, Yu Zhang, Z. H. Zhang, Z. P. Zhang, Z. Y. Zhang, G. Zhao, J. W. Zhao, J. Y. Zhao, J. Z. Zhao, Lei Zhao, Ling Zhao, M. G. Zhao, Q. Zhao, S. J. Zhao, T. C. Zhao, Y. B. Zhao, Z. G. Zhao, A. Zhemchugov, B. Zheng, J. P. Zheng, Y. Zheng, Y. H. Zheng, B. Zhong, L. Zhou, L. P. Zhou, Q. Zhou, X. Zhou, X. K. Zhou, X. R. Zhou, PHYSICAL REVIEW D 100, 072008 (2019)


I. INTRODUCTION
Hadronic decays of mesons with charm are an important tool for understanding the dynamics of the strong interaction in the low energy regime. The amplitudes describing D meson weak decays into four-body final states are dominated by (quasi-) two-body processes, such as D → VP, D → SP, D → VV, and D → AP, where P, V, S, and A denote pseudoscalar, vector, scalar, and axial-vector mesons, respectively. Final-state interactions can cause significant changes in decay rates and shifts in the phases of decay amplitudes. Experimental measurements can help to refine theoretical models of these phenomena [1][2][3]. Many measurements on D → PP and D → VP decays have been performed [4]. However, there are only a few studies focusing on D → AP decays [4]. We have therefore measured D → AP decays via an amplitude analysis of the decay D þ → K 0 S π þ π þ π − (the inclusion of charge conjugate reaction is implied throughout the paper), which is expected to be dominated by D þ → K 0 S a 1 ð1260Þ þ . In addition, the measurements of the intermediate processes containing K 1 ð1270Þ and K 1 ð1400Þ are helpful for understanding the mixture between these two axial-vector kaons [3].
In this paper, we present an amplitude analysis of the decay D þ → K 0 S π þ π þ π − to study the resonant substructures and nonresonant components, where the a Also amplitude model is constructed using the covariant tensor formalism [5].

II. DETECTOR AND DATA SETS
The data used in this analysis were accumulated with the BESIII detector [6]. The event sample is based on 2.93 fb −1 of e þ e − collisions at the ψð3770Þ mass [7,8]. At this energy, D meson pairs are produced without any additional hadrons. To suppress backgrounds from other charmed meson decays and continuum (QED process and light quark productions), only the decay mode D − → K þ π − π − is used to tag the D þ D − pairs. This provides a clean environment for selecting the decay D þ → K 0 S π þ π þ π − (the signal side) by requiring the D − → K þ π − π − decay to be observed (the tag side).
The BESIII detector located at Beijing Electron Positron Collider [9] is described in Ref. [6]. The geometrical acceptance of the BESIII detector is 93% of the full solid angle. Starting from the interaction point (IP), it consists of a main drift chamber (MDC), a time-of-flight (TOF) system, and a CsI(Tl) electromagnetic calorimeter, 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 momentum resolution for charged tracks in the MDC is 0.5% at a transverse momentum of 1 GeV=c. The energy resolution for the photon in electromagnetic calorimeter measurement is 2.5% (5%) in the barrel (end caps) region at 1 GeV. The time resolution of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps.
Monte Carlo (MC) simulations of the BESIII detector are based on GEANT4 [10]. The production of ψð3770Þ is simulated with the KKMC [11] package, taking into account the beam energy spread and the initial-state radiation (ISR). The PHOTOS [12] package is used to simulate the final-state radiation of charged particles. The EVTGEN [13] package is used to simulate the known decay modes with branching fractions (BFs) taken from the Particle Data Group (PDG) [4], and the remaining unknown decays are generated with the LUNDCHARM model [14]. The MC sample referred to as "generic MC," including the processes of ψð3770Þ decays to DD, non-DD, ISR production of low mass charmonium states and continuum processes, is used to study the background contribution. The effective luminosities of the generic MC samples correspond to at least five times the data sample luminosity. Two kinds of MC samples with the decay chain ψð3770Þ → D þ D − with D þ → K 0 S π þ π þ π − and D − → K þ π − π − using different decay models are generated for the amplitude analysis. One sample, "PHSP MC," is generated with a uniform distribution in phase space for the D þ → K 0 S π þ π þ π − decay, which is used to calculate the MC integrations. The other sample,"signal MC," is generated according to the results obtained in this analysis for the D þ → K 0 S π þ π þ π − decay. It is used to validate the fit performance, calculate the goodness of fit and estimate the detector efficiency.

III. EVENT SELECTION
Good charged tracks other than K 0 S daughters are required to have a point of closest approach to the IP within 10 cm along the beam axis and within 1 cm in the plane perpendicular to the beam. The polar angle θ between the track and the e þ beam direction is required to satisfy j cos θj < 0.93. Separation of charged kaons from charged pions is implemented by combining the energy loss (dE=dx) in the MDC and the time-of-fight information from the TOF. We calculate the probabilities PðKÞ and PðπÞ with the hypothesis of K or π, and require that K candidates have PðKÞ > PðπÞ, while π candidates have PðπÞ > PðKÞ. Tracks without particle identification (PID) information are rejected. Furthermore, a vertex fit with the hypothesis that all tracks originate from the IP is performed, and the χ 2 of the fit is required to be less than 100.
The K 0 S candidates are reconstructed from a pair of oppositely charged tracks which satisfy j cos θj < 0.93 and whose distances to the IP along the beam direction are within 20 cm. The two charged tracks are assumed to be a π þ π − pair without PID. In order to improve the signal-tobackground ratio, the decay vertex of the π þ π − pair is required to be more than two standard deviations away from the IP [15], and their invariant mass is required to be in the region ½467.6; 527.6 MeV=c 2 .
The D þ D − pair with D þ → K 0 S π þ π þ π − and D − → K þ π − π − is reconstructed with the requirement that they do not have any tracks in common. If there are multiple D þ D − candidates reconstructed in an event, the one with the average invariant mass closest to the nominal D AE mass [4] is selected. To characterize the D candidates, two variables, M BC and ΔE, defined as and In order to suppress the peaking background of D þ → K 0 S K 0 S π þ with an additional K 0 S → π þ π − , which has the same final state as our signal decay, we perform a decay vertex constrained fit on any remaining π þ π − pair with invariant mass within AE30 MeV=c 2 of the mass of the K 0 S . The events are removed if the obtained decay length is greater than twice its uncertainty. After applying all selection criteria, the expected yield from the background D þ → K 0 S K 0 S π þ is estimated to be 72.9 AE 8.5 by using the generic MC sample. In the amplitude analysis, it is subtracted by giving negative weights to the background events, as discussed in Sec. IVA. Self cross-feed events with misreconstructed signal decays are estimated from signal MC samples to be ∼0.1%. This effect is considered as a systematic uncertainty.
To estimate the contribution from the nonpeaking background, a 2D unbinned maximum likelihood fit is performed to the M BC ðD tag Þ versus M BC ðD signal Þ distribution in Fig. 1(c). The signal shape is modeled with the shape extracted from MC-simulated events. The function used to describe the diagonal background band is the product of an ARGUS function [16] in the M BC ðD tag Þ þ M BC ðD signal Þ plane and a Gaussian in the M BC ðD tag Þ − M BC ðD signal Þ plane. The background with only the tag candidate (signal candidate) properly reconstructed peaks at the charged D mass and spreads out on the other axis, which is parametrized as the product of a MC-simulated shape in and an ARGUS function on the other axis. The number of background events within the signal region extracted from the fit is 37.5 AE 7.5. The projection on M BC ðD signal Þ from the 2D fit is shown in Fig. 1(d). The small background bump under the signal is from the events with the D signal properly reconstructed but the D tag improperly reconstructed. In the amplitude analysis, the general background is ignored and its effect is considered as a systematic uncertainty.
To improve the momentum resolution and ensure that all events fall within the phase space boundary, the selected candidate events are further subjected to a six-constraint (6C) kinematic fit. It constrains the total four-momentum of all final state particles to the initial four-momentum of the e þ e − system, the invariant mass of signal side D þ → K 0 S π þ π þ π − constrains to the D þ nominal mass, and the K 0 S invariant mass constrains to the K 0 S nominal mass. We discard events with a χ 2 of 6C kinematic fit larger than 100. After applying all selection criteria, 4559 candidate events are obtained with a purity of 97.5%.

IV. AMPLITUDE ANALYSIS
The goal of this analysis is to determine the intermediate components in the four-body D þ → K 0 S π þ π þ π − decay. The decay modes that may contribute to the D þ → K 0 S π þ π þ π − decay are listed in Table I. The letters S, D in square brackets refer to the relative angular momentum between the daughter particles. The amplitudes and the relative phases between the different decay modes are determined with a maximum likelihood fit.

A. Likelihood function construction
The unbinned maximum likelihood fit is performed by minimizing the negative log likelihood (NLL) of the observed events (N data ) and the MC-simulated background events (N bkg ), where the indices k and k 0 refer to the kth event of the data sample and the k 0 th background event, respectively. The index j refers to the jth particle in the final state, f S ðp j Þ is the signal probability density function (PDF) in terms of the final four-momentum p j , and w bkg k 0 is the weight of the k 0 th background event. The contribution from the background is subtracted by assigning a negative weight to the background events.
The signal PDF f S ðp j Þ is given by where Mðp j Þ is the total decay amplitude describing the dynamics of the D þ decays, ϵðp j Þ is the detection efficiency parametrized in terms of the final four-momentum p j . R 4 ðp j Þdp j is the standard element of four-body phase space, which is given by The ϵðp j Þ in the numerator of Eq. (4) is independent of the fitted variables, leading to a constant term in minimizing the likelihood and can be ignored in the fit. The normalization integral of Eq. (4) is performed with a MC technique, which is then given by where k MC is the index of the kth event of the MC sample and N MC is the number of the selected MC events. M gen ðp j Þ is the PDF function used to generate the MC sample for the integration. This analysis uses an isobar model formulation in which the total decay amplitude Mðp j Þ is given by the coherent sum over all contributing amplitudes, where ρ n and ϕ n are the magnitude and phase of the nth amplitude, respectively. The nth amplitude A n ðp j Þ is given by where the indices 1 and 2 correspond to the two intermediate resonances. S n ðp j Þ is the spin factor, P α n ðm α Þ and B α n ðp j Þ (α ¼ 1, 2) are the propagator and the Blatt-Weisskopf barrier factor [17], respectively, and B D n ðp j Þ is the Blatt-Weisskopf barrier factor of the D þ decay. The parameters m 1 and m 2 in the propagators are the invariant masses of the corresponding resonances. For nonresonant contributions with orbital angular momentum between the daughters, we set the propagator to unity. This means that the amplitude has negligible m dependence. Since the D þ → K 0 S π þ π þ π − decay contains two identical π þ s in the final state, A n ðp j Þ is symmetrized by exchanging the two π þ s to take into account the Bose symmetry.

Spin factor
The spin factor S n ðp j Þ is constructed with the covariant tensor formalism [5]. The amplitudes with angular momenta larger than 2 are not considered due to the limited phase space. For a specific process a → bc, the covariant tensorst L μ 1 ÁÁÁμ l for the final states of pure orbital angular momentum L are constructed from the relevant momenta p a , p b , p c [5], μ 1 ÁÁÁμ L ν 1 ÁÁÁν L is the spin projection operator and is defined as μ ðVÞ for spin 1, and for spin 2.
The spin factors of the decay modes used in the analysis are listed in Table I. We useT ðLÞ μ 1 ÁÁÁμ L to represent the decay of the D þ meson andt ðLÞ μ 1 ÁÁÁμ L to represent the decay of the intermediate state.

Blatt-Weisskopf barrier factors
For the process a → bc, the Blatt-Weisskopf barrier factor [17] B L ðp j Þ is parametrized 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 for the intermediate resonances and 5.0 GeV −1 for the D þ meson. X L ðqÞ is given by With the invariant mass squared s a=b=c of the particle a=b=c, q is

Resonance line shapes
The propagator PðmÞ describes the line shape of the intermediate resonance. The resonances K Ã− ,K 1 ð1400Þ 0 , a 1 ð1260Þ þ andKð1460Þ 0 are parametrized with a relativistic Breit-Wigner (RBW) line shape where m 0 is the mass of resonance and ΓðmÞ is the massdependent width. The latter is expressed as where Γ 0 is the width of resonance and q 0 denotes the value of q at m ¼ m 0 . The ω and K 1 ð1270Þ − are parametrized as a RBW with a constant width ΓðmÞ ¼ Γ 0 .
The resonance ρ 0 is described by the Gounaris-Sakurai (GS) function P GS ρ ðmÞ with the ρ − ω interference taken into account [18,19], where ρ ω and ϕ ω are the relative magnitude and phase, respectively. P GS ρ ðmÞ is given by where and the function hðmÞ is defined as with dh dðm 2 Þ where m π is the charged pion mass [4]. The normalization condition at P GS ð0Þ fixes the parameter d ¼ fð0Þ=ðΓ 0 m 0 Þ. It is found to be [18] d The resonance f 0 ð500Þ is parametrized with the formula given in Ref. [20], which is identical to Eq. (17) with ΓðmÞ being decomposed into two parts, and Here, ρ ππ is the phase space of the π þ π − system and ρ 4π is the 4π phase space approximated by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 − 16m 2 π =m 2 p = ½1 þ e ð2.8−m 2 Þ=3.5 [20]. The parameters b 1 , b 2 , and a are fixed to the values given in Ref. [21].

AMPLITUDE ANALYSIS OF …
PHYS. REV. D 100, 072008 (2019) 072008-7 The resonance K Ã 0 ð1430Þ − is considered in a Kπ S-wave [denoted as ðK 0 S π − Þ S-wave ] parametrization extracted from the scattering data [22]. The same parametrization was used in Ref. [23], where a and r denote the scattering length and effective interaction length, respectively. Fðϕ F Þ and Rðϕ R Þ are the relative magnitudes (phases) for the nonresonant and resonant terms, and q and Γðm Kπ Þ are defined as in Eqs. (16) and (18), respectively. In the fit, the parameters M, Γ, F, ϕ F , R, ϕ R , a and r are fixed to the values obtained from the fit to the D 0 → K 0 S π þ π − Dalitz plot in Ref. [23], given in Table II.

B. Fit fraction
The fit fraction (FF) for an amplitude or a component (a certain subset of amplitudes) is calculated using a large set of generation-level PHSP MC samples by whereÃ n ðp k j Þ is either the nth amplitude (Ã n ðp k j Þ ¼ ρ n e iϕ n A n ðp k j Þ) or the nth component of a coherent sum of amplitudes (Ã n ðp k j Þ ¼ P ρ n i e iϕ n i A n i ðp k j Þ) and N gen is the number of the PHSP MC events. Note that the sum of the FFs is not necessarily equal to unity due to the interferences among the contributing amplitudes.
To obtain the statistical uncertainties of the FFs, the FFs are calculated 500 times by randomly varying the floated parameters according to the full covariance matrix.
The distribution for each amplitude or each component is fitted with a Gaussian function. The width of the Gaussian function is the statistical uncertainty of the corresponding FF.

V. RESULTS
We start the fit of the data by considering the amplitudes containing K Ã− , ρ 0 ,K 1 ð1270Þ 0 ,K 1 ð1400Þ 0 , a 1 ð1260Þ þ resonances, as these resonances are clearly observed in the corresponding invariant mass spectra. We then add amplitudes with resonances listed in the PDG [4] and nonresonant components until no additional amplitude has a significance larger than 5σ. To avoid either artificial or missing components, the total FF of each fit in the procedure is required to be less than 1.5. The cases of high correlation are also avoided, which is discussed in the next paragraph. In addition, in the iteration of adding amplitudes by comparing with the previous step, a better fit quality is required. The statistical significance for any new amplitude is calculated from the change of the log-likelihood value ΔðNLLÞ and the change of the degrees of freedom Δν. In the fits, the amplitude and phase of D þ → K 0 S a 1 ð1260Þ þ ðρ 0 π þ ½SÞ are fixed to 1 and 0 as the reference, while the magnitudes and phases of the other amplitudes are floating. Here, [S] means the angular momentum of the ρ 0 π þ combination is 0 (S wave). The corresponding D-wave amplitude D þ → K 0 S a 1 ð1260Þ þ ðρ 0 π þ ½DÞ is found to have a FF of about 1% of the S wave, which is consistent with both BESIII and LHCb amplitude analyses on D 0 → K − π þ π þ π − [24,25]. We consider therefore this D-wave amplitude in the nominal fit although its significance is 4.3σ.
The resonant term D þ → K 0 S a 1 ð1260Þ þ ðρ 0 π þ ½SÞ and its nonresonant partner D þ → K 0 S ðρ 0 π þ ½SÞ A (the subscript A represents the axial-vector nonresonant state for the ρ 0 π þ combination) are both found with significances greater than 10σ, while they are highly correlated because of the same angular distribution and large common region in phase space. For the resonant term in the fit model with the nonresonant partner, its FF becomes highly uncertain and is significantly different to the one in the fit model without the nonresonant partner. However the combined FF of these two amplitudes is almost unchanged. We, therefore, only consider the resonant term. Similar cases are also found with the amplitude pairs of D þ → Kð1460Þ 0 ðK 0 S ρ 0 Þπ þ and D þ → ðK 0 S ρ 0 Þ P π þ , D þ → Kð1460Þ 0 ðK Ã− π þ Þπ þ and D þ → ðK Ã− π þ Þ P π þ , as well as D þ →K 1 ð1650Þ 0 ðK Ã− π þ ½SÞπ þ and D þ →ðK Ã− π þ ½SÞ A π þ . Throughout this paper, we denote K Ã− → K 0 S π − and ρ 0 → π þ π − , which is also included in the FFs and BFs of corresponding submodes. In the nominal fit, we only use the resonant terms, as done in the analysis of Mark III [26]. S π − Þ S-wave parameters, obtained from the fit to the D 0 → K 0 S π þ π − Dalitz plot in Ref. [23]. The uncertainties are statistical. The masses and widths of ρ 0 , ω, K Ã− ,K 1 ð1270Þ 0 , where the uncertainties are statistical only. In the nominal fit, these parameters are set to be the values determined by likelihood scans. The scan results are shown in Fig. 2. In Fig. 2(a), three scan points at the right of the minimum point are higher than smooth scan expectations due to the correlation between the states with resonances a 1 ð1260Þ þ orKð1460Þ 0 involved.
Finally, our nominal fit model includes 13 amplitudes (labeled as I; II; III; …; XIII), in which eight of them can be summarized into four different components. To quantify the fit quality for this unbinned likelihood fit, an unbinned "mixed-sample method" is performed, which is described in Refs. [27,28]. With this method, the p-value is 25.5%. The projections of the invariant mass spectra and the distribution of χ are shown in Fig. 3. All the amplitudes and the corresponding significances and phases, as well as the FFs of amplitudes and components are listed in Table III, where the last row of each box is the coherent sum of the preceding amplitudes (components). For the phases and FFs, the first and second uncertainties are statistical and systematic, respectively. The systematic uncertainties are discussed below. Other tested amplitudes when determining the nominal fit model, but finally not used, are listed in Appendix A. The interference fit fractions between each amplitude are given in Appendix B.  To estimate the systematic uncertainties, the fit is altered to investigate the effect from each source. For the masses and widths of the intermediate resonances given by the PDG [4], they are shifted within the uncertainties from the PDG [4]. The masses and widths of a 1 ð1260Þ þ and Kð1460Þ, as well as the relative magnitude and phase of ω in ρ − ω parametrization are shifted within the uncertainties given by the likelihood scans. The barrier effective radius R is varied within AE1 GeV −1 . The input parameters of K 0 S π þ S-wave model are varied within their uncertainties given by Ref. [23]. For the resonance f 0 ð500Þ, the propagator is replaced by the RBW function with mass and width fixed at 526 MeV=c 2 and 535 MeV [21], respectively. For the resonance a 1 ð1260Þ þ , a constant width Breit-Wigner with mass and width determined by the fit is used to estimate the effect from the a 1 ð1260Þ þ line shape. For the effects from different line shapes, only the changes in the fit fractions are given. Since different propagators have different normalization factors, for the amplitude with f 0 ð500Þ involved, the shift effects on the FF are only considered. The effect from the peaking background D þ → K 0 S K 0 S π þ is estimated by altering the number of background events to be half of that in the nominal fit. The uncertainty from general background is studied by taking the background events into account, which are estimated from the average M BC (ðM BC ðD tag Þ þ M BC ðD signal ÞÞ=2) sideband region of ½1.830; 1.858 GeV=c 2 . Individual changes of the results with respect to the nominal one are taken as the corresponding systematic uncertainties.

VI. SYSTEMATIC UNCERTAINTIES
To evaluate the uncertainty from the fit procedure, we generate 300 sets of signal MC samples according to the nominal results in this analysis. Each sample, which has equivalent size as the data, is analyzed with the same S π þ 2 π − , and (h) π þ π þ π − invariant mass spectra, where the dots with error are data, and the curves are the fit projections. The small red histogram in each projection shows the D þ → K 0 S K 0 S π þ peaking background. The dip around the K 0 S peak comes from the used requirement to suppress the D þ → K 0 S K 0 S π þ peaking background. For the identical pions, the one resulting in a lower π þ π − invariant mass is denoted as π þ 1 ; the other is denoted as π þ 2 .
method as data analysis. We fit the resulting pull distributions, where V input is the input value in the generator, and V fit and σ fit are the output value and the corresponding statistical uncertainty, respectively. Fits to the pull distributions with Gaussian functions show no obvious biases and under-or overestimations on statistical uncertainties. We add in quadrature the mean and the mean error of the pull and multiply this number with the statistical error to get the systematic error. The results are given in Table VI, in which the corresponding uncertainties are the statistical uncertainties of the respective fits.
The effects from tracking/PID efficiency and the kinematic fit arising due to the imperfect modeling of the data by the simulation, as well as the resolution, are also III. Significances and phases for different amplitudes, labeled as I; II; III; …; XIII, respectively, as well as FFs for amplitudes and components (the last row of each box), where the first and second uncertainties are statistical and systematic, respectively. The f 0 ð500Þ and ρ 0 resonances decay to π þ π − , and the K Ã− resonance decays to K 0 S π − .
Amplitude Significance (σ) Phase FF       )   I  II  III  IV  V  VI  VII  VIII  IX     investigated. For tracking/PID efficiency and kinematic fit, a factor related to the correction is considered when calculating the normalization integral of Eq. (4). The difference between the alternative fit and the nominal fit is found to be negligible. The effect from the resolution is estimated from the difference of the pull distribution obtained from these 300 sets of signal MC samples using the generated and reconstructed four-momenta, which is also found to be negligible.

VII. CONCLUSION
We have determined the intermediate state contributions to the decay D þ → K 0 S π þ π þ π − from an amplitude analysis. With the fit fraction of the nth component FFðnÞ obtained from this analysis, we calculate the corresponding BF: 11Þ% is the total inclusive BF quoted from the PDG [4]. The results on the BFs are shown in Table VII.
Compared with the previous measurements [26], the precisions of the subdecay modes are significantly improved. The dominant intermediate process is D þ → K 0 S a 1 ð1260Þ þ ðρ 0 π þ Þ, which agrees with the measurement of Mark III [26]. We also extract the BFs of D þ → K 0 S a 1 ð1260Þ þ ðf 0 ð500Þπ þ Þ, D þ →K 1 ð1400Þ 0 ðK Ã− π þ Þπ þ , and D þ →K 1 ð1270Þ 0 ðK 0 S ρ 0 Þπ þ decays for the first time. Comparing with the decay of D 0 → K − π þ π þ π − [24,25], the decay mode D → Ka 1 ð1260Þ is found to be the dominant substructure in both D 0 and D þ decays. For the two K 1 states, the contributions from D → K 1 ð1270Þπ are at the same level for both D þ and D 0 decays. For D → K 1 ð1400Þπ, the related BF in D þ decays is found to be greater than that in D 0 decay by 1 order of magnitude. These results provide criteria to further investigate the mixture between these two axial-vector kaon states [1][2][3]. TABLE VII. The results of BFs for different components. The first, second and third errors are statistical, systematical and the uncertainty related to BðD þ → K 0 S π þ π þ π − Þ [4], respectively. The f 0 ð500Þ and ρ 0 resonances decay to π þ π − , and the K Ã− resonance decays to K 0 S π − .

APPENDIX B: INTERFERENCE OF FIT FRACTION
The interference between each amplitude is listed in Table VIII.