Measurement of the relative branching fractions of $B^+ \to h^+h^{\prime +}h^{\prime -}$ decays

The relative branching fractions of $B^+ \to h^+h^{\prime +}h^{\prime -}$ decays, where $h^{(\prime)}$ is a pion or kaon, are measured. The analysis is performed with a data sample, collected with the LHCb detector, corresponding to an integrated luminosity of $3.0 {\rm fb}^{-1}$ of $pp$ collisions. The results obtained improve significantly on previous measurements of these quantities, and are important for the interpretation of Dalitz plot analyses of three-body charmless hadronic decays of $B^+$ mesons.


Introduction
Three-body hadronic B meson decays to final states without any charm or charmonium hadrons are of great interest since they can be mediated by both tree and loop (socalled penguin) diagrams, and consequently CP -violation effects can manifest. Such charmless three-body decays can proceed through a number of different intermediate resonances, which increases the range of ways in which CP -violation effects can occur. Model-independent studies of the B + → K + K + K − , π + K + K − , K + π + π − and π + π + π − decays, collectively referred to as B + → h + h + h − decays, have revealed large CP -violation effects in certain regions of their Dalitz plots [1][2][3], with these results confirmed for B + → π + K + K − and π + π + π − decays by model-dependent Dalitz-plot analyses [4][5][6]. 1 It is as yet unclear whether the observed effects can be explained within the Standard Model or if new dynamics are involved.
The results of Dalitz-plot analyses typically include fit fractions of contributing resonances. These can be converted to quasi-two-body branching fractions, which can be predicted theoretically (see, for example, Refs. [7][8][9][10][11][12][13][14]), by multiplication with the branching fraction for the three-body decay. Interpretation of the data requires both branching fractions and CP asymmetries to be considered. Consequently, precise measurements of the branching fractions of charmless hadronic three-body B + decays are needed. Current knowledge of the B + → h + h + h − branching fractions, as tabulated by the Particle Data Group (PDG) [15], is summarised in Table 1. The precision ranges from 4 % to 9 %, which is not sufficient given the sensitivity of the most recent Dalitz-plot analyses. Improved knowledge of these quantities is therefore required.
The relative size of the branching fractions of B + → h + h + h − decays, as given in Table 1, can be understood to first approximation through the Cabibbo-Kobayashi-Maskawa matrix elements that enter the relevant Feynman diagrams. Examples of such diagrams are shown in Fig. 1. Interference between different amplitudes contributing to the same process can cause CP violation.
In this paper, the relative branching fractions of the B + → h + h + h − decays are determined. The analysis is based on a data sample corresponding to an integrated luminosity of 3.0 fb −1 of pp collision data collected with the LHCb detector, of which 1.0 fb −1 was collected in 2011 when the centre-of-mass energy, √ s, was equal to 7 TeV and the remaining 2.0 fb −1 was collected in 2012 at √ s = 8 TeV. Since currently B(B + → K + K + K − ) is known most precisely, results are presented primarily as ratios with this mode as the denominator. However, determinations of all ratios of one mode to another are also presented, as are the correlations between the results, in order to profit from future Decay PDG average (10 −6 ) References B + → K + K + K − 34.0 ± 1.4 [16,17] B + → π + K + K − 5.2 ± 0.4 [18,19] B + → K + π + π − 51.0 ± 2.9 [20,21] B + → π + π + π − 15.2 ± 1.4 [22]  Where final-state particles other than π + and K + are given, it should be understood that a range of resonances are possible, and where these are unflavoured in many cases decays to both π + π − and K + K − are possible. Other types of Feynman diagrams that can also contribute, such as internal W emission and annihilation amplitudes as well as rescattering processes, are not shown.
improvements of any of the individual branching fraction measurements. The analysis presented here does not include study of the suppressed three-body charmless hadronic decays B + → K + K + π − and B + → π + π + K − , which require dedicated measurements [23][24][25]. Previous measurements have used slightly different definitions of the three-body branching fractions, B(B + → h + h + h − ), and given the desired precision it is important to have a clear definition. In the work presented here, any B + → h + h + h − decay where the three final-state particles originate from the same vertex is considered to be part of the signal. This definition thus includes all charmonium resonances, since all have negligible lifetimes, and excludes all contributions from weakly decaying charm mesons. This choice differs from that used in some Dalitz-plot analyses, where contributions from the J/ψ resonance are often vetoed to avoid the need to account for resolution effects, which are negligible for other, broader, resonances. Existing knowledge of B(B + → J/ψh + ) and B(J/ψ → h + h − ) [15] is sufficient to correct for such differences in definition, which have an impact not larger than 1%.
To determine the relative branching fraction of two modes, it is necessary to know the relative signal yields and efficiency of each. By considering only ratios of these quantities, many sources of potentially large systematic uncertainty are rendered negligible. However, the efficiency of each mode depends on its Dalitz-plot distribution, and for B + → K + K + K − and K + π + π − decays the most recent Dalitz-plot models [16,17,20,21] were obtained from analyses of significantly smaller samples than those in the current analysis. To avoid a dominant systematic uncertainty due to lack of knowledge of the Dalitz-plot distributions, a model-independent approach is pursued whereby an efficiency correction is applied to each candidate depending on its Dalitz-plot position. The remainder of the paper is organised as follows. In Sec. 2 the detector and simulation software is described. The selection of signal candidates is discussed in Sec. 3, with the efficiency of these requirements, including the variations of the efficiency across the Dalitz plot of each of the final states, presented in Sec. 4. In Sec. 5 the simultaneous fit to the invariant-mass distributions of selected candidates is described, with emphasis on the various constraints that are imposed. A detailed discussion of the evaluation of systematic uncertainties is presented in Sec. 6, with the results and their correlations given in Sec. 7. A summary concludes the paper in Sec. 8.

Detector and simulation
The LHCb detector [26,27] 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 siliconstrip vertex detector surrounding the pp interaction region [28], 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 [29] placed downstream of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with 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 pp collision 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 [30]. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [31].
The online event selection is performed by a trigger [32], 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 2011 (2012) data. 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 deposited in the calorimeters. For hadrons, the transverse energy threshold is 3.5 GeV. The software trigger requires a two-, three-or four-track vertex with significant displacement from any PV. At least one charged particle must have transverse momentum p T > 1.6 GeV/c and be inconsistent with originating from a PV. A multivariate algorithm [33] is used for the identification of displaced vertices consistent with the decay of a b hadron.
In the offline selection, trigger signals are associated with reconstructed particles.
Selection requirements can therefore be made on the trigger output and on whether the decision was due to the signal candidate, other particles produced in the pp collision, or a combination of both. In this analysis it is required that the hardware trigger decision is due to either clusters in the hadronic calorimeter created by one or more of the final-state particles, or only by particles produced in the pp bunch crossing not involved in forming the B candidate.
Simulation is used to model the effects of the detector acceptance and the selection requirements. In the simulation, pp collisions are generated using Pythia [34] with a specific LHCb configuration [35]. Decays of unstable 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].

Selection of signal candidates
The procedure to select signal candidates is similar to those used in previous LHCb analyses of B + → h + h + h − decays [1][2][3][4][5][6], but is optimised for the set of relative branching fraction measurements of this analysis. A loose set of initial requirements is applied, and particle identification (PID) requirements are imposed to reject background with misidentified final-state particles. A multivariate algorithm (MVA) is used to distinguish signal from combinatorial background. Further specific requirements are applied to remove potentially large background sources from candidates where two of the final-state particles originate from a charm-or beauty-meson decay.
The initial selection includes requirements on the quality of each of the three tracks comprising the signal candidate. They are required to be displaced from all PVs, as quantified through the variable χ 2 IP , which is the difference in the vertex-fit χ 2 of a given PV reconstructed with and without the particle under consideration. The three tracks must form a common, good-quality vertex, and have invariant mass within a broad window of the known B + mass [15]. The B candidate is associated to the PV with which it forms the minimum χ 2 IP value, which must be below a certain threshold, and the B-candidate momentum must be aligned with the vector between its production and decay vertices. The B decay vertex must be displaced significantly from its associated PV. Requirements are also imposed on the p and p T of the B candidate and of the individual tracks. Variables used subsequently in the analysis are obtained from a kinematic fit to the decay [40] in which the tracks are constrained to a common vertex. For the computation of Dalitz-plot variables, the B candidate is additionally constrained to have the known B + mass [15].
Information from the ring-imaging Cherenkov detectors is combined with tracking information to obtain variables that quantify how likely a given track is to be caused by either a pion or a kaon [30]. Disjoint regions in the plane formed by these two variables are used to separate tracks that are likely to originate from kaons and unlikely to come from pions and vice versa. For each of the four final states, requirements on these PID variables are imposed to reduce the potential cross-feed background from misidentification of the other modes. Optimal requirements are evaluated by considering the figure of merit N S / √ N S + N B cf , where N S and N B cf are the expected signal and cross-feed background yields for each case. The relative sizes of N S and N B cf depend on the branching fractions of the four signal modes, which are taken from previous measurements [15], as well as Table 2: Probability (%), due to the particle identification requirements, for each of the four signal modes to be correctly identified, or to form a cross-feed background to one of the other final states. Empty entries correspond to values below 0.05 %. The decays B + → π + K + K − and B + → K + π + π − can, through both K + → π + and π + → K + misidentification, appear as a cross-feed background in the correct final state with probabilities of below 0.05 % and 0.4 %, respectively.

Decay
Reconstructed final state efficiencies and misidentification rates. These are determined from data control samples of D * + → D 0 (→ K − π + )π + decays [30], weighted to reproduce the p and η distributions of signal tracks, and -since the PID performance depends on the detector occupancy -the number of reconstructed tracks in the pp bunch crossing. Requirements on the ranges of these three variables are applied to ensure reliable performance of the PID calibration procedure. Tracks are also required to not have any associated signal in the muon detectors. For the B + → π + K + K − channel, the expected significant cross-feed background from partially reconstructed B → K + π + π − X and K + K + K − X decays, where X denotes any additional particles, is accounted for by doubling the value of N B cf from that obtained considering the three-body B + decays only. A baseline set of PID requirements is applied, in the cases where the optimisation procedure returns loose values, to ensure that no candidate can be selected in more than one of the final states under consideration. The outcome of this procedure is a set of requirements that, after further tightening in certain regions of phase space as described below, corresponds to the efficiencies and misidentification rates given in Table 2.
Variables that provide good discrimination between signal and combinatorial background without introducing significant distortions into the B-candidate mass or Dalitz-plot distributions, are identified for inclusion in the MVA. In order of discriminating power, these are: the pointing angle, which characterises how well the B-candidate momentum aligns with the vector from the associated PV to the B decay vertex; the p T asymmetry, which quantifies the isolation of the B candidate through the p T asymmetry between itself and other tracks within a cone around its flight direction [41]; the distance between the B-candidate production and decay vertices, divided by its uncertainty; the χ 2 of the B-candidate vertex; the χ 2 IP of the track with the largest p T out of the three that form the B candidate; the p of the same track; the χ 2 IP of the B candidate. These variables are distributed almost identically for all signal modes, justifying the use of a single MVA. The distributions of all input variables, and the MVA output, are confirmed to agree well between data and simulation, where the data distributions are obtained from the B + → K + K + K − sample with background subtracted using weights obtained from a fit to the B-candidate mass distribution [42].
The combination of variables into the MVA is implemented with the NeuroBayes package [43]. The MVA is trained to discriminate between a signal sample, taken from simulation, and a background sample obtained from the data sideband with B-candidate mass values significantly above the B + mass. Since the decay B + → π + K + K − is the most challenging of the four modes to separate from background, the training is performed with both signal and background samples corresponding to that mode, with initial selection and PID requirements applied. A requirement on the output of the MVA is optimised by considering the figure of merit N S / √ N S + N B cb , where N B cb is the expected combinatorial background yield in the signal region [5240, 5320] MeV/c 2 .
Background from B + → D 0 h + decays, with D 0 → K + π − , K + K − or π + π − , passes the selection requirements for the correctly reconstructed final state in large numbers, since the D 0 lifetime is sufficiently small that the three tracks can still form a good B-candidate vertex. This background is vetoed by removing any candidate with one of the corresponding two-body invariant masses in the region [1830, 1890] MeV/c 2 . Such decays can still cause background when final-state particles are misidentified. Tighter PID requirements are imposed when one of the two-body invariant masses of oppositely charged final-state particles is in the range [1890, 2000] MeV/c 2 for π → K misidentification or [1700, 1850] MeV/c 2 for K → π misidentification. These requirements reduce most misidentified charm background components to negligible levels with minimal impact on the signal efficiency.
The so-called partially combinatorial background, where a two-body B-meson decay is combined with a random track, can populate the B-candidate invariant-mass region at values above the signal peaks. The shape of such background can be hard to model in the B-candidate invariant-mass fit, introducing a potential source of systematic uncertainty on the signal yield. Therefore, candidates that may contain  [44] is negligible. After all selection requirements are imposed, a small fraction of selected pp bunch crossings, ranging from 0.2 % for the K + K + K − final state to 2.4 % for π + π + π − , contain more than one B candidate. In such cases, only the candidate with the highest MVA output value is retained. The systematic uncertainty associated with this procedure is negligible.

Signal efficiency
The total signal efficiency, tot , can be expressed in terms of factorising components, where sel+geom includes the effects of the geometrical efficiency of the LHCb detector and of both online and offline selection requirements, and PID is the PID efficiency for candidates that have passed the selection requirements. The former can be evaluated quite reliably from simulation, although small data-driven corrections are applied, while the latter is obtained from control samples. As explained in Sec. 1, the variation of the efficiency across the phase space, or Dalitz plot, of each decay, must be accounted for. It is convenient to do so using the so-called square Dalitz plot (SDP) representation of the phase space, since this provides greater granularity in regions close to the edges of the regular Dalitz plot where resonances tend to populate and where the efficiency variation tends to be larger. Moreover, the SDP definition in terms of two variables m and θ , each of which is bounded in the range [0, 1], aligns a rectangular grid with the edges of the phase space, avoiding edge effects associated with rectangular binning of the regular Dalitz plot. The variable m is a transformation of the invariant mass of two of the three final-state particles, while θ is a transformation of the helicity angle associated with that pair, i.e. the angle between the momentum of one of the pair and the third particle in the rest frame of the pair. The explicit definitions are [45] where the ordering of the particles used in the analysis is given in Table 3, m α is the mass of the particle labelled α and m αβ is the two-body invariant mass of particles α and β.
For decays with two identical particles, i.e. B + → K + K + K − and B + → π + π + π − , the SDP is folded along the line θ = 0.5, making the initial ordering, i.e. which of the two identical particles is i and which is j, irrelevant. The simulated samples of signal decays used in the analysis to determine sel+geom are generated with uniform density in these SDP coordinates. The impact of the hardware trigger is a potentially significant source of discrepancy between data and simulation in the evaluation of sel+geom . Corrections to the simulation are applied for two mutually exclusive subsamples of the selected candidates. The first includes candidates that are triggered at hardware level by clusters in the hadronic calorimeter created by one or more of the final-state signal particles, and the second contains those triggered only by other particles produced in the pp bunch crossing. For the first subsample, a correction is calculated from the probability of an energy deposit in the hadronic calorimeter to fire the trigger, evaluated from calibration data samples as a function of particle type (kaon or pion), charge, dipole magnet polarity, transverse energy and position in the calorimeter. In the second subsample, the simulation is weighted so that the rates of the different categories of hardware trigger (hadron, muon, dimuon, electron, photon) match those observed in data. As described in Sec. 6, the former of these corrections has a non-negligible impact on the results, while the effect of the latter is smaller. Additional small corrections are applied to the simulation to ensure that the tracking efficiency [46], and the kinematic (p T , η) distributions of selected B mesons match those of data.
The PID efficiency is calculated, in the same way as described above for the optimisation of the PID requirements, from calibration samples. The efficiencies for each final-state particle are parameterised in terms of their total and transverse momentum, and the number of tracks in the event, and these are multiplied to form the overall efficiency PID .
The total efficiency, tot , is shown in Fig. 2 as a function of SDP position for the four signal modes, with all selection requirements except the charm vetoes applied. Bands in the phase space are nevertheless visible around the charm-meson mass due to the tighter PID requirements applied in these regions. For example, the depleted region in tot for B + → K + K + K − decays is due to tightened PID requirements to remove B + → D 0 (→ K + π − )K + decays with π − → K − misidentification. The choice of 30 × 30 bins in these efficiency maps is made so that the minimum bin content remains above 10 and hence the efficiency in each bin is determined with reasonably small uncertainty, although some fluctuations are visible at the edges, and particularly the corners, of the SDP. These fluctuations occur where the Jacobian of the transformation from conventional to SDP coordinates takes extreme values, and hence affect modes with final-state pions more than kaons.
Since candidate-by-candidate efficiency corrections are applied in the evaluation of the relative branching fractions, the impact of charm vetoes that completely remove regions of phase space are accounted for separately. The veto efficiencies are determined by generating ensembles of samples according to the most recent Dalitz-plot models of the Decay veto (%) 98.05 ± 0.10 signal modes [4-6, 17, 21], and evaluating the impact of the veto. Each sample contains a number of decays sampled from a Poisson distribution with mean corresponding to the signal yield in the analysis where the model was determined, and the corresponding uncertainties are estimated from the spread of veto efficiency values in the ensemble. The efficiencies obtained for each channel, veto , are given in Table 4.

B-candidate invariant-mass fit
A simultaneous unbinned extended maximum-likelihood fit is performed to the four Bcandidate invariant-mass distributions, in the range [5100, 5500] MeV/c 2 , to determine the yields of the signal components. The fit model includes components for signal, crossfeed from misidentified three-body B decays, partially reconstructed background and combinatorial background. The signal mass distributions are modelled as the sum of two Crystal Ball functions [47], with a common peak position and width, and tails to opposite sides of the peak. The shape parameters of the double Crystal Ball function are determined from fits to simulation and then fixed in the data fit, with the exception of an offset to the peak position and a scaling factor of the width. These two parameters, shared by all four modes, are both left free to vary in the fit to data to account for small differences between data and simulation.
All possible cross-feed background contributions from one B + → h + h + h − decay to another, or to itself, with single or double misidentification are accounted for in the fit. The shapes are described empirically with the sum of two Crystal Ball functions, with parameters obtained from simulated samples weighted to reproduce the underlying Dalitz-plot distributions [4-6, 17, 21] and with per-track data-calibrated PID efficiencies applied. The peak positions and widths of these shapes are adjusted, in the fit to data, by the same offset and scale factor as the signal functions. Other potential sources of similar background, involving misidentified three-body b-hadron decays such as Ξ + b → h + h + p [48] are found to have negligible contribution.
The sources of partially reconstructed background differ between the four final states considered. All include a component from four-body charmless B + and B 0 decays with an additional soft neutral or charged pion that is not reconstructed. The shapes of these, and all partially reconstructed background components, are modelled with ARGUS functions [49], where the threshold is fixed to the known difference between the B-meson and pion masses [15], convolved with a Gaussian resolution function with width of the corresponding signal mode. The shape parameters are fixed to the values obtained from fitting simulated samples of the background.
For all modes except B + → K + K + K − , there is significant background from B 0 s → D − s π + decays, with subsequent D − s decay to the corresponding pair of particles plus an additional soft pion that is not reconstructed. The shapes of these components differ from those of the corresponding charmless four-body decays because of differences in the momentum distributions of the missing pion. The same parametric functions are used as for the charmless four-body decays, but with parameters determined independently from appropriate simulation samples.
The π + K + K − final state has a further source of partially reconstructed background through B 0 s → π + K + K − π − decays. The latest study of this process [50] reveals that it is composed of a mixture of Kπ resonances, rather than being dominated by the B 0 s → K * (892) 0 K * (892) 0 decay, so a data-driven approach is used to determine the shape of this component.
The K + π + π − final state contains background from B + → η K + with η → π + π − γ decays. The ARGUS function shape parameter is determined from a fit to a sample of simulation weighted to reproduce the appropriate π + π − invariant-mass shape [51]. The threshold parameter is fixed to the peak value of the B + → K + π + π − signal decay including, in the fit to data, the offset.
Background to the B + → π + π + π − decay from misidentified B + → D 0 (→ K + π − )π + decays remains at non-negligible level after the PID requirements. This is modelled in the fit with an ARGUS function convolved with a Gaussian resolution, with parameters determined from a fit to simulation, in a similar way as for the partially reconstructed background. Misidentified B + → D 0 (→ K + π − )π + decays are also a source of background in the π + K + K − final state, but this is found to be readily absorbed by other fit components and is therefore not included explicitly. The combinatorial background in each final state is described by an exponential function.
The free parameters of the fit are the four signal yields, the common offset and scale factor of the signal shape functions, the four combinatorial background yields and their associated exponential shape parameters, one partially reconstructed background yield for each of the K + K + K − , π + K + K − and π + π + π − final states and two for the K + π + π − channel. All misidentified background yields are constrained, within uncertainty, to their expected levels based on the signal yields in the corresponding correctly identified final states and the known misidentification probabilities, as given in the off-diagonal elements of Table 2. For background from misidentified B + → D 0 (→ K + π − )π + decays, the known branching fraction, relative to those of the signal channels, also enters the calculation of the constraint. Similarly, the relative yields of the different sources of partially reconstructed background in the π + K + K − and π + π + π − final states, and of the B + → η K + background to the K + π + π − final state, are constrained to their expected values.
The invariant-mass distributions m(h + h + h − ) for selected candidates in all four signal modes together with the fit projections are shown in Fig. 3 for the K + K + K − and π + K + K − final states and in Fig. 4 for the K + π + π − and π + π + π − final states. The signal yields are given in Table 5. There is good agreement of the fit model with the data in all four final states, with some potential small residual discrepancies accounted for as sources of systematic uncertainty. The stability of the fit is investigated with pseudoexperiments, and the signal yields are found to be unbiased within the statistical precision of the ensemble.

Decay
Fit yield 6 Systematic uncertainties Systematic uncertainties are minimised by measuring the ratios of the B + → h + h + h − branching fractions relative to one another, but given the statistical precision of the results several sources of significant uncertainty remain. These originate from possible imperfections in the fit model used to determine the signal yields and the precision with which the relative efficiencies are known. A summary of the uncertainties assigned to each ratio of branching fractions is given in Table 6. Pseudoexperiments are used to determine the effect on the signal yields of using alternative shapes to describe the different fit components. Three variants of the fit model are constructed where in each an alternative shape is used for a particular category of fit component. In Model I, the signal and cross-feed components are changed to double Hypatia functions [52]. In Model II, a set of Chebyshev polynomials up to second order is used to describe the combinatorial background shape. In Model III, the partially reconstructed background shapes are replaced with non-parametric functions. The pseudoexperiments are generated according to the alternative model, then fitted with both the baseline and alternative model. The mean of the distribution of the difference between the results with the two models is taken as the corresponding systematic uncertainty. Overall, the Model II and III uncertainties are the dominant sources of systematic uncertainty for all measured branching fraction ratios. Uncertainty from possible bias on the fitted yields is also investigated using pseudoexperiments, generated and fitted using the nominal fit model. The effect of the fixed parameters in the fit model is estimated by evaluating the impact of varying these parameters within their uncertainties.
Uncertainties on the signal efficiencies originate from residual differences in the be-haviour of data and simulation, as well as the limited size of the simulation and control samples. Data-driven corrections are applied in the determination of the signal efficiency related to the performance of the hardware trigger (denoted L0 TOS and L0 TIS in Table 6 for cases where the trigger is associated to the tracks that comprise the B candidate and to other particles in the event, respectively), the reconstruction of tracks, and the B-meson production kinematics. The L0 TOS uncertainty is determined from the difference between results with and without the correction applied; this is a more conservative approach compared to those used for other uncertainties, reflecting the fact that the method used to obtain the correction does not account for all possible variables that the efficiency may depend upon. Effects associated with the reweighting of L0 TIS categories, and with the correction to the track reconstruction efficiency, are both determined by varying the correction within its uncertainties. The systematic uncertainty associated with the production kinematics correction is estimated by determining the correction factors from an alternative background-subtracted data sample. Possible small differences between data and simulation in the distribution of the variables included in the MVA are accounted for by weighting the simulated events to match the distributions observed in data. The changes in results when this weighting is applied are assigned as the associated systematic uncertainties. Uncertainty in the efficiency of the charm vetoes is obtained by propagating the corresponding values, given in Table 4. Effects related to the choice of binning of the efficiency maps are estimated by changing the granularity, while those due to the finite size of the simulated signal samples (denoted "MC stats" in Table 6) are evaluated by varying the efficiency maps according to the uncertainties in each SDP bin. The determination of the PID efficiency from control samples is also a source of uncertainty. Effects related to the differing kinematic distributions of tracks in the signal modes and the control samples, to the finite size of the control samples, and to the background-subtraction procedure are determined.
The stability of the results is cross-checked by determining the relative branching fraction ratios in various subsets of the data. The data are subdivided by year of datataking and (separately) by magnet polarity, with consistent results obtained. When comparing results obtained in subsamples separated by hardware trigger decision, by B-meson pseudorapidity and by detector occupancy some discrepancies can be seen if considering statistical uncertainties alone. These, however, are compatible with the size of relevant systematic uncertainties.

Results
The relative branching fractions of the signal modes are determined, for example with B + → K + K + K − as denominator, as where N corr is, for the mode indicated in the subscript, the efficiency-corrected signal yield accounting both for the variation of the total efficiency across the SDP and for the charm vetoes that completely remove certain regions of the phase space. These   B ratio Value 3.488 ± 0.035 ± 0.053 efficiency-corrected yields are [42] N corr = 1 veto N bins where the index j runs over the N bins bins of the SDP, tot j is the corresponding efficiency in bin j (as given in Fig. 2), and for each value of j the index i runs over the candidates in that bin. The per-candidate signal sWeights w i , which implement the background subtraction, are obtained from individual fits to the B-candidate mass distribution of each mode in which all nuisance parameters are fixed to the values obtained in the simultaneous fit. In these fits the only varying parameters are the yields of the signal and all background components except those of the cross-feed background contributions, which are fixed. The term cM j accounts for these fixed components, where the coefficient c is determined from the fit [42] and M j is the fraction of the cross-feed background in SDP bin j. The statistical uncertainty on each N corr value is calculated as described in Ref. [53], accounting for the reduction in the uncertainties of the yields, compared to the baseline fit, due to the nuisance parameters being fixed.
The complete set of results for twelve relative branching fractions of B + → h + h + h − decays is shown in Table 7. Six of these are the inverse of the other six. Moreover, since there are only three independent measurements, correlations between the ratios must Table 8: Statistical correlations between the measured branching fraction ratios.
−0.01 −0.34 0.78 −0.08 0.32 also be taken into account. The statistical and systematic correlations are presented in Tables 8 and 9, respectively. The statistical correlations are determined from ensembles of pseudoexperiments. In each experiment, the ratios are calculated and the correlation is obtained from the distribution of one ratio against another in the ensemble. Large statistical correlations are observed between the two ratios that share a decay with a yield that is small compared to that of the other decay channel in the ratios; this affects in particular pairs of ratios that have B + → π + K + K − as a common channel. Ratios which do not have any mode in common have smaller correlations, which can however be non-zero due to the nature of the simultaneous fit from which the yields are obtained. Correlations related to systematic uncertainties obtained from ensembles of pseudoexperiments, as described in Sec. 6 are evaluated with the same method as the statistical correlations. For those that are determined from the difference between the results obtained when a single variation is made and those in the baseline analysis, 100% correlation or anticorrelation (depending on the relative sign of the shift) is assumed. For each source of systematic uncertainty, these correlations are converted into a covariance matrix. These are summed, and the total systematic covariance matrix thus obtained is converted back into the total systematic correlation matrix. The size of the systematic correlations is related to whether two ratios share dominant sources of systematic uncertainty. In particular, for pairs of ratios with B + → π + K + K − as a common channel, the uncertainty due to limited knowledge of the background shapes induces significant correlations.

Summary
Data collected by the LHCb experiment in 2011 and 2012, corresponding to an integrated luminosity of 3.0 fb −1 , has been used to determine the relative branching fractions of the B + → h + h + h − decays. The measured ratios relative to the B + → K + K + K − channel are B(B + → π + K + K − )/B(B + → K + K + K − ) = 0.151 ± 0.004 (stat) ± 0.008 (syst) , B(B + → K + π + π − )/B(B + → K + K + K − ) = 1.703 ± 0.011 (stat) ± 0.022 (syst) , B(B + → π + π + π − )/B(B + → K + K + K − ) = 0.488 ± 0.005 (stat) ± 0.009 (syst) . Table 9: Systematic correlations between the measured branching fraction ratios. π + K + K − K + K + K − K + π + π − K + K + K − π + π + π − K + K + K − K + π + π − π + K + K − π + π + π − π + K + K − π + π + π − K + π + π − B(B + →π +  The dominant systematic uncertainties are related to knowledge of the background shapes in the invariant-mass fit, and are reducible if knowledge of the various sources of background can be improved or if the background can be suppressed in future analyses. Several other sources of systematic uncertainty are, however, not negligible compared to the statistical uncertainty of these results, so that further significant reduction in uncertainty will be challenging. Comparisons with the current world averages are given, for the three measurements above, in Fig. 5. All measurements are in good agreement with the previous world-average results and, furthermore, significant improvement in the precision of all measured ratios is obtained.  Figure 5: Comparisons of the measured branching fraction ratios, with B(B + → K + K + K − ) as denominator, with the current world averages [15]. Light (dark) bands associated with the branching fraction ratio correspond to the ±1σ total (statistical) uncertainty intervals. For horizontal and vertical bands taken from the PDG only the total uncertainty is shown.