Updated branching fraction measurements of $B^0_{(s)} \to K_{\mathrm{\scriptscriptstyle S}}^0 h^+ h^{\prime -}$ decays

The charmless three-body decays $B^0_{(s)} \to K_{\mathrm{\scriptscriptstyle S}}^0 h^+ h^{\prime -}$ (where $h^{(\prime)} = \pi, K$) are analysed using a sample of $pp$ collision data recorded by the LHCb experiment, corresponding to an integrated luminosity of $3\mbox{fb}^{-1}$. The branching fractions are measured relative to that of the $B^0 \to K_{\mathrm{\scriptscriptstyle S}}^0 \pi^{+} \pi^{-}$ decay, and are determined to be: \begin{eqnarray*} \frac{{\cal B}(B^0\rightarrow K^0_{\rm S}K^{\pm}\pi^{\mp})}{{\cal B}(B^0\rightarrow K^0_{\rm S}\pi^{+}\pi^{-})} = {}&0.123 \pm 0.009 \; \mathrm{\,(stat)}\; \pm 0.015 \; \mathrm{\,(syst)}\,, \frac{{\cal B}(B^0\rightarrow K^0_{\rm S}K^{+}K^{-})} {{\cal B}(B^0\rightarrow K^0_{\rm S}\pi^{+}\pi^{-})} = {}&0.549 \pm 0.018 \; \mathrm{\,(stat)}\; \pm 0.033 \; \mathrm{\,(syst)}\,, \frac{{\cal B}(B_{s}^0\rightarrow K^0_{\rm S}\pi^{+}\pi^{-})}{{\cal B}(B^0\rightarrow K^0_{\rm S}\pi^{+}\pi^{-})} = {}&0.191 \pm 0.027 \; \mathrm{\,(stat)}\; \pm 0.031 \; \mathrm{\,(syst)}\; \pm 0.011 \; (f_s/f_d) \,, \frac{{\cal B}(B_{s}^0\rightarrow K^0_{\rm S}K^{\pm}\pi^{\mp})} {{\cal B}(B^0\rightarrow K^0_{\rm S}\pi^{+}\pi^{-})} = {}&1.70\phantom{0} \pm 0.07\phantom{0} \; \mathrm{\,(stat)} \; \pm 0.11\phantom{0} \; \mathrm{\,(syst)}\; \pm 0.10\phantom{0} \; (f_s/f_d) \,, \frac{{\cal B}(B_{s}^0\rightarrow K^0_{\rm S}K^{+}K^{-})}{{\cal B}(B^0\rightarrow K^0_{\rm S}\pi^{+}\pi^{-})} \in {}&[0.008 - 0.051] \rm ~at~90\%~confidence~level, \end{eqnarray*} where $f_s/f_d$ represents the ratio of hadronisation fractions of the $B^0_s$ and $B^0$ mesons.


Updated branching fraction measurements
The LHCb collaboration †

Introduction
The measurement of CP -violation observables in the decays B 0 → K 0 S π + π − and B 0 → K 0 S K + K − , which are dominated by b → qqs (q = u, d, s) loop transitions, are of great theoretical interest. 1 In particular, the mixing-induced CP asymmetries in these decays are predicted by the Standard Model (SM) Cabibbo-Kobayashi-Maskawa mechanism [1,2] to be approximately equal to those governed by b → ccs transitions, such as B 0 → J/ψ K 0 S . Within the SM the weak phase measurements in b → qqs decays are expected to deviate from the values determined in b → ccs decays but for certain contributions to these decays, such as B 0 → φK 0 S and B 0 → ρ 0 K 0 S , this deviation is either expected to be small or can be controlled using flavour symmetries [3][4][5]. The existence of new particles predicted in several extensions of the SM could introduce additional weak phases that contribute along with the SM mixing phase to the amplitudes of these loop-dominated charmless decays, potentially leading to much greater deviations from the b → ccs values [6][7][8]. The mixing-induced CP -violating phase can be measured by means of a flavour-tagged time-dependent analysis of the three-body Dalitz plot of these decays [9][10][11][12]. The current experimental measurements of this phase in b → qqs decays [13] show a generally good agreement with the results for the weak phase β from b → ccs decays for each of the CP eigenstates studied. The experimental uncertainties are, however, currently rather larger than the size of the expected deviations, both in the SM and beyond-the-SM scenarios, and so there is a need for more precise measurements of these quantities. A similar determination of the mixing-induced CP -violating phase in the B 0 s system is possible with, among others, the B 0 s → K 0 S K ± π ∓ decays [14]. It is also possible to determine the CKM angle γ by combining information from several B → Khh decays, using either the methods originally proposed in Refs. [15,16] and recently developed further in Ref. [17], or those proposed in Refs. [18][19][20]. The existing experimental results, which come from the BaBar collaboration [21,22], demonstrate the feasibility of the measurement, albeit with large statistical uncertainties. The decay B 0 s → K 0 S π + π − is dominated by tree-level processes and as such is of particular interest for this effort, with the potential to yield a theoretically clean determination of γ [23].
The measurements of the branching fractions themselves are of great importance in order to confront theoretical predictions. These predictions are based on various approaches to modelling the hadronisation processes, such as QCD factorisation or PQCD, see for example Refs. [24][25][26][27][28][29]. Comparison of the different approaches with the experimental data will allow further refinement of the theoretical models, which in turn will yield improved predictions of branching fractions and CP asymmetries of these and many other charmless decay modes. In addition, these results can be used to test the level of breaking of the flavour symmetries: isospin, U-spin and SU(3), see for example Ref. [30].
Of the decays of neutral B mesons to K 0 S π + π − , K 0 S K ± π ∓ and K 0 S K + K − final states, only the decay B 0 s → K 0 S K + K − remains to be observed [10,12,[31][32][33][34]. Most recently, a search for the three B 0 s decays was reported by the LHCb experiment using the 1 fb −1 data sample recorded in 2011 [34]. While first observations were made for the B 0 s → K 0 S π + π − and B 0 s → K 0 S K ± π ∓ modes, no evidence for the decay B 0 s → K 0 S K + K − was found. In this work, all the aforementioned charmless three-body decays of the B 0 and B 0 s mesons are studied using the pp collision data recorded by the LHCb detector, corresponding to an integrated luminosity of 1.0 fb −1 at a centre-of-mass energy of 7 TeV in 2011 and 2.0 fb −1 at a centre-of-mass energy of 8 TeV in 2012. This sample is three times larger than that used in Ref. [34]. The measurements of the time-integrated branching fractions [35] relative to that of B 0 → K 0 S π + π − are presented. The notation B(B 0 → K 0 S K ± π ∓ ) is used throughout the document to indicate the sum of the branching fractions B(B 0 → K 0 S K + π − ) and B(B 0 → K 0 S K − π + ), and similarly for the corresponding B 0 s decays.

Detector and simulation
The LHCb detector [36,37] 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 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 to 1.0% at 200 GeV/c. 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/c. 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 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. Simulated data samples are used to investigate backgrounds from other b-hadron decays and also to study the detection and reconstruction efficiency of the signal. In the simulation, pp collisions are generated using Pythia [38] with a specific LHCb configuration [39]. Decays of hadronic particles are described by EvtGen [40], in which final-state radiation is generated using Photos [41]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [42] as described in Ref. [43].

Trigger and event selection
The online event selection is performed by a trigger [44], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, in which all charged particles with p T > 500 (300) MeV/c are reconstructed for data collected in 2011 (2012). At the hardware trigger stage, events are required to have a muon with high p T or a hadron, photon or electron with high transverse energy in the calorimeters. The software trigger requires a two-, three-or four-track secondary vertex with a significant displacement from all primary pp interaction vertices. At least one charged particle must have transverse momentum p T > 1.7 (1.6) GeV/c in the 2011 (2012) data and be inconsistent with originating from a PV. A multivariate algorithm [45] is used for the identification of secondary vertices consistent with the decay of a b hadron.
It is required that the software trigger decision must have been caused entirely by tracks from the decay of the signal B candidate.
To suppress the 'combinatorial' background formed by combinations of unrelated tracks, the events satisfying the trigger requirements are filtered in two stages: a preselection based on loose requirements, followed by a multivariate selection. In order to minimise the variation of the selection efficiency over the Dalitz plot, the selection procedure uses only loose requirements on the momenta of the B-meson decay products and relies mainly on topological features such as the flight distance of the B candidate. These features depend on whether the B candidate or the K 0 S , h + , h − candidates are consistent with having originated from a particular PV. It is therefore necessary to 'associate' each candidate with a single PV -that from which it is most consistent with having originated. The association is defined in terms of the χ 2 IP quantity, which is the difference in fit χ 2 of the given PV reconstructed with and without the track or tracks from the particle in question. In events that contain more than one PV, each candidate is associated with the PV that has the smallest χ 2 IP . Decays of K 0 S → π + π − are reconstructed in two different categories: the first involving K 0 S mesons that decay early enough for the resulting pions to be reconstructed in the VELO; and the second containing those K 0 S mesons that decay later, such that track segments of the pions cannot be formed in the VELO. These K 0 S reconstruction categories are referred to as long and downstream, respectively. The long category has better mass, momentum and vertex resolution than the downstream category. There are however approximately twice as many K 0 S candidates reconstructed as downstream than as long, simply due to the lifetime of the K 0 S meson and the geometry of the detector. In the following, B candidates reconstructed from either a long or downstream K 0 S candidate, in addition to two oppositely charged tracks, are also referred to using these category names. During the 2012 data taking, a significant improvement of the trigger efficiency for long-lived particles, specifically for downstream candidates, was obtained following an update of the software trigger algorithms. To take into account the differences in trigger efficiencies and the different data-taking conditions, the data sample is divided into 2011, 2012a, and 2012b data-taking periods, and each period is divided in two sub-samples according to the K 0 S reconstruction category. The 2012b sample is the largest, corresponding to 1.4 fb −1 , and also has the highest trigger efficiency.
The two charged pions that form the K 0 S candidates are both required to have momentum p > 2 GeV/c and have χ 2 IP with respect to their associated PV greater than 9 (4) for long (downstream) candidates. They are then required to form a vertex with good fit quality (quantified by the fit χ 2 , χ 2 vtx < 12) and to have invariant mass within 20 MeV/c 2 (30 MeV/c 2 ) of the nominal K 0 S mass [46] for long (downstream) candidates. A requirement on the square of the ratio of the separation of the K 0 S vertex from its associated PV and the corresponding uncertainty, χ 2 VS > 80 (50) for long (downstream) candidates, ensures a significant vertex separation. Downstream K 0 S candidates are required in addition to have a momentum p > 6 GeV/c.
The B candidates are formed from a K 0 S candidate and two oppositely charged tracks (initially reconstructed under the pion mass hypothesis). Each of these two tracks is required to have p < 100 GeV/c, a value beyond which there is little pion-kaon discrimination. The scalar sum of the transverse momenta of the K 0 S and the two h + h − candidates must be greater than 3.0 GeV/c (4.2 GeV/c), for long (downstream) candidates, and at least two of the three decay products must have p T > 0.8 GeV/c. The IP of the B-meson decay product with the largest p T is required to be greater than 0.05 mm relative to the PV associated to the B candidate. The B candidate decay products are then required to form a vertex that has χ 2 vtx < 12 and which is separated from any PV by at least 1.7 mm. The difference in χ 2 vtx when adding another track must be greater than 4. The B candidates must have p T > 1.5 GeV/c and invariant mass within the range 4000 < m K 0 S π + π − < 6200 MeV/c 2 . They are further obliged to be consistent with originating from a PV, quantified by requiring, for long (downstream) candidates, both that χ 2 IP < 8 (6) and that the cosine of the angle θ DIR between the reconstructed momentum of the B candidate and the vector between the associated PV and the decay vertex be greater than 0.9999 (0.999). Finally, the decay vertex of the K 0 S candidate is required to be at least 30 mm downstream, along the beam direction, from that of the B candidate.
Multivariate discriminants based on a boosted decision tree (BDT) algorithm [47,48] are used to further reduce combinatorial backgrounds. Simulated B 0 → K 0 S π + π − events and data from upper mass sidebands, 5425 < m K 0 S π + π − < 6200 MeV/c 2 , are used as the signal and background training samples, respectively. Contributions from muons and protons are removed from these samples using particle identification (PID) variables. Each of the six samples (resulting from the division by the three data-taking periods and the two K 0 S reconstruction categories) is further subdivided into two equally-sized subsamples. Each subsample is then used to train an independent discriminant. In the subsequent analysis the BDT trained on one subsample of a given category is used to select events from the other subsample, in order to avoid bias. The input quantities for the BDTs are: the p T , η, χ 2 IP , χ 2 VS , cos θ DIR and χ 2 vtx values of the B candidate; the smallest change in the B-candidate χ 2 vtx value when adding another track from the event; the sum of the χ 2 IP values of the h + and h − candidates; the χ 2 IP , χ 2 VS and χ 2 vtx values of the K 0 S candidate; and the p T asymmetry where p T cone is the transverse component of the sum of all particle momenta inside a cone around the B-candidate direction, of radius R ≡ δη 2 + δφ 2 = 1.5, where δη and δφ are the difference in pseudorapidity and azimuthal angle (in radians) around the beam direction, between the momentum vector of the track under consideration and that of the B candidate.
The selection requirement placed on the output of the BDTs is independently optimised for each data sample. For all signal decay modes that have previously been observed, the following figure of merit is used where N sig (N bg ) represents the number of expected signal (combinatorial background) events for a given selection. The value of N sig is estimated based on the known branching fractions and efficiencies, while N bg is calculated by fitting the sideband above the signal region and extrapolating into the signal region, defined as the invariant-mass window of five times the typical resolution around the B 0 and the B 0 s masses. For the yet unobserved B 0 s → K 0 S K + K − mode, an alternative figure of merit [49] is where the signal efficiency (ε sig ) is estimated from the signal simulation. The optimisation is performed separately for each of the six categories. As each final state contains both B 0 and B 0 s signals, one of which is favoured and the other suppressed, this procedure results in applying two differently optimised selections on each final state.
Particle identification requirements are subsequently applied in order to reduce backgrounds from decays such as where, respectively, the proton and muons are misidentified as pions or kaons. PID information is also used to assign each candidate exclusively to one out of four possible final states: The PID requirements are optimised to reduce the cross-feed between the different signal decay modes using the same figures of merit introduced for the BDT optimisation.
Fully reconstructed B-meson decays into two-body Dh or (cc)K 0 S combinations, where (cc) indicates a charmonium resonance, may result in a K 0 S h ± h ∓ final state that satisfies the selection criteria and has the same B-candidate invariant mass distribution as the signal candidates. The decays of Λ 0 b baryons to Λ + c h with Λ + c → pK 0 S also peak under the signal when the proton is misidentified. Therefore, the following D and Λ + c decays are explicitly reconstructed under the relevant particle hypotheses and vetoed in all the spectra: Additional vetoes on charmonium resonances, J/ψ → π + π − , K + K − and χ c0 → π + π − , K + K − , are applied to remove the small number of fully reconstructed and well identified peaking B 0 (s) → (J/ψ , χ c0 ) K 0 S decays. The vetoed region for each reconstructed charm (charmonium) state is an invariant-mass window of 30 (48) MeV/c 2 around the world average mass value of that state [46]. This range reflects the typical mass resolution obtained at LHCb.
The fraction of selected events containing more than one B candidate is at the percent level. The candidate to be retained in each event is chosen randomly, but reproducibly.

Fit model
The signal yields corresponding to each of the BDT optimisations are determined by means of a simultaneous unbinned extended maximum likelihood fit to the B-candidate invariant mass distributions of all final states in the six categories. Four types of components contribute to each invariant mass distribution: signal decays, backgrounds resulting from cross-feeds, partially reconstructed decays, and random combinations of unrelated tracks.
Signal B 0 (s) → K 0 S h ± h ∓ decays with correct identification of the final-state particles are modelled with the sum of two Crystal Ball (CB) functions [50] that share common values for the peak position and width but have independent power law tails on opposite sides of the peak. The B 0 and B 0 s masses (peak positions of the CB functions) are free parameters in the fit and are allowed to take different values in the different data-taking periods in order to allow for small differences in momentum calibration. Seven parameters related to the widths of the CB functions are also free parameters of the fit: the width of the downstream B 0 → K 0 S π + π − signal in each of the three data-taking periods; the ratio of the widths of the B 0 s and B 0 decay modes; the relative widths of K 0 S K ± π ∓ and K 0 S K + K − to K 0 S π + π − ; and the ratio of the widths in the long and downstream categories. The dependence of the width on each of these divisions is assumed to factorise; for example the width σ of the long B 0 s → K 0 S K + K − signal in the 2011 data-taking period is related to that of the downstream B 0 → K 0 S π + π − signal in the same data-taking period by σ 2011 long where r x/y indicates the ratio of the widths of categories x and y. These assumptions are made necessary by the otherwise poor determination of the width of the suppressed mode in each spectrum. The other parameters of the CB components are obtained by a simultaneous fit to simulated samples.
Cross-feed contributions from misidentified signal decays are modelled empirically by the sum of two CB functions using simulated events. Only contributions from the Other potential misidentified decays are neglected, as their contributions have been checked to be below one event. The relative yield of each misidentified decay is constrained with respect to the yield of the corresponding correctly identified decay. The constraints are implemented using Gaussian prior probability distributions included in the likelihood. The mean values are obtained from the ratio of selection efficiencies and the widths include uncertainties originating from the finite size of the simulated event samples and the systematic uncertainties related to the determination of the PID efficiencies.
Backgrounds from partially reconstructed decays such as where the neutral pion is not reconstructed, are also modelled. Four categories are included in each of the final state spectra, where the background results from either charmed or charmless decays of B 0,+ or B 0 s mesons. These decays are modelled by means of generalised ARGUS functions [51] convolved with a Gaussian resolution function. Their parameters are determined from simulated samples of the expected dominant decays in each category. Radiative decays and those from B 0 → η K 0 S are considered separately and included only in the K 0 S π + π − final state. The normalisation of all such contributions is constrained with respect to the signal in the relevant final state using Gaussian prior probability distributions based on the ratio of efficiencies and the ratio of branching fractions from world averages [46]. The relative uncertainties on these ratios vary between 20% and 100%.
The combinatorial background is modelled by a linear function. The variations of the slope parameter between data-taking periods, K 0 S reconstruction categories and the different final states are assumed to factorise (in an analogous way to the widths of the signal distributions), leaving six free parameters. This assumption, as well as the choice of the linear model, are considered as sources of systematic uncertainties.
The fit results for each BDT optimisation, combining all data-taking periods, are displayed in Figs. 1 and 2. The separate plots for the individual data-taking periods are shown in Figs. 3-14 in Appendix A. Table 1 shows the signal yields for each mode summed over all data-taking periods and K 0 S reconstruction categories, along with a weighted sum of efficiencies. The fitted yields of each decay mode for each of the three data-taking periods and two K 0 S reconstruction categories are given in Appendix A. Statistical correlations between the signal yields are below 10% in all cases and are neglected. For the suppressed modes, the combinatorial background is negligible in the high invariant-mass region for Table 1: Signal yields obtained from the simultaneous fit to the data. The yields are the sum of those obtained in the three data-taking periods when fitting the data sample selected using the BDT optimisation chosen for the given decay mode. The uncertainties are statistical only. The average selection efficiencies, described in Sec. 5, are also shown for each decay mode together with the corresponding total uncertainty due to the limited simulation sample size and systematic effects in their determination.

downstream long Decay
Yield leading to a small systematic uncertainty related to the assumptions used to fit this component. In order to determine the significance of in each fit category, taking into account systematic uncertainties. The profiles are constructed from fits where the shape parameters of the B 0 s → K 0 S K + K − signal are fixed to the values obtained from the nominal fit, which allows the change in the fit likelihood to be interpreted using Wilks' theorem [52]. Combining these profiles yields a significance of 2.5 σ.

Determination of the efficiencies
The measurements of the branching fractions of the B 0 (s) where ε sel is the selection efficiency (which includes geometrical acceptance, reconstruction, selection, trigger and particle identification components), N is the fitted signal yield, and f d and f s are the hadronisation fractions of a b quark into a B 0 and B 0 s meson, respectively. The ratio f s /f d has been precisely determined by the LHCb experiment from hadronic and semileptonic measurements to be f s /f d = 0.259 ± 0.015 [53]. Since the CP content of the three B 0 s decays is currently unknown, the calculation of the corresponding efficiencies assumes an effective lifetime of 1/Γ s , where Γ s is the average width of the two CP -eigenstates of the B 0 s meson. The effect of varying the decay width by ±∆Γ s /2, where ∆Γ s is the width difference between the two B 0 s CP -eigenstates, results in relative changes to the average efficiency of ∓4%.
Three-body decays are, in general, composed of several quasi-two-body decays and nonresonant contributions, all of them possibly interfering. The signal reconstruction,  selection and trigger efficiencies also vary over the phase space. Hence, both the distribution of the signal events and the variation of the efficiency over the Dalitz plot [54] must be determined in order to calculate the efficiency-corrected yield. In this analysis, efficiencies candidates, summing the three periods of data taking, with the selection optimisation chosen for the suppressed decay modes for (left) downstream and (right) long K 0 S reconstruction categories. Colours and line styles follow the same conventions as in Fig. 1 are determined for each decay mode from simulated signal samples in bins of the "square Dalitz plot" [55], where the usual Dalitz-plot coordinates have been transformed in order to map a rectangular space. The edges of the phase space are spread out in the square Dalitz plot, which permits a more precise modelling of the efficiency in the regions where it varies the most and where most of the signal candidates are expected. The square Dalitz-plot distribution of each signal mode is determined from the data using the sPlot technique [56], using the K 0 S h + h − invariant mass as the discriminating variable. The distributions of signal events on the Dalitz plot, as obtained using this technique, are shown in Appendix C. The efficiency-corrected yields are calculated as: where w i is the sPlot weight and ε i is the efficiency for event i, and the sum includes all events in the fitted data sample. The average efficiency for each decay mode: where is the fitted signal yield, is given for each signal decay in Table 1. They are presented for each data-taking period and K 0 S reconstruction category in Tables 3-5 in Appendix A. Their relative uncertainties due to the finite size of the simulated event samples vary from 1% to 20%.
The imperfections of the simulation are corrected for in several respects. Inaccuracies of the tracking efficiency in the simulation are mitigated by weighting the h + and h − tracks by a correction factor obtained from a data calibration sample [57]. An analogous correction is applied for the K 0 S tracking and vertex reconstruction efficiency. A control data sample of D * + → D 0 (→ K − π + )π + s decays, where π + s indicates a slow pion, is used to quantify the differences between the data and simulation hardware trigger stage for pions and kaons, independently for positive and negative hadron charges, as a function of p T [44]. Corrections to account for differences between data and simulation related to tracking efficiency are O(1%), while those related to trigger efficiency can be O(10%), depending on the position on the Dalitz plot. The uncertainties attached to these various corrections are propagated to the branching fraction measurements as systematic uncertainties and are further discussed in Sec. 6.
The particle identification efficiencies and misidentification rates are determined from simulated signal decays on an event-by-event basis using a data-driven calibration to match the kinematic properties of the tracks in the decay of interest. A weighting procedure is performed in bins of p, η and event multiplicity, accounting for kinematic correlations between the tracks. Calibration tracks are taken from D * + → D 0 π + s decays where the D 0 decays to the Cabibbo-favoured K − π + final state. In this case, the charge of the soft pion π + s provides the kaon or pion identity of the tracks. The dependence of the PID efficiency over the Dalitz plot induced by the variations of PID performance with the track kinematics is included in the procedure described above. This calibration is performed using samples from the same data-taking period, accounting for the variation in the performance of the Cherenkov detectors over time.

Systematic uncertainties
Most of the systematic uncertainties are eliminated or greatly reduced by normalising the branching fraction measurements to the B 0 → K 0 S π + π − mode. A summary of the contributions, expressed as relative uncertainties, is given in Table 2, including the uncertainty in the measurement of f s /f d [53]. A detailed breakdown of systematic uncertainties per data-taking period and K 0 S reconstruction category is provided in Table 2: Summary of the systematic uncertainties on the ratio of branching fractions. A weighted average of the two K 0 S reconstruction categories and three data-taking periods is performed. The values quoted for the individual contributions, which are illustrative of the hierarchy between sources of systematic uncertainty, each result from a weighted average in which the other systematic uncertainty contributions are disregarded. The total uncertainty is the weighted average including all contributions. All uncertainties are relative and are quoted as percentages. Tables 6 and 7 in Appendix B. The dominant contributions arise from the modelling of the combinatorial background shape in the invariant mass fit and from the determination of the efficiency of the hardware trigger.

Fit model
The fit model relies on a number of assumptions, both in terms of the values of parameters being taken from simulation and in terms of the choice of the functional forms describing the various contributions. In both cases, the uncertainties are evaluated using pseudoexperiments that are generated from the alternative parameterisation and are fitted using both the nominal and the alternative fit models. The distribution of the difference in the value of a given parameter determined in the two fit model is subsequently fitted with a Gaussian function and the corresponding systematic uncertainty is assigned as the sum in quadrature of the mean and the resolution of the Gaussian. This procedure is employed for the fixed parameters of the signal, partiallyreconstructed and cross-feed backgrounds and for the functional forms used for the signal and combinatorial background. Due to the limited sizes of the simulated event samples used to parameterise both the partially-reconstructed and cross-feed background shapes, the uncertainty related to the fixed parameters also covers any reasonable variation of the shape. For the combinatorial background, the ratios of the slopes in different K 0 S reconstruction categories and in different data-taking periods are constrained to be the same for all final states. Two alternative models are considered: allowing independent ratios for each of the final states (testing the assumption of the factorisation of the slope ratios) and using an exponential model instead of the nominal linear one (testing the functional form of the combinatorial shape).
Finally, in order to evaluate the impact of residual contributions from Λ 0 b decays that survive the proton PID veto described in Sec. 3, fits to data are performed including a model for this contribution. As these fits show negligible difference to the nominal model, no systematic uncertainty is assigned. The total fit model systematic uncertainty is given by the sum in quadrature of all the contributions and is mostly dominated by the combinatorial background model uncertainty. Some uncertainties are fully correlated among the different data samples that are averaged and are treated as such in the uncertainty propagation. The correlated fit model systematics include uncertainties due to the fit biases and combinatorial and signal shapes. The combination of these contributions is shown in Table 2 as "Fit model", while they are referred to in Tables 6 and 7 as "Fit model (corr.)" and "Fit model (uncorr.)", respectively.

Selection and trigger efficiencies
The accuracy of the efficiency determination is limited in most cases by the finite size of the samples of simulated signal events, duly propagated as a systematic uncertainty. In addition, the effect related to the choice of binning for the square Dalitz plot is estimated from the spread of the average efficiencies determined from several alternative binning schemes. These two sources of uncertainties are detailed in Tables 6 and 7, and are labelled "Selection (statistics)" and "Selection (binning)", respectively.
As introduced in Sec. 5, the sources of uncertainties related to the imperfections of the tracking simulation are two-fold: the reconstruction of both long and downstream tracks and the reconstruction of the K 0 S decay vertex (in particular for the downstream category). In both cases, the reconstruction efficiencies are determined from data calibration samples and the simulated events are weighted to match the performance measured in the data. The uncertainties arising from the finite size of the calibration samples are propagated to assign a systematic uncertainty.
Possible sources of systematic uncertainty related to the efficiency estimation of the hardware trigger have been studied. Two additional systematic uncertainties are assigned: one related to the imperfect simulation of the rate of overlapping tracks in the hadron calorimeter forming a single hadron trigger candidate and one related to the choice of the data calibration sample itself. These two sources of uncertainties are detailed in Tables 6 and 7 in Appendix B, labelled "Trigger (overlap)" and "Trigger (calib. sample)", respectively. For the first source, the systematic uncertainty is estimated as the difference between the trigger efficiency correction with and without the overlapping cluster corrections. For the latter source of uncertainty, the correction factors have been determined from a sample of reconstructed B 0 → J/ψ K + π − events. Twice the difference between the correction factors determined from the two calibration samples is taken as the estimate of the associated systematic uncertainty.
The uncertainties due to the choice of the binning of the square Dalitz plot used to produce the efficiency maps and those related to the hardware trigger calibration samples are treated as correlated among the different data samples split by year of data taking.

Particle identification efficiencies
The procedure to evaluate the efficiencies of the PID selections uses calibration tracks that differ from the signal tracks in terms of their kinematic distributions. While the binning procedure attempts to mitigate these differences, there could be some remaining systematic effect. This is addressed by considering different ensembles of kinematical binnings to determine the efficiency. An overall 1% systematic uncertainty is assigned to quantify any bias due to the procedure. The statistical uncertainties originating from the finite sample sizes are added in quadrature.

Results and conclusion
The decay modes B 0 (s) → K 0 S h ± h ∓ have been analysed using a dataset, corresponding to an integrated luminosity of 3.0 fb −1 recorded by the LHCb detector at a centre-of-mass energy of 7 TeV and 8 TeV. The branching fraction of each decay is measured relative to that of B 0 → K 0 S π + π − . The ratios of branching fractions are determined independently for the two K 0 S reconstruction categories and three data-taking periods and then combined by performing a χ 2 fit. The corresponding covariance matrix includes the statistical and systematic uncertainties. A 100% linear correlation factor is assumed for the correlated systematic uncertainties. Good agreement is obtained among all determinations from each data-taking period and K 0 S reconstruction category. The results obtained from the combination are = 0.026 ± 0.011 (stat) ± 0.007 (syst) ± 0.002 (f s /f d ) .
All measurements of branching fractions are in good agreement with the earlier LHCb determinations [34], which they supersede. The measurement of the relative branching fractions of B 0 → K 0 S K ± π ∓ and B 0 → K 0 S K + K − are consistent with the world average results [13,46].
The significance of the measured signal yield for the decay B 0 s → K 0 S K + K − , including systematic uncertainties, is 2.5 standard deviations. A 90% confidence level (C.L.) interval for the corresponding branching fraction relative to that of B 0 → K 0 S π + π − is derived, following the approach of Feldman-Cousins [58] ∈ [0.008 − 0.051] at 90% C.L.
The first Dalitz-plot analyses by the LHCb experiment of the dominant decays are the next step of the physics programme introduced in this work. These studies will follow and benefit from the selection methods developed for this analysis.

Appendices A Fit results by category
Signal yields and efficiencies for the different decays, data-taking periods and K 0 S reconstruction categories are shown for each of the two BDT optimisation points in Tables 3-5. Fit results for the different data-taking periods and K 0 S reconstruction categories are shown for each of the two BDT optimisation points in Figs. 3-14. Table 3: Signal yields obtained for the 2011 category from the simultaneous fit to the data. The yields shown are those obtained when fitting the data sample selected using the BDT optimisation chosen for the given decay mode. The uncertainties are statistical only. The average selection efficiencies, described in Sec. 5, are also shown for each decay mode together with the corresponding total uncertainty due to the limited simulation sample size and systematic effects in their determination.

downstream long Decay
Yield 0.0244 ± 0.0052 4 ± 3 0.0129 ± 0.0029 Table 4: Signal yields obtained for the 2012a category from the simultaneous fit to the data. The yields shown are those obtained when fitting the data sample selected using the BDT optimisation chosen for the given decay mode. The uncertainties are statistical only. The average selection efficiencies, described in Sec. 5, are also shown for each decay mode together with the corresponding total uncertainty due to the limited simulation sample size and systematic effects in their determination.

downstream long Decay
Yield S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.  Figure 4: Results of the simultaneous fit to data (downstream, 2012a) with the BDT optimisation corresponding to the favoured decay modes. The modes K 0 S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.  Figure 5: Results of the simultaneous fit to data (downstream, 2012b) with the BDT optimisation corresponding to the favoured decay modes. The modes K 0 S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1. LHCb Long Figure 6: Results of the simultaneous fit to data (long, 2011) with the BDT optimisation corresponding to the favoured decay modes. The modes K 0 S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.  Figure 7: Results of the simultaneous fit to data (long, 2012a) with the BDT optimisation corresponding to the favoured decay modes. The modes K 0 S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.
S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.
S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.
S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.  Figure 11: Results of the simultaneous fit to data (downstream, 2012b) with the BDT optimisation corresponding to the suppressed modes. The modes K 0 S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.  Figure 12: Results of the simultaneous fit to data (long, 2011) with the BDT optimisation corresponding to the suppressed modes. The modes K 0 S K + K − , K 0 S K + π − , K 0 S π + K − and K 0 S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.  Figure 13: Results of the simultaneous fit to data (long, 2012a) with the BDT optimisation corresponding to the suppressed modes. The modes K 0 S K + K − , K 0 S K + π − , K 0 S π + K − and K 0 S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.  Figure 14: Results of the simultaneous fit to data (long, 2012b) with the BDT optimisation corresponding to the suppressed modes. The modes K 0 S K + K − , K 0 S K + π − , K 0 S π + K − and K 0 S π + π − are shown from top to bottom. The left-hand side plots show the results with a logarithmic scale and the right-hand side with a linear scale. Legend is similar to that of plots shown in Fig. 1.

B Breakdown of systematic uncertainties
The full breakdown of the systematic uncertainties for each data-taking period is given in Tables 6 and 7 for, respectively, the downstream and long K 0 S reconstruction categories. Table 6: Systematic uncertainties on the ratios of branching fractions for downstream K 0 S reconstruction. All uncertainties are relative and are quoted as percentages.

C Dalitz-plot distributions of signal events
Dalitz-plot distributions of signal events as extracted using the sPlot technique are shown for each signal mode are shown in Fig. 15.  Figure 15: Distribution of sPlot weights in data. The K 0 S K + K − , K 0 S K + π − , K 0 S π + K − , and K 0 S π + π − final states are shown from top to bottom, with a B 0 parent on the left, and a B 0 s on the right. All data-taking periods and K 0 S reconstruction categories are added.