Study of $B_c^+$ meson decays to charmonia plus multihadron final states

Four decay modes of the $B_c^+$ meson into a $J/\psi$ meson and multiple charged kaons or pions are studied using proton-proton collision data, collected with the~LHCb detector at centre-of-mass energies of 7, 8, and 13~TeV and corresponding to an integrated luminosity of $9$~fb$^{-1}$. The decay $B_c^+\to J/\psi K^+ K^- \pi^+ \pi^+ \pi^-$ is observed for the first time, and evidence for the $B_c^+\to J/\psi 4\pi^+ 3\pi^-$ decay is found. The decay $B_c^+\to J/\psi 3\pi^+ 2\pi^-$ is observedand and the previous observation of the $B_c^+\to\psi(2S) \pi^+ \pi^+ \pi^-$ decay is confirmed using the $\psi(2S) \to J/\psi \pi^+ \pi^-$ decay mode. Ratios of the branching fractions of these four $B_c^+$ decay channels are measured.


Introduction
The B + c meson, discovered in 1998 by the CDF collaboration [1,2] at the Tevatron pp collider, is the only known meson that contains two different heavy-flavour quarks, charm and beauty.The high b-quark production cross-section at the Large Hadron

Detector and simulation
The LHCb detector [53,54] 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 [55], 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 [56,57] placed downstream of the magnet.
The tracking system provides a measurement of the momentum of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c.The momentum scale is calibrated using samples of J/ψ → µ + µ − and B + → J/ψK + decays collected concurrently with the data sample used for this analysis [58,59].The relative accuracy of this procedure is estimated to be 3 × 10 −4 using samples of other fully reconstructed b hadrons, Υ and K 0 S mesons.The minimum distance between a track and a primary pp-collision vertex (PV) [60,61], the impact parameter, 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/c.Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors (RICH) [62].Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter.Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [63].
The online event selection is performed by a trigger [64], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which performs a full event reconstruction.The hardware trigger selects muon candidates with high transverse momentum or dimuon candidates with a high value of the product of the transverse momenta of the two muons.In the software trigger, two oppositely-charged muons are required to form a good-quality vertex that is significantly displaced from any PV, and the mass of the µ + µ − pair is required to exceed 2.7 GeV/c 2 .
Simulated events are used to model the signal mass shapes and to compute the efficiencies needed to determine the branching fraction ratios.In the simulation, pp collisions are generated using Pythia [65] with a specific LHCb configuration [66].Decays of unstable particles are described by the EvtGen package [67], in which final-state radiation is generated using Photos [68].The decay channels in this study are simulated using the BLL model [49,69].The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [70] as described in Ref. [71].To account for imperfections in the simulation of charged-particle reconstruction, the track-reconstruction efficiency determined from simulation is corrected using calibration samples [72].

Decay
Yield power-law tails on both sides of the distribution [79,80].The tail parameters are fixed to the values obtained from simulation; 2. contributions from non-resonant B + c → (J/ψπ + π − ) NR π + π + π − decays, not proceeding through the intermediate ψ(2S) state, but falling into the 3.67 < m J/ψπ + π − < 3.70 GeV/c 2 region, parameterised as the product of the B + c signal function and a phase-space function describing a three-body out of the six-body final state [81], modified by a positive linear function of the J/ψπ + π − mass; 3. random combinations for ψ(2S) and π + π + π − candidates, parameterised as the product of the ψ(2S) signal function and a positive linear function of the mass of the J/ψ3π + 2π − system; 4. random J/ψ3π + 2π − combinations, described by a two-dimensional positive-definite second-order polynomial function.
For all B + c signal functions, the peak-position parameter is shared by all decays and allowed to vary in the fit.The ratio of the mass resolutions of the B + c decays in data and simulation,  s B + c , is shared by all decay modes and is allowed to vary in the fit, to account for a discrepancy in the mass resolution between data and simulation [77,78,82].The ratio of the mass resolution of the ψ(2S) → J/ψπ + π − decays in data and simulation, s ψ(2S) = 1.048 ± 0.004, and the peak-position parameter for the ψ(2S) signal component are Gaussian constrained to the values obtained from a previous LHCb study [77].The projections of the fit are overlaid in Fig. 1 for B The signal yields obtained from the fit are listed in Table 1, along with the statistical significance estimated using Wilks' theorem [83].The resolution correction factors are found to be s B + c = 1.00 ± 0.06 and s ψ(2S) = 1.048 ± 0.004.For all previously unobserved modes, the significance is confirmed by simulating a large number of pseudoexperiments according to the background distribution observed in data.
The background-subtracted mass spectra for the light-hadron system for the observed decays of the B + c mesons are obtained using the sPlot technique [84], based on the results of the fit described above.The distributions are shown in Figs. 3 and 4    The background-subtracted π + π + π − mass distribution from the B + c → J/ψ3π + 2π − decays is shown in Fig. 4 (left).The observed spectrum is in good agreement with the expectations from the BLL model.The background-subtracted π + π − mass spectra from the B + c → J/ψ3π + 2π − and B + c → ψ(2S)π + π + π − decays are shown in Fig. 5. Figures 4 and 5 contain all possible π + π + π − and π + π − combinations from a single B + c candidate.The fits to the π + π − mass distributions are performed using a function that contains two terms: a component corresponding to decays via the intermediate ρ 0 → π + π − resonance and a smooth function describing the π + π − mass spectrum without a ρ 0 → π + π − signal, labelled as "non-resonant" in Fig. 5.The resonance component is parameterised with a relativistic P-wave Breit-Wigner function with a Blatt-Weisskopf form factor with a meson radius of 3.5 GeV −1 [85].The non-resonant component is parameterised with the product of the phase-space function describing a two-body combination from a six-body combination in the B + c → J/ψ3π + 2π − case and a two-body combination from a four-body ) combination in the B + c → ψ(2S)π + π + π − case [81], and a positive first-order polynomial function that accounts for the unknown decay dynamics.The results of the fits, overlaid in Fig. 5, are consistent with a large fraction of the decays proceeding via an intermediate ρ 0 → π + π − resonance, as expected within the BLL model.Making a more quantitative statement would require a more complicated treatment of the multihadron system, which is beyond the scope of this paper.
The background-subtracted K + π − and K − π + mass spectra and the low-mass part of the K + K − mass spectrum from the B + c → J/ψK + K − π + π + π − decays are shown in Fig. 6.A fit to the K ± π ∓ mass spectrum is performed using a two-component function, similar to the function described above, and consisting of a component corresponding to decays via the intermediate K * 0 or K * 0 resonance and a smooth function describing decays without a K * 0 or K * 0 resonance.The resonance component is parameterised with a relativistic P-wave Breit-Wigner function.Fit results are overlaid in Fig. 6 (left) and indicate a presence of decays via intermediate K * 0 and K * 0 mesons.The K + K − mass spectrum, shown in Fig. 6 (right), exhibits no sign of the ϕ resonance, in agreement both with the expected suppression of the ϕ meson production due to the Okubo-Zweig-Iizuka rule [86][87][88][89][90] and with expectations from the BLL model.A similar suppression has been observed for the B + c → J/ψK + K − π + decays [14,35].
Each ratio of branching fractions for the decays of B + c mesons into the final states X and Y is calculated as R where N is the signal yield reported in Table 1 and ε denotes the corresponding efficiency.The efficiency is defined as the product of geometric acceptance and of reconstruction, selection, hadron-identification and trigger efficiencies.All of these contributions, except that of the hadron-identification efficiency, are determined using simulated samples, corrected as described in Sec. 2. The hadron-identification efficiency is calculated separately for each hadron track [62], determined from large calibration samples of [91].The measured ratios of branching fractions are where uncertainties are statistical only and correlation coefficients are listed in Table A.1.

Systematic uncertainties
The decay channels under study have similar kinematics and topologies, therefore, many sources of systematic uncertainty cancel in the branching fraction ratios, R X Y .The remaining contributions to the systematic uncertainty are summarised in Table 2 and are discussed below.
An important source of systematic uncertainty on the ratios is the imperfect knowledge of the shapes of signal and background components used in the fits.To estimate this uncertainty, several alternative models are tested.For the B + c and ψ(2S) signal shapes, a generalized Student's t-distribution [92,93] and a modified Apollonios function [94] are employed as an alternative model.For the background components, the degree of the polynomials used in the fits is increased by one.Also, the product of an exponential function and a first-order polynomial function is considered as an alternative background shape.The systematic uncertainty related to the fit model is estimated with large ensembles of pseudoexperiments.For each alternative model an ensemble of pseudoexperiments is generated and each pseudoexperiment is fitted with the baseline model.The maximal deviations in the ratios of the mean values of signal yields over the ensemble with respect to the baseline model do not exceed 2.5% for the variations of the signal model and 1.0% for the variations of background model, and are taken as systematic uncertainties.The sample of B + c → J/ψ3π + 2π − decays is used to assess the systematic uncertainty due to the procedure of multiple candidate exclusion, if two or more B + c candidates are found from the same pp collision.A large set of pseudoexperiments is performed with a random rejection of multiple candidates.The variation of the signal yield for the B + c → J/ψ3π + 2π − channel between the pseudoexperiments is found to be of 1.3% and this value is assigned as the corresponding systematic uncertainty.
To assess the systematic uncertainty related to the B + c decay model used in the simulation [49,69], the reconstructed mass distributions of the light-hadron systems in simulation are adjusted to reproduce the distributions observed in data.The uncertainty associated with the low yield of the target data distributions is accounted for by varying them within their uncertainties.The changes in the ratios R X Y do not exceed 5.1% and are taken as systematic uncertainties related to the B + c decay model.An additional uncertainty arises from the difference between data and simulation in the reconstruction efficiency of charged-particle tracks.The track-finding efficiencies obtained from simulation are corrected using data calibration samples [72].The uncertainties related to the correction factors, together with the uncertainty in the hadron-identification efficiency due to the finite size of the calibration samples [62,91], are propagated to the ratio of total efficiencies using pseudoexperiments.The obtained systematic uncertainty for the R X Y ratios does not exceed 1.1%.The hadronic interaction length of the detector is known with 10% uncertainty [95].It corresponds to an additional uncertainty for the track-finding efficiency of 1.1% (1.4%) per charged kaon (pion) track [72,95,96].This uncertainty is assumed to be totally correlated and partly cancels for the ratios.The systematic uncertainty of 1.1% related to the trigger efficiency is estimated by comparing the ratios of trigger efficiencies in data and simulation using large samples of B + → J/ψK + and B + → ψ(2S)K + decays [97].Another source of uncertainty is a potential disagreement between data and simulation in the estimation of efficiencies, due to possible effects not explicitly considered above.This is studied by varying the selection criteria of the high yield B + c → J/ψ3π + 2π − data sample in ranges that lead up to a ±20% change in the measured signal yields.The resulting difference between the efficiencies estimated using data and simulation does not exceed 2.3%, which is taken as a systematic uncertainty for the ratios R X Y .The last systematic uncertainty considered is due to the finite size of the simulated samples, and it varies between 1.5% and 2.4%.The total systematic uncertainty is estimated as the quadratic sum of individual contributions.For each choice of the alternative fit model the statistical significance for the channels under study is recalculated from data using Wilks' theorem [83].The smallest signifi-cances found are 9.0, 5.2 and 4.7 standard deviations for the B + c → J/ψK + K − π + π + π − , B + c → (ψ(2S) → J/ψπ + π − ) π + π + π − and B + c → J/ψ4π + 3π − decays, respectively.These values are taken as the significance including systematic uncertainty.
Three ratios of branching fractions, defined in Eqs.(1), are measured as where the first uncertainty is statistical and the second systematic.Correlation coefficients for statistical and systematic uncertainties for the measured ratios of branching fractions are given in Appendix A. The mass spectra for the light-hadron system, as well as the mass spectra for the intermediate combinations of light hadrons agree with the phenomenological model by Berezhnoy, Likhoded and Luchinsky based on QCD factorisation [43][44][45][46][47][48][49].The ratio R PIC (Spain), GridPP (United Kingdom), CSCS (Switzerland), IFIN-HH (Romania), CBPF (Brazil), Polish WLCG (Poland) and NERSC (USA).We are indebted to the communities behind the multiple open-source software packages on which we depend.Individual groups or members have received support from ARC and ARDC (Australia); Minciencias (Colombia); AvH Foundation (Germany); EPLANET, Marie Sk lodowska

Figure 6 :
Figure 6: Background-subtracted (left) K ± π ∓ (3 entries per B + c candidate) mass and (right) low-mass part of the K + K − mass distribution from the B + c → J/ψK + K − π + π + π − decays.The results of the fit described in the text are overlaid on the left plot, while expectations from the BLL model are overlaid on the right plot.
(right) together with the expectations from the BLL model.For all cases, good agreement with the BLL model is observed.No D + s → 3π + 2π − , D + s → 4π + 3π − or D + s → π + π + π − signals are observed in the studied spectra.

Table 2 :
Ranges of relative systematic uncertainties for the various ratios of branching fractions, R X Y .The total systematic uncertainty is the quadratic sum of individual contributions.