First observation of the decays $\bar{B}^0_{(s)}\to D_s^+K^-\pi^+\pi^-$ and $\bar{B}^0_s\to D_{s1}(2536)^+\pi^-$

The first observation of the decays $\bar{B}^0_{s}\to D_s^+K^-\pi^+\pi^-$ and $\bar{B}^0\to D_s^+K^-\pi^+\pi^-$ are reported using an integrated luminosity of 1.0 fb$^{-1}$ recorded by the LHCb experiment. The branching fractions, normalized with respect to $\bar{B}^0_{s}\to D_s^+\pi^-\pi^+\pi^-$ and $\bar{B}^0_{s}\to D_s^+K^-\pi^+\pi^-$, respectively, are measured to be {\br(\bar{B}^0_{s}\to D_s^+K^-\pi^+\pi^-)\over\br(\bar{B}^0_{s}\to D_s^+\pi^-\pi^+\pi^-)}&= (5.2\pm0.5\pm0.3)\times10^{-2}, {\br(\bar{B}^0\to D_s^+K^-\pi^+\pi^-)\over\br(\bar{B}^0_{s}\to D_s^+K^-^\pi^+\pi^-)}&= 0.54\pm0.07\pm0.07, where the first uncertainty is statistical and the second is systematic. The $\bar{B}^0_{s}\to D_s^+K^-\pi^+\pi^-$ decay is of particular interest as it can be used to measure the weak phase $\gamma$. First observation of the $\bar{B}^0_s\to D_{s1}(2536)^+\pi^-, D_{s1}^+\to D_s^+\pi^-\pi^+$ decay is also presented, and its branching fraction relative to $\bar{B}^0_{s}\to D_s^+\pi^-\pi^+\pi^-$ is found to be {\br(\bar{B}^0_s\to D_{s1}(2536)^+\pi^-, D_{s1}^+\to D_s^+\pi^-\pi^+)\over\br(\bar{B}^0_{s}\to D_s^+\pi^-\pi^+\pi^-)}&= (4.0\pm1.0\pm0.4)\times10^{-3}.

First observation of the decays B 0 (s) → D + s K − π + π − and B 0 s → D s1 (2536) + π − The LHCb collaboration † Abstract The first observation of the decays B 0 s → D + s K − π + π − and B 0 → D + s K − π + π − are reported using an integrated luminosity of 1.0 fb −1 recorded by the LHCb experiment. The branching fractions, normalized with respect to B 0 s → D + s π − π + π − and B 0 s → D + s K − π + π − , respectively, are measured to be B(B 0 s → D + s K − π + π − ) B(B 0 s → D + s π − π + π − ) = (5.2 ± 0.5 ± 0.3) × 10 −2 , where the first uncertainty is statistical and the second is systematic. The B 0 s → D + s K − π + π − decay is of particular interest as it can be used to measure the weak phase γ. First observation of the B 0 s → D s1 (2536) + π − , D + s1 → D + s π − π + decay is also presented, and its branching fraction relative to B 0 s → D + s π − π + π − is found to be

Introduction
In the Standard Model (SM), the amplitudes associated with flavor-changing processes depend on four Cabibbo-Kobayashi-Maskawa (CKM) [1, 2] matrix parameters. Contributions from physics beyond the Standard Model (BSM) add coherently to these amplitudes, leading to potential deviations in rates and CP-violating asymmetries when compared to the SM contributions alone. Since the SM does not predict the CKM parameters, it is important to make precise measurements of their values in processes that are expected to be insensitive to BSM contributions. Their values then provide a benchmark to which BSM-sensitive measurements can be compared. The least well-determined of the CKM parameters is the weak phase γ ≡ arg − V * ub V ud V * cb V cd , which, through direct measurements, is known to a precision of ∼ 10 o − 12 o [3,4]. It may be probed using time-independent rates of decays such as B − → DK − [5,6,7], or by analyzing the time-dependent decay rates of processes such as B 0 s → D ∓ s K ± [8,9,10,11]. Sensitivity to the weak phase γ results from the interference between b → c and b → u transitions, as indicated in Figs. 1(a-c). Such measurements may be extended to multibody decay modes, such as B − → DK − π + π − [12] for a timeindependent measurement, or B 0 s → D + s K − π + π − in the case of a time-dependent analysis. The B 0 → D + K − π + π − decay, while having the same final state as B 0 s → D + s K − π + π − , receives contributions not only from the W -exchange process ( Fig. 1(d)), but also from b → c transitions in association with the production of an extra ss pair (Figs. 1(e-f)). The decay may also proceed through mixing followed by a b → u, W -exchange process (not shown). However, this amplitude is Cabibbo-, helicity-and color-suppressed, and is therefore negligible compared to the b → c amplitude.
This paper reports the first observation of B 0 s → D + s K − π + π − and B 0 → D + s K − π + π − , and measurements of their branching fractions relative to B 0 s → D + s π − π + π − and B 0 s → D + s K − π + π − , respectively. The data sample is based on an integrated luminosity of 1.0 fb −1 of pp collisions at √ s = 7 TeV collected by the LHCb experiment. The same data sample is also used to observe the B 0 s → D s1 (2536) + π − , D + s1 → D + s π − π + decay for the first time, and measure its branching fraction relative to B 0 s → D + s π − π + π − . The inclusion of charge-conjugated modes is implied throughout this paper.

Detector and simulation
The LHCb detector [13] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift-tubes placed downstream. The combined tracking system has a momentum resolution (∆p/p) that varies from 0.4% at 5 GeV/c to 0.6% at 100 GeV/c, and an impact parameter (IP) resolution of 20 µm for tracks with high transverse mo-  Figure 1: Diagrams contributing to the B 0 s , B 0 s → D + s K − π + π − (a-c) and B 0 → D + s K − π + π − (d-f) decays, as described in the text. In (a-d), the additional (π + π − ) indicates that the K − π + π − may be produced either through an excited strange kaon resonance decay, or through fragmentation. mentum (p T ). Charged hadrons are identified using two ring-imaging Cherenkov (RICH) detectors. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and pre-shower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
The trigger consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. The software trigger requires a two-, three-or four-track secondary vertex with a high p T sum of the tracks and a significant displacement from the primary pp interaction vertices (PVs). At least one track should have p T > 1.7 GeV/c, an IP χ 2 greater than 16 with respect to all PVs, and a track fit χ 2 /ndf < 2, where ndf is the number of degrees of freedom. The IP χ 2 is defined as the difference between the χ 2 of the PV reconstructed with and without the considered particle. A multivariate algorithm is used for the identification of secondary vertices [14].
For the simulation, pp collisions are generated using Pythia 6.4 [15] with a specific LHCb configuration [16]. Decays of hadronic particles are described by EvtGen [17] in which final state radiation is generated using Photos [18]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [19] as described in Ref. [20].

Signal selection
Signal B 0 (s) decay candidates are formed by pairing a D + s → K + K − π + candidate with either a π − π + π − (hereafter referred to as X d ) or a K − π + π − combination (hereafter referred to as X s ). Tracks used to form the D + s and X d,s are required to be identified as either a pion or a kaon using information from the RICH detectors, have p T in excess of 100 MeV/c, and be significantly detached from any reconstructed PV in the event.
Signal D + s candidates are required to have good vertex fit quality, be significantly displaced from the nearest PV, and have invariant mass, M(K + K − π + ), within 20 MeV/c 2 of the D + s mass [21]. To suppress combinatorial and charmless backgrounds, only those D + s candidates that are consistent with decaying through either the φ (M( [21]). The remaining charmless background yields are determined using the D + s mass sidebands. For about 20% of candidates, when the K + is assumed to be a π + , the corresponding K − π + π + invariant mass is consistent with the D + mass. To suppress cross-feed from B 0 → D + X decays, a tighter particle identification (PID) requirement is applied to the K + in the D + s → K + K − π + candidates when |M(K − π + π + ) − m D + | < 20 MeV/c 2 (m D + is the D + mass [21]). Similarly, if the invariant mass of the particles forming the D + s candidate, after replacing the K + mass with the proton mass, falls within 15 MeV/c 2 of the Λ + c mass, tighter PID selection is applied. The sizes of these mass windows are about 2.5 times the invariant mass resolution, and are sufficient to render these cross-feed backgrounds negligible.
Candidate X d and X s are formed from π − π + π − or K − π + π − combinations, where all invariant mass values up to 3 GeV/c 2 are accepted. To reduce the level of combinatorial background, we demand that the X d,s vertex is displaced from the nearest PV by more than 100 µm in the direction transverse to the beam axis and that at least two of the daughter tracks have are suppressed by applying more stringent PID requirements to the K − and π + in X s . The PID requirements have an efficiency of about 65% for selecting X s , while rejecting about 97% of the favored threepion background. To suppress peaking backgrounds from B 0 is more than 20 MeV/c 2 away from the D + s mass. Signal B meson candidates are then formed by combining a D + s with either an X d or X s . The reconstructed B candidate is required to be well separated from the nearest PV with a decay time larger than 0.2 ps and have a good quality vertex fit. To suppress remaining charmless backgrounds, which appear primarily in B 0 → D + s K − π + π − , the vertex separation (VS) χ 2 between the D + s and B decay vertices is required to be greater than 9. Candidates passing all selection requirements are refit with both D + s mass and vertex constraints to improve the mass resolution [22].
To further suppress combinatorial background, a boosted decision tree (BDT) selection [23] with the AdaBoost algorithm [24] is employed. The BDT is trained using simulated B 0 s → D + s K − π + π − decays for the signal distributions, and the high B mass sideband in data are used to model the backgrounds. The following thirteen variables are used: • B candidate: IP χ 2 , VS χ 2 , vertex fit χ 2 , and p T ; • X d,s candidate: IP χ 2 , maximum of the distances of closest approach between any pair of tracks in the decay; • X d,s daughters: min(IP χ 2 ), max(IP χ 2 ), min(p T ); where min and max denote the minimum and maximum of the indicated values amongst the daughter particles. The flight distance significance is the separation between the D + s and B vertices, normalized by the uncertainty. The training produces a single variable, x, that provides discrimination between signal decays and background contributions. The cut value is chosen by optimizing S(x cut )/ S(x cut ) + B(x cut ), where S(x cut ) and B(x cut ) are the expected signal and background yields, respectively, after requiring x > x cut . At the optimal point, a signal efficiency of ∼90% is expected while rejecting about 85% of the combinatorial background (after the previously discussed selections are applied). After all selections, about 3% of events have more than one signal candidate in both data and simulation. All candidates are kept for further analysis.

Fits to data
The B 0 s → D + s π − π + π − and B 0 (s) → D + s K − π + π − invariant mass spectra are each modeled by the sum of a signal and several background components. The signal shapes are obtained from simulation, and are each described by the sum of a Crystal Ball (CB) [25] shape and a Gaussian function. The CB shape parameter that describes the tail toward low mass is fixed based on simulated decays. A common, freely varying scale factor multiplies the width parameters in the CB and Gaussian functions to account for slightly larger resolution in data than in simulation. For the B 0 (s) → D + s K − π + π − mass fit, the difference between the mean B 0 s and B 0 masses is fixed to 87.35 MeV/c 2 [21]. Several non-signal b-hadron decays produce broad peaking structures in the D + s π − π + π − and D + s K − π + π − invariant mass spectra. For B 0 s → D + s π − π + π − , the only significant source of peaking background is from B 0 s → D * + s π − π + π − , where the photon or π 0 from the D * + s decay is not included in the reconstructed decay. Since the full decay amplitude for B 0 s → D * + s π − π + π − is not known, the simulation may not adequately model the decay. Simulation is therefore used to provide an estimate for the shape, but the parameters are allowed to vary within one standard deviation about the fitted values.
mass fit. This same shape is assumed for both B 0 and B 0 s , where for the former, a shift by the B 0 − B 0 s mass difference is included. For the B 0 s → D + s π − π + π − and B 0 s → D * + s π − π + π − cross-feed, simulated decays and kaon misidentification rates taken from D * + calibration data are used to obtain their expected yields and invariant mass shapes. The cross-feed contribution is about 3% of the B 0 s → D + s π − π + π − and B 0 s → D * + s π − π + π − yields; the corresponding cross-feed yields are fixed in the B 0 (s) → D + s K − π + π − fit. The shape is obtained by parameterizing the invariant mass spectrum obtained from the simulation after replacing the appropriate π − mass in X d with the kaon mass. The combinatorial background is described by an exponential function whose slope is allowed to vary independently for both mass fits. Figure 2 shows the invariant mass distribution for B 0 s → D + s π − π + π − candidates passing all selection criteria. The fitted number of B 0 s → D + s π − π + π − signal events is 5683 ± 83. While it is expected that most of the low mass background emanates from B 0 s → D * + s π − π + π − decays, contributions from other sources such as B 0 s → D + s π − π + π − π 0 are also possibly absorbed into this background component. Figure 3 shows the invariant mass distribution for B 0 (s) → D + s K − π + π − candidates. The fitted signal yields are 402±33 B 0 → D + s K − π + π − and 216 ± 21 B 0 s → D + s K − π + π − events. The D + s mass sidebands, defined to be from 35 to 55 MeV/c 2 on either side of the nominal D + s mass, are used to estimate the residual charmless background that may contribute to the observed signals. The numbers of B 0 s decays in the D + s sidebands are 61 ± 16 , 0 +5 −0 , and 9 ± 5 for the B 0 s → D + s π − π + π − , B 0 s → D + s K − π + π − and B 0 → D + s K − π + π − decays, respectively; they are subtracted from the observed signal yields to obtain the corrected number of signal decays. The yields in the signal and sideband regions are summarized in Table 1.

Mass distributions of X d,s and two-body masses
In order to investigate the properties of these B 0 (s) decays, sWeights [26] obtained from the mass fits are used to determine the underlying X d,s invariant mass spectra as well as the two-body invariant masses amongst the three daughter particles. Figure 4 shows (a) the π − π + π − mass, (b) the smaller π + π − mass and (c) the larger π + π − mass in B 0 s → D + s π − π + π − data and simulated decays. A prominent peak, consistent with the a 1 (1260) − → π − π + π − is observed, along with structures consistent with the ρ 0 in the two-body masses. There appears to be an offset in the peak position of the a 1 (1260) −  between data and simulation. Since the mean and width of the a 1 (1260) − resonance are not well known, and their values may even be process dependent, this level of agreement is reasonable. A number of other spectra have been compared between data and simulation, such as the p T spectra of the D + s , X d and the daughter particles, and excellent agreement is found. Figure 5 shows the corresponding distributions for the B 0 s → D + s K − π + π − decay. A  peaked structure at low K − π + π − mass, consistent with contributions from the lower-lying excited strange mesons, such as the K 1 (1270) − and K 1 (1400) − , is observed. As many of these states decay through K * 0 and ρ 0 mesons, significant contributions from these resonances are observed in the K − π + and π + π − invariant mass spectra, respectively. The simulation provides a reasonable description of the distributions in the data. Figure 6 shows the same distributions for B 0 → D + s K − π + π − . The K − π + π − invariant Table 1: Summary of event yields from data in the D + s signal and sidebands regions, and the background corrected yield. The signal and sideband regions require D + s candidates to have invariant mass |M(K + K − π + ) − m D + s | < 20 MeV/c 2 and 35 < |M(K + K − π + ) − m D + s | < 55 MeV/c 2 , respectively, where m D + s is the D + s mass [21].  Figure 4: Invariant mass distributions for (a) X d , (b) smaller π + π − mass in X d and (c) the larger π + π − mass in X d , from B 0 s → D + s π − π + π − decays using sWeights. The points are the data and the solid line is the simulation. The simulated distribution is normalized to have the same yield as the data. mass is quite broad, with little indication of any narrow structures. There are indications of K * 0 and ρ 0 contributions in the K − π + and π − π + invariant mass spectra, respectively, but the contribution from resonances such as the K 1 (1270) − or K 1 (1400) − appear to be small or absent. In the K − π + invariant mass spectrum, there may be an indication of a K * 0 (1430) 0 contribution. The simulation, which models the K − π + π − final state as 10% K 1 (1270) − , 10% K 1 (1400) − , 40% K * 0 π − and 40% K − ρ 0 , provides a reasonable description of the data, which suggests that processes such as those in Figs. 1(e-f) constitute a large portion of the total width for this decay.  Figure 5: Invariant mass distributions for (a) X s , (b) π + π − in X s and (c) the K − π + in X s , from B 0 s → D + s K − π + π − data using sWeights. The points are data and the solid line is the simulation. The simulated distribution is normalized to have the same yield as the data.  Figure 6: Invariant mass distributions for (a) X s , (b) π + π − in X s and (c) the K − π + in X s , from B 0 → D + s K − π + π − data using sWeights. The points are data and the solid line is the simulation. The simulated distribution is normalized to have the same yield as the data.

Decay Signal Region Sideband Region Corrected Yield
B 0 s mass are selected, and from them the invariant mass difference, ∆M = M(D + s π − π + )− M(D + s ) is formed, where both π + π − combinations are included. The ∆M distribution for candidates in the B 0 s signal window is shown in Fig. 7. A peak corresponding to the D s1 (2536) + is observed, whereas no significant structures are observed in the upper B 0 s mass sideband (5450−5590 MeV/c 2 ). The distribution is fitted to the sum of a signal Breit-Wigner shape convolved with a Gaussian resolution function, and a second order polynomial to describe the background contribution. The Breit-Wigner width is set to 0.92 MeV/c 2 [21], and the Gaussian resolution is fixed to 3.8 MeV/c 2 based on simulation. A signal yield of 20.0 ± 5.1 signal events is observed at a mass difference of 565.1 ± 1.0 MeV/c 2 , which is consistent with the known D s1 (2536) + − D + s mass difference of 566.63 ± 0.35 MeV/c 2 [21]. The significance of the signal is 5.9, obtained by fitting the invariant mass distribution with the mean mass difference fixed to 566.63 MeV/c 2 [21], and computing −2ln(L 0 /L max ). Here, L max and L 0 are the fit likelihoods with the signal yields left free and fixed to zero, respectively. Several variations in the background shape were investigated, and in all cases the signal significance exceeded 5.5. This decay is therefore observed for the first time. To obtain the yield in the normalization mode (B 0 s → D + s π − π + π − ), the signal function is integrated from 40 MeV/c 2 below to 40 MeV/c 2 above the nominal B 0 s mass. A yield of 5505 ± 85 events is found in this restricted mass interval.

Selection efficiencies
The ratios of branching fractions can be written as and where Y are the measured yields, ǫ s are the relative selection efficiencies (including trigger), and f s /f d = 0.267±0.021 [27] is the B 0 s fragmentation fraction relative to B 0 . The ratios of selection efficiencies are obtained from simulation, except for the PID requirements, which are obtained from a dedicated D * + calibration sample, weighted to match the momentum spectrum of the particles that form X d and X s . The selection efficiencies for each decay are given in Table 2. The efficiency of the B 0 s → D + s π − π + π − decay is about 35% larger than the values obtained in either the B 0 s → D + s K − π + π − or B 0 → D + s K − π + π − decay; the efficiencies of the latter two are consistent with each other. The lower efficiency is due almost entirely to the tighter PID requirements on the K − and π + in X s . Two additional multiplicative correction factors, also shown in Table 2, are applied to the measured ratio of branching fractions in Eqs. 1 and 2. The first is a correction for the D + s mass veto on M(X d,s ), and the second is due to the requirement that M(X s,d ) < 3 GeV/c 2 . The former, which represents a small correction, is estimated from the sWeight-ed distributions of M(X d,s ) shown previously. For the latter, the fraction of events with M(X d,s ) > 3 GeV/c 2 is obtained from simulation, and scaled by the ratio of yields in data relative to simulation for the mass region 2.6 < M(X s,d ) < 3.0 GeV/c 2 . A 50% uncertainty is assigned to the estimated correction. Based on the qualitative agreement between data and simulation in the M(X d,s ) distributions (see Sect. 5), and the fact that the phase space approaches zero as M(X d,s ) → 3.5 GeV/c 2 , this uncertainty is conservative. The relative efficiency between B 0 s → D s1 (2536) + π − , D + s1 → D + s π − π + and B 0 s → D + s π − π + π − is estimated from simulation, and is found to be 0.90 ± 0.05.

Systematic uncertainties
Several uncertainties contribute to the ratio of branching fractions. The sources and their values are listed in Table 3. The largest uncertainty, which applies only to the ratio , is from the b hadronization fraction, f s /f d = 0.267 ± 0.021 [27], which 4.97 ± 0.08 3.67 ± 0.10 3.59 ± 0.10 D + s veto corr.
1.02 ± 0.01 1.04 ± 0.02 1.14 ± 0.07 is 7.9%. Another large uncertainty results from the required correction factor to account for the signal with M(X s,d ) > 3 GeV/c 2 . Those corrections are described in Sect. 7.
The selection efficiency depends slightly on the modeling of the X d,s decay. The momentum spectra of the B, D + s , X d,s and the X d,s daughters have been compared to simulation, and excellent agreement is found. The selection efficiency is consistent with being flat as a function of M(X d,s ) at the level of two standard deviations or less. To assess a potential systematic uncertainty due to a possible M(X d,s )-dependent efficiency, the relative differences between the nominal selection efficiencies and the ones obtained by reweighting the measured efficiencies by the X d,s mass spectra in data are computed. The relative deviations of 0.5%, 1.1% and 1.2% for B 0 s → D + s K − π + π − , B 0 s → D + s π − π + π − and B 0 → D + s K − π + π − , respectively, are the assigned uncertainties. The systematic uncertainty on the BDT efficiency is determined by fitting the B 0 s → D + s π − π + π − mass distribution in data with and without the BDT requirement. The efficiency is found to agree with simulation to better than the 1% uncertainty assigned to this source. In total, the simulated efficiencies have uncertainties of 1.6% and 1.9% in the two ratios of branching fractions. The PID efficiency uncertainty is dominated by the usage of the D * + calibration sample to determine the efficiencies of a given PID requirement [28]. This uncertainty is assessed by comparing the PID efficiencies obtained directly from simulated signal decays with the values obtained using a simulated D * + calibration sample that is re-weighted to match the kinematics of the signal decay particles. Using this technique, an uncertainty of 2% each on the B 0 s → D + s K − π + π − and B 0 → D + s K − π + π − PID efficiencies is obtained, which is 100% correlated, and a 1% uncertainty for B 0 s → D + s π − π + π − . The trigger is fully simulated, and given the identical number of tracks and the well-modeled p T spectra, the associated uncertainty cancels to first order. Based on previous studies [12], a 2% uncertainty is assigned.
The uncertainties in the signal yield determinations have contributions from both the background and signal modeling. The signal shape uncertainty was estimated by varying all the fixed signal shape parameters one at a time by one standard deviation, and adding the changes in yield in quadrature (0.5%). A double Gaussian signal shape model was also tried, and the difference was negligible. For the combinatorial background, the shape was modified from a single exponential to either the sum of two exponentials, or a linear function. For B 0 s → D + s π − π + π − , the difference in yield was 0.4%. For B 0 s → D + s K − π + π − ,  the maximum change was 4%, and for B 0 → D + s K − π + π − , the maximum shift was 1%. In the B 0 (s) → D + s K − π + π − mass fit, the B 0 (s) → D * + s K − π + π − contribution was modeled using the shape from the B 0 s → D + s π − π + π − mass fit. To estimate an uncertainty from this assumption, the data were fitted with the shape obtained from B 0 s → D * + s K − π + π − simulation. A deviation of 5.5% in the fitted B 0 → D + s K − π + π − yield is found, with almost no change in the B 0 s → D + s K − π + π − yield. The larger sensitivity on the B 0 yield than the B 0 s yield arises because these background contributions have a rising edge in the vicinity of the B 0 mass peak, which is far enough below the B 0 s mass peak to have negligible impact. These yield uncertainties are added in quadrature to obtain the values shown in Table 3. The uncertainties due to the finite simulation sample sizes are 3.0%.
The major source of systematic uncertainty on the branching fraction for B 0 s → D s1 (2536) + π − , D + s1 → D + s π − π + , is from the relative efficiency (5%), and on the fraction of events with M > 3 GeV/c 2 (10%). This 10% uncertainty is conservatively estimated by assuming a flat distribution in M(X d ) up to 3 GeV/c 2 , and then a linear decrease to zero at the phase space limit of ∼3.5 GeV/c 2 . Other systematic uncertainties related to the fit model are negligible. Thus in total, a systematic uncertainty of 11% is assigned to the ratio B(B 0 s → D s1 (2536) + π − , D + s1 → D + s π − π + )/B(B 0 s → D + s π − π + π − ).