Amplitude analysis of $B^{-} \to D^{+} \pi^{-} \pi^{-}$ decays

The Dalitz plot analysis technique is used to study the resonant substructures of $B^{-} \to D^{+} \pi^{-} \pi^{-}$ decays in a data sample corresponding to 3.0 ${\rm fb}^{-1}$ of $pp$ collision data recorded by the LHCb experiment during 2011 and 2012. A model-independent analysis of the angular moments demonstrates the presence of resonances with spins 1, 2 and 3 at high $D^{+}\pi^{-}$ mass. The data are fitted with an amplitude model composed of a quasi-model-independent function to describe the $D^{+}\pi^{-}$ S-wave together with virtual contributions from the $D^{*}(2007)^{0}$ and $B^{*0}$ states, and components corresponding to the $D^{*}_{2}(2460)^{0}$, $D^{*}_{1}(2680)^{0}$, $D^{*}_{3}(2760)^{0}$ and $D^{*}_{2}(3000)^{0}$ resonances. The masses and widths of these resonances are determined together with the branching fractions for their production in $B^{-} \to D^{+} \pi^{-} \pi^{-}$ decays. The $D^{+}\pi^{-}$ S-wave has phase motion consistent with that expected due to the presence of the $D^{*}_{0}(2400)^{0}$ state. These results constitute the first observations of the $D^{*}_{3}(2760)^{0}$ and $D^{*}_{2}(3000)^{0}$ resonances.

The Dalitz plot analysis technique is used to study the resonant substructures of B − → D þ π − π − decays in a data sample corresponding to 3.0 fb −1 of pp collision data recorded by the LHCb experiment during 2011 and 2012. A model-independent analysis of the angular moments demonstrates the presence of resonances with spins 1, 2 and 3 at high D þ π − mass. The data are fitted with an amplitude model composed of a quasi-model-independent function to describe the D þ π − S wave together with virtual contributions from the D Ã ð2007Þ 0 and B Ã0 states, and components corresponding to the D Ã 2 ð2460Þ 0 , D Ã 1 ð2680Þ 0 , D Ã 3 ð2760Þ 0 and D Ã 2 ð3000Þ 0 resonances. The masses and widths of these resonances are determined together with the branching fractions for their production in B − → D þ π − π − decays. The D þ π − S wave has phase motion consistent with that expected due to the presence of the D Ã 0 ð2400Þ 0 state. These results constitute the first observations of the D Ã 3 ð2760Þ 0 and D Ã 2 ð3000Þ 0 resonances, with significances of 10σ and 6.6σ, respectively. DOI: 10.1103/PhysRevD.94.072001

I. INTRODUCTION
There is strong theoretical and experimental interest in charm meson spectroscopy because it provides opportunities to study QCD predictions within the context of different models [1][2][3][4][5]. Experimental knowledge of the masses, widths and spins of the charged and neutral orbitally excited (1P) charm meson states has been gained through analyses of both prompt production [6,7] and three-body decays of B mesons [8][9][10][11][12][13]. Progress has been equally strong for excited charm-strange (cs) mesons [14][15][16][17][18]. These studies have in addition revealed several new states at higher masses, most of which have not yet been confirmed by analyses of independent data samples. Moreover, quantum numbers are only known for states studied in amplitude analyses of multibody B meson decays, since analyses of promptly produced excited charm states only determine whether the spin-parity is natural (i.e. J P ¼ 0 þ ; 1 − ; 2 þ ; …) or unnatural (i.e. J P ¼ 0 − ; 1 þ ; 2 − ; …), not the resonance spin. The experimental status of the neutral excited charm states is summarized in Table I (here and throughout the paper, natural units with ℏ ¼ c ¼ 1 are used). The D Ã 0 ð2400Þ 0 , D 1 ð2420Þ 0 , D 0 1 ð2430Þ 0 and D Ã 2 ð2460Þ 0 mesons are generally understood to be the four 1P states. The spectroscopic identification for heavier states is not clear.
The B − → D þ π − π − decay mode has been previously studied in Refs. [8,9]. The inclusion of charge-conjugate processes is implied throughout the paper. The Dalitz plot (DP) models that were used contained components for two excited charm states, the D Ã 0 ð2400Þ 0 and D Ã 2 ð2460Þ 0 resonances, together with nonresonant amplitudes. More recently, a DP analysis of B − → D þ K − π − decays [12] included, in addition, a contribution from the D Ã 1 ð2760Þ 0 state. The properties of this state indicate that it belongs to the 1D family [20,21]. The D Ã 1 ð2760Þ 0 width is found to be larger than in previous measurements based on prompt production, which may be due to a contribution from an additional resonance, as would be expected if both 2S and 1D states with spin-parity J P ¼ 1 − are present in this TABLE I. Measured properties of neutral excited charm states. World averages are given for the 1P resonances (top part), while all measurements are listed for the heavier states (bottom part). Where two uncertainties are given, the first is statistical and second systematic; where a third is given, it is due to model uncertainty. The uncertainties on the averages for the D Ã 0 ð2400Þ 0 mass and the D 1 ð2420Þ 0 and D Ã 2 ð2460Þ 0 masses and widths are inflated by scale factors to account for inconsistencies between measurements. The quoted D Ã 2 ð2460Þ 0 averages do not include the recent result from Ref. [12].

Resonance
Mass (MeV) Width (MeV) J P Ref.  [12] region. There should also be a 1D state with J P ¼ 3 − at similar mass, as seen in the charm-strange system [15,16].
As yet there is no evidence for such a neutral charm state, but a DP analysis ofB 0 → D 0 π þ π − decays [11] led to the first observation of the D Ã 3 ð2760Þ þ state. One challenge for DP analyses with large data samples is the modeling of broad resonances that interfere with nonresonant amplitudes in the same partial wave. Inclusion of both contributions in an amplitude fit can violate unitarity in the decay matrix element, and also gives results that are difficult to interpret due to large interference effects. In the case of B − → D þ π − π − decays this is particularly relevant for the D þ π − S wave, where both the D Ã 0 ð2400Þ 0 resonance and a nonresonant contribution are expected. In the π þ π − and K þ π − systems such effects can be handled with a K-matrix approach or specific models such as the LASS function [22] inspired by low-energy scattering data, respectively. In the absence of any D þ π − scattering data, a viable alternative approach is to use a quasi-modelindependent description, in which the partial wave is fitted using splines to describe the magnitude and phase as a function of mðD þ π − Þ. Determination of the phase depends on interference of the S wave with another partial wave, so that some model dependence remains due to the description of the other amplitudes in the decay. This approach was first applied to the Kπ S wave using D þ → K − π þ π þ decays [23]. Subsequent uses include further studies of the Kπ S wave [24][25][26][27] as well as the K þ K − [28] and π þ π − [29] S waves, in various processes. Similar methods have been used to determine the phase motion of exotic hadron candidates [30,31]. Quasi-model-independent information on the D þ π − S wave could be used to develop better models of the dynamics in the D þ π − system [32][33][34][35].
In this paper, the DP analysis technique is employed to study the contributing amplitudes in B − → D þ π − π − decays, where the charm meson is reconstructed through D þ → K − π þ π þ decays. The analysis is based on a data sample corresponding to an integrated luminosity of 3.0 fb −1 of data collected with the LHCb detector during 2011 when the pp collision centerof-mass energy was ffiffi ffi s p ¼ 7 TeV, and 2012 with ffiffi ffi The paper is organized as follows. Section II provides a brief description of the LHCb detector and the event reconstruction and simulation software. The selection of signal candidates is described in Sec. III and the determination of signal and background yields is presented in Sec. IV. The angular moments of B − → D þ π − π − decays are studied in Sec. V and are used to guide the amplitude analysis. The DP analysis formalism is reviewed briefly in Sec. VI, and implementation of the amplitude fit is given in Sec. VII. Experimental and modeldependent systematic uncertainties are evaluated in Sec. VIII, and the results and a summary are presented in Sec. IX.

II. LHCb DETECTOR
The LHCb detector [36,37] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector 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 polarity of the dipole magnet is reversed periodically throughout data taking. The tracking system provides a measurement of momentum, p, of charged particles with 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 vertex, the impact parameter, 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. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
The trigger consists of a hardware stage based on information from the calorimeter and muon systems followed by a software stage, in which all tracks with p T > 500ð300Þ MeV are reconstructed for data collected in 2011 (2012). The software trigger line used in the analysis reported in this paper requires a two-, threeor four-track secondary vertex with significant displacement from the primary pp interaction vertices (PVs). At least one charged particle must have p T > 1.7 GeV and be inconsistent with originating from the PV. A multivariate algorithm [38] is used for the identification of secondary vertices consistent with the decay of a b hadron.
In the off-line selection, the objects that fired the trigger are associated with reconstructed particles. Selection requirements can therefore be made not only on the trigger line that fired, but on whether the decision was due to the signal candidate, other particles produced in the pp collision, or a combination of both. Signal candidates are accepted off-line if one of the final state particles created a cluster in the hadronic calorimeter with sufficient transverse energy to fire the hardware trigger.
Simulated events are used to characterize the detector response to signal and certain types of background events. In the simulation, pp collisions are generated using PYTHIA [39] with a specific LHCb configuration [40]. Decays of hadronic particles are described by EVTGEN [41], in which final state radiation is generated using PHOTOS [42]. The interaction of the generated particles with the detector and its response are implemented using the GEANT4 toolkit [43] as described in Ref. [44].

III. SELECTION REQUIREMENTS
The selection criteria are the same as those used in Ref. [12], where a detailed description is given, with the exception that only candidates that are triggered by at least one of the signal tracks are retained in order to minimize the uncertainty on the efficiency. First, loose requirements are applied in order to obtain a visible peak in the B candidate invariant mass distribution. These criteria are found to be 91% efficient on simulated signal decays. The remaining data are then used to train two artificial neural networks [45] that separate signal from different categories of background. The first is designed to distinguish candidates that contain real D þ → K − π þ π þ decays from those that do not; the second separates signal B − → D þ π − π − decays from background combinations. The SPLOT technique [46] is used to statistically separate signal decays from background combinations using the D (B) candidate mass as the discriminating variable for the first (second) network. The first network takes as input properties of the D candidate and its decay product tracks, including information about kinematics, track and vertex quality. The second uses a total of 27 input variables, including the output of the first network, as described in Ref. [12]. The neural network input quantities depend only weakly on the position in the DP, so that training the networks with the same data sample used for the analysis does not bias the results. A requirement that reduces the combinatorial background by an order of magnitude, while retaining about 75% of the signal, is imposed on the second neural network output.
Particle identification (PID) requirements are applied to all five final state tracks to select pions or kaons as necessary. Background from D þ s → K − K þ π þ decays, where the K þ is misidentified as a π þ meson, are suppressed using a tight PID criterion on the higher momentum π þ from the D þ decay. The combined efficiency of the PID requirements on the five final state tracks is determined using D Ãþ → D 0 π þ , D 0 → K − π þ calibration data [47] and found to be around 70%.
Potential background from Λ þ c → pK − π þ decays, misreconstructed as D þ candidates, is removed if the invariant mass lies in the range 2280-2300 MeV when the proton mass hypothesis is applied to the low momentum pion track. Decays of B − mesons to the K − π þ π þ π − π − final state that do not proceed via an intermediate charm state are removed by requiring that the D and B candidate decay vertices are separated by at least 1 mm. The signal efficiency of this requirement is approximately 85%.
To improve mass resolution, the momenta of the final state tracks are rescaled [48,49] using weights obtained from a sample of J=ψ → μ þ μ − decays where the measured mass peak is matched to the known value [19].
Additionally, a kinematic fit [50] is performed to candidates in which the invariant mass of the D decay products is constrained to equal the world average D mass [19]. A B mass constraint is added in the calculation of the variables that are used in the DP fit.
Candidate B mesons with invariant mass in the range 5100-5800 MeV are retained for further analysis. Following all selection requirements, multiple candidates are found in approximately 0.4% of events. All candidates are retained and treated in the same way.

IV. DETERMINATION OF SIGNAL AND BACKGROUND YIELDS
The signal and background yields are measured using an extended unbinned maximum likelihood fit to the D þ π − π − invariant mass distribution. The candidates are comprised of true signal decays and several sources of background. Partially reconstructed backgrounds come from b hadron decays where one or more final state particles are not reconstructed. Combinatorial background originates from random combinations of tracks, potentially including a real D þ → K − π þ π þ decay. Misidentified background arises from b hadron decays in which one of the final state particles is not correctly identified. Potential residual background from charmless B decays is reduced to a negligible level by the requirement that the flight distance of the D candidate be greater than 1 mm.
Signal candidates are modeled by the sum of two Crystal Ball (CB) functions [51] with a common peak position of the Gaussian core and tails on opposite sides. The relative normalization of the narrower CB shape and the ratio of widths of the CB functions are constrained, by including a Gaussian penalty term in the likelihood, to the values found in fits to simulated samples. The tail parameters of the CB shapes are fixed to those found in simulation.
The main source of partially reconstructed background is the B − → D Ãþ π − π − channel with subsequent D Ãþ → D þ γ or D Ãþ → D þ π 0 decay, where the neutral particle is not reconstructed. A nonparametric shape derived from simulation is used to model this contribution. The shape is characterized by an edge around 100 MeV below the B peak, where the exact position of the edge depends on properties of the decay, including the D Ãþ polarization. As in previous studies of similar processes [12,52], the fit quality improves when the shape is allowed to be offset by a small shift (≈3.5 MeV) that is determined from the data.
The combinatorial background is modeled with a linear function, where the slope is free to vary. Many sources of misidentified background have broad D þ π − π − invariant mass distributions that can be absorbed into the combinatorial background component. The exceptions are B − → D ðÃÞþ K − π − decays that produce distinctive shapes in the B candidate invariant mass distribution. These backgrounds are combined into a single nonparametric shape determined from simulated samples that are weighted to account for the known DP distribution for B − → D þ K − π − decays [12]. The ratio of D þ and D Ãþ components in the There are ten parameters in the fit that are free to vary: the yields for signal and combinatorial B − → D ðÃÞþ K − π − and B − → D Ãþ π − π − backgrounds, the combinatorial background slope, the shared mean of the double CB shape, the width and relative normalization of the narrower CB and the ratio of CB widths, and the shift parameter of the B − → D Ãþ π − π − shape. The result of the fit is shown in Fig. 1 and gives a signal yield of approximately 29 000 decays. The χ 2 per degree of freedom for this projection of the fit is 1.16 calculated with statistical uncertainties only. Component yields are shown in Table II for both the full fit range and the signal region defined as AE2.5σ around the B peak, where σ is the width parameter of the dominant CB function in the signal shape; this corresponds to 5235.3 < mðD þ π − π − Þ < 5320.8 MeV.
A Dalitz plot [53] is a two-dimensional representation of the phase space for a three-body decay in terms of two of the three possible two-body invariant mass squared combinations. In B − → D þ π − π − decays there are two indistinguishable pions in the final state, so the two m 2 ðD þ π − Þ combinations are ordered by value and the DP axes are defined as m 2 ðD þ π − Þ min and m 2 ðD þ π − Þ max . The ordering causes a "folding" of the DP from the minimum value of The DP distribution of the candidates in the signal region that are used in the DP fit is shown in Fig. 2 (left). The same data are shown in the square Dalitz plot (SDP) in Fig. 2 (right). The SDP is defined by the variables m 0 and θ 0 , which are given by where m max π − π − ¼ m B − − m D þ and m min π − π − ¼ 2m π − are the kinematic boundaries of mðπ − π − Þ and θðπ − π − Þ is the helicity angle of the π − π − system (the angle between the momenta of the D meson and one of the pions, evaluated in the π − π − rest frame). With m 0 and θ 0 defined in terms of the π − π − mass and helicity angle in this way, only the region of the SDP with θ 0 ≤ 0.5 is populated due to the symmetry of the two pions in the final state. The SDP is used to describe the signal efficiency variation and distribution of background candidates, as described in Sec. VII.

V. STUDY OF ANGULAR MOMENTS
The angular moments of the B − → D þ π − π − decays are studied to investigate which amplitudes to include in the DP fit model. Angular moments are determined by weighting the data by the Legendre polynomial P L ðcos θðD þ π − ÞÞ, where θðD þ π − Þ is the helicity angle of the D þ π − system, i.e. the angle between the momenta of the pion in the D þ π − system and the other pion from the B − decay, evaluated in the D þ π − rest frame. The moment hP L i is the sum of the weighted data in a bin of D þ π − mass with background contributions subtracted using sideband data and efficiency corrections, determined as in Sec. VII A, applied. Each of TABLE II. Yields of the various components in the fit to B − → D þ π − π − candidate invariant mass distribution. Note that the yields in the signal region are scaled from the full mass range.

Component
Full mass range Signal region the moments contains contributions from certain partial waves and interference terms. For the S-, P-, D-and F-wave amplitudes denoted by h j e iδ j (j ¼ 0, 1, 2, 3 respectively), These expressions assume that there are no contributions from partial waves higher than F wave. Thus, they are valid only in regions of the DP unaffected by the folding, i.e. for mðD þ π − Þ ≲ 3.2 GeV, where the full range of the D þ π − helicity angle distribution is available. Above this mass, the orthogonality of the Legendre polynomials does not hold and a straightforward interpretation of the angular moments in terms of the contributing partial waves is not possible. Nevertheless, the angular moments provide a useful way to judge the agreement of the fit result with the data, complementary to the projections onto the invariant masses.
The unnormalized angular moments hP 0 i-hP 6 i are shown in Fig. 3 for the D þ π − invariant mass range 2.0-4.0 GeV. The D Ã 2 ð2460Þ 0 resonance is clearly seen in the hP 4 i distribution of Fig. 3(e). From Eqs. (3) and (5) it can be inferred that the structures in the distributions of hP 1 i and hP 3 i below 3 GeV suggest that there is interference both between the S-and P-wave amplitudes and between the P-and D-wave amplitudes. Therefore broad spin 0 and spin 1 components are required in the DP model. In addition, structure in hP 2 i around 2.76 GeV implies the possible presence of a spin 1 resonance in that region. The angular moments hP 7 i and hP 8 i shown in Fig. 4, show no structure, consistent with the assumption that contributions from higher partial waves and from the isospin-2 dipion channel are small.
Zoomed views of the fourth and sixth moments in the region around mðD þ π − Þ ¼ 3 GeV are shown in Fig. 5. A wide bump is visible in the distribution of hP 4 i at mðD þ π − Þ ≈ 3 GeV. Although close to the point where the DP folding affects the interpretation of the moments, this enhancement suggests that an additional spin 2 resonance could be contributing in this region. A peak is also seen at mðD þ π − Þ ≈ 2.76 GeV in the hP 6 i distribution, suggesting that a spin 3 resonance should be included in the DP model. As discussed in Sec. I, other recent analyses [6,7,11,12,15,16] suggest that both spin 1 and spin 3 states could be expected in this region.

VI. DALITZ PLOT ANALYSIS FORMALISM
The isobar approach [54][55][56] is used to describe the complex decay amplitude as the coherent sum of amplitudes for intermediate resonant and nonresonant decays. The total amplitude is given by where the complex coefficients c j describe the relative contribution of each intermediate process.
Here, and for the remainder of this section, m 2 ðD þ π − Þ min and m 2 ðD þ π − Þ max are referred to as s and t, respectively. The resonant dynamics are encoded in the F j ðs; tÞ terms, each of which is normalized such that the integral of the magnitude squared across the DP is unity. The amplitude is explicitly symmetrized to take account of the Bose symmetry of the final state due to the identical pions, i.e.
Aðs; tÞ↦Aðs; tÞ þ Aðt; sÞ: This substitution is implied throughout this section. For a D þ π − resonance wherep andq are the momenta calculated in the D þ π − rest frame of the particle not involved in the resonance and one of the resonance decay products, respectively. The functions X, T and R are described below. The XðzÞ terms are Blatt-Weisskopf barrier factors [57], where z ¼ jqjr BW or jpjr BW and r BW is the barrier radius, and are given by where L is the spin of the resonance and z 0 is defined as the value of z where the invariant mass is equal to the mass of the resonance. Since the B − meson has zero spin, L is also the orbital angular momentum between the resonance and the other pion. The barrier radius r BW is taken to be 4.0 GeV −1 ≈ 0.8 fm [16,58] for all resonances.
The Tðp;qÞ functions describe the angular distribution and are given in the Zemach tensor formalism ] 59,60 ], L ¼ 0∶ Tðp;qÞ ¼ 1; These are proportional to the Legendre polynomials, P L ðxÞ, where x is the cosine of the helicity angle betweeñ p andq. The function RðsÞ of Eq. (11) describes the resonance line shape. Resonant contributions to the total amplitude are modeled by relativistic Breit-Wigner (RBW) functions given by with a mass-dependent decay width defined as where q 0 is the value of q ≡ jqj when m ¼ m 0 and Γ 0 is the full width. Virtual contributions, from resonances with pole masses outside the kinematically allowed region, can be described by RBW functions with one modification: the pole mass m 0 is replaced with an effective mass, m eff 0 , in the allowed region of s, when the parameter q 0 is calculated. The term m eff 0 is given by the ad hoc formula [16] m eff where m max and m min are the upper and lower thresholds of s. Note that m eff 0 is only used in the calculation of q 0 , so only the tail of such virtual contributions enters the DP.
A quasi-model-independent approach is used to describe the entire D þ π − spin 0 partial wave. The total D þ π − S wave is fitted using cubic splines to describe the magnitude and phase variation of the spin 0 amplitude. Knots are defined at fixed values of mðD þ π − Þ and splines give a smooth interpolation of the magnitude and phase of the S wave between these points. The S-wave magnitude and phase are both fixed to zero at the highest mass knot in order to ensure sensible behavior at the kinematic limit. For the knot at mðD þ π − Þ ¼ 2.4 GeV, close to the peak of the D Ã 0 ð2400Þ 0 resonance, the magnitude and phase values are fixed to 0.5 and 0, respectively, as a reference. The magnitude and phase values at every other knot position are determined from the fit.
The folding of the Dalitz plot has implications for the choice of knot positions. Since the S-wave amplitude varies with mðD þ π − Þ, its reflection onto the other DP axis gives a helicity angle distribution that corresponds to higher partial waves. Equally, if knots are included at high mðD þ π − Þ, the quasi-model-independent D þ π − S-wave amplitude can absorb resonant contributions with nonzero spin due to their reflections. To avoid this problem, only a single knot with floated parameters is used above the minimum value of m 2 ðD þ π − Þ max , specifically at 4.1 GeV (as mentioned above, the amplitude is fixed to zero at the highest mass knot at 5.1 GeV). At lower mðD þ π − Þ, knots are spaced every 0.1 GeV from 2.0 GeV up to 3.1 GeV, except that the knot at 3.0 GeV is removed in order to stabilize the fit.
Neglecting reconstruction effects, the DP probability density function would be P phys ðs; tÞ ¼ jAðs; tÞj 2 R R DP jAðs; tÞj 2 dsdt : ð17Þ The effects of nonuniform signal efficiency and of background contributions are accounted for as described in Sec. VII. The probability density function depends on the complex coefficients introduced in Eq. (9), as well as the masses and widths of the resonant contributions and the parameters describing the D þ π − S wave. These parameters are allowed to vary freely in the fit. Results for the complex coefficients are dependent on the amplitude formalism, normalization and phase convention, and consequently may be difficult to compare between different analyses. It is therefore useful to define fit fractions and interference fit fractions to provide convention-independent results. Fit fractions are defined as the integral over the DP for a single contributing amplitude squared divided by that of the total amplitude squared, The sum of the fit fractions is not required to be unity due to the potential presence of net constructive or destructive interference. Interference fit fractions are defined, for i < j only, as VII. DALITZ PLOT FIT

A. Signal efficiency
Variation of the efficiency across the phase space of B − → D þ π − π − decays is studied in terms of the SDP, since the efficiency variation is typically greatest close to the kinematic boundaries of the conventional DP. The causes of efficiency variation across the SDP are the detector acceptance and trigger, selection and PID requirements. Simulated samples generated uniformly over the SDP are used to evaluate the efficiency variation. Data-driven corrections are applied to correct the simulation for known discrepancies with the data, for the tracking, trigger and PID efficiencies, using identical methods to those described in Ref. [16]. The efficiency distributions are fitted with two-dimensional cubic splines to smooth out statistical fluctuations due to limited sample size. Figure 6 shows the efficiency variation over the SDP.

B. Background studies
The yields presented in Table II show that the important background components in the signal region are from combinatorial background and B − → D ðÃÞþ K − π − decays. The SDP distribution of B − → D ðÃÞþ K − π − decays is obtained from simulated samples using the same procedures as described in Sec. IV to apply weights and combine the D þ and D Ãþ contributions. The distribution of combinatorial background events is obtained from D þ π − π − candidates in the high-mass sideband defined to be 5500-5800 MeV. Figure 7 shows the SDP distributions of these backgrounds, which are used in the Dalitz plot fit.
The DP fit is performed using the LAURA++ [61] package, and the likelihood function is given by where the index i runs over n c candidates, while k sums over the probability density functions P k with a yield of N k candidates in each component. For signal events P k ≡ P sig is similar to Eq. (17), but is modified such that the jAðs; tÞj 2 terms are multiplied by the efficiency function described in Sec. VII A. The mass resolution is approximately 2.4 MeV, which is much less than the width of the narrowest contribution to the Dalitz plot (∼50 MeV); therefore, this has negligible effect on the likelihood. Its effect on the measurement of masses and widths of resonances is, however, considered as a systematic uncertainty.
Using the results of the moments analysis presented in Sec. V as a guide, a B − → D þ π − π − DP model is constructed by including various resonant, nonresonant and virtual amplitudes. Only intermediate states with natural spin-parity are included because unnatural spin-parity states do not decay to two pseudoscalars. Amplitudes that do not contribute significantly and cause the fit to become unstable are discarded. Alternative and additional contributions that have been considered include an isobar description of the D þ π − S wave including the D Ã 0 ð2400Þ 0 resonance and a nonresonant amplitude, a nonresonant P-wave component, an isospin-2 ππ interaction described by a unitary model as in Refs. [24,62] (see also Refs. [63][64][65]), and quasi-model-independent descriptions of partial waves other than the D þ π − S wave.
The resulting baseline signal model consists of the seven components listed in Table III: four resonances, two virtual resonances and a quasi-model-independent description of the D þ π − S wave. There are 42 free parameters in this model. The broad P-wave structure indicated by the angular moments is adequately described by the virtual D Ã ð2007Þ 0   Table I), while the parameters of the D Ã 3 ð2760Þ 0 resonance seem to be consistent with earlier measurements. An excess at mðD þ π − Þ ≈ 3000 MeV was reported in Ref. [7], but the parameters of this state were not reported with systematic uncertainties. The baseline model provides a better quality fit than the alternative models that are discussed in Sec. VIII. The inclusion of all components of the model is necessary to obtain a good description of the data, as described in Sec. IX.
The real and imaginary parts of the complex coefficients for each of the components are free parameters of the fit, except for the D Ã 2 ð2460Þ 0 contribution that is taken to be a reference amplitude with real and imaginary parts of its complex coefficient c k fixed to 1 and 0, respectively. Parameters such as magnitudes and phases for each amplitude, the fit fractions and interference fit fractions are calculated from these quantities. The statistical uncertainties are determined using large samples of pseudoexperiments to ensure that correlations between parameters are accounted for.

D. Dalitz plot fit results
The masses and widths of the D Ã 2 ð2460Þ 0 , D Ã 1 ð2680Þ 0 , D Ã 3 ð2760Þ 0 and D Ã 2 ð3000Þ 0 resonances are determined from the fit and are given in Table IV. The floated complex coefficients at each knot position and the splines describing the total D þ π − S wave are shown in Fig. 8. The phase motion at low mðD þ π − Þ is consistent with that expected due to the presence of the D Ã 0 ð2400Þ 0 state. There is, however, an ambiguous solution with the opposite phase motion in this region, which occurs since there are significant contributions only from S and P waves and thus only cosðδ 0 − δ 1 Þ can be determined as seen in Eq. (3). Since the P wave in this region is described by the D Ã v ð2007Þ 0 amplitude, and hence has slowly varying phase, the entire D þ π − S wave has a sign ambiguity. Similar ambiguities have been observed previously [23]. Only results consistent with the expected phase motion are reported. Table V shows the values of the complex coefficients and fit fractions for each amplitude. The interference fit fractions are given in the Appendix.
Given the complexity of the DP fit, the minimization procedure may find local minima in the likelihood function. To try to ensure that the global minimum is found, the fit is performed many times with randomized initial values for the c j terms. No other minima are found with negative  [19]. States labeled with subscript v are virtual contributions. The model "MIPW" refers to the quasi-model-independent partial wave approach.

Resonance
Spin Model Parameters D Ã 2 ð2460Þ 0 2 RBW Determined from data (see Table IV)  log-likelihood values close to that of the global minimum so they are not considered further. The consistency of the fit model and the data is evaluated in several ways. Numerous one-dimensional projections comparing the data and fit model (including several shown below and those from the moments study in Sec. V) show good agreement. Additionally, a two-dimensional χ 2 value is calculated by comparing the data and the fit model distributions across the SDP in 484 equally populated bins. Figure 9 shows the normalized residual in each bin. The distribution of the z-axis values from Fig. 9 is consistent with a unit Gaussian centered on zero. Further checks using unbinned fit quality tests [66] show satisfactory agreement between the data and the fit model.
One-dimensional projections of the baseline fit model and data onto mðD þ π − Þ min , mðD þ π − Þ max and mðπ − π − Þ are shown in Fig. 10. The model is seen to give a good description of the data sample, with the most evident discrepancy at low values of mðD þ π − Þ max , a region of the DP [that corresponds to high values of mðπ − π − Þ and mðD þ π − Þ min ≈ 3.2 GeV] in which many different amplitudes contribute. In Fig. 11, zoomed views of the mðD þ π − Þ min invariant mass projection are provided for regions at threshold and around the D Ã 2 ð2460Þ 0 , D Ã 1 ð2680Þ 0 -D Ã 3 ð2760Þ 0 and D Ã 2 ð3000Þ 0 resonances. Projections of the cosine of the D þ π − helicity angle in the same regions of mðD þ π − Þ min are also shown in Fig. 11. Good agreement is seen in all these projections, suggesting that the model gives an acceptable description of the data and the spin assignments of the D Ã 1 ð2680Þ 0 , D Ã 3 ð2760Þ 0 and D Ã 2 ð3000Þ 0 states are correct.

VIII. SYSTEMATIC UNCERTAINTIES
Sources of systematic uncertainty are divided into two categories: experimental and model uncertainties. The sources of experimental systematic uncertainty are the signal and background yields in the signal region, the SDP distributions of the background components, the efficiency variation across the SDP, and possible fit bias. Model uncertainties arise due to the fixed parameters in the amplitude model, the addition of amplitudes not included in the baseline fit, the modeling of the amplitudes from virtual resonances, and the effect of removing the least wellmodeled part of the phase space. The systematic uncertainties from each source are combined in quadrature.
The signal and background yields in the signal region are determined from the fit to the B candidate invariant mass distribution, as described in Sec. IV. The total uncertainty on each yield, including systematic effects due to the modeling of the components in the B candidate mass fit, is calculated, and the yields varied accordingly in the DP fit. The deviations from the baseline DP fit result are assigned as systematic uncertainties.
The effect of imperfect knowledge of the background distributions over the SDP is tested by varying the bin contents of the histograms used to model the shapes within their statistical uncertainties. For B − → D ðÃÞþ K − π − decays the ratio of the D Ãþ and D þ contributions is varied. Where applicable, the reweighting of the SDP distribution of the simulated samples is removed. Changes in the results compared to the baseline DP fit result are again assigned as systematic uncertainties.
The uncertainty related to the knowledge of the variation of efficiency across the SDP is determined by varying the  efficiency histograms before the spline fit is performed. The central bin in each 3 × 3 cluster is varied by its statistical uncertainty and the surrounding bins in the cluster are varied by interpolation. This procedure accounts for possible correlations between the bins, since a systematic effect on a given bin is likely also to affect neighboring bins. An ensemble of DP fits is performed, each with a unique efficiency histogram, and the effects on the results are assigned as systematic uncertainties. An additional systematic uncertainty is assigned by varying the binning scheme of the control sample used to determine the PID efficiencies. Systematic uncertainties related to possible intrinsic fit bias are investigated using an ensemble of pseudoexperiments. Differences between the input and fitted values from the ensemble for the fit parameters are found to be small. Systematic uncertainties are assigned as the sum in quadrature of the difference between the input and output values and the uncertainty on the mean of the output value determined from a fit to the ensemble. The only fixed parameter in the line shapes of resonant amplitudes is the Blatt-Weisskopf barrier radius, r BW . To account for potential systematic effects, this is varied between 3 and 5 GeV −1 [16], and the difference compared to the baseline fit model is assigned as an uncertainty. The choice of knot positions in the quasi-model-independent description of the D þ π − S wave is another source of possible systematic uncertainty. This is evaluated from the change in the fit results when more knots are added at low mðD þ π − Þ. As discussed in Sec. VI, it is not possible to add more knots at high mðD þ π − Þ without destabilizing the fit. As discussed in Sec. I, it is possible that there is more than one spin 1 resonance in the range 2.6 < mðD þ π − Þ < 2.8 GeV. The measured parameters of the D Ã 1 ð2680Þ 0 resonance are most consistent with those given for the D Ã ð2650Þ state in Table I; therefore the effect of including an additional D Ã ð2760Þ contribution is considered as a source of systematic uncertainty. Separate fits are performed with the parameters of the D Ã ð2760Þ state fixed to the values determined by BABAR [6] and LHCb [7] and the larger of the deviations from the baseline results is taken as the associated uncertainty. Additional fits are performed with the value of the D Ã v ð2007Þ 0 width given in Table III, which corresponds to the current experimental upper limit [19] replaced by the measured central value for the D Ã ð2010Þ þ (83.4 keV); the associated systematic uncertainty is negligible. The dependence of the results on the effective pole mass description of Eq. (16) that is used for the virtual resonance contributions is found by using a fixed width in Eq. (14), removing the dependence on m eff 0 . A discrepancy between the model and the data is seen in the low mðD þ π − Þ max region, as discussed in Sec. VII D. Since this may not be accounted for by the other sources of systematic uncertainty, the effect on the results is determined by performing fits where this region of the DP is vetoed by removing separately candidates with either mðD þ π − Þ max < 3.3 GeV or mðπ − π − Þ > 3.05 GeV. Systematic uncertainties are assigned as the difference in the fitted parameters compared to the baseline fit.
Contributions to the experimental and model systematic uncertainties for the fit fractions, masses and widths are 0 resonance, and some changing when the poorly modeled region of phase space is vetoed. The effect of the finite mass resolution described in Sec. VII C on the measurements of the masses and widths of resonances is found to be negligible.
Several cross-checks are performed to confirm the stability of the results. The data sample is divided into two parts depending on the charge of the B candidate, the polarity of the magnet and the year of data taking. All fits give consistent results.

IX. RESULTS AND SUMMARY
Results for the complex coefficients multiplying each amplitude are reported in Table VIII, and those that describe the D þ π − S wave amplitude are shown in Table IX. These complex numbers are reported in terms of real and imaginary parts and also in terms of magnitude and phase as, due to correlations, the propagation of uncertainties from one form to the other may not be trivial. Results for the interference fit fractions are given in the Appendix.
The fit fractions summarized in Table X for resonant contributions are converted into quasi-two-body product branching fractions by multiplying by the B − → D þ π − π − branching fraction. This value is taken from the world average after a correction for the relative branching fractions of B þ B − and B 0B0 pairs at the ϒð4SÞ resonance, Γðϒð4SÞ → B þ B − Þ=Γðϒð4SÞ → B 0B0 Þ ¼ 1.055 AE 0.025 [19], giving BðB − → D þ π − π − Þ ¼ ð1.014 AE 0.054Þ × 10 −3 . The product branching fractions are shown in Table XI; they cannot be converted into absolute branching fractions because the branching fractions for the resonance decays to D þ π − are unknown.
The masses and widths of the where the three quoted errors are statistical, experimental systematic and model uncertainties. The results for the D Ã 2 ð2460Þ 0 are consistent with the PDG averages [19] given in Table I. The D Ã 1 ð2680Þ 0 state has parameters close to those measured for the D Ã ð2650Þ resonance observed by LHCb in prompt production in pp collisions [7]. As discussed in Sec. I, both 2S and 1D states with spin-parity J P ¼ 1 − are expected in this region. Similarly, the D Ã 3 ð2760Þ 0 state has parameters close to those for the D Ã ð2760Þ states reported in Refs. [6,7] and for the charged D Ã 3 ð2760Þ þ state [11]. It appears likely to be a member of the 1D family. The D Ã 2 ð3000Þ 0 state has parameters that are not consistent with any previously observed resonance, although due to the large uncertainties it cannot be ruled out that it has a common origin with the D Ã ð3000Þ state that was reported, without evaluation of systematic uncertainties, in Ref. [7]. It could potentially be a member of the 2P or 1F family.
Removal of any of the D Ã 1 ð2680Þ 0 , D Ã 3 ð2760Þ 0 and D Ã 2 ð3000Þ 0 states from the baseline fit model results in large changes of the likelihood value. To investigate the effect of the systematic uncertainties, a similar likelihood ratio test is performed in the alternative models that give the largest uncertainties on the parameters of these resonances. Accounting for the 4 degrees of freedom associated with each resonance, the significances of the D Ã 1 ð2680Þ 0 and D Ã 3 ð2760Þ 0 states including systematic uncertainties are found to be above 10σ, while that for the D Ã 2 ð3000Þ 0 state is 6.6σ. Assigning alternative spin  hypotheses to these states results in similarly large changes in likelihood.
In summary, an analysis of the amplitudes contributing to B − → D þ π − π − decays has been performed using a data sample corresponding to 3.0 fb −1 of pp collision data recorded by the LHCb experiment. The Dalitz plot fit model containing resonant contributions from the D Ã 2 ð2460Þ 0 , D Ã 1 ð2680Þ 0 , D Ã 3 ð2760Þ 0 and D Ã 2 ð3000Þ 0 states, virtual D Ã v ð2007Þ 0 and B Ã0 v resonances and a quasi-modelindependent description of the full D þ π − S wave has been found to give a good description of the data. These results constitute the first observations of the D Ã 3 ð2760Þ 0 and D Ã 2 ð3000Þ 0 resonances and may be useful to develop improved models of the dynamics in the D þ π − system.

ACKNOWLEDGMENTS
We express our gratitude to our colleagues in the CERN accelerator departments for the excellent performance of the LHC. We thank the technical and administrative staff at the LHCb institutes. We acknowledge support from CERN and from the national agencies:

APPENDIX: RESULTS FOR INTERFERENCE FIT FRACTIONS
The central values and statistical errors for the interference fit fractions are shown in Table XII. The experimental systematic and model uncertainties are given in Table XIII.     0.08