The $e^+ e^-\to K^+ K^- \pi^+\pi^-$, $K^+ K^- \pi^0\pi^0$ and $K^+ K^- K^+ K^-$ Cross Sections Measured with Initial-State Radiation

We study the processes $e^+ e^-\to K^+ K^- \pi^+\pi^-\gamma$, $K^+K^-\pi^0\pi^0\gamma$ and $K^+ K^- K^+ K^-\gamma$, where the photon is radiated from the initial state. About 34600, 4400 and 2300 fully reconstructed events, respectively, are selected from 232 \invfb of \babar data. The invariant mass of the hadronic final state defines the effective \epem center-of-mass energy, so that the $K^+ K^- \pi^+\pi^-\gamma$ data can be compared with direct measurements of the $e^+ e^-\to K^+K^- \pipi$ reaction; no direct measurements exist for the $e^+ e^-\to K^+ K^- \pi^0\pi^0$ or $\epem\to K^+ K^- K^+ K^-$ reactions. Studying the structure of these events, we find contributions from a number of intermediate states, and we extract their cross sections where possible. In particular, we isolate the contribution from $e^+ e^-\to\phi(1020) f_{0}(980)$ and study its structure near threshold. In the charmonium region, we observe the $J/\psi$ in all three final states and several intermediate states, as well as the $\psi(2S)$ in some modes, and measure the corresponding branching fractions. We see no signal for the Y(4260) and obtain an upper limit of $\BR_{Y(4260)\to\phi\pi^+\pi^-}\cdot\Gamma^{Y}_{ee}<0.4 \ev$ at 90% C.L.


I. INTRODUCTION
Electron-positron annihilation at fixed center-of-mass (c.m.) energies has long been a mainstay of research in elementary particle physics. The idea of utilizing initialstate radiation (ISR) to explore e + e − reactions below the nominal c.m. energies was outlined in Ref. [1], and discussed in the context of high-luminosity φ and B factories in Refs. [2][3][4]. At high energies, e + e − annihilation is dominated by quark-level processes producing two or more hadronic jets. However, low-multiplicity exclusive processes dominate at energies below about 2 GeV, and the region near charm threshold, 3.0-4.5 GeV, features a number of resonances [5]. These allow us to probe a wealth of physics parameters, including cross sections, spectroscopy and form factors.
Of particular current interest are the recently observed states in the charmonium region, such as the Y (4260) [6], and a possible discrepancy between the measured value of the anomalous magnetic moment of the muon, g µ − 2, and that predicted by the Standard Model [7]. Charmonium and other states with J P C = 1 −− can be observed as resonances in the cross section, and intermediate states may be present in the hadronic system. Measurements of the decay modes and their branching fractions are important in understanding the nature of these states. For example, the glue-ball model [8] predicts a large branching fraction for Y (4260) into φππ. The prediction for g µ − 2 is based on hadronic-loop corrections measured from low-energy e + e − → hadrons data, and these dominate the uncertainty on the prediction. Improving this prediction requires not only more precise measurements, but also measurements over the entire energy range and inclusion of all the important subprocesses in order to understand possible acceptance effects. ISR events at B factories provide independent and contiguous measurements of hadronic cross sections from the production threshold to about 5 GeV.
The cross section for the radiation of a photon of energy E γ followed by the production of a particular hadronic final state f is related to the corresponding direct e + e − → f cross section σ f (s) by where √ s is the initial e + e − c.m. energy, x = 2E γ / √ s is the fractional energy of the ISR photon and E c.m. ≡ * Deceased † Also with Università di Perugia, Dipartimento di Fisica, Perugia, Italy ‡ Also with Università della Basilicata, Potenza, Italy § Also with IPPP, Physics Department, Durham University, Durham DH1 3LE, United Kingdom s(1 − x) is the effective c.m. energy at which the final state f is produced. The probability density function W (s, x) for ISR photon emission has been calculated with better than 1% precision (see e.g. Ref. [4]). It falls rapidly as E γ increases from zero, but has a long tail, which combines with the increasing σ f (s(1−x)) to produce a sizable cross section at very low E c.m. . The angular distribution of the ISR photon peaks along the beam directions, but 10-15% [4] of the photons are within a typical detector acceptance.
Experimentally, the measured invariant mass of the hadronic final state defines E c.m. . An important feature of ISR data is that a wide range of energies is scanned simultaneously in one experiment, so that no structure is missed and the relative normalization uncertainties in data from different experiments or accelerator parameters are avoided. Furthermore, for large values of x the hadronic system is collimated, reducing acceptance issues and allowing measurements at energies down to production threshold. The mass resolution is not as good as a typical beam energy spread used in direct measurements, but the resolution and absolute energy scale can be monitored by the width and mass of well known resonances, such as the J/ψ produced in the reaction e + e − → J/ψγ. Backgrounds from e + e − → hadrons events at the nominal √ s and from other ISR processes can be suppressed by a combination of particle identification and kinematic fitting techniques. Studies of e + e − → µ + µ − γ and several multi-hadron ISR processes using BABAR data have been reported [9][10][11][12], demonstrating the viability of such measurements.
The K + K − π + π − final state has been measured directly by the DM1 collaboration [13] for √ s < 2.2 GeV, and we have previously published ISR measurements of the K + K − π + π − and K + K − K + K − final states [11] for E c.m. < 4.5 GeV. We recently reported [14] an updated measurement of the K + K − π + π − final state with a larger data sample, along with the first measurement of the K + K − π 0 π 0 final state, in which we observed a structure near threshold in the φf 0 intermediate state. In this paper we present a more detailed study of these two final states along with an updated measurement of the K + K − K + K − final state. In all cases we require detection of the ISR photon and perform a set of kinematic fits. We are able to suppress backgrounds sufficiently to study these final states from their respective production thresholds up to 5 GeV. In addition to measuring the overall cross sections, we study the internal structure of the events and measure cross sections for a number of intermediate states. We study the charmonium region, measure several J/ψ and ψ(2S) branching fractions, and set limits on other states.

II. THE BABAR DETECTOR AND DATASET
The data used in this analysis were collected with the BABAR detector at the PEP-II asymmetric energy e + e − storage rings. The total integrated luminosity used is 232 fb −1 , which includes 211 fb −1 collected at the Υ (4S) peak, √ s = 10.58 GeV, and 21 fb −1 collected below the resonance, at √ s = 10.54 GeV.
The BABAR detector is described elsewhere [15]. Here we use charged particles reconstructed in the tracking system, which comprises the five-layer silicon vertex tracker (SVT) and the 40-layer drift chamber (DCH) in a 1.5 T axial magnetic field. Separation of charged pions, kaons and protons uses a combination of Cherenkov angles measured in the detector of internally reflected Cherenkov light (DIRC) and specific ionization measured in the SVT and DCH. For the present study we use a kaon identification algorithm that provides 90-95% efficiency, depending on momentum, and pion and proton rejection factors in the 20-100 range. Photon and electron energies are measured in the CsI(Tl) electromagnetic calorimeter (EMC). We use muon identification provided by the instrumented flux return (IFR) to select the µ + µ − γ final state.
To study the detector acceptance and efficiency, we use a simulation package developed for radiative processes. The simulation of hadronic final states, including K + K − π + π − γ, K + K − π 0 π 0 γ and K + K − K + K − γ, is based on the approach suggested by Czyż and Kühn [16]. Multiple soft-photon emission from the initialstate charged particles is implemented with a structurefunction technique [17,18], and photon radiation from the final-state particles is simulated by the PHOTOS package [19]. The accuracy of the radiative corrections is about 1%.
We simulate the K + K − ππ final states both according to phase space and with models that include the φ(1020) → K + K − and/or f 0 (980) → ππ channels, and the K + K − K + K − final state both according to phase space and including the φ → K + K − channel. The generated events go through a detailed detector simulation [20], and we reconstruct them with the same software chain as the experimental data. Variations in detector and background conditions are taken into account.
We also generate a large number of background processes, including the ISR channels e + e − → π + π − π + π − γ and π + π − π 0 π 0 γ, which can contribute due to particle misidentification, and φηγ, φπ 0 γ, π + π − π 0 γ, which have larger cross sections and can contribute via missing or spurious tracks or photons. In addition, we study the non-ISR backgrounds e + e − → qq (q = u, d, s, c) generated by JETSET [21] and e + e − → τ + τ − by KORALB [22]. The contribution from the Υ (4S) decays is found to be negligible. The cross sections for these processes are known with about 10% accuracy or better, which is sufficient for these measurements.

III. EVENT SELECTION AND KINEMATIC FIT
In the initial selection of candidate events, we consider photon candidates in the EMC with energy above 0.03 GeV and charged tracks reconstructed in the DCH or SVT or both that extrapolate within 0.25 cm of the beam axis in the transverse plane and within 3 cm of the nominal collision point along the axis. These criteria are looser than in our previous analysis [11], and have been chosen to maximize efficiency. We require a highenergy photon in the event with an energy in the initial e + e − c.m. frame of E γ > 3 GeV, and either exactly four charged tracks with zero net charge and total momentum roughly opposite to the photon direction, or exactly two oppositely charged tracks that combine with a set of other photons to roughly balance the highest-energy photon momentum. We fit a vertex to the set of charged tracks and use it as the point of origin to calculate the photon direction. Most events contain additional soft photons due to machine background or interactions in the detector material.
We subject each of these candidate events to a set of constrained kinematic fits, and use the fit results, along with charged-particle identification, both to select the final states of interest and to measure backgrounds from other processes. We assume the photon with the highest E γ in the c.m. frame is the ISR photon, and the kinematic fits use its direction along with the four-momenta and covariance matrices of the initial e + e − and the set of selected tracks and photons. Because of excellent resolution for the momenta in the DCH and good angular resolution for the photons in the EMC, the ISR photon energy is determined with better resolution through fourmomentum conservation than through measurement in the EMC. Therefore we do not use its measured energy in the fits, eliminating the systematic uncertainty due to the EMC calibration for high energy photons. The fitted three-momenta for each charged track and photon are used in further kinematical calculations.
For the four-track candidates, the fits have three constraints (3C). We first fit to the π + π − π + π − hypothesis, obtaining a χ 2 4π . If the four tracks include one identified K + and one K − , we fit to the K + K − π + π − hypothesis and retain the event as a K + K − π + π − candidate. For events with one identified kaon, we perform fits with each of the two oppositely charged tracks given the kaon hypothesis, and the combination with the lower χ 2 4π . If the event contains three or four identified K ± , we fit to the K + K − K + K − hypothesis and retain the event as a K + K − K + K − candidate.
For the events with two charged tracks and five or more photon candidates, we require both tracks to be identified as kaons to suppress background from ISR π + π − π 0 π 0 and K ± K 0 S π ∓ events. We then pair all non-ISR photon candidates and consider combinations with invariant mass within ±30 MeV/c 2 of the π 0 mass as π 0 candidates. We perform a six-constraint (6C) fit to each set of Events/unit χ 2 FIG. 1: Distribution of χ 2 from the three-constraint fit for K + K − π + π − candidates in the data (points). The open histogram is the distribution for simulated signal events, normalized as described in the text. The cross-hatched (hatched) histogram represents the background from non-ISR events (plus that from ISR 4π events), estimated as described in the text.
two non-overlapping π 0 candidates plus the ISR photon direction, the two tracks and the beam particles. Both π 0 candidates are constrained to the π 0 mass, and we retain the combination with the lowest χ 2 KKπ 0 π 0 .

A. Final Selection and Backgrounds
The experimental χ 2 KKπ + π − distribution for the K + K − π + π − candidates is shown in Fig. 1 as points, and the open histogram is the distribution for the simulated K + K − π + π − events. The simulated distribution is normalized to the data in the region χ 2 KKπ + π − < 10 where the backgrounds and radiative corrections are insignificant. The experimental distribution has contributions from background processes, but the simulated distribution is also broader than the expected 3C χ 2 distribution. This is due to multiple soft-photon emission from the initial state and radiation from the final-state charged particles, which are not taken into account by the fit, but are present in both data and simulation. The shape of the χ 2 distribution at high values was studied in detail [11,12] using specific ISR processes for which a very clean sample can be obtained without any limit on the χ 2 value.
The cross-hatched histogram in Fig. 1 represents the background from e + e − → qq events, which is based on the JETSET simulation. It is dominated by events with a hard π 0 producing a fake ISR photon, and the similar kinematics cause it to peak at low values of χ 2 KKπ + π − .
We evaluate this background in a number of E c.m. ranges by combining the ISR photon candidate with another photon candidate in both data and simulated events, and comparing the π 0 signals in the resulting γγ invariant mass distributions. The simulation gives an E c.m.dependence consistent with the data, so we normalize it by an overall factor. The hatched histogram represents the sum of this background and that from ISR e + e − → π + π − π + π − events with one or two misidentified π ± , which also contributes at low χ 2 values. We estimate the contribution as a function of E c.m. from a simulation using the known cross section [11]. All remaining background sources are either negligible or give a χ 2 KKπ + π − distribution that is nearly uniform over the range shown in Fig. 1. We therefore define a signal region χ 2 KKπ + π − < 30, and estimate the sum of the remaining backgrounds from the difference between the number of data and simulated entries in a control region, 30 < χ 2 KKπ + π − < 60. This difference is normalized to the corresponding difference in the signal region, as described in detail in Refs. [11,12]. The signal region contains 34635 data and 14077 simulated events, and the control region contains 4634 data and 723 simulated events. The invariant mass distribution for K + K − π + π − candidates in the data (points): the cross-hatched, hatched and open histograms represent, cumulatively, the non-ISR background, the contribution from ISR π + π − π + π − events, and the ISR background from the control region of Fig. 1. Figure 2 shows the K + K − π + π − invariant mass distribution from threshold up to 5.0 GeV/c 2 for events in the signal region. Narrow peaks are apparent at the J/ψ and ψ(2S) masses. The cross-hatched histogram represents the qq background, which is negligible at low mass but becomes large at higher masses. The hatched region represents the ISR π + π − π + π − contribution, which we estimate to be 2.4% of the selected events on average. The open histogram represents the sum of all backgrounds, including those estimated from the control region. They total 6-8% at low mass but account for 20-25% of the observed data near 4 GeV/c 2 and become the largest contribution near 5 GeV/c 2 .
We subtract the sum of backgrounds in each mass bin to obtain a number of signal events. Considering uncertainties in the cross sections for the background processes, the normalization of events in the control region and the simulation statistics, we estimate a systematic uncertainty on the signal yield that is less than 3% in the 1.6-3 GeV/c 2 mass region, but increases to 3-5% in the region above 3 GeV/c 2 .

B. Selection Efficiency
The selection procedures applied to the data are also applied to the simulated signal samples. The resulting K + K − π + π − invariant-mass distributions in the signal and control regions are shown in Fig. 3(a) for the phase space simulation. The broad, smooth mass distribution is chosen to facilitate the estimation of the efficiency as a function of mass, and this model reproduces the observed distributions of kaon and pion momenta and polar angles. We divide the number of reconstructed simulated events in each mass interval by the number generated in that interval to obtain the efficiency shown as the points in Fig. 3(b). The 3 rd order polynomial fit to the points is used for further calculations. We simulate events with the ISR photon confined to the angular range 20-160 • with respect to the electron beam in the e + e − c.m. frame, which is about 30% wider than the EMC acceptance. This efficiency is for this fiducial region, but includes the acceptance for the final-state hadrons, the inefficiencies of the detector subsystems, and event loss due to additional soft-photon emission.
The simulations including the φ(1020)π + π − and/or K + K − f 0 (980) channels have very different mass and angular distributions in the K + K − π + π − rest frame. However, the angular acceptance is quite uniform for ISR events, and the efficiencies are consistent with those from the phase space simulation within 3%. To study possible mis-modeling of the acceptance, we repeat the analysis with the tighter requirements that all charged tracks be within the DIRC acceptance, 0.45 < θ ch < 2.4 radians, and the ISR photon be well away from the edges of the EMC, 0.35 < θ ISR < 2.4 radians. The fraction of selected data events satisfying the tighter requirements differs from the simulated ratio by 3.7%. We conservatively take the sum in quadrature of this variation and the 3% model variation (5% total) as a systematic uncertainty due to acceptance and model dependence.
We correct for mis-modeling of the shape of the χ 2 KKπ + π − distribution by (3.0±2.0)% and the track finding efficiency following the procedures described in detail in Ref. [11]. We use a comparison of data and simulated χ 2 4π distributions in the much larger samples of ISR π + π − π + π − events. We consider data and simulated  events that contain a high-energy photon plus exactly three charged tracks and satisfy a set of kinematical criteria, including a good χ 2 from a kinematic fit under the hypothesis that there is exactly one missing track in the event. We find that the simulated track-finding efficiency is overestimated by (0.8 ± 0.5)% per track, so we apply a correction of +(3 ± 2)% to the signal yield. We correct the simulated kaon identification efficiency using e + e − → φ(1020)γ → K + K − γ events. Events with a hard ISR photon and two charged tracks, one of which is identified as a kaon, with a K + K − invariant mass near the φ mass provide a very clean sample, and we compare the fractions of data and simulated events with the other track also identified as a kaon, as a function of momentum. The data-simulation efficiency ratio averages 0.990 ± 0.001 in the 1-5 GeV/c momentum range with variations at the 0.01 level. We conservatively apply a correction of +(1.0 ± 1.0)% per kaon, or +(2.0 ± 2.0)% to the signal yield.
C. Cross Section for e + e − → K + K − π + π − We calculate the e + e − → K + K − π + π − cross section as a function of the effective c.m. energy from where E c.m. ≡ m KKπ + π − c 2 , m KKπ + π − is the measured invariant mass of the K + K − π + π − system, dN KKπ + π − γ is the number of selected events after background subtraction in the interval dE c.m. , and ǫ KKπ + π − (E c.m. ) is the corrected detection efficiency. We calculate the differential luminosity, dL(E c.m. ), in each interval dE c.m. from ISR µ + µ − γ events with the photon in the same fiducial range used for the simulation; the procedure is described in Refs. [11,12]. From data-simulation comparison we conservatively estimate a systematic uncertainty on dL of 3%. This dL has been corrected for vacuum polar-TABLE II: Summary of corrections and systematic uncertainties on the e + e − → K + K − π + π − cross section. The total correction is the linear sum of the components and the total uncertainty is the sum in quadrature.
For the cross section measurement we use the tighter angular criteria on the charged tracks and the ISR photon, discussed in Sec. IV B, to exclude possible errors from incorrect simulation of the EMC and DCH edge effects. We show the cross section as a function of E c.m. in Fig. 4, with statistical errors only, and provide a list of our results in Table I. The result is consistent with the direct measurement by DM1 [13], and with our previous measurement of this channel [11] but has much better statistical precision. The systematic uncertainties, summarized in Table II, affect the normalization, but have little effect on the energy dependence.
The cross section rises from threshold to a peak value of about 4.7 nb near 1.85 GeV, then generally decreases with increasing energy. In addition to narrow peaks at the J/ψ and ψ(2S) masses, there are several possible wider structures in the 1.8-2.8 GeV region. Such structures might be due to thresholds for intermediate resonant states, such as φf 0 (980) near 2 GeV. Gaussian fits to the simulated line shapes give a resolution on the measured K + K − π + π − mass that varies between 4.2 MeV/c 2 in the 1.5-2.5 GeV/c 2 region and 5.5 MeV/c 2 in the 2.5-3.5 GeV/c 2 region. The resolution function is not purely Gaussian due to soft-photon radiation, but less than 10% of the signal is outside the 25 MeV/c 2 mass bin. Since the cross section has no sharp structure other than the J/ψ and ψ(2S) peaks discussed in Sec. VIII below, we apply no correction for resolution.
Our previous study [11] showed many intermediate resonances in the K + K − π + π − final state. With the larger data sample used here, they can be seen more clearly and, in some cases, studied in detail. Figure 5(a) shows a scatter plot of the invariant mass of the K − π + pair versus that of the K + π − pair, and Fig. 5(b) shows the sum of the two projections. Here we have suppressed the contributions from φπ + π − and K + K − ρ(770) by requir- where m(φ) and m(ρ) values are taken from the Particle Data Group (PDG) tables [5]. Bands and peaks corresponding to the K * 0 (892) and K * 0 2 (1430) are visible. In Fig. 5(c) we show the sum of projections of the K * 0 (892) bands, defined by lines in Fig. 5(a), with events in the overlap region plotted only once. No K * 0 (892) signal is seen, confirming that the e + e − → K * 0 (892)K * 0 (892) cross section is small. We observe associated K * 0 (892)K * 0 2 (1430) production, but it is mostly from J/ψ decays (see Sec. VIII).
We combine K * 0 /K * 0 candidates within the lines in Fig. 5(a) with the remaining pion and kaon to obtain the K * 0 π +− invariant mass distribution shown in Fig. 6(a), and the K * 0 π +− vs. K * 0 K −+ mass scatter plot in Fig. 6(b). The bulk of Fig. 6(b) shows a strong positive correlation, characteristic of K * 0 Kπ final states with no higher resonances. The horizontal band in Fig. 6(b) corresponds to the peak region in Fig. 6(a), and is consistent with contributions from the K 1 (1270) and K 1 (1400) resonances. There is also an indication of a vertical band in Fig. 6(b), perhaps corresponding to a K * 0 K resonance at ∼1.5 GeV/c 2 . We now suppress K * 0 Kπ by considering only events outside the lines in Fig. 5(a). In Fig. 7 we show the K ± π + π − invariant mass (two entries per event) vs. that of the π + π − pair, along with its two projections. There is a strong ρ(770) → π + π − signal, and the K ± π + π − mass projection shows further indications of the K 1 (1270) and K 1 (1400) resonances, both of which decay into Kρ(770). There are suggestions of additional structure in the π + π − mass distribution, including possible f 0 (980) shoulder and a possible enhancement near the f 2 (1270), however the current statistics do not allow us to make definitive statements.
ate states involving relatively wide resonances requires a partial wave analysis. This is beyond the scope of this paper. Here we present the cross section for the sum of all states including a K * 0 (892), and study intermediate states that include a narrow φ or f 0 resonance.
Signals for the K * 0 (892) and K * 0 2 (1430) are clearly visible in the K ± π ∓ mass distributions in Fig. 5(b) and, with a different bin size, in Fig. 8(a). We perform a fit to this distribution using P-wave Breit-Wigner (BW) functions for the K * 0 and K * 0 2 signals and a third-order polynomial function for the remainder of the distribution taking into account the Kπ threshold. The result is shown in Fig. 8(a). The fit yields a K * 0 signal of 19738 ± 266 events with m(Kπ) = 896.2 ± 0.3 MeV/c 2 and Γ(Kπ) = 50.6 ± 0.9 MeV, and a K * 0 2 signal of 1786±127 events with m(Kπ) = 1428.5±3.9 MeV/c 2 and Γ(Kπ) = 113.7 ± 9.2 MeV. These values are consistent with current world averages [5], and the fit describes the data well, indicating that contributions from any other resonances decaying into K ± π ∓ are small.
We perform a similar fit to the data in bins of the K + K − π + π − invariant mass, with the resonance masses and widths fixed to the values obtained by the overall fit. Since there is at most one K * 0 per event, we convert the resulting K * 0 yield in each bin into an "inclusive" e + e − → K * 0 Kπ cross section, following the procedure described in Sec. IV C. This cross section is shown in Fig. 8(b) and listed in Table III for the effective c.m. energies from threshold up to 3.5 GeV. At higher energies the signals are small and contain an unknown, but possibly large, contribution from e + e − → qq events. There is a rapid rise from threshold to a peak value of about 4 nb at 1.84 GeV, followed by a very rapid decrease with increasing energy. There are suggestions of narrow structure in the peak region, but the only statistically significant structure is the J/ψ peak, which is discussed below.
The e + e − → K * 0 Kπ contribution is a large fraction of the total K + K − π + π − cross section at all energies above its threshold, and dominates in the 1.8-2.0 GeV region. We are unable to extract a meaningful measurement of the K * 0 2 Kπ cross section from this data sample because it is more than ten times smaller. The K + K − ρ 0 (770) intermediate state makes up the majority of the remainder of the cross section and it can be estimated as a difference of the values in Table I and Table III for the K + K − π + π − and K * 0 Kπ final states.
Intermediate states containing relatively narrow resonances can be studied more easily. Figure 9(a) shows a scatter plot of the invariant mass of the π + π − pair versus that of the K + K − pair. Horizontal and vertical bands corresponding to the ρ 0 (770) and φ, respectively, are visible, and there is a concentration of entries on the φ band corresponding to the correlated production of φ and f 0 (980). The φ signal is also visible in the K + K − mass projection, Fig. 9(c). The large contribution from the ρ(770), coming from the K 1 decay, is nearly uniform in the K + K − mass, and the cross-hatched histogram shows the non-K + K − π + π − background estimated from the control region in χ 2 KKπ + π − . The cross-hatched his-  togram also shows a φ peak, but this is a small fraction of the events. Subtracting this background and fitting the remaining data gives 1706±56 events produced via the φπ + π − intermediate state.
To study the φπ + π − channel, we select candidate events with a K + K − invariant mass within 10 MeV/c 2 of the φ mass, indicated by the inner vertical lines in Figs. 9(a,c), and estimate the non-φ contribution from the mass sidebands between the inner and outer vertical lines. In Fig. 9(b) we show the π + π − invariant mass distributions for φ candidate events, sideband events and χ 2 control region events as the open, hatched and crosshatched histograms, respectively, and in Fig. 9(d) we show the numbers of entries from the candidate events minus those from the sideband and control region. There is a clear f 0 peak over a broad mass distribution, with no indication of associated ρ 0 production.
A coherent sum of two Breit-Wigner functions is sufficient to describe the invariant mass distribution of the π + π − pair recoiling against a φ in Fig. 9(d). We fit with the function: where m is the π + π − invariant mass, m i and Γ i are the parameters of the i th resonance, ψ is their relative phase and N i are normalization parameters, corresponding to the number of events under each BW. One BW corresponds to the f 0 (980), but a wide range of values of the other parameters can describe the data. Fixing the relative phase to ψ = π and the parameters of the first BW to m 1 = 0.6 GeV/c 2 and Γ 1 = 0.45 GeV (which could be interpreted as the f 0 (600) [5]), we obtain the fit shown in Fig. 9(d). It describes the data well and gives an f 0 (980) signal of 262±30 events, with m 2 = 0.973±0.003 GeV/c 2 and Γ 2 = 0.065 ± 0.013 GeV, consistent with the PDG values [5]. There is a suggestion of an f 2 (1270) peak in the data, but it is much smaller than the f 0 peak and we do not consider it further.
We obtain the number of e + e − → φπ + π − events in bins of φπ + π − invariant mass by fitting the K + K − in-variant mass projection in that bin after subtracting non-K + K − π + π − background. Each projection is a subset of Fig. 9(c), where the curve represent a fit to the full sample. In each mass bin, all parameters are fixed to the values obtained from the overall fit except the numbers of events in the φ peak and the non-φ component.
The efficiency may depend on the details of the production mechanism. Using the two-pion mass distribution in Fig. 9(d) as input, we simulate the π + π − system as an S-wave comprising two scalar resonances, with parameters set to the values given above. To describe the φπ + π − mass distribution we use a simple model with one resonance, the φ(1680), of mass 1.68 GeV/c 2 and width 0.2 GeV, decaying to φf 0 . The simulated reconstructed spectrum is shown in Fig. 10(a). There is a sharp increase at about 2 GeV/c 2 due to the φf 0 (980) threshold. All other structure is determined by phase space and a m −2 falloff with increasing mass. Dividing the number of reconstructed events in each bin by the number of generated ones, we obtain the efficiency as a function of φπ + π − mass shown in Fig. 10(b). The solid line represents a fit to a third order polynomial, and the dashed line the corresponding fit to the phase space model from Fig. 3. The model dependence is weak, giving confidence in the efficiency calculation. We calculate the e + e − → φπ + π − cross section as described in Sec. IV C but using the efficiency from the fit to Fig. 10(b) and dividing by the φ → K + K − branching fraction of 0.491 [5]. We show our results as a function of energy in Fig. 11 and list them in Table IV. The cross section has a peak value of about 0.6 nb at about 1.7 GeV, then decreases with increasing energy until φ(1020)f 0 (980) threshold, around 2.0 GeV. From this point it rises, falls sharply at about 2.2 GeV, and then decreases slowly. Except in the charmonium region, the results at energies above 2.9 GeV are not meaningful due to small signals and potentially large backgrounds, and are omitted from Table IV. Figure 11 displays the cross-section up to 4.5 GeV to show the signals from the J/ψ and ψ(2S) decays. They are discussed in Sec. VIII. There are no previous measurements of this cross section, and our results are consistent with the upper limits given in Ref. [13]. We perform a study of the angular distributions in the φ(1020)π + π − final state by considering all K + K − π + π − candidate events with mass below 3 GeV/c 2 , binning them in terms of the cosine of the angles defined below, and fitting the background-subtracted K + K − mass projections. The efficiency is nearly uniform in these angles, so we study the number of events in each bin. We define the φ production angle, Θ φ as the angle between the φ momentum and the e − beam direction in the rest frame of the φπ + π − system. The distribution of cos Θ φ , shown in Fig. 12(a), is consistent with the uniform distribution expected if S-wave two-body channels φX, X → π + π − dominate the φπ + π − system. We define the pion and kaon helicity angles, Θ π + and Θ K + as those between the π + and the π + π − -system momenta in the π + π − rest frame and between the K + and ISR photon momenta in the φ rest frame, respectively. The distributions of cos Θ π + and cos Θ K + , shown in Figs. 12(b) and 12(c), respectively, are consistent with those expected from scalar and vector meson decays. IV: Measurements of the e + e − → φ(1020)π + π − cross section (errors are statistical only).  The narrow f 0 (980) peak seen in Fig. 9(d) allows the selection of a fairly clean sample of φf 0 events. We repeat the analysis just described with the additional requirement that the π + π − invariant mass be in the range 0.85-1.10 GeV/c 2 . The fit to the full sample yields about 700 events; all of these contain a true φ, but about 10% are from e + e − → φπ + π − events where the pion pair is not produced through the f 0 (980).
We convert the numbers of fitted events in each mass bin into a measurement of the e + e − → φ(1020)f 0 (980) cross section as described above and dividing by the f 0 → π + π − branching fraction of two-thirds. The cross section is shown in Fig. 13 as a function of the effective c.m. energy and is listed in Table V. Its behavior near threshold does not appear to be smooth, but is more consistent with a steep rise to a value of about 0.3 nb at 1.95 GeV followed by a slow decrease that is interrupted by a structure around 2.175 GeV. Possible interpretations of this structure are discussed in Sec. VII. Again, the values are not meaningful for the effective c.m. above about 2.9 GeV/c 2 , except for the J/ψ and ψ(2S) signals, discussed in Sec. VIII.

A. Final Selection and Backgrounds
The K + K − π 0 π 0 sample contains background from the ISR processes e + e − → K + K − π 0 γ and K + K − ηγ, in which two soft photon candidates from machine-or detector-related background combine with the relatively energetic photons from the π 0 or η to form two fake π 0 candidates. We reduce this background using the helicity angle between each reconstructed π 0 direction and the direction of its higher-energy photon daughter calculated in its rest frame. If the cosines of both helicity angles are higher than 0.85, we remove the event. Figure 14 shows the distribution of χ 2 KKπ 0 π 0 for the remaining candidates together with the simulated FIG. 13: The e + e − → φ(1020)f0(980) cross section as a function of the effective e + e − c.m. energy obtained from the K + K − π + π − final state.
K + K − π 0 π 0 events. Again, the distributions are broader than those for a typical 6C χ 2 due to higher order ISR, and we normalize the histogram to the data in the region χ 2 KKπ 0 π 0 < 10. The cross-hatched histogram in Fig. 14 represents background from e + e − → qq events, evaluated in the same way as for the K + K − π + π − final state. The hatched histogram represents the sum of this background and that from ISR π + π − π 0 π 0 events with both charged pions misidentified as kaons, evaluated using the simulation.
The dominant background in this case is from residual ISR K + K − π 0 and K + K − η events, as well as ISR K + K − π 0 π 0 π 0 events. Their simulated contribution, shown as the dashed histogram in Fig. 14, is consistent with the data in the high χ 2 KKπ 0 π 0 region. All other backgrounds are either negligible or distributed uniformly in χ 2 KKπ 0 π 0 . We define a signal region, χ 2 KKπ 0 π 0 < 50, containing 4425 data and 6948 simulated events, and a control region, 50 < χ 2 KKπ 0 π 0 < 100, containing 1751 data FIG. 14: Distribution of χ 2 from the six-constraint fit for K + K − π 0 π 0 candidates in the data (points). The open histogram is the distribution for simulated signal events, normalized as described in the text. The cross-hatched, hatched and dashed histograms represent, cumulatively, the backgrounds from non-ISR events, ISR π + π − π 0 π 0 events, and ISR K + K − π 0 , K + K − η and K + K − π 0 π 0 π 0 events. and 848 simulated events. Figure 15 shows the K + K − π 0 π 0 invariant mass distribution from threshold up to 5 GeV/c 2 for events in the signal region. The qq background (cross-hatched histogram) is negligible at low masses but forms a large fraction of the selected events above about 4 GeV/c 2 . The ISR π + π − π 0 π 0 contribution (hatched region) is negligible except in the 1.5-2.5 GeV/c 2 region. The sum of all other backgrounds, estimated from the control region, is the dominant contribution below 1.6 GeV/c 2 and non negligible everywhere. The total background in the 1.6-2.5 GeV/c 2 region is 15-20% (open histogram).
We subtract the sum of backgrounds from the number of selected events in each mass bin to obtain a number of signal events. Considering uncertainties in the cross sections for the background processes, the normalization open histograms represent, cumulatively, the non-ISR background, the contribution from ISR π + π − π 0 π 0 events, and the ISR background from the control region of Fig. 14. of events in the control region and the simulation statistics, we estimate a systematic uncertainty on the signal yield after background subtraction as less than 5% in the 1.6-3.0 GeV/c 2 region, but increases to 10% in the region above 3 GeV/c 2 .

B. Selection Efficiency
The detection efficiency is determined in the same manner as in Sec. IV B. Figure 16(a) shows the simulated K + K − π 0 π 0 invariant mass distributions in the signal and control regions from the phase space model. We divide the number of reconstructed events in each 40 MeV/c 2 mass interval by the number generated ones in that interval to obtain the efficiency shown as the points in Fig. 16(b); a third order polynomial fit to the efficiency is used to calculate the cross section. Again, the simulation of the ISR photon covers a limited angular range, about 30% wider than EMC acceptance, and shown efficiency is factor 0.7 lower than for the hadronic system alone. Simulations assuming dominance of the φ → K + K − and/or f 0 → π 0 π 0 channels give consistent results, and we apply the same 5% systematic uncertainty for possible model dependence as in Sec. IV B.
We correct for mis-modeling of the track finding and kaon identification efficiencies as in Sec. IV B, and for the shape of the χ 2 KKπ 0 π 0 distribution analogously, using the result in Ref. [12], (0 ± 6)%. We correct the π 0finding efficiency using the procedure described in detail in Ref. [12]. From ISR e + e − → ωπ 0 γ→ π + π − π 0 π 0 γ events selected with and without the π 0 from the ω decay, we find that the simulated efficiency for one π 0 is too high by (2.8±1.4)%. Conservatively we apply a correction of +(5.6 ± 2.8)% for two π 0 in the event.
C. Cross Section for e + e − → K + K − π 0 π 0 We calculate the cross section for e + e − → K + K − π 0 π 0 in 40 MeV E c.m. intervals from the analog of Eq. 2, using the invariant mass of the K + K − π 0 π 0 system to determine the effective c.m. energy. We show the first measurement of this cross section in Fig. 17 and list the results obtained in Table VI. The cross section rises to a peak value near 1 nb at 2 GeV, falls sharply at 2.2 GeV, then decreases slowly. The only statistically significant structure is the J/ψ peak. The drop at 2.2 GeV is similar to that seen in the K + K − π + π − mode. Again, dL includes corrections for vacuum polarization that should be omitted from calculations of g µ −2.
The simulated K + K − π 0 π 0 invariant mass resolution is 8.8 MeV/c 2 in the 1.5-2.5 GeV/c 2 mass range, and increases with mass to 11.2 MeV/c 2 in the 2.5-3.5 GeV/c 2 range. Since less than 20% of the events in a 40 MeV/c 2 bin are reconstructed outside that bin and the cross section has no sharp structure other than the J/ψ peak, we again make no correction for resolution. The point-topoint systematic errors are much smaller than statistical ones, and the errors on the normalization are summarized in Table VII, along with the corrections that were applied to the measurements. The total correction is +9.2%, and the total systematic uncertainty is 10% at low mass, increasing to 14% above 3 GeV/c 2 .  VII: Summary of corrections and systematic uncertainties on the e + e − → K + K − π 0 π 0 cross section. The total correction is the linear sum of the components and the total uncertainty is the sum in quadrature.

Source Correction Uncertainty
Rad. Corrections -1% Backgrounds -5%, m KKπ 0 π 0 < 3 GeV/c 2 10%, m KKπ 0 π 0 > 3 GeV/c 2 Model Dependence -5% χ 2 KKπ 0 π 0 Distn. 0% 6% Tracking Efficiency +1.6% 0.8% Kaon ID Efficiency +2% 2% π 0 Efficiency +5.6% 2.8% ISR Luminosity -3% Total +9.2% 10%, m KKπ 0 π 0 < 3 GeV/c 2 14%, m KKπ 0 π 0 > 3 GeV/c 2 D. Substructure in the K + K − π 0 π 0 Final State A scatter plot of the invariant mass of the K − π 0 versus that of the K + π 0 pair is shown in Fig. 18(a) with two entries per event selected in the χ 2 signal region. Horizontal and vertical bands corresponding to the K * + (892) and K * − (892), respectively, are visible. Figure 18(b) shows as points the sum of the two projections of Fig. 18(a); a large K * ± (892) signal is evident. Fitting this distribution with the function discussed in Sec. IV E gives a good χ 2 and the curve shown on Fig. 18(b). The K * ± (1430):K * ± (892) ratio is consistent with that for neutral K * seen in the K + K − π + π − channel, and the number of K * ± (892) in the peak is consistent with one per selected event. The hatched histogram in Fig. 18(b) represents the K ± π 0 mass in events with the other K ∓ π 0 mass within the lines in Fig. 18(a), but with events in the overlap region used only once, and shows no K * ± (892) signal. These results indicate that the e + e − → K * ± K * ∓ cross section is small and that the K * ± (892)K ∓ π 0 chan-  18: (a) Invariant mass of the K − π 0 pair versus that of the K + π 0 pair in selected K + K − π 0 π 0 events (two entries per event); (b) sum of projections of (a) (dots, four entries per event). The curve represents the result of the fit described in the text. The hatched histogram is the K ± π 0 distribution for events in which the other K ∓ π 0 combination is within the K * ± (892) bands indicated in (a), with events in the overlap region taken only once. nels dominate the overall cross section.
We find no signals for resonances in the K + K − π 0 or K ± π 0 π 0 decay modes. Since the K * ± (892)K ∓ π 0 channels dominate and the statistics are low in any mass bin, we do not attempt to extract a separate K * ± (892)K ∓ π 0 cross section. The total K + K − π 0 π 0 cross section is roughly a factor of four lower than the K * 0 (892)K ± π ∓ cross section observed in the K + K − π + π − final state. This is consistent with what one might expect from isospin and the charged vs. neutral K * branching fractions into charged kaons.
E. The φ(1020)π 0 π 0 Intermediate State The selection of events containing a φ(1020) → K + K − decay follows that in Section IV F. Figure 19(a) shows a scatter plot of the invariant mass of the π 0 π 0 pair versus that of the K + K − pair. A vertical band corresponding to the φ is visible, whose intensity decreases with increasing π 0 π 0 mass except for an enhancement in the f 0 (980) region. The φ signal is also visible in the K + K − invariant mass projection shown in Fig. 19(c). The relative non-φ background is smaller than in the K + K − π + π − mode, but there is a large background from ISR φπ 0 , φη and/or φπ 0 π 0 π 0 events, as indicated by the control region histogram (hatched) in Fig. 19(c). The contributions from non-ISR and ISR π + π − π 0 π 0 events are negligible. Selecting φ candidate and side band events as for the K + K − π + π − mode (vertical lines in Figs. 19(a,c)), we obtain the π 0 π 0 mass projections shown as the open and cross-hatched histograms, respectively, in Fig. 19(b). Control region events (hatched histogram) are concentrated at low masses. A peak corresponding to the f 0 (980) is visible over a relatively low background.
In Fig. 19(d) we show the numbers of entries from the candidate events minus those from the sideband and control regions. A sum of two Breit-Wigner functions is again sufficient to describe the data. Fitting Eq. 3 with the parameters of one BW fixed to the values given in Sec. IV F, corresponding to the f 0 (600), we obtain a good fit, shown as the curve in Fig. 19(d). This fit yields a f 0 (980) signal of 54 ± 9 events with a mass m = 0.970 ± 0.007 GeV/c 2 and width Γ = 0.081 ± 0.021 GeV consistent with PDG values [5]. Due to low statistics and high backgrounds, we do not extract an e + e − → φ(1020)π 0 π 0 cross section. Since the background under the f 0 (980) peak in Figs. 19(b,d) is relatively low we are able to extract the φ(1020)f 0 (980) contribution. As in Sec. IV G, we require the dipion mass to be in the range 0.85-1.10 GeV/c 2 and fit the background-subtracted K + K − mass projection in each bin of K + K − π 0 π 0 mass to obtain a number of φf 0 events. Again, about 10% of these are φπ 0 π 0 events in which the π 0 π 0 pair is not produced through the f 0 , but this does not affect the conclusions.
We convert the number of fitted events in each mass bin into a measurement of the e + e − → φ(1020)f 0 (980) cross section as described above and dividing by the f 0 (980) → π 0 π 0 branching fraction of one-third. The cross section is shown in Fig. 20 as a function of E c.m. and is listed in Table VIII. Due to smaller number of events, we have used larger bins at higher energies. The overall shape is consistent with that obtained in the K + K − π + π − mode (see Fig. 13), and there is a sharp   Figure 21 shows the distribution of χ 2 4K for the K + K − K + K − candidates as points, and the open histogram is the distribution for simulated K + K − K + K − events, normalized to the data in the region χ 2 4K < 5 where the backgrounds and radiative corrections are small. The hatched histogram represents the background from e + e − → qq events, evaluated as for the other modes. The cross-hatched histogram represents the background from simulated ISR K + K − π + π − events with both charged pions misidentified as kaons.
We define signal and control regions of χ 2 4K < 20 and 20 < χ 2 4K < 40, respectively. The signal region contains 2,305 data and 20,616 simulated events, and the control region contains 463 data and 1,601 simulated events. . The open histogram is the distribution for simulated signal events, normalized as described in the text. The hatched histogram represents the background from non-ISR events, estimated as described in the text. The cross-hatched histograms is for simulated ISR K + K − π + π − events. Figure 22 shows the K + K − K + K − invariant mass distribution from threshold up to 5 GeV/c 2 for events in the signal region as points with errors. The qq background (hatched histogram) is small at low masses, but dominant above about 4.5 GeV/c 2 . Since the ISR K + K − π + π − background does not peak at low χ 2 4K values, we include it in the background evaluated from the control region, according to the method explained in Sec. IV A. It dominates this background, which is 10% or lower at all masses. The total background is shown as the open histogram in Fig. 22.
We subtract the sum of backgrounds from the number of selected events in each mass bin to obtain a number of signal events. Considering uncertainties in the cross sections for the background processes, the normalization of events in the control region, and the simulation statistics, we estimate a systematic uncertainty on the signal yield of less than 5% in the 2-3 GeV/c 2 region, increasing to 10% in the region above 3 GeV/c 2 .

B. Selection Efficiency
The detection efficiency is determined as for the other two final states. Figure 23(a) shows the simulated K + K − K + K − invariant-mass distributions in the signal and control regions from the phase space model. We divide the number of reconstructed events in each mass interval by the number of generated ones in that interval to obtain the efficiency shown as the points in Fig. 23(b). It is quite uniform, and we fit a third order polynomial, which we use to extract the cross section. A factor of 0.7 is again applicable for only hadronic system efficiency due to the limited angular coverage of the ISR photon simulation. A simulation assuming dominance of the φK + K − channel, with the K + K − pair in an S-wave, gives consistent results, and we apply the same 5% systematic uncertainty as for the other final states. We correct for mis-modeling of the track finding and kaon identification efficiencies as in Sec. IV B, and for the shape of the χ 2 4K distribution analogously, using the result in Ref. [11], (3.0 ± 2.0)%.
C. Cross Section for e + e − → K + K − K + K − We calculate the e + e − → K + K − K + K − cross section in 40 MeV E c.m. intervals from the analog of Eq. 2, using the invariant mass of the K + K − K + K − system to determine the effective c.m. energy. We show this cross section in Fig. 24 and list it in Table IX. It rises to a peak value near 0.1 nb in the 2.3-2.7 GeV region, then decreases slowly with increasing energy. The only statistically significant narrow structure is the J/ψ peak. Again, dL includes corrections for vacuum polarization that should be omitted from calculations of g µ −2. This supersedes our previous result [11]. The simulated K + K − K + K − invariant mass resolution is 3.0 MeV/c 2 in the 2.0-2.5 GeV/c 2 range, increasing with mass to 4.7 MeV/c 2 in the 2.5-3.5 GeV/c 2 range. Since the cross section has no sharp structure except for the J/ψ peak, we again make no correction for resolution. The errors shown in Fig. 24 and Table IX are statistical only. The point-to-point systematic errors are much smaller than this, and the errors on the normalization are summarized in Table X, along with the corrections applied to the measurement. The total correction is +10%, and the total systematic uncertainty is 9% at low mass, increasing to 13% above 3 GeV/c 2 .  Figure 25 shows the invariant mass distribution for all K + K − pairs in the selected K + K − K + K − events (4 entries per event) as the open histogram. A prominent φ peak is visible, along with possible peaks at 1.5, 1.7 and 2.0 GeV/c 2 . The hatched histogram is for the pair in each event with mass closest to the nominal φ mass, and indicates that the φK + K − channel dominates the K + K − K + K − final state. Our previous finding of very little φ signal [11] was incorrect due to an error in the analysis algorithm. If the pair mass closest to the φ mass is within 10 MeV/c 2 of the φ mass, then we include the invariant mass of the other K + K − combination in Fig. 26. The contribution from events in the J/ψ peak is shown as the hatched histogram which is in agreement with the BES experiment [24] which studied the structures around 1.5, 1.7 and 2.0 GeV/c 2 in detail. There is no evidence for the φf 0 channel, but there is an enhancement at threshold that can be interpreted as the tail of the f 0 (980). This is expected in light of the φf 0 cross sections measured above in the K + K − π + π − and K + K − π 0 π 0 final states. However the statistics and uncertainties in the f 0 (980) → K + K − lineshape do not allow a meaningful extraction of the cross section in this final state. We observe no significant structure in the K + K − K ± mass distribution. We use these events to study the possibility that part of our φπ + π − signal is due to φK + K − events with the two kaons not from the φ taken as pions. No structure is present in the reconstructed K + K − π + π − invariant mass distribution from these events.
VII. e + e − → φf0 NEAR THRESHOLD The behavior of the e + e − → φf 0 cross section near threshold shows a structure near 2150 MeV/c 2 , and we have published this result in Ref. [14]. Here we provide a more detailed study of this cross section in the 1.8-3 GeV region. In Fig. 27 we superimpose the cross sections measured in the K + K − π + π − and K + K − π 0 π 0 final states (shown in Figs. 13 and 20); they are consistent with each other. The K + K − K + K − cross section (Sec. VI D) is also consistent with the presence of a structure near 2150 MeV/c 2 and shows a contribution from the φf 0 channel, but since we cannot extract a meaningful φf 0 cross section, we do not discuss this final state further.
First, we attempt to reproduce this spectrum with a smooth threshold function. In the absence of resonances, the only theoretical constraint on the cross section well above threshold is that it should decrease smoothly with increasing E c.m. . However the form of the cutoff at threshold is determined by the properties of the interme- The e + e − → φ(1020)f0(980) cross section measured in the K + K − π + π − (circles) and K + K − π 0 π 0 (squares) final states. The hatched histogram shows the simulated cross section, assuming no resonant structure. The solid (dashed) line represents the result of the one-resonance (no-resonance) fit described in the text.
diate resonances and the final state particle spins, phase space and detector resolution. The model discussed in Sec. IV F takes the φ and f 0 (980) lineshapes, the spins of all particles and their phase space into account, and postulates a simple E −4 c.m. dependence of the cross section. For the e + e − → φf 0 reaction, it predicts the cross section shown as the hatched histogram in Fig. 27, normalized to the same total area as the data. It shows a sharp rise from the threshold with a peak near 2070 MeV and is inconsistent with the data.
To account for uncertainties in the f 0 width and the shape of the cross section well above threshold, we seek a functional form that describes the simulation and whose parameters can be varied to cover a reasonable range of possibilities. This can be achieved by the product of a phase space term, an exponential rise and a second order polynomial: where the a i are free parameters, σ 0 is a normalization factor, and P (µ) is a good approximation of the two-body phase space for particles with similar masses. Both the φ(1020) and  Table XI, and the latter is shown as the dashed curve on Fig. 27; all fits are inconsistent with the data. We now add a resonance and fit the data with the function where m 1 and Γ 1 are the mass and width of the resonance, σ 1 is its peak cross section, and ψ 1 is its phase relative to the non-resonant component. We obtain good fits both assuming no interference between the two components, ψ 1 = π, and with ψ 1 floating. The result of the latter fit is shown as the solid curve on Fig. 27. The data are somewhat above this curve near 2.4 GeV/c 2 and a fit with two resonances can also describe the data. Due to the sharp drop near 2.2 GeV/c 2 , the single-resonance fit with interference gives a resonance mass about 30 MeV/c 2 higher than the other two fits. All these fits, with or without resonances, give a peak non-resonant cross section in the range 0.3-0.4 nb, which is of independent theoretical interest, because it can be related to the φ → f 0 (980)γ decay studied at the φ-factory [25].
Under the hypothesis of one resonance interfering with the non-resonant component, the fit gives the resonance parameters σ 1 = 0.13 ± 0.04 nb, m 1 = 2.175 ± 0.010 GeV/c 2 , Γ 1 = 0.058 ± 0.016 GeV, and ψ 1 = −0.57 ± 0.30 radians, along with χ 2 /n.d.f.= 37.6/(56 − 9) (C.L. 0.84). We can estimate the product of its electronic width and branching fraction to φf 0 as where we fit the product Γ 1 σ 1 to reduce correlations, and the conversion constant C = 0.389 mb( GeV/c 2 ) 2 . The second error is systematic and corresponds to the normalization errors on the cross section. The significance of the structure calculated from the change in χ 2 between the best fit and the null hypothesis is 6.2 standard deviations. Since this calculation can be unreliable in the case of low statistics and functions that vary rapidly on the scale of the bin size, we perform a set of simulations in which we generate a number of events according to a Poisson distribution about the number observed in the data and with a mass distribution given by either the simulation or fitted function in Fig. 27 without resonant structure. On each sample, we perform fits to Eqs. 4 and 5 and calculate the difference in χ 2 . The fraction of trials giving a χ 2 difference larger than that seen in the data corresponds to a significance of approximately 5 standard deviations.
We search for this structure in other submodes with different and/or fewer intermediate resonances. The total cross sections are dominated by K * Kπ channels, and the K * 0 K + π − cross section is shown in Fig. 8. There is no significant structure in the 2.1-2.5 GeV region, but the point-to-point statistical uncertainties are large. If we remove events within the bands in Figs. 5 and 18, then most of the events containing a K * are eliminated and we obtain the raw mass distributions shown as the points with errors in Figs. 28 and 29, respectively. Both distribution show evidence of a structure around 2.15 GeV/c 2 and the K + K − π + π − distribution also shows a structure near 2.4 GeV/c 2 . We cannot exclude the presence of these structures in events with a K * , but we can conclude that they do not dominate those events, whereas they comprise a substantial fraction of the remaining events in that mass region.
Applying the further requirement that the dipion mass be in the range 0.85-1.10 GeV/c 2 , we remove most of the events without an f 0 , and obtain the mass distributions shown as the hatched histograms in Figs. 28 and 29. Peaks are visible at both 2.15 GeV/c 2 and 2.4 GeV/c 2 in both distributions, and they contain enough events to account for the corresponding structures in the distributions for all non-K * events. These peaks contain at least as many events as are present in the φf 0 samples, but the non-resonant components are higher and there is a substantial kinematic overlap between K + K − f 0 events and K * Kπ events in this mass range.
Since this f 0 (980) band appears to contain a large fraction of the events within the structure, we now consider all selected events with a dipion mass inside or outside this range. Figure 31 shows the mass distribution for all selected K + K − π 0 π 0 events as the open histogram, and the subsets of events with π 0 π 0 mass inside and outside the range 0.85-1.10 GeV/c 2 as the hatched and crosshatched histograms, respectively. It is evident that the K + K − f 0 channel contains the majority of the structure in the 2.0-2.6 GeV/c 2 range.
We show the corresponding distributions for the distribution for all selected K + K − π 0 π 0 events lying outside the K * 0 (892) bands of Fig. 18 (points), and the subset of these events with 0.85 < m(π 0 π 0 ) < 1.10 GeV/c 2 (hatched). Fig. 30. Due to the presence of the ρ 0 , the relative f 0 contribution is much smaller in this final state, but the events in the f 0 band show clear indications of structure in the 2.0-2.4 GeV/c 2 region. The remaining events may also have structure in this region, but the statistical significance is marginal and it could be due to other sources, such as the φf 2 (1270) threshold at 2.3 GeV/c 2 . Figures 32 and 33 show enlarged views of the mass distributions within the f 0 bands from Figs. 30 and 31, respectively. The two-peak structure is more evident here than in the φf 0 events. The 0.85 < m(ππ) < 1.10 GeV/c 2 requirement gives enough phase space for K + K − invariant mass to cover the region from threshold to ∼1.3 GeV/c 2 for m(K + K − ππ) ≈ 2.15 GeV/c 2 . From the measured kaon form factor we expect to find only about two-thirds of K + K − P-wave in our fitted φ peak. Since the non-ISR and ISR ππππ backgrounds have not been subtracted and the samples contain an unknown mixture of intermediate states, we fit them with a modi- XII: Summary of parameters obtained from the fits described in the text to the K + K − π + π − and K + K − π 0 π 0 events with dipion mass in the f0(980) band. An asterisk denotes a value that was fixed in that fit.

No Resonance
One Resonance Two Resonances Fit Here, the normalization is in terms of events rather than cross section (σ i → N i ) and a fraction a 4 of the nonresonant component does not interfere with the resonances. We first fit the distribution with no resonances (and a 4 = 1). The results are shown as the dashed lines in Figs. 32 and 33 and listed in Table XII; both are inconsistent with the data.
We next include one resonance in the fit. The parameter a 4 is not well constrained by the data and its value has a small influence on all other fit parameters except for the number of events assigned to the resonance, so we present results with a 4 fixed to the reasonable values of 0.75 and 0.50 for the K + K − π + π − and K + K − π 0 π 0 data, respectively. The results are shown as the solid lines in Figs. 32 and 33 and listed in Table XII. The fit quality is good in both cases, the fitted resonance parameters are consistent with those from the φf 0 study, and the calculated significance of the structure for the K + K − π + π − data is similar, 5.2 standard deviations. The K + K − π 0 π 0 data show much more pronounced structure than in the φf 0 study, allowing a full fit to this sample with a significance of 5.0 standard deviations.
We then add a second resonance to the fit, keeping a 4 fixed and floating all other parameters. The results are shown as the dotted lines in Figs. 32 and 33, and listed in Table XII. These fits are also of good quality, but do not change the χ 2 CL or the parameters of the first resonance significantly. We also perform fits with no interference between the non-resonant component and any resonance (a 4 = 1), obtaining good quality fits for both one resonance and two resonances with relative phase π/2. The fitted resonance parameters are consistent in all cases, except that the mass of the first resonance is lower by about 50 MeV/c 2 , similar to the 30 MeV/c 2 shift seen in the φf 0 study.
From these studies we conclude that we have observed a new vector structure at a mass of about 2150 MeV/c 2 with a significance of over six standard deviations. It decays into K + K − f 0 (980), with the K + K − pair produced predominantly via the φ(1020). There is an additional structure at about 2400 MeV/c 2 , and the two structures can be described by either two resonances or a single resonance that interferes with the non-resonant K + K − f 0 (980) process. More data and searches in other final states are needed to understand the nature of these structures.
If the main structure is due to a resonance, then it is relatively narrow and might be interpreted as the strange analog of the recently observed charmed Y(4260) state [6], which decays to J/ψπ + π − . The value of B φf0 · Γ ee = (2.5 ± 0.8 ± 0.4) eV measured here is similar to the value of B Y →J/ψπ + π − · Γ Y ee = (5.5 ± 1.0 ± 0.8) eV reported in Ref. [6].

VIII. THE CHARMONIUM REGION
The data at masses above 3 GeV/c 2 can be used to measure or set limits for the branching fractions of narrow resonances, such as charmonia, and the narrow J/ψ and ψ(2S) peaks allow measurements of our mass scale and resolution. Figures 34, 35 and 36 show the invariant mass distributions for the selected K + K − π + π − , K + K − π 0 π 0 and K + K − K + K − events, respectively, in this region, with finer binning than in the corresponding Figs. 2, 15 and 22. We do not subtract any background FIG. 32: The K + K − π + π − invariant mass distribution in the K + K − f0(980) threshold region for events with a π + π − mass inside the f0 band. The lines represent the results of the fits including no (dashed), one (solid) and two (dotted) resonances described in the text.
FIG. 33: The K + K − π 0 π 0 invariant mass distribution in the K + K − f0(980) threshold region for events with a π 0 π 0 mass inside the f0 band. The lines represent the results of the fits including no (dashed), one (solid) and two (dotted) resonances described in the text. from the K + K − π + π − or K + K − K + K − data, since it is small and nearly uniformly distributed, but we use the χ 2 KKπ 0 π 0 control region to subtract part of the ISR background from the K + K − π 0 π 0 data. Signals from the J/ψ are visible in all three distributions, and the ψ(2S) is visible in the K + K − π + π − mode.
We fit each of these distributions using a sum of two Gaussian functions to describe the J/ψ and ψ(2S) signals plus a polynomial to describe the remainder of the dis- tribution. We take the signal function parameters from the simulation, but let the overall mean and width float in the fit, along with the amplitude and the coefficients of the polynomial. The fits are of good quality and are shown as the curves on Figs. 34, 35 and 36. In all cases, the fitted mean value is within 1 MeV/c 2 of the PDG [5] J/ψ or ψ(2S) mass, and the width is consistent within 10% with the simulated resolution discussed in Sec. IV C, V C or VI C. The fits yield 1586 ± 58 events in the J/ψ peak for the K + K − π + π − final state, 203 ± 16 events for K + K − π 0 π 0 and 156 ± 15 events for K + K − K + K − . From these numbers of observed events in each final state f , N J/ψ→f , we calculate the product of the J/ψ branching fraction to f and the J/ψ electronic width: where dL/dE = 89.8 nb −1 / MeV and ǫ f (m J/ψ ) are the ISR luminosity and corrected selection efficiency, respectively, at the J/ψ mass and C is the conversion constant. We estimate ǫ K + K − π + π − = 0.202, ǫ K + K − π 0 π 0 = 0.069 and ǫ K + K − K + K − = 0.176. Using Γ J/ψ ee = 5.40 ± 0.18 keV [5], we obtain the branching fractions listed in Table XIII, along with the measured products and the current PDG values. The systematic errors include a 3% uncertainty on Γ J/ψ ee . The branching fractions to K + K − π + π − and K + K − K + K − are more precise than the current PDG values, which were dominated by our previous results of (6.25±0.80)×10 −3 and (7.4±1.8)×10 −4 , respectively [11]. This is the first measurement of the K + K − π 0 π 0 branching fraction.  These fits also yield 91±15 K + K − π + π − events in the ψ(2S) peak, but no other significant signals. We expect 6.3 events from ψ(2S) → J/ψπ + π − → K + K − π + π − from the relevant branching fractions [5], which is less than the statistical error. Subtracting this contribution and using a calculation analogous to Eq. 7, with dL/dE = 115.3 nb −1 / MeV, we obtain the product of the ψ(2S)→ K + K − π + π − branching fraction and its electronic width. Dividing by the world average value of Γ ψ(2S) ee [5], we obtain the branching fraction listed in Table XIII; it is consistent with the current PDG value [5].
As noted in Sec. IV D and shown in Fig. 5, the K + K − π + π − final state is dominated by the K * 0 (892)Kπ channels, with a small fraction seen in the K * 0 (892)K * 0 2 (1430) + c.c. channels. Figure 37 shows a scatter plot of the invariant mass of a K ± π ∓ pair versus that of the K + K − π + π − system in events with the other K ∓ π ± pair near the K * 0 (892) mass, i.e. within the bands in Fig. 5(a) with overlapped region taken only once. There is a large concentration of entries in the J/ψ band with K ± π ∓ masses near 1430 MeV/c 2 , but no solid evidence for a horizontal band corresponding to K * 0 2 (1430) production other than in J/ψ decays. We show the K ± π ∓ mass projection for the subset of events with a K + K − π + π − mass within 50 MeV/c 2 of the J/ψ mass in Fig. 38 as the open histogram. The hatched histogram is the projection for events with a K + K − π + π − mass between 50 and 100 MeV/c 2 below the J/ψ mass. FIG. 37: The K ± π ∓ invariant mass versus K + K − π + π − invariant mass for events with the other K ∓ π ± combination in the K * 0 (892) bands of Fig. 5(a). The overlapped region is taken only once.
The J/ψ component appears to be dominated by the K * 0 2 (1430). Also seen is a small signal from K * 0 (892) indicating the K * 0 (892)K * 0 (892) decay of J/ψ: this is also seen as an enhancement in the vertical J/ψ band in Fig. 37. The enhancement at 1.8 GeV/c 2 of Fig. 38 can be explained by the J/ψ decay into K * 0 (892)K 2 (1770)+c.c. (or K * 0 (892)K 2 (1820) + c.c.), a mode which has not previously been reported. Subtracting the number of sideband events from the number in the J/ψ mass window, we obtain 317±23 events with a K ± π ∓ mass in the range 1200-1700 MeV/c 2 , which we take as a measure of J/ψ decays into K * 0 (892)K * 0 2 (1430), 25 ± 8 events in the 0.8-1.0 GeV/c 2 window for the K * 0 (892)K * 0 (892) decay and 110 ± 14 events for the K * 0 (892)K 2 (1770) or K * 0 (892)K 2 (1830) final state in the 1.7-2.0 GeV/c 2 region. We convert these to branching fractions using Eq. 7 and dividing by the known branching fractions of excited kaons [5]. The results are listed in Table XIII: they are considerably more precise than the PDG values. We cannot calculate B J/ψ→K * 0 K2(1770) because no branching fractions of K 2 (1770) or K 2 (1830) to Kπ are reported.
We study decays into φπ + π − and φπ 0 π 0 using the mass distributions shown in Figs. 39 and 40, respectively. The open histograms are for the events with a K + K − mass within the φ bands of Figs. 9(c) and 19(c). The cross-hatched histogram in Fig. 39 is from the φ sidebands of Fig. 9(c) and represents the dominant background in the φπ + π − mode. The hatched histogram in Fig. 40 is from the χ 2 KKπ 0 π 0 control region and represents the dominant background in the φπ 0 π 0 mode. Subtracting these backgrounds, we find 103±12 J/ψ → φπ + π − events, 23±6 J/ψ → φπ 0 π 0 events, and 10±4 ψ(2S) → φπ + π − events. We convert these to branching fractions and list them in Table XIII. This is the first measurement of the J/ψ → φπ 0 π 0 branching fraction, and the other two are consistent with current PDG values.
If the Y (4260) has a substantial branching fraction into φπ + π − , then we would expect to see a signal in Fig. 39. In the mass range |m(φπ + π − ) − m(Y )| < 0.1 GeV/c 2 , we find 10 events, and assuming a uniform distribution we estimate 9.2 background events from the 3.8-5.0 GeV/c 2 region. This corresponds to a signal of 0.8 ± 3.3 events or a limit of < 5 events at the 90% C.L. Using dL/dE = 147.7 nb −1 / MeV at the Y (4260) mass, we calculate B Y →φπ + π − · Γ Y ee < 0.4 eV which is well below the value of B Y →J/ψπ + π − · Γ Y ee = (5.5 ± 1.0 ± 0.8) eV [6]. No Y (4260) signal is seen in any other mode studied here. both cases, but φf 0 is not the dominant mode of the J/ψ → φπ + π − decay. Figure 41(b) shows the π + π − invariant mass distribution for events in the J/ψ peak of Fig. 39, 3.05 < m(K + K − π + π − ) < 3.15 GeV/c 2 . A twopeak structure is visible that can interpreted as due to the f 0 (980) and f 2 (1270) resonances. Fitting the distribution in Fig. 41(b) with a sum of two Breit-Wigner functions with parameters fixed to PDG values [5], we find 19.5 ± 4.5 J/ψ → φf 0 events and 44 ± 7 J/ψ → φf 2 events. From Fig. 42 we estimate 7.0 ± 2.8 φf 0 events in the π 0 π 0 mode.  events (open histogram) and events in the φ side bands (cross-hatched) in the charmonium region; (b) the π + π − invariant mass distribution for φπ + π − events from the J/ψ peak of Fig. 39. The line represents the result of the fit described in the text.
Using Eq. 7 and dividing by the appropriate branching fractions, we obtain the J/ψ branching fractions listed in Table XIII. The measurements of B J/ψ→φf0 in the π + π − and π 0 π 0 decay modes of the f 0 are consistent with each other and with the PDG value, and combined they have roughly the same precision as listed in the PDG [5]. This is the first measurement of B J/ψ→φf2 , and the value is consistent with the previous upper limit [5]. We also KKπ 0 π 0 control region (hatched) in the charmonium region. observe 6 ± 3 ψ(2S) → φf 0 , f 0 → π + π − events, which we convert to the branching fraction listed in Table XIII; it is consistent with the PDG value [5], assuming B f0→π + π − = 2/3.