Analysis of the resonant components in B0->J/\psi pi+pi-

Interpretation of CP violation measurements using charmonium decays, in both the B0 and Bs systems, can be subject to changes due to"penguin"type diagrams. These effects can be investigated using measurements of the Cabibbo-suppressed B0->J/\psi pi+pi- decays. The final state composition of this channel is investigated using a 1.0/fb sample of data produced in 7 TeV pp collisions at the LHC and collected by the LHCb experiment. A modified Dalitz plot analysis is performed using both the invariant mass spectra and the decay angular distributions. An improved measurement of the B0->J/\psi pi+pi- branching fraction of (3.97 +/-0.09+/- 0.11 +/- 0.16)x10^{-5} is reported where the first uncertainty is statistical, the second is systematic and the third is due to the uncertainty of the branching fraction of the decay B- ->J/\psi K- used as a normalization channel. In the J/\psi pi+pi- final state significant production of f0(500) and rho(770) resonances is found, both of which can be used for CP violation studies. In contrast evidence for the f0(980) resonance is not found, thus establishing the first upper limit on the branching fraction product B(B0->J/\psi f0(980) x B(f0(980)->pi+ pi-)<1.1x10^{-6}, leading to an upper limit on the absolute value of the mixing angle of the f0(980)$ with the f0(500) of<31 degrees, both at 90% confidence level.


Introduction
CP violation measurements using neutral B meson decays into J/ψ mesons are of prime importance both for determinations of Standard Model (SM) parameters and searching for physics beyond the SM. In the case of B 0 decays, the final state J/ψ K 0 S is the most important for measuring sin 2β [1], while in the case of B 0 s decays, used to measure φ s , only the final states J/ψ φ [2][3][4], and J/ψ π + π − [5] have been used so far, where the largest component of the latter is J/ψ f 0 (980) [6]. The decay rate for these J/ψ modes is dominated by the color-suppressed tree level diagram, an example of which is shown for B 0 decays in Fig. 1(a), while penguin processes, an example of which is shown in Fig. 1(b), are expected to be suppressed. Theoretical predictions on the effects of such "penguin pollution" vary widely for both B 0 and B 0 s decays [7], so it is incumbent upon experimentalists to limit possible changes in the value of the CP violating angles measured using other decay modes. The decay B 0 → J/ψ π + π − can occur via a Cabibbo suppressed tree level diagram, shown in Fig. 2(a), or via several penguin diagrams. An example is shown in Fig. 2(b), while others are illustrated in Ref. [8]. These decays are interesting because they can also be used to measure or limit the amount of penguin pollution. The advantage in using the decay B 0 → J/ψ π + π − arises because the relative amount of pollution is larger. In the allowed decays, e.g. B 0 → J/ψ K 0 S , the penguin amplitude is multiplied by a factor of λ 2 Re iφ , where λ is the sine of the Cabibbo angle (≈ 0.22), while in the suppressed decays the factor becomes R e iφ , where R and R , and φ and φ are expected to be similar in size [8]. A similar study uses the decay B 0 s → J/ψ K 0 S [9]. CP violation measurements in the J/ψ π + π − mode utilizing B 0 − B 0 mixing determine sin 2β eff which can be compared to the well measured sin 2β. Differences can be used to estimate the magnitude of penguin effects. Knowledge of the final state structure is the first step in this program. Such measurements on sin 2β eff have been attempted in the system by using the J/ψ π 0 final state [10].
In order to ascertain the viability of such CP violation measurements we perform a full "Dalitz like" analysis of the final state. Regions in π + π − mass that correspond to spin-0 final states would be CP eigenstates. Final states containing vector resonances, such as the ρ(770) can be analyzed in a similar manner as was done for the decay B 0 s → J/ψ φ [2][3][4]. It is also of interest to search for the f 0 (980) contribution and to obtain information concerning the mixing angle between the f 0 (980) and the f 0 (500) 1 partners in the scalar nonet, as the latter should couple strongly to the dd system. Branching fractions for B 0 → J/ψ π + π − and J/ψ ρ 0 have previously been measured by the BaBar collaboration [11].
In this paper the J/ψ π + and π + π − mass spectra and decay angular distributions are used to determine the resonant and non-resonant components. This differs from a classical Dalitz plot analysis [12] because one of the particles in the final state, the J/ψ meson, has spin-1 and its three decay amplitudes must be considered. We first show that there are no evident structures in the J/ψ π + invariant mass, and then model the π + π − invariant mass with a series of resonant and non-resonant amplitudes. The data are then fitted with the coherent sum of these amplitudes. We report on the resonant structure and the CP content of the final state.

Data sample and selection requirements
The data sample consists of 1.0 fb −1 of integrated luminosity collected with the LHCb detector [13] using pp collisions at a center-of-mass energy of 7 TeV. The detector is a singlearm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. Components include a high precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift-tubes placed downstream. The combined tracking system has a momentum 2 resolution ∆p/p that varies from 0.4% at 5 GeV to 0.6% at 100 GeV, and an impact parameter resolution of 20 µm for tracks with large transverse momentum (p T ) with respect to the proton beam direction. Charged hadrons are identified using two ring-imaging Cherenkov (RICH) detectors. Photon, electron and hadron candidates are identified by a calorimeter system consisting of scintillating-pad and 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 that applies a full event reconstruction [14].
Events are triggered by a J/ψ → µ + µ − decay, requiring two identified muons with opposite charge, p T (µ ± ) greater than 500 MeV, an invariant mass within 120 MeV of the J/ψ mass [15], and form a vertex with a fit χ 2 less than 16. After applying these requirements, there is a large J/ψ signal over a small background [16]. Only candidates with dimuon invariant mass between −48 MeV and +43 MeV relative to the observed J/ψ mass peak are selected, corresponding a window of about ±3σ. The requirement is asymmetric because of final state electromagnetic radiation. The two muons subsequently are kinematically constrained to the known J/ψ mass.
Other requirements are imposed to isolate B 0 candidates with high signal yield and minimum background. This is accomplished by combining the J/ψ → µ + µ − candidate with a pair of pion candidates of opposite charge, and then testing if all four tracks form a common decay vertex. Pion candidates are each required to have p T greater than 250 MeV, and the scalar sum of the two transverse momenta, p T (π + ) + p T (π − ), must be larger than 900 MeV. The impact parameter (IP) is the distance of closest approach of a track to the primary vertex (PV). To test for inconsistency with production at the PV, the IP χ 2 is computed as the difference between the χ 2 of the PV reconstructed with and without the considered track. Each pion must have an IP χ 2 greater than 9. Both pions must also come from a common vertex with an acceptable χ 2 and form a vertex with the J/ψ with a χ 2 per number of degrees of freedom (ndf) less than 10 (here ndf equals five). Pion and kaon candidates are positively identified using the RICH system. Cherenkov photons are matched to tracks, the emission angles of the photons compared with those expected if the particle is an electron, pion, kaon or proton, and a likelihood is then computed. The particle identification makes use of the logarithm of the likelihood ratio comparing two particle hypotheses (DLL). For pion selection we require DLL(π − K) > −10.
The four-track B 0 candidate must have a flight distance of more than 1.5 mm, where the average decay length resolution is 0.17 mm. The angle between the combined momentum vector of the decay products and the vector formed from the positions of the PV and the decay vertex (pointing angle) is required to be less than 2.5 • .
Events satisfying this preselection are then further filtered using a multivariate analyzer based on a Boosted Decision Tree (BDT) technique [17]. The BDT uses six variable that are chosen in a manner that does not introduce an asymmetry between either the two muons or the two pions. They are the minimum DLL(µ − π) of the µ + and µ − , the minimum p T of the π + and π − , the minimum of the IP χ 2 of the π + and π − , the B 0 vertex χ 2 , the B 0 pointing angle, and the B 0 flight distance. There is discrimination power between signal and background in all of these variables, especially the B 0 vertex χ 2 .
The background sample used to train the BDT consists of the events in the B 0 mass sideband having 5566 < m(J/ψ π + π − ) < 5616 MeV. The signal sample consists of two million B 0 → J/ψ (→ µ + µ − )π + π − Monte Carlo simulated events that are generated uniformly in phase space, using Pythia [18] with a special LHCb parameter tune [19], and the LHCb detector simulation based on Geant4 [20] described in Ref. [21]. Separate samples are used to train and test the BDT. The distributions of the BDT classifier for signal and background are shown in Fig. 3. To minimize a possible bias on the signal acceptance due to the BDT, we choose a relatively loose requirement of the BDT classifier > 0.05 which has a 96% signal efficiency and a 92% background rejection rate.
The invariant mass of the selected J/ψ π + π − combinations, where the dimuon pair is constrained to have the J/ψ mass, is shown in Fig. 4 Figure 4: Invariant mass of J/ψ π + π − combinations. The data are fitted with a double-Gaussian signal and several background functions. The (red) solid double-Gaussian function centered at 5280 MeV is the B 0 signal, the (brown) dotted line shows the combinatorial background, the (green) short-dashed shows the B − background, the (purple) dot-dashed line shows the contribution of B 0 s → J/ψ π + π − decays, the (black) dot-long dashed is the sum of B 0 s → J/ψ η (→ ργ) and B 0 s → J/ψ φ(→ π + π − π 0 ) backgrounds, the (light blue) long-dashed is the B 0 → J/ψ K − π + reflection, and the (blue) solid line is the total. the B 0 s and B 0 masses on top of the background. Double-Gaussian functions are used to fit both signal peaks. They differ only in their mean values, which are determined by the data. The core Gaussian width is also allowed to vary, while the fraction and width ratio of the second Gaussian is fixed to that obtained in the fit of B 0 s → J/ψ φ events. (The details of the fit are given in Ref. [6].) Other components in the fit model take into account background contributions. One source is from B − → J/ψ K − decays, which contributes when the K − is misidentifed as a π − and then combined with a random π + ; the smaller J/ψ π − mode contributes when it is combined with a random π + . The next source contains B 0 s → J/ψ η (→ ργ) and B 0 s → J/ψ φ(→ π + π − π 0 ) decays where the γ and the π 0 are ignored respectively. Finally there is a B 0 → J/ψ K − π + reflection where the K − is misidentified as π − . Here and elsewhere charged conjugated modes are included when appropriate. The exponential combinatorial background shape is taken from same-sign combinations, that are the sum of J/ψ π + π + and J/ψ π − π − candidates. The shapes of the other components are taken from the simulation with their normalizations allowed to vary. The fit gives 5287 ± 112 signal and 3212 ± 80 background candidates within ±20 MeV of the B 0 mass peak, where a K 0 S veto, discussed later, is applied. We use the well measured B − → J/ψ K − mode as a normalization channel to determine the branching fractions. To minimize the systematic uncertainty from the BDT selection, we employ a similar selection on B − → J/ψ K − decays after requiring the same preselection except for particle identification criteria on the K − candidates. Similar variables are used for the BDT except that the variables describing the combination of π + and π − in the J/ψ π + π − final state are replaced by ones describing the K − meson. For BDT training, the signal sample uses simulated events and the background sample consists of the data events in the region 5400 < m(J/ψ K − ) < 5450 MeV. The resulting invariant mass distribution of the candidates satisfying BDT classifier > 0.05 is shown in Fig. 5. Fitting the distribution with a double-Gaussian function for the signal and linear function for the background gives 350,727 ± 633 signal and 4756 ± 103 background candidates within ±20 MeV of the B − mass peak.

Analysis formalism
We apply a formalism similar to that used in Belle's analysis [22] of B 0 → K − π + χ c1 decays and later used in LHCb's analysis of B 0 s → J/ψ π + π − decays [6]. The decay B 0 → J/ψ π + π − , with J/ψ → µ + µ − , can be described by four variables. These are taken to be the invariant mass squared of J/ψ π + (s 12 ≡ m 2 (J/ψ π + )), the invariant mass squared of π + π − (s 23 ≡ m 2 (π + π − )), where we use label 1 for J/ψ , 2 for π + and 3 for π − , the J/ψ helicity angle (θ J/ψ ), which is the angle of the µ + in the J/ψ rest frame with respect to the J/ψ direction in the B 0 rest frame, and the angle between the J/ψ and π + π − decay planes (χ) in the B 0 rest frame. To improve the resolution of these variables we perform a kinematic fit constraining the B 0 and J/ψ masses to their nominal values [15], and recompute the final state momenta. To simplify the probability density function, we analyze the decay process after integrating over χ, which eliminates several interference terms.

The decay model for
The overall probability density function (PDF) given by the sum of signal, S, and background, B, functions is where f sig is the fraction of the signal in the fitted region and ε is the detection efficiency. The fraction of the signal is obtained from the mass fit and is fixed for the subsequent analysis. The normalization factors are given by N sig = ε(s 12 , s 23 , θ J/ψ )S(s 12 , s 23 , θ J/ψ ) ds 12 ds 23 d cos θ J/ψ , The event distribution for m 2 (π + π − ) versus m 2 (J/ψ π + ) in Fig. 6 shows obvious structure in m 2 (π + π − ). To investigate if there are visible exotic structures in the J/ψ π + system as claimed in similar decays [23], we examine the J/ψ π + mass distribution shown in Fig. 7 (a). No resonant effects are evident. Figure 7 (b) shows the π + π − mass distribution. There is a clear peak at the ρ(770) region, a small bump around 1250 MeV, but no evidence for the f 0 (980) resonance. The favored B 0 → J/ψ K 0 S decay is mostly rejected by the B 0 vertex χ 2 selection, but about 150 such events remain. We eliminate them by excluding the candidates that have |m(  Figure 7: Distribution of (a) m(J/ψ π + ) and (b) m(π + π − ) for B 0 → J/ψ π + π − candidate decays within ±20 MeV of B 0 mass shown with the solid line. The (red) points with error bars show the background contribution determined from m(J/ψ π + π − ) fits performed in each bin.

The signal function
The signal function for B 0 is taken to be the coherent sum over resonant states that can decay into π + π − , plus a possible non-resonant S-wave contribution 3 where A R i λ (s 12 , s 23 , θ J/ψ ) is the amplitude of the decay via an intermediate resonance R i with helicity λ. Each R i has an associated amplitude strength a R i λ for each helicity state λ and a phase φ R i λ . Note that the spin-0 component can only have a λ = 0 term. The amplitudes for each i are defined as where P B is the J/ψ momentum in the B 0 rest frame and P R is the momentum of either of the two pions in the dipion rest frame, m B is the B 0 mass, F are the B 0 meson and R resonance Blatt-Weisskopf barrier factors [24], L B is the orbital angular momentum between the J/ψ and π + π − system, and L R is the orbital angular momentum in the π + π − decay and is equal to the spin of resonance R because pions have spin-0. Since the parent B 0 has spin-0 and the J/ψ is a vector, when the π + π − system forms a spin-0 resonance, L B = 1 and L R = 0. For π + π − resonances with non-zero spin, L B can be 0, 1 or 2 (1, 2 or 3) for L R = 1(2) and so on. We take the lowest L B as the default and consider the other possibilities in the systematic uncertainty.
The Blatt-Weisskopf barrier factors F For the B meson z = r 2 P 2 B , where the hadron scale r is taken as 5.0 GeV −1 , and for the R resonance z = r 2 P 2 R with r taken as 1.5 GeV −1 [25]. In both cases z 0 = r 2 P 2 0 where P 0 is the decay daughter momentum calculated at the resonance pole mass.
The angular term, T λ , is obtained using the helicity formalism and is defined as where d is the Wigner d-function, J is the resonance spin, θ ππ is the π + π − resonance helicity angle which is defined as the angle of the π + in the π + π − rest frame with respect to the π + π − direction in the B 0 rest frame and calculated from the other variables as The J/ψ helicity dependent term Θ λ (θ J/ψ ) is defined as The function A R (s 23 ) describes the mass squared shape of the resonance R, that in most cases is a Breit-Wigner (BW) amplitude. Complications arise, however, when a new decay channel opens close to the resonant mass. The proximity of a second threshold distorts the line shape of the amplitude. This happens for the f 0 (980) resonance because the K + K − decay channel opens. Here we use a Flatté model [26] which is described below.
The BW amplitude for a resonance decaying into two spin-0 particles, labeled as 2 and 3, is where m R is the resonance pole mass, Γ(s 23 ) is its energy-dependent width that is parametrized as Here Γ 0 is the decay width when the invariant mass of the daughter combinations is equal to m R . The Flatté model is parametrized as The constants g ππ and g KK are the f 0 (980) couplings to ππ and KK final states respectively. The ρ factors account for the Lorentz-invariant phase space and are given as For non-resonant processes, the amplitude A(s 12 , s 23 , θ J/ψ ) is derived from Eq. 4, considering that the π + π − system is S-wave (i.e. L R = 0, L B = 1) and A R (s 23 ) is constant over the phase space s 12 and s 23 . Thus, it is parametrized as

Detection efficiency
The detection efficiency is determined from a sample of two million B 0 → J/ψ (→ µ + µ − )π + π − simulated events that are generated uniformly in phase space. Both s 12 and s 13 are centered at about 18.4 GeV 2 . We model the detection efficiency using the symmetric dimensionless Dalitz plot observables These variables are related to s 23 since The acceptance in cos θ J/ψ is not uniform, but depends on s 23 , as shown in Fig. 8. If the efficiency was independent of s 23 , then the curves would have the same shape. On the other hand, no clear dependence on s 12 is seen. Thus the efficiency model can be expressed as To study the cos θ J/ψ acceptance, we fit the cos θ J/ψ distributions from simulation in 24 bins of m 2 (π + π − ) with the function giving 24 values of a as a function of m 2 (π + π − ). The resultant distribution shown in Fig. 9 can be described by an exponential function a(s 23 ) = exp(a 1 + a 2 s 23 ), with a 1 = −1.48 ± 0.20 and a 2 = (−1.45 ± 0.33) GeV −2 .    Figure 9: Exponential fit to the acceptance parameter a(s 12 ) used in Eq. 18.

Equation 18
is normalized with respect to cos θ J/ψ . Thus, after integrating over cos θ J/ψ , Eq. 17 becomes This term of the efficiency is parametrized as a symmetric fourth order polynomial function given by where the i are the fit parameters. Figure 10 shows the polynomial function obtained from a fit to the Dalitz-plot distributions of simulated events. The projections of the fit are shown in Fig. 11 and the resulting parameters are given in Table 1.   Candidates/ (0.1 GeV ) Figure 11: Projections onto (a) m 2 (J/ψ π + ) and (b) m 2 (π + π − ) of the simulated Dalitz plot used to determine the efficiency parameters. The points represent the simulated event distributions and the curves the projections of the polynomial fits.

Background composition
Backgrounds from B decays into J/ψ final states have already been discussed in Section 2.
The main background source is combinatorial and its shape can be determined from the same-sign π ± π ± combinations within ±20 MeV of the B 0 mass peak; this region also contains the small B − background. In addition, there is background arising from partially reconstructed B 0 , and a B 0 → J/ψ K − π + reflection, which cannot be present in same-sign combinations. We use simulated samples of inclusive B 0 s decays, and exclusive B 0 → J/ψ K * 0 (892) and B 0 → J/ψ K * 0 2 (1430) decays to model the additional backgrounds. The background fraction of each source is studied by fitting the J/ψ π + π − candidate invariant mass distributions in bins of m 2 (π + π − ). The resulting background distribution in the ±20 MeV B 0 signal region is shown in Fig. 12. It is fit by histograms from the same-sign combinations and two additional simulations, giving a partially reconstructed B 0 s background of 12.8%, and a reflection background that is 5.2% of the total background. LHCb π π Figure 12: The m 2 (ππ) distribution of background. The (black) histogram with error bars shows the same-sign data combinations with additional background from simulation, the (blue) points with error bars show the background obtained from the mass fits, the (black) dashed line is the partially reconstructed B 0 s background, and the (red) dotted is the misidentified B 0 → J/ψ K − π + contribution.
The background is parametrized as where the first part m(π + π − ) 2P R P B m B converts phase space from s 12 to cos θ ππ , and The variable ζ = 2(s 23 −s min )/(s max −s min )−1, where s min and s max give the fit boundaries, B 2 (ζ) is a fifth-order Chebychev polynomial with parameters b i (i = 1-5), and q(ζ) and p(ζ) are both second-order Chebychev polynomials with parameters c i (i=2, 3, 5, 6), and c 1 , and c 4 are free parameters. In order to better approximate the real background in the B 0 s signal region, the J/ψ π ± π ∓ candidates are kinematically constrained to the B 0 s mass. A fit to the same-sign sample, with additional background from simulation, determines b i , c i , m 0 and Γ 0 . Figure 13 shows the mass squared projections from the fit. The fitted background parameters are shown in Table 2.
The 1 + α cos 2 θ J/ψ term is a function of the J/ψ helicity angle. The cos θ J/ψ distribution of background is shown in Fig. 14, and is fit with the function 1 + α cos 2 θ J/ψ that determines the parameter α = −0.38 ± 0.04. We have verified that α is independent of s 23 .   Figure 14: distribution of the background in cos θ J/ψ resulting from J/ψ π + π − candidate mass fits in each bin of cos θ J/ψ . The curve represents the fitted function 1 + α cos 2 θ J/ψ .

Fit fractions
While a complete description of the decay is given in terms of the fitted amplitudes and phases, the knowledge of the contribution of each component can be summarized by defining a fit fraction, F R λ , as the integration of the squared amplitude of R over the Dalitz plot divided by the integration of the entire signal function, 0.42 ± 0.03 c 5 1.7 ± 0.8 c 6 2.5 ± 0.8 Note that the sum of the fit fractions over all λ and R is not necessarily unity due to the potential presence of interference between two resonances. If the Dalitz plot has more destructive interference than constructive interference, the total fit fraction will be greater than one. Interference term fractions are given by and the sum of the two is Note that interference terms between different spin-J states vanish, because the d J λ0 angular functions in A R λ are orthogonal. The statistical errors of the fit fractions depend on the statistical errors of every fitted magnitude and phase, and their correlations. Therefore, to determine the uncertainties the covariance matrix and parameter values from the fit are used to generate 500 sample parameter sets. For each set, the fit fractions are calculated. The distributions of the obtained fit fractions are described by bifurcated Gaussian functions. The widths of the Gaussians are taken as the statistical errors on the corresponding parameters. The correlations of fitted parameters are also taken into account.

Final state composition 4.1 Resonance models
To study the resonant structures of the decay B 0 → J/ψ π + π − we use those combinations with an invariant mass within ±20 MeV of the B 0 mass peak and apply a J/ψ K 0 S veto. The total number of remaining candidates is 8483, of which 3212 ± 80 are attributed to background. Possible resonances in the decay B 0 → J/ψ π + π − are listed in Table 3. In addition, there could be some contribution from non-resonant B 0 → J/ψ π + π − decays. Table 3: Possible resonances in the B 0 → J/ψ π + π − decay mode.
Resonance Spin Helicity Resonance formalism The masses and widths of the BW resonances are listed in Table 4. When used in the fit they are fixed to these values except for the parameters of the f 0 (500) resonance which are constrained by their uncertainties. Besides the mass and width, the Flatté resonance shape has two additional parameters g ππ and g KK , which are also fixed in the fit to values obtained in our previous Dalitz analysis of B 0 s → J/ψ π + π − [6], where a large fraction of B 0 s decays are to J/ψ f 0 (980). The parameters are taken to be m 0 = 939.9 ± 6.3 MeV, g ππ = 199 ± 30 MeV and g KK /g ππ = 3.0 ± 0.3. All background and efficiency parameters are fixed in the fit.
To determine the complex amplitudes in a specific model, the data are fitted maximizing the unbinned likelihood given as where N is the total number of candidates, and F is the total PDF defined in Eq. 1. The PDF is constructed from the signal fraction f sig , the efficiency model ε(s 12 , s 23 , θ J/ψ ), the background model B(s 12 , s 23 , θ J/ψ ), and the signal model S(s 12 , s 23 , θ J/ψ ). In order to ensure proper convergence using the maximum likelihood method, the PDF needs to be normalized. This is accomplished by first normalizing the J/ψ helicity dependent part ε(s 23 , θ J/ψ )Θ λ (θ J/ψ ) over cos θ J/ψ by analytical integration. This integration results in additional factors as a function of s 23 . We then normalize the mass dependent part multiplied by the additional factors using numerical integration over 500×500 bins. The fit determines the relative amplitude magnitudes a R i λ and phases φ R i λ defined in Eq. 3; we choose to fix a ρ(770) 0 to 1. As only relative phases are physically meaningful, one phase in each helicity grouping has to be fixed; we choose to fix those of the f 0 (500) and the ρ(770) (|λ| = 1) to 0. In addition, since the final state J/ψ π + π − is a self-charge-conjugate mode and as we do not determine the B flavor, the signal function is an average of B 0 and B 0 decays. If we do not consider π + π − partial waves of a higher order than D-wave, then we can express the differential decay rate derived from Eqs. 3, 4 and 8 in terms of S-, P-, and D-waves including helicity 0 and ±1 dΓ dm ππ d cos θ ππ d cos θ J/ψ for B 0 decays, where A s k λ and φ s k λ are the sum of amplitudes and reference phase for the spin-k resonance group, respectively. The B 0 function for decays is similar, but θ π + π − and θ J/ψ are changed to π − θ π + π − and π − θ J/ψ respectively, as a result of using π − and µ − to define the helicity angles, yielding Summing Eqs. 28 and 29 results in cancellation of the interference involving the λ = 0 terms for spin-1, and the λ = ±1 terms for spin-2, as they appear with opposite signs for B 0 and B 0 decays. Therefore, we have to fix one phase in spin-1 (λ = 0) group (φ s P 0 ) and one in spin-2 (λ = ±1) group (φ s D ±1 ); the phases of ρ(770) (λ = 0) and f 2 (1270) (λ = ±1) are fixed to zero. The other phases in each corresponding group are relative to that of the fixed resonance.

Fit results
To find the best model, we proceed by fitting with all the possible resonances and a nonresonance (NR) component, then subsequently remove the most insignificant component one at a time. We repeat this procedure until each remaining contribution has more than 3 statistical standard deviation (σ) significance. The significance is estimated from the fit fraction divided by its statistical uncertainty. The best fit model contains six resonances, the f 0 (500), f 0 (980), f 2 (1270), ρ(770), ρ(1450), and ω(782).
In order to compare the different models quantitatively an estimate of the goodness of fit is calculated from three-dimensional partitions of the one angular and two mass squared variables. We use the Poisson likelihood χ 2 [28] defined as where n i is the number of events in the three dimensional bin i and x i is the expected number of events in that bin according to the fitted likelihood function. A total of 1021 bins (N bin ) are used to calculate the χ 2 , based on the variables m 2 (J/ψ π + ), m 2 (π + π − ), and cos θ J/ψ . The χ 2 /ndf and the negative of the logarithm of the likelihood, −lnL, of the fits are given in Table 5; ndf is equal to N bin − 1 − N par , where N par is the number of fitting parameters. The difference between the best fit results and fits with one additional component is taken as a systematic uncertainty. Figure 15 shows the best fit model projections of m 2 (π + π − ), m 2 (J/ψ π + ), cos θ J/ψ and m(π + π − ). We calculate the fit fraction of each component using Eq. 24. For a P-or D-wave resonance, we report its total fit fraction by summing all the helicity components, and the fraction of the helicity λ = 0 component. The results are listed in Table 6. Systematic uncertainties will be discussed in Section 6. Two interesting ratios of fit fractions are (0.93 +0.37+0.47 −0.22−0.23 )% for ω(782) to ρ(770), and (9.5 +6.7 −3.4 ± 3.0)% for f 0 (980) to f 0 (500). The fit fractions of the interference terms are computed using Eq. 25 and listed in Table 7. Table 8 shows the resonant phases from the best fit. For the systematic uncertainty study, Table 9 shows the fit fractions of components for the best model with one additional resonance.   Figure 15: Dalitz fit projections of (a) m 2 (π + π − ), (b) m 2 (J/ψ π + ), (c) cos θ J/ψ and (d) m(π + π − ) for the best model. The points with error bars are data, the signal fit is shown with a (red) dashed line, the background with a (black) dotted line, and the (blue) solid line represents the total. In (a) and (d), the shape variations near the ρ(770) mass is due to ρ(770) − ω(782) interference, and the dip at the K 0 S mass [15] is due to the K 0 S veto.    Components 160 ± 48 f 0 (500) 0 (fixed)

Helicity angle distributions
We show the helicity angle distributions in the ρ(770) mass region defined within one full width of the ρ(770) resonance (the width values are given in Table 4) in Fig. 16. The cos θ J/ψ and cos θ ππ background subtracted and efficiency corrected distributions for this mass region are presented in Fig. 17. The distributions are in good agreement with the best fit model.

Branching fractions
Branching fractions are measured by normalizing to the well measured decay mode B − → J/ψ K − , which has two muons in the final state and has the same triggers as the B 0 → J/ψ π + π − decays. Assuming equal production of charged and neutral B mesons at the LHC due to isospin symmetry, the branching fraction is calculated as where N and denote the yield and total efficiency of the decay of interest. The branching fraction B(B − → J/ψ K − ) = (10.18 ± 0.42) × 10 −4 is determined from an average of recent Belle [29] and BaBar [30] measurements that are corrected with respect to the reported values, which assume equal production of charged and neutral B mesons at the Υ(4S), using the measured value of Γ(B + B − ) Γ(B 0 B 0 ) = 1.055 ± 0.025 [31]. Signal efficiencies are derived from simulations including trigger, reconstruction, and event selection components. Since the efficiency to detect the J/ψ π + π − final state is not uniform across the Dalitz plane, the efficiency is averaged according to the Dalitz model, where the best fit model is used. The K 0 S veto efficiency is also taken into account. Small corrections are applied to account for differences between the simulation and the data. We measure the kaon and pion identification efficiencies with respect to the simulation using D * + → π + D 0 (→ K − π + ) events selected from data. The efficiencies are measured in bins of p T and η and the averages are weighted using the signal event distributions in the data. Furthermore, to ensure that the p and p T distributions of the generated B mesons are correct we weight the B − and B 0 simulation samples using B − → J/ψ K − and B 0 → J/ψ K * 0 data, respectively. Finally, the simulation samples are weighted with the charged tracking efficiency ratio between data and simulation in bins of p and p T of the track. The average of the weights is the correction factor. The total correction factors are below 1.04 and largely cancel between the signal and normalization channels. Multiplying the simulation efficiencies and correction factors gives the total efficiency (1.163±0.003±0.017)% for B 0 → J/ψ π + π − and (3.092±0.012±0.038)% for B − → J/ψ K − , where the first uncertainty is statistical and the second is systematic.
Using N B − = 350,727 ± 633 and N B 0 = 5287 ± 112, we measure where the first uncertainty is statistical, the second is systematic and the third is due to the uncertainty of B(B − → J/ψ K − ). The systematic uncertainties are discussed in Section 6. Our measured value is consistent with and more precise than the previous BaBar measurement of (4.6 ± 0.7 ± 0.6) × 10 −5 [11]. Table 10 shows the branching fractions of resonant modes calculated by multiplying the fit fraction and the total branching fraction of B 0 → J/ψ π + π − . Since the f 0 (980) contribution has a significance of less than 3σ we quote also an upper limit of B B 0 → J/ψ f 0 (980) × B (f 0 (980) → π + π − ) < 1.1 × 10 −6 at 90% confidence level (CL); this is the first such limit. The limit is calculated assuming a Gaussian distribution as the central value plus 1.28 times the addition in quadrature of the statistical and systematic uncertainties. This branching ratio is predicted to be in the range (1 − 3) × 10 −6 if the f 0 (980) resonance is formed of tetra-quarks, but can be much smaller if the f 0 (980) is a standard quark anti-quark resonance [8]. Our limit is at the lower boundary of the tetra-quark prediction, and is consistent with a quark anti-quark resonance with a small mixing angle. In Section 7.2, we show that the mixing angle, describing the admixture of ss and light quarks, is less than 31 • at 90% CL.

Systematic uncertainties
The contributions to the systematic uncertainties on the branching fractions are listed in Table 11. Since the branching fractions are measured with respect to the B − → J/ψ K − mode, which has a different number of charged tracks than the decays of interest, a 1% systematic uncertainty is assigned due to differences in the tracking performance between data and simulation. Another 2% uncertainty is assigned because of the difference between two pions and one kaon in the final states, due to decay in flight, multiple scattering, Table 10: Branching fractions for each channel. The upper limit at 90% CL is also quoted for the f 0 (980) resonance which has a significance smaller than 3σ. The first uncertainty is statistical and the second the total systematic.
and hadronic interactions. Small uncertainties are introduced if the simulation does not have the correct B meson kinematic distributions. We are relatively insensitive to any differences in the B meson p and p T distributions since we are measuring the relative rates. By varying the p and p T distributions we see at most a change of 0.5%. There is a 1.0% systematic uncertainty assigned for the relative particle identification efficiencies (0.5% per particle). These efficiencies have been corrected from those predicted in the simulation by using the data from D * + → π + D 0 (→ K − π + ). A 0.6% uncertainty is included for the J/ψ π − π + efficiency, estimated by changing the best model to that including all possible resonances. The B 0 signal yield is changed by 0.5% when the shape of the combinatorial background is changed from an exponential to a linear function. The total systematic uncertainty is obtained by adding each source of systematic uncertainty in quadrature as they are uncorrelated. In addition, the largest source is 4.1% due to the uncertainty of B(B − → J/ψ K − ) which is quoted separately. The sources of the systematic uncertainties on the results of the Dalitz plot analysis are summarized in Table 12. For the uncertainties due to the acceptance or background modeling, we repeat the data fit 100 times where the parameters of acceptance or background modeling are generated according to the corresponding covariance matrix. We also study the acceptance function by changing the minimum IP χ 2 requirement from 9 to 12.5 on both of the pion candidates. As shown previously [6], this increases the χ 2 of the fit to the angular distributions by one unit. The acceptance function is then applied to the data with the original minimum IP χ 2 selection of 9, and the likelihood fit is redone and the uncertainties are estimated by comparing the results with the best fit model. The larger of the two variations is taken as uncertainty due to the acceptance.
We study the effect of ignoring the experimental mass resolution in the fit by comparing fits between different pseudo-experiments with and without the resolution included. As the widths of the resonances we consider are much larger than the mass resolution, we find that the effects are negligible except for the ω(782) resonance whose fit fraction is underestimated by (0.09 ± 0.08)%. Thus, we apply a 0.09% correction to the ω(782)

±3.0
fraction and assign an additional ±0.08% in the acceptance systematic uncertainty. The results shown in the previous sections already include this correction.
In the default fit, the signal fraction f sig = 0.621 ± 0.009, defined in Eq. 1 is fixed; we vary its value within its error to estimate the systematic uncertainty. The change is added in quadrature with the background modeling uncertainties.
The uncertainties due to the fit model include adding each resonance that is listed in Table 4 but not used in the best model, changing the default values of L B in P-and D-wave cases, varying the hadron scale r parameters for the B meson and R resonance to 3.0 GeV −1 for both, replacing the f 0 (500) model by a Zhou and Bugg function [33,34] and using the alternate Gounaris and Sakurai model [35] for ρ resonances. Then the largest variations among those changes are assigned as the systematic uncertainties for modeling (see Table 12).
Finally, we repeat the data fit by varying the mass and width of resonances (see Table 4) within their errors one at a time, and add the changes in quadrature.  [15], except for the f 0 (500) [27].
Isospin Vector particle Vector mass ( MeV) Scalar particle Scalar mass ( MeV) 0 Thus the f 0 (500) state is firmly established in B 0 → J/ψ π + π − decays. As discussed in the introduction, a region with only S-and P-waves is preferred for measuring sin 2β eff . The best fit model demonstrates that the mass region within ±149 MeV (one full width) of the ρ(770) mass contains only (0.72 ± 0.09)% D-wave contribution, thus this region can be used for a clean CP measurement. The S-wave in this region is (11.9±1.7)%, where the fraction is the sum of individual fit fractions and the interference.

7.2
Mixing angle between f 0 (980) and f 0 (500) The scalar nonet is quite an enigma. The mysteries are summarized in Ref. [36], and in the "Note on scalar mesons" in the PDG [15]. Let us contrast the masses of the lightest vector mesons with those of the scalars, listed in Table 13. For the vector particles, the ω and ρ masses are nearly degenerate and the masses increase as the s-quark content increases. For the scalar particles, however, the mass dependence differs in several ways which requires an explanation. Some authors introduce the concept of qqqq states or superpositions of the four-quark state with the qq state. In either case, the I = 0 f 0 (500) and the f 0 (980) are thought to be mixtures of the underlying states whose mixing angle has been estimated previously (see Ref. [8] and references contained therein).
Assuming that the ππ and KK decays are dominant we obtain B f 0 (980) → π + π − = (46 ± 6) %, where we have assumed that the only other decays are to π 0 π 0 , half of the π + π − rate, and to neutral kaons, taken equal to charged kaons. We use B (f 0 (500) → π + π − ) = 2 3 , which results from isospin Clebsch-Gordon coefficients, and assuming that the only decays are into two pions. Since we have only an upper limit on the J/ψ f 0 (980) final state, we will only find an upper limit on the mixing angle, so if any other decay modes of the f 0 (500) (f 0 (980)) exist, they would make the limit more (less) stringent. Our limit then is B B 0 → J/ψ f 0 (500) Φ(500) Φ(980) < 0.35 at 90% confidence level, which translates into a limit |ϕ m | < 31 • at 90% confidence level.

Conclusions
We have studied the resonance structure of B 0 → J/ψ π + π − using a modified Dalitz plot analysis where we also include the decay angle of the J/ψ meson. The decay distributions are formed from a series of final states described by individual π + π − interfering decay amplitudes. The largest component is the ρ(770) resonance. The data are best described by adding the f 2 (1270), f 0 (500), ω(782), ρ(1450) and f 0 (980) resonances, where the f 0 (980) resonance contributes less than 3σ significance. The results are listed in Table 6.
We set an upper limit B B 0 → J/ψ f 0 (980) × B (f 0 (980) → π + π − ) < 1.1 × 10 −6 at 90% confidence level that favors somewhat a quark anti-quark interpretation of the f 0 (980) resonance. We also have firmly established the existence of the J/ψ f 0 (500) intermediate resonant state in B 0 decays, and limit the absolute value of the mixing angle between the two lightest scalar states to be less than 31 • at 90% confidence level.
Our six-resonance best fit shows that the mass region within one full width of the ρ(770) contains mostly P-wave, (11.9 ± 1.7)% S-wave, and only (0.72 ± 0.09)% D-wave. Thus this region can be used to perform CP violation measurements, as the S-and P-wave components can be treated in the same manner as in the analysis of B 0 s → J/ψ φ [2][3][4]. The measured value of the asymmetry can be compared to that found in other modes such as B 0 → J/ψ K 0 in order to ascertain the possible effects due to penguin amplitudes.