Amplitude analysis and branching fraction measurement of Bs->J/\psi K+K-

An amplitude analysis of the final state structure in the Bs->J/\psi K+K- decay mode is performed using 1.0/fb of data collected by the LHCb experiment in 7 TeV center-of-mass energy pp collisions produced by the LHC. A modified Dalitz plot analysis of the final state is performed using both the invariant mass spectra and the decay angular distributions. Resonant structures are observed in the K+K- mass spectrum as well as a significant non-resonant S-wave contribution. The largest resonant component is the \phi(1020), accompanied by f0(980), f'2(1525), and four additional resonances. The overall branching fraction is measured to be B(Bs->J/\psi K+K-)=(7.70 +/-0.08 +/- 0.39 +/- 0.60)x 10^(-4), where the first uncertainty is statistical, the second systematic, and the third due to the ratio of the number of Bs to B- mesons produced. The mass and width of the f'2(1525) are measured to be 1522.2 +/- 2.8^{+5.3}_{-2.0} MeV and 84 +/- 6^{+10}_{-5} MeV, respectively. The final state fractions of the other resonant states are also reported.


Introduction
The study of B 0 s decays to J/ψh + h − , where h is either a pion or kaon, has been used to measure mixing-induced CP violation in B 0 s decays [1][2][3][4][5][6][7].‡ In order to best exploit these decays a better understanding of the final state composition is necessary.This study has been reported for the B 0 s → J/ψ π + π − channel [8].Here we perform a similar analysis for B 0 s → J/ψ K + K − .While a large φ(1020) contribution is well known [9] and the f 2 (1525) component has been recently observed [10] and confirmed [11], other components have not heretofore been identified including the source of S-wave contributions [12].The tree-level Feynman diagram for the process is shown in Fig. 1.In this paper the J/ψK + and K + K − mass spectra and decay angular distributions are used to study resonant and non-resonant structures.This differs from a classical "Dalitz plot" analysis [13] since the J/ψ meson has spin-1, and its three helicity amplitudes must be considered.

Data sample and detector
The event sample is obtained using 1.0 fb −1 of integrated luminosity collected with the LHCb detector [14] 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 momentum § resolution ∆p/p that varies from 0.4% at 5 GeV to 0.6% at 100 GeV.The impact parameter (IP) is defined as the minimum distance of approach of the track with respect to the primary vertex.For tracks with large transverse momentum with respect to the proton beam direction, the IP resolution is approximately 20 µm.Charged hadrons are identified using two ring-imaging Cherenkov detectors.Photon, electron and hadron candidates are ‡ Mention of a particular mode implies use of its charge conjugate throughout this paper.
§ We work in units where c = 1.
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 [15] 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.Events selected for this analysis are triggered by a J/ψ → µ + µ − decay, where the J/ψ is required at the software level to be consistent with coming from the decay of a B 0 s meson by use either of IP requirements or detachment of the J/ψ from the primary vertex.Monte Carlo simulations are performed using Pythia [16] with the specific tuning given in Ref. [17], and the LHCb detector description based on Geant4 [18] described in Ref. [19].Decays of B mesons are based on EvtGen [20].

Signal selection and backgrounds
We select B 0 s → J/ψK + K − candidates trying to simultaneously maximize the signal yield and reduce the background.Candidate J/ψ → µ + µ − decays are combined with a pair of kaon candidates of opposite charge, and then requiring that all four tracks are consistent with coming from a common decay point.To be considered a J/ψ → µ + µ − candidate, particles identified as muons of opposite charge are required to have transverse momentum, p T , greater than 500 MeV, and form a vertex with fit χ 2 per number of degrees of freedom (ndf) less than 11.These requirements give rise to a large J/ψ signal over a small background [21].Only candidates with a dimuon invariant mass between −48 MeV to +43 MeV relative to the observed J/ψ mass peak are selected.The asymmetric requirement is due to final-state electromagnetic radiation.The two muons are subsequently kinematically constrained to the known J/ψ mass [9].
Our ring-imaging Cherenkov system allows for the possibility of positively identifying kaon candidates.Charged tracks produce Cherenkov photons whose emission angles are compared with those expected for electrons, pions, kaons or protons, and a likelihood for each species is then computed.To identify a particular species, the difference between the logarithm of the likelihoods for two particle hypotheses (DLL) is computed.There are two criteria used: loose corresponds to DLL(K −π) > 0, while tight has DLL(K −π) > 10 and DLL(K − p) > −3.Unless stated otherwise, we require the tight criterion for kaon selection.
We select candidate K + K − combinations if each particle is inconsistent with having been produced at the primary vertex.For this test we require that the χ 2 formed by using the hypothesis that the IP is zero be greater than 9 for each track.Furthermore, each kaon must have p T > 250 MeV and the scalar sum of the p T of the kaon candidates must be greater than 900 MeV.To select B 0 s candidates we further require that the two kaon candidates form a vertex with χ 2 < 10, and that they form a candidate B 0 s vertex with the J/ψ where the vertex fit χ 2 /ndf < 5. We require that this B 0 s vertex be more 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.
The B 0 s candidate invariant mass distribution is shown in Fig. 2. The vertical lines indicate the signal and sideband regions, where the signal region extends to ±20 MeV around the nominal B 0 s mass [9] and the sidebands extend from 35 MeV to 60 MeV on either side of the peak.The small peak near 5280 MeV results from B 0 decays, and will be subject to future investigation.The background consist of combinations of tracks, which have a smooth mass shape through the J/ψ K + K − region, and peaking contributions caused by the reflection of specific decay modes where a pion is misidentified as a kaon.The reflection background that arises from the decay B 0 → J/ψK − π + , where the π + is misidentified as a K + , is determined from the number of B 0 candidates in the control region 25 − 200 MeV above the B 0 s mass peak.For each of the candidates in the J/ψK + K − control region, we reassign each of the two kaons in turn to the pion mass hypothesis.The resulting J/ψKπ invariant mass distribution is shown in Fig. 3.The peak at the B 0 mass has 906 ± 51 candidates, determined by fitting the data to a Gaussian function for the signal, and a polynomial function for the background.From these events we estimate the number in the B 0 s signal region, based on a simulation of the shape of the reflected distribution as a function of J/ψK − K + mass.Using simulated B 0 → J/ψK * 0 (892) and B 0 → J/ψK * 2 (1430) samples, we calculate 309 ± 17 reflection candidates within ±20 MeV of the B 0 s peak.This number is used as a constraint in the mass fit described below.
To determine the number of B 0 s signal candidates we perform a fit to the candidate J/ψK + K − invariant mass spectrum shown in Fig. 4. The fit function is the sum of the B 0 s signal component, combinatorial background, and the contribution from the B 0 → J/ψK − π + reflections.The signal is modeled by a double-Gaussian function with a common mean.The combinatorial background is described by a linear function.The   reflection background is constrained as described above.The mass fit gives 19,195±150 signal together with 894 ± 24 combinatorial background candidates within ±20 MeV of the B 0 s mass peak.We use the decay B − → J/ψ K − as the normalization channel for branching fraction determinations.The selection criteria are similar to those used for J/ψ K + K − , except for particle identification as here a loose kaon identification criterion is used.Figure 5 shows the J/ψK − mass distribution.The signal is fit with a double-Gaussian function and a linear function is used to fit the combinatorial background.There are 342,786±661 signal and 10,195±134 background candidates within ±20 MeV of the B − peak.

Analysis formalism
One of the goals of this analysis is to determine the intermediate states in B 0 s → J/ψK + K − decay within the context of an isobar model [22,23], where we sum the resonant and non-resonant components testing if they explain the invariant mass squared and angular distributions.We also determine the absolute branching fractions of B 0 s → J/ψ φ(1020) and B 0 s → J/ψ f 2 (1525) final states and the mass and width of the f 2 (1525) resonance.Another important goal is to understand the S-wave content in the φ(1020) mass region.
Four variables completely describe the decay of B 0 s → J/ψK + K − with J/ψ → µ + µ − .Two are the invariant mass squared of J/ψK + , s 12 ≡ m 2 (J/ψK + ), and the invariant mass squared of K + K − , s 23 ≡ m 2 (K + K − ).The other two are 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 K + K − decay planes, χ, in the B 0 s rest frame.To simplify the probability density function (PDF), we analyze the decay process after integrating over the angular variable χ, which eliminates several interference terms.

The model for
In order to perform an amplitude analysis a PDF must be constructed that models correctly the dynamical and kinematic properties of the decay.The PDF is separated into two components, one describing signal, S, and the other background, B. The overall PDF given by the sum is where ε is the detection efficiency.The background is described by the sum of combinatorial background, C, and reflection, R, functions where f com and f refl are the fractions of the combinatorial background and reflection, respectively, in the fitted region.The fractions f com and f refl obtained from the mass fit are fixed for the subsequent analysis.
The normalization factors are given by This formalism similar to that used by Belle in their analysis of B 0 → K − π + χ c1 [24], and later used by LHCb for the analysis of B 0 s → J/ψπ + π − [8].The invariant mass squared of J/ψK + versus K + K − is shown in Fig. 6 for B 0 s → J/ψK + K − candidates.No structure is seen in m 2 (J/ψ K + ).There are however visible horizontal bands in the K + K − mass squared spectrum, the most prominent of which correspond to the φ(1020) and f 2 (1525) resonances.These and other structures in m 2 (K + K − ) are now examined.
The signal function is given by the coherent sum over resonant states that decay into where A R i λ (s 12 , s 23 , θ J/ψ ) describes the decay amplitude via an intermediate resonance state R i with helicity λ.Note that the J/ψ has the same helicity as the intermediate K + K − resonance.Each R i has an associated amplitude strength a R i λ and a phase φ R i λ for each helicity state λ.The amplitude for resonance R, for each i, is given by where P R is the momentum of either of the two kaons in the di-kaon rest frame, m B is the B 0 s mass, P B is the magnitude of the J/ψ three-momentum in the B 0 s rest frame, and are the B 0 s meson and R i resonance decay form factors.The orbital angular momenta between the J/ψ and K + K − system is given by L B , and the orbital angular momentum in the K + K − decay is given by L R ; the latter is the same as the spin of the K + K − system.Since the parent B 0 s has spin-0 and the J/ψ is a vector, when the K + K − system forms a spin-0 resonance, L B = 1 and L R = 0.For K + K − 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 value and consider the other possibilities in the systematic uncertainty.
The Blatt-Weisskopf barrier factors and F (L R ) R [25] 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 [26].In both cases z 0 = r 2 P 2 0 where P 0 is the decay daughter momentum at the pole mass; for the B 0 s decay the J/ψ momentum is used, while for the R resonances the kaon momentum is used.
In the helicity formalism, the angular term, T λ (θ KK ) is defined as where d is the Wigner d-function, J is the resonance spin, θ KK is the helicity angle of the K + in the K + K − rest frame with respect to the K + K − direction in the B 0 s rest frame, and may be calculated directly from the other variables as The J/ψ helicity dependent term Θ λ (θ J/ψ ) is defined as The mass squared shape of each resonance, R is described by the function A R (s 23 ).In most cases this is a Breit-Wigner (BW) amplitude.When a decay channel opens close to the resonant mass, complications arise, since the proximity of the second threshold distorts the line shape of the amplitude.The f 0 (980) can decay to either ππ or KK.While the ππ channel opens at much lower masses, the K + K − decay channel opens near the resonance mass.Thus, for the f 0 (980) we use a Flatté model [27] that takes into account these coupled channels.
We describe the BW amplitude for a resonance decaying into two spin-0 particles, labeled as 2 and 3, as 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é mass shape is parametrized as where 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 For non-resonant processes, the amplitude A(s 12 , s 23 , θ J/ψ ) is constant over the variables s 12 and s 23 , but has an angular dependence due to the J/ψ decay.The amplitude is derived from Eq. ( 5), assuming that the non-resonant K + K − contribution is an S-wave (i.e.L R = 0, L B = 1) and is uniform in phase space (i.e.A R = 1),

Detection efficiency
The detection efficiency is determined from a phase space simulation sample containing 3.4 × 10 6 B 0 s → J/ψK + K − events with J/ψ → µ + µ − .We also use a separate sample of 1.3 × 10 6 B 0 s → J/ψφ events.The p and p T distributions of the generated B 0 s mesons are weighted to match the distributions found using J/ψφ data.The simulation is also corrected by weighting for difference between the simulated kaon detection efficiencies and the measured ones determined by using a sample of Next we describe the efficiency in terms of the analysis variables.Both s 12 and s 13 range from 12.5 GeV 2 to 24.0 GeV 2 , where s 13 is defined below, and thus are centered at s 0 = 18.25 GeV 2 .We model the detection efficiency using the dimensionless symmetric Dalitz plot observables and the angular variable θ J/ψ .The observables s 12 and s 13 are related to s 23 as To parametrize this efficiency, we fit the cos θ J/ψ distributions of the J/ψK + K − and J/ψφ simulation samples in bins of m 2 (K + K − ) with the function giving values of a as a function of m 2 (K + K − ).The resulting distribution, shown in Fig. 7, is described by an exponential function with a 1 = −0.76± 0.18 and a 2 = (−1.02± 0.15) GeV −2 .Equation ( 18) is normalized with respect to cos θ J/ψ .The efficiency in cos θ J/ψ depends on s 23 , and is observed to be independent of s 12 .Thus the detection efficiency can be expressed as After integrating over cos θ J/ψ , Eq. ( 20) becomes and is modeled by a symmetric fifth order polynomial function given by where i are the fit parameters.The B 0 s → J/ψK + K − phase space simulation sample is modeled with the polynomial function.The fitted function is shown in Fig. 8, and the projections of the fit are shown in Fig. 9.The efficiency is well described by the  parametrization.
For the region within ±20 MeV of the φ(1020) mass, the cos θ KK acceptance is used separately, due to the large number of signal events.Here the cos θ KK distribution shows a variation in efficiency, which can be parametrized using the efficiency function where the parameter 12 is measured from a fit to the simulated J/ψφ sample with ε 1 (x, y) × A(θ KK ), giving 12 = −0.099± 0.010, as shown in Fig. 10.The mass resolution is ∼ 0.7 MeV at the φ(1020) mass peak, which is added to the fit model by increasing the Breit-Wigner width of the φ(1020) to 4.59 MeV.

Background composition
The shape of the combinatorial background is modeled as where C 1 (s 12 , s 23 ) is parametrized as with c i , m 0 , Γ 0 and α as the fit parameters.The variables x and y are defined in Eq. ( 16).
Figure 11 shows the mass squared projections from the B 0 s mass sidebands with the fit projections overlaid.The χ 2 /ndf of the fit is 291/305.The value of α is determined by fitting the cos θ J/ψ distribution of background, as shown in Fig. 12, with a function of the form 1 + α cos 2 θ J/ψ , yielding α = −0.14 ± 0.08.
5 Final state composition

Resonance models
The resonances that are likely to contribute are produced from the ss system in Fig. 1, and thus are isoscalar (I = 0).The K + K − system in the decay B 0 s → J/ψK + K − can, in principle, have zero or any positive integer angular momentum.Both the Pparity and C-parity of K + K − pair in a state of relative angular momentum L are given by (−1) L .Therefore the allowed resonances decaying to K + K − are limited to J P C = 0 ++ , 1 −− , 2 ++ , ..., with isospin I = 0.In the kinematically accessible mass range up to 2 GeV, resonances with J P C = 3 −− or higher are not expected and thus the subsequent analysis only uses spins up to J = 2. Possible resonance candidates are listed in Table 1.
There could also be a contribution from non-resonant events which we assume to be S-wave and evenly distributed over the available phase space.To study the resonant structures of the decay B 0 s → J/ψK + K − we use 20, 425 candidates with an invariant mass within ±20 MeV of the observed B 0 s mass peak.This includes both signal and background, with 94% signal purity.We begin our analysis considering only the resonance components φ(1020), f 2 (1525) and a non-resonant component, established in our earlier measurement [10], and add resonances until no others are found with more than two standard deviation statistical significance (2σ).The significance is estimated from the fit fraction divided by its statistical uncertainty.Our best fit model includes a non-resonant component and 8 resonance states: φ(1020), f 0 (980), f 0 (1370), f 2 (1525), f 2 (1640) (|λ| = 1), φ(1680) (|λ| = 1), f 2 (1750), and f 2 (1950).Most of the resonances considered here are well established except for the modes f 2 (1640), f 2 (1750), and f 2 (1950).Although the existence of f 2 (1640) is not confirmed yet [9], the right shoulder The f 2 (1640) (λ = 0) and φ(1680) (λ = 0) components have less than two standard deviation significance when added separately to the fit, and therefore are not included in the best fit model. of f 2 (1525) fits better when we add this state.The presence of multiple broad overlapping resonances in this region may indicate a failure of the isobar model used in this analysis, but with the present data sample alternative descriptions are not feasible.Indeed, the situation is not clear for the resonance states in the vicinity of 1750 MeV.The PDG lists a spin-0 resonance, f 0 (1710), around 1.72 GeV of K + K − invariant mass [9].The Belle collaboration observed a resonance in the vicinity of 1.75 GeV with J P C = (even) ++ in their study of γγ → K + K − [28], but could not establish its spin.A state of mass 1767±14 MeV was seen by the L3 collaboration decaying into K 0 S K 0 S with J = 2 [29].We find that our data are better fit including the f 2 (1750) mode.If we substitute either the f 0 (1710) or f 0 (1750) resonance the fit is worsened, as the −lnL increase by 59 and 7 units, respectively.
In the same analysis of γγ → K + K − , Belle also observed the f 2 (1950) [28] resonance.We include this state in our best fit model.Furthermore, we do not expect significant contributions from the f 2 (1270) and f 0 (1500) resonances, since the PDG branching fractions are much larger in the π + π − final state than in K + K − [9] and we did not see significant contributions from these two resonances in the B 0 s → J/ψπ + π − final state [8].Therefore, these two resonances are not considered in the best fit model.However, we add these states, in turn, to the best fit model in order to test for their possible presence.
The masses and widths of the BW resonances are listed in Table 2.When used in the fit they are fixed to the central values, except for the f 2 (1525), whose mass and width are allowed to vary.
Belle [28] The f 0 (980) is described by a Flatté resonance shape, see Eq. ( 12).The parameters describing the function are the mass, and the couplings g ππ and g KK , which are fixed in the fit from the previous analysis of B 0 s → J/ψπ + π − [8].The parameters are m 0 = 939.9±6.3MeV, 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 maxi-mizing 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 normalization is accomplished by first normalizing the J/ψ helicity dependent part by analytical integration, and then for the mass dependent part using numerical integration over 400×800 bins.
The fit determines the relative values of the amplitude strengths, a R i λ , and phases, φ R i λ , defined in Eq. ( 4).We choose to fix a φ(1020) 0 = 1.As only relative phases are physically meaningful, one phase in each helicity grouping must be fixed.In addition, because J/ψK + K − is a self-charge-conjugate mode and does not determine the initial B flavor, the signal function is an average of B 0 s and B 0 s .If we consider no K + K − partial-waves of a higher order than D-wave, then we can express the differential decay rate (dΓ/dm KK d cos θ KK d cos θ J/ψ ) derived from Eq. ( 4) in terms of S-, P-, and D-waves including helicity 0 and ±1 components.The differential decay rates for B 0 s and B 0 s , respectively are where A s k λ and φ s k λ are the sum of amplitudes and reference phases, for the spin-k resonance group, respectively.The decay rate for B 0 s is similar to that of B 0 s , except θ K + K − and θ J/ψ are now changed to π − θ K + K − and π − θ J/ψ respectively, as a result of using K − and µ − to define the helicity angles and hence the signs change in front of the A s P 0 and A s D ±1 terms.Summing Eqs. ( 28) and ( 29) results in cancellation of the interference involving λ = 0 terms for spin-1, and the λ = ±1 terms for spin-2, as they appear with opposite signs for B 0 s and B 0 s decays.Therefore we have to fix one phase in the spin-1 (λ = 0) group (φ s P 0 ) and one in the spin-2 (λ = ±1) group (φ s D ±1 ).The other phases in each corresponding group are determined relative to that of the fixed resonance.

Fit results
The goodness of fit is calculated from 3D partitions of s 12 , s 23 and cos θ J/ψ .We use the Poisson likelihood χ 2 [30] defined as where n i is the number of candidates in the three dimensional bin i and x i is the expected number of candidates in that bin according to the fitted likelihood function.An adaptive binning algorithm is used, requiring a minimum of 25 entries in each bin.The associated number of degrees of freedom (ndf) is where k is the number of free parameters in the likelihood function.The χ 2 /ndf and the negative of the logarithm of the likelihood, −lnL, of the fits are given in Table 3. Starting values of parameters are varied in order to ensure that global likelihood minimums are found rather than local minimums.Attempts to add one more resonance such as f 2 (1270) and f 0 (1500) improve the −lnL marginally, but the χ 2 /ndf are worse than the best fit model.We retain only those resonances that are more than 2σ significant, except for the f 2 (1750) where we allow the |λ| = 1 component, since the λ = 0 component is significant.For models with one more resonance, the additional components never have more than 2σ significance.Figure 15 shows the projection of m 2 (K + K − ) for the best fit model, the m 2 (J/ψK + ) and cos θ J/ψ projections are displayed in Fig. 16.The projection of the K + K − invariant mass spectrum is shown in Fig. 17.
While a complete description of the B 0 s → J/ψ K + K − 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 If the Dalitz plot has more destructive interference than constructive interference, the sum of the fit fractions will be greater than unity.Conversely, the sum will be less than one if the Dalitz plot exhibits constructive interference.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 uncertainties of the fit fractions is difficult because they depend on the statistical uncertainty of every fitted magnitude and phase.Therefore we determine the uncertainties from simulated experiments.We perform 500 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 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 Gaussian functions are taken as the statistical uncertainties on the corresponding parameters.The fit fractions and phases of the contributing components are given in Table 4, while the fit fractions of the interference terms are quoted in Table 5.
Table 4: Fit fractions (%) and phases of contributing components.For P-and D-waves λ represents the helicity.
The off-diagonal elements give the fit fractions of the interference.The null values originate from the fact that any interference contribution between different spin-J state integrates to zero.Here the resonances are labeled by their masses in MeV and the subscripts denote the helicities.
9.9 ± 0.7 9.8 ± 0.7 9.5 ± 0.7 f 2 (1525), |λ| = 1 5.1 ± 0.9 5.1 ± 0.9 4.9 ± 0.9 the second largest contribution is the f 2 (1525), and the third the f 0 (980) resonance.There are also significant contributions from the f 0 (1370), f 2 (1640), φ(1680), f 2 (1750), f 2 (1950) resonances, and non-resonant final states.The amount of f 0 (980) is strongly parametrization dependent, so we treat these three models separately and do not assign any systematic uncertainty based on the use of these different f 0 (980) shapes.Therefore we refrain from quoting a branching fraction measurement for the decay B 0 s → J/ψf 0 (980).The determination of the parameters of the f 2 (1525) resonance are not dependent on the f 0 (980) parametrization.The parameters of the f 2 (1525) are determined to be: Whenever two or more uncertainties are quoted, the first is the statistical and the second systematic.The latter will be discussed in Section 5.6.These values are the most accurate determinations of the f 2 (1525) resonant parameters [9].Note that our determination of the mass has the same uncertainty as the current PDG average.

K + K − S-wave in the φ(1020) mass region
It was claimed by Stone and Zhang [12] that in the decay of B 0 s → J/ψφ, the K + K − system can have S-wave contributions under the φ(1020) peak of order 7% of the total yield.In order to investigate this possibility we calculate the S-wave fractions as given by the fit in 4 MeV mass intervals between 990 < m(K + K − ) < 1050 MeV.The resulting behavior is shown in Fig. 18.Here we show the result from our preferred model and also from the alternative f 0 (980) parameterizations discussed above.The observation of significant S-wave fractions in this region means that this contribution must be taken into account when measuring CP violation in the φ mass region.The total S-wave fraction as a function of the mass interval around the φ mass is also shown in Fig. 19.Using a time dependent analysis of B 0 s → J/ψφ(1020), LHCb reported (2.2 ± 1.2 ± 0.07)% [3] of S-wave within ±12 MeV of the φ(1020) mass peak.We measure the S-wave fraction within the same mass window as a consistent, and more precise (1.1 ± 0.1 +0.2 −0.1 )%.CDF measured the S-wave fraction as (0.8 ± 0.2)% for m(K + K − ) within about ±9.5 MeV of the φ mass [6], while ATLAS quotes (2±2)% for an 11 MeV interval [7].These results are consistent with ours.The D0 collaboration, however, claimed a (14.7±3.5)%S-wave fraction within approximately ±10 MeV of the φ meson mass [5], in disagreement with all of the other results.The squares (blue), triangles (red), and circles (green) represent the LHCb, BES and BaBar parameterizations of f 0 (980), respectively.The experimental statistical uncertainties are only shown for the LHCb model; they are almost identical for the other cases.The experimental mass resolution is not unfolded.

Helicity angle distributions
The decay angular distributions or the helicity angle distributions are already included in the signal model via Eqs.( 8) and (9).In order to test the fit model we examine the cos θ J/ψ and cos θ KK distributions in two different K + K − mass regions: one is the φ(1020) region defined within ±12 MeV of the φ(1020) mass peak and the other is defined within one full width of the f 2 (1525) mass.The background-subtracted efficiency-corrected distributions are shown in Figs.20 and 21.The distributions are in good agreement with the fit model.

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 [8,34].We define the angular moments Y 0 l as the efficiency-corrected and backgroundsubtracted K + K − invariant mass distributions, weighted by orthogonal and normalized spherical harmonic functions Y 0 l (cos θ KK ), If we assume that no K + K − partial-waves of a higher order than D-wave contribute, then we can express the differential decay rate, derived from Eq. ( 4) in terms of S-, P-, and D-waves including helicity 0 and ±1 components as where S λ , P λ , D λ and Φ k λ are real-valued functions of m KK , and we have factored out the S-wave phase.We can then calculate the angular moments The angular moments for l > 4 vanish.Figures 22 and 23 show the distributions of the angular moments for the fit model around ±30 MeV of the φ(1020) mass peak and above the φ(1020), respectively.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 sum of the interference between S-wave and P-wave, and P-wave and D-wave amplitudes, Y 0 2 the sum of 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 amplitudes, and Y 0 4 the D-wave.As discussed in Section 5.1, the average of B 0 s and B 0 s cancels the interference terms that involve P 0 and D ± .This causes the angular moments Y 0 1 and Y 0 3 to be zero when averaging over B 0 s and B 0 s decays.We observe that the fit results well describe the moment distributions, except for the Y 0 1 and Y 0 4 values below 1.2 GeV.This may be the result of statistical fluctuations or imperfect modeling.

Systematic uncertainties
The sources of the systematic uncertainties on the results of the Dalitz plot analysis are summarized in Table 7.The uncertainties due to the background parametrization are estimated by comparing the results from the best fit model with those when the background shape parameters are obtained from a fit to the lower sideband region only.The uncertainties in the efficiency are estimated by comparing the fit results when the efficiency parameters are changed by their statistical uncertainties and are added in quadrature.The effect on the fit fractions of changing the efficiency function is evaluated using a similar method to that used previously [8].Briefly, we change the efficiency model by increasing the minimum IP χ 2 requirement from 9 to 12.5 on both of the kaon candidates.This has the effect of increasing the χ 2 of the fit to the angular distributions of B 0 s → J/ψ φ data by 1 unit.The new efficiency function is then applied to the data with the original minimum IP χ 2 selection of 9, the likelihood is re-evaluated and the uncertainties are estimated by comparing the results with the best fit model.The largest variations among these two efficiency categories are included in the uncertainty.We estimate additional uncertainties by comparing the results when one more resonance is added to the best fit model.The uncertainties due to the line shape of the contributing resonances with fixed mass and width parameters are estimated by varying them individually in the fit according to their combined statistical and systematic uncertainties added in quadrature.We compare the results with the best fit and add them in quadrature to estimate the uncertainties due to the line shape.
Another source of systematic uncertainty is the value we choose for L B , the orbital angular momentum in the B 0 s decay.If L R equals zero then L B equals zero.If, however, L R is 1 then L B can either be 0 or 1, and if L R is 2, L B can be 1, 2 or 3.For our best fit we don't allow multiple values for L B , but choose the lowest allowed value.To estimate the systematic uncertainties due to the choice of L B , we repeat the fit changing the default value of L B , in turn, to each higher allowed value and compare the fit results with the best fit.The differences are grouped into the fit model category, and we assign the largest variations as the systematic uncertainties.These later two categories often give in asymmetric uncertainties.

Absolute branching fractions
Branching fractions are measured from ratios of the decay rates of interest normalized to the well established decay mode B − → J/ψK − .This decay mode, in addition to having a well measured branching fraction, has the advantage of having two muons in the final state and hence the same triggers as the B 0 s decay.However, we require knowledge of the B 0 s /B − production ratio.For this we assume isospin invariance and use the B 0 s /B 0 production ratio f s /f d = 0.256 ± 0.020, given in Ref. [35].The branching fractions are calculated using where X indicates a specific K + K − state, N represents the yield of the decay of interest, and corresponds to the overall efficiency.We form an average of B(B − → J/ψ K − ) = (10.18± 0.42) × 10 −4 using the recent Belle [36] and BaBar [37] measurements, corrected to take into account different rates of B + B − and B 0 B 0 pair production from Υ(4S) using The detection efficiency is obtained from simulation and is a product of the geometrical acceptance of the detector, the combined reconstruction and selection efficiency and the trigger efficiency.The efficiency also includes the efficiency of the Dalitz plot model for the case of B 0 s → J/ψK + K − , where the best fit model is used.The detection efficiencies and their various correction factors are given in Table 8.To ensure that the p and p T distributions of the generated B meson are correct we weight the B 0 s simulations using B 0 s → J/ψφ(1020) data and the B − simulations using B − → J/ψK − data.Since the control channel has a different number of charged tracks than the decay channel, we weight the simulations with the tracking efficiency ratio by comparing the data and simulations in bins of the track's p and p T .we further weight the B 0 s → J/ψ K + K − simulation, using the PDG value of B 0 s lifetime, (1.497 ± 0.015) × 10 −12 s [9], as input.where the branching fractions B(φ(1020) → K + K − ) = (48.9± 0.5)% and B(f 2 (1525) → K + K − ) = (44.4± 1.1)% are used [9].Here the first uncertainty in each case is statistical, the second is systematic and the third reflects the uncertainty due to f s /f d .Note that these are the time-integrated branching fractions.Results on the polarization fractions of B 0 s → J/ψ φ(1020) from a time-dependent analysis will be forthcoming in a separate publication [38].The ratio of B(B 0 s → J/ψf 2 (1525))/B(B 0 s → J/ψφ(1020)) is consistent with our previous result [10], D0 [11], and the Belle result [39].The current PDG value of B(B 0 s → J/ψφ(1020)) = (1.4 ± 0.5) × 10 −3 is dominated by the CDF measurement [40].Our measured value is in good agreement with this measurement and also the most recent yet unpublished values measured by CDF [41] and Belle [39].The Belle collaboration has also recently reported the branching fraction of B(B 0 s → J/ψK + K − ) [39], where B 0 s → J/ψφ(1020) is excluded.
The systematic uncertainty on the branching fraction has several contributions listed in Table 9.Since the branching fractions are measured with respect to B − → J/ψK − 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 Table 9: Relative systematic uncertainties on branching fractions (%).data and simulation.Another 2% uncertainty is assigned for the additional kaon which is due to decay in flight, large multiple scatterings and hadronic interactions along the track.Using the PDG value for the B 0 s lifetime [9] as input gives rise to an additional 1.5% systematic uncertainty.Small uncertainties are introduced if the simulation does not have the correct B meson kinematic distributions.We are relatively insensitive to any of these 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% systematic uncertainty assigned for the relative particle identification efficiencies.An uncertainty of 0.02% is included due to the change of the efficiency function Eq. (20).Three additional uncertainties are considered in the branching fractions of B(B 0 s → J/ψφ(1020)) and B(B 0 s → J/ψf 2 (1525)) as these are measured from the fit fractions of the Dalitz plot analysis.The total systematic uncertainty is obtained by adding each source of systematic uncertainty in quadrature as they are uncorrelated.

Conclusions
We have determined the final state composition of the B 0 s → J/ψK + K − decay channel using a modified Dalitz plot analysis where we include the decay angle of the J/ψ.The largest contribution is the φ(1020) resonance, along with other S-, P-and D-wave K + K − states, and a non-resonant K + K − contribution.All of the components are listed in Table 4.The mass and width of the where the first uncertainty in each case is statistical, the second is systematic and the third due to f s /f d .These results provide a good understanding of the J/ψ K + K − final state in B 0 s decays over the entire kinematically allowed region.The J/ψf 2 (1525) results supersede those of Ref. [10].This decay mode offers the opportunity for additional measurements of CP violation [42].

Figure 3 :
Figure 3: Invariant mass distribution for J/ψK + K − candidates 25 − 200 MeV above the B 0 smass, reinterpreted as B 0 → J/ψK ∓ π ± events.The fit is to a signal Gaussian whose mass and width are allowed to vary as well as the polynomial background.

Figure 4 :
Figure 4: Fit to the invariant mass spectrum of J/ψK + K − combinations.The dotted (black) line is the combinatorial background, the dashed (red) shape shows the misidentified B 0 → J/ψK − π + decays, and the solid (blue) curve shows the total.The vertical dashed lines indicate the signal region.

Figure 5 :
Figure 5: Fit to the invariant mass spectrum of J/ψK − candidates.The dotted line shows the combinatorial background and the solid (blue) curve is the total.

Figure 7 :
Figure 7: Exponential fit to the efficiency parameter a(s 23 ).The point near the φ(1020) meson mass is determined more precisely due to the use of a large simulation sample.

Figure 9 :
Figure 9: Projections of the invariant mass squared (a) K + K − and (b) J/ψK + from the simulation used to measure the efficiency parameters.The points represent the generated event distributions and the curves the polynomial fit.

Figure 11 :Figure 12 :Figure 13 :
Figure 11: Invariant mass squared projections of (a) K + K − and (b) J/ψK + from the background Dalitz plot of candidates in the B 0 s mass sidebands.

Figure 15 :
Figure 15: Dalitz plot fit projection of m 2 (K + K − ) using a logarithmic scale.The points with error bars are data, the (black) dotted curve shows the combinatorial background, the (red) dashed curve indicates the reflection from the misidentified B 0 → J/ψK − π + decays, and the (blue) solid line represents the total.

Figure 16 :Figure 17 :
Figure 16: Dalitz plot fit projections of (a) m 2 (J/ψK + ) and (b) cos θ J/ψ .The points with error bars are data, the (black) dotted curve shows the combinatorial background, the (red) dashed curve indicates the reflection from the misidentified B 0 → J/ψK − π + decays, and the (blue) solid line represents the total fit results.

Table 6 :
Comparison of the fit fractions (%) with the LHCb, BES and BaBar f 0 (980) parameterizations described in the text.For P-and D-waves, λ represents the helicity.

Figure 18 :Figure 19 :
Figure 18: S-wave fraction as a function of m(K + K − ) starting from 990 MeV up to 1050 MeV in 4 MeV mass intervals.The squares (blue), triangles (red), and circles (green) represent the LHCb, BES and BaBar parameterizations of f 0 (980), respectively.The experimental statistical uncertainties are only shown for the LHCb model; they are almost identical for the other cases.The experimental mass resolution is not unfolded.

Figure 22 :
Figure 22: Dependence of the spherical harmonic moments of cos θ KK as a function of the K + K − mass around the φ(1020) mass peak after efficiency corrections and background subtraction.The points with error bars are the data and the solid curves are derived from the fit model.

Figure 23 :
Figure 23: Dependence of the spherical harmonic moments of cos θ KK as a function of the K + K − mass above 1050 MeV, after efficiency corrections and background subtraction.The points with error bars are the data and the solid curves are derived from the fit model.

Table 1 :
Possible resonance candidates in the B 0 s → J/ψK + K − decay mode.

Table 7 :
Absolute systematic uncertainties on the fit results.

Table 8 :
Detector efficiencies determined from simulation and the correction factors.