Study of B^{+-} ->K^{+-}(K_S K pi)^0 Decay and Determination of eta_c and eta_c(2S) Parameters

We report the results of a study of $B^{\pm}\to K^{\pm}\eta_c$ and $B^{\pm}\to K^{\pm}\eta_c(2S)$ decays followed by $\eta_c$ and $\eta_c(2S)$ decays to $(K_SK\pi)^0$. The results are obtained from a data sample containing 535 million $B\bar{B}$-meson pairs collected by the Belle experiment at the KEKB $e^+e^-$ collider. We measure the products of the branching fractions ${\mathcal B}(B^{\pm}\to K^{\pm}\eta_c){\mathcal B}(\eta_c\to K_S K^{\pm}\pi^{\mp})=(26.7\pm 1.4(stat)^{+2.9}_{-2.6}(syst)\pm 4.9(model))\times 10^{-6}$ and ${\mathcal B}(B^{\pm}\to K^{\pm}\eta_c(2S)){\mathcal B}(\eta_c(2S)\to K_S K^{\pm}\pi^{\mp})=(3.4^{+2.2}_{-1.5}(stat+model)^{+0.5}_{-0.4} syst))\times 10^{-6}$. Interference with the non-resonant component leads to significant model uncertainty in the measurement of these product branching fractions. Our analysis accounts for this interference and allows the model uncertainty to be reduced. We also obtain the following charmonia masses and widths: $M(\eta_c)=(2985.4\pm 1.5(stat)^{+0.5}_{-2.0}(syst))$ MeV/$c^2$, $\Gamma(\eta_c)=(35.1\pm 3.1(stat)^{+1.0}_{-1.6}(syst))$ MeV/$c^2$, $M(\eta_c(2S))=(3636.1^{+3.9}_{-4.2}(stat+model)^{+0.7}_{-2.0}(syst))$ MeV/$c^2$, $\Gamma(\eta_c(2S))=(6.6^{+8.4}_{-5.1}(stat+model)^{+2.6}_{-0.9}(syst))$ MeV/$c^2$.

10 −6 . Interference with the non-resonant component leads to significant model uncertainty in the measurement of these product branching fractions. Our analysis accounts for this interference and allows the model uncertainty to be reduced. We also obtain the following charmonia masses and widths: M (η c ) = (2985.4 ± 1.5(stat) +0. 5 −2.0 (syst)) MeV/c 2 , Γ(η c ) = (35.1±3.1(stat) +1.0 −1.6 (syst)) MeV/c 2 , M (η c (2S)) = (3636.1 +3. 9 −4.2 (stat+model) +0.7 −2.0 (syst)) MeV/c 2 , Γ(η c (2S)) = (6.6 +8.4 −5.1 (stat+model) +2. 6 1 Introduction Charmonium states consist of a heavy charm-anticharm quark pair, which allows the prediction of some of the parameters of these states using nonrelativistic and relativistic potential models [1], lattice QCD [2], non-relativistic effective field theory (NRQCD) [3], and sum rules [4] (see the recent review in [5]). The comparison of these predictions with experimental results provides an opportunity to tune the parameters of theoretical models and, therefore, improve the accuracy of other values predicted by these models. We have to measure the charmonium masses, widths, and product branching fractions with enough accuracy to compare them with theoretical predictions. Parameters of (cc) states such as the η c and η c (2S) mesons have been studied in various experiments using a variety of decay channels [6]. Tables 1 and 2, which list a selection of most precise mass and width determinations, show that there is quite a large spread of measured masses and widths of the η c and, especially, η c (2S) mesons resulting in large scale factors for the world average values [6]. Moreover, our knowledge of hadronic decays of these charmonia is rather poor.
In the Belle and BaBar B factory experiments, charmonia are produced in various ways: from fragmentation in electron-positron annihilation, from twophoton processes, and in B decays. The advantages of B ± → K ± cc decay are the relatively large reconstruction efficiency, small background, and the fixed quantum numbers (J P = 0 − ) of the initial state. Here we consider the following decays of charged B mesons: A 492 fb −1 data sample provides an opportunity to determine the corresponding products of the branching fractions as well as masses and widths of the η c and η c (2S) mesons. Table 1 Previously measured η c parameters.

Event selection
The results are based on a data sample that contains 535 × 10 6 BB pairs, collected with the Belle detector at the KEKB asymmetric-energy e + e − collider [18] operating at the Υ(4S) resonance.
The Belle detector [19] is a large-solid-angle magnetic spectrometer that consists of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC) for charged particle tracking and specific ionization measurement (dE/dx), an array of aerogel threshold Cherenkov counters (ACC), time-of-flight scintillation counters (TOF), and an array of 8736 CsI(Tl) crystals for electromag-netic calorimetry (ECL) located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux return located outside the coil is instrumented to detect K 0 L mesons and identify muons (KLM). We use a GEANT-based Monte Carlo (MC) simulation to model the response of the detector and determine its acceptance [20].
Pions and kaons are separated by combining the responses of the ACC and the TOF with dE/dx measurements in the CDC to form a likelihood L(h) where h = π or K. Charged particles are identified as pions or kaons using the likelihood ratio R: Charged tracks are selected with requirements based on the χ 2 of the track fits and the impact parameters relative to the interaction point. We require that the polar angle of each track be in the angular range 18 • -152 • and that the track momentum perpendicular to the positron beamline be greater than 100 MeV/c.
Charged kaon candidates are identified by the requirement R(K) > 0.6, which has an efficiency of 90% and a pion misidentification probability of 3 -10% depending on momentum. For pion candidates we require R(π) > 0.2. K S candidates are reconstructed via the π + π − mode. We apply the following cut on the π + π − invariant mass: 0.489 GeV/c 2 < M(π + π − ) < 0.505 GeV/c 2 . The flight length of the K S is required to lie within the interval [0.1 : 20] mm. The condition on the K S direction angle ϕ is cos(ϕ) > 0.95.
B meson candidates are identified by their center-of-mass (c.m.) energy differ- √ s/2 is the beam energy in the Υ(4S) c.m. frame, and p i and E i are the c.m. three-momenta and energies, respectively, of the B meson candidate decay products. The signal region is defined as: |∆E| < 0.03 GeV, 5.273 GeV/c 2 < M bc < 5.285 GeV/c 2 . The ∆E sideband region is defined as ||∆E| − 0.06(GeV)| < 0.03 GeV.
To suppress the large continuum background (e + e − → qq, where q = u, d, s, c), topological variables are used. Since the produced B mesons are nearly at rest in the c.m. frame, the signal tends to be isotropic while continuum qq background tends to have a two-jet structure. We use the angle between the thrust axis of the B candidate and that of the rest of the event (Θ thrust ) to discriminate between these two cases. The distribution of | cos Θ thrust | is strongly peaked near | cos Θ thrust | = 1 for qq events and is nearly uniform for Υ(4S) → BB events. We require | cos Θ thrust | < 0.8.
If there is more than one combination that satisfies the selection criteria, we select the B candidate with the minimum difference |M(K S ) − M(π + π − )| and the minimum difference between the vertex z-coordinates of the kaon and pion from the charmonium decay, and the kaon from the B decay. If the final state includes two kaons of the same charge, as in K ± (K ± K S π ∓ ), we must choose the one from the charmonium decay. For this purpose we select the candidate with the minimum difference |M(η c /η c (2S)) − M(K S K ± π ∓ )|. The mean number of multiple candidates per event is 1.6.

Interference study
The data sample can contain signal, which has the same final state as a resonant decay but does not include a charmonium resonance. The contribution of these events is referred to as the non-resonant amplitude. Since the final state is the same, this amplitude interferes with the signal. However, if the final particles form narrow resonances such as D, D S , and φ mesons, the interference effect cancels after integration over mass. Thus, we can reject some background channels by applying appropriate cuts on the mass combinations of the final state particles. If the intermediate resonances have a substantial width or the B meson decays directly into the final state particles, the effect of interference must be taken into account.
The non-resonant contribution can be seen as a peak in the ∆E distribution in the charmonium sideband regions shown in Figs. 1 and 2 (plots on the right) 1 . We fit the ∆E distributions with the sum of a Gaussian distribution and a second-order polynomial function. From these fits we obtain the number of events in the signal region N obs and in the sideband region N sb , which can be rescaled to obtain the number of non-resonant events N non−res . In the η c case N obs = 889 ± 37(stat) and N non−res = 87 ± 11(stat). In the η c (2S) case N obs = 279 ± 29(stat) and N non−res = 156 ± 13(stat).
Different values of the interference phase can give significant variations in the number of signal events while the total number of observed events remains the same. Using the hypotheses of maximal constructive and destructive interference, we would obtain 410 and 1550 signal events, respectively. Therefore, the model uncertainty in the number of η c signal events is rather large: N signal = 980 ± 570(model). It would be even larger for the η c (2S) decay. A dedicated study of the interference effect allows this uncertainty to be reduced.
In the B → K (1) K S K (2) π decay (see Fig. 3) there are four particles in the final state, which gives 4 × 3 measured parameters. Taking into account the four constraints of energy-momentum conservation and integrating over the three angles that characterize the B decay (it is a pseudoscalar and there should be no dependence on these angles), we have 5 independent variables to describe the amplitude of the process. We chose the following parameters: K S K (2) π invariant mass, two Dalitz variables for the η c (or η c (2S)) decay -q 2 1 and q 2 2 (for example, M(K (2) π) 2 and M(K S π) 2 ), the angle between the K S and K (1) in the rest frame of the K (2) K S π system (θ), and the angle between the planes of the K (1) − π and K (1) − K S in the same system (φ).
The M(K S Kπ) distribution has four peaks corresponding to η c , J/ψ, χ c1 and η c (2S) production (see Fig. 4). In addition to these peaks, there is a non-resonant signal, which interferes with the η c (or η c (2S)) signal 2 . Unfortunately, the shape of the one dimensional (1-D) mass distribution alone does not allow the interference contribution to be obtained, so other variables should be used.  Fig. 4. The signal distribution of (K S K ± π ∓ ) invariant mass in the B ± → K ± (K S Kπ) 0 decay. The charmonium states η c , J/ψ, χ c1 , and η c (2S) (in order of mass) can be seen. The solid histogram is the combinatorial background determined from the ∆E sideband region.
The Dalitz plots of 3-body η c and η c (2S) decays are shown in Fig. 5. In the η c case the distribution is not uniform and has a peaking structure around a Kπ mass of 1.4 GeV/c 2 , which can be a combination of the K * (1410), K * 0 (1430), and K * 2 (1430) states. However, the statistics are rather low, so it is impossible to determine with acceptable accuracy either the mass and width of these states, or the product branching fractions of η c decay modes involving these states. The small number of events does not allow the Dalitz analysis to be efficiently performed and makes it difficult to use the Dalitz plot variables to distinguish the η c signal and non-resonant amplitudes. The same conclusion holds for the η c (2S) decay.
Another variable that can be used for the amplitude separation is cos θ. Since the η c and η c (2S) are pseudoscalars (J P = 0 − ), we expect a uniform distribution in cos θ. In Figs. 6 and 7 the cos θ distributions are shown for the η c and η c (2S) signal and sideband regions, while the combinatorial background is subtracted. One can see that the sideband distribution has contributions from higher angular waves. A good fit can be obtained with the sum of S-, P-, and D-waves. The signal region also contains non-resonant background but mostly consists of signal events, so the S-wave contribution here is dominant. Separation of the P-and D-waves from the S-wave in the non-resonant background allows the uncertainty from interference to be reduced. Thus, we analyze a 2-D M(KK S π)-cos θ histogram assuming that the nonresonant signal amplitude is constant within the (2.5 -3.46) and (3.14 -4.06) GeV/c 2 mass ranges. The number of events in a single bin is N bin = N ∆Esignal − k · N ∆Esideband , where the coefficient k is used to normalize the number of events in the ∆E sideband region. We minimize the likelihood function, assuming that the events in the signal and sideband regions are described by the Poisson statistics. In the η c analysis the bin size along the cos θ axis is 0. We exclude the J/ψ region (3.07 -3.13 GeV/c 2 ) from the fit because the interference of the J/ψ and non-resonant signal is negligible due to the small width of the former and inclusion of this region does not constrain the η c interference contribution. Moreover, the J/ψ angular distribution has contributions from several amplitudes that are not well determined. The same arguments apply to the exclusion of the χ c1 mass region (3.48 -3.54 GeV/c 2 ). We perform separate fits to the J/ψ and χ c1 with Gaussian functions using 1-D K S Kπ invariant mass distributions from signal MC and data. After comparing the obtained widths we determine the degradation of the resolution in data. Taking this into account, we recalculate the detector resolution obtained from signal MC in the η c and η c (2S) regions. We obtain σ(η c ) = (6.2 ± 1.1) MeV and σ(η c (2S)) = (9.8 ± 1.7) MeV.
The fitting function can be represented as the square of the absolute value of the sum of the signal and non-resonant amplitudes integrated over all variables except M(K S Kπ) and cos θ: where x = cos θ, s = M 2 (K S Kπ); q 2 1 and q 2 2 are Dalitz plot variables; ε 1 and ε 2 are constants that characterize the efficiency dependence on x and are determined from MC; δ and ∆ are the bin widths in cos θ and M(K S Kπ) invariant mass, respectively; M and Γ are mass and width of the η c (η c (2S)) meson; N is the η c (η c (2S)) signal yield; α, β, γ are the relative fractions of the S-, P-, and D-waves, respectively; S = 1 √ 2 , P = 3 2 x, D = 3 2 5 2 (x 2 − 1 3 ) are the functions characterizing the angular dependence of the S-, P-, and D-waves, respectively; A η is the signal S-wave amplitude, A S,P,D are the background S-, P-, and D-wave amplitudes, respectively. The absolute values of the amplitudes squared are normalized to unity: To account for the momentum resolution, Eq. (1) is convolved with a Gaussian detector resolution function that is determined from the MC and calibrated from the J/ψ (χ c1 ) width in data.
Since the function F (s, x) is a sum of the η c (η c (2S)) Breit-Wigner and S-, P-, and D-waves, it can be represented as a rational function of s and x: In its most general form, such a function has 15 independent coefficients (C ij ) in the numerator and two (M and Γ) in the denominator, however, in our case some coefficients are not independent: Thus, we have (15 + 2 − 4 = 13) independent terms only, which is not enough to determine all 15 parameters of the function F (s, x). Two (15 − 13 = 2) of the parameters must be either obtained from other measurements or allowed to vary over the full allowed range. Since we have no additional information on these parameters, we scan over them. The result does not depend on the choice of the two scanned parameters. Since the interference is more significant in the S-wave, we choose α and ℑ ηS . To perform this scan, we randomly sample α and ℑ ηS in a reasonable range 3 with the remaining 13 parameters free. After fitting the distributions we obtain a set of parameters and a χ 2 . No additional local minima are found.
The dependences of the signal yields on χ 2 are shown in Fig. 8. In the η c case, one can see that this distribution has a "plateau", which consists of fits with different N signal (and other fit parameters) and the same χ 2 . This feature arises because our system of equations for the fit parameters is underdetermined. The variation of parameters within this plateau will be referred to as the model uncertainty of our analysis and the average of their statistical errors as the statistical uncertainty. In the η c (2S) case, the minimum χ 2 plateau is not reached, because the parameters ℜ ηS,ηP,ηD , ℑ ηS,ηP,ηD , and Π SP,SD,P D tend to their bounds 4 . These bounds make our system of equations fully determined, so we can float the parameters α and ℑ ηS . In that case, the model and statistical errors cannot be separated. Thus we obtain N signal = 920±50(stat)±170(model) for the η c decay and N signal = 128 +83 −58 (stat+model) for the η c (2S) decay. Projections of the fits using the function F (s, x) are shown in Figs. 9 and 10.
The fit procedure described above was also applied to MC signal samples. The obtained number of signal events was used to determine the detection efficiency ((9.32 ± 0.10)% and (10.18 ± 0.10)% for η c and η c (2S) decays, respectively) and hence to calculate the product branching fractions. The decay B ± → K ± (K S Kπ) 0 has two possible final states: K ± K S K ± π ∓ and K ± K S K ∓ π ± . We assume that the η c decay signal amplitudes are the same, but the non-resonant contributions (α, β, γ) could be different in each decay channel. Because of limited statistics we do not treat these final states separately. A single distribution (1) effectively describes the incoherent sum of two distributions for the two decay channels. This does not affect the parameters of η c (η c (2S)) states, but the parameters describing the non-resonant part obtained in the fit take effective values that depend on both non-resonant amplitudes.
The method described above was checked using toy MC, which showed that the described procedure gives parameter values consistent with the generated ones. Moreover, a generic MC test was performed that included a full simulation of all b → c decays at the Υ(4S) without interference. This test verified that the signal determination is not biased and gave an interference value consistent with zero.

Systematic uncertainties
We evaluated possible sources of systematic uncertainties in the product branching fractions. The number of BB pairs is calculated from the difference of the number of hadronic events on resonance and the scaled number of those offresonance. The systematic error is dominated by the uncertainty in the scale factor and is equal to ∼1.3%. We assume that the combinatorial background can be parameterized with a first-order polynomial. To obtain the background shape uncertainty, we describe the background by a second-order polynomial and compare the results. The uncertainty on the K S decay branching fraction is taken from [6]. The contribution of the K S reconstruction uncertainty was estimated in the Belle experiment to be 4.4% [21]. In our fitting procedure we take into account the efficiency dependence on cos θ and assume that it does not depend on the K S Kπ invariant mass. By adding a linear dependence on M(K S Kπ) we estimate the corresponding systematic error. Moreover, we take into account the dependence of the efficiency on the η c and η c (2S) decay models, such as KK * , KK * 0 (1430), and KK * 2 (1430). The corresponding contribution to the systematic uncertainty is estimated by varying the efficiency obtained using these models and taking the difference in the results. An analysis of the charged track reconstruction uncertainty as a function of particle momenta has been performed in Belle data and gave an estimate of 1% per charged track. To determine the errors due to K and π meson identification, data from analysis of the process D * + → D 0 π + followed by the decay D 0 → K − π + were used. The uncertainty in K ± identification is 0.8% per K meson and the corresponding value for π ± identification is 0.5% per π meson. We also take into account the deviation of MC from the data by applying a correction to the efficiency: ε Data ε M C is 0.9996 for each kaon and 0.9756 for each pion. We vary the bin size along the x axis from 0.15 to 0.225, the ∆E window from 20 MeV to 40 MeV, and the detector resolution within the limits of its statistical uncertainty.
Sources of systematic uncertainties for masses and widths include the background parameterization, bin size, and detector resolution as described above, as well as a scale uncertainty and the effect of identical kaons in the final state. The scale uncertainty is determined from a comparison of masses of the J/ψ (χ c1 ) resonances, which were obtained by fitting the K S Kπ invariant mass distribution, with the world average values [6]. In case two kaons of the same charge are present in the final state, the charmonium amplitude is a sum of two amplitudes corresponding to K S K (1) π and K S K (2) π combinations. This can lead to the deformation of the K S Kπ invariant mass distribution in the region of phase space where the K S K (1) π and K S K (2) π invariant mass values overlap. We use toy MC to estimate this effect and take it into account as an additional systematic error.
All the contributions to the systematic uncertainties are listed in Table 3 for the product branching fractions and in Table 4 for the masses and widths. Table 3 Systematic uncertainties of the product branching fractions (in %).  Table 4 Systematic uncertainties of masses and widths of the η c and η c (2S) mesons (in MeV/c 2 ).  Table 5 shows a comparison of the results obtained assuming no interference (1-D fits to the ∆E and to the K S Kπ invariant mass distributions) and those obtained using the analysis described above. One can see that taking interference into account leads to the introduction of a model error for the product branching fractions B(B ± → K ± η c )B(η c → K S K ± π ∓ ) (for the η c mass and width this error turns out to be negligibly small). In the η c (2S) decay analysis the model error is not listed separately, but the results differ noticeably from those that assume no interference. Table 5 Comparison of the results obtained under the assumption of no interference between the signal and the non-resonant contribution and those obtained with interference.

No interference
Taking interference into account  Table 2 shows that there is a large spread in the η c (2S) width values. A possible explanation of this spread is that the previous studies did not take interference into account. For each of the studied processes the interference could have a different effect on the results and shift the η c (2S) mass value significantly. Thus it is important to take interference into account.
In addition to affecting the value of the branching fraction, interference changes the Breit-Wigner shape. This effect can allow the improvement of the statistical accuracy with which the Breit-Wigner width is determined. In particular, the η c (2S) width, obtained in the present work, has a rather good accuracy, despite limited statistics and a detector resolution broader than the intrinsic width. The interference deforms the Breit-Wigner, lengthening its tail and thus improves the fit to the width (see Fig. 10).

Conclusion
We report a study of the decay B ± → K ± (cc), where the (cc) state decays to (K S Kπ) 0 and includes the η c and η c (2S) charmonia states. Both decay channels contain B ± → K ± (K S Kπ) 0 decays without intermediate charmonia that interfere with the signal. For the first time, the analysis takes interference into account with no assumptions on the phase or absolute value of the interference.
As a result, we obtain an estimate of the model error for We  Table 5).
These results are consistent with those obtained in the most accurate existing measurements. Our errors are comparable to those in other experiments despite the fact that they include additional uncertainty related to interference effects.
Thus, the function F is determined by 15 parameters listed in Table 6. Parameters α and ℑ ηS are fixed in the η c case, so they do not have a statistical error. In case the model error is not quoted, it is negligibly small compared to the statistical one. Table 6 Fit results. Parameter Function F (s, x) can be represented as a rational function of s and x: