Amplitude analysis of the $D^+\to\pi^- \pi^+\pi^+$ decay and measurement of the $\pi^-\pi^+$ S-wave amplitude

An amplitude analysis of the $D^+ \to \pi^- \pi^+ \pi^+$ decay is performed with a sample corresponding to $1.5\rm fb^{-1}$ of integrated luminosity of $pp$ collisions at a centre-of-mass energy $\sqrt{s}=8$ TeV collected by the LHCb detector in 2012. The sample contains approximately six hundred thousand candidates with a signal purity of $95\%$. The resonant structure is studied through a fit to the Dalitz plot where the $\pi^- \pi^+$ S-wave amplitude is extracted as a function of $\pi^-\pi^+$ mass, and spin-1 and spin-2 resonances are included coherently through an isobar model. The S-wave component is found to be dominant, followed by the $\rho(770)^0\pi^+$ and $f_2(1270)\pi^+$ components. A small contribution from the $\omega(782)\to\pi^-\pi^+$ decay is seen for the first time in the $D^+ \to \pi^- \pi^+ \pi^+$ decay.


Introduction
determined by a fit to the data, while spin-1 and spin-2 states are included through an isobar model.

LHCb detector and simulation
The LHCb detector [22,23] 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 (VELO) 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 of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV. 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. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillatingpad 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.
The online event selection is performed by a trigger system [24], which consists of a hardware stage, followed by a software stage, which applies a full event reconstruction. At the hardware-trigger stage, events are selected based on information from particles that are not related to the signal by requiring a hadron, photon or electron with high transverse energy in the calorimeter, or a muon with high p T in the muon system. The software trigger is divided into two parts. The first part employs a partial reconstruction of the tracks, and a requirement on p T and IP is applied to, at least, one of the final-state particle forming the D + candidate. In the second part a full event reconstruction is performed and dedicated algorithms are used to select D + candidates decaying into three charged hadrons.
Simulation is used to model the effects of the detector acceptance and the selection requirements, to validate the fit models and to evaluate efficiencies. In the simulation, pp collisions are generated using Pythia [25,26] with a specific LHCb configuration [27]. Decays of unstable particles are described by EvtGen [28], in which final-state radiation is generated using Photos [29]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [30,31] as described in Ref. [32]. For part of the total simulation sample, the underlying pp interaction is reused multiple times, with an independently generated signal decay for each [33].
3 D + → π − π + π + selection In the offline selection, combinations of three charged particles with net charge ±1 must form a good-quality vertex well detached from any PV. The PV associated to the D + candidate is chosen as that with the smallest value of χ 2 IP , where χ 2 IP is defined as the difference in the vertex-fit χ 2 of the PV reconstructed with and without the particle under consideration, in this case the D + candidate. Requirements are placed on the following variables: the significance of the distance between the PV and the D + decay vertex; the angle between the reconstructed D + momentum vector and the vector connecting the PV to the decay vertex; the χ 2 of the D + decay vertex fit; the distance of closest approach between any two final-state tracks; the momentum, the transverse momentum and the χ 2 IP of the D + candidate and of its decay products. The invariant mass of the D + candidates, calculated using the pion mass hypothesis for the three tracks, is required to be within the interval [1800,1940] MeV.
Particle identification (PID) is also used to separate pions from kaons and muons. The requirements, placed on all tracks, reduce to the percent level the cross-feed from decays such as D 0 → K − π + plus an unrelated track, D + → K − π + π + and D + → µ − π + π + ν µ . The D + → K 0 S π + decay is removed by vetoing candidates with π + π − invariant mass in the interval [0.485,0.500] GeV.
A multivariate analysis (MVA) [34] is used to further reduce the remaining backgrounds, mostly from random association of three charged tracks, but also from contamination of D + (s) → (η ( ) → π + π − γ)π + decays, where the photon is not detected. The MVA uses a Gradient Boosted Decision Tree (BDTG) [35] as a classifier and is based on the quantities (or combinations of them) described above, except PID and those that can potentially cause large efficiency variation across the Dalitz plot, such as the χ 2 IP and p T of the decay products. For training, simulated D + → π − π + π + decays are used to represent the signal, whereas data from the invariant-mass sidebands ([1810,1830] and [1910,1930] MeV) of the signal peak are used for the background. Before performing the MVA, a weighting procedure [36] is used to account for small differences between the simulation and data, including differences in the momentum and transverse momentum distributions of the decay products. The data are divided into two independent sets in a pseudorandom manner. The π − π + π + mass spectrum of one of the datasets is fitted using the sPlot technique [37] to obtain signal and background weights needed to determine the distributions of the quantities to be used in the MVA, as well as the kinematic distributions, needed for a weighting procedure explained above. The second dataset, corresponding to an integrated luminosity of about 0.75 fb −1 , is used to select the final candidates for the Dalitz plot analysis.
The requirement on the BDTG output is chosen to yield a sample of D + → π − π + π + decays with 95% purity in order to reduce the impact of systematic effects related to the background modelling in the Dalitz plot fit. The efficiency drops very rapidly for more stringent requirements, with only a modest gain in purity. The invariant mass of D + → π − π + π + candidates after all requirements is shown in Fig. 2. An extended binned maximum-likelihood fit to this distribution is performed. The probability distribution function (PDF) of the signal is represented by a sum of a Gaussian function and two Crystal Ball (CB) [38] functions, while the background is modelled by an exponential function. The signal PDF is where µ and σ G are the mean value and the width of the Gaussian function G. The two Crystal Ball functions, CB 1 and CB 2 , have widths R 1 σ G and R 2 σ G , and tail parameters α 1 , N 1 and α 2 , N 2 . A common parameter, µ, describes the most probable mass value of the two Crystal Balls and the mean of the Gaussian function. The fractions of each PDF component are f G for the Gaussian function, The parameters α i , N i , R i , f CB and f G defining the signal PDF are fixed to the values obtained from a fit to the simulated sample, while the parameters µ and σ G are allowed to float freely. For the Dalitz Plot analysis, candidates are selected within a 2σ eff mass window around the mean µ, where the effective mass resolution, σ eff , is defined by MeV. The final sample has 601 171 candidates with a signal purity of (95.2 ± 0.1)% in the corresponding mass interval of [1854.1,1889.0] MeV. Multiple candidates in the same event correspond to only 0.15% of the final sample and are retained.
The Dalitz plot distribution of the selected candidates is shown in Fig. 3. The axes are the Dalitz variables s 12 ≡ (p 1 + p 2 ) 2 and s 13 ≡ (p 1 + p 3 ) 2 , where p i (i = 1, 2, 3) are the 4-momenta of the three pions in the final state and the particle ordering is such that the pion with charge opposite to that of the D + meson is always particle 1, and the same-sign pions are randomly assigned particles 2 and 3. These Lorentz-invariant quantities are computed after refitting tracks with the constraint on the invariant mass of the candidate to the known D + mass [1].

Background and efficiency models
The remaining background in the D + → π − π + π + sample, which amounts to 4.8%, needs to be parameterised to be included in the Dalitz plot fit. The background distribution is inferred from the sidebands of the D + → π − π + π + invariant-mass signal region, specifically, the intervals [1810,1830] MeV and [1910,1930] MeV. The background within the signal region is assumed to be the average composition of both sidebands. The candidates in these regions are projected to the Dalitz plot and a two-dimensional (2D) cubic spline procedure [39] is used to smooth the distribution in order to avoid binning discontinuities. The resulting model, shown in Fig. 4 (left), is used in the Dalitz plot fit.
The signal distribution in the Dalitz plot includes efficiency effects, which need to be corrected for. The efficiency model across the Dalitz plot includes the effects of the geometrical acceptance of the detector, as well as reconstruction, trigger, selection, and  PID requirements. All these effects apart from those associated with PID are quantified using simulation. The PID efficiencies for the pions are evaluated from calibration samples of D * + → D 0 (→ K − π + )π + decays [40] and depend on the particle momentum, pseudorapidity and charged-particle multiplicity. The final efficiency model is constructed from a two-dimensional histogram with 15 × 15 uniform bins which is then smoothed by a 2D cubic spline, as shown in Fig. 4 (right).

The QMIPWA formalism
The Dalitz plot of D + → π − π + π + candidates shown in Fig. 3 shows a rich resonant structure: contributions from the subchannels ρ(770) 0 π + (with an evident interference with ω(782)π + ), f 0 (500)π + , f 0 (980)π + and possible high-mass vector and tensor states are seen. To disentangle these contributions, a full amplitude analysis is needed. The approach of this work is to describe the total D + → π − π + π + amplitude as a coherent sum of an S-wave contribution and higher-spin waves, where the total amplitude is Bose-symmetrised with respect to s 12 and s 13 due to the two identical pions. The first term is the S-wave amplitude, where the real functions a 0 (s 12 ) and δ 0 (s 12 ) are to be determined by the Dalitz plot fit: the π − π + invariant-mass spectrum is divided in 50 intervals (knots) where interpolation to obtain a continuous S-wave complex function is attained through a linear spline. The intervals are not uniformly distributed to allow better determination of the magnitude and phase variations in regions where they are expected to vary most -such as near the f 0 (980) resonance. Variations in the knot spacing and number of knots are addressed for systematic uncertainties. The P-and D-wave components are included through an isobar model, represented by the terms in the sum in Eq. 2, where A i (s 12 , s 13 ) is the complex amplitude of resonance R i , with magnitude a i and phase δ i as free parameters. This approach to the total amplitude is referred to as quasi-model-independent since any limitation of the isobar model to describe the higher-spin components may reflect in the description of the S-wave amplitude which has 100 free parameters. Within the isobar model, the individual resonant amplitude for a process of type

Spin
Blatt-Weisskopf factor Since the D + meson and the final-state pions are spinless particles, the spin of the resonance, J, is equal to the orbital angular momentum in both the decays D → R i c and The Blatt-Weisskopf barrier factors [41], F J D and F J R , take into account the finite size of the D meson and the R i resonance in the decay processes D → R i c and R i → a b, respectively. They are functions of z = rp * , as shown in Table 1, where r is the effective radius of the decaying particle, p * is the modulus of the momentum of the decay products measured in the decaying particle rest frame, and z 0 is the value of z calculated at the known mass of the decaying particle. In this analysis, the effective radii of the D + meson and the intermediate resonances are set to 5 GeV −1 and 1.5 GeV −1 , respectively [15,[42][43][44].
The function M J describes the angular distribution of the decay particles in the Zemach formalism [45] with explicit forms for J = 1 and J = 2 given by and respectively, where m D is the known mass of the D meson and m a , m b and m c those of the decay products. It is clear that these expressions take a simpler form for the particular case when there are three charged pions in the final state (m a = m b = m c = m π ). The dynamical function T R i (s ab ) in Eq. 4 represents the resonance lineshape, usually a relativistic Breit-Wigner (RBW), where m 0 is the nominal mass of the resonance and Γ(s ab ) is the mass-dependent width which is given by with Γ 0 being the nominal width of the resonance. A Gounaris-Sakurai (GS) function [46] is a modification of the RBW lineshape, commonly used to describe the pion electromagnetic form factor in the parameterisation of spin-1 ρ-type resonances, where The function h(s ab ) is given by where m π is the pion mass and the derivative is given by The In this analysis, the GS lineshape is used for the spin-1 ρ-type resonances, while the RBW lineshape is used for other resonances such as the ω(782) and f 2 (1270) states. The values m 0 and Γ 0 for all resonances are fixed in the fit to their known values [1].
The ω(782) → π − π + decay violates isospin, and it is not clear whether this process occurs through direct decay or through mixing with the ρ(770) 0 state (or both). As an alternative to representing the ρ(770)π + and ω(782)π + amplitudes as a sum of isobars, their combined contribution is parametrised through a ρ − ω mixing lineshape given by [47] T where T ρ and T ω are the GS and RBW lineshapes for the ρ(770) 0 and ω(782) resonances, respectively. The magnitude |B| and the phase φ B quantify the relative contribution of the ω(782) and the ρ(770) resonances, and are free parameters in the fit. The factor ∆ = δ(m ρ 0 + m ω 0 ) governs the electromagnetic mixing of these states, where the value of δ is fixed to 2.15 MeV [47] and m ρ 0 and m ω 0 are the known masses [1]. This parameterisation is equivalent to that used in Ref. [48] given that ∆ 2 is small and therefore the term where it appears in the denominator can be neglected.

The Dalitz plot fit methodology
Given the large data sample, the large number of parameters used in the decay amplitude, and the need to normalise the total PDF at each iteration of the minimisation process, the GooFit framework [49] for maximum-likelihood fits is used. GooFit is based on GPU acceleration with parallel processing.
An unbinned maximum-likelihood fit to the data distribution in the Dalitz plot is performed. The likelihood function is written as a combination of the signal and background PDFs given by where f sig is the signal fraction and the product runs over the candidates in the final data sample. The background PDF, P bkg (s 12 , s 13 ), is the normalised background model, and is provided as the high definition histogram shown in Fig. 4. The normalised signal PDF, P sig (s 12 , s 13 ), is given by where (s 12 , s 13 ) is the efficiency model function included as the smoothed histogram shown in Figure 4; the denominator is the integral of the numerator over the Dalitz plot (DP) to guarantee that P sig (s 12 , s 13 ) is normalised at each iteration of the minimisation process. The fit parameters are the 50 pairs of magnitudes and phases of the S-wave amplitude, and the magnitudes and phases of the higher-spin components, except the ρ(770) 0 π + channel which is taken as the reference mode, with magnitude fixed to 1 and phase fixed to zero. The set of optimal parameters is determined by minimising the quantity −2 log L using the MINUIT [50] package. The fit fraction for the i th intermediate channel is defined as .
Due to interference, the sum of fit fractions can be less than or greater than 100%.
Interference fit fractions can also be defined, quantifying the level of interference between any pair i, j (i = j) of amplitude components, FF ij = DP 2 Re a i a j e i(δ i −δ j ) A i (s 12 , s 13 ) A * j (s 12 , s 13 ) ds 12 ds 13 DP | k a k e iδ k A k (s 12 , s 13 )| 2 ds 12 ds 13 .
By construction, the sum of fit fractions and interference terms is 100%. The fit quality is measured through the statistical quantity χ 2 defined as where the Dalitz plot is divided in N b bins and, for each bin, the number of observed candidates, N obs i , the number of candidates estimated from the fit model, N est i , and the uncertainty on their difference, σ i , are obtained. For unbinned maximum-likelihood fits, the number of degrees of freedom (ndof) range as [N b − q − 1, N b − 1] [51], where q is the number of free parameters, and it is used to calculate the corresponding range of χ 2 /ndof. This is done using the folded Dalitz plot -due to the symmetry of the D + → π − π + π + Dalitz plot with respect to the axis s 13 = s 12 , the variables s high and s low are defined, respectively, as the higher and the lower values of each pair (s 12 , s 13 ). The folded Dalitz plot is divided in N b = 625 bins using an adaptive binning algorithm, such that all bins have the same population. Besides the χ 2 /ndof, the value of −2 log L is also used to compare models. In addition, the distribution of residuals (N obs i − N est i )/σ i across the folded Dalitz plot is used for visual inspection of any local discrepancy between fit model and data, which are also compared through the projections of s high , s low , the sum of these projections, denoted s π − π + , and s 23 ≡ (p 2 + p 3 ) 2 .
As anticipated by visual inspection, there is a small but clear contribution from the resonant state ω(782), with an interference pattern with ρ(770) 0 , which can also be seen in the s π − π + projection in Fig. 5. This contribution was not observed in previous analyses due to their limited datasets, thus this effect is observed for the first time in this analysis. The alternative parameterisation of the ρ − ω mixing lineshape shown in Eq. 14 was tested, with |B| and φ B as free parameters. The outcome of the fit represents essentially the same solution: a difference in −2 log L of only +1 unit, no significant differences in all the other fit parameters, and the resulting ρ − ω amplitude being almost indistinguishable to that of the sum of ρ(770) 0 and ω(782) isobars. The values found for ρ − ω mixing parameters are |B| = 0.52 ± 0.02 ± 0.05 ± 0.01 and φ B = (158.8 ± 2.1 ± 2.6 ± 0.4) • , where the quoted uncertainties are, in order, statistical, experimental systematics and model systematic uncertainties.
While the contribution of the ρ(1450) 0 resonance has been reported previously in the D + → π − π + π + decay [15], this state alone is not enough to describe the π − π + P-wave amplitude at high mass. The inclusion of the ρ(1700) 0 resonance in the model improves the fit significantly, resulting in a difference in −2 log L of −488 units compared to that of the model without it. Its contribution is robust, with its stability in the fit being tested through many fit variations (as part of the systematic studies discussed in Sec. 8). The ρ(1700) 0 is a broad state (nominal width of 250 ± 100 MeV [1]), and its inclusion affects the whole Dalitz plot, with an interesting interference pattern with the other vector states, in particular with the ρ(1450) 0 resonance. The results presented here provide input for the debate of ρ(1450) 0 − ρ(1700) 0 interference, supporting the need for these two overlapping states to exist (see for instance the entry for ρ(1700) 0 under "Particle listings" in [1]). Note that if these states are excitations of the ρ 0 (770), their production should be also favoured in the D + decay.
The S-wave amplitude obtained from the final-model fit is shown in Fig. 6. The values of the magnitude and phases for each knot are shown in Table 4. The magnitude a 0 (m π − π + ) is large close to threshold and decreases until ∼ 0.9 GeV, with a steady phase increase, as expected from the dominant f 0 (500) contribution reported in previous analyses, and consistent with a dd source of the D + → π − π + π + weak decay. Starting at m π − π + ∼ 0.9 GeV the f 0 (980) signature is observed both as a sharp increase of the magnitude and a rapid variation (decrease and increase) in the phase, enhanced possibly by the opening of the KK channel. Starting at about m π − π + ∼ 1.4 GeV and peaking near m π − π + ∼ 1.5 GeV, another structure in the amplitude with a corresponding phase movement is observed, indicating the presence of at least one more scalar resonance, possibly f 0 (1370) or f 0 (1500), or a combination of the two.

Systematic uncertainties
The systematic uncertainties on the fit parameters are divided into two categories: first those coming from the impact of experimental aspects; and second, referred to as model  systematics, those corresponding to the uncertainties in resonance lineshape parameters such as masses, widths and form factors.
In the first category, effects due to the efficiency modelling, background contributions and intrinsic fit biases are considered. For the efficiency map, uncertainties arise from the finite size of the simulation samples, the effect of the binning scheme of the efficiency histogram prior to the 2D spline smoothing, and the procedure for obtaining PID efficiencies from the calibration samples. The first effect is estimated by generating a set of alternative efficiency histograms where the bin contents are varied according to a Poisson distribution, and performing the Dalitz fit with each alternative map. The root mean square (rms) of the distribution of each of the fit parameters is assigned as the corresponding systematic uncertainty. The systematic uncertainty due to the binning scheme of the efficiency map is assessed by performing the fits with smoothed efficiency maps obtained from varying the binning grid of the efficiency histogram from 15 × 15 (default) to 12 × 12 and 20 × 20. The largest variation in each fit parameter is assigned as the systematic uncertainty. The dominant uncertainties on the PID efficiency are due to the finite size of the calibration samples; the effect is assessed by varying the efficiencies from the calibration tables according to a Gaussian distribution centered at the nominal values and with widths equal to the statistical uncertainties, and obtaining new smoothed efficiency maps to refit the Dalitz plot and assess the impact on the fit parameters.
The effect of the signal-to-background ratio is estimated by varying the signal purity obtained from the invariant-mass fit within one standard deviation, repeating the Dalitz fit and assigning the difference in each parameter as the systematic uncertainty. The impact of the background model is addressed by determining it from either just the lower or just the higher mass sideband, and assigning the largest deviation in each fit parameter as the corresponding systematic uncertainty. The S-wave interval definition (knots in m π − π + ) is varied both in the number of knots (from 45 to 55) and position of the knots, by using an adaptive binning where the division is made to equalise the event population.
The resulting S-wave amplitude in each case after the fit is used to calculate the values of magnitude and phase in the original knot scheme and the rms in each fit parameter (including those from spin-1 and spin-2 states) is assigned as systematic uncertainty.
Finally, uncertainties due to biases from the fitting algorithm, including the ability to reproduce a given (continuous) S-wave amplitude, are assessed by generating a large number of pseudoexperiments given by the full PDF (signal and background). The PDF for the signal model has its parameters set to those obtained from a QMIPWA fit with an alternative 50-knot choice using the adaptive binning. The pseudoexperiments are fitted with the baseline 50-knot model, and the rms of the distribution of the difference in each parameter is assigned as the systematic uncertainty. Overall, the largest systematic uncertainties come from the finite size of the simulation samples, the efficiency map binning, and the fit bias. The systematic uncertainties are all added in quadrature and comprise the second uncertainty in Tables 2 and 4. Other studies are performed as crosschecks and lead to variations within statistical uncertainties: the use of the efficiency map without correcting for mismatches between data and simulation; splitting the sample in terms of magnet orientation during data taking; and studying the effect of mass resolution around the region of the narrow ω(782) state.
In the second category, isobar model systematics, two studies are performed. First, the mass and the width of all the spin-1 and spin-2 states are varied within their quoted uncertainties [1], and the rms of the difference in each fit parameter is assigned as a systematic uncertainty. The second study addresses the impact of changing the effective radii of both the resonances and the D + meson used in the Blatt-Weisskopf barrier factors, to 1.0 and 2.0 GeV −1 and to 4.0 and 6.0 GeV −1 , respectively, and the largest variation for each parameter is assigned as systematic uncertainty. These two effects are added in quadrature and comprise the third uncertainties presented in Tables 2 and 4.

Summary and conclusions
This paper presents a Dalitz plot analysis of the D + → π − π + π + decay. Using a sample containing more than six hundred thousand candidates, with a purity of 95%, the resonant structure of the decay is studied using the QMIPWA method, where the magnitude and the phase of the S-wave amplitude is obtained as a function of m π − π + , while the spin-1 and spin-2 contributions are included with an isobar model. This approach is motivated by the presence of broad and overlapping light scalar resonances below 2 GeV, with poorly known masses and widths.
The result of the Dalitz plot fit shows that the decay is dominated by the π − π + S-wave component corresponding to nearly 62% of the D + → π − π + π + decay rate, which is consistent with results from previous analyses. The P-wave amplitude is the second largest component, led by the ρ(770) 0 contribution at the level of 26% but including also the ω(782) resonance and the high-mass states ρ(1450) 0 and ρ(1700) 0 . The D-wave amplitude, consisting only of the f 2 (1270) state, accounts for about 14% of the D + → π − π + π + decay rate.
The contribution of D + → ω(782)π + in the D + → π − π + π + decay is observed for the first time, with a fit fraction of (0.103 ± 0.016)%. The ω(782) → π − π + decay is isospin violating and has been observed in association with the ρ(770) → π − π + in a few processes with very different interference patterns [44,48,[52][53][54][55][56][57]. Table 4: Fitted magnitude and phase of the S-wave amplitude at each π − π + mass knot, relative to the ρ(770) 0 π + channel. The uncertainties quoted are, in order, statistical, experimental systematics and model systematics. Interesting structures are observed in the S-wave amplitude shown in Fig. 6. The broad structure near threshold is associated with the f 0 (500) resonance in model-dependent analyses, the f 0 (980) state is clearly visible with an asymmetric peak lineshape, and a further peak with corresponding phase variation is observed around 1.5 GeV indicating the presence of a further high-mass scalar state.
While this analysis is based on the concept of a well-isolated companion pion, that is, a 2+1 approximation, the capacity of the QMIPWA approach to absorb some three-body final-state interaction effects should be considered in the interpretation of the resulting (2-body) π − π + S-wave amplitude obtained from the fit. This is the first time that the π − π + S-wave amplitude is extracted through a quasimodel-independent approach for the D + → π − π + π + decay, from threshold up to 1.7 GeV. Together with the companion analysis of the D + s → π − π + π + channel [20] this provide an important input to phenomenological studies.