Amplitude analysis of $D^{+} \rightarrow K_{S}^{0} \pi^{+} \pi^{+} \pi^{-}$

The decay $D^{+} \rightarrow K_{S}^{0} \pi^{+} \pi^{+} \pi^{-}$ is studied with an amplitude analysis using a data set of 2.93${\mbox{\,fb}^{-1}}$ of $e^+e^+$ collisions at the $\psi(3770)$ peak accumulated by the BESIII detector. Intermediate states and non-resonant components, and their relative fractions and phases have been determined. The significant amplitudes, which contribute to the model that best fits the data, are composed of five quasi-two-body decays $ K_{S}^{0} a_{1}(1260)^{+}$, $ \bar{K}_{1}(1270)^{0} \pi^{+}$ $ \bar{K}_{1}(1400)^{0} \pi^{+}$, $ \bar{K}_{1}(1650)^{0} \pi^{+}$, and $ \bar{K}(1460)^{0} \pi^{+}$, a three-body decays $K_{S}^{0}\pi^{+}\rho^{0}$, as well as a non-resonant component $ K_{S}^{0}\pi^{+}\pi^{+}\pi^{-}$. The dominant amplitude is $ K_{S}^{0} a_{1}(1260)^{+}$, with a fit fraction of $(40.3\pm2.1\pm2.9)\%$, where the first and second uncertainties are statistical and systematic, respectively.


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 → V P , D → SP , D → V V , and D → AP , where P , V , S, and A denote pseudoscalar, vector, scalar, and axial-vector mesons, respectively. Finalstate 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 → P P and D → V P 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) will be 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 non-resonant components, where the 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 [7]. The event sample is based on 2.93 fb −1 of e + e − collisions at the ψ(3770) mass [8,9]. At this energy, D meson pairs are produced without any additional hadrons. To suppress backgrounds from other charmed meson decays and continuum (QED and qq) processes, 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 [6] is described in Ref. [7]. 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, 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 photon in EMC 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 sim-ulate 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 5 times the data sample luminosity. Two kind of MC samples with the decay chain of ψ(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 an 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 | cos θ| < 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 | cos θ| < 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 signalto-background 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 ± mass [4] is selected. To characterize the D candidates, two variables, M BC and ∆E, defined as and are calculated, where (E D , p D ) is the reconstructed fourmomentum of D candidate, and E beam is the calibrated beam energy. The signal events form a peak around zero in the ∆E distribution and around the charged D mass in the M BC distribution. In order to suppress the 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 ±30 MeV/c 2 of the mass of the K 0 S . The events are removed if the obtained decay length greater than twice of 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±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. IV A. Self crossfeed events with mis-reconstructed 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 general 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 MC-simulated shape. The diagonal background band is described by an ARGUS function [16] multiplied by a Gaussian in the anti-diagonal axis. 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 parameterized as the product of a MC-simulated shape in M BC (D tag ) (M BC (D signal )) and an ARGUS function on the other axis. The number of background events within the signal region extracted from the fit is 37.5±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 sixconstraint (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 which 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 ′ refer to the k th event of the data sample and the k ′th background event, respectively. The index j refers to the j th 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 ′ is the weight of the k ′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 detec-tion efficiency parameterized in terms of the final fourmomentum p j . R 4 (p j )dp j is the standard element of fourbody 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 k th 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 n th amplitude, respectively. The n th amplitude A n (p j ) is given by where the indexes 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 two 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] where µ1···µLν1···νL is the spin projection operator and is defined as 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

Blatt-Weisskopf barrier factors
For the process a → bc, the Blatt-Weisskopf barrier factor [17] B 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 the 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 X L=2 (q) = 13 With the invariant mass squared s a/b/c of the particle a/b/c, the 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 parameterized 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 parameterized 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 [19,20]: where ρ ω and φ ω are the relative magnitude and phase, respectively. P GS ρ (m) is given by where and the function h(m) is defined as 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 [19] The resonance f 0 (500) is parameterized with the formula given in Ref. [18]: where Γ tot (m) is decomposed into two parts: and Here is ρ ππ the phase space of the π + π − system and ρ 4π is the 4π phase space approximated by 1 − 16m 2 π /m 2 /[1 + e (2.8−m 2 )/3.5 ]. The parameters are fixed to the values given in Ref. [21].
The resonance K * (1430) − is considered in a Kπ Swave (denoted as (K 0 S π − ) S−wave ) parameterization extracted from the scattering data [22]. The same parameterization was used by BABAR in Ref. [23]: with 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 non-resonant and resonant terms, and q and Γ(m Kπ ) are defined as in Eq. (16) and Eq. (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 performed by BABAR [23] and are 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 sample by whereÃ n (p k j ) is either the n th amplitude (Ã n (p k j ) = ρ n e iφn A n (p k j )) or the n th component of a coherent sum of amplitudes (Ã n (p k j ) = ρ ni e iφn i A ni (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.

C. Goodness of fit
To examine the performance of the fit process, the goodness of fit is defined as follows. Since the D + and all four final-state particles have spin zero, the phase space of the decay can be completely described by five linearly independent Lorentz invariant variables. For convenience, one of the two identical pions which results in a lower π + π − invariant mass is denoted as π + 1 , while the other as π + 2 . The four final-state particles K 0 S , π + 1 , π + 2 , π − are marked with the indices 1, 2, 3, 4, respectively. Then the five invariant masses m 24 , m 34 , m 124 , m 134 , and m 234 are chosen to map out the phase space. To calculate the goodness of fit, the five-dimensional phase space is first divided into cells with equal size. Then, adjacent cells are combined until the number of events in each cell is larger than 20. The fit residual χ in p th cell is calcu- The goodness of fit is quantified as χ 2 = n p=1 χ 2 p , where N p and N exp p are the number of observed and expected number of events from the fit in the p th cell, respectively, and n is the total number of cells. The number of degrees of freedom (NDF) ν is given by ν = (n − 1) − n par , where n par is the number of the free parameters in the fit.

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 non-resonant components until no additional amplitude has a significance larger than 5σ. The statistical significance for any new amplitude is calculated with the change of the log-likelihood value ∆(NLL) after taking the change of the degrees of freedom ∆ν into account. 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 ρ 0 π + combination is zero (S-wave). The corresponding D-wave amplitude ) 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 non-resonant partner D + → K 0 S (ρ 0 π + [S]) A (the subscript A represents the axial-vector non-resonant 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 non-resonant partner, its FF becomes highly uncertain and is significantly different with the one in the fit model without the non-resonant 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 sub modes. In the nominal fit, we only use the resonant terms, as done in the analysis of Mark III [26].
The masses and widths of ρ 0 , ω, K * − ,K 1 (1270) 0 , K 1 (1400) 0 , andK 1 (1650) 0 are fixed at the values from PDG [4]. Since there are no world average values for the masses and widths of a 1 (1260) + andK(1460) 0 and the resonances lie on the upper boundary of the corresponding invariant mass spectrum, their values are determined by likelihood scans. The values of the parameters related to ρ − ω mixing are also determined by likelihood scans. The scan results are: 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, in which 8 of them can be summarized into four different components. The nominal fit yields a goodness of fit value of χ 2 /ν = 275.15/204 = 1.35. 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 earlier published 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.

VI. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties are categorized into the following sources: (I) masses and widths of the intermediate resonances, (II) effective radius of intermediate resonances and D + , (III) parameters in K 0 S π + Swave parameterization, (IV) parameters in ρ − ω mixing parameterization, (V) line shape of f 0 (500), (VI) effect from peaking background, (VII) effect from general background, and (VIII) fit procedure. The systematic uncertainties of the phases of amplitudes and the FFs of amplitudes and components due to different contributions are given in Tables IV and V, respectively. These uncertainties are given in units of standard deviations σ stat and are added in quadrature to obtain the total systematic uncertainties, as they are uncorrelated.
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) + andK(1460), as well as the relative magnitude and phase of ω in ρ − ω parameterization are shifted within the uncertainties given by the likelihood scans. The barrier effective radius R is varied within ±1 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 RBW function with mass and width fixed at 526 MeV/c 2 and 535 MeV [21], respectively. 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 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 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 over-estimations 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.

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 n th component FF(n) obtained from this analysis, we calculate the corresponding BF: B(n) = B(D + → K 0 S π + π + π − ) × FF(n), where B(D + → K 0 S π + π + π − ) = (2.97 ± 0.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 sub decay 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)π is at the same level for both D + and D 0 decays. For D → K 1 (1400)π, the re- S π + 2 π − , and (h) π + π + π − invariant mass spectra, where the dots with error are data, and the curves are the fit projections. The small red (colors online) 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. The identical pions are sorted with same method mentioned in Sec. IV C. Figure ( lated BF in D + decays is found to be greater than that in D 0 decay by one order of magnitude. These results provide criteria to further investigate the mixture between these two axial-vector kaon states [1][2][3]. S π + S-wave parameterization, (IV) parameters in the ρ − ω mixing parameterization, (V) line shape of the f0(500), (VI) effect from peaking background, (VII) effect from general background, and (VIII) fit procedure.