Study of charmonium production in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${b} $$\end{document}b-hadron decays and first evidence for the decay \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{B}} ^0_{{s}}} \!\rightarrow \phi \phi \phi $$\end{document}Bs0→ϕϕϕ

Using decays to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document}ϕ-meson pairs, the inclusive production of charmonium states in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${b} $$\end{document}b-hadron decays is studied with pp collision data corresponding to an integrated luminosity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3.0 {\,\mathrm{fb}}^{-1} $$\end{document}3.0fb-1, collected by the LHCb experiment at centre-of-mass energies of 7 and 8 TeV. Denoting by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {B}} _C \equiv {\mathcal {B}} ( {{b}} \!\rightarrow C X ) \times {\mathcal {B}} ( C\!\rightarrow \phi \phi )$$\end{document}BC≡B(b→CX)×B(C→ϕϕ) the inclusive branching fraction of a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{b}} $$\end{document}b hadron to a charmonium state C that decays into a pair of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document}ϕ mesons, ratios \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{C_1}_{C_2}\equiv {\mathcal {B}} _{C_1} / {\mathcal {B}} _{C_2}$$\end{document}RC2C1≡BC1/BC2 are determined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{{\upchi _{{{c}} 0}}}_{{\eta _{{c}}} (1S)} = 0.147 \pm 0.023 \pm 0.011$$\end{document}Rηc(1S)χc0=0.147±0.023±0.011, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{{\upchi _{{{c}} 1}}}_{{\eta _{{c}}} (1S)} = 0.073 \pm 0.016 \pm 0.006$$\end{document}Rηc(1S)χc1=0.073±0.016±0.006, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{{\upchi _{{{c}} 2}}}_{{\eta _{{c}}} (1S)} = 0.081 \pm 0.013 \pm 0.005$$\end{document}Rηc(1S)χc2=0.081±0.013±0.005, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{{\upchi _{{{c}} 1}}}_{{\upchi _{{{c}} 0}}} = 0.50 \pm 0.11 \pm 0.01$$\end{document}Rχc0χc1=0.50±0.11±0.01, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{{\upchi _{{{c}} 2}}}_{{\upchi _{{{c}} 0}}} = 0.56 \pm 0.10 \pm 0.01$$\end{document}Rχc0χc2=0.56±0.10±0.01 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{{\eta _{{c}}} (2S)}_{{\eta _{{c}}} (1S)} = 0.040 \pm 0.011 \pm 0.004$$\end{document}Rηc(1S)ηc(2S)=0.040±0.011±0.004. Here and below the first uncertainties are statistical and the second systematic. Upper limits at 90% confidence level for the inclusive production of X(3872), X(3915) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\upchi _{{{c}} 2}} (2P)$$\end{document}χc2(2P) states are obtained as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{X(3872)}_{{\upchi _{{{c}} 1}}} < 0.34$$\end{document}Rχc1X(3872)<0.34, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{X(3915)}_{{\upchi _{{{c}} 0}}} < 0.12$$\end{document}Rχc0X(3915)<0.12 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^{{\upchi _{{{c}} 2}} (2P)}_{{\upchi _{{{c}} 2}}} < 0.16$$\end{document}Rχc2χc2(2P)<0.16. Differential cross-sections as a function of transverse momentum are measured for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\eta _{{c}}} (1S)$$\end{document}ηc(1S) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi _c$$\end{document}χc states. The branching fraction of the decay \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{B}} ^0_{{s}}} \!\rightarrow \phi \phi \phi $$\end{document}Bs0→ϕϕϕ is measured for the first time, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {B}} ( {{{B}} ^0_{{s}}} \!\rightarrow \phi \phi \phi ) = (2.15 \pm 0.54 \pm 0.28 \pm 0.21_{{\mathcal {B}}}) \times 10^{-6}$$\end{document}B(Bs0→ϕϕϕ)=(2.15±0.54±0.28±0.21B)×10-6. Here the third uncertainty is due to the branching fraction of the decay \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{B}} ^0_{{s}}} \!\rightarrow \phi \phi $$\end{document}Bs0→ϕϕ, which is used for normalization. No evidence for intermediate resonances is seen. A preferentially transverse \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document}ϕ polarization is observed. The measurements allow the determination of the ratio of the branching fractions for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\eta _{{c}}} (1S)$$\end{document}ηc(1S) decays to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi \phi $$\end{document}ϕϕ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{p}} {\overline{{{p}}}} $$\end{document}pp¯ as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {B}} ( {\eta _{{c}}} (1S)\!\rightarrow \phi \phi )/{\mathcal {B}} ( {\eta _{{c}}} (1S)\!\rightarrow {{p}} {\overline{{{p}}}} ) = 1.79 \pm 0.14\pm 0.32$$\end{document}B(ηc(1S)→ϕϕ)/B(ηc(1S)→pp¯)=1.79±0.14±0.32.

Here the third uncertainty is due to the branching fraction of the decay B 0 s → φφ, which is used for normalization. No evidence for intermediate resonances is seen. A preferentially transverse φ polarization is observed. The measurements allow the determination of the ratio of the branching fractions for the η c (1S) decays to φφ and p p as B(η c (1S) → φφ)/B(η c (1S) → p p) = 1.79 ± 0.14 ± 0.32.

Introduction
The production of the J PC = 1 −− charmonium states has been extensively studied using decays to clean dilepton final states. Other states such as those from the χ c family can be accessed via the radiative transition to J/ψ. Studies of the production of the non-1 −− charmonium states can be performed by reconstructing their decays to fully hadronic final e-mail: sergey.barsuk@cern.ch states [1]. This paper reports a measurement of the inclusive production rates of the η c and χ c states in b-hadron decays, b → η c X and b → χ c X , using charmonia decays to a pair of φ mesons. In addition, the first evidence for the decay B 0 s → φφφ is reported.
Results on inclusive charmonium production in b-hadron decays are available from e + e − experiments operating at centre-of-mass energies around the ϒ(4S) and ϒ(5S) resonances, studying mixtures of B + and B 0 mesons 1 (light mixture) or B + , B 0 and B 0 s mesons, respectively. Mixtures of all b-hadrons (B + , B 0 , B 0 s , B + c and b-baryons) have been studied at LEP, the Tevatron and the LHC. The world average values for charmonium branching fractions from the light mixture are dominated by results from the CLEO [2,3], Belle [4] and BaBar [5] collaborations. For the J/ψ, ψ(2S) and χ c1 states the measured branching fractions are consistent within uncertainties. The new Belle result for the b → χ c2 X branching fraction [4], which supersedes the previous measurement [6], is below the BaBar result [5] by more than 2.5 standard deviations, while the CLEO collaboration does not observe a statistically significant b → χ c2 X signal [3]. An upper limit on the inclusive production rate of η c (1S) mesons in the light mixture, B(B → η c (1S)X ) < 9 × 10 −3 at 90% confidence level (CL), was reported by CLEO [7].
The branching fractions of b-hadron decays to final states including a J/ψ or ψ(2S) charmonium state, where all b-hadron species are involved, are known with uncertainties of around 10%, with the world averages dominated by the measurements performed at LEP [8][9][10]. The ratio of b → ψ(2S)X and b → J/ψX yields has been measured at the LHC by the LHCb, CMS and ATLAS collaborations with a precision of around 5% [11][12][13]. The only available results for the χ c family are the χ c1 inclusive production rates in b-hadron decays measured by the DELPHI and L3 collaborations [8,9], with an average value of B(b → χ c1 X ) = (14 ± 4) × 10 −3 [14]. Recently, LHCb measured the η c (1S) production rate, B(b → η c (1S)X ) = (4.88 ± 0.64 ± 0.29 ± 0.67 B ) × 10 −3 , where the third uncertainty is due to uncertainties on the J/ψ inclusive branching fraction from b-hadron decays and the branching fractions of the decays of J/ψ and η c (1S) to the p p final state [15].
While experimentally the reconstruction of charmonia from b-hadron decays allows an efficient control of combinatorial background with respect to charmonium candidates produced in the p p collision vertex via hadroproduction or in the decays of heavier resonances (prompt charmonium), inclusive b-hadron decays to charmonia are theoretically less clean. Since a description of the strong interaction dynamics in b-hadron inclusive decays improves with respect to exclusive decays due to consideration of more final states, and the formation of charmonium proceeds through a shortdistance process, a factorization of a cc pair production and its hadronization in a given charmonium state becomes a reasonable assumption [16]. The relative inclusive production of χ c states in b-hadron decays provides a clean test of charmonia production models. For example, the colour evaporation model predicts a χ c2 /χ c1 production ratio of 5/3 [17], while the perturbative QCD-based computation predicts that the V-A current, which is responsible for the b decays, forbids the χ c2 and χ c0 production at leading order. In the non-relativistic QCD (NRQCD) framework [18][19][20], the colour-octet contributions have to be included, predicting the rates to be proportional to (2J +1) for the χ c J states. The NRQCD framework can be applied to both prompt charmonium production and secondary production from b-hadron decays and the comparison between these two production mechanisms can provide a valuable test of this theoretical framework.
In this paper we report the first measurements of inclusive χ c and η c (2S) production rates in b-hadron decays using charmonium decays to hadronic final states in the high-multiplicity environment of a hadron collider. Experimentally, charmonium candidates from b-hadron decays are distinguished from prompt charmonia by exploiting the bhadron decay time and reconstructing a b-hadron (and charmonium) decay vertex well separated from the primary vertex where the b-hadron candidate was produced. The charmonium states are reconstructed via their decays to a φφ final state. The η c (1S) production followed by the decay η c (1S) → φφ is used for normalization, so that systematic uncertainties partially cancel in the ratios. As a by-product of the production rate measurements, the masses of the η c (1S), χ c0 , χ c1 , χ c2 and η c (2S) charmonium states and the natural width of the η c (1S) meson are determined.
The B 0 s decay to the φφ final state has been observed by the CDF collaboration [21] and recently precisely measured by the LHCb collaboration [22], where it was also used to search for CP-violating asymmetries [23]. In the Standard Model (SM) the amplitude for the decay B 0 s → φφ is dominated by a loop diagram. Experimental verification of the partial width, polarization amplitudes and triple-product asymmetries of the B 0 s → φφ decay probes the QCD contribution to the weak processes described by nonfactorizable penguin diagrams [24,25], and contributions from particles beyond the SM to the penguin loops [26][27][28][29][30]. A three-body B 0 s → φφφ decay leads to a final state with six strange quarks. In the SM it is described by the penguin diagram of the B 0 s → φφ decay with the creation of an additional ss quark pair. The B 0 s → φφφ decay can also receive contributions from an intermediate charmonium state decaying to a φφ state. Here we report first evidence for the B 0 s → φφφ decay and study its resonance structure. The branching fraction of this decay is determined relative to the branching fraction B(B 0 s → φφ) [22]. To cross-check the technique exploited in this paper, the value of B(B 0 s → φφ) is also determined relative to the η c (1S) production rate. Finally, the ratio of the branching fractions for the decays η c (1S) → φφ and η c (1S) → p p is determined, using additional external information.
The LHCb detector and data sample used for the analysis are presented in Sect. 2. Section 3 explains the selection details and the signal extraction technique. Inclusive production of charmonium states in b-hadron decays is discussed in Sect. 4. In Sect. 5 measurements of the η c (1S) mass and natural width are described. First evidence for the B 0 s → φφφ decay is reported in Sect. 6. The main results of the paper are summarized in Sect. 7.

LHCb detector and data sample
The LHCb detector [31,32] 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 includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum 2 to 1.0% at 200 GeV. The minimum distance of a track to a primary vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/ p T ) μm, where p T is the component of the momentum transverse to the beam, in GeV. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillatingpad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers. The online event selection is performed by a trigger, 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.
The analysis is based on p p collision data recorded by the LHCb experiment at a centre-of-mass energy √ s = 7 TeV, corresponding to an integrated luminosity of 1.0 fb −1 , and at √ s = 8 TeV, corresponding to an integrated luminosity of 2.0 fb −1 . Events enriched in signal decays are selected by the hardware trigger, based on the presence of a single deposit of high transverse energy in the calorimeter. The subsequent software trigger selects events with displaced vertices formed by charged particles having a good track-fit quality, transverse momentum larger than 0.5 GeV, and that are incompatible with originating from any PV [23]. Charged kaon candidates are identified using the information from the Cherenkov and tracking detectors. Two oppositely charged kaon candidates having an invariant mass within ±11 MeV of the known mass of the φ meson are required to form a good quality vertex.
Precise mass measurements require a momentum-scale calibration. The procedure [33] uses J/ψ → μ + μ − decays to cross-calibrate a relative momentum scale between different data-taking periods. The absolute scale is determined using B + → J/ψK + decays with known particle masses as input [14]. The final calibration is checked with a variety of fully reconstructed quarkonium, B + and K 0 S decays. No residual bias is observed within the experimental resolution.
In the simulation, pp collisions are generated using Pythia [34,35] with a specific LHCb configuration [36]. Decays of hadronic particles are described by EvtGen [37], in which final-state radiation is generated using Photos [38]. The interaction of the generated particles with the detector material and its response are implemented using the Geant4 toolkit [39,40] as described in Ref. [41]. Simulated samples of η c and χ c mesons decaying to the φφ final state, and B 0 s decaying to two and three φ mesons, are used to estimate efficiency ratios and to evaluate systematic uncertainties. Charmonium and B 0 s → φφφ decays are generated with uniform phase-space density, while B 0 s → φφ decays are generated according to the amplitudes from Ref. [21].

Selection criteria and signal extraction
The signal selection is largely performed at the trigger level. The offline analysis selects combinations of two or three φ candidates that are required to form a good-quality common vertex displaced from the primary vertex. A good separation between the two vertices (χ 2 > 100) is required, reducing the contribution from charmonia produced directly at the primary vertex to below 1%. Pairs of φ mesons origi-nating from different b-hadrons produced in the same beam crossing event are suppressed by the requirement of a goodquality common vertex. Detector acceptance and selection requirements, and in particular the requirement of the kaon p T to exceed 0.5 GeV, significantly suppress φφ combinations with p T below 4 GeV.
Two-dimensional (2D) or three-dimensional (3D) extended unbinned maximum likelihood fits, corresponding to the two or three randomly ordered K + K − combinations, are performed in bins of the invariant mass of the four-kaon and six-kaon combinations, denoted as 2(K + K − ) or 3(K + K − ), respectively, to determine the numbers of φφ and φφφ combinations. The 2D fit accounts for the signal, φφ, and background, φ (K + K − ) and 2(K + K − ), components, while the 3D fit accounts for the signal, φφφ, and background, φφ (K + K − ), φ 2(K + K − ) and 3(K + K − ), components. A φ signal is described by the convolution of a Breit-Wigner function and a sum of two Gaussian functions with a common mean. The ratio of the two Gaussian widths and the fraction of the narrower Gaussian are taken from simulation. The contribution from combinatorial background, due to K + K − pairs not originating from the decay of a φ meson, is assumed to be flat. In addition, a threshold function to account for the available phase-space in the K + K − system is introduced for both signal and combinatorial background. While no visible contribution from the f 0 (980) resonance decaying into a K + K − pair is observed in the 2(K + K − ) or 3(K + K − ) combinations, a potential effect due to contributions from such decays is considered as a source of systematic uncertainty. Figures 1 and 2 show the results of the 2D fits to the 2(K + K − ) invariant mass distributions along with the projections to the K + K − invariant mass axes in the η c (1S) and B 0 s signal regions, 2.91 − 3.06 GeV and 5.30 − 5.43 GeV. Figure 3 shows the projections to the K + K − invariant mass axes of the 3D fit to the 3(K + K − ) invariant mass distribution in the B 0 s signal region. While using the known value for the natural width of the φ resonance [14], the φ mass and the remaining resolution parameter are determined from the fits in the enlarged signal φφ and φφφ invariant mass regions. In the 2D and 3D fits in the bins of φφ and φφφ invariant mass the φ mass and the resolution parameter are then fixed to the values determined in the enlarged signal regions.
Unless they are extracted from the 2D or 3D fits, throughout the paper the error bars shown in the histograms are determined as follows: the upper (lower) error bar assigned to a given bin content N is defined by the expectation value λ of the Poissonian distribution giving 16% probability to observe the number of events N or less (more). When obtained from the 2D or 3D fits the histogram bin contents are constrained to be positive, with error bars defined by the range in the allowed region where the best fit negative-log likelihood value is within half a unit from the global minimum. Candidates/(1 MeV) Candidates/(1 MeV)  In the following, production ratios are determined from the signal yields obtained from the fits of the φφ or φφφ invariant mass spectra. The relative efficiencies are taken into account to determine the ratio of the branching fractions of the corresponding decays. Signal yields corresponding to the data samples accumulated at √ s = 7 and 8 TeV are found to be compatible. Unless otherwise specified, the results below are based on the analysis of the combined data sample.
from the η c (1S) and η c (2S) mesons, and the χ c0 , χ c1 and χ c2 mesons. The charmonium-like states X (3872), X (3915) and χ c2 (2P) with masses and natural widths from Ref. [14] are taken into account in alternative fits in order to evaluate systematic uncertainties, as well as to obtain upper limits on the inclusive production of these states in b-hadron decays. Each resonance is described by the convolution of a relativistic Breit-Wigner function and a sum of two Gaussian functions with a common mean. The combinatorial background, i.e. contributions due to random combinations of two true φ mesons, is described by the product of a first-order polynomial with an exponential function and a threshold factor. The natural width of the η c (1S) state is a free parameter in the fit, while the natural widths of the η c (2S) and the χ c states, which have lower signal yields, are fixed to their known values [14]. Possible variations of the η c (2S) production rate depending on its natural width, η c (2S) , are explored by providing the result as a function of the assumed natural width. The ratio of the two Gaussian widths and the fraction of the narrow Gaussian are fixed from the simulation. The mass resolution for different charmonium resonances is scaled according to the energy release, so that a single free parameter in the φφ invariant mass fit accounts for the detector resolution. This scaling of the mass resolution for different charmonium states has been validated using simulation. The signal yields are given in Table 1. The ratios of the resonance yields from the fit are given in Table 2, both for the Table 2 The ratio of charmonium signal yields with respect to the η c (1S) yield and between pairs of χ c states. The first uncertainties are statistical and the second systematic

Resonances
Signal yield ratio ratios with respect to the η c (1S) yield and between pairs of χ c states; the systematic uncertainties are discussed below. The statistical significance for the N η c (2S) signal is estimated from the χ 2 -profile to be 3.7 standard deviations. Systematic uncertainties in the ratios of the charmonium yields are estimated by considering potential contributions from other states, from imperfect modelling of detector resolution, signal resonances and background, and from the 2D fit technique. In order to evaluate the systematic uncertainty related to potential contributions from other states, signal shapes for the X (3872), X (3915), and χ c2 (2P) states are included in the fit. Systematic uncertainties related to detector resolution are estimated by fixing the η c (1S) mass resolution to the value determined from the simulation. In addition, systematic uncertainties associated to the impact of the detector resolution on the signal shapes are estimated by comparing the nominal fit results to those obtained using a single instead of a double Gaussian function. An uncertainty associated with the description of the mass resolution of the φ meson is estimated by fixing the resolution in the 2D fits to the value determined from simulation. The uncertainty associated with the description using the relativistic Breit-Wigner line shape [42] is estimated by varying the radial parameter of the Blatt-Weisskopf barrier factor [43] between 0.5 and 3 GeV −1 . In order to estimate the uncertainty related to the natural width of the η c (2S) meson, the value of η c (2S) is varied within the uncertainties of the world average [14]. The uncertainty in the description of the χ c signal peaks is estimated by fixing the χ c masses to their known values. A reduced fit range, covering only the χ c and η c (2S) region (3.15−3.95 GeV), is used to estimate a systematic uncertainty associated to the choice of the fit range. An alternative background parametrization using a parabolic instead of a linear function is used to estimate the systematic uncertainty due to the choice of the background parametrization. A systematic uncertainty associated to the background parametrization in the 2D fits is estimated by adding slope parameters to the description of the non-φ K + K − combinations in the φK + K − and the 2 × (K + K − ) components. The effect of a potential contribution from the f 0 (980) state in the 2D fits is estimated by including the f 0 (980) contribution following the analysis described in Ref. [44], and varying the f 0 (980) mass and natural width within the uncertainties quoted in Ref. [14]. Contributions from multiple candidates are below 2% per event and are not considered in the evaluation of systematic uncertainties. The uncertainty related to the momentum-scale calibration is negligible. The total systematic uncertainty is obtained as the quadratic sum of the individual systematic contributions. The systematic uncertainties are shown in Table 3. The description of the background and the potential contributions from other resonances dominate the total systematic uncertainties. In the yield ratios the systematic uncertainty is smaller than or comparable to the statistical uncertainty.
Complementary cross-checks are performed by comparing the distributions of kinematic variables in simulation and data. The stability of the results is checked by using an alternative binning of the φφ invariant mass distribution (shifted by half a bin) and by using the weighted signal events from the sPlot [45] instead of the nominal 2D fit technique. The same check is performed for the determination of the masses and widths of the charmonium states. In all cases no significant changes are observed and no additional contributions to the systematic uncertainties are assigned.

Production of η c and χ c in b-hadron decays
The production ratios of charmonium C with respect to the η c (1S) yield and between pairs of χ c states where, here and throughout, the first uncertainties are statistical and the second are systematic. The dominant contributions to the systematic uncertainty on the relative χ c production rates arise from accounting for possible other resonances and from uncertainties on the χ c masses [14]. The systematic uncertainties are smaller than the statistical uncertainties, so that the overall uncertainty on the measurements will be reduced with a larger dataset. The systematic uncertainty on the relative production rate of η c (2S) mesons is dominated by possible contributions from other resonances and variation of the η c (2S) natural width.
In the following, the η c (1S) production rate in b-hadron decays and the branching fractions of the charmonium decays to φφ are used to obtain single ratios of charmonium production rates in b-hadron decays and the branching fractions of inclusive b-hadron transitions to charmonium states. The η c (1S) inclusive production rate in b-hadron decays was measured by LHCb as B(b → η c (1S)X ) = (4.88 ± 0.97) × 10 −3 [15] using decays to p p. Branching fractions of the charmonia decays to φφ from Ref. [14] are used. However, Ref. [14] indicates a possible discrepancy for the B(η c (1S) → φφ) value when comparing a direct determination and a fit including all available measurements. Therefore, an average of the results from Belle [46] and BaBar [47] using B + decays to φφK + , B(η c (1S) → φφ) = (3.21 ± 0.72) × 10 −3 , is used below. The uncertainty of this average dominates the majority of the following results in this section, and an improved knowledge of the B(η c (1S) → φφ) is critical to reduce the uncertainties of the derived results. The values B(χ c0 → φφ) = (7.7 ± 0.7) × 10 −4 , B(χ c1 → φφ) = (4.2 ± 0.5) × 10 −4 , and B(χ c2 → φφ) = (1.12 ± 0.10) × 10 −3 , are used for the χ c decays [14]. The relative branching fractions of b-hadron inclusive decays into χ c states are then derived to be where the third uncertainty is due to those on the branching fractions B(χ c → φφ). The result for the relative χ c1 and χ c2 production in inclusive b-hadron decays is close to the values measured in B 0 and B + production [14]. The branching fractions of b-hadron decays into χ c states relative to the η c (1S) meson are where the dominating uncertainty is due to the uncertainty of the branching fractions B(η c (1S) → φφ) and B(χ c → φφ).
The absolute branching fractions are determined to be where the third uncertainty is due to the uncertainties on the branching fractions of the b-hadron decays to the η c (1S) meson, B(b → η c (1S)X ), and of η c (1S) and χ c decays to φφ. The branching fraction of b-hadron decays into χ c0 is measured for the first time, and is found to be larger than the values predicted in Ref. [18]. Throughout the paper comparisons of the results on the production of charmonium states to theory predictions neglect the fact that the measured branching fractions also contain decays via intermediate higher-mass charmonium resonances, whereas theory calculations consider only direct b-hadron transitions to the charmonium state considered. Among the contributions that can be quantified the most sizeable comes from the ψ(2S) state decaying to the χ c states. With the branching fraction B(b → ψ(2S)X ) recently measured [11] by LHCb the branching fractions B(b → χ c X ), measured in this paper, are influenced by about 10%. The branching fractions B(b → χ c0 X ) and B(b → χ c2 X ) remain different from the predictions in Ref. [18].
The branching fraction measurement for b-hadron decays into χ c1 is most precise in mixtures of B 0 , B + , B 0 s , B + c and b baryons. The central value is lower than the central values measured by the DELPHI [8] and L3 [9] experiments at LEP, (11.3 +5.8 −5.0 ± 0.4) × 10 −3 and (19 ± 7 ± 1) × 10 −3 , respectively. For the measurements with different b-hadron content, the LHCb result is consistent with measurements by CLEO [2], Belle [4], and BaBar [5]. Finally, the LHCb result for the inclusive b-hadron decays into χ c1 is consistent with the prediction in Ref. [18].
The branching fraction of b-hadron decays into χ c2 is measured for the first time with a mixture of B 0 , B + , B 0 s , B + c and b-baryons. The result is consistent with the world average [14] measured with the B 0 and B + mixture, and with individual results from CLEO [3], Belle [4] and BaBar [5]. The value obtained is below the range predicted in Ref. [18]. A deviation of the η c (2S) natural width from the world average value [14] would affect the measured ratio of η c (2S) and η c (1S) production rates in b-hadron inclusive decays, as shown in Fig. 5. The decay η c (2S) → φφ has not been observed so far. Hence the product of the branching fraction of b-hadron decays to η c (2S) and the branching fraction of the η c (2S) → φφ decay mode is determined as where the systematic uncertainty is dominated by the uncertainty on the η c (1S) production rate in b-hadron decays. This is the first evidence for η c (2S) production in b-hadron decays, and for the decay of the η c (2S) meson into a pair of φ mesons.

Transverse momentum dependence of the differential cross-sections for η c (1S) and χ c production
The shapes of the differential production cross-sections as a function of transverse momentum are studied in the LHCb acceptance (2 < η < 5) and for 3 < p T < 17 GeV and 2 < p T < 19 GeV for the η c (1S) and χ c states, respectively. Each differential production cross-section is normalized to the production cross-section integrated over the studied p T region. Figure 6 shows the normalized differential crosssections of η c (1S), χ c0 , χ c1 and χ c2 production at √ s = 7 and 8 TeV. An exponential function proportional to exp(−α p T ) is fitted to the integral of the each bin of the distributions. No significant difference is observed between the √ s = 7 TeV and 8 TeV data. The results for the slope parameters α are given in Table 4. For χ c1 and χ c2 production in b-hadron     [48] in p T and rapidity.

Search for production of X (3872), X (3915) and χ c2 (2P)
The observation of the X (3915) and χ c2 (2P) states in bhadron decays or the X (3872) decaying to a pair of φ mesons would provide interesting information on the properties of these states. The invariant mass spectrum of φφ combinations in Fig. 4 shows no evidence for a signal from the X (3872), X (3915), or χ c2 (2P) states. Bayesian upper limits assuming a uniform prior in the event yields are obtained on the branching fractions relative to those involving decays to the states with similar quantum numbers. For the states with similar quantum numbers, in the efficiency ratios systematic uncertainties largely cancel. Using the efficiency ratios from the simulation, the upper limits at 95% (90%) CL on the ratios of inclusive branching fractions are < 0.14 (0.12), Using the measured production rates of the χ c states in bhadron decays and branching fractions for the χ c decays to the φφ final state [14], the upper limits at 95% (90%) CL on the production rates of the X (3872), X (3915), and χ c2 (2P) states in b-hadron decays are

Masses and natural widths of charmonium states
The majority of the η c (1S) mass measurements, used in the fit of Ref. [14], were performed with two-photon pro-  [49][50][51], do not provide consistent results. In 2009, the CLEO collaboration observed a significant asymmetry in the line shapes of radiative J/ψ → γ η c (1S) and ψ(2S) → γ η c (1S) transitions [52], which, when ignored, could lead to significant bias in the mass and width measurement via J/ψ or ψ(2S) radiative decays. Recent BES III results [53,54] obtained using radiative decays of ψ(2S), shifted the world average value by more than two standard deviations. Therefore precise η c (1S) mass measurements using a different technique are needed. LHCb measured M η c (1S) = 2982.2 ± 1.5 ± 0.1 MeV [15] using η c (1S) from b-hadron decays and reconstructing η c (1S) via decays to p p. A similar situation occurs with the η c (1S) natural width determination, where recent BES III results obtained using radiative decays of ψ(2S) shifted the world average from 29.7 ± 1.0 MeV to 31.8 ± 0.8 MeV.
The properties of the η c (2S) state are less well studied. Measurements at the CLEO [55], BaBar [56,57], Belle [51,58] and BES III [59,60] experiments, using γ γ → η c (2S) → hadrons, double charmonium production in e + e − annihilation, exclusive B decays and radiative transitions of ψ(2S), yield the world averages [14] of 3639.4 ± 1.3 MeV for the η c (2S) mass, and 11.3 +3.2 −2.9 MeV for its natural width. Table 5 presents measurements of the masses of the η c and χ c states and of the natural width of the η c (1S) from the fit of the φφ invariant mass spectrum in Fig. 4. For the determination of the systematic uncertainties, except for the test of the impact of the f 0 (980) meson, the same variations of the analysis are performed as for the determination of the charmonium yields. In addition, the effect of excluding the η c (2S) mass region (2.8−3.7 GeV) is studied, and the uncertainties related to the momentum-scale calibration are estimated by varying the calibration parameter by ±3×10 −4 [33]. The resulting total systematic uncertainty is obtained as the quadratic sum of the individual contributions. The uncertainty related to the momentum-scale calibration dominates the mass determination for all η c and χ c states. The uncertainty of the η c (1S) measurement is dominated by the background description.
The measured charmonium masses agree with the world averages [14]. The measured η c (1S) mass is in agreement with the previous LHCb measurement using decays to the p p final states [15] and has a better precision. The precision obtained for the η c (1S) mass is comparable to the precision of the world average value. The value of the η c (1S) natural width is consistent with the world average [14].
The charmonium mass differences M χ c1 − M χ c0 , M χ c2 − M χ c0 , and M η c (2S) − M η c (1S) are obtained (Table 6) as a    Figure 7 shows the η c (1S) , M η c (1S) contour plot, obtained from the analysis of b-hadron decays into η c mesons, where the η c candidates are reconstructed via the decay η c (1S) → φφ. The measurements of the η c (1S) mass and natural width using η c (1S) meson decays to φφ are consistent with the studies using decays to p p [15] and with the world average [14]. The measured η c (1S) mass is below the result in Ref. [61]. The precision obtained on the η c (1S) mass is comparable to the precision of the world average.

First evidence of the B 0 s → φφφ decay
In order to extract φφφ combinations a 3D extended unbinned maximum likelihood fit is used, as described in Sect. 3. Figure 8 shows the invariant mass distribution for φφφ combinations. The fit to the invariant φφφ mass spectrum is performed using a sum of two Gaussian functions with a common mean to describe the B 0 s signal, and an exponential function to describe combinatorial background. The ratio of the two Gaussian widths and the fraction of the nar- This measurement World average [14] Inclusive pp [15] Exclusive [60] Γ Only statistical uncertainties are shown. The red cross, black square and blue triangle with error bars indicate the world average [14], the result from Ref. [15], and the result from Ref. [61], respectively row Gaussian are taken from simulation so that a single free parameter in the φφφ invariant mass fit accounts for the detector resolution. A signal of 41 ± 10 ± 5 B 0 s decays over a low background of about 3 events is obtained. Uncertainties related to the background description in the 3D fit and to the decay model defining the φ polarization dominate the systematic uncertainty in the B 0 s signal yield determination. The significance of the signal is estimated from the distributions of the difference in the logarithm of the best-fit χ 2 with and without including the signal shape in toy simulation samples. This leads to a signal significance of 4.9 standard deviations.
The B 0 s → φφ decay mode is chosen as a normalization mode for the B(B 0 s → φφφ) measurement. The invariant mass spectrum obtained from 2D fits in bins of the φφ invariant mass in the region of the B 0 s mass is shown in Fig. 9. A sum of two Gaussian functions with a common mean is used to describe the B 0 s signal shape, while an exponential function models the combinatorial background. The ratio of the two Gaussian widths and the fraction of the narrow Gaussian function are taken from simulation. In total 2701 ± 114 ± 84 B 0 s decays are found. The uncertainties related to the description of the resolution in the 2D fits and the description of the φφ invariant mass resolution dominate the systematic uncertainty in the B 0 s signal yield determination. The ratio of the B 0 s → φφφ and B 0 s → φφ branching fractions is obtained from the relative B 0 s → φφφ and B 0 s → φφ signal yields and their efficiencies as In the above expression, the event yields are determined from the fits. The efficiency ratio, ε B 0 s →φφφ /ε B 0 s →φφ = 0.26 ± 0.01, is obtained from simulation and corrected to account for different B 0 s transverse momentum spectra in data and simulation. The B 0 s → φφφ transition is assumed to proceed via a three-body decay with uniform phase-space density. This assumption is supported below by comparing the φφ invariant mass distribution in data and simulation. The systematic uncertainty is dominated by the uncertainty in polarization of the φ mesons in the decay B 0 s → φφφ, as discussed at the end of this section. Using the branching fraction of the B 0 s → φφ decay, B(B 0 s → φφ) = (1.84 ± 0.05 ± 0.07 ± 0.11 f s / f d ± 0.12 norm ) × 10 −5 [22], the branching fraction for the B 0 s meson decay to three φ mesons is determined to be where the last uncertainty is due to the branching fraction B(B 0 s → φφ). The A phase-space distribution as obtained from simulation is overlaid for comparison. No indication of significant contributions from η c , χ c , f 2 (2300) or f 2 (2340) states is seen. A symmetrized Dalitz plot constructed following the approach described in Ref. [62] shows no evidence for resonant contributions either.
The polarization of the φ mesons is studied by means of the angle θ between the direction of flight of a φ meson in the B 0 s rest frame and the B 0 s direction in the laboratory frame. With the limited sample of B 0 s → φφφ candidates the 3D fit technique to remove contributions from K + K − combinations that are not from φ decays cannot be used for this measurement. Instead, all φ mesons contributing in the mass range of the B 0 s are used, with an estimated signal purity of 71%. Figure 11 compares the cos(θ ) distribution for the B 0 s → φφφ signal candidates in data with expectations from simulation using different assumptions for the polarization. The purely longitudinal polarization clearly does not describe the data. The difference between the expectations for no polarization and purely transverse polarization is used to estimate the corresponding systematic uncertainty in the B(B 0 s → φφφ) measurement. The most probable value for the fraction of transverse polarization, f T , is found to be f T = 0.86. Assuming a uniform prior in the physically allowed range, a Bayesian lower limit of f T > 0.28 at 95% CL is found.

Summary and discussion
Charmonium production in b-hadron inclusive decays is studied in p p collisions collected at √ s = 7 and 8 TeV corresponding to an integrated luminosity of 3.0 fb −1 , using charmonium decays to φ-meson pairs. The masses and natural widths of the η c and χ c states are determined. In addition, the first evidence of B 0 s → φφφ decay is obtained. Ratios of charmonium C production rates, where the first uncertainties are statistical and the second ones are systematic. Using the branching fractions of χ c decays to φφ from Ref. [14], relative branching fractions of b hadrons decaying inclusively to χ c states are derived, = 0.92 ± 0.20 ± 0.02 ± 0.14 B , where the third uncertainty is due to the branching fractions B(χ c → φφ). These results are consistent with the ratio of the χ c1 and χ c2 production rates measured in B 0 and B + decays [14]. Inclusive production rates of the χ c states in b-hadron decays are derived using branching fractions of the χ c decays to φφ from Ref. [14], an average of the results from Belle [46] and BaBar [47] B(η c (1S) → φφ) = (3.21 ± 0.72) × 10 −3 , and the η c (1S) inclusive production rate measured using decays to p p, B(b → η c (1S)X ) = (4.88 ± 0.97) × 10 −3 [15]. They are where the third uncertainty is due to the uncertainties on the branching fraction of the b-hadron decays to the η c (1S) meson, B(b → η c (1S)X ), and η c (1S) and χ c decays to φφ. No indirect contribution to the production rate is subtracted. However, since contributions from ψ(2S) decays to the χ c states are limited, the results disfavour dominance of either colour-octet or colour-singlet contributions. The observed relations between the χ c branching fractions are not consistent with those predicted in Ref. [18]. The branching fraction B(b → χ c0 X ) is measured for the first time. The result for b-hadron decays into χ c1 is the most precise measurement for the mixture of B 0 , B + , B 0 s , B + c and b-baryons. The central value of the result for b-hadron decays into χ c1 is lower than the central values measured by the DELPHI [8] and L3 [9] experiments at LEP. The value obtained is consistent with the branching fraction of b-hadron decays into χ c1 measured by CLEO [2], Belle [4] and BaBar [5] with the light mixture of B 0 and B + . The branching fraction of b-hadron decays into χ c2 is measured for the first time with the B 0 , B + , B 0 s and b-baryons mixture. The result is consistent with the world average corresponding to the B 0 , B + mixture [14] and with individual measurements from CLEO [3], Belle [4], and BaBar [5].
Scaled differential charmonium production cross-sections as a function of p T are presented for the η c (1S) and χ c states in the LHCb acceptance and for p T > 4 GeV. Nextto-leading-order calculations of the p T dependence of the η c and χ c production rates in b-hadron decays will help to relate the results to conclusions on production mechanisms.
The production rate of the η c (2S) state in b-hadron decays is determined to be This is the first measurement for inclusive η c (2S) production rate in b-hadron decays and the first evidence for the decay η c (2S) → φφ. The production rate as a function of the assumed natural width is given in Fig. 5. These are the first χ c and η c (2S) inclusive production measurements, using charmonium decays to a hadronic final state, in the highmultiplicity environment of a hadron machine. In addition, upper limits at 95% (90%) CL on the production rates of the X (3872), X (3915), and χ c2 (2P) states in b-hadron decays are obtained, Masses and natural widths of the η c and χ c states agree with the world averages. The precision of the η c (1S) mass is comparable to the precision of the world average value. The measured η c (1S) mass is in agreement with the LHCb measurement using decays to the p p final states [15].
First evidence for the transition B 0 s → φφφ is reported with a significance of 4.9 standard deviations. Its branching fraction is measured to be B(B 0 s → φφφ) = (2.15 ± 0.54 ± 0.28 ± 0.21 B ) × 10 −6 . No resonant structure is observed in the φφ invariant mass distribution. In the B 0 s → φφφ decay, transverse polarization is preferred for the φ mesons, with an estimate of f T > 0.28 at 95% CL and the most probable value of f T = 0.86 for the fraction of transverse polarization.
As a by-product of the analysis, the branching fraction B(B 0 s → φφ) is determined to be B(B 0 s → φφ) = (2.18±0.17±0.11±0.14 f s ±0.65 B )×10 −5 with a different technique with respect to the previous results [21,22,63,64]. This technique is based on relation of B 0 s production in p p collisions and η c (1S) inclusive production rate in b-hadron decays, and reconstruction of B 0 s and η c (1S) via decays to φφ. The measurement is consistent with the recent LHCb result [22] and the current world average [14], as well as with theoretical calculations [29,30,65].
Finally, using the measurements presented and external input, the ratio of the branching fractions for the η c (1S) decays to φφ and to p p is determined. The measured B 0 s and η c (1S) yields and efficiency ratio, the branching fraction B(B 0 s → φφ) = (1.84 ± 0.05 ± 0.07 ± 0.11 f s / f d ± 0.12 norm ) × 10 −5 [22], the J/ψ production rate in b-hadron decays B(b → J/ψX ) = (1.16 ± 0.10)% [14], the relative production rates of η c (1S) and J/ψ in b-hadron decays B(b→η c (1S)X )×B(η c (1S)→ p p) B(b→J/ψX)×B(J/ψ→ p p) = 0.302 ± 0.042 [15], the branching fraction B(J/ψ → p p) = (2.120 ± 0.029) × 10 −3 [14], the ratio of fragmentation fractions f s / f d = 0.259 ± 0.015 [66], and the 0 b fragmentation fraction f 0 b momentum dependence from Ref. [67] are used. The ratio of the branching fractions for the η c (1S) decays to φφ and to p p is determined as B(η c (1S) → φφ) B(η c (1S) → p p) = 1.79 ± 0.14 ± 0.09 ± 0.10 f s / f d where the third uncertainty is related to f s / f d , the fourth uncertainty is related to f 0 b , and the fifth uncertainty is related to uncertainties of the production rates and decay branching fractions involved. This value is larger than the value computed from the world average branching fractions given in Ref. [14].