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

The decay Bs->J/\psi\pi+\pi- can be exploited to study CP violation. A detailed understanding of its structure is imperative in order to optimize its usefulness. An analysis of this three-body final state is performed 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 of the final state is performed using both the invariant mass spectra and the decay angular distributions. The \pi+\pi- system is shown to be dominantly in an S-wave state, and the CP-odd fraction in this Bs decay is shown to be greater than 0.977 at 95% confidence level. In addition, we report the first measurement of the J/\psi\pi+\pi- branching fraction relative to J/\psi\phi\ of (19.79 +/- 0.47 +/- 0.52)%.


Introduction
Measurement of mixing-induced CP violation in B 0 s decays is of prime importance in probing physics beyond the Standard Model. Final states that are CP eigenstates with large rates and high detection efficiencies are very useful for such studies. The B 0 s → J/ψf 0 (980), f 0 (980) → π + π − decay mode, a CP -odd eigenstate, was discovered by the LHCb collaboration [1] and subsequently confirmed by several experiments [2]. As we use the J/ψ → µ + µ − decay, the final state has four charged tracks, and has high detection efficiency. LHCb has used this mode to measure the CP violating phase φ s [3], which complements measurements in the J/ψφ final state [4,5]. It is possible that a larger π + π − mass range could also be used for such studies. Therefore, to fully exploit the J/ψπ + π − final state for measuring CP violation, it is important to determine its resonant and CP content. The tree-level Feynman diagram for the process is shown in Fig. 1. In this paper the J/ψπ + and π + π − mass spectra, and decay angular distributions are used to study the resonant and non-resonant structures. This differs from a classical "Dalitz plot" analysis [6] because one of the particles in the final state, the J/ψ, 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 analysis requirements
The data sample contains 1.0 fb −1 of integrated luminosity collected with the LHCb detector [7] using pp collisions at a center-of-mass energy of 7 TeV. The detector is a single-arm 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 resolution ∆p/p that varies from 0.4% at 5 GeV to 0.6% at 100 GeV (we work in units where c = 1), and an impact parameter resolution of 20 µm for tracks with large transverse momentum 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 pre-shower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a muon 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 which applies a full event reconstruction.
Events selected for this analysis are triggered by a J/ψ → µ + µ − decay. Muon candidates are selected at the hardware level using their penetration through iron and detection in a series of tracking chambers. They are also required in the software level to be consistent with coming from the decay of a B 0 s meson into a J/ψ. Only J/ψ decays that are triggered on are used.

Selection requirements
The selection requirements discussed here are imposed to isolate B 0 s candidates with high signal yield and minimum background. This is accomplished by first selecting candidate J/ψ → µ + µ − decays, selecting a pair of pion candidates of opposite charge, and then testing if all four tracks form a common decay vertex. To be considered a J/ψ → µ + µ − candidate particles of opposite charge are required to have transverse momentum, p T , greater than 500 MeV, be identified as muons, and form a vertex with fit χ 2 per number of degrees of freedom (ndf) less than 11. After applying these requirements, there is a large J/ψ signal over a small background [1]. Only candidates with dimuon invariant mass between −48 MeV to +43 MeV relative to the observed J/ψ mass peak are selected. The requirement is asymmetric because of final state electromagnetic radiation. The two muons subsequently are kinematically constrained to the known J/ψ mass [8].
Pion and kaon candidates are positively identified using the RICH system. Cherenkov photons are matched to charged 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 is done by using the logarithm of the likelihood ratio comparing two particle hypotheses (DLL). For pion selection we require DLL(π − K) > −10.
Candidate π + π − combinations are selected if each particle is inconsistent with having been produced at the primary vertex. This is done by use of the impact parameter (IP) defined as the minimum distance of approach of the track with respect to the primary vertex. We require that the χ 2 formed by using the hypothesis that the IP is zero be greater than 9 for each track. Furthermore, each pion candidate must have p T > 250 MeV and the scalar sum of the two pion candidate momentum, p T (π + ) + p T (π − ), must The dashed line (barely visible along the x-axis) shows the S-wave contribution and the solid curve is the sum of the S-wave and a P-wave Breit-Wigner functions, fitted to the data. be greater than 900 MeV. To select B 0 s candidates we further require that the two pion candidates form a vertex with a χ 2 < 10, that they form a candidate B 0 s vertex with the J/ψ where the vertex fit χ 2 /ndf < 5, that this vertex is greater than 1.5 mm from the primary vertex and the angle between the B 0 s momentum vector and the vector from the primary vertex to the B 0 s vertex must be less than 11.8 mrad We use the decay B 0 s → J/ψφ, φ → K + K − as a normalization and control channel in this paper. The selection criteria are identical to the ones used for J/ψπ + π − except for the particle identification requirement. Kaon candidates are selected requiring that DLL(K − π) > 0. Figure 2(a) shows the J/ψK + K − mass for all events with m(K + K − ) < 1050 MeV. The K + K − combination is not, however, pure φ due to the presence of an S-wave contribution [9]. We determine the φ yield by fitting the data to a relativistic Pwave Breit-Wigner function that is convolved with a Gaussian function to account for the experimental mass resolution and a straight line for the S-wave. We use the S P lot method to subtract the background [10]. This involves fitting the J/ψK + K − mass spectrum, determining the signal and background weights and then plotting the resulting weighted mass spectrum, shown in Fig. 2(b). There is a large peak at the φ meson mass with a small S-wave component. The mass fit gives 20,934±150 events of which (95.5 ± 0.3)% are φ and the remainder is the S-wave contribution.
The invariant mass of the selected J/ψπ + π − combinations, where the dimuon candidate pair is constrained to have the J/ψ mass, is shown in Fig. 3. There is a large peak at the B 0 s mass and a smaller one at the B 0 mass on top of a background. A double-Gaussian function is used to fit the signal, the core Gaussian mean and width are allowed to vary, and the fraction and width ratio for the second Gaussian are fixed to that obtained in the  Figure 3: Invariant mass of J/ψπ + π − candidate combinations. The data have been fitted with double-Gaussian signal and several background functions. The (red) solid line shows the B 0 s signal, the (brown) dotted line shows the combinatorial background, the (green) short-dashed shows the B − background, the (purple) dot-dashed is B 0 → J/ψπ + π − , the (black) dot-long dashed is the sum of B 0 s → J/ψη and B 0 s → J/ψφ when φ → π + π − π 0 backgrounds, the (light blue) long-dashed is the B 0 → J/ψK − π + reflection, and the (blue) solid line is the total. fit of B 0 s → J/ψφ. Other components in the fit model take into account contributions from B − → J/ψK − (π − ), B 0 s → J/ψη , η → ργ, B 0 s → J/ψφ, φ → π + π − π 0 , B 0 → J/ψπ + π − backgrounds and a B 0 → J/ψK − π + reflection. Here and elsewhere charged conjugated modes are used when appropriate. The shape of the B 0 → J/ψπ + π − signal is taken to be the same as that of the B 0 s . The exponential combinatorial background shape is taken from wrong-sign combinations, that are the sum of π + π + and π − π − candidates. The shapes of the other components are taken from the Monte Carlo simulation with their normalizations allowed to vary (see Sect. 4.2). The mass fit gives 7598 ± 120 signal and 5825 ± 54 background candidates within ±20 MeV of the B 0 s mass peak.

Analysis formalism
The decay of B 0 s → J/ψπ + π − with the 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 (π + π − )), 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 s rest frame, and the angle between the J/ψ and π + π − decay planes (χ) in the B 0 s rest frame. To improve the resolution of these variables we perform a kinematic fit constraining the B 0 s and J/ψ masses to their PDG mass values [8], and recompute the final state momenta. To simplify the probability density function (PDF), we analyze the decay process after integrating over χ, that eliminates several interference terms. The χ distribution is shown in Fig. 4 after background subtraction using wrong-sign events. The distribution has little structure, and thus the χ acceptance can be integrated over without biasing the other variables.

The decay model for
One of the main challenges in performing a Dalitz plot angular analysis is to construct a realistic probability density function (PDF), where both the kinematic and dynamical properties are modeled accurately. The overall 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 normalization factors are given by In this analysis we apply a formalism similar to that used in Belle's analysis of B 0 → K − π + χ c1 decays [11].
To investigate if there are visible exotic structures in the J/ψπ + system as claimed in similar decays [12], we examine the J/ψπ + mass distribution shown in Fig Fig. 6 shows obvious structure in m 2 (π + π − ) that we wish to understand.

The signal function
The signal function is taken to be the sum over resonant states that can decay into π + π − , plus a possible non-resonant S-wave contribution 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 λ . The amplitudes are defined as where P B is the J/ψ momentum in the B 0 s 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 s mass, F are the B 0 s meson and R i resonance decay form factors, L B is the orbital angular momentum between the J/ψ and π + π − system, and L R the orbital angular momentum in the π + π − decay, and thus is the same as the spin of the π + π − . Since the parent B 0 s 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.
The Blatt-Weisskopf barrier factors F and F (L R ) R [13] are For the B meson z = r 2 P 2 B , where r, the hadron scale, is taken as 5.0 GeV −1 ; for the R resonance z = r 2 P 2 R , and r is taken as 1.5 GeV −1 . In both cases z 0 = r 2 P 2 0 where P 0 is the decay daughter momentum at the pole mass, different for the B 0 and the resonance decay.
The angular term, T λ , is obtained using the helicity formalism and is defined as where d is the Wigner d-function [8], J is the resonance spin, θ ππ is the π + π − resonance helicity angle which is defined as the angle of π + in the π + π − rest frame with respect to the π + π − direction in the B 0 s 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) because the K + K − decay channel opens. Here we use a Flatté model [14]. For non-resonant processes, the amplitude A R (s 23 ) is constant over the variables s 12 and s 23 , and has an angular dependence due to the J/ψ decay.
The BW amplitude for a resonance decaying into two spin-0 particles, labeled as 2 and 3, is where m R is the resonance 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 K + K − final states respectively. The ρ factors are given by Lorentz-invariant phase space The non-resonant amplitude is parametrized as

Detection efficiency
The detection efficiency is determined from a sample of one million B 0 s → J/ψπ + π − Monte Carlo (MC) events that are generated flat in phase space with J/ψ → µ + µ − , using Pythia [15] with a special LHCb parameter tune [16], and the LHCb detector simulation based on Geant4 [17] described in Ref [18]. After the final selections the MC has 78,470 signal events, reflecting an overall efficiency of 7.8%. The acceptance in cos θ J/ψ is uniform.
Next we describe the acceptance in terms of the mass squared variables. Both s 12 and s 13 range from 10.2 GeV 2 to 27.6 GeV 2 , where s 13 is defined below, and thus are centered at 18.9 GeV 2 . We model the detection efficiency using the symmetric Dalitz plot observables x = s 12 − 18.9 GeV 2 , and y = s 13 − 18.9 GeV 2 .
These variables are related to s 23 as The detection efficiency is parametrized as a symmetric 4th order polynomial function given by where the ε i are the fit parameters. The fitted polynomial function is shown in Fig. 7. The projections of the fit used to measure the efficiency parameters are shown in Fig. 8. The efficiency shapes are well described by the parametrization.
To check the detection efficiency we compare our simulated J/ψφ events with our measured J/ψφ helicity distributions. The events are generated in the same manner as for J/ψπ + π − . Here we use the measured helicity amplitudes of A || (0) 2 = 0.231 and |A 0 (0)| 2 = 0.524 [5]. The background subtracted J/ψφ angular distributions, cos θ J/ψ and cos θ KK , defined in the same manner as for the J/ψπ + π − decay, are compared in Fig. 9 with the MC simulation. The χ 2 /ndf =389/400 is determined by binning the angular distributions in two dimensions. The p-value is 64.1%. The excellent agreement gives us confidence that the simulation accurately predicts the acceptance.

Background composition
The main background source is taken from the wrong-sign combinations within ±20 MeV of the B 0 s mass peak. In addition, an extra 4.5% contribution from combinatorial background formed by J/ψ and random ρ(770), which cannot be present in wrong-sign combinations, is included using a MC sample. The level is determined by measuring the background yield as a function of π + π − mass. The background model is parametrized as where the first part B 1 (s 12 , s 23 ) is modeled using the technique of multiquadric radial basis functions [19]. These functions provide a useful way to parametrize multi-dimensional data giving sensible non-erratic behaviour and yet they follow significant variations in a smooth and faithful way. They are useful in this analysis in providing a modeling of the decay angular distributions in the resonance regions. Figure 10 shows the mass squared projections from the fit. The χ 2 /ndf of the fit is 182/145. We also used such functions with half the number of parameters and the changes were insignificant. The second part 1 + α cos θ J/ψ + β cos 2 θ J/ψ is a function of J/ψ helicity angle. The cos θ J/ψ distribution of background is shown in Fig. 11, fit with the function 1 + α cos θ J/ψ + β cos 2 θ J/ψ that determines the parameters α = −0.0050 ± 0.0201 and β = −0.2308 ± 0.0036.

Resonance models
To study the resonant structures of the decay B 0 s → J/ψπ + π − we use 13,424 candidates with invariant mass within ±20 MeV of the B 0 s mass peak. This includes both signal and background. Possible resonance candidates in the decay B 0 s → J/ψπ + π − are listed in Table 1. To understand what resonances are likely to contribute, it is important to realize that the ss system in Fig. 1 is isoscalar (I = 0) so when it produces a single meson it must have zero isospin, resulting in a symmetric isospin wavefunction for the two-pion system. Since the two-pions must be in an overall symmetric state, they must have even total angular momentum. In fact we only need to consider spin-0 and spin-2 particles as there are no known spin-4 particles in the kinematically accessible mass range below 1600 MeV. The particles that could appear are spin-0 f 0 (600), spin-0 f 0 (980), spin-2 f 2 (1270), spin-0 f 0 (1370) and spin-0 f 0 (1500). Diagrams of higher order than the one shown in Fig. 1 could result in the production of isospin-one π + π − resonances, thus we use the ρ(770) as a test of the presence of these higher order processes.   We proceed by fitting with a single f 0 (980), established from earlier measurements [1], and adding single resonant components until acceptable fits are found. Subsequently, we try the addition of other resonances. The models used are listed in Table 2.
The masses and widths of the BW resonances are listed in Table 3. When used in the fit they are fixed to these values, except for the f 0 (1370), for which they are not well measured, and thus are allowed to vary using their quoted errors as constraints in the fits, taking the errors as being Gaussian.
Besides the mass and width, the Flatté resonance shape has two additional parameters g ππ and g KK , which are also allowed to vary in the fit. Parameters of the non-resonant amplitude are also allowed to vary. One magnitude and one phase in each helicity grouping have to be fixed, since the overall normalization is related to the signal yield, and only relative phases are physically meaningful. The normalization and phase of f 0 (980) are fixed to 1 and 0 respectively. The phase of f 2 (1270), with helicity = ±1 is also fixed to Components zero when it is included. All background and efficiency parameters are held static 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 events, and F is the total PDF defined in Eq. 1. The PDF is constructed from the signal fraction f sig , efficiency model ε(s 12 , s 23 ), background model B(s 12 , s 23 , θ J/ψ ) and the signal model S(s 12 , s 23 , θ J/ψ ). The PDF needs to be normalized. This is accomplished by first normalizing the J/ψ helicity dependent part by analytical integration, and then for the mass dependent part using numerical integration over 500×500 bins.

Fit results
In order to compare the different models quantitatively an estimate of the goodness of fit is calculated from 3D partitions of the one angular and two mass-squared variables. We use the Poisson likelihood χ 2 [22] 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 N bin = 1356 bins are used to calculate the χ 2 , using 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 4. There are two solutions of almost equal likelihood for the 3R+NR model. Based on a detailed study of angular distributions (see Section 5.3) we choose one of these solutions and label it as "preferred". The other solution is called "alternate." We will use the differences between these to assign systematic uncertainties to the resonance fractions. The probability is improved noticeably adding components up to 3R+NR. Figure 12 shows the preferred model projections of m 2 (π + π − ) for the preferred model including only the 3R+NR components. The projections for the other considered models are indiscernible. The preferred model projections of m 2 (J/ψπ + ) and cos θ J/ψ are shown in Fig. 13 for the preferred model 3R+NR fit. The projections of the other preferred model fits including the additional resonances are almost identical. While a complete description of the decay is given in terms of the fitted amplitudes and phases, knowledge of the contribution of each component can be summarized by defining a fit fraction, F R λ . To determine F R λ we integrate the squared amplitude of R over the Dalitz plot. The yield is then normalized by integrating the entire signal function over the same area. Specifically, Note that the sum of the fit fractions is not necessarily unity due to the potential presence of interference between two resonances. Interference term fractions are given by and If the Dalitz plot has more destructive interference than constructive interference, the total fit fraction will be greater than one. Note that, interference between different spin-J states vanishes because the d J λ0 angular functions in A R λ are orthogonal. The determination of the statistical errors of the fit fractions is difficult because they depend on the statistical errors of every fitted magnitude and phase. A toy Monte Carlo approach is used. We perform 500 toy experiments: each sample is generated according to the model PDF, input parameters are taken from the fit to the data. The correlations of fitted parameters are also taken into account. For each toy experiment the fit fractions are calculated. The distributions of the obtained fit fractions are described by Gaussian functions. The r.m.s. widths of the Gaussians are taken as the statistical errors on the corresponding parameters. The fit fractions are listed in Table 5.
The 3R+NR fit describes the data well. For models adding more resonances, the additional components never have more than 3 standard deviation (σ) significance, and the fit likelihoods are only slightly improved. In the 3R+NR solution all the components have more than 3σ significance, except for the f 2 (1270) where we allow the helicity ±1 components since the helicity 0 component is significant. In all cases, we find the dominant contribution is S-wave which agrees with our previous less sophisticated analysis [3]. The D-wave contribution is small. The P-wave contribution is consistent with zero, as expected. The fit fractions from the alternate model are listed in Table 6. There are only small changes in the f 2 (1270) and ρ(770) components.
The fit fractions of the interference terms for the preferred and alternate models are computed using Eq. 22 and listed in Table 7.

Helicity distributions
Only S and D waves contribute to the B 0 s → J/ψπ + π − final state in the m(π + π − ) region below 1550 MeV. Helicity information is already included in the signal model via Eqs. 7  and 8. For a spin-0 π + π − system cos θ J/ψ should be distributed as 1−cos 2 θ J/ψ and cos θ ππ should be flat. To test our fits we examine the cos θ J/ψ and cos θ ππ distribution in different regions of π + π − mass. The decay rate with respect to the cosine of the helicity angles is where A 00 is the S-wave amplitude, A 2i , i = −1, 0, 1, the three D-wave amplitudes, and φ is the strong phase between A 00 and A 20 amplitudes. Non-flat distributions in cos θ ππ would indicate interference between the S-wave and D-wave amplitudes.
To investigate the angular structure we then split the helicity distributions into three different π + π − mass regions: one is the f 0 (980) region defined within ±90 MeV of the f 0 (980) mass and the others are defined within one full width of the f 2 (1270) and f 0 (1370) masses, respectively (the width values are given in Table 3). The cos θ J/ψ and cos θ ππ background-subtracted efficiency corrected distributions for these three different mass regions are presented in Figs. 14 and 15. The distributions are in good agreement with the 3R+NR preferred signal model. Furthermore, splitting into two bins, [−90, 0] and [0,90] MeV, we see different shapes, because across the pole mass of f 0 (980), the f 0 (980)'s phase changes by π. Hence the relative phase between f 0 (980) and the small D-wave in the two regions changes very sharply. This feature is reproduced well by the "preferred" model and shown in Fig. 16. The "alternate" model gives an acceptable, but poorer description.

Resonance parameters
The fit results from the four-component best fit are listed in Table 8 for both the preferred and alternate solutions. The table summarizes the f 0 (980) mass, the Flatté resonances parameters g ππ , g KK /g ππ , f 0 (1370) mass and width and the phases of the contributing resonances.
The mass and resonance parameters depend strongly on the final state in which they are measured, and the form of the resonance fitting function. Thus we do not quote systematic errors on these values. The value found for the f 0 (980) mass in the Flatté function 939.9 ± 6.3 MeV is lower than most determinations, although the observed peak value is close to 980 MeV, the estimated PDG value [8]. This is due to the interference from other resonances. The BES collaboration using the same functional form found a   [8]. Our result is within both of these ranges.

Angular moments
The angular moment distributions provide an additional way of visualizing the effects of different resonances and their interferences, similar to a partial wave analysis. This technique has been used in previous studies [24]. We define the angular moments Y 0 l as the efficiency corrected and background sub- 939.9 ± 6.3 939.2 ± 6.5 113 ± 11 108 ± 11 φ 980 0 (fixed) 0 (fixed) φ 1370 241.5 ± 6.3 181.7 ± 8.4 φ NR 217.0 ± 3.7 232.2 ± 3.7 φ 1270 , λ = 0 165 ± 15 118 ± 15 φ 1270 , |λ| = 1 0 (fixed) 0 (fixed) tracted π + π − invariant mass distributions, weighted by spherical harmonic functions The spherical harmonic functions satisfy If we assume that no π + π − partial-waves of a higher order than D-wave contribute, then we can express the differential decay rate (dΓ) derived from Eq. (3) in terms of S-, P-, and D-waves including helcity 0 and ±1 components as where A k λ and φ k λ are real-valued functions of m ππ , and we have factored out the S-wave phase. We then calculate the angular moments √ Figure 17 shows the distributions of the angular moments for the preferred solution. In general the interpretation of these moments is that Y 0 0 is the efficiency corrected and background subtracted event distribution, Y 0 1 the interference of the sum of S-wave and P-wave and P-wave and D-wave amplitudes, Y 0 2 the sum of the P-wave, D-wave and the interference of S-wave and D-wave amplitudes, Y 0 3 the interference between P-wave and D-wave, and Y 0 4 the D-wave. In our data the Y 0 1 distribution is consistent with zero, confirming the absence of any P-wave. We do observe the effects of the f 2 (1270) in the Y 0 2 distribution including the interferences with the S-waves. The other moments are consistent with the absence of any structure, as expected.

CP content
The main result in this paper is that CP -odd final states dominate. The f 2 (1270) helicity ±1 yield is (0.21 ± 0.65)%. As this represents a mixed CP state, the upper limit on the CP -even fraction due to this state is < 1.3 % at 95% confidence level (CL). Adding the ρ(770) amplitude and repeating the fit shows that only an insignificant amount of ρ(770) can be tolerated; in fact, the isospin violating J/ψρ(770) final state is limited to < 1.5% at 95% CL. The sum of f 2 (1270) helicity ±1 and ρ(770) is limited to < 2.3% at 95% CL. In the π + π − mass region within ±90 MeV of 980 MeV, this limit improves to < 0.6% at 95% CL.
The points with error bars are the data points and the solid curves are derived from the 3R+NR preferred model.

Total branching fraction ratio
To avoid the uncertainties associated with absolute branching fraction measurements, we quote branching fractions relative to the B 0 s → J/ψφ channel. The detection efficiency for this channel from Monte Carlo simulation is (1.07 ± 0.01)%, where the error is due to the limited Monte Carlo sample size.
The simulated detection efficiency for B 0 s → J/ψπ + π − as a function of the m 2 (π + π − ) is shown in Fig. 18. The simulation does not model the pion and kaon identification efficiencies with sufficient accuracy for our purposes. Therefore, we measure the kaon identification efficiency with respect to the Monte Carlo simulation. We use samples of D * + → π + D 0 , D 0 → K − π + events selected without kaon identification to measure the kaon and pion efficiencies with respect to the simulation, and an additional sample of K 0 s → π + π − decay for pions. The identification efficiency is measured in bins of p T and η and then the averages are weighted using the event distributions in the data. We find the correction to the J/ψφ efficiency is 0.970 (two kaons) and that to the J/ψf 0 efficiency is 0.973 (two pions). The additional correction due to particle identification then is 0.997±0.010. In addition, we re-weight the B 0 s p and p T distributions in the simulation which lowers the π + π − efficiency by 1.01% with respect to the K + K − efficiency.
Dividing the number of the J/ψπ + π − signal events by the J/ψK + K − yield, applying the additional corrections as described above, and taking into account B (φ → K + K − ) = (48.9 ± 0.5)% [8], we find Whenever two uncertainties are quoted the first is statistical and the second systematic. The latter will be discussed later in Section 7. This branching fraction ratio has not been previously measured.

Relative resonance yields
Next we evaluate the relative yields for the 3R+NR fit to the J/ψπ + π − final state from the preferred solution. We normalize the individual fit fractions reported in Table 5 by the sum. These normalized fit fractions are listed in Table 9 along with the branching fraction relative to J/ψφ, φ → K + K − , defined as R r , where r refers to the particular final state under consideration. Thus We use the difference between the preferred and alternate solutions found for the 3R+NR fit to assign a systematic uncertainty. Other systematic uncertainties are described in Section 7. The value found for R r for the f 0 (980), 0.139 ± 0.006 +0.025 −0.012 , is consistent with the prediction of Ref. [9], and consistent with the our first observation using 33 pb −1 of integrated luminosity [1], after multiplying by B (φ → K + K − ). The decay B 0 s → J/ψf 0 (1370) is now established. Previously both LHCb [1] and Belle [2] had seen evidence for this final state. The normalized f 2 (1270) helicity zero rate is (0.49±0.16)% in the preferred model and (0.42±0.11)% for the alternate solution.

Systematic uncertainties
Systematic uncertainties on the CP -odd fraction are negligible. In fact, using any of the alternate fits with different additional components does not introduce any significant fractions of CP -odd final states.
The systematic uncertainties on the branching fraction ratios have several contributions listed in Table 10. Since R r is measured relative to J/ψφ there is no systematic uncertainty due to differences in the tracking performance between data and simulation. The J/ψφ P-wave yield is fully correlated with the S-wave yield whose uncertainty we estimate as 0.7% by changing the signal PDF, and the background shape. By far the largest uncertainty in every rate, except the total, is caused by our choice of the preferred versus the alternate solutions. Using the difference between these fit results for the systematic uncertainty causes relatively large and asymmetric values. We also include systematic uncertainties due to the possible presence of the ρ(770), the f 0 (1500), or the f 0 (600) resonances by taking the maximum difference between the fit including one of these resonances and our preferred solution, if the difference is larger than the one between the preferred and alternate 3R+NR fit. In the case of the f 0 (1500) the preferred solution is pathological in that it produces an unacceptably large f 0 (1370) component along with a 214% component sum; therefore here we use the alternate solution that is much better behaved.
The uncertainty from Monte Carlo sample size for the mass dependent π + π − efficiencies are accounted for in the statistical errors, a residual systematic uncertainty is included that results from allowed changes in the shape due to the distribution of the events. The size of these differences depends on the mass range for the particular component multiplied by the possible efficiency variation across this mass range. This is estimated as 1% for the entire mass range and is smaller for individual resonances. Small uncertainties are introduced if the simulation does not have the correct B 0 s kinematic distributions. We are relatively insensitive to any these differences in the B 0 s p and p T distributions since we are measuring relative rates. These distributions are varied by changing the weights in each bin by plus and minus the statistical error in that bin. We see at most a 0.5% change. There is a 2% systematic uncertainty assigned for the relative particle identification efficiencies. These efficiencies have been corrected from those predicted in the simulation by using pion data from K 0 s → π + π − decays and kaon and pion data from D * ± → π ± D 0 (D 0 ), D 0 (D 0 ) → K ∓ π ± decays. The uncertainty on the corrections is 0.5% per track. The background modeling was changed by using a second-order polynomial shape in the J/ψπ + π − mass fit giving a 0.6% change in the signal yield. Since the input f 0 (1370) mass and width parameters were allowed to vary within Gaussian constraints, there is no additional uncertainty to account for. The effect on the fit fractions of changing the acceptance function is also evaluated. Since the acceptance model was tested by its agreement with the B 0 s → J/ψK + K − data in Fig. 9, we vary the data so that the model does not fit as well. This is accomplished by increasing the minimum IP χ 2 requirement from 9 to 12.25 on both of the kaon candidates, which has the effect of increasing the χ 2 /ndf of the fit to angular distributions by 1 unit. The Monte Carlo simulation of B 0 s → J/ψπ + π − with the changed requirement is then fitted to get an acceptance function. This acceptance function is then applied to the data with the original minimum IP χ 2 cut of 9, and the likelihood fit is redone. The resulting fitted values from the preferred solution are compared with the original values in Table 11.
The changes are small and well within the statistical uncertainties.

Conclusions
We have studied the resonance structure of B 0 s → J/ψπ + π − using a modified Dalitz plot analysis where we also include the decay angle of the J/ψ. The decay distributions are formed from a series of final states described by individual π + π − interfering decay amplitudes. The largest component is the f 0 (980) that is described by a Flatté function. The data are best described by adding Breit-Wigner amplitudes for the f 0 (1370), the f 2 (1270) resonances and a non-resonance contribution. Adding a ρ(770) into the fit does not improve the overall likelihood. Inclusion of f 0 (600) or f 0 (1500) does not result in significant signals for these resonances.
Our three resonance plus non-resonance best fit is dominantly CP -odd S-wave over the entire signal region. We also have a D-wave component arising from the f 2 (1270) resonance. Part of this corresponds to the A 20 amplitude which is also pure CP -odd and is (0.49 ± 0.16 +0.02 −0.08 )% of the total rate. A mixed CP part corresponding to the A 2±1 amplitude is (0.2 ± 0.7)% of the total. Adding this to the amount of allowed ρ(770), less than 1.5% at 95% CL, we find that the CP -odd fraction is greater than 0.977 at 95% CL. Thus, the entire mass range can be used to study CP violation with this almost pure CP -odd final state.
The measured relative branching ratio is where the first uncertainty is statistical and the second systematic. The largest component is the f 0 (980) resonance. We also determine This state was predicted to exist and have a branching fraction about 10% that of J/ψφ [9].
Our new measurement is consistent with and somewhat larger than this prediction. Other models give somewhat higher rates [25]. We also have firmly established the existence of the J/ψf 0 (1370) final state in B 0 s decay.