Measurement of the absolute branching fraction of the inclusive decay $D_s^+\to \pi^+\pi^+\pi^- X$

Using an $e^+ e^-$ collision data sample with a total integrated luminosity of $3.19$ fb$^{-1}$ collected with the BESIII detector at a center-of-mass energy of 4.178 GeV, the branching fraction of the inclusive decay of the $D_s^+$ meson to final states including at least three charged pions is measured for the first time to be ${\cal B}(D_s^+\to\pi^+ \pi^+ \pi^- X) = (32.81 \pm 0.35_{\rm stat} \pm {0.63_{\rm syst}})\%$. In this measurement the charged pions from $K_S^0$ meson decays are excluded. The partial branching fractions of $D_s^+\to\pi^+ \pi^+ \pi^- X$ are also measured as a function of the $\pi^+ \pi^+ \pi^-$ invariant mass.

Using an e + e − collision data sample with a total integrated luminosity of 3.19 fb −1 collected with the BESIII detector at a center-of-mass energy of 4.178 GeV, the branching fraction of the inclusive decay of the D + s meson to final states including at least three charged pions is measured for the first time to be B(D + s → π + π + π − X) = (32.81 ± 0.35stat ± 0.63syst)%. In this measurement the charged pions from K 0 S meson decays are excluded. The partial branching fractions of D + s → π + π + π − X are also measured as a function of the π + π + π − invariant mass.

I. INTRODUCTION
Lepton flavor universality (LFU) is a very hot topic in the field of particle physics. Recent reports from several high energy experiments studying B meson semileptonic decays of the form B → D ( * ) ℓ + ν ℓ (ℓ = e, µ, or τ ) suggest tensions with Standard Model based predictions [1], and have drawn attention from both experimentalists and theorists. Among these results, the LHCb measurement of the ratio of the branching fractions (BFs) of the B 0 meson semileptonic decays, R(D * − ) ≡ B(B 0 → D * − τ + ν τ )/B(B 0 → D * − µ + ν µ ), has attracted a lot of attention. The R(D * − ) measurement with the τ lepton reconstructed from three prongs, i.e. charged pions (τ + → π + π + π − X, where X stands for all possible particle combinations) [2], based on a 3 fb −1 data sample of pp collisions, has one of the smallest statistical uncertainties. However, it is systematically limited and one of the limiting systematic uncertainty is associated with the background from B 0 → D * − D + s (X), followed by D + s → π + π + π − X decays. In the LHCb measurement, the total systematic uncertainty is 9.1%, out of which 2.5% is attributed to the "D + s → 3πX decay model" for these background events [2]. Here, 3π ≡ π + π + π − . A better understanding of the D + s → π + π + π − X decay dynamics and of the related π + π + π − -system kinematics would directly enhance the sensitivity of the R(D * − ) measurement at LHCb, thus becoming more and more important in view of the steady increase of the integrated luminosity collected by the experiment. Similarly, the measurement at LHCb would strongly benefit from the knowledge on the D + s → π + π + π − X decay as input [3].
The inclusive decay D + s → π + π + π − X has never been studied directly in experiments. Table I lists the BFs of major exclusive D + s channels with at least three charged pions of π + π + π − in the final states that have been measured with relatively good precision to date [4,[6][7][8][9][10]. The decays D + s → η (′) π + (X) make up about half of the contributions. Measurements on the inclusive decay rate of D + s into into final states with at least one combination of π + π + π − , B(D + s → π + π + π − X), will shed light on yet unobserved decay modes with rich hadronic contents, such as D + s → η ′ π + π + π − . Charge-conjugate states are implied throughout this paper.
In this paper, the first measurement of the BF of the inclusive decay of D + s → π + π + π − X based on 3.19 fb −1 of e + e − collision data collected at a center-of-mass energy (E cm ) of 4.178 GeV with the BESIII detector in 2016 is presented. The paper is organized as follows. Section II introduces the BESIII detector as well as the data and Monte Carlo (MC) simulated samples used in this analysis. Section III gives an overview of the analysis technique. Event selection requirements to fully reconstruct D − s decays are described in Sec. IV, while the selection requirements for the 3π system at the signal D + s side and the further analysis based on the selected signal candidates are presented in Sec. V. The systematic uncertainties on our measurements are evaluated in Sec. VI, and the final results are summarized in Sec. VII.

II. DETECTOR AND DATA SETS
The BESIII detector [11] records the final state particles of symmetric e + e − collisions provided by the BEPCII storage ring [12] in the E cm range from 2.00 to 4.95 GeV. BESIII has collected large data samples in this energy region [13], while this analysis uses the entirety of the 3.19 fb −1 data sample collected at E cm = 4.178 GeV in 2016. The cylindrical core of the BESIII detector covers 93% of the full solid angle and consists of a heliumbased multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field [14]. The solenoid is supported by an octagonal flux-return yoke instrumented with resistiveplate-counter muon identification modules interleaved with steel. The MDC reconstructs charged-particle momenta with a resolution of 0.5% at 1 GeV/c, measuring the characteristic energy loss (dE/dx) with 6% precision. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end-cap) region. The time resolution in the TOF barrel region is 68 ps. The end-cap TOF system was upgraded in 2015 using multigap resistive plate chambers, providing a time resolution of 60 ps [15].
Simulated data samples are produced with a geant4 [16] based MC simulation toolkit, which includes the geometric description of the BESIII detector and the detector response. The simulation also exploits the kkmc [17] generator to take into account the beam energy spread and the initial-state radiation in the e + e − annihilations. The MC simulation samples, referred to as "inclusive MC simulation", include the production of D + s particles via the process e + e − → D , and other open charm processes, as well as the production of vector charmonium(like) states, and the continuum processes incorporated in kkmc. The known decay modes are modeled with evtgen [18] using the BFs taken from the Particle Data Group [4], and the remaining unknown decays from the charmonium states with lundcharm [19]. Final state radiation from charged final state particles is incorporated using photos [20]. The effective luminosities of the inclusive MC simulation samples correspond to 40 times the data luminosity. The inclusive MC simulation samples are used to estimate background contributions. A subset of the D s inclusive MC simulation sample, with one D + s meson decaying into final state modes including at least three charged pions, is also used as the signal MC simulation sample of D + s → π + π + π − X to determine detection efficiencies.

III. ANALYSIS TECHNIQUE
At E cm = 4.178 GeV, D ± s mesons are produced in pairs predominantly through the process of e + e − → D * ± s D ∓ s . The D * ± s mesons then primarily decay into D ± s γ (93.5%) and D ± s π 0 (5.8%) [4]. The double-tag (DT) technique, first employed by the Mark III collaboration [21], is used to determine the BF for the inclusive decay mode of the D + s meson. We first fully reconstruct a D − s meson, named "single tag (ST)", in one of the hadron tag modes TABLE I. Major exclusive decay modes that contribute to the inclusive decay of D + s → π + π + π − X, along with the known BFs. The modes with BFs of large uncertainties (K + η ′ , π + π + π + π − π − π 0 (without contributions from η ′ /ω → π + π − π 0 ), etc.) are not shown. None of the charged pions in the 3π system comes from K 0 S meson decays. The values of B(M → π + π − X), M = η, ω, η ′ , ϕ and B(τ + → π + π + π − X) are estimated by combining all known exclusive decays that contribute [4]. listed in Sec. IV [22]. Then in the ST D − s sample, we reconstruct the D + s → π + π + π − X signal in the side recoiling against the D − s candidate, referred to as the DT. The number of ST candidates for a specific tag mode α (N α ST ) determined from the fit shown in Fig. 1 is used to obtain the total number of e + e − → D ( * )± s D ∓ s events in the data (N tot ) by where B α is the BF for the D − s tag mode α, and ϵ α ST is the ST detection efficiency.
Then the signal decay which contains at least three charged pions is selected at the recoil (signal) side. The tag side and signal side selection criteria are described in detail in Sec. IV and Sec. V, respectively. To ensure one 3π system per signal decay, in the case of more than one π − (two π + s) at the signal side, we select the π − (two π + s) with the largest momentum (largest and second largest momenta), in the laboratory frame. The true (reconstructed) 3π invariant mass M (π + π + π − ) is calculated based on true (reconstructed) pion momenta.
Since the reconstruction efficiency is dependent on the kinematics of the 3π system at the signal side, which is not well known, the partial BFs of the inclusive D + s → π + π + π − X decay are measured in intervals of M (π + π + π − ). The efficiency is parametrized as a function of M (π + π + π − ), and its dependence on other kinematic variables is later examined in Sec. VI G. For tag mode α, the number of produced DT events and the numbers of observed DT events are related in M (π + π + π − ) intervals through a detector response matrix that accounts for detector efficiency and detector resolution for M (π + π + π − ), where N α prod,j (N α obs,i ) is the number of signal events produced (observed) in the j th (i th ) M (π + π + π − ) interval. The ij th element of the efficiency matrix, ϵ α ij , describes the efficiency and migration across the M (π + π + π − ) intervals. It is calculated as ϵ α ij ≡ N α obs,ij /N α prod,j using simulated signal D + s → π + π + π − X decays with D − s decaying into the tag mode α. Here, N α obs,ij is the number of signal MC simulated events produced in the j th interval and reconstructed in the i th interval of the M (π + π + π − ) distribution.
The number of observed D + s → π + π + π − X signal candidates in each M (π + π + π − ) interval is determined by fitting the invariant mass distribution of ST D − s candidates to extract the raw signal yield N α raw,i . The yields of expected K 0 S background contributions, N α K 0 S ,i , and contributions from particle misidentification (misID), N α misID,i , estimated from MC simulation, are subtracted from the raw yield N α Here, N α misID,i = N α e→π,i + N α µ→π,i + N α K→π,i includes contributions from the three misID scenarios described in Sec. V B.
For tag mode α, the number of D + s → π + π + π − X signal candidates produced in the i th interval of the M (π + π + π − ) distribution is then obtained by solving Eq. (2), which gives and its statistical uncertainty is given by where σ stat N α obs,j is the statistical uncertainty of N α obs,j . The statistical uncertainties of ϵ α ij due to the limited size of the signal MC simulation sample are later considered as a source of systematic uncertainties, as discussed in Sec. VI H.
Taking N tot from Eq. (1), the corresponding partial BF of D + s → π + π + π − X is given by In each tag mode α, the partial BFs of the D + s → π + π + π − X decays are summed to obtain the total BF B α (D + s → π + π + π − X) = i ∆B α 3πX,i . The total BF and its associated statistical uncertainty are and respectively, where σ B α is the statistical uncertainty of B α (D + s → π + π + π − X).

IV. SINGLE TAG ANALYSIS
ST D − s candidates are selected in two hadronic tag modes, Other tag modes are not considered due to their relatively high background level.
All charged track candidates detected in the MDC must be within a polar angle (θ) range of |cos θ| < 0.93, where θ is defined with respect to the z axis, which is the symmetry axis of the MDC. For charged tracks not originating from K 0 S decays, the distance of closest approach to the interaction point (IP) is required to be less than 10 cm along the z axis, and less than 1 cm in the transverse plane.
Charged tracks are identified as pions or kaons with the particle identification algorithm (PID), which combines measurements of dE/dx in the MDC and the time of flight in the TOF to form likelihoods L(h)(h = K, π) for different hadron h hypotheses. Charged kaons and pions are identified by comparing the likelihoods for the kaon and pion hypotheses, L(K) > L(π) and L(π) > L(K), respectively. Furthermore, charged pion candidates with momenta below 100 MeV/c are rejected to suppress soft pions from D * + → D 0 π + .
The K 0 S candidates are reconstructed from two oppositely charged tracks each with the distance of closest approach to the IP less than 20 cm along the z axis. The two charged tracks are assigned as π + π − without imposing further PID criteria. They are constrained to originate from a common vertex and are required to have an invariant mass satisfying S are the invariant mass of the pion pair and the known K 0 S mass [4], respectively. For each D − s candidate, the reconstructed invariant mass M (D − s ) is required to be in the range 1.90 and 2.05 GeV/c 2 . The recoil mass of the D − s candidate is defined as where ⃗ p Ds is the reconstructed momentum of the D − s candidate in the e + e − center-of-mass frame, and m Ds is the known D − s mass [4].
Finally, for either tag mode, in about 6% of the cases, more than one D − s candidate per event passes the tag side event selection requirements. In these cases, only the D − s candidate with M rec closest to the known D * + s mass per event per tag mode is retained for further study. The systematic uncertainties related to this single-candidate requirement are discussed in Sec. VI K.
For each tag mode, the number of correctly reconstructed ST D − s candidates is determined by performing an unbinned maximum likelihood fit to the distribution of M (D − s ), as shown in Fig. 1. The combinatorial background is modeled using a linear probability density function (PDF). The shapes of the D − background con- The yields of the D − meson background contributions are estimated with the simulation and fixed in the fits. The signal lineshapes are modeled by the sum of a Gaussian PDF and a double-sided Crystal-ball (DSCB) PDF [23] with a common mean value m 0 . The Gaussian PDF has a width of σ 0 that is allowed to float along with m 0 . The DSCB PDF tail parameters and the width ratio between the DSCB and Gaussian PDFs are determined from signal MC simulated events and fixed in the fits to data, as in Ref. [24]. The fitted yields are summarized in Table II, along with the ST efficiencies.
. Data (black points) are shown overlaid with the fit results including the total (solid blue), signal PDF (dashed magenta), D − background PDF (dotted red), and combinatorial background PDF (long-dashed blue) components.

V. DOUBLE TAG ANALYSIS
In each event, once an ST D − s candidate is found, three charged pions with the sum of charges opposite to the charge of the D − s tag in the rest of the event (recoil side) are searched for. The tracking and PID requirements of charged pions are the same as those in the tag side as described in Sec. IV.
There are two major background sources that contribute to the 3π inclusive final states. The first source of background contribution is due to events with three correctly identified charged pions where at least one of them comes from K 0 S → π + π − decays. The other type of background comes from events for which at least one charged pion is misidentified. These two types of background sources are studied separately using dedicated MC simulations of inclusive D + s decays to determine additional selection criteria, and to estimate the remaining contribution after imposing all the requirements. The background subtracted pion momentum spectra and mass spectra from the 3π system for the tag mode D − s → K − K + π − are shown in Figs. 2 and 3, while the spectra for the tag mode D − s → K 0 S K − follow very similar patterns with lower yields. There is reasonable agreement between data and the simulated sample, and any differences in the pion momentum spectra are later considered as a source of systematic uncertainties.
To suppress the background contribution due to K 0 S decays, each of the three charged pions of the 3π system is checked against a list of K 0 S candidates in an event. Each of the K 0 S candidates is reconstructed with the criteria defined in Sec. IV . If at least one pion is part of one reconstructed K 0 S candidate that has L/σ L > 2 (L is defined as the distance between the K 0 S vertex and the IP, and σ L is the uncertainty of L), or within the mass window |M (π + π − ) − m K 0 S | < 10 MeV/c 2 , this signal candidate is rejected. With this requirement, about 88% of the K 0 S meson related background contributions are rejected, while losing about 11% of the signal events.

B. Particle misidentification background
For one charged pion, there are three possible misID scenarios: • e ± → π ± , Normalized momentum spectra after background subtraction from the 3π system of the inclusive decay D + s → π + π + π − X for the tag mode D − s → K − K + π − . The points are obtained from data and the solid line histograms are from inclusive MC simulation samples. The subscripts for π + s indicate the order p(π + 1 ) > p(π + 2 ).
The points are obtained from data and the solid line histograms are from inclusive MC simulation samples. The subscripts for the two π + π − combinations indicate the order M high (π + π − ) > M low (π + π − ). The dips in the M (π + π − ) spectra are due to the 10 MeV/c 2 K 0 S mass veto defined in Sec. V A.
As the charged pions are identified against kaons, the main misID background contributions are from the semileptonic D + s decays, and the misidentified pion in the 3π system has the same charge as the signal D + s side. Electron background can be largely removed by rejecting signal candidates for which at least one e + has been identified at the signal D + s side. The electron PID is based on the likelihoods for the electron, pion, and kaon hypotheses L(h)(h = e, π, K) using information from the MDC, TOF, and EMC, as in Ref. [27]. With this, about 90% of the e + → π + misID background contributions are removed at the cost of mere ∼5% of the signal events. The residual misID backgrounds are then dominated by the semimuonic decays D + s → Xµ + ν µ , and are estimated from the D s inclusive MC simulation sample (∼70%), , and the related systematic uncertainties are discussed in Sec. V B.

C. Partial branching fractions
The D + s → π + π + π − X partial BFs are measured in eleven M (π + π + π − ) intervals. The intervals are populated with approximately equal numbers of DT candidates, except for the last interval, where the signal candidates are predominantly from the exclusive decay of D + s → π + π + π − . The interval boundaries in the M (π + π + π − ) distribution are chosen as The number of D + s → π + π + π − X signal candidates in each M (π + π + π − ) interval is determined by fitting the corresponding M (D − s ) distribution of the tag side. From this fit, the raw signal yield N α raw,i is measured for each interval i of the M (π + π + π − ) distribution, for tag mode α. The signal shape for each tag mode is determined in the corresponding ST fit (results shown in Fig. 1). The same shape is used in all M (π + π + π − ) intervals, and the linear background shape parameter is left unconstrained in the fit. The fitted raw signal yields (N α raw,i ) are presented in Tables III and IV for the tag modes . Two examples of the 22 fits are shown in Fig. 4, and the full set of fit results is made available as Supplemental Material [28]. Tables III and IV summarize the observed signal yields, after subtracting the K 0 S -decay related and misID background contributions from the raw signal yields as described earlier in this section.
Using the observed signal yields and the efficiency matrix elements ϵ α ij presented in Tables V and VI, the partial BFs in various M (π + π + π − ) intervals are calculated using Eqs. 4 and 6 and shown in Table VII after combining the measurements from both tag modes. III. Raw signal yields and estimated background contributions for D − s → K 0 S K − versus D + s → π + π + π − X events. The associated uncertainties are statistical only. The uncertainties from different background contributions are due to the limited size of our inclusive MC simulation samples, whose integrated luminosity is 40 times of data.
M (π + π + π − ) interval Raw signal yields and estimated background contributions for D − s → K − K + π − versus D + s → π + π + π − X events. The associated uncertainties are statistical only. The uncertainties from different background contributions are due to the limited size of our inclusive MC simulation samples, whose integrated luminosity is 40 times of data.

VI. SYSTEMATIC UNCERTAINTIES
The methods to determine the systematic uncertainties in the measured total and partial BFs of inclusive D + s → π + π + π − X decays are described below. The signal M (D − s ) PDF is changed from the DSCB shape to a shape modeled with signal MC simulated events convolved with a Gaussian function. The width and mean parameters of the Gaussian function are floating in the fits to account for the resolution difference between data and simulation. Alternatively, the fits are performed by allowing the two signal shape parameters TABLE VI. Efficiency matrix ϵij in percentage for D − s → K − K + π − versus D + s → π + π + π − X determined from signal MC simulated events. Each column gives the true M (π + π + π − ) interval j, while each row gives the reconstructed M (π + π + π − ) interval i. m 0 and σ 0 in different M (π + π + π − ) intervals to differ. The largest variations are taken as the relative systematic uncertainties.
The background shape is changed from the linear PDF to an exponential or a histogram PDF based on MC simulation. Also the fits are performed with the yields TABLE VII. Partial BFs (in percentage) of the decay D + s → π + π + π − X combined from the two tag modes. The first and second uncertainties are statistical and systematic, respectively.

C. Pion tracking and PID efficiency
The π ± tracking efficiency is studied with the control sample of e + e − → K + K − π + π − decay. The data-MC tracking efficiency ratios of π + and π − are measured in bins of pion transverse momentum. The signal MC samples are then weighted by these ratios to recalculate signal efficiency matrices. This leads to 0.98% variation in total BF, which is taken as the systematic uncertainty associated with tracking efficiency.
The π ± PID efficiency is studied with the control samples of e + e − → K + K − π + π − (π 0 ) and π + π − π + π − (π 0 ) decays. The data-MC tracking efficiency ratios of π ± are measured in bins of pion momentum. The signal MC samples are then weighted by these ratios to recalculate signal efficiency matrices. This leads to 0.99% variation in the total BF, which is taken as the systematic uncertainty associated with PID efficiency. (3) is determined by the D s inclusive MC simulation sample, the data and simulation differences on the D + s → K 0 S π + X contributions are studied by using control samples with dedicated stringent requirements on the 3π system: one π + π − pair forms the reconstructed K 0 S that has L/σ L > 2 and |M (π + π − ) − m K 0 S | < 10 MeV/c 2 . These requirements lead to a data sample with more than 99% of the D + s candidates having a π + π − pair from K 0 S . The agreement on the observed number of D + s candidates in data and the number from simulation is at ∼10% level across different M (π + π + π − ) intervals. Therefore, the N α K 0 S ,i values are all scaled by ±10% in turn, and the BFs are recalculated according to altered N α obs,i values based on Eq. (3). The larger variations from the two measurements are taken as the relative systematic uncertainties on the K 0 S background.
Furthermore, the systematic uncertainty related to the K 0 S veto of |M (π + π − ) − m K 0 S | < 10 MeV/c 2 is examined by varying the mass window by ±5 MeV/c 2 . The resulting changes in the BFs are smaller than the statistical uncertainties on the differences, so no uncertainty is assigned for this source according to Ref. [25].

E. Particle misidentification background
The D s inclusive MC simulation sample is also used to estimate N α misID,i defined in Eq. (3). The uncertainties on N α misID,i come from limited knowledge on the BFs of different decays that contribute to the misID backgrounds, and differences in misID rates between data and simulation, that are considered separately in the three misID scenarios as listed in Sec. V B.
The D + s → π + π − ℓ + X decays form the background where the leptons ℓ = e, µ are misidentified as pions, and the decays D + s → K − π + π + X and D + s → K + π + π − X are relevant for the background when a kaon is misidentified instead. The D + s → π + π − ℓ + X decays are dominated by the semileptonic decays of D + s → ϕℓ + ν ℓ and D + s → η (′) ℓ + ν ℓ . By assuming LFU in semileptonic charm decays, the data and simulation differences on modeling the D + s → π + π − ℓ + X decays are studied using samples of D + s → π + π − e + X. This LFU assumption has been found to hold experimentally, e.g. within a precision of 1.4% in D 0 → K − ℓ + ν ℓ decays [26]. The electron PID requirement used on one of the positively charged track in the 3π system is the same as that in Ref. [27]. As the relative difference between data and MC simulation is below 5%, the sums of N α e→π, i + N α µ→π, i are all scaled by ±5% in turn. The larger variations from the two measurements are taken as the systematic uncertainties related to the lepton misID background. Similarly, the dedicated simulation samples of D + s → K − π + π + X and D + s → K + π + π − X decays can be compared to data by imposing kaon identification requirements on one track from the 3π system. This study indicates a disagreement between data and simulation at a level of about 15%. This is significantly larger than in the lepton case due to the poor knowledge of hadronic D + s decays with at least one charged kaon. The N α K→π, i values are all scaled by ±15% in turn, and the larger variations from the two measurements are taken as the systematic uncertainties related to the kaon misID. Separate studies on the probabilities of tracks being misidentified as pions show the effects due to data and simulation differences are negligible. Therefore, this is not included as part of the systematic uncertainties.

F. Pion multiplicity
The number of charged pions in D + s → π + π + π − X decays may be inaccurately modeled in the simulation mainly due to the limited knowledge of hadronic D + s decays with five or more charged pions in the final states such as D + s → π + π + π − π + π − π 0 (π 0 ) [4]. The signal MC simulated events are weighted so that the simulated distribution of the number of charged pion candidates at the signal side matches with that in the data. This causes changes in signal efficiencies, and the resulting BF variations using the weighted MC simulated events are taken as the relative systematic uncertainties.

G. Pion momenta
The signal MC simulated events are weighted to improve the agreement of the pion momentum distributions between data and simulation. The multivariate algorithm [29] to determine the weights is trained on the distributions from unweighted signal MC simulation and background subtracted data samples as shown in Fig. 2. Using the weighted MC candidates based on the threedimensional weighting scheme to recalculate signal efficiency matrices for both tag modes, the variations in the BF central values are taken as the relative systematic uncertainties.
As the data-MC differences in the angular distributions of pion momenta may result in differences in the mass spectra of the 3π system, we also adopt a five-dimensional weighting scheme with the addition of M high (π + π − ) and M low (π + π − ) to account for the data-MC differences shown in Fig. 3. The resulting variations in BFs are smaller than those from the three-dimensional weighting scheme, and thus no additional uncertainty is assigned.

H. Monte Carlo simulation sample size
The BFs are recalculated 1000 times by randomly perturbing the efficiency matrix elements shown in Tables V and VI according to their uncertainties. The distributions of the recalculated BF central values are fitted with a Gaussian function and the fitted Gaussian widths are taken as the relative systematic uncertainties.

I. Measurement bias
The uncertainties due to the entire analysis procedure to extract the BFs are estimated by using 40 inclusive MC simulation samples, each with a total candidate number statistically matched to that in the data. While good agreement between the input and measured BF values is found, the mean biases are taken as the relative systematic uncertainties.

J. Signal MC simulation model
The purpose of measuring B(D + s → π + π + π − X) in intervals of M (π + π + π − ) is to be independent of the D + s decay model used in the signal simulation. A number of exclusive D + s decay modes such as D + s → η ′ ρ + , D + s → ϕπ + π + π − , and D + s → ϕρ + are used to perform the tests. The yield of the respective signal MC simulation events is weighted by about 40%. This in effect changes the related input BF by about 40% (up or down). The signal efficiency matrices are then updated based on the weighted MC simulation samples. The obtained BF results show good agreement with the baseline results well within one standard deviation. Therefore, a systematic uncertainty related to the signal MC simulation model is not assigned.

K. Single-candidate requirement
The single-candidate requirement described in Sec. IV may bias the measurement due to potentially different candidate multiplicities and M rec distributions in data and MC simulation. The systematic uncertainties from the single-candidate requirement are examined by rerunning the measurement with this requirement removed. The resulting BF variations are smaller than the statistical uncertainties on the differences, so no uncertainty is assigned for this source.
L. Summary of systematic uncertainties Table VIII summarizes contributions from the different systematic sources in the measurement of B(D + s → π + π + π − X). The largest sources of systematic uncertainties come from tracking and PID of the 3π system. These contributions are combined in quadrature to determine the total systematic uncertainty in the BF measurement. The summary table for the systematic uncertainties of the partial BFs of D + s → π + π + π − X is given in the Supplemental Material [28].
TABLE VIII. Sources of systematic uncertainties in the measurement of B(D + s → π + π + π − X).
Based on 3.19 fb −1 of e + e − collision data recorded with the BESIII detector at E cm = 4.178 GeV, the BF of the inclusive D + s → π + π + π − X decay is measured for the first time. It is found to be B(D + s → π + π + π − X) = (32.81 ± 0.35 stat ± 0.63 syst ) %. This is larger than the sum of all observed exclusive BFs of ∼25% based on the PDG and recent measurements as summarized in Table I, hinting at potentially unobserved decay modes with at least three charged pions in the final states.
Furthermore, the partial BFs of the D + s → π + π + π − X decay as a function of M (π + π + π − ) are also measured, as presented in Table VII and Fig. 5. As the feed-up contribution from lower intervals is expected to be negligible (∼1% of exclusive signals), the partial BF in the last M (π + π + π − ) interval is consistent with the exclusive BF of B(D + s → π + π + π − ) = (1.08 ± 0.04)% from the PDG [4], with somewhat lower precision. The measured total and partial BFs of D + s → π + π + π − X offer important inputs to constrain the systematic uncertainties in future LHCb measurements on R(D * − ) and R(Λ + c ) with much larger data samples [30].

Supplementary material
In addition to the results presented in the main body of the paper, we provide a supplementary collection of figures on the fits used in determining the double-tag yields in the two tag modes, and a table summarizing the systematic uncertainties of the partial BFs of D + s → π + π + π − X in intervals of M (π + π + π − ). Table 1: Summary of all different sources of systematic uncertainties for partial branching fractions of D + s → π + π + π − X. The uncertainties are summed up in quadrature to make the total systematic uncertainties.