Measurements of the branching fractions of $\eta_c\to K^+K^-\pi^0$, $K^0_S K^{\pm}\pi^{\mp}$,$2(\pi^+\pi^-\pi^0)$, and $p \bar{p}$

Using data samples collected with the BESIII detector at center-of-mass energies $\sqrt{s} = 4.23, 4.26, 4.36,$ and $4.42$~\rm{GeV}, we measure the branching fractions of $\eta_c\to K^+K^-\pi^0$, $K^0_S K^{\pm}\pi^{\mp}$, $2(\pi^+\pi^-\pi^0)$, and $p \bar{p}$, via the process $e^+e^-\to\pi^+\pi^-h_c$, $h_c\to\gamma\eta_c$. The corresponding results are $(1.15\pm0.12\pm0.10)\%$, $(2.60\pm0.21\pm0.20)\%$, $(15.2\pm1.8\pm1.7)\%$, and $(0.120\pm0.026\pm0.015)\%$, respectively. Here the first uncertainties are statistical, and the second ones systematic. Additionally, the charged track multiplicity of $\eta_c$ decays is measured for the first time.

Using data samples collected with the BESIII detector at center-of-mass energies √ s = 4. 23, 4.26, 4.36, and 4.42 GeV, we measure the branching fractions of ηc → K + K − π 0 , K 0 S K ± π ∓ , 2(π + π − π 0 ), and pp, via the process e + e − → π + π − hc, hc → γηc. The corresponding results are (1.15 ± 0.12 ± 0.10)%, (2.60 ± 0.21 ± 0.20)%, (15.2 ± 1.8 ± 1.7)%, and (0.120 ± 0.026 ± 0.015)%, respectively. Here the first uncertainties are statistical, and the second ones systematic. Additionally, the charged track multiplicity of ηc decays is measured for the first time. Many new charmonium or charmonium-like states have been discovered recently [1], which broaden our horizon on understanding the charmonium family. These states have led to a revived interest in improving the quarkmodel picture of hadrons. However, the knowledge of the lowest lying charmonium state, η c , is relatively poor compared to the other charmonium states. The reason is that most of the measurements involving η c were performed using M1 transitions from J/ψ or hindered M1 transitions from ψ(3686). In these decays, the interference between η c and non-η c amplitudes affects the η c lineshape [2]. The branching fraction (BF) of η c decays and the M1 transition rate are entangled. The insufficient understanding of the η c properties has so far prevented precise studies of η c decays themselves or of decays involving the η c . For example, in 2002, the Belle Collaboration release the measurements on the total cross section of the exclusive production of J/ψ + η c via the e + e − annihilation at the center-of-mass collision energy √ s = 10.58 GeV [3] with the result of σ[e + e − → J/ψ + η c ] × BF(η c →≥ 4 charged) = 33 +7 −6 ± 9 fb. These measurements were improved as σ[e + e − → J/ψη c (γ)]×BF(η c →≥ 2 charged) = 25.6 ± 2.8 ± 3.4 fb [4]. In 2005, the BABAR Collaboration independently measured the total cross section as 17.6 ± 2.8 +1. 5 −2.1 fb [5]. As the number of charged tracks is required in these measurements, the results will be improved if the charged tracks multiplicity is fully studied.
Recently, the E1 transition h c → γη c was found to be a perfect process to measure both η c resonant parameters and its decay BFs [6]. In addition, the h c production proceeds via ψ(3686) → π 0 h c , where the interference effect between η c and non-η c is much less than that in J/ψ, ψ(3686) radiative transition. One can draw such a conclusion according to the following calculation. The E1 transition rate, BF(h c → γη c ) = 50%, is about 2 orders of magnitude larger than that of the M 1 transition BF(ψ(3686) → γη c ) = 0.3% [7]. On the other hand, the background that can interfere with the signal comes from charmonium radiative decays, e.g. h c , ψ(3686) → γ + hadrons. If we assume the radiative decay rates of h c and ψ(3686) to be at the same level, therefore, this kind of background in the process h c → γη c should be 1 to 2 orders of magnitude less than in ψ(3686) → γη c .
BESIII has collected sizable data samples between 4.009 and 4.600 GeV (called "XYZ data" hereafter) since 2013 to study the XYZ states [9]. A large production rate of e + e − → π + π − h c has been found [10]. The total number of h c events in all these data samples combined is comparable to that from ψ(3686) → π 0 h c decays in BE-SIII data, according to the measured cross section and the corresponding integrated luminosity at each energy point. The h c is tagged by the recoil mass (RM ) of π + π − in XYZ data, while it is tagged by the recoil mass of π 0 in ψ(3686) data. Generally, the two-charged-pion mode has lower background and higher detection efficiency than the neutral pion mode.
In this paper, we report a measurement of the BFs of four η c exclusive decays via the process e + e − → π + π − h c , h c → γη c . These exclusive decays are η c → K + K − π 0 , K 0 S K ± π ∓ , 2(π + π − π 0 ), and pp, respectively. Apart from the BF measurement mentioned above, we also measure the charged tracks multiplicities in inclusive η c decays by using an unfolding method [11].

II. METHODOLOGY
The BFs of η c exclusive decays are obtained by a simultaneous fit to the RM spectrum of π + π − γ for both inclusive and exclusive modes. The BFs are common parameters independent of the center of mass energy. The numbers of the η c signal events of the exclusive and inclusive decay modes can be calculated by the following formulas, and where the subscript i denotes the different center-of-mass energy points. L and σ denote the luminosity and cross section, respectively. X denotes a certain η c exclusive decay mode, Y denotes the possible π 0 or K 0 S final state from X decay. ǫ denotes the detection efficiency determined by Monte Carlo (MC) simulations.
By comparing Eq. (1) and Eq. (2), BF(η c → X) can be extracted as In the simultaneous fit, the total number of free parameters is less than in the fits taken individually, due to common parameters such as the η c mass and width, etc. In addition, some parameters, for example, σ(e + e − → π + π − h c ), L, are not necessary in the measurement according to Eq. (3), resulting in reduced statistical uncertainties. In addition, systematic uncertainties from the same sources, e.g., the tracking efficiency of two pions from e + e − → π + π − h c , can be canceled.

III. DETECTOR AND DATA SAMPLES
The BESIII detector is a magnetic spectrometer [12] located at the Beijing Electron Positron Collider (BEPCII) [13]. The cylindrical core of the BE-SIII detector consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over 4π solid angle. The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the specific energy loss (dE/dx) resolution is 6% for the electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region. The time resolution of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps.
The data samples collected at 4 center-of-mass energies, i.e. √ s = 4.23, 4.26, 4.36, and 4.42 GeV [9], are used for our studies. Simulated samples produced with the geant4-based [14] MC package which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the backgrounds. The simulation includes the beam energy spread and initial state radiation (ISR) in the e + e − annihilations modeled with the generator kkmc [15].
The inclusive MC samples with equivalent luminosities the same as the data samples consist of the production of open charm processes, the ISR production of vector charmonium(-like) states, and the continuum processes incorporated in kkmc [15]. The known decay modes are modeled with evtgen [16] using branching fractions taken from PDG [7], and the remaining unknown decays from the charmonium states with lundcharm [17]. The final state radiations (FSR) from charged final state particles are incorporated with the photos package [18].
Signal MC samples with 200 000 events each are generated for each η c decay mode (inclusive and exclusive decays) at each center-of-mass energy. ISR is simulated using kkmc with a maximum energy for the ISR photon corresponding to the π + π − h c mass threshold. The E1 transition h c → γη c is generated with an angular distribution of 1 + cos 2 θ, where θ is the angle of the E1 photon with respect to the h c helicity direction in the h c rest frame. The inclusive decays of η c are produced similarly to the inclusive MC samples.

IV. EVENT SELECTIONS
In this analysis, the η c signal is tagged with RM (π + π − γ) by requiring RM (π + π − ) in h c signal region. For the inclusive mode, at least two charged tracks and one photon is required. For the exclusive modes, the requirements on charged tracks and photon candidates depend on their respective final state.
Charged tracks at BESIII are reconstructed from MDC hits within a polar-angle (θ) acceptance range of | cos θ| < 0.93. We require that these tracks pass within 10 cm of the interaction point in the beam direction and within 1 cm in the plane perpendicular to the beam. Tracks used in reconstructing K 0 S decays are exempted from these requirements.
A vertex fit constrains charged tracks to a common production vertex, which is updated on a run-by-run basis. For each charged track, TOF and dE/dx information is combined to compute particle identification (PID) confidence levels for the pion, kaon, and proton hypotheses.
Electromagnetic showers are reconstructed by clustering EMC crystal energies. Efficiency and energy resolution are improved by including energy deposits in nearby TOF counters. A photon candidate is defined as an isolated shower with an energy deposit of at least 25 MeV in the barrel region (| cos θ| < 0.8), or of at least 50 MeV in the end-cap region (0.86 < | cos θ| < 0.92). Showers in the transition region between the barrel and the end-cap are not well measured and are rejected. An additional requirement on the EMC hit timing suppresses electronic noise and energy deposits unrelated to the event.
A candidate π 0 is reconstructed from pairs of photons with an invariant mass in the range |M γγ − m π 0 | < 15 MeV/c 2 [7]. A one-constraint (1-C) kinematic fit is performed to improve the energy resolution, with the M γγ constrained to the known π 0 mass. We reconstruct K 0 S → π + π − candidates using pairs of oppositely charged tracks with an invariant mass in the range [7]. To reject random π + π − combinations, a secondary-vertex fitting algorithm is employed to impose the kinematic constraint between the production and decay vertices [19]. Accepted K 0 S candidates are required to have a decay length of at least twice the vertex resolution. If there is more than one π + π − combinations in an events, the one with the smallest χ 2 of the secondary vertex fit is retained.
For the selection of exclusive η c decays, the requirements on the number of photons and charged tracks are listed in Table I. A four-constraint (4C) kinematic fit imposing overall energy-momentum conservation is per-TABLE I. Requirements of the number of photons, charged tracks, π 0 , and K 0 S candidates in exclusive ηc decay modes, denoted as N charge , Nγ , N π 0 , and N K 0 S , respectively.
To determine the species of final state particles and to select the best combination when additional photons (or π 0 candidates) are found in an event, the combination with the minimum value of Vertex is selected for further analysis, where χ 2 4C is the χ 2 from the four-momentum conservation kinematic fit and χ 2 1C is the sum of the 1C (mass constraint of the two daughter photons) χ 2 of the π 0 in the final state. χ 2 PID is the χ 2 from PID of different particle hypothesis, using the energy loss in the MDC and the time measured with the TOF system, N charge is the number of the charged tracks in the final states.
is required to be not more than 50 depending on the η c decay modes, which is optimized using the figure of merit where N S is the number of signal events obtained from MC simulation (normalized to data luminosity), while N B is the number of background events obtained from the sidebands of h c in data. The require-   Table II. In addition, we require the same h c mass windows on the RM (π + π − ) spectra for both inclusive and exclusive modes.

V. NUMERICAL RESULTS OF BF(ηc → X)
A simultaneous unbinned maximum likelihood fit to the RM(π + π − γ) spectrum of the exclusive decays and the inclusive decay of η c at the four center-of-mass energies is performed to obtain the branching fractions BF(η c → X). The fit function is parameterized as follows: where the signal function is described by a Breit-Wigner function, BW (M ), convolved with the detection resolution, σ. The mass and width of BW (M ) are fixed to the η c nominal values taken from the PDG [7]. M represents the recoil mass RM (π + π − γ). The detection resolution is described by a double Gaussian function, whose parameters are obtained from MC simulations. ǫ(M ) is the efficiency curve, obtained from a fit of the efficiencies along the RM (π + π − γ) spectrum with a polynomial function and fixed in the fit to data. Figure 2 shows the efficiencies along the RM (π + π − γ) spectrum for the inclusive η c decay and the exclusive decay η c → K + K − π 0 at √ s = 4.23 GeV.
E γ = (m 2 hc − M 2 )/2m hc is the energy of the transition photon, where m hc is the h c mass [7].
is the damping factor [20], where E 0 = E γ (m ηc ) is the most probable transition energy.
B(M ) denotes the function which is used to describe the background shape. For an exclusive decay mode, a polynomial function is used. For the inclusive decay mode, it is a combination of the distribution from h c sidebands and a polynomial function. Figure 3 shows the simultaneous fit results. The fitted BFs are summarized in Table III, together with the detection efficiencies and signal yields at each energy point.

VI. CHARGED TRACK MULTIPLICITY OF ηc INCLUSIVE DECAYS
The MC simulation for the inclusive η c decay has been introduced in section III. The performance of the inclusive simulation, to some extent, can be investigated by the consistency of the charged track multiplicity [11,21,22]. Below, we introduce how to obtain the true charged track multiplicity of η c inclusive decay. An even number of charged tracks is generated in an event due to the charge conservation, while any number of charged tracks can be observed due to the detector acceptance and reconstruction efficiency. The observed charged track multiplicity of η c can be obtained by fitting for the η c signal in the π + π − γ recoil mass with the number of extra candidate tracks required to be 0, 1, 2, 3, · · · , respectively. To obtain the charged track multiplicity at the production level, an unfolding method is employed based on an efficiency matrix, whose matrix elements, ǫ ij , represent the probabilities of an event generated with j tracks being observed with i tracks. The efficiency matrix is determined from the inclusive η c MC samples. The unfolding of data is achieved by minimizing a χ 2 value, defined as where the values N obs i (i = 0, 1, 2, · · · ) are the observed multiplicities of charged tracks in the data sample, σ obs i are the corresponding uncertainties, while N j (j = 0, 2, 4, · · · ) are the true multiplicities of charged tracks at the production level in the data sample. For simplicity, the events with eight or more tracks are considered in a single value, N ≥8 , so are the efficiencies, ǫ ≥8 .   (d) show the 4 exclusive decay modes of ηc, namely, ηc → K + K − π 0 , K 0 S K ± π ∓ , 2(π + π − π 0 ), and pp. Row (e) shows the fit to the inclusive ηc decay, while (f) denotes the background-subtracted RM (π + π − γ) spectrum with the signal shape overlaid.   The systematic uncertainties on the BF measurements for exclusive η c decays from different sources are described below and listed in Table V. The total systematic uncertainty is determined by the sum in quadrature of the individual values, assuming all sources to be independent.

MDC tracking and PID
The uncertainty from the tracking efficiency and PID for the two soft pions in the process e + e − → π + π − h c cancels since the BFs are measured by a relative method, as mentioned in the introduction. We only consider the uncertainty from tracking efficiency and PID of the η c decay products. The involved charged tracks are pions (not including the pions from K 0 S decay), kaons, and protons. Their uncertainties are studied with different control samples, e + e − → π + π − K + K − for pions and kaons, e + e − → pπ −p π + (e + e − → pπ −p π + π + π − ) for protons, The uncertainties from tracking efficiency are 1% for each pion, and 2% for each kaon or proton. The uncertainties for PID are 1% for each pion, kaon or proton.

Kinematic fit
The systematic uncertainty from the kinematic fit is estimated by correcting the helix parameters of the charged tracks in the MC simulation [23]. The differences in the detection efficiency between the MC samples with and without the corrections are taken as the uncertainties due to the kinematic fit.

K 0 S reconstruction
The K 0 S reconstruction is studied with two control samples, J/ψ → K * ±K ∓ and J/ψ → φK 0 S K ± π ∓ . The difference in the K 0 S reconstruction efficiency between the MC simulation and the data is 1.2% [24], which is taken as the uncertainty due to K 0 S reconstruction.

MC model
In the MC simulation, the process e + e − → π + π − h c is modeled with a phase space (phsp) distribution. In fact, there is a confirmed intermediate state Z c (4020) and a potential intermediate state Z c (3900), in the π + π − h c final state. The uncertainty caused by the intermediate states is estimated by mixing the MC events including Z c (4020)/Z c (3900) component according to the measured fractions [10,25]. The difference in the detection efficiency is taken as the uncertainty. The uncertainty due to the inconsistency between data and MC simulation on the charged track multiplicity in inclusive η c decays is estimated based on the multiplicity obtained by the unfolding method mentioned in Sec. VI. The detection efficiency for inclusive decay can also be re-calculated with the following formula, where N i are the normalized multiplicities in data, listed in Table IV, and ǫ ij are the elements of the efficiency matrix in Eq. (5). The differences between this result and the original one are taken into account in the simultaneous fit. It is found that the influence on BF(η c → X) is negligible.

hc mass window
The uncertainty from the h c mass window is estimated by randomly changing the low and high boundaries of the h c signal region in the ranges of [3.512, 3.518] GeV/c 2 and [3.532, 3.538] GeV/c 2 and fitting the spectrum with efficiencies estimated in the corresponding intervals. The procedure is repeated for 800 times, and the distributions of the fitted BFs follow Gaussian functions. The obtained standard deviations are taken as the uncertainties due to the h c mass window selection.

Fit procedure
This uncertainty arises from the fit range, the background shape, the mass resolution, the parameters of the η c resonance, the efficiency curves, and the damping factor.
The uncertainty from the fit range is estimated by randomly changing the lower side in the range of [2.540, 2.555] GeV/c 2 and higher side in [3.200, 3.215] GeV/c 2 and repeating the fit for 800 times. The root mean square (RMS) of the resulting distributions are taken as the systematic uncertainties from the fit range.
The uncertainty due to the assumed background shape in the exclusive modes is estimated by changing the order of the Chebychev polynomial functions. For the inclusive decay mode, the h c sidebands need to be considered as well, whose systematic uncertainty is estimated by randomly changing the left and right margins of the lower and upper sidebands and repeating the fit. The procedure is performed 800 times. The left and right margins of the sidebands are changed in the ranges of [3.496, 3.450], [3.503, 3.507] GeV/c 2 and [3.543, 3.547], [3.548, 3.552] GeV/c 2 for the lower and upper sideband regions, respectively. The distributions of the fitted results follow Gaussian functions, and the standard deviations are taken as the uncertainties from the h c sidebands selection. The uncertainty from the polynomial is estimated by changing the order of the polynomial.
The discrepancy between data and MC simulation on detection resolution is estimated by a control sample, ψ(2S) → π + π − J/ψ, J/ψ → γη ′ , η ′ → γπ + π − . By fitting the η ′ signals, we can obtain the mass resolution for both data and MC. We change the mass resolutions according to the result obtained from control sample to re-fit the RM (γπ + π − ). The differences on the BFs with and without changing the mass resolution are taken as the systematic uncertainties.
The η c resonance parameters are fixed to the world average values in the fit. We change these values by ±1σ, and the larger difference is taken as the uncertainty.
The efficiency curves, as shown in Fig. 2, change slowly with RM (π + π − γ). We find only a very small change in results when constant efficiencies are used. Therefore, the uncertainties due to efficiencies can be neglected.
The uncertainty from the damping factor is estimated by using an alternative form of the damping factor, which is used in the CLEO's published paper [26]. The differences between the results with the two forms of damping factor are taken as the systematic uncertainty.

B. Charged track multiplicity
The systematic uncertainties on the charged track multiplicity in η c inclusive decay from different sources are described below and listed in Table VI. They are estimated in a similar way as introduced in Sec. VII A. The total systematic uncertainty is determined by the sum in quadrature of the individual values, assuming that all the sources are independent.

MDC tracking and PID
The uncertainties from MDC tracking and PID are the same as those in the measurement of BF(η c → X).

hc mass window
The uncertainties are estimated by changing the h c mass window from [3.515, 3.535] GeV/c 2 to [3.518, 3.532] GeV/c 2 and [3.512, 3.538] GeV/c 2 . The largest changes on the multiplicity are taken as the uncertainty.

MC model
Similar to that in measurement of BF(η c → X), the uncertainty due to MC model mainly comes from the potential intermediate states and the inclusive η c decay. The uncertainty from the former is estimated as before, while the latter is estimated by removing the unknown modes simulated by lundcharm model, and only considering the known η c decay modes.

Fit
The uncertainties due to the fit to the recoil mass spectra of π + π − γ are evaluated by varying the fit range, sideband ranges, mass resolution, resonant parameters of η c , and damping factors used in the fit, in similar ways as introduced in Sec. VII A. The spreads of the results obtained with the alternative assumptions are used to assign the systematic uncertainties.

VIII. SUMMARY
In summary, with the data samples collected at √ s =4.23, 4.26, 4.36, and 4.42 GeV, by comparing the exclusive and inclusive decays of η c , we determine the BFs for η c → K + K − π 0 , K 0 S K ± π ∓ , 2(π + π − π 0 ), and pp via e + e − → π + π − h c , h c → γη c . The results are presented in Table VII; they agree with previous measurements by BESIII [6] within uncertainties, while the accuracy of these BFs is improved. With this improved accuracy, the measurements of the M1 transitions of J/ψ → γη c and ψ(3686) → γη c can be more precise, since such measurements provide combined results of BF(J/ψ(ψ(3686)) → γη c ) × BF(η c → X).
Moreover, the charged track multiplicity of η c inclusive decay at production level is quantitatively presented for the first time in Table IV. The good consistency between data and MC simulation for this charged track multiplicity indicates that the current MC simulation works generally well. With this charged track multiplicity, many studies with η c in the final state [27] are possible with higher precision than previously.

IX. ACKNOWLEDGEMENT
The BESIII collaboration thanks the staff of BEPCII and the IHEP computing center for their strong support. This work is supported in part by National