Determination of quantum numbers for several excited charmed mesons observed in $B^- \to D^{*+} \pi^- \pi^-$ decays

A four-body amplitude analysis of the $B^- \to D^{*+} \pi^- \pi^-$ decay is performed, where fractions and relative phases of the various resonances contributing to the decay are measured. Several quasi-model-independent analyses are performed aimed at searching for the presence of new states and establishing the quantum numbers of previously observed charmed meson resonances. In particular the resonance parameters and quantum numbers are determined for the $D_1(2420)$, $D_1(2430)$, $D_0(2550)$, $D^*_1(2600)$, $D_2(2740)$ and $D^*_3(2750)$ states. The mixing between the $D_1(2420)$ and $D_1(2430)$ resonances is studied and the mixing parameters are measured. The dataset corresponds to an integrated luminosity of 4.7 $fb^{-1}$, collected in proton-proton collisions at center-of-mass energies of 7, 8 and 13 TeV with the LHCb detector.

1 Introduction details on the fits to the data. The measurements of the partial branching fractions are given in Sec. 9 and results are summarized in Sec. 10.

LHCb detector
The LHCb detector [15,16] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the pp interaction region, 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 the momentum of charged particles with relative uncertainty that varies from 0.5 % at low momentum (3 GeV) to 1.0 % at 200 GeV. The minimum distance of a track to a primary vertex (PV) (defined as the location of a reconstructed pp collision) 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. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
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 at √ s = 7 TeV (8)(9)(10)(11)(12)(13). The software trigger used in this analysis requires a two-, three-or four-track secondary vertex with significant displacement from the primary pp interaction vertices. At least one charged particle must have p T > 1.7 GeV and be inconsistent with originating from any PV. A multivariate algorithm [17] is used for the identification of secondary vertices consistent with the decay of a b hadron.
In the offline selection, the objects that fired the trigger are associated with reconstructed particles. Selection requirements can therefore be made not only on the software trigger that fired, but also on whether the decision was due to the signal candidate, other particles produced in the pp collision, or a combination of both. Both cases are retained for further analysis.
Simulated samples are used to characterize the detector response to signal and certain types of backgrounds. In the simulation, pp collisions are generated using Pythia [18] with a specific LHCb configuration [19]. Decays of hadronic particles are described by EvtGen [20], in which final state radiation is generated using Photos [21]. The interaction of the generated particles with the detector and its response are implemented using the Geant4 toolkit [22] as described in Ref. [23].

Selection requirements
The selection of the B − meson is performed using the decay chain where X represents a system composed of any collection of charged or neutral particles. After applying selections on the quality of the B − candidate tracks, further requirements are applied on their momenta, p, and transverse momenta, p T . The D 0 meson is reconstructed through its K − π + decay, applying loose particle identification criteria on both particles and good vertex quality requirements. The remaining tracks associated to the B − final state form a π + π − π − system which defines the B − decay vertex. Very loose particleidentification criteria are applied to the three pions together with good vertex-quality and impact-parameter constraints. The invariant mass of the above π + π − π − system is required to be below the physical boundary m(π + π − π − ) < 3.6 GeV. In the data collected at √ s = 7 and 8 TeV (48.5% of the total dataset), the requirement is m(π + π − π − ) < 3.0 GeV, which also removes 1.2% of the signal. Although the loss in the Run 1 data is rather small, it produces a nonnegligible distortion in the B − Dalitz plot and in the π − π − invariant-mass distribution.
The momentum scale is calibrated using samples of J/ψ → µ + µ − and B + → J/ψ K + decays collected concurrently with the data sample used for this analysis [24,25]. The relative accuracy of this procedure is estimated to be 3 × 10 −4 using samples of other fully reconstructed b hadrons, Υ and K 0 S mesons. Figure 1(a) shows the D 0 π + mass spectrum, computed as m(K − π + π + s ) − m(K − π + ) + m PDG D 0 (here π + s indicates the "slow pion" from the D * + decay and m PDG D 0 indicates the known D 0 mass value), where a clean D * + signal can be observed.
The D * + candidate is selected within 3.5σ of the fitted D * + mass value, where σ = 0.45 MeV is the effective mass resolution obtained from a fit to the D 0 π + mass spectrum with a sum of two Gaussian functions. The B − → D * + π − π − decay is affected by background from B 0 → D * + π − decays combined with a random π − candidate in the event. This contribution populates the high-mass sideband of the B − signal and is removed if either D * + π − have a mass within 2σ the B 0 known value [3], where σ = 18.5 MeV is obtained from a Gaussian fit to the B 0 mass distribution. The resulting D * + π − π − mass and (b) D * + π − π − invariant masses for candidates after the selection on the χ 2 /ndf from the fit to the B − decay tree. distribution is shown in Fig 1(b), where the B − signal can be observed over a significant background.
A significant source of background is due to B − → D * + π − π − π 0 or B 0 → D * + π − π − π + decays, where a pion is not reconstructed. However the B − → D * + π − π − mass combination from these final states populate the low-mass sideband of the B − signal and do not extend into the signal region. A coherent source of background which could affect the signal region is due to B 0 → D * + π − π 0 decays, where the π 0 meson is not reconstructed but is replaced by a random π − candidate from the event. In this case the D * + π − system could have a definite spin-parity configuration, however this contribution is found to be negligible. A possible source of background comes from the B − → D * + K − π − decay, where the K − candidate is misidentified as a pion. This background is kinematically confined in the lower sideband of the B − signal. Its contribution relative to B − → D * + π − π − has been measured in Ref. [26] and found to be negligible. Therefore the background under the B − signal is dominated by combinatorial background.
To reduce the combinatorial background while keeping enough signal for an amplitude analysis, a multivariate selection is employed, in the form of a likelihood ratio defined, for each event, as where i runs over a set of 6 variables and P s (i) and P b (i) are probability density functions (PDFs) of the signal and background contributions, respectively. The signal PDFs are obtained from simulated signal samples, while background PDFs are obtained from the B − sideband regions, defined within 4.5 − 6.6σ on either side of the B − mass peak, where σ is obtained from the fit to the D * + π − π − mass spectrum defined below. The variables used are: the B − decay length significance, defined as the ratio between the decay length and its uncertainty; the B − transverse momentum; the χ 2 of the primary vertex associated to the B − meson; the B − and D * + impact parameters with respect to the primary vertex; and the χ 2 /ndf from the fit to the B − → (D 0 (→ K − π + )π + s )π − π − decay tree. The choice of the selection value on the variable R is performed using an optimization procedure where the D * + π − π − mass spectrum of candidates selected with increasing cut on R is fitted. The fits are performed using two Gaussian functions with a common mean to describe the B − signal and a quadratic function for the background. Defining σ as the weighted mean of the widths of the two Gaussian functions, the signal region is defined within 2σ, where σ = 17.7 MeV. For each selection the fit estimates the signal and background yields, N sig and N bkg . From these quantities the purity p = N sig N sig +N bkg and the significance s = N sig √ N sig +N bkg are evaluated. To obtain both the largest purity and significance, the figure of merit s · p is evaluated. It has its maximum at R > 0.5 which is taken as the default selection. For this selection, Fig. 2 shows the resulting D * + π − π − mass spectrum where the B − signal is observed over small background.
For the above selection the signal purity is p = 0.9, while the efficiency is 81.9%. The yield within the signal region is 79 120, of which 48.5% and 51.5% are from Run 1 and Run 2, respectively. The purities of the two data sets are found to be the same. The number of events with multiple B − candidate combinations is negligible.
The B − → D * + π − π − decay mode contains two indistinguishable π − mesons, giving two D * + π − mass combinations. In the following, m(D * + π − ) low and m(D * + π − ) high indicate the lower and higher values of the two mass combinations, respectively. The B − Dalitz plot, described as a function of m 2 (D * + π − ) high and m 2 (D * + π − ) low , is shown in Fig. 3 for Run 2 data.   The Dalitz plot contains clear vertical bands in the 6 GeV 2 mass region, due to the presence of the well-known D 1 (2420) and D * 2 (2460) resonances. The presence of further weaker bands can be observed in the higher mass region. The prominent presence of the above two resonances can be observed in the m(D * + π − ) low projection, shown in Fig. 4 for the total dataset. On the other hand, the presence of additional states is rather weak in the mass projection and therefore an angular analysis is needed to separate the different contributions.
The following angles are useful in discriminating between different J P contributions: θ H , the helicity angle defined as the angle between the π + s direction in the D 0 π + s rest frame and the D 0 π + s direction in the D 0 π + s π − rest frame (see Fig. 5(a)); θ, the helicity angle defined as the angle formed by the π − 1 direction in the D * + π − 1 rest frame and the D * + π − 1 direction in the D * + π − 1 π − 2 rest frame (see Fig. 5(b)); γ, the angle in the D * + π − π − rest frame formed by the π + s direction in the D 0 π + s rest frame and the normal to the D * + π − π − plane (see Fig. 5(c)).
The angle θ H is useful to discriminate between natural and unnatural spin-parity resonances for which the expected angular distributions are sin 2 θ H and 1 − h cos 2 θ H (where h < 1 depends on the properties of the decay), respectively, except for J P = 0 − where a cos 2 θ H distribution is expected. Figure 6 shows the two-dimensional distribution of cos θ H vs. m(D * + π − ) low . The two vertical bands are due to the D 1 (2420) and D * 2 (2460) states which exhibit the expected 1 − h cos 2 θ H and sin 2 θ H distributions, respectively.

Background model
The background model is obtained from the data in the signal region using the method of signal subtraction. Using the R variable defined in Eq. (2), two datasets are extracted, (a) with low purity (R > 0.0, p a = 0.865 and signal yield N a = 77 644), (b) with high purity (R > 2.5, p b = 0.949 and signal yield N b = 34 019). The background distribution for a given variable is then obtained by subtracting the high-purity distributions, scaled by the factor N a /N b , from the low-purity distributions. The variables m(D * + π − ) low and cos θ are found to be independent and different for signal and background, therefore the resulting background model is obtained by the product of the PDFs of these distributions. Figure 7(a) shows the m(D * + π − ) low distribution for the low-purity (filled points) and high-purity (open points) selections. Figure 7(b) shows the signal-subtracted distribution, where no significant resonant structure is seen. The superimposed curve is the result from a fit performed using the sum of two exponential functions multiplied by the B − phase-space factor. Similarly, Fig. 7(c) shows the cos θ distribution for the low-purity and high-purity samples, and Fig. 7(d) shows the signal-subtracted distribution. The curve is the result of a fit using a 6 th order polynomial function.

Efficiency
The efficiency is computed using simulated samples of signal decays analyzed using the same procedure as for data. Due to the different trigger and reconstruction methods, the efficiency is computed separately for Run 1 and Run 2 data. It is found that the efficiency mainly depends on the variables, m(D * + π − ) and cos θ. Weak or no dependence is found on other variables. The efficiency model is obtained by dividing the simulated sample into 22 slices in log(m(D * + π − )/ MeV) and fitting the cos θ distribution for each slice using 5 th order polynomial functions. The efficiency for a given value of log(m(D * + π − )/ MeV) is then obtained by linear interpolation between two adjacent bins. Figure 8 shows the interpolated efficiency maps in the (log(m(D * + π − )/ MeV), cos θ) plane, separately for Run 1 and Run 2. The empty region in Run 1 data is caused by the requirement m(π + π − π − ) < 3.0 GeV. Although this region is populated by a small fraction of signal, estimated using Run 2 data, this introduces some uncertainty in the description of the Run 1 data.
The mass resolution is studied as function of m(D * + π − ) using simulation. For each slice in m(D * + π − ) the difference between the generated and reconstructed mass is computed and the resulting distributions are fitted using the sum of two Gaussian functions. The effective resolution σ, increases almost linearly from σ = 4 MeV at 2.4 GeV to σ = 7 MeV at 3.0 GeV. This value of the mass resolution is much smaller than the minimum width of the known resonances present in the D * + π − mass spectrum (for D 1 (2420), Γ = 31 MeV), therefore, in the following, resolution effects are ignored.

Amplitude analysis
An amplitude analysis of the four-body decay where π − 1 and π 2 indicate the two indistinguishable pions, is performed to extract the fractions and the phases of the charmed resonances contributing to the decay and to measure their parameters and quantum numbers. All the amplitudes are symmetrized according to the exchange of the π − 1 with π − 2 mesons.

Description of the amplitudes
The amplitudes contributing to the decay are parameterized using the nonrelativistic Zemach tensors formalism [27][28][29]. It is assumed that reaction (3) proceeds as where R 0 is an intermediate charmed meson resonance which decays as Reaction (4) is a weak decay and does not conserve parity while reaction (5) is a strong decay and conserves both angular momentum and parity. The four particles in the final state are labeled as In the description of the decay R 0 → D 0 π + s π − 1 , the 3-vectors p i (i=1,2,3) indicate the momenta of the particles in the D 0 π + s π − 1 rest frame and L indicates the angular momentum between the D 0 π + s system and the π − 1 meson. For the resonance D * + decaying as D * + → D 0 π + s , the decay products D 0 and π + s having 3-momenta p 1 , p 2 and masses m 1 and m 2 , the 3-vector t 3 is defined as with m 12 indicating the D 0 π + s invariant mass. To describe the decay B − → R 0 π − 2 , the 3-vector q 4 indicates the momentum of π − 2 in the B − rest frame. The amplitudes are obtained as follows: • a symmetric and traceless tensor of rank L, P , constructed with p 3 is used to describe the orbital angular momentum L; • a symmetric and traceless tensor of rank S, T , constructed with t 3 is used to describe the spin of intermediate resonances; • the tensors P and T are then combined into a tensor J of rank J to obtain the total spin J of the D 0 π + π − 1 system; • a symmetric and traceless tensor, Q, of rank J constructed with q 4 is used to describe the orbital angular momentum between R 0 and π − 2 ; • the scalar product of the two tensors J and Q gives the scalar which represents the 0 spin of the B − meson.
The resonance R 0 , having a given J P , decays into a 1 − resonance (D * + ) and a 0 − particle with a given orbital angular momentum L. In a first approach, the resonance lineshape is described by a complex relativistic Breit-Wigner function, BW(m 123 ), with appropriate Blatt-Weisskopf centrifugal barrier factors [14,30] which are computed assuming a radius 4.5 GeV −1 . In a second approach the resonance lineshapes are described by the quasi-model-independent method (QMI) [31,32], described later.
The list of the amplitudes used in the present analysis is given in Table 1. The J P = 0 − nonresonant contribution term is omitted because it is found to be negligible.

Amplitude analysis fit
The amplitude analysis of the D * + π − π − candidates in the B − mass region is performed using unbinned maximum-likelihood fits. The likelihood function is written as where • N is the number of candidates in the signal region; • for the n th event, x n is the set of variables describing the 4-body B meson decay; • (x n ) is the efficiency function; • A i (x n ) represents the i th complex signal-amplitude contribution; • c i is the complex intensity of the i th signal component; the c i parameters are allowed to vary during the fit process; numerical integration is performed on phase-space-generated decays with the B − signal lineshape generated according to the experimental distribution; • p is the signal purity obtained from the fit to the D * + π − π − mass spectrum; • B(z n ) is the normalized background contribution, parameterized as a function of the two variables described in Sec. 5.2.
The efficiency-corrected fraction f i due to a resonant or nonresonant contribution i is defined as follows The f i fractions do not necessarily sum to 1 because of interference effects. The uncertainty for each f i fraction is evaluated by propagating the full covariance matrix obtained from the fit. Similarly, the efficiency-corrected interference fractional contribution f ij , for i < j is defined as The amplitude analysis is started by including, one by one, all the possible charmed resonance contributions with masses and widths listed in Ref. [3]. Resonances are kept if a significant likelihood increase (∆(2 log L) > 3) is observed. The list of the states giving significant contributions at the end of the process is given in the upper section of Table 2. The fit procedure is tested on pseudo-experiments using different combinations of amplitudes, input fractions and phases, obtaining a good agreement between generated and fitted values.
The quality of the description of the data is tested by the χ 2 /ndf, defined as the sum of two χ 2 values, calculated from the two two-dimensional distributions (m , cos θ H ) and (cos θ, cos γ) as Here N model i and N data i are the fit predictions and observed yields in each bin of the two-dimensional distributions. The variable ndf is defined as ndf = N cells − N par , where N cells is the number of bins having at least 6 entries and N par is the number of free parameters in the fit. The variable m , defined in the range 0-1, is computed as

Fits to the data using quasi-model-independent amplitudes
It is found in several analyses that the mass terms of some amplitudes may not be well described by Breit-Wigner functions, because they are broad or because additional contributions may be present at higher mass. Therefore, for a given value of J P , a quasi-model-independent method is tested to describe the amplitude, while leaving all the other resonances described by Breit-Wigner functions. The method is also used to perform a scan of the mass spectrum to search for additional resonances. The D * + π − mass spectrum is divided into 31 slices with nonuniform bin widths and, for a given contributing resonance, the complex Breit-Wigner term is replaced by a set of 31 complex coefficients (magnitude anf phase) which are free to float. These values are fixed to arbitrary values in one bin, at a mass value in the 2.42 − 2.60 GeV range, depending on the amplitude and therefore the set of additional free parameters is reduced to 60.
The largest amplitude, usually the 1 + D amplitude, is taken as the reference wave. Due to the large number of fit parameters, QMI amplitudes can only be introduced one by one. The fit is performed using as free parameters the real and imaginary parts of the amplitude in each bin of the D * + π − mass spectrum. The search for the QMI parameters is performed using a random search, starting from zero in each mass bin. The fitted solution is then given as input to a second iteration modifying the value for the fixed bin to the average value obtained from the two adjacent bins. Obvious spikes are smoothed in the input of the second iteration. Normally the second iteration converges and is able to compute the full covariance matrix. The fitted QMI amplitude is then modeled through a cubic-spline interpolation function.
The method is tested using different initial values for the first iteration. In all cases the fit converges to the same solution. It is also tested with simulation obtaining good agreement between input and fitted values of the amplitudes.
The process starts with a QMI fit to the J P = 1 + S amplitude, including all the amplitudes listed in the upper part of Table 2 and described by Breit-Wigner functions with initial parameter values fixed to those reported in Ref. [3]. In this fit, due to significant Table 2: Resonance parameters from the amplitude analysis. The first uncertainty is statistical, the second systematic. The upper part reports the resonance parameters obtained from the amplitude analysis described in Sec. 7, the lower part those obtained from the mixing analysis described in Sec. 8.3. The labels indicating the spin-parity of D 0 (2550), D * 1 (2600), and D 2 (2740) resonances are updated, with respect to those reported in Ref. [3], according to the results from the amplitude analysis reported in this work. The D * 2 (2460) parameters are fixed to the world averages.
interference effects between the 1 + D, the narrow 1 + S and the broad 1 + S amplitudes, the narrow 1 + D/1 + S D 1 (2420) parameters, described by a Breit-Wigner function, are left free as well. The resulting parameters for D 1 (2420) resonance are given in Table 2. Statistical significances are computed as the fitted fraction divided by its statistical uncertainty. The QMI J P = 1 + S amplitude is then fixed and a QMI analysis of the J P = 0 − is performed. The process continues by fixing the J P = 0 − QMI amplitude and leaving, one by one, free Breit-Wigner parameters for all the resonances listed in the upper part of Table 2. The parameters of the D * 2 (2460) resonances are fixed to the world averages because they are well determined. The process is iterative, with QMI analyses of the J P = 1 + S and J P = 0 − amplitudes and free parameters for the resonances described by Breit-Wigner functions repeated several times, until the process converges and no significant variation of the free parameters is observed. The resulting fitted parameters of D * 1 (2600), D 2 (2740), and D * 3 (2750) amplitudes are listed in Table 2. To obtain the parameters of the broad D 1 (2430) resonance, a fit is performed with the QMI model for the J P = 1 + S amplitude replaced by the Breit-Wigner function model. Similarly, to obtain the D 0 (2550) parameters, the QMI model for the J P = 0 − amplitude is replaced by the Breit-Wigner function. The presence of a broad J P = 1 + D D 1 (2430) contribution has been tested but excluded from the final fit. Its effect, due to the presence the broad J P = 1 + S resonance, is to produce large interference effects so that the total fraction increases to large and rather unphysical values without significantly improving the fit quality. Figure 9 shows the fitted magnitude and phase of the 1 + S amplitude. The presence of a broad structure can be noted close to threshold with a corresponding phase motion as expected for a resonance. The magnitude and phase show further activity in the 2.8 GeV mass region, suggesting the presence of an additional 1 + S resonance. However, the introduction of a new 1 + S Breit-Wigner resonance with floating parameters in that mass region does not produce a significant contribution. The high-mass enhancement in the amplitude, on the other hand, is due to symmetrization effects due to the presence of two identical pions.
The J P = 0 − QMI magnitude and phase are shown in Fig. 10. In addition to the D 0 (2550) resonance, further activity can be seen in the 2.8 GeV mass region both in amplitude and phase, suggesting the presence of a possible new excited D 0 resonance. The J P = 0 − amplitude and phase distributions are fitted using the model where ps(m) is the D * + momentum in the D * + π − center-of-mass frame and a, c 1 , c 2 , α and β are free parameters. The parameters of the D 0 (2550) resonance (m 0 , Γ 0 ) are fixed to the values extracted from the amplitude analysis (see Table 2), while the parameters of the D 0 resonance, (m 1 , Γ 1 ) are free. The first term in the above equation represents a threshold J P = 0 − nonresonant term. The fit is performed in terms of real and imaginary parts of the amplitude and then converted into amplitude and phase when projected on the data in Fig. 10. The D 0 fitted parameters are m(D 0 ) = 2782 ± 13 MeV, Γ(D 0 ) = 146 ± 23 MeV (14) and the significance, computed as the ratio between the fitted fraction divided by the statistical uncertainty, is 3.2σ. However, an attempt to include this new possible resonance in the amplitude analysis gives a fraction consistent with zero.
To search for additional states, the QMI method is used for the most significant amplitudes, i.e. those with J P = 1 − (Fig. 11), J P = 1 + D (Fig. 12), and J P = 2 + (Fig. 13). In mass regions where the amplitude is consistent with zero the phase is not well measured and therefore statistical uncertainties are large. Superimposed on the QMI amplitudes are the Breit-Wigner functions, with arbitrary normalizations, describing the D * 1 (2600) (Fig. 11) and D 1 (2420) resonances (Fig. 12), respectively, using the fitted parameters given in Table 2. Similarly, the J P = 2 + amplitude is shown in Fig. 13 with D * 2 (2460) resonance parameters fixed to the values reported in Ref. [3].
A good agreement between the results from the QMI method and the expected lineshape of the Breit-Wigner description of the resonances is found. In the case of the J P = 1 − amplitude no additional structure can be seen, and the enhancement at high  mass can be associated to the reflection due to the presence of two identical π − mesons. Some amplitudes as J P = 1 + D, J P = 2 + and J P = 1 − evidence some points off from the Breit-Wigner behavior in the threshold region. Since in these regions phase space is limited, these effects can be due to cross-feeds from other partial waves.

Fit results
The data are fitted using three different models described below.
• The J P = 1 + S and J P = 0 − are described by QMI amplitudes. This model gives the best description of the data and is used to search for new states and obtain the Breit-Wigner parameters for several resonances.   • Mixing is allowed between the 1 + amplitudes. This model allows to measure the D 1 and D 1 Breit-Wigner parameters and their mixing angle and phase.

Results from the QMI model
In this fitting model the J P = 1 + S and J P = 0 − are described by QMI, while all the other amplitudes are described by relativistic Breit-Wigner functions with parameters given in Table 2. The results from the fit are given in Table 3. The dominance of the D 1 (2420) (1 + D) resonance can be noted, with important contributions from 1 + S QMI and D * 2 (2460) amplitudes. The sum of fractions is larger than 100%, indicating important interference effects.
The fit projections for Run 2 data (not biased by the π + π − π − mass cut) are shown in Fig. 14. Figure 15 shows the fit projections on m(D * + π − ) low using the total dataset together with all the contributing amplitudes and the significant interference contributions. Using statistical uncertainties only, the separate fits for Run 1 and Run 2 give χ 2 /ndf = 2348/1748 = 1.34 and χ 2 /ndf = 2111/1780 = 1.19, respectively. For a fit to the Table 3: Fit results from the amplitude analysis for the model where the J P = 1 + S and J P = 0 − amplitudes are described by QMI. The first uncertainty is statistical, the second systematic.

Resonance
J P fraction (%) phase (rad) D 1 (2420) 1 + D 59.8 ± 0.3 ± 2.9 0 128.2 ± 0.6 ± 3.8 total dataset χ 2 /ndf = 2551/1784 = 1.43. However it has to be taken into account that in this fit the total sample size is double and therefore statistical uncertainties are smaller. These χ 2 /ndf values indicate a good description of Run 2 data, but a worse description of the Run 1 data indicating some limitation in the handling of the efficiency for this data set.

Systematic uncertainties
Systematic uncertainties reported in Table 2 and alternative fit models (described later), whose results are given in Tables 9 and 10, are evaluated as follows. When multiple contributions are needed to describe a given effect, the average value of the absolute deviations from the reference fit is taken as a systematic uncertainty. Table 4 gives details on the contributions to the systematic uncertainties on the fractions and phases for the model where the J P = 1 + S and J P = 0 − amplitudes are described by QMI. The effect of the background (labeled as Purity) is studied by changing the selection cut corresponding to lower (with R > 1.1, p = 0.92, 66 064 candidates) or higher (with R > 0. The effect of the uncertainty of the background size (labeled as Bkg size) is estimated by modifying the value of the fixed purity value in the fit (90%) by ±0.5%. The effect of the small discrepancy between the data and fit projections on the cos θ and cos θ H distribution (labeled as Data/sim) is evaluated by weighting the efficiency distribution to match the data. The effect of the limited simulation sample (labeled as Sim) is evaluated by fitting the data using 100 binned 2-dimensional efficiency tables obtained from the  [14] (labeled as Mod) are included and excluded in the fit. The root-mean-square value of the deviations of the fraction from the reference fit are taken as systematic uncertainties. All the different contributions are added in quadrature. The dominant sources of systematic uncertainties are found to be due to the Blatt-Weisskopf radius. Table 5 gives details on the contributions to the systematic uncertainties for the measured masses and widths of the resonances contributing to the B − decay. In this case only the most relevant contributions are listed. From the study of large control samples, a systematic uncertainty of 0.0015 Q on the mass scale is added, where Q is the Q-value involved in the resonance decay.
The consistency between the Run 1 and Run 2 datasets is tested performing separate fits to the data and good agreement is obtained, within the uncertainties, on fractions and relative phases. Separate fits are performed to subsamples of the data where the D * + π − π − final state is directly (69%) or undirectly (31%) selected by the trigger conditions. The fitted fractions and phases are found consistent within the statistical uncertainties. A test is performed by weighting the simulated p T distribution to match the data and recomputing the efficiencies. The impact on the fitted fractions and phases is found to be negligible.

Legendre polynomial moments projections
A more detailed understanding of the resonant structures present in the D * + π − mass spectrum and of the agreement between data and fitting model is obtained by looking at the angular distributions as functions of cos θ, cos θ H and cos γ. This is obtained by weighting the D * + π − mass spectrum by the Legendre polynomial moments computed as functions of the above three angles. The D * + π − mass spectrum weighted by Legendre polynomial moments expressed as functions of cos θ is shown in Fig. 16 for L between 1 and 6 and reveals a rich structure. Higher moments are consistent with zero. Equations (15) relate the moments with orbital angular momentum between the D * + and π − mesons, assuming only partial waves between L = 0 and L = 3. Here S, P , D and F indicate the magnitudes of the amplitudes with angular momenta L = 0, 1, 2, 3 and φ denotes their relative phases.
A comparison with Table 1 allows for the identification of the resonant contributions to each distribution, listed in Table 6.
Significant interference effects between 1 + amplitudes can be observed in the Y 0 2 distribution and a clean 2 + signal due to D * 2 (2460) can be seen in the Y 0 4 distribution. Other moments show rather complex structures. An overall good description of the data is obtained, although some small discrepancy can be seen in Y 0 3 and Y 0 5 . This is expected, given the large number of physical contributions (see Table 6) to the shape of the Y 0 L moments which are also sensitive to efficiency effects.
Additional information can be obtained from the fit projections on the D * + π − mass spectrum weighted by Legendre polynomial moments computed as functions of cos θ H (labeled as Y H L ) and cos γ (labeled as Y γ L ), shown in Fig. 17(a) and Fig. 17(b) respectively. The two Y H 2 (cos θ H ) and Y γ 2 (cos γ) projections show interference effects between the D 1 (2420) and D * 2 (2460) resonances. The Y γ 2 distribution also shows an enhancement at Table 6: Relationship between the Legendre polynomial moments Y 0 L and spin amplitudes. In the column describing the interfering amplitudes, the left side amplitude is intended to interfere with any of the amplitudes listed on the right side.

Moment Squared amplitudes
Interfering amplitudes Y 0 the position of the D * 1 (2600) resonance. Other moments are consistent with zero.

Search for additional contributions and spin-parity determination
The presence of additional contributions is tested by adding them to the reference fit using the total dataset. The significance of each contribution is computed as its fitted fraction divided by its statistical uncertainty. No evidence is found for the D * 1 (2760) or D * 2 (3000) contributions, previously observed in the B − → D + π − π − decay [14]. Their statistical significance is found to be 2.4 and 0.0, respectively. Virtual contributions, as described in Ref. [14], are found to be small with a statistical significance of 4.4σ but ignored because they have a small fraction (0.12 ± 0.03)%, an uncertain physical meaning [14] and do not significantly improve the fit χ 2 .
The presence of a D 0 π + π − π − nonresonant contribution has been tested but excluded from the final fit. Its effect, due to the presence of broad J P = 1 + resonances, is to produce large interference effects so that the total fraction increases to large and rather unphysical values without significantly improving the fit quality.
It has been noted that the QMI 1 + S amplitude (Fig. 9) shows activity both in amplitude and phase in the mass region around 2.8 GeV which could correspond to the presence of an additional D 1 resonance. A test is performed including an additional 1 + S Breit-Wigner resonance in this mass region with free parameters. However, no significant contribution for this additional state is found.
The QMI approach is used for the most significant amplitudes and Breit-Wigner behavior is obtained for J P = 1 + , J P = 2 + , J P = 0 − and J P = 1 − resonances. For other contributions, such as the J P = 2 − D 2 (2740) or J P = 3 − D * 3 (2750) resonances, this is not possible due to the weakness of these contributions. For these two states, a spin analysis is performed. For each state additional fits are performed where the masses and widths are fixed to the results given in Table 2 but where the angular distributions are replaced by those from other possible spin assignments. For the D 2 (2740) resonance, J P = 0 − , 1 + D, 1 + S, 1 − , 2 + are tried but the likelihood and χ 2 variations exclude all the alternative hypotheses with significances greater than 8σ. The estimate of the significance is obtained using ∆χ 2 , where ∆χ 2 is the variation of the fit χ 2 for the given spin hypothesis. Similarly, for the D * 3 (2750) resonance, values of J P = 0 − , 1 + D, 1 + S, 1 − and 2 + are tried but excluded with significances greater than 8σ. In conclusion, the  Table 7. Note that different methods have been used to extract the resonances parameters. The results from the BaBar [11] and LHCb [12] collaborations come from inclusive studies of the D + π − , D 0 π + and D * + π − combinations where signals are fitted directly on the mass spectra. In the case of the D * + π − mass spectrum, resonance production is enhanced by the use of selections on the helicity angle θ H . Cross feeds from the resonance production in the D * + π − system are present in the D + π − and D 0 π + mass spectra. The LHCb results from Ref. [13] and [14], on the other hand, come from Dalitz plot analyses of B decays. Table 2 gives the resonance parameters for the D 1 (2430) and D 0 (2550) states when they are described by relativistic Breit-Wigner functions. An amplitude analysis performed using this model gives the results shown in Table 8. In this case χ 2 /ndf = 2348/1748 = 1.34 and χ 2 /ndf = 2205/1780 = 1.24 for Run 1 and Run 2 data, respectively. Somewhat reduced fractional contributions from the D 1 (2430) and D 0 (2550) resonances with respect to the QMI approach can be seen. This effect can be understood since in this model the J P = 1 + S and J P = 0 − contributions do not include possible additional contributions from higher mass resonances. Systematic uncertainties are evaluated as described in Sec. 7.

8.3
Results from the J P = 1 + mixing model A consequence of the heavy-quark symmetry is that, in the infinite-mass heavy quark limit, heavy-light Qq mesons can be classified in doublets labeled by the value of the total angular momentum j q of the light degrees of freedom with respect to the heavy quark Q [33]. In the quark model j q would be given by j q = s q + L, where L is the light quark orbital angular momentum. The Heavy Quark Effective Theory predicts that the two  J P = 1 + mesons, with j q = 1 2 and j q = 3 2 , decay into the D * π final state via the S-and D-wave, respectively. Due to the finite c-quark mass, the observed physical 1 + states can be a mixture of such pure states. The mixing can occur for instance via the common D * π decay channel and the resulting D 1 and D 1 amplitudes are a superposition of the S-and D-wave amplitudes where ω is the mixing angle and ψ is a complex phase. In this model the J P = 1 + D 1 amplitude is taken as reference. The J P = 1 + S amplitudes are described by relativistic Breit-Wigner functions with free parameters, while the J P = 0 − amplitude is described by the QMI model. All the other resonances are described by relativistic Breit-Wigner functions with parameters fixed to the values reported in Table 2. Table 9 gives details on the fractions and relative phases.
Combining statistical and systematic uncertainties the mixing angle deviates from zero by 2.3σ. The χ 2 /ndf for the fit to the total dataset is χ 2 /ndf = 2739/1780 = 1.54. Systematic uncertainties on the mixing parameters, fractional contributions and relative phases are computed as described in Sec. 7. The measured D 1 and D 1 masses and widths are reported in Table 2.

Measurement of the branching fractions
The known branching fraction of the B − → D * + π − π − decay mode is B(B − → D * + π − π − ) = (1.35 ± 0.22) × 10 −3 [3]. Table 10 reports the partial branching fractions for the resonances contributing to the total branching fraction. They are obtained multiplying the B − branching fraction by the fractional contributions obtained from the amplitude analysis performed using the Breit-Wigner model for all the resonances and reported in Table 8. For the D 1 and D 1 branching fractions the fractional contributions obtained from the mixing model and reported in Table 9 are used. Since the uncertainty on the absolute branching fraction is large, it has been separated from the other sources of systematic uncertainty. The D 1 (2420) resonance decays to D-and S-wave states and therefore the two contributions are added; a similar procedure is followed for the D 2 (2740) resonance, which decays to P -and F -wave states. Table 10: Summary of the measurements of the branching fractions. The first uncertainty is statistical, the second systematic, the third is due to the uncertainty on the measurement of the B − → D * + π − π − absolute branching fraction. A comparison with measurements obtained by the Belle collaboration [4] is shown.

Summary
A four-body amplitude analysis of the B − → D * + π − π − decay is performed using pp collision data, corresponding to an integrated luminosity of 4.7 fb −1 , collected at centerof-mass energies of 7, 8 and 13 TeV with the LHCb detector.