Observation of the $B^{0}_{s}\rightarrow J/\psi K_{{\rm S}}^{0} K^{\pm} \pi^{\mp}$ decay

Decays of the form $B^{0}_{(s)}\rightarrow J/\psi K_{{\rm S}}^{0} h^+ h^{\left(\prime\right) -}$ ($h^{(\prime)} = K, \pi$) are searched for in proton-proton collision data corresponding to an integrated luminosity of $1.0 \, {\rm fb}^{-1}$ recorded with the LHCb detector. The first observation of the $B^{0}_{s}\rightarrow J/\psi K_{{\rm S}}^{0} K^{\pm} \pi^{\mp}$ decay is reported, with significance in excess of 10 standard deviations. The $B^{0}\rightarrow J/\psi K_{{\rm S}}^{0} K^{+} K^{-}$ decay is also observed for the first time. The branching fraction of $B^{0}\rightarrow J/\psi K_{{\rm S}}^{0} \pi^{+} \pi^{-}$ is determined, to significantly better precision than previous measurements, using $B^0 \rightarrow J/\psi K_{{\rm S}}^{0}$ as a normalisation channel. Branching fractions and upper limits of the other $B^{0}_{(s)}\rightarrow J/\psi K_{{\rm S}}^{0} h^+ h^{\left(\prime\right) -}$ modes are determined relative to that of the $B^{0}\rightarrow J/\psi K_{{\rm S}}^{0} \pi^{+} \pi^{-}$ decay.


Introduction
All current experimental measurements of CP violation in the quark sector are well described by the Cabibbo-Kobayashi-Maskawa mechanism [1,2], which is embedded in the framework of the Standard Model (SM). However, it is known that the size of CP violation in the SM is not sufficient to account for the asymmetry between matter and antimatter observed in the Universe; hence, additional sources of CP violation are being searched for as manifestations of non-SM physics.
The measurement of the phase φ s ≡ −2 arg (−V ts V * tb /V cs V * cb ) associated with B 0 s -B 0 s mixing is of fundamental interest (see, e.g., Ref. [3] and references therein). To date, only the decays B 0 s → J/ψ φ [4][5][6][7][8], B 0 s → J/ψ π + π − [9, 10] and B 0 s → φφ [11] have been used to measure φ s . To maximise the sensitivity to all possible effects of non-SM physics, which might affect preferentially states with certain quantum numbers, it would be useful to study more decay processes. Decay channels involving J/ψ mesons are wellsuited for such studies since the J/ψ → µ + µ − decay provides a distinctive experimental signature and allows a good measurement of the secondary vertex position. Observation of the decay B 0 s → J/ψ π + π − π + π − , with a significant contribution from the J/ψ f 1 (1285) component, has recently been reported by LHCb [12]. There are several unflavoured mesons, including a 1 (1260), f 1 (1285), η(1405), f 1 (1420) and η(1475), that are known to decay to K 0 S K ± π ∓ [13], and that could in principle be produced in B 0 s decays together with a J/ψ meson. If such decays are observed, they could be used in future analyses to search for CP violation.
No measurements exist of the branching fractions of B 0 (s) → J/ψ K 0 S K ± π ∓ decays. The decays B 0 → J/ψ K 0 S π + π − [14][15][16] and B 0 → J/ψ K 0 S K + K − [17,18] have been previously studied, though the measurements of their branching fractions have large uncertainties. In addition to being potential sources of "feed-across" background to B 0 (s) → J/ψ K 0 S K ± π ∓ , these decays allow studies of potential exotic charmonia states. For example, the decay chain B + → X(3872)K + with X(3872) → J/ψ π + π − has been observed by several experiments [19][20][21], and it is of interest to investigate if production of the X(3872) state in B 0 decays follows the expectation from isospin symmetry. Another reported state, dubbed the X(4140), has been seen in the decay chain B + → X(4140)K + , X(4140) → J/ψ φ by some experiments [22][23][24] but not by others [25], and further experimental studies are needed to understand if the structures in the J/ψ φ system in B + → J/ψ φK + decays are of resonant nature. In addition, the relative production of an isoscalar meson in association with a J/ψ particle in B 0 and B 0 s decays can provide a measurement of the mixing angle between the 1 √ 2 uū + dd and |ss components of the meson's wavefunction [26][27][28]. Therefore studies of B 0 (s) → J/ψ K 0 S K ± π ∓ decays may provide further insights into light meson spectroscopy.
In this paper, the first measurements of B 0 and B 0 s meson decays to J/ψ K 0 S K ± π ∓ final states are reported. All J/ψ K 0 S h + h ( )− final states are included in the analysis, where h ( ) = K, π. The inclusion of charge-conjugate processes is implied throughout the paper. The J/ψ and K 0 S mesons are reconstructed through decays to µ + µ − and π + π − final states, respectively. The analysis strategy is to reconstruct the B meson decays with minimal bias on their phase-space to retain all possible resonant contributions in the relevant invariant mass distributions. If contributions from broad resonances overlap, an amplitude analysis will be necessary to resolve them. Such a study would require a dedicated analysis to follow the exploratory work reported here.
This paper is organised as follows. An introduction to the LHCb detector and the data sample used in the analysis is given in Sec. 2, followed by an overview of the analysis procedure in Sec. 3. The selection algorithms and fit procedure are described in Secs. 4 and 5, respectively. In Sec. 6 the phase-space distributions of the decay modes with significant signals are shown. Sources of systematic uncertainty are discussed in Sec. 7 and the results are presented together with a summary in Sec. 8.

The LHCb detector
The analysis is based on a data sample corresponding to an integrated luminosity of 1.0 fb −1 of pp collisions at centre-of-mass energy √ s = 7 TeV recorded with the LHCb detector at CERN. The LHCb detector [29] 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 (VELO) surrounding the pp interaction region, a largearea 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 [30] placed downstream of the magnet. The combined tracking system provides a momentum measurement with relative uncertainty that varies from 0.4% at low momentum to 0.6% at 100 GeV/c, and impact parameter resolution of 20 µm for tracks with large transverse momentum, p T . Different types of charged hadrons are distinguished by information from two ring-imaging Cherenkov detectors [31]. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad 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 [32].
The trigger [33] consists of hardware and software stages. The events selected for this analysis are triggered at the hardware stage by a single muon candidate with p T > 1.48 GeV/c or a pair of muon candidates with p T product greater than (1.296 GeV/c) 2 . In the software trigger, events are initially required to have either two oppositely charged muon candidates with combined mass above 2.7 GeV/c 2 , or at least one muon candidate or one track with p T > 1.8 GeV/c with impact parameter greater than 100 µm with respect to all pp interaction vertices (PVs). In the subsequent stage, only events containing J/ψ → µ + µ − decays that are significantly displaced from the PVs are retained.
Simulated events are used to study the detector response to signal decays and to investigate potential sources of background. In the simulation, pp collisions are generated using Pythia [34] with a specific LHCb configuration [35]. Decays of hadronic particles are described by EvtGen [36], in which final state radiation is generated using Photos [37]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [38] as described in Ref. [39].

Analysis overview
The main objective of the analysis is to measure the relative branching fractions of the Since the most precise previous measurement of any of these branching fractions is B(B 0 → J/ψ K 0 π + π − ) = (10.3±3.3±1.5)×10 −4 [14], where the first uncertainty is statistical and the second is systematic, conversion of relative to absolute branching fractions would introduce large uncertainties. To alleviate this, a measurement of the branching fraction of B 0 → J/ψ K 0 S π + π − relative to that of B 0 → J/ψ K 0 S is also performed. For this measurement the optimisation of the selection criteria is performed based on simulation, whereas for the B 0 (s) → J/ψ K 0 S h + h ( )− relative branching fraction measurements, the optimisation procedure uses data. The two sets of requirements are referred to as "simulation-based" and "data-based" throughout the paper. The use of two sets of requirements is to avoid bias in the measurements, since the selection requirements for the yield of the numerator in each branching fraction ratio are optimised on independent samples. Furthermore, the regions of the invariant mass distributions potentially containing previously unobserved decays were not inspected until after all analysis procedures were established.
The relative branching fractions are determined from where represents the total efficiency, including effects from acceptance, trigger, reconstruction, and selection and particle identification requirements. The relative efficiencies are determined from samples of simulated events, generated with either a phase-space distribution for previously unobserved decay modes, or including known contributions from resonant structures. The relevant ratio of fragmentation fractions, denoted f d /f q , is either trivially equal to unity or is taken from previous measurements, f s /f d = 0.259 ± 0.015 [40][41][42].
The measured numbers N of decays for each channel are determined from fits to the appropriate invariant mass spectra. To determine the ratios in Eq. (2), a simultaneous fit to all J/ψ K 0 S h + h ( )− final states is used to account for possible feed-across coming from kaon-pion misidentification. The contribution from ψ(2S) decays to the J/ψ K 0 S π + π − final state is vetoed, and the veto is inverted to determine the relative branching fraction for B 0 → ψ(2S)K 0 S using a relation similar to that of Eq. (1). In Eq. (2) effects due to the width difference between mass eigenstates in the B 0 s system [43] are neglected, since the final states in B 0 s → J/ψ K 0 S h + h ( )− decays are expected to be CP mixtures. (The quantity determined using Eq. (2) is the time-integrated branching fraction.) The long lifetime of K 0 S mesons and the large boost of particles produced in LHC pp collisions causes a significant fraction of K 0 S decays to occur outside the VELO detector.
Following Refs. [44][45][46][47][48], two categories are considered: "long", where both tracks from the K 0 S → π + π − decay products contain hits in the VELO, and "downstream", where neither does. The long candidates have better momentum and vertex resolution, so different selection requirements are imposed for candidates in the two K 0 S decay categories, and the ratios given in Eqs. (1) and (2) are determined independently for each. These are then combined and the absolute branching fractions obtained by multiplying by the relevant normalisation factor. Upper limits are set for modes where no significant signal is observed In addition, the phase-space is inspected for resonant contributions from either exotic or conventional states in channels where significant signals are seen. The presence or absence of resonances could guide future analyses. However, no attempt is made to determine the relative production rates of the different possible contributions.

Selection requirements
After a set of preselection requirements to allow B candidates to be formed, additional criteria are imposed based on the output of a recursive algorithm designed to optimise the signal significance for B 0 For the measurement of the ratio of B 0 → J/ψ K 0 S π + π − and B 0 → J/ψ K 0 S branching fractions, the same requirements are also applied to B 0 → J/ψ K 0 S candidates, with the exception of those on variables that are related to the two extra pions in the numerator final state.
To optimise the simulation-based selection, used only for the determination of the relative branching fraction of B 0 → J/ψ K 0 S π + π − and B 0 → J/ψ K 0 S decays, the algorithm is applied to simulated signal events and to background events in the data. These background events are taken from invariant mass sideband regions that are not otherwise used in the analysis. For the tuning of the data-based selection, used for the relative branching fraction measurements of B 0 (s) → J/ψ K 0 S h + h ( )− decays, the properties of the B 0 → J/ψ K 0 S π + π − decays in data are used instead of simulation, since a highly pure signal can be isolated with loose requirements. Since the amount of background varies depending on whether each of h and h is a pion or kaon, different requirements are imposed for each final state. For both simulation-and data-based selections, different sets of requirements are obtained for long and downstream categories.
In the preselection, the J/ψ → µ + µ − decay is reconstructed from two oppositely charged tracks with hits in the VELO, the tracking stations and the muon chambers. The tracks are required to have p T > 500 MeV/c, to be positively identified as muons [49], to form a common vertex with χ 2 < 16, and to have an invariant mass within ±80 MeV/c 2 of the known J/ψ mass [13].
The K 0 S → π + π − decay is reconstructed from pairs of tracks with opposite charge, each with momentum greater than 2 GeV/c, that form a common vertex with χ 2 < 20. The mass of the pion pair must be within ±30 MeV/c 2 of the known K 0 S mass [13]. When considering the pair under the hypothesis that one of the tracks is a misidentified proton, the invariant mass for candidates in the long (downstream) K 0 S category must differ by more than 10 MeV/c 2 (25 MeV/c 2 ) from the known Λ baryon mass [13].
Candidates for the pions and kaons coming directly from the B decay (referred to as "bachelor" tracks) are selected if they have impact parameter χ 2 IP , defined as the difference in χ 2 of the primary pp interaction vertex reconstructed with and without the considered particle, greater than 4 and p T > 250 MeV/c. They must have momentum less than 100 GeV/c to allow reliable particle identification, and must not be identified as muons. Kaons, pions and protons are distinguished using the difference in the natural logarithm of the likelihoods (DLL) obtained from the particle identification subdetectors under the different mass hypotheses for each track [31]. Bachelor pions are selected with the requirements DLL Kπ < 0 and DLL pπ < 10, while kaons must satisfy DLL Kπ > 2 and DLL pK < 10. The particle identification efficiencies, determined from control samples of D 0 → K − π + decays reweighted to match the kinematic properties of the signal, are found to range from around 73% for B 0 The bachelor candidates are required to form a vertex with χ 2 < 10.
The B candidates are reconstructed using a kinematic fit [50] to their decay products, including the requirements that the B meson is produced at a PV and that the J/ψ and K 0 S decay products combine to the known masses of those mesons [13]. Candidates with invariant mass values between 5180 and 5500 MeV/c 2 are retained for the fits to determine the signal yields described in Sec. 5.
The recursive algorithm tunes requirements on a number of variables that are found to discriminate between signal and background and that are not strongly correlated. The most powerful variables are found to be the significance of the separation of the K 0 S vertex from the PV for the long category and the B candidate χ 2 IP . The other variables are the following: the B, J/ψ and K 0 S candidates' vertex χ 2 p-values; the J/ψ and K 0 S candidates' and the bachelor tracks' χ 2 IP values; the separation of the J/ψ vertex from the PV divided by its uncertainty; the angle between the B momentum vector and the vector from the PV to the B decay vertex; and the B candidate p T . These variables are found to be only weakly correlated with the B candidate mass or the position in the phase-space of the decay. For the simulation-based selection, the efficiency of the requirements relative to those made during preselection is around 50%. For the data-based selection the corresponding value is between around 40% for B 0 (s) → J/ψ K 0 S π + π − and around 55% for B 0 where the background is low due to the particle identification requirements and the narrow signal peak. The efficiency of the requirement that the B meson decay products lie within the detector acceptance also depends on the final state, ranging from around 10% for B 0 Backgrounds may arise from decays of b baryons. In addition to decay modes where the K 0 S meson is replaced by a Λ baryon, which are removed by the veto described above, there may be decays such as Λ 0 b → J/ψ K 0 S ph − , which have the same final state as the signal under consideration except that a kaon or pion is replaced by a proton. There is currently no measurement of such decays that could enable the level of potential background to be assessed, though the yields observed in the Λ 0 b → J/ψ pK − channel [51,52] suggest that it may not be negligible. Therefore, this background is vetoed by recalculating the candidate mass under the appropriate mass hypothesis for the final state particles and removing candidates that lie within ±25 MeV/c 2 of the known Λ 0 b mass [13].
In the J/ψ K 0 S π + π − final state, the π + π − system could potentially arise from a K 0 S meson that decays close to the B candidate vertex. This background is removed by requiring that the π + π − invariant mass is more than 25 MeV/c 2 from the known K 0 S mass [13]. In addition, There could potentially be a similar contribution in the B 0 s decay to the same final state. Such decays are removed from the sample by vetoing candidates with invariant masses of the J/ψ π + π − system within ±15 MeV/c 2 of the known ψ(2S) mass [13].
In around 2% of events retained after all criteria are applied, more than one candidate is selected. A pseudorandom algorithm is used to select only a single candidate from these events.

Determination of signal yields
After all selection requirements are applied, the only sources of candidates in the selected invariant mass ranges are expected to be signal decays, feed-across from B 0 (s) → J/ψ K 0 S h + h ( )− decays with kaon-pion misidentification, and combinatorial background. The suppression to negligible levels of other potential sources of background, such as b baryon decays, is confirmed with simulation. For each mode, the ratios of yields under the correct particle identification hypothesis and as feed-across are found to be at the few percent level from the kaon and pion control samples from D 0 → K − π + decays reweighted to the appropriate kinematic distributions. The feed-across contribution can therefore be neglected in the fit to the J/ψ K 0 S π + π − final state, as is done in the fit to the candidates passing the simulation-based selection, shown in Fig. 1.
The signal shape is parametrised in the same way for all B 0 and follows the approach used in Ref. [45]. Namely, the signal is described with the sum of two Crystal Ball functions [53] with common mean and independent tails on opposite sides of the peak. This shape is found to give an accurate description of simulated signal decays. In the fit to data, the tail parameters are fixed according to values determined from simulation. The mean and the widths as well as the relative normalisation of the two Crystal Ball functions are allowed to vary freely in the fit to data. The B 0 s region is excluded from the fit to candidates passing the simulation-based selection in the J/ψ K 0 S π + π − final state. In the fit to the J/ψ K 0 S candidates, shown in Fig. 2, a B 0 s component is included with shape identical to that for the B 0 decays except with mean value shifted by the known value of the B 0 s -B 0 mass difference [13]. The signal yields are obtained from extended unbinned maximum likelihood fits to the mass distributions of the reconstructed candidates. Independent fits are carried out for candidates in the long and downstream categories. In addition to the signal components, an exponential function is included to describe the combinatorial background with both yield and slope parameter allowed to vary freely. The results of the fits to the J/ψ K 0 S π + π − and J/ψ K 0 S invariant mass distributions are summarised in Table 1. The ratio of B 0 s → J/ψ K 0 S and B 0 → J/ψ K 0 S yields is consistent with that found in a dedicated study of those  channels [45]. Also included in Table 1 are the results of fits to the B 0 Fig. 3. These fits provide a consistency check of the analysis procedures, since the measured ratio of the B 0 → ψ(2S)K 0 S and B 0 → J/ψ K 0 S branching fractions can be compared to its known value [13]. For consistency with the fit to B 0 (s) → J/ψ K 0 S π + π − candidates, the B 0 s region is not examined in these fits. The fit to the sample selected with data-based criteria is similar to that for the sample selected with simulation-based criteria, but with some important differences. Signal shapes are included for both B 0 and B 0 s decays to each of the final states considered. The signal components are described with the same sum of two Crystal Ball functions as used in the fits  to the sample selected with simulation-based criteria, with tail parameters fixed according to values determined from simulation. For each final state, the shape of the B 0 s component is identical to that for the B 0 decays, with mean value shifted by the known value of the B 0 s -B 0 mass difference [13]. To reduce the number of freely varying parameters in the fit, the relative widths of the signal shapes in the final states with long and downstream K 0 S candidates are constrained to be identical for all signal components. The combinatorial background is modelled as a linear function, rather than the exponential model used in the fits to the samples obtained from the simulation-based selection. The use of the linear shape is found to make the fit more stable in channels with low background yields, such as B 0 (s) → J/ψ K 0 S K + K − , and it is preferable to use the same shape for all channels in the simultaneous fit. The linear function has independent parameters in each final state. A single extended unbinned maximum likelihood fit is performed for the long and downstream categories, with all final states fitted simultaneously. This procedure allows the amount of each feed-across contribution to be constrained according to the observed yields and known misidentification rates. The shapes of the feed-across contributions are described with kernel functions [54] obtained from simulation. All correlations between fitted yields are found to be less than 10% and are neglected when determining the branching fraction ratios.
The results of the fit to the samples obtained with the data-based selection are shown in Fig. 4 for the J/ψ K 0 S π + π − hypothesis, in Fig. 5 for the J/ψ K 0 S K ± π ∓ hypothesis and in Fig. 6 for the J/ψ K 0 S K + K − hypothesis. A summary of the fitted yields is given in Table 2.

Phase-space distributions of signal decays
Clear signals are seen The significance of each of the signals is discussed in Sec. 8. The distributions of the two-and three-body invariant mass combinations of the signal decay products are examined using the sPlot technique [55] with the B candidate invariant mass as the discriminating variable.
None of the channels show significant structures in any invariant mass combinations involving the J/ψ meson. In B 0 → J/ψ K 0 S π + π − decays the ψ(2S) contribution is vetoed and therefore does not appear in m(J/ψ π + π − ); there is also a small but not significant excess around the X(3872) mass. In the same channel, excesses from K * (892) and ρ(770) mesons are seen in m(K 0 S π ± ) and m(π + π − ) respectively, and there is an enhancement from the K 1 (1400) state in m(K 0 S π + π − ), as shown in Figs. 7 and 8. In B 0 s → J/ψ K 0 S K ± π ∓ decays (Figs. 9 and 10), excesses from K * (892) resonances are seen in m(K 0 S π ± ) and m(K ± π ∓ ), but no significant narrow structures are seen in m(K 0 S K ± π ∓ ).

Systematic uncertainties
Systematic uncertainties arise from possible inaccuracies in the determination of the yields, and imprecision of the knowledge of the efficiencies and fragmentation fractions that enter Eq. (1) and Eq. (2). These contributions are summarised in Tables 3 and 4 for measurements with the simulation-based and data-based selection, respectively. Total systematic uncertainties are obtained by addition in quadrature.
The systematic uncertainties on the yields are estimated by (i) varying all fixed fit parameters within their uncertainties; (ii) replacing the double Crystal Ball shape that describes the signal with a double Gaussian function; (iii) scaling the relative width of the B 0 s and B 0 peaks according to the available phase-space for the decays; (iv) replacing  Figure 7: Background-subtracted distributions of the possible two-body invariant mass combinations in B 0 → J/ψ K 0 S π + π − decays. Contributions from the ρ(770) 0 and K * (892) ± mesons are seen in the m(π + π − ) and m(K 0 S π ± ) distributions, respectively. Table 3: Systematic uncertainties (%) for the relative branching fraction measurements with B 0 → J/ψ K 0 S as normalisation channel, given separately for long and downstream categories. The total systematic uncertainty is the sum in quadrature of all contributions.

Source
Total Normalisation Yield Efficiency systematic sample size long the function that describes the combinatorial background with a second-order polynomial shape. The changes in the fitted yields are assigned as the corresponding uncertainties.
In addition, for channels where both signal and background yields are low, a small bias (less than 20% of the statistical uncertainty) on the signal yield is observed in samples of pseudoexperiments. To have a coherent treatment of all channels, each fitted yield is corrected for the bias, and the uncertainty on the bias combined in quadrature with half the correction is assigned as a systematic uncertainty. One source of systematic uncertainty that affects the relative efficiencies arises from the particle identification requirements. This is estimated by applying the method to determine the efficiency from control samples of D 0 → K − π + decays to simulated signal events, and comparing the result to the true value. The systematic uncertainty due to the variation of the efficiency over the phase-space is evaluated by reweighting the simulated samples for each signal decay to match the main features of the distributions seen in data (see Sec. 6). However, this method can only be applied for the channels where significant signals are observed: For the other decay channels the root-mean-square variation of the efficiency over the phase-space is obtained by binning the simulated events in each invariant mass combination, and this value is assigned as the associated uncertainty. There is also a small uncertainty arising from the limited simulation sample sizes.
The effective lifetimes of B 0 s meson decays depend on the CP -admixture of the final state [56]. Since the selection efficiency depends on decay time, this in principle leads to a source of uncertainty in the measurement. The scale of the efficiency variation is ±4% for the extreme ranges of possible effective lifetime distributions. Although knowledge of the exact composition of CP -even and CP -odd states requires either a detailed amplitude analysis or a measurement of the effective lifetime, all of the J/ψ K 0 S h + h ( )− final states are expected to be approximately equal admixtures. Therefore, no systematic uncertainty is assigned due to the assumption that the effective lifetime is given by 1/Γ s , where Γ s is the mean width of the two B 0 s mass eigenstates [57]. The decay time distribution observed in data is consistent with that obtained in a sample simulated with lifetime 1/Γ s , verifying that this is a reasonable assumption.
For the relative branching fraction measurement of B 0 → J/ψ K 0 S π + π − to B 0 → J/ψ K 0 S , there are two more tracks in the former channel than the latter. Therefore additional small systematic uncertainties arise due to the limited knowledge of the track reconstruction and trigger efficiencies. Uncertainty on the ratio of fragmentation fractions f s /f d affects the measurement of the B 0 s decay branching fractions relative to that of B 0 → J/ψ K 0 S π + π − . Finally, for each relative branching fraction measurement the statistical uncertainty on the normalisation channel also contributes. To allow a straightforward evaluation of the absolute branching fractions of the modes studied with the data-based selection, this source is treated separately.

Results and conclusions
Results are obtained separately for the relative branching fractions in the long and downstream categories and then combined. The combinations are performed using the full likelihood functions, though the uncertainties are symmetrised for presentation of the results. Possible correlations between systematic uncertainties in the different categories, due to the fit model, particle identification efficiencies and f s /f d , are accounted for in the combinations. All pairs of results in long and downstream categories are consistent within 2.5 standard deviations (σ). The signal significances are obtained from the change in likelihood when the signal yields are fixed to zero. Systematic uncertainties that affect the yield are accounted for in the calculation by smearing the likelihood with a Gaussian function of appropriate width. The significances are summarised in Table 5. Since the significances of the B 0 → J/ψ K 0 S K + K − and B 0 s → J/ψ K 0 S K ± π ∓ signals exceed 5 σ, these results constitute the first observations of those decays.
The results from the simulation-based selection are and where the first uncertainties are statistical and the second systematic. The measurement of B(B 0 → J/ψ K 0 S π + π − ) excludes the contribution from ψ(2S) → J/ψ π + π − decays. These results are converted to absolute branching fraction measurements using known values of the other branching fractions involved in the ratios [13]. For consistency with the standard convention, the results for the absolute branching fractions are multiplied by a factor of two to obtain values corresponding to a K 0 , instead of K 0 S , meson in the final state, giving B(B 0 → J/ψ K 0 π + π − ) = (43.0 ± 3.0 (stat) ± 3.3 (syst) ± 1.6 (PDG)) × 10 −5 , B(B 0 → ψ(2S)K 0 ) = (4.7 ± 0.7 (stat) ± 0.4 (syst) ± 0.6 (PDG)) × 10 −4 , where the last uncertainty is from knowledge of the normalisation branching fractions. These results are consistent with previous measurements [13] and, in the case of the former, significantly more precise. Table 4: Systematic uncertainties (%) for the relative branching fraction measurements with B 0 → J/ψ K 0 S π + π − as normalisation channel, given separately for long and downstream categories. The total systematic uncertainty is the sum in quadrature of all contributions.
These results are converted to absolute branching fraction measurements by multiplying by the value of the normalisation channel branching fraction determined with the simulationbased selection. In this process, the statistical uncertainty of the B 0 → J/ψ K 0 S π + π − yield is taken to be 100% correlated between the samples with simulation-based and data-based selection, since differences are small enough to be neglected. The results are B(B 0 → J/ψ ( ) K 0 K ± π ∓ ) = (11 ± 5 (stat) ± 3 (syst) ± 1 (PDG)) × 10 −6 , < 21 × 10 −6 at 90% CL , < 24 × 10 −6 at 95% CL , B(B 0 → J/ψ K 0 K + K − ) = (20.2 ± 4.3 (stat) ± 1.7 (syst) ± 0.8 (PDG)) × 10 K 0 K ± π ∓ ) = (91 ± 6 (stat) ± 6 (syst) ± 3 (f s /f d ) ± 3 (PDG)) × 10 −5 , B(B 0 s → J/ψ K 0 K + K − ) = (5 ± 9 (stat) ± 2 (syst) ± 1 (f s /f d )) × 10 −6 , < 12 × 10 −6 at 90% CL , < 14 × 10 −6 at 95% CL , where the contribution from the PDG uncertainty to the last result is negligible. The expression B(B 0 (s) → J/ψ ( ) K 0 K ± π ∓ ) denotes the sum of the branching fractions for B 0 (s) → J/ψ K 0 K − π + and B 0 (s) → J/ψ K 0 K + π − decays. In all results the strangeness of the produced neutral kaon is assumed to be that which is least suppressed in the SM. In summary, using a data sample corresponding to an integrated luminosity of 1.0 fb −1 of pp collisions at centre-of-mass energy √ s = 7 TeV recorded with the LHCb detector at CERN, searches for the decay modes B 0 (s) → J/ψ K 0 S h + h ( )− have been performed. The most precise measurement to date of the B 0 → J/ψ K 0 S π + π − branching fraction and the first observations of the B 0 → J/ψ K 0 S K + K − and B 0 s → J/ψ K 0 S K ± π ∓ decays are reported. The first limits on the branching fractions of B 0 s → J/ψ K 0 S π + π − , B 0 → J/ψ K 0 S K ± π ∓ and B 0 s → J/ψ K 0 S K + K − decays are set. Inspection of the phase-space distributions of the decays with significant signals does not reveal any potentially exotic narrow structure, nor is any significant excess from a narrow resonance seen in the K 0 S K ± π ∓ invariant mass distribution in B 0 s → J/ψ K 0 S K ± π ∓ decays. Further studies will be needed to investigate the underlying dynamics of these channels, and to understand whether they can in future be used for CP violation studies.