1 Introduction

The \(B^{+} \to p \bar{p} K^{+}\) decayFootnote 1 offers a clean environment to study \(c\bar{c}\) states and charmonium-like mesons that decay to \(p\bar{p}\) and excited \(\bar{ \varLambda }\) baryons that decay to \(\bar{p} K^{+}\), and to search for glueballs or exotic states. The presence of \(p \bar{p}\) in the final state allows intermediate states of any quantum numbers to be studied and the existence of the charged kaon in the final state significantly enhances the signal to background ratio in the selection procedure. Measurements of intermediate charmonium-like states, such as the X(3872), are important to clarify their nature [1, 2] and to determine their partial width to \(p\bar{p}\), which is crucial to predict the production rate of these states in dedicated experiments [3]. BaBar and Belle have previously measured the \(B^{+} \to p \bar{p} K^{+}\) branching fraction, including contributions from the J/ψ and η c (1S) intermediate states [4, 5]. The data sample, corresponding to an integrated luminosity of 1.0 fb−1, collected by LHCb at \(\sqrt{s}=7~\text {TeV}\) allows the study of substructures in the \(B^{+}\to p\bar{p} K^{+}\) decays with a sample ten times larger than those available at previous experiments.

In this paper we report measurements of the ratios of branching fractions

$$ \mathcal{R}({\rm mode}) =\frac{\mathcal{B}(B^{+} \to{\rm mode}\to p\bar{p} K^{+})}{\mathcal{B}(B^{+}\to J/\psi K^{+}\to p\bar{p} K^{+})}, $$
(1)

where “mode” corresponds to the intermediate η c (1S), ψ(2S), η c (2S), χ c0(1P), h c (1P), X(3872) or X(3915) states, together with a kaon.

2 Detector and software

The LHCb detector [6] 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{\rm\,Tm}\), and three stations of silicon-strip detectors and straw drift-tubes placed downstream. The combined tracking system has momentum (p) resolution Δp/p that varies from 0.4 % at 5 GeV/c to 0.6 % at 100 GeV/c, and impact parameter resolution of 20 μm for tracks with high transverse momentum (\(p_{\rm T}\)). Charged hadrons are identified using two ring-imaging Cherenkov (RICH) detectors. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and pre-shower 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 trigger [7] consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage where candidates are fully reconstructed. The hardware trigger selects hadrons with high transverse energy in the calorimeter. The software trigger requires a two-, three- or four-track secondary vertex with a high \(p_{\rm T}\) sum of the tracks and a significant displacement from the primary pp interaction vertices (PVs). At least one track should have \(p_{\rm T}> 1.7~\text{GeV/}c\) and impact parameter (IP) χ 2 with respect to the primary interaction greater than 16. The IP χ 2 is defined as the difference between the χ 2 of the PV reconstructed with and without the considered track. A multivariate algorithm is used for the identification of secondary vertices consistent with the decay of a b hadron.

Simulated \(B^{+} \to p \bar{p} K^{+}\) decays, generated uniformly in phase space, are used to optimize the signal selection and to evaluate the ratio of the efficiencies for each considered channel with respect to the J/ψ channel. Separate samples of \(B^{+} \to J/\psi K^{+} \to p \bar{p} K^{+}\) and \(B^{+} \to\eta_{c}(1S) K^{+} \to p \bar{p} K^{+}\) decays, generated with the known angular distributions, are used to check the dependence of the efficiency ratio on the angular distribution. In the simulation, pp collisions are generated using Pythia 6.4 [8] with a specific LHCb configuration [9]. Decays of hadronic particles are described by EvtGen [10] in which final state radiation is generated by Photos [11]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [12, 13] as described in Ref. [14].

3 Candidate selection

Candidate \(B^{+}\to p\bar{p} K^{+}\) decays are reconstructed from any combination of three charged tracks with total charge of +1. The final state particles are required to have a track fit with a \(\chi ^{2}/{\rm ndf} < 3\) where ndf is the number of degrees of freedom. They must also have p>1500 MeV/c, \(p_{\rm T} > 100~\text {MeV}/c\), and IP χ 2>1 with respect to any primary vertex in the event. Particle identification (PID) requirements, based on the RICH detector information, are applied to p and \(\bar{p}\) candidates. The discriminating variables between different particle hypotheses (π, K, p) are the differences between log-likelihood values \(\Delta\ln\mathcal{L}_{\alpha\beta}\) under particle hypotheses α and β, respectively. The p and \(\bar{p}\) candidates are required to have \(\Delta \ln\mathcal{L}_{p\pi}>-5\). The reconstructed B + candidates are required to have an invariant mass in the range 5079–5579 MeV/c 2. The asymmetric invariant mass range around the nominal B + mass is designed to select also \(B^{+} \to p \bar{p} \pi^{+}\) candidates without any requirement on the PID of the kaon. The PV associated to each B + candidate is defined to be the one for which the B + candidate has the smallest IP χ 2. The B + candidate is required to have a vertex fit with a \(\chi^{2}/{\rm ndf}<12\) and a distance greater than 3 mm, a χ 2 for the flight distance greater than 500, and an IP χ 2<10 with respect to the associated PV. The maximum distance of closest approach between daughter tracks has to be less than 0.2 mm. The angle between the reconstructed momentum of the B + candidate and the B + flight direction (\(\theta_{\rm fl}\)) is required to have \(\cos\theta _{\rm fl}>0.99998\).

The reconstructed candidates that meet the above criteria are filtered using a boosted decision tree (BDT) algorithm [15]. The BDT is trained with a sample of simulated \(B^{+} \to p \bar{p} K^{+}\) signal candidates and a background sample of data candidates taken from the invariant mass sidebands in the ranges 5080–5220 MeV/c 2 and 5340–5480 MeV/c 2. The variables used by the BDT to discriminate between signal and background candidates are: the \(p_{\rm T}\) of each reconstructed track; the sum of the daughters’ p T ; the sum of the IP χ 2 of the three daughter tracks with respect to the primary vertex; the IP of the daughter, with the highest \(p_{\rm T}\), with respect to the primary vertex; the number of daughters with \(p_{\rm T} > 900~\text {GeV}/c\); the maximum distance of closest approach between any two of the B + daughter particles; the IP of the B + candidate with respect to the primary vertex; the distance between primary and secondary vertices; the \(\theta_{\rm fl}\) angle; the \(\chi^{2}/{\rm ndf}\) of the secondary vertex; a pointing variable defined as \(\frac{P\sin\theta}{P\sin\theta+ \sum_{i} p_{\rm T,i}}\), where P is the total momentum of the three-particle final state, θ is the angle between the direction of the sum of the daughter’s momentum and the direction of the flight distance of the B + and \(\sum_{i} p_{{\rm T},i}\) is the sum of the transverse momenta of the daughters; and the log likelihood difference for each daughter between the assumed PID hypothesis and the pion hypothesis. The selection criterion on the BDT response (Fig. 1) is chosen in order to have a signal to background ratio of the order of unity. This corresponds to a BDT response value of −0.11. The efficiency of the BDT selection is greater than 92 % with a background rejection greater than 86 %.

Fig. 1
figure 1

Distribution of the BDT algorithm response evaluated for background candidates from the data sidebands (red hatched area), and signal candidates from simulation (blue filled area). The dotted line (black) indicates the chosen BDT response value (Color figure online)

4 Signal yield determination

The signal yield is determined from an unbinned extended maximum likelihood fit to the invariant mass of selected \(B^{+} \to p \bar{p} K^{+}\) candidates, shown in Fig. 2(a). The signal component is parametrized as the sum of two Gaussian functions with the same mean and different widths. The background component is parametrized as a linear function. The signal yield of the charmless component is determined by performing the same fit described above to the sample of \(B^{+} \to p \bar{p} K^{+}\) candidates with \(M_{p\bar{p}} < 2.85~\text {GeV}/c^{2}\), shown in Fig. 2(b). The B + mass and widths, evaluated with the invariant mass fits to all of the \(B^{+} \to p \bar{p} K^{+}\) candidates, are compatible with the values obtained for the charmless component.

Fig. 2
figure 2

Invariant mass distribution of (a) all selected \(B^{+} \to p \bar{p} K^{+}\) candidates and (b) candidates having \(M_{p \bar {p}} < 2.85~\text {GeV}/c^{2}\). The points with error bars are the data and the solid lines are the result of the fit. The dotted lines represent the two Gaussian functions (red) and the dashed line the linear function (green) used to parametrize the signal and the background, respectively. The vertical lines (black) indicate the signal region. The two plots below the mass distributions show the pulls (Color figure online)

The signal yields for the charmonium contributions, \(B^{+}\to(c\bar{c}) K^{+} \to p\bar{p} K^{+}\), are determined by fitting the \(p\bar{p}\) invariant mass distribution of \(B^{+} \to p\bar{p} K^{+}\) candidates within the B + mass signal window, \(\vert M_{p\bar{p} K^{+}} - M_{B^{+}}\vert< 50~\text {MeV}/c^{2}\). Simulations show that no narrow structures are induced in the \(p \bar{p}\) spectrum as kinematic reflections of possible \(B^{+} \to p \bar{\varLambda} \to p \bar{p} K^{+}\) intermediate states.

An unbinned extended maximum likelihood fit to the \(p\bar{p}\) invariant mass distribution, shown in Fig. 3, is performed over the mass range 2400–4500 MeV/c 2. The signal components of the narrow resonances J/ψ, ψ(2S), h c (1P), and X(3872), whose natural widths are much smaller than the \(p\bar{p}\) invariant mass resolution, are parametrized by Gaussian functions. The signal components for the η c (1S), χ c0(1P), η c (2S), and X(3915) are parametrized by Voigtian functions.Footnote 2 Since the \(p\bar{p}\) invariant mass resolution is approximately constant in the explored range, the resolution parameters for all resonances, except the ψ(2S), are fixed to the J/ψ value (σ J/ψ =8.9±0.2 MeV/c 2). The background shape is parametrized as \(f(M)=e^{c_{1}M+ c_{2}M^{2}}\) where c 1 and c 2 are fit parameters. The J/ψ and ψ(2S) resolution parameters, the mass values of the η c (1S), J/ψ, and ψ(2S) states, and the η c (1S) natural width are left free in the fit. The masses and widths for the other signal components are fixed to the corresponding world averages [16]. The \(p\bar{p}\) invariant mass resolution, determined by the fit to the ψ(2S) is σ ψ(2S)=7.9±1.7 MeV/c 2.

Fig. 3
figure 3

Invariant mass distribution of the \(p \bar{p}\) system for \(B^{+} \to p\bar{p} K^{+}\) candidates within the B + mass signal window, \(\vert M(p\bar{p} K^{+}) - M_{B^{+}}\vert< 50~\text {MeV}/c^{2}\). The dotted lines represent the Gaussian and Voigtian functions (red) and the dashed line the smooth function (green) used to parametrize the signal and the background, respectively. The bottom plot shows the pulls (Color figure online)

The fit result is shown in Fig. 3. Figures 4 and 5 show the details of the fit result in the regions around the η c (1S) and J/ψ, η c (2S) and ψ(2S), χ c0(1P) and h c (1P), and X(3872) and X(3915) resonances. Any bias introduced by the inaccurate description of the tails of the η c (1S), J/ψ and ψ(2S) resonances is taken into account in the systematic uncertainty evaluation.

Fig. 4
figure 4

Invariant mass distribution of the \(p \bar{p}\) system in the regions around (a) the η c (1S) and J/ψ and (b) the η c (2S) and ψ(2S) states. The dotted lines represent the Gaussian and the Voigtian functions (red) and the dashed line the smooth function (green) used to parametrize the signal and the background, respectively. The two plots below the mass distribution show the pulls (Color figure online)

Fig. 5
figure 5

Invariant mass distribution of the \(p \bar{p}\) system in the regions around (a) the χ c0(1P) and h c and (b) the X(3872) and X(3915) states. The dotted lines represent the Gaussian and Voigitian functions (red) and the dashed line the smooth function (green) used to parametrize the signal and the background, respectively. The two plots below the mass distribution show the pulls (Color figure online)

The contribution of \(c\bar{c}\to p\bar{p}\) from processes other than \(B^{+} \to p \bar{p} K^{+}\) decays, denoted as “non-signal”, is estimated from a fit to the \(p \bar{p}\) mass in the B + mass sidebands 5130–5180 and 5380–5430 MeV/c 2. Except for the J/ψ mode, no evidence of a non-signal contribution is found. The non-signal contribution to the J/ψ signal yield in the B + mass window is 43±11 candidates and is subtracted from the number of J/ψ signal candidates.

The signal yields, corrected for the non-signal contribution, are reported in Table 1. For the intermediate charmonium states η c (2S), χ c0(1P), h c (1P), X(3872) and X(3915), there is no evidence of signal. The 95 % CL upper limits on the number of candidates are shown in Table 1 and are determined from the likelihood profile integrating over the nuisance parameters. Since for the X(3872) the fitted signal yield is negative, the upper limit has been calculated integrating the likelihood only in the physical region of a signal yield greater than zero.

Table 1 Signal yields for the different channels and corresponding 95 % CL upper limits for modes with less than 3σ statistical significance. For the J/ψ mode, the non-signal yield is subtracted. Uncertainties are statistical only

5 Efficiency determination

The ratio of branching fractions is calculated using

$$\begin{aligned} \mathcal{R}({\rm mode}) =&\frac{\mathcal{B}(B^{+} \to{\rm mode}\to p\bar{p} K^{+})}{\mathcal{B}(B^{+}\to J/\psi K^{+}\to p\bar{p} K^{+})} \\=& \frac {N_{\rm mode}}{N_{ J/\psi }}\times \frac{\epsilon_{ J/\psi }}{\epsilon_{\rm mode}}, \end{aligned}$$
(2)

where \(N_{\rm mode}\) and N J/ψ are the signal yields for the given mode and the reference mode, \(B^{+}\to J/\psi K^{+}\to p\bar{p} K^{+}\), and \(\epsilon_{\rm mode}/\epsilon_{ J/\psi }\) is the corresponding ratio of efficiencies. The efficiency is the product of the reconstruction, trigger, and selection efficiencies, and is estimated using simulated data samples.

Since the track multiplicity distribution for simulated events differs from that observed in data, simulated candidates are assigned a weight so that the weighted distribution reproduces the observed multiplicity distribution. The distributions of \(\Delta\ln\mathcal{L}_{K\pi}\) and \(\Delta\ln\mathcal{L}_{p\pi}\) for kaons and protons in data are obtained in bins of momentum, pseudorapidity and number of tracks from control samples of D ∗+D 0(→K π +)π + decays for kaons and Λ decays for protons, which are then used on a track-by-track basis to correct the simulation. The efficiency as a function of \(M_{p \bar{p}}\) is shown in Fig. 6. A linear fit to the efficiency distribution is performed and the efficiency ratios are determined based on the fit result.

Fig. 6
figure 6

Efficiency as a function of \(M_{p \bar{p}}\) for \(B^{+} \to p \bar{p} K^{+}\) decays. The solid line represents the linear fit to the efficiency distribution; the dashed line is the point-by-point interpolation used to estimate the systematic uncertainty

6 Systematic uncertainties

The measurements of the relative branching fractions depend on the ratios of signal yields and efficiencies with respect to the reference mode. Since the final state is the same in all cases, most of the systematic uncertainties cancel. The systematic uncertainty on the efficiency ratio, in each region of \(p \bar{p}\) invariant mass, is determined from the difference between the efficiency ratios calculated using the solid fitted line and the dashed point-by-point interpolation shown in Fig. 6. The uncertainty associated with the evaluation of the B + signal yield has been determined by varying the fit range by ±30 MeV/c 2, using a single Gaussian instead of a double Gaussian function to model the signal PDF, and using an exponential function to model the background. For each charmonium resonance the systematic uncertainty on the signal yield has been investigated by varying the B mass signal window by ±10 MeV/c 2, the signal and background shape parametrization and the subtraction of the \(c\bar{c}\) contribution from the continuum. The systematic uncertainty associated with the parametrization of the signal tails of the J/ψ, η c (1S) and ψ(2S) resonances is taken into account by taking the difference between the number of candidates in the observed distribution and the number of candidates calculated from the integral of the fit function in the range −6σ to −2.5σ. The systematic uncertainty associated with the selection procedure is estimated by changing the value of the BDT selection to −0.03, which retains 85 % of the signal with a 30 % background, and is found to be negligible. The contributions to the systematic uncertainties from the different sources are listed in Table 2. The total systematic uncertainty is determined by adding the individual contributions in quadrature.

Table 2 Relative systematic uncertainties (in %) on the relative branching fractions from different sources. The total systematic uncertainty is determined by adding the individual contributions in quadrature

7 Results

The results are summarized in Table 3 and the values of the product of branching fractions derived from our measurement using the world average values \(\mathcal{B}(B^{+} \to J/\psi K^{+}) =(1.013\pm0.034)\times10^{-3}\) and \(\mathcal{B}(J/\psi\to p\bar{p}) =(2.17\pm0.07)\times10^{-3}\) [16] are listed in Table 4. The branching fractions obtained are compatible with the world average values [16]. The upper limit on \(\mathcal{B}(B^{+} \to\chi_{c0}(1P) K^{+} \to p \bar{p} K^{+})\) is compatible with the world average \(\mathcal{B}(B^{+} \to\chi _{c0}(1P) K^{+}) \times\mathcal{B}(\chi_{c0}(1P) \to p \bar{p}) = (0.030 \pm0.004) \times 10^{-6}\) [16]. We combine our upper limit for X(3872) with the known value for \(\mathcal{B} (B^{+} \to X(3872) K^{+} ) \times\mathcal{B} (X(3872) \to J/\psi\pi^{+} \pi^{-})= (8.6 \pm0.8) \times10^{-6}\) [16] to obtain the limit

$$ \frac{\mathcal{B} (X(3872) \to p \bar{p})}{\mathcal{B} (X(3872) \to J/\psi\pi^{+} \pi^{-})}< 2.0\times10^{-3}. $$

This limit challenges some of the predictions for the molecular interpretations of the X(3872) state and is approaching the range of predictions for a conventional χ c1(2P) state [17, 18]. Using our result and the η c (2S) branching fraction \(\mathcal{B} (B^{+} \to\eta_{c}(2S) K^{+})\times \mathcal{B} (\eta_{c}(2S) \to K \bar{K} \pi) = (3.4\, ^{+2.3}_{-1.6}) \times10^{-6}\) [16], a limit of

$$ \frac{\mathcal{B} (\eta_{c}(2S) \to p \bar{p})}{\mathcal{B} (\eta_{c}(2S) \to K \bar{K} \pi)} < 3.1 \times10^{-2} $$

is obtained.

Table 3 Signal yields, efficiency ratios, ratios of branching fractions and corresponding upper limits
Table 4 Branching fractions for \(B^{+}\to({\rm mode})\to p\bar{p} K^{+}\) derived using the world average value of the \(\mathcal{B}(B^{+}\to J/\psi K^{+})\) and \(\mathcal{B}( J/\psi \to p\bar{p})\) branching fractions [16]. For the charmonium modes we compare our values to the product of the independently measured branching fractions. The first uncertainties are statistical, the second systematic in the present measurement, and the third systematic from the uncertainty on the J/ψ branching fraction

8 Summary

Based on a sample of 6951±176 \(B^{+} \to p \bar{p} K^{+}\) decays reconstructed in a data sample, corresponding to an integrated luminosity of 1.0 fb−1, collected with the LHCb detector, the following relative branching fractions are measured

An upper limit on the ratio \(\frac{\mathcal{B} (B^{+} \to X(3872) K^{+} \to p \bar{p} K^{+})}{\mathcal{B}(B^{+} \to J/\psi K^{+} \to p \bar{p} K^{+})} < 0.017\) is obtained, from which a limit of

$$ \frac{\mathcal{B} (X(3872) \to p \bar{p})}{\mathcal{B} (X(3872) \to J/\psi\pi^{+} \pi^{-})}< 2.0\times10^{-3} $$

is derived.