Study of charmonium decays to $K^0_S K \pi$ in the $B \to (K^0_S K \pi) K$ channels

A study of the $B^+\to K^0_SK^+K^-\pi^+$ and $B^+\to K^0_SK^+K^+\pi^-$ decays is performed using proton-proton collisions at center-of-mass energies of 7, 8 and 13 TeV at the LHCb experiment. The $K^0_SK \pi$ invariant mass spectra from both decay modes reveal a rich content of charmonium resonances. New precise measurements of the $\eta_c$ and $\eta_c(2S)$ resonance parameters are performed and branching fraction measurements are obtained for $B^+$ decays to $\eta_c$, $J/\psi$, $\eta_c(2S)$ and $\chi_{c1}$ resonances. In particular, the first observation and branching fraction measurement of $B^+ \to \chi_{c0} K^0 \pi^+$ is reported as well as first measurements of the $B^+\to K^0K^+K^-\pi^+$ and $B^+\to K^0K^+K^+\pi^-$ branching fractions. Dalitz plot analyses of $\eta_c \to K^0_SK\pi$ and $\eta_c(2S) \to K^0_SK\pi$ decays are performed. A new measurement of the amplitude and phase of the $K \pi$ $S$-wave as functions of the $K \pi$ mass is performed, together with measurements of the $K^*_0(1430)$, $K^*_0(1950)$ and $a_0(1700)$ parameters. Finally, the branching fractions of $\chi_{c1}$ decays to $K^*$ resonances are also measured.


Introduction
Understanding strong interaction effects in exclusive weak decays of heavy hadrons is of importance to gain information on fundamental aspects of the phenomenology of strong and weak interactions. In two-body nonleptonic decays such as B → (cc)K, a simple factorization method has been adopted to compute nonleptonic decay amplitudes [1]. The method involves expressing the hadronic matrix elements of four-quark operators in the effective Hamiltonian inducing the decay as the product of two matrix elements of quark currents. It has been successful, such as in the description of the B → η c K and J/ψK decays [2]. However, it clearly misses important effects, since it fails to describe the B → χ c0 K mode. Indeed, in this mode the factorized amplitude involves the matrix element of the (cc) V,A vector (V) and axial (A) currents between the vacuum and χ c0 resonance, which vanishes due to charmed vector current conservation and parity conservation. In contrast, the measured B → χ c0 K branching fraction is sizable [3], clearly indicating that the nonfactorizable part of the amplitude plays an important role [4]. Additional observations of new B → χ c0 X decay modes are therefore of interest.
The only established strange scalar meson is the K * 0 (1430) resonance, whose parameters are yet to be precisely measured [3]. Scalar resonances decaying to Kπ are particularly interesting, since many amplitude analyses of heavy-flavor decays involve a Kπ system [5], whose theoretical description is a source of large systematic uncertainty. A widely used method of modelling the Kπ S-wave relies on the results from Ref. [6], which consists of a large threshold enhancement described by a scattering length term and a relativistic Breit-Wigner (BW) function describing the K * 0 (1430) resonance. A similar behaviour is observed in the Kπ S-wave measured in D + decays [7][8][9]. Still unresolved is the possible existence of a broad scalar resonance, κ/K * 0 (700), claimed by several experiments [3]; its existence would suggest the possible presence of tetraquark states in the light meson system [10].
Further information has been obtained from the Dalitz plot analysis of η c decays to KKπ with the η c meson produced in two-photon interactions, and in an extended Kπ mass region [11]. Due to its large width, the description of the lineshape of the K * 0 (1430) resonance could be complicated by the effects of the opening of the Kη and Kη ′ thresholds. The decay of the K * 0 (1430) resonance to Kη has been observed in a Dalitz-plot analysis of η c → K + K − η [12], and its branching fraction has been found to be small. Its decay to Kη ′ has been observed in Refs. [13,14]. Another resonance, K * 0 (1950), seen in the Kπ decay mode [6], is still to be confirmed. Evidence for the K * 0 (1950) → Kη ′ decay mode has been found in a Dalitz-plot analysis of η c → K + K − η ′ decays [14].
In the Dalitz-plot analysis of the η c → ηπ + π − decay, a new a 0 (1700) resonance has been observed in the ηπ mass spectrum [14] and recently confirmed in the a 0 (1700) → K 0 S K decay mode [15]. The a 0 (1700) decay to K 0 S K is therefore expected to contribute to η c /η c (2S) → K 0 S Kπ decays. To date, no Dalitz-plot analysis of η c (2S) has been performed. The χ c1 → K 0 S Kπ decay has been studied in Ref. [16] in ψ(2S) → γK 0 S Kπ decays with 220 ± 16 events and a low background. By fitting to the Kπ mass projections, partial branching fractions to K * resonances have been measured. Given the small dataset, only upper limits have been obtained for some χ c1 → K * X decay modes, and therefore further measurements of these branching fractions are useful.
The B + → K 0 S K + K − π + and B + → K 0 S K + K + π − decays have been previously studied in Refs. [17,18], but their branching fractions are yet to be measured. New large datasets may therefore help in clarifying several of the above-listed issues related to light-meson spectroscopy and B to charmonium decays. This paper is organized as follows: Sec. 2 describes the LHCb detector; Sec. 3, the signal candidate selection procedure; Sec. 4, the study of various mass spectra; Sec. 5, the measurement of the charmonium-resonance parameters; Sec. 6, the efficiency evaluation; Secs. 7 and 8, the Dalitz plot analysis of the η c and η c (2S) mesons, respectively; Sec. 9, the study of χ c1 decays; and Sec. 10, the measurements of various branching fractions. Finally, Sec. 11 summarizes the results.

Detector, simulation and analysis
The LHCb detector [19,20] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector elements that are particularly relevant to this analysis are: a silicon-strip vertex detector (VELO) [21] surrounding the pp interaction region that allows c and b hadrons to be identified from their characteristically long flight distance; a tracking system that provides a measurement of the momentum, p, of charged particles; and two ring-imaging Cherenkov detectors that are able to discriminate between different species of charged hadrons. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The entire dataset collected with the LHCb experiment during Runs 1 and 2 is used, corresponding to center-of-mass energies √ s=7, 8 and 13 TeV and comprising an integrated luminosity of 9 fb −1 . The online event selection starts with a trigger [22], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. During offline selection, trigger signatures are associated with reconstructed particles. Since the trigger system uses the transverse momentum of the charged particles with respect to the beam axis, p T , the phase-space and time acceptance is different for events where signal tracks were involved in the trigger decision (called trigger-on-signal or TOS throughout) and those where the trigger decision was made using information from the rest of the event only (noTOS). Data from both trigger conditions are used and studied separately for consistency tests and the evaluation of systematic uncertainties.
Simulation is required to model the effects of the detector acceptance and the imposed selection requirements. In the simulation, pp collisions are generated using Pythia [23] with a specific LHCb configuration [24]. Decays of unstable particles are described by EvtGen [25], in which final-state radiation is generated using Photos [26]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [27] as described in Ref. [28].
Two types of simulations are performed: (a) where the B + is decayed according to a phase-space model and (b) where the B + decays as B + → (cc)K, where (cc) indicates a charmonium resonance decaying to K 0 S Kπ by phase space.

Event selection
The present work reports a study of the two B + decays 1 and with K 0 S → π + π − . Decays of the K 0 S are reconstructed in two categories: the first involving K 0 S mesons that decay early enough for the pions to be reconstructed inside the VELO; and the second containing K 0 S mesons that decay later such that track segments from the pions are outside the VELO. These categories are referred to as long K 0 S (indicated in the following with K 0 SLL ) and downstream K 0 S (indicated in the following with K 0 SDD ), respectively. While the K 0 SLL category has better mass, momentum and vertex resolution, there are approximately twice as many K 0 SDD candidates. Candidate B + particles are formed by combining the K 0 S candidate with three other charged tracks, having a total charge of one, performing a global vertex fit to the decay tree and requiring the B + candidate to originate from one of the primary pp collision vertices in the event. Selection criteria and efficiency measurements are performed separately for each K 0 S category. To suppress backgrounds, in particular combinatorial background formed from random combinations of unrelated tracks, the events satisfying the trigger requirements are filtered by a loose preselection, followed by a multivariate selection optimized separately for each data sample. Selection requirements are tuned to minimize correlation of the signal efficiency with kinematic variables, resulting in better control of the corresponding systematic uncertainties. Consequently, the selection relies minimally on the kinematics of the final-state particles and instead exploits the topological features that arise from the detached vertex of the B + candidate. These include: the impact parameters of the B + candidate and its decay products, the quality of the decay vertices of the B + and K 0 S candidates and the separation of these vertices from each other and from the primary vertex.
The preselection of K 0 S and B + candidates requires, for each track, the presence of appropriate particle identification (PID) information and imposes invariant mass selections around the known K 0 S and B + particle masses. The separation of signal from combinatorial background is achieved by means of a boosted decision tree (BDT) classifier [29,30], implemented using the TMVA toolkit [31]. The multivariate classifier chosen for this analysis is a BDT with a gradient boosting algorithm [32], separated for K 0 SLL and K 0 SDD data. The classifiers are trained using simulated signal decays, composed of samples (a) and (b) described in Sec. 2 with proportions corresponding to the resonance composition observed in the data (see Sec. 4). The simulation also matches the relative proportions of the dataset at the various center-of-mass energies. It is assumed that the efficiencies for the reconstruction of the decays displayed in Eq. (1) and Eq. (2) are the same. For the background input to the training, data in the lower and upper sidebands of the B + signal region are used. The composition of the background sample reflects the data-taking conditions and the decays in Eq. (1) and Eq. (2) are used in equal proportions. The optimization of the BDT classifier working point is performed by scanning the where N sig (N bkg ) represents the signal (combinatorial background) yield in the signal region. The yields are evaluated in a ±2.5σ window around the B + mass, where σ is the mass resolution, through fits to the K 0 S K + K − π mass spectra using two Gaussian functions sharing the same mean for the signal and a linear function for the background. Figure 1 shows the π + π − invariant-mass distributions at the K 0 S candidate vertex, separated for K 0 SLL and K 0 SDD , after selecting candidates using the optimized figure of merit S defined in Eq. (3).
The two π + π − invariant mass distributions are fitted using the sum of two Gaussian functions sharing the same mean, with σ 1 and σ 2 resolutions, and a linear function for the background. An effective resolution is computed as σ = f σ 1 + (1 − f )σ 2 where f is the fraction of the first Gaussian contribution. The resulting effective resolutions for the LL and DD categories are σ LL = 2.53 MeV and σ DD = 6.46 MeV. The K 0 S signals are selected within 3.0σ of the fitted K 0 S mass of 497.8 MeV. 2 In order to facilitate the extraction of the B + and K 0 S signal and combinatorial background components from these invariant-mass spectra, no kinematic mass constraint is applied to the K 0 S and B + signals. To improve the mass resolution (see Sec. 4.1), the energy of the selected candidate K 0 S is computed as where p K 0 S is the reconstructed K 0 S momentum and m K 0 S the known K 0 S mass. Compared with the resolution obtained from the use of the K 0 S mass constraint, this method gives the same resolution for the K 0 S KKπ invariant mass and a slightly worse resolution, by ≈ 6%, for the K 0 S Kπ invariant mass. Particle identification of the three charged hadrons is performed using the output of a probabilistic neural network (NN) trained on the output of all the subdetectors. The figures of merit are expressed as P K = NN K (1 − NN π ) for kaon identification and  P π = NN π (1 − NN K ) for pion identification, where NN π and NN K are the NN probabilities for pion and kaon identification, respectively. Thresholds are applied to these quantities to maximize the significance of the B + -candidate invariant-mass peak as a function of P K or P π .
Open charm production in the B + decays is significant, with the presence of several signals of D + , D 0 and D + s in two-body and three-body mass combinations. The largest contributions are due to D 0 → K 0 S K + K − decays in B + → K 0 S K + K − π + decays (17.4 ± 0.2)% and D 0 → K + π − decays (6.4 ± 0.1)% in B + → K 0 S K + K + π − decays. As both D 0 signals have large signal to background ratios, they are removed from the B + -candidate samples.
mass, with m 0 = 1864.6 MeV and σ = 4.5 MeV for K 0 SLL data and σ = 6.3 MeV for K 0 SDD data. The background from D 0 → K + π − decays is removed by requiring |m − m 0 | > 3.5σ, where m indicates the K + π − mass, with m 0 = 1864.5 MeV and σ = 8.4 MeV. In both cases the D 0 parameters are extracted from a fit to the data, using a Gaussian function and a linear polynomial for signal and background, respectively. Figure 2 shows the K 0 S K + K − π + and K 0 S K + K + π − invariant mass spectra, after the optimized offline selection requirements, separated by K 0 S category. The distributions are fitted to a sum of two Gaussian functions sharing the same mean values and a linear background function. The fits give an average B + -mass value of 5280 MeV and an average width of σ = 17.7 MeV. Signal candidates are selected in a window of ±2σ Table 1: Fitted B + signal yield and purity for B + → K 0 S K + K − π + and B + → K 0 S K + K + π − final states separated by K 0 S type.

Final state
79.9 ± 0.2 of the fitted B + mass, common to the four datasets. Table 1 lists the fitted yields and purities (P ) in the B + signal region for the different datasets, where the purity is defined as P = N sig /(N sig + N bkg ). The dependence of P on the collision energy and data-taking conditions is approximately uniform for all the four datasets, simplifying the Dalitz-plot analyses reported in the following. Approximately 0.02% of events contain multiple B + decay candidates, as selected with the above procedure, and therefore their effect is considered negligible. An inspection of the π + π − mass spectra for K 0 SLL and K 0 SDD candidates after all selections, and in the B + signal region, shows K 0 S signals with negligible background.

Mass spectra
The K 0 S Kπ invariant-mass spectra for events in the B + signal region, summed over the K 0 SLL and K 0 SDD datasets, are shown in Fig. 3 for B + → K 0 S K + K − π + and B + → K 0 S K + K + π − final states. The lower and upper mass sidebands around the B + signal peak are defined in the ranges [−6σ, −4σ] and [4σ, 6σ], respectively, with σ = 17.7 MeV. The corresponding K 0 S Kπ invariant mass spectra, representative of background candidates, are superimposed onto the K 0 S Kπ invariant mass spectrum from the B + signal region. The B + → K 0 S K + K + π − final state has two kaons with the same charge and therefore both combinations are included in the K 0 S Kπ invariant-mass spectrum. There are 119,198 entries and 159,694 combinations in the K 0 S Kπ invariant mass spectra from the two B + decay modes, respectively. The K 0 S Kπ invariant-mass distributions show a signal at threshold in the position of the f 1 (1285) resonance and a broad complex structure in the 1.5 GeV mass region. A D 0 → K 0 S Kπ signal can be observed, due to the open-charm final state B + → D 0 K + . Prominent signals of η c , J/ψ, χ c1 and η c (2S) can be observed in both invariant-mass spectra. A broad enhancement in the η c mass region is present in the K 0 S Kπ mass spectrum from B + sidebands. The effect can be understood as due to the presence of prompt η c → K 0 S Kπ decays [33] reconstructed using an incorrect decay chain. The structure above 1.9 GeV that appears mostly in the sidebands is due to the reflection from D 0 → K 0 S π + π − decays, where one pion is misidentified as a kaon. The strong η c signal present in the data allows for a test of the agreement between data and simulation for the PID. Particle identification is removed for each final-state kaon and pion in turn on both data and simulation, and the η c event losses due to the kaon (3%) and the pion (0.4%) identification are compared. The results for data and simulation agree within 2σ and this effect is therefore ignored.

Mass resolution
The mass resolution is obtained from simulated data. In this analysis it is used to perform fits to the invariant-mass spectra and obtain resonance parameters in which the effects of the mass resolution are expected to be significant, in particular in fitting narrow charmonium states. As the resolution functions are mass dependent, they are evaluated in specific mass intervals, namely the η c -J/ψ mass region, defined in the [2.90-3.15] GeV region, and the χ c1 -η c (2S) mass region, defined in the [3.46-3.70] GeV region. Due to the presence of two types of reconstructed K 0 S with different resolutions, the mass resolution is computed separately for K 0 SLL and K 0 SDD data. The resulting mass-difference distributions are well described by the sum of a Gaussian and a Crystal Ball function. The width (σ) of the dominant Gaussian contribution is 8.8 MeV (9.9 MeV) in the η c -J/ψ mass region and 10.3 MeV (12.0 MeV) in the χ c1 -η c (2S) mass region for the K 0 SLL (K 0 SDD ) data.

Measurement of charmonium-resonance parameters
The measurements of charmonium-resonance parameters are performed with binned fits to the K 0 S Kπ invariant-mass spectra separately in the η c -J/ψ and the χ c1 -η c (2S) mass regions in the ranges shown in Figs. 4 and 5, respectively. In both invariant-mass regions, the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − data are fitted separately. Since the experimental resolutions for K 0 SLL and K 0 SDD data are different, a simultaneous fit to the two invariant-mass spectra is performed, sharing only the resonance parameters. In the fit to the η c -J/ψ mass region, the J/ψ width is fixed to the known value [3] and the K 0 S Kπ mass value is shifted by the J/ψ mass: m ′ = m − m J/ψ , where m J/ψ is fixed to the known value [3]. Here, all resonances are described by simple non-relativistic BW functions convolved with the appropriate mass resolution functions. The backgrounds, which are due to a combination of B → K 0 S KKπ decays and incoherent K 0 S Kπ production, are represented by first-order polynomials. Table 2 lists the fitted resonance parameters and yields, together with fit p-values. A comparison with PDG [3] averages shows improvements both in the mass and width of the η c resonance, with the J/ψ mass in good agreement with the known value. A small negative shift, consistent with zero, can be seen on the fitted J/ψ mass. The good description of the J/ψ lineshape, whose observed width is dominated by the mass resolution, demonstrates the good agreement of the mass resolution in data and simulation.
Systematic uncertainties include the following sources. The bin width is varied from 2.5 to 3.0 MeV and the background shape is changed from linear to a second-order polynomial. The η c component is allowed to interfere with the background that has a significant contribution from the B + decay. The interference is parameterized as where A nres is the non-resonant amplitude and |A nres | 2 is described by a linear function.
The resonant contribution is constructed as A res = α · BW(m) · exp(iϕ), where α and ϕ are free parameters and BW(m) is the Breit-Wigner function of Eq. (5) describing the η c lineshape convolved with the experimental resolution. The coherence factor, c, is a free parameter. It is found that the interference model produces a small improvement in the description of the data with fitted phases of ϕ = 1.596 ± 0.009 rad and ϕ = 1.631 ± 0.011 rad, both close to π/2, for B + → K 0 S K + K − π + and B + → K 0 S K + K + π − data, respectively. The deviations of the fitted resonance parameters from those obtained without interference are included as systematic uncertainties.
The systematic uncertainty associated with the background model is evaluated by varying the BDT classifier selection working point, resulting in a variation of the B + purity of around ±5%. The average value of the absolute values of the two resulting variations of the resonance parameters is used to quantify the systematic uncertainty associated to the background level. Starting with the fitted functions from the reference fit, 400 pseudoexperiments are generated and fitted. The average value of the deviations of the fitted parameters is included as a systematic uncertainty.
The momentum-scale uncertainty is evaluated as 0.03Q [34], where Q is evaluated as the difference between the η c mass and the sum of the masses of the decay particles. The uncertainties on the resonance parameters arising from the limited sizes of the simulation Table 2: Fitted η c , J/ψ, η c (2S) and χ c1 parameters. For the J/ψ the m ′ = m − m J/ψ value is reported. The first uncertainty is statistical, the second systematic.  Table 3: Summary of systematic uncertainties on J/ψ and η c resonance parameters.
Contribution samples, used to obtain the resolution functions, are obtained by fitting the K 0 S Kπ invariant-mass spectrum from 400 pseudoexperiments where, in each fit, all the parameters describing the resolution functions are varied randomly from a Gaussian distribution defined by their statistical uncertainties. Table 3 gives the resulting contributions to the systematic uncertainties which are then added in quadrature.
A similar model is used to fit the χ c1 -η c (2S) mass region, shown in Fig. 5, with the fit results summarized in Table 2 together with the inverse-variance-weighted averages  of the fitted parameters from the two B + decay modes. In this case, the χ c1 width is fixed to the known value [3]. The inverse-variance method is also used, here and in the following, to evaluate average systematic uncertainties. A comparison of the results listed in Table 2 with PDG [3] measurements shows improvements both in the mass and width of the η c (2S) resonance, with the χ c1 mass in good agreement with the known value. Systematic uncertainties on the fit parameters are evaluated in a similar way as for the η c -J/ψ fit except now the η c (2S) is allowed to interfere with the background. The fit with interference returns relative phases of ϕ = 1.52 ± 0.05 rad and ϕ = 1.82 ± 0.09 rad for the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − final states, respectively (similarly close to π/2). The presence of signals corresponding to the h c (1P ) and χ c2 (1P ) resonances is explored by adding additional components to the fit function and resonance parameters fixed to the known values [3]; their yields are found to be consistent with zero. A summary of the systematic uncertainties, together with their quadratic sum, is given in Table 4. Table 4: List of the systematic uncertainties on the χ c1 and η c (2S) parameters.   Figure 6 shows the K + K − invariant-mass spectrum for B + → K 0 S K + K − π + candidates in the χ c0 -χ c2 mass region, where a prominent χ c0 signal can be seen, together with a weaker signal at the χ c2 mass. Also superimposed is the K + K − invariant-mass spectrum from the B + sidebands, where no χ c0 signal can be seen. For this B + decay mode, the PDG [3] reports a branching-fraction upper limit of B(B + → χ c0 K * ) < 2.1 × 10 −4 [35].

First observation of
Using a similar method as for the fit to the other charmonium resonances, a simultaneous binned fit of the K + K − mass spectrum for the K 0 SLL and K 0 SDD data is performed. The background is parameterized by a second-order polynomial; the resonance parameters for the χ c0 state are unconstrained, while those for the χ c2 component are fixed to their known values. The χ c0 signal is described by the relativistic spin-0 Breit-Wigner function where m 0 is the resonance mass. The mass dependent width Γ is written as where Γ 0 is the resonance width and q (q 0 ) is the momentum of either decay particle in the two-body (resonance) rest frame. The χ c2 signal, due to its narrow width, is described by a simple Breit-Wigner function (Eq. 5); both distributions are convolved with the experimental resolution function modelled, as described in Sec. 4.1, by the sum of a Crystal Ball function and a Gaussian function having σ = 11.9 MeV. In this fit, the measured χ c0 mass and width are both shifted by approximately 5 MeV with respect to their known values [3].
The results of the fit are overlayed. The curves also include interference terms and therefore the χ c0 lineshape takes slightly negative values.
Allowing for interference of the χ c0 resonance with a non-resonant K + K − component, as described in Eq. 6, results in the data description shown in Fig. 6, with a p-value=27.9% and χ c0 parameters consistent with PDG averages [3]. The χ c0 yield is N χ c0 = 1920 ± 90 and the relative phase ϕ = −1.290 ± 0.073 rad. The fit also returns a χ c2 yield of N χ c2 = 190 ± 30. To evaluate the significance of the χ c2 signal, a fit without its contribution is performed. This results in a χ 2 variation of ∆χ 2 = 22.5 for the difference of one parameter, which gives a statistical significance of 4.6σ.

Efficiency
Two types of efficiencies are evaluated, total and local. The total efficiency describes the effects of the reconstruction on the B + decay to the 4-body final state, needed to evaluate the relative charmonium branching fractions. The local efficiencies are evaluated in specific K 0 S Kπ mass regions, i.e. the η c -J/ψ and χ c1 -η c (2S) mass regions, where Dalitzplot analyses are performed and descriptions of the K 0 S Kπ detection efficiency are needed. The efficiencies are evaluated by generating simulation samples that undergo the same reconstruction and analysis selections as the data. The efficiency is evaluated as the ratio of selected over generated distributions projected over the relevant kinematic variables.
The choice of the phase-space variables which describe the efficiency is somewhat arbitrary, and a mixture of two-body and three-body invariant-mass projections is used, together with variables related to the angular distributions. A comparison between the Figure 7: Definition of the angles (a) θ π , (b) θ K 1 and (c) ϕ K .
p T distributions of the B + candidates in simulated samples and in data shows a small disagreement, which is corrected by weighting the former to match the latter.

Total efficiency
To study efficiencies, simulated events are generated according to a four-body phase-space model. Since the physics of the K 0 S Kπ system is of interest, the charged kaon participating in the decay under study is labeled K 1 , while the second charged kaon behaves as a spectator and is labeled K 2 . Therefore, the reaction can be written Only B + → K 0 S K + K − π + decays are generated because the efficiency for B + → K 0 S K + K + π − is expected to be the same. The generated angular distributions are all uniform; therefore, any observed variation is due to inefficiency. The efficiencies are evaluated separately for the K 0 SLL and K 0 SDD samples. In both cases the resulting distributions have only a mild dependence on the K 0 S Kπ invariant mass. The kinematics of a four-body decay are fully described by five independent variables. Different invariant-mass combinations have different kinematic bounds, so mass-reduced variables are used instead because they always range between 0 and 1. They are defined as [36] where m, m min and m max indicate the invariant mass and its minimum and maximum kinematically allowed values, respectively. As an example in the above equation, m x (K 0 S K) for the three-body K 0 S Kπ system is computed as and m min = m K 0 S + m K and m max = m(K 0 S Kπ) − m π , similarly for other two-body mass combinations. The angular distributions make use of helicity angles defined as follows. In the K 0 S Kπ rest frame, shown in Fig. 7(a) and (b), θ π (θ K 1 ) indicates the angle formed by the K 1 (π) with respect to the K 0 S K 1 (K 0 S π) direction in the K 0 S K 1 (K 0 S π) rest frame. Similarly, the angle θ K 2 is defined by exchanging K 1 with K 2 . Figure 7(c) shows the definition of ϕ K , the angle formed by the spectator K 2 with the normal to the K 0 S K 1 π plane.
The curves are the results of fits made using seventh-order polynomials.
The model used to describe the efficiency is obtained in an iterative manner. First, the variables showing the strongest deviation from uniformity in the simulation are chosen. The efficiency projection as a function of the first chosen variable, m x (K 0 S K 1 ), is fitted using a seventh-order polynomial labeled as ϵ 1 (m x (K 0 S K 1 )). Figure 8 shows the efficiency projected on m x (K 0 S K 1 ) separately for the K 0 SLL and K 0 SDD samples. The events are then weighted by the inverse of the efficiency 1/ϵ 1 (m x (K 0 S K 1 )) and a second variable (m x (K 0 S π)) is chosen which is itself fitted by a seventh-order polynomial labeled as ϵ 2 (m x (K 0 S π)). The events are then weighted by ). The process continues in this fashion, terminating when, after weighting, the efficiency is consistent with being uniform across all the nine considered variables, both on one-and two-dimensional projections. The total efficiency for each K 0 S category, ϵ LL and ϵ DD , is found to be well described by

Efficiency for different trigger conditions
As described in Sec. 2, the reconstructed events belong to two trigger categories: TOS and noTOS. Separate efficiencies are needed for each trigger condition, further divided into K 0 SLL and K 0 SDD categories. Figure 9 shows the efficiency distributions as functions of the K 0 S Kπ invariant mass, separated by trigger condition and K 0 S type. The efficiency evaluations are performed using the same sample of generated events and therefore the scale of the distributions also gives the fraction of the simulation samples belonging to each category. It is found that the efficiency has a weak dependence on m(K 0 S K 1 π) for all the considered samples. The efficiency in the η c -J/ψ mass region is evaluated from a dedicated sample of simulated B + → η c K + decays, where the η c decays to a K 0 S Kπ state, in which the η c is generated according to a BW function and decays uniformly in its three-body phase space. These simulations are used to evaluate the efficiency across the Dalitz plot, which can be described in terms of two independent variables, chosen to be m x (K 0 S K 1 ) and cos θ π . Labeling x = m x (K 0 S K 1 ) and y = cos θ π , the efficiency map is smoothed by fitting with a two-dimensional polynomial separately for the K 0 SLL and K 0 SDD samples. Figure 10 shows the fitted two-dimensional efficiency distributions with the fit projections on m x (K 0 S K 1 ) and cos θ π . The efficiency in the χ c1 -η c (2S) mass region is computed from the sample of simulated B → K 0 S KKπ decays by selecting events in the [3.46-3.70] GeV mass region. Due to the size of the simulated sample and because of the weak dependence of the efficiency on cos θ π , the efficiency model is obtained from fits with seventh-order polynomials to the efficiency projections on m x (K 0 S K 1 ) and cos θ π , as shown in Fig. 11. The efficiency is then parameterized as ϵ = ϵ 1 (m x (K 0 S K 1 )) · ϵ 2 (cos θ π ).
7 Dalitz plot analysis of the η c decay to K 0 S Kπ This section is devoted to the study of the decay where the K + meson is considered a spectator. For simplicity here and in the following, kaons are labeled using their charge. Possible backgrounds originating from charm decays to particles amongst the η c → K 0 S K − π + decay products and the spectator K + are investigated, and no open charm-production is observed in any two-body or three-body mass combinations. The absence of significant structures ensures that the η c signal is decoupled from the K + and the η c analysis can be considered as a simple study of the three-body η c decay. The π − K + invariant-mass spectrum, where the π − comes from the K 0 S decay, shows no evidence of a D 0 background. Possible backgrounds from pions misidentified as kaons are considered by assigning the pion mass to each of the two kaons in the final state. Except for the small D 0 → K 0 S π + π − signal observed in the B + sidebands, discussed in Sec. 4, no charm signal is observed in the resulting two-body or three-body mass combinations. Also, no evidence for the B + decay to the K 0 S K + π + π − state is observed. The η c signal region is defined as [2.935-3.035] GeV while the sideband regions are defined as [2.83-2.88] GeV and [3.15-3.20] GeV as shown in Fig. 12. The background under the η c signal peak has two components: (a) from B + background, estimated from B + sidebands, labeled in the following as incoherent and shown as dark-gray areas in Fig. 12; (b) from B + signal, labeled as coherent and indicated by light-gray shading in Fig. 12. The relative fraction of these two components is estimated by integrating the contents of the dark-and light-gray areas in the η c lower and higher sidebands. The resulting incoherent backgrounds fractions, labeled as f B , are listed in Table 5 together with the η c event yields in the signal region and purities, separated for the K 0 SLL and K 0 SDD data.

Data selection for the
This section is devoted to the selection of the decay where two identical K + mesons are present. Because of this, the two K + mesons are alternately considered as part of the η c signal decay or as a spectator to it, and every The dark-gray area represents the K 0 S Kπ invariant-mass spectrum from the B + sideband; the light-gray areas indicate the η c signal and sideband regions. Table 5: Candidate events and purities in the η c signal region and fractions of B + sideband contributions in the η c sideband regions for B + → K 0 S K + K − π + and B + → K 0 S K + K + π − decays. Low and high indicate the lower and higher sideband candidates. Because of the limited statistics, the values of f B are summed for the K 0 SLL and K 0 SDD data.
candidate decay therefore appears twice in the sample under study. This final state is affected by a significant background from D 0 → K + π − decays, which is removed as discussed in Sec. 3. Labeling the spectator K + as K + 2 , the K 0 S K + 2 π − invariant-mass spectrum shows a small D 0 signal, and events are removed if they lie within ±26 MeV of the D 0 mass. No other open-charm signal is observed in other two-body or three-body mass combinations. The K 0 S Kπ mass spectrum is shown in Fig. 12(b) with indications of background, signal and sideband regions. Table 5 gives information about the event yields, purities and background composition.

Dalitz plot analysis
The η c Dalitz plot is shown in Fig. 13, for (a It is dominated by horizontal and vertical bands due to the presence of the K * 0 (1430) +,0 resonances.

Dalitz plot analysis method
An amplitude analysis of the K 0 S Kπ final state in the η c mass region is performed, using unbinned maximum-likelihood fits. The likelihood function is written as where N is the number of events in the signal region. For the n-th event, x n = m 2 (K 0 π + ) and y n = m 2 (K − π + ). The signal purity (P ) is obtained from the fit to the K 0 S Kπ invariant-mass spectrum (see Fig. 4) and listed in Table 5. The efficiency ϵ(x ′ n , y ′ n ) is parameterized as a function of x ′ n = m(K 0 S K) and y ′ n = cos θ π and described in Sec. 6. The complex signal-amplitude contribution for the i-th signal component is indicated as A i (x n , y n ). The corresponding complex parameter c i is allowed to vary in the fit, except for the largest amplitude, the (Kπ S-wave)K (labeled in the following as (Kπ) S K) in the quasi-model-independent (QMI) analysis (Sec. 7.4) or the K * 0 (1430)K in the isobar model analysis (Sec. 7.5), which is taken as the reference, setting the modulus |c 1 | = 1 and phase ϕ 1 = 0. The background probability density function for the i-th background component is indicated as B i (x n , y n ), and it is assumed that interference between signal and background amplitudes can be ignored. The magnitude of the i-th background component is indicated by k i and is obtained by fitting to the sideband regions and interpolating linearly in the signal region. The expressions y) dxdy represent normalization integrals for signal and background. Numerical integration is performed on phase-space generated events [37] with η c signal and background generated according to their experimental distributions. For resonances with free parameters, integrals are re-computed at each minimization step. Integrals corresponding to the background components, and fits involving amplitudes with fixed resonance parameters, are computed only once. Both the K 0 SLL and K 0 SDD data are included in the likelihood function, each dataset with its own efficiency and purity. The A i amplitudes are taken from Ref. [38,39]. The standard Blatt-Weisskopf [40] form factors are used at the resonance vertices with radius r fixed to 1.5 GeV −1 . The decay of the η c to KKπ is governed by strong interactions with two allowed processes: η c → a L π and η c → K * L K, where a L indicates an isospin-1 spin-L resonance and K * L indicates an isospin-1/2 spin-L resonance. Conservation of G-parity for the η c → a L π case implies L to be even. For the decay η c → K * L K, G-parity relates the final states K * L K and K * Since G-parity is positive, it follows that in the Dalitz plot, the K * L K and K * L K bands will interfere constructively and all values of L are in principle allowed. In particular, the a 0 (980) resonance is described by a coupled-channel Breit-Wigner as measured in Ref. [41]. In the following, all resonances listed in Ref. [3] and satisfying the above criteria are included or tested in the Dalitz-plot analysis. The twobody invariant-mass resolution varies with invariant-mass and is approximately 7 MeV in the Kπ mass region between 1 and 2 GeV, much smaller than the width of the resonances present in this analysis. Therefore, resolution effects are ignored. The efficiency-corrected fractional contribution f i , due to resonant or non-resonant contribution i, is defined as follows: where x and y indicate the Dalitz-plot variables and the integrals in the numerator and denominator are computed on a large sample of simulated events generated uniformly in the decay phase space. The f i do not necessarily sum to 100% because of interference effects. The uncertainty for each f i is evaluated by propagating the full covariance matrix obtained from the fit. Interference fractions between amplitudes i and j are evaluated as: After the fit, a large simulated sample is prepared, where events are generated uniformly in the decay phase space [37]. These events are weighted by the fitted likelihood function, normalized to the data-event yields and compared to the data on several invariant-mass and angular projections. To assess the fit quality, the Dalitz plot is divided into a grid of n × n cells (n depends on the event yields and Dalitz-plot size) and only those cells containing at least two simulated events are considered. A χ 2 estimator is used, defined as where N i obs and N i exp are event yields from data and simulation, respectively. Here ε = N i exp for cells containing more than 9 entries, while it is approximated as the average of the lower and upper Poisson 68% CL errors for lower statistics cells. The figure of merit is defined as χ 2 /ndf where ndf is computed as N cells − n par and n par is the number of free parameters.
Two amplitude fits are then performed. The first employs a QMI description of the Kπ S-wave (Sec. 7.4) and the second uses an isobar model (Sec. 7.5).

Fits to the η c sidebands
The background model is obtained by fitting the η c sidebands. An unbinned maximumlikelihood fit is performed by inserting, one-by-one, all possible a 0 and a 2 resonances decaying to a K 0 S K state and all possible K * resonances decaying to a Kπ state [3] which could contribute to the K 0 S Kπ final state. To describe the fraction of background candidates associated to the B + decay (see Table 5), resonances are described by relativistic Breit-Wigner functions with appropriate Blatt-Weisskopf [40] form factors using angular terms just as for the η c decay but, for K * resonances, there is no symmetrization requirement for charged and neutral contributions. Interference is also allowed between all the contributing amplitudes. The fraction of background not associated to the B + decay is assumed to be uniform, except for a small contribution from K * (892) resonances, described only by relativistic Breit-Wigner functions. Small open-charm contributions are described by simple Gaussian functions with parameters fitted to the data. No efficiency is included in the description of the sidebands. The amplitudes are numerically integrated using a phase-space simulation generated according to the K 0 S Kπ invariant-mass distributions in the sidebands shown in Fig. 12. Contributions which are uniform in the phase space of the decay, having fixed fractions f B (given in Table 5), are included. Due to the limited sample size in the sideband regions, the K 0 SLL and K 0 SDD data are combined.

Fit using a QMI description of the Kπ S-wave
To measure the I = 1/2 Kπ S-wave, the QMI technique described in Refs. [11,42] is used. In this fit, the Kπ S-wave, being the largest contribution, is taken as the reference amplitude. The Kπ invariant-mass spectrum is divided into 37 equally-spaced mass intervals each 50 MeV wide, and two new fit parameters are added for each interval: the amplitude and the phase of the Kπ S-wave. The amplitude is fixed to 1.0 and its phase to π/2 at an arbitrary point in the mass spectrum, chosen to be interval 16, corresponding to a mass of 1.45 GeV. The number of free parameters associated to the description of the QMI Kπ S-wave is therefore 72. The Kπ S-wave amplitude in bin j is written as where the amplitude a Kπ j (m) = a K 0 π j (m) and the phase ϕ Kπ j (m) = ϕ K 0 π j (m) at the Kπ mass m. All the K * L (1430)K contributions with L = 0, 2 are symmetrized in the same way as for the S-wave amplitude.
The fit is initiated by performing a randomized scan for the Kπ S-wave solution with arbitrary start values of parameters. In addition, expected contributions from known resonances such as a 0 (980), a 2 (1320), K * 2 (1430), a 0 (1450), a 2 (1750) and K * 2 (1900) are included. The a 0 (1700) resonance, recently observed in the Dalitz-plot analysis of η c → ηπ + π − [14] is also included and found to be significant. The presence of an a 0 (1950) resonance, for which some evidence was found in the η c Dalitz-plot analysis Ref. [11], is tested, but its contribution is found to be consistent with zero. In the search for the best solution, the fit is iterated from the first found solution and additional resonances are added one-by-one; the process is completed when additional contributions give fit fractions consistent with zero. Spin-1 K * resonances (K * (892), K * (1410) and K * (1680)) are also tested but their contributions are found to be consistent with zero. The inclusion of additional spin-0 and spin-2 resonance components with unconstrained resonance parameters leads to no improvement of the fit quality. Figures 14 and 15 show the Dalitz-plot projections with superimposed projections of the fit function for B + → K 0 S K + K − π + and B + → K 0 S K + K + π − data, respectively. The distributions of the fitted background, interpolated from the η c and B + sidebands, are also included. The fractional contributions and relative phases from the fit are given in Table 6. There is a significant contribution from interference terms, evidenced by a sum of fractions much larger than 100%. This is a common feature of Dalitz-plot analyses dealing with several broad resonances having the same quantum numbers [43]. Interferences between amplitudes are listed in Table 7 for absolute values above 3%.
To test the fit quality, the efficiency-uncorrected Legendre-polynomial moments P L are computed in each K − π + , K 0 S π + , and K 0 S K − mass interval by weighting each event by a P L (cos θ) function, where θ is the corresponding helicity angle [44]. These distributions are shown in Figs. 26, 27 and 28 of Appendix A as functions of the K − π + , K 0 S π + , and K 0 S K − invariant masses, respectively. The extracted moment distributions are compared with the expected Legendre polynomial moment distributions obtained from simulated data weighted by the fit results. Good agreement for all the distributions is observed, which indicates that the fit is able to reproduce the local structures apparent in the Dalitz plot.
Systematic uncertainties, evaluated separately for the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − final states, are summarized in Tables 8 and 9, respectively. The effect of the efficiency correction, labeled as "Eff" in Table 8, is evaluated by replacing the two-dimensional fitted efficiency ϵ(m x (K 0 S K), cos θ π )) with the product of the fits to the efficiency projections ϵ = ϵ 1 (m x (K 0 S K)) · ϵ 2 (cos θ π ). The uncertainty related to trigger effects (Trig) is estimated by performing separate Dalitz-plot analyses of the TOS and no- Table 6: Results from the Dalitz-plot analysis of the η c decay in (top) B + → K 0 S K + K − π + , (center) B + → K 0 S K + K + π − and (bottom) inverse-variance-weighted averages. The QMI model is used for the Kπ S-wave.
motion not continuous with the rest of the distribution is present in the threshold region of Fig. 16(b). This behavior may be attributed to the inability of the algorithm to obtain a correct phase motion in this region of the phase space due to the absence, in this mass region, of additional significant resonant contributions in addition to that of the K * 0 (1430) resonance.
The inverse-variance-weighted average of the Kπ S-wave from the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − data is shown in terms of the real and imaginary parts of its amplitude in the Argand diagram of Fig. 17. The plot shows the expected anticlockwise phase motion due to the K * 0 (1430) resonance and a second loop due to the presence of the K * 0 (1950) resonance, observed here for the first time. Note that it was not possible to observe this behavior in Ref. [6] due to the presence of two solutions above the mass of 1.6 GeV.  Argand diagram for the Kπ S-wave averaged over the QMI results from Systematic uncertainties on the QMI Kπ S-wave amplitudes are evaluated as follows. The B + → K 0 S K + K − π + and B + → K 0 S K + K + π − datasets are each fitted independently using QMI for lower and higher purity selections, averaging the results; this is also done separately for TOS and noTOS selections. The four solutions are compared with the reference solution and the absolute values of the deviations are added in quadrature. The effect on the QMI Kπ S-wave due to the variation of the resonance parameters fixed to known values, in particular those for the a 0 (980) and K * 2 (1430) resonances, is tested and found to be negligible. The numerical values of the resulting QMI Kπ S-wave amplitudes and phases are listed in Appendix A, and this topic is further discussed in Sec. 7.5.

Isobar model η c Dalitz plot analysis
The QMI method used in the previous section has allowed information to be obtained about the Kπ S-wave which is produced by a combination of several scalar resonances. However, the procedure does not give insight into the parameters of the contributing resonances. An additional description of the data is therefore performed, known as the isobar model, which fits the data with a superposition of known resonances and additional (unknown or poorly known) states. Both methods are expected to give similar descriptions of the data.
The fitting method is the same as that described in Sec. 7.4, but in this case all the resonances are described by relativistic Breit-Wigner functions with parameters initially fixed to their default values [3]. The search for the best solution is performed by adding resonances one-by-one, using the K * 0 (1430) resonance as the reference amplitude, and considering as figures of merit, the likelihood and χ 2 /ndf behavior, as discussed in Sec. 7.4. It is found that, using only the list of the known resonances [3], a poor description of the data is obtained. A large improvement (increased likelihood ∆(2 log L) = 1338 and decreased χ 2 /ndf) is obtained by adding an additional scalar Kπ resonance with free parameters, labeled in the following as κ(2600) and included in the fit using the simple nonrelativistic BW function of Eq. (5). Resonances having poorly known parameters, namely the K * 0 (1430), K * 0 (1950) and a 0 (1700) resonances, are allowed to have free parameters. The significance of the K * 2 (1980) is found to be consistent with zero and its contribution is removed from the fit. The procedure is developed on the B + → K 0 S K + K − π + data, because of the better quality of this final state due to reduced combinatorial background, and is then tested on the B + → K 0 S K + K + π − final state. The fit projections are shown in Fig. 18 for the B + → K 0 S K + K − π + data; they are similar for the B + → K 0 S K + K + π − data and are shown in Fig. 29 of Appendix A.
The resulting fitted parameters are listed in Table 10 together with the estimated significances, evaluated using Wilks' theorem [45]. A test is made for the presence of the low-mass κ/K * 0 (700) resonance, with parameters fixed to average values [3]. The resulting likelihood variation is ∆(2 log L) = 6.6 for the difference of two parameters which corresponds to a significance of 1.8σ. The same is true for the B + → K 0 S K + K + π − decay where ∆(2 log L) = 12.6, corresponding to a significance of 2.9σ. Replacing the broad κ(2600) resonance with a uniform non-resonant contribution results in a decrease of the likelihood function by ∆(2 log L) = 254.
The fitted resonance contributions and relative phases are given in Table 11. Interferences between amplitudes are listed in Table 12 and shown in Fig. 18 for absolute values above 3%.
Systematic uncertainties on fractions and relative phases are evaluated as in Sec. 7.4 with the addition of an alternative model for the K * 0 (1430) resonance: a simplified coupledchannel Breit-Wigner function, which ignores the small Kη coupling. Its lineshape is parameterized as where m 0 is the resonance mass, g Kπ and g Kη ′ are the couplings to the Kπ and Kη ′ final states, and ρ j (m) = 2p/m are the respective Lorentz-invariant phase-space factors, with p the decay-particle momentum in the K * 0 (1430) rest frame. The ρ 2 (m) function becomes imaginary below the Kη ′ threshold. The Dalitz-plot analysis is insensitive to the g Kη ′ parameter, which has been fixed to g 2 Kη ′ = 0.197 GeV 2 /c 4 [14]. A similar method is used for the evaluation of the systematic uncertainties on the resonance parameters.
The isobar model allows the separation of the resonance composition of the Kπ S-wave  Figure 18: Fit projections on the K − π + , K 0 S π + and K 0 S K − invariant-mass distributions from the Dalitz-plot analysis of the η c decay using the isobar model for the B + → K 0 S K + K − π + data. Panels (a, c, e) show the most important resonant contributions. To simplify the plots, only resonant contributions relative to that mass projection are shown. Panels (b, d, f) show the interference terms for contributions greater than 3%. The legend in (a) also applies to (c) and the legend in (b) also applies to (d) and (f).
obtained from the QMI method. A fit to the Kπ S-wave amplitude and phase is performed using the sum of three interfering spin-0 relativistic Breit-Wigner amplitudes: where m indicates the Kπ mass. In this fit, the mass and width of the three resonances Table 11: Results from the Dalitz-plot analysis of the η c using the isobar model for (top) B + → K 0 S K + K − π + decays, B + → K 0 S K + K + π − decays (center) and their inverse-varianceweighted averages (bottom).  are fixed to the results obtained in the isobar model analysis while their coefficients and relative phases are left as free parameters. The fitted Kπ S-wave used in the fit is obtained from the inverse-variance-weighted averages from the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − data. The QMI amplitude and phase, together with the squared amplitude, are shown in Fig. 19, where the statistical and systematic uncertainties are added in quadrature and the phase measurements in the Kπ mass region below 0.90 GeV Table 12: Fractional interference contributions from the Dalitz-plot analysis of the η c decay in B + → K 0 S K + K − π + decays using the isobar model. Absolute values less than 3% are not listed.
is described. Figure 5 shows the K 0 S Kπ mass spectra for the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − data. The η c (2S) signal is accompanied by a significant background contribution which is especially large in the B + → K 0 S K + K + π − data by virtue of the two indistinguishable particles in the final state. Therefore, only the B + → K 0 S K + K − π + data are used for the η c (2S) Dalitz-plot analysis. In addition to the background from D 0 → K 0 S K + K − (see Sec. 3), additional backgrounds involving the spectator K + are removed. A significant contribution from ϕ(1020) → K + K − decays is removed by selecting events outside the 1.01 < m(K + K − ) < 1.03 GeV mass region, and a small contribution from D + s → K + K − π + decays is removed in the ±2σ (σ = 7.8 MeV) mass region close to the known D + s mass. Finally, a significant contribution from D − s → K 0 S K − decays is removed within ±2.5σ (σ = 7.8 MeV) from the known D + s mass [3].    Table 13 shows the η c (2S) event yields in the signal region, purities and Table 13: Candidate events in the signal region, purities, and B + background contributions to the η c (2S) → K 0 S Kπ Dalitz-plot analysis separated for K 0 S types. The incoherent background fractions f B in the low and high sidebands refer to the sum of the K 0 SLL and K 0 SDD data.   background compositions. Figure 21 shows the η c (2S) Dalitz plot. The distribution is dominated by the K * 0 (1430) resonance with bands due to the K * (892) resonance originating from background.
The decay of the η c (2S) to K 0 S Kπ is expected to be very similar to that of the η c , except for the different size of the available phase space. Therefore, the Dalitz-plot analysis follows closely the method used for the η c analysis, except that, due to the limited statistics and the significant background, only the isobar model is used. In addition, all resonance parameters are fixed to those extracted from the η c isobar model. The efficiency model used for the η c (2S) Dalitz-plot analysis has been described in Sec. 6.2.1. The resulting fitted amplitudes and phases are listed in Table 15. The Dalitz-plot projections, together with the fitted background interpolated from the sidebands and the superimposed fit, are shown in Fig. 22. Large interference effects are observed, evidenced by the high value of the sum of the fractions. Systematic uncertainties are evaluated in a similar way as for the η c Dalitz-plot analysis. In particular, the effect of the efficiency model is evaluated by replacing the fitted efficiencies by two-dimensional binned numerical maps. The effects of the uncertainties on the resonance parameters are evaluated by varying masses and widths by their statistical uncertainties according to a Gaussian distribution and averaging the absolute values of the deviations from the reference fit. The effects of the different sources of systematic uncertainties on the fitted fractions and phases are listed in Table 14.
The K 0 S Kπ invariant-mass spectrum in the χ c1 -η c (2S) mass region, summed over the K 0 SLL and K 0 SDD data, is shown in Fig. 23 for the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − final states. The signal region is [3.48-3.54] GeV and the lower and In light gray are indicated the χ c1 signal region and sidebands. In dark gray is indicated the incoherent background estimated from the B + sidebands. Table 17: Candidate events and purities of χ c1 for the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − final states separated for the K 0 SLL and K 0 SDD data. The incoherent background fractions f B values for the lower and higher sideband regions are computed for the sum of the K 0 SLL and K 0 SDD data.  Table 17 reports the χ c1 event yields, purities and background compositions, separated for the K 0 SLL and K 0 SDD data, for the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − decays. The χ c1 Dalitz plot for B + → K 0 S K + K − π + data, shown in Fig. 24, is dominated by horizontal and vertical bands due to the presence of K * (892) and K * 2 (1430) resonances. A Dalitz-plot analysis of the χ c1 decay to the K 0 S Kπ state is not feasible with the present dataset due to its limited sample size and the high level of background; therefore a simplified approach, similar to that used in Ref [16], is taken where fits are applied to the Kπ invariant-mass projections. The efficiency distribution as a function of the Kπ and K 0 S π invariant masses is computed using the method described in Sec. 6.2.1 and expressed in terms of m x (K 0 S K) and cos θ π . Each event in the χ c1 mass region is weighted by the inverse of the efficiency. The ratios of the unweighted over weighted Kπ and K 0 S π invariant-mass distributions are then computed, yielding the efficiency as functions of the Kπ and K 0 S π invariant masses. These efficiencies are consistent with being uniform, so no efficiency correction is used in fitting to the Kπ or K 0 S π mass projections. The background contribution is evaluated using the sidebands shown in Fig. 23. For each K 0 S category, the Kπ and K 0 S π invariant-mass distributions from the χ c1 signal region and from the sidebands are obtained. The Kπ and K 0 S π invariant-mass distributions in the sidebands are normalized to the χ c1 fitted purity, and subtracted from the invariant-mass spectra in the χ c1 signal region, obtaining the distributions shown in Fig. 25. Prominent K * (892) and K * 2 (1430) resonances are observed. A fit to the Kπ and K 0 S π invariant-mass distributions is performed using an empirical and flexible background shape in addition to two relativistic Breit-Wigner functions describing the K * (892) and K * 2 (1430) resonances with parameters fixed to their known values [3]. The background is parameterized as B(m) = P (m)e a 1 m+a 2 m 2 for m < m 0 , and B(m) = P (m)e b 0 +b 1 m+b 2 m 2 for m > m 0 , where P (m) is the two-body phase space and m 0 , a 1 , a 2 , b 0 , b 1 , b 2 are free parameters. The two functions and their first derivatives are required to be continuous at m = m 0 , and this constraint reduces the number of freely varying parameters to four. Figure 25 shows the fits to the background-subtracted Kπ and K 0 S π invariant-mass distributions in the χ c1 signal region for both the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − data. Table 18 reports candidate yields and fractional contributions obtained from the fits.
The inverse-variance-weighted averages of the fractional K * (892) and K * 2 (1430) contributions obtained from the fits to the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − data are listed in Table 19. The fractional K * (892) and K * 2 (1430) contributions are converted into χ c1 → K * K branching fractions (listed in Table 19) by multiplying the fractional contributions by the known χ c1 branching fraction [3], The results of the fits are overlaid. and correcting for unseen K * decay modes. A comparison with Ref. [16] shows good agreement, while improving the statistical uncertainties and significances. Systematic uncertainties are evaluated as follows. The deviations using different BDT classifier selections are evaluated by fitting the lower-and higher-purity datasets and averaging the absolute values of the deviations of the fitted fractions in comparing to the values obtained from the default fit. The differences resulting from the use of the Table 19: Inverse-variance-weighted averages of the fractional contributions for χ c1 decays to K * K resonances from the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − final states. The branching fractions of the χ c1 → K 0 S Kπ decays to K * (892) and the K * 2 (1430) resonances are also included. The reported uncertainties are statistical, systematic and from the uncertainty on the χ c1 branching fraction (see text) [3].
The entire dataset, selected as described in Sec. 3 is used in this section. Efficiencycorrected invariant-mass distributions are obtained by weighting each event by the inverse of the total efficiency described in Sec. 6. Fits are performed to the resulting efficiencycorrected invariant-mass spectra, and the extracted yields of B + decays and of various charmonium resonances are used to obtain efficiency-corrected ratios with respect to the η c and J/ψ resonances. The procedure is performed separately for K 0 SLL and K 0 SDD data, which are subsequently combined by means of inverse-variance-weighted averages.
Prior to evaluating the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − candidate yields, it is necessary to subtract all contributions from open-charm production. These are observed in the two-or three-body invariant-mass spectra of the B → D (s) X decays, where X represents one or two additional particles. The charmed-meson vetoes from the previous sections are not applied. Since open-charm resonances may also be present in the B + -meson background, a background subtraction is performed using the normalized B +candidate invariant-mass sidebands. This is achieved by plotting the different two-body or three-body invariant-mass combinations in each of the two sideband regions, with both halves as wide as the B + signal region. The resulting invariant-mass distributions are Table 20: Resulting yields and total charm fraction from the fits to the four-body, three-body and two-body invariant-mass spectra uncorrected for efficiency for the (top) B + → K 0 S K + K − π + and (bottom) B + → K 0 S K + K + π − data, separated by K 0 S category. The quoted uncertainties are statistical only.

Resonance
LL DD  then subtracted from the corresponding ones in the signal region. In performing fits to the two-or three-body mass spectra, the charmed mesons are described by two Gaussian functions, sharing the same mean, with first-or second-order polynomials used to model the background. Table 20 summarizes the fit results for the B + → K 0 S K + K − π + and B + → K 0 S K + K + π − data, without efficiency corrections. The efficiency-corrected K 0 S Kπ and K + K − invariant mass spectra are fitted as described in Sec. 5 to obtain efficiency-corrected yields for η c , J/ψ (see Fig. 4), χ c1 , η c (2S) (see Fig. 5), χ c0 and χ c2 (see Fig. 6) resonances. The efficiency corrected K 0 S K + K − π + and K 0 S K + K + π − invariant-mass spectra are fitted as described in Sec. 3 (see Fig. 2). These yields are used to compute the ratios with respect to the η c and J/ψ resonance yields, labeled as and respectively, where X labels the charmonium resonance. Similar expressions are used to calculate the B + four-body decays ratios R 1 (B + ) and R 2 (B + ). Table 21 lists the measurements of the efficiency-corrected relative ratios.
The following systematic uncertainties are considered. In fitting for B + -candidate and charmonium-resonance yields, the fit model is varied to use Gaussian functions sharing or not sharing the same mean on a quadratic, linear or cubic polynomial background shape. The maximum deviation from the default fit is taken as the systematic uncertainty for each parameter. The open-charm fraction is re-evaluated for efficiency-corrected invariant-mass Resonance 4.14 ± 0.04 ± 0.25 3.40 ± 0.04 ± 0.09 R 1 (J/ψ) 0.201 ± 0.005 ± 0.008 0.200 ± 0.005 ± 0.013 0.082 ± 0.004 ± 0.003 0.073 ± 0.004 ± 0.007 R 1 (η c (2S)) 0.113 ± 0.005 ± 0.004 0.109 ± 0.007 ± 0.016 4.95 ± 0.12 ± 0.23 5.26 ± 0.24 ± 0.29 0.56 ± 0.03 ± 0.02 0.55 ± 0.04 ± 0.08 0.06 ± 0.01 ± 0.01 distributions and the difference with respect to the default fit is considered as a systematic uncertainty. In fitting the K 0 S Kπ invariant-mass spectra, the bin size is changed and the background shape is modified from a linear to quadratic function. To evaluate the size of the fit bias on the resonance yields, starting from the reference fits, 400 randomized pseudoexperiments are generated according to Poissionian statistics, having the same sample size as the original ones, and then fitted. The resonance yields are compared with the reference fits. In all cases, the deviations from the reference fits are included as systematic uncertainties. The entire procedure of fitting the invariant-mass spectra and evaluating the relative ratios is repeated dividing the data into the different trigger conditions, TOS and noTOS, whose results are then averaged and compared with the reference fit. All the systematic uncertainties are added in quadrature to calculate the total systematic uncertainty.
Using the listed ratios it is possible to evaluate the branching fractions as for the η c mode and for the J/ψ mode, where the f i are correction factors for unseen decay modes, including K 0 modes. The resulting measured branching fractions are listed in Table 22 (B 1 ) and Table 23 (B 2 ). For the modes involving χ c0 and χ c2 resonances, listed in Table 24, the f i parameters also include the corrections for their known branching fractions to the K + K − state, (6.05 ± 0.31) × 10 −3 and (1.01 ± 0.06) × 10 −3 , respectively [3]. The B + → K 0 K + K − π + and B + → K 0 K + K + π − branching fractions are measured here for the first time. Tension is found for the B + → χ c2 K 0 π + branching fraction with respect to the PDG value [3], measured only in Ref. [47] to be (0.116±0.022±0.012)×10 −3 using the χ c2 → J/ψγ decay mode as a reference with 76.4 ± 14.7 events and a significance of 4.6σ. This measurement differs from that reported here by 2.6σ (3.2σ), using the η c Table 22: Measured branching fractions with the η c resonance as a reference using the (top) B + → K 0 K + K − π + and (bottom) B + → K 0 K + K + π − data. The first uncertainty is statistical, the second systematic and the third due to the PDG uncertainty on the B + → η c K + branching fraction.

Summary
A study is presented of B + → K 0 S K + K − π + and B + → K 0 S K + K + π − decays at protonproton collision energies of 7, 8 and 13 TeV using the LHCb detector. The K 0 S Kπ invariant-mass spectra from both B + decay modes show a rich spectrum of charmonium resonances. New measurements of the η c and η c (2S) masses and widths are performed resulting in the best determinations from a single measurement. Branching fractions are measured for B + decays to η c , J/ψ, η c (2S) and χ c1 resonances. Evidence is also found for the B + → χ c2 K 0 π + decay. Furthermore, the first observation and branching fraction measurement of B + → χ c0 K 0 π + decays is reported. The first measurements of B + → K 0 K + K − π + and B + → K 0 K + K + π − branching fractions are also given. A Dalitzplot analysis of η c → K 0 S Kπ decay is undertaken using both a quasi-model-independent approach for the Kπ S-wave and an isobar-model approach. The K * 0 (1950) → Kπ resonance is established and its parameters are measured. A new parametrization of the Kπ S-wave is obtained, which includes K * 0 (1430), K * 0 (1950) and broad κ(2600) resonances. No evidence for the κ(700) resonance is found in η c decays. The a 0 (1700) resonance as well as its decay to K 0K are confirmed, and its parameters are measured. The first Dalitz-plot analysis of the η c (2S) → K 0 S Kπ decay is performed. Finally, branching fractions of χ c1 → K * (892) 0K 0 , χ c1 → K * 2 (1430) 0K 0 , χ c1 → K * (892) + K − and χ c1 → K * 2 (1430) + K − decays are measured; these measurements are more precise than previous evaluations and correspond to the first observations of the two χ c1 → K * 2 (1430)K decay modes.

A Appendix
The inverse-variance-weighted averages of the numerical values of the Kπ S-wave as a function of the Kπ mass from the Dalitz plot analysis of η c mesons in B + → K 0 S K + K − π + and B + → K 0 S K + K + π − decays are listed in Table 25.   Figure 26: Projected K − π + invariant mass spectra weighted by Legendre polynomial moments P L (cos θ K 0 S ) for η c → K 0 S K − π + from B + → K 0 S K + K − π + data. The drawn curves result from the Dalitz plot fit using the QMI method described in the text. Figures 26,27, and 28 show the K − π + , K 0 S π + and K 0 S K − invariant mass distributions weighted by Legendre polynomial moments computed as functions of cos θ K 0 S , cos θ K − and cos θ π + , respectively, from the Dalitz plot analysis of the η c decay to K 0 S Kπ using the QMI approach to the fit to B + → K 0 S K + K − π + data. Figure 29 shows the mass projections from the η c Dalitz plot analysis using the QMI model from B + → K 0 S K + K + π − data.  Figure 27: Projected K 0 S π + invariant mass spectra weighted by Legendre polynomial moments P L (cos θ K − ) for η c → K 0 S K − π + from B + → K 0 S K + K − π + data. The drawn curves result from the Dalitz plot fit using the QMI method described in the text.  Figure 28: Projected K 0 S K − invariant mass spectra weighted by Legendre polynomial moments P L (cos θ π + ) for η c → K 0 S K − π + from B + → K 0 S K + K − π + data. The drawn curves result from the Dalitz plot fit using the QMI method described in the text.