Study of the reactions $e^+e^-\to2(\pi^+\pi^-)\pi^0\pi^0\pi^0$ and $e^+e^-\to2(\pi^+\pi^-)\pi^0\pi^0\eta$ at center-of-mass energies from threshold to 4.5 GeV using initial-state radiation

We study the processes $e^+e^-\to 2(\pi^+\pi^-)\pi^0\pi^0\pi^0\gamma$ and $2(\pi^+\pi^-)\pi^0\pi^0\eta\gamma$ in which an energetic photon is radiated from the initial state. The data were collected with the \babar~ detector at SLAC. About 14\,000 and 4700 events, respectively, are selected from a data sample corresponding to an integrated luminosity of 469 fb^{-1}. The invariant mass of the hadronic final state defines the effective \epem center-of-mass energy. The center-of-mass energies range from threshold to 4.5 GeV. From the mass spectra, the first ever measurements of the $e^+e^-\to 2(\pi^+\pi^-)\pi^0\pi^0\pi^0$ and the $e^+e^-\to2(\pi^+\pi^-)\pi^0\pi^0\eta$ cross sections are performed. The contributions from $\omega\pi^+\pi^-\pi^0\pi^0$, $\eta2(\pi^+\pi^-)$, and other intermediate states are presented. We observe the $J/\psi$ and $\psi(2S)$ in most of these final states and measure the corresponding branching fractions, many of them for the first time.


I. INTRODUCTION
The Standard Model (SM) calculation of the muon anomalous magnetic moment (g µ −2) requires input from experimental e + e − hadronic cross section data in order to account for hadronic vacuum polarization (HVP) terms.In particular, the calculation is most sensitive to the lowenergy region, from the hadronic threshold to about 2 GeV, where the inclusive hadronic cross section cannot be measured reliably and a sum of exclusive states must be used.Despite the large data set accumulated in the past years and the analysis studies performed, there is still a ∼3.5 sigma discrepancy between the SM calculation and the experimental value [2].Not all exclusive states have yet been measured, and new measurements will improve the reliability of the calculation.Finally, these studies provide information on the resonant spectroscopy.
In this paper, we report the first measurements of the 2(π + π − )3π 0 and 2(π + π − )2π 0 η channels.The final states are produced in conjunction with a hard photon, assumed to result from ISR.To reduce background from Υ (4S) decays, the analysis is restricted to the c.m. energy below 4.5 GeV.As part of the analysis, we search for and observe intermediate states, including the η, ω, and ρ resonances.In the charmonium region, we observe J/ψ and ψ(2S) signals in the studied final states and the corresponding branching fractions are measured.
The data used in this analysis were collected with the BABAR detector at the PEP-II asymmetric-energy e + e − storage ring.The total integrated luminosity used is 468.6 fb −1 [19], which includes data collected at the Υ (4S) resonance (424.7 fb −1 ) and at a c.m. energy 40 MeV below this resonance (43.9 fb −1 ).
The BABAR detector is described in detail elsewhere [20].Charged particles are reconstructed using a BABAR tracking system, which is comprised of a silicon vertex tracker (SVT) and a drift chamber (DCH), both located inside a 1.5 T solenoid.Separation of pions and kaons is accomplished by means of a detector of internally reflected Cherenkov light (DIRC) and energy-loss measurements in the SVT and DCH.Photons are detected in an electromagnetic calorimeter (EMC).Muon identification is provided by an instrumented flux return.
To evaluate the detector acceptance and efficiency, we have developed a special package of Monte Carlo (MC) simulation programs for radiative processes based on the approach of Kühn and Czyż [21].Multiple collinear softphoton emission from the initial e + e − state is implemented with a structure function technique [22,23], while additional photon radiation from final-state particles is simulated using the PHOTOS package [24].The precision of the radiative simulation is such that it contributes less than 1% to the uncertainty in the measured hadronic cross sections.
A sample of 100-200k simulated events is generated for each signal reaction and processed through the detector response simulation, based on the GEANT4 package [26].These events are reconstructed using the same software chain as the data.Variations in the detector conditions are taken into account.The simulation includes random trigger events to account for the observed distributions of the background tracks and photons.Most of the experimental events contain additional soft photons due to machine background or interactions in the detector material, which are properly modeled in the simulation.
For the purpose of background estimation, large samples of events from the main relevant ISR processes (4πγ, 5πγ, ωηγ, and 2(π + π − )π 0 π 0 γ) are simulated.The background from the relevant non-ISR processes, namely e + e − → qq (q = u, d, s) and e + e − → τ + τ − , are generated using the jetset [27] and koralb [28] programs, respectively.The cross sections for the above processes are known with an accuracy about or better than 10%, which is sufficient for the present purpose.

III. EVENT SELECTION AND KINEMATIC FIT
Candidates for the 2(π + π − )3π 0 γ and 2(π + π − )2π 0 ηγ events are selected by requiring that there be four well measured tracks and seven or more detected photons, with an energy above 0.02 GeV in the EMC.We assume that the photon with the highest energy is the ISR photon, and we require its c.m. energy to be larger than 3 GeV.
The four tracks must have zero total charge and extrapolate to within 0.25 cm of the beam axis and 3.0 cm of the nominal collision point along that axis.In order to recover a relatively small fraction of signal events that contain a background track from secondary decay or interaction, we allow for the presence of a fifth track in the event, which however must not fulfill the above condition.The four tracks that satisfy the extrapolation criteria are fit to a vertex to determine the collision point, which is used in the calculation of the photon directions.
We subject each candidate event to a set of constrained kinematic fits and use the fit results, along with chargedparticle identification, to select the final states of interest and evaluate backgrounds from other processes.The kinematic fits use the four-momenta and covariance matrices of the colliding electrons and selected tracks and photons.The fitted three-momenta of each track and photon are then used in further kinematic calculations.We exclude the photon with the highest c.m. energy, which is assumed to arise from ISR, and consider each independent set of six other photons, and combine them into three pairs.For each set of six photons, there are 15 independent combinations of photon pairs.We retain those combinations in which the diphoton mass of at least two pairs lies within 35 MeV/c 2 of the π 0 mass m π 0 .The selected combinations are subjected to a fit in which the diphoton masses of the two pairs with |m(γγ) − m π 0 | < 35 MeV/c 2 are constrained to m π 0 .In combination with the constraints due to four-momentum conservation, there are thus six constraints (6C) in the fit.The photons in the remaining ("third") pair are treated as being independent.If all three photon pairs in the combination satisfy |m(γγ) − m π 0 | < 35 MeV/c 2 , then we test all possible combinations, allowing each of the three diphoton pairs in turn to be the third pair, i.e., the pair without the m π 0 constraint.
The above procedure allows us not only to search for events with π 0 → γγ in the third photon pair, but also for events with η → γγ.
The 6C fit is performed under the signal hypothesis e + e − → 2(π + π − )π 0 π 0 γγγ ISR .The combination with the smallest χ 2 is retained, along with the obtained χ 2 4π2π 0 γγ value and the fitted three-momenta of each track and photon.Each selected event is also subjected to a 6C fit under the e + e − → 2(π + π − )π 0 π 0 γ ISR background hypothesis, and the χ 2 4π2π 0 value is retained.The 2(π + π − )π 0 π 0 process has a larger cross section than the 2(π + π − )3π 0 signal process and can contribute to the background when two background photons are present.

IV. ADDITIONAL SELECTION CRITERIA
We require the tracks to lie within the fiducial region of the DCH (0.45-2.40 radians) and to be inconsistent with being a kaon or muon.The photon candidates are required to lie within the fiducial region of the EMC (0.35-2.40 radians) and to have an energy larger than 0.035 GeV.The angular distance between the ISR photon and the closest track must be greater than 1 radian; this requirement significantly suppresses the non-ISR background, in particular reducing the background from e + e − → τ + τ − to a negligible level.A requirement that any extra photons in an event must have an energy below 0.7 GeV reduces the multi-photon background by 10-20%.Finally, the background from the ISR process e + e − → 2(π + π − )2π 0 γ is reduced from 30% to about 1-2%, with a loss of only 5% of signal events, by requiring χ 2 4π2π 0 > 30.

FIG. 3:
The MC-simulated distribution for e + e − → ηπ + π − π 0 π 0 events of (a) the third-photon-pair invariant mass m(γγ), and (b) m(γγ) vs m(2(π + π − )2π 0 γγ).shows the m(γγ) distribution for events in the signal region after the above requirements have been applied.The dip in this distribution at the π 0 mass value is a consequence of the kinematic fit constraint of the best two photon pairs to the π 0 mass, so the third photon pair is always formed from photon candidates that are less well measured.
Figure 2 shows the m(γγ) distribution vs the invariant mass m(2(π + π − )2π 0 γγ) for events in the signal (a) and control (b) region.Events from the e + e − → 2(π + π − )3π 0 and 2(π + π − )2π 0 η processes are clearly seen in the signal region, as well as J/ψ decays to these final states.In the control region no significant structures are seen; we use these events to evaluate background.
Our strategy to extract the signals for the e + e − → 2(π + π − )3π 0 and 2(π + π − )2π 0 η processes is to perform a fit for the π 0 and η yields in intervals of 0.05 GeV/c 2 in the distribution of the 2(π + π − )2π 0 γγ invariant mass.This mass interval is about three times wider than the experimental resolution.

A. Number of signal events in simulation
As mentioned in Sec.II, the model used in the MC simulation assumes that the seven-pion final state results from ωπ 0 η and ηπ + π − π 0 π 0 production, with ω decays to three pions and η decays to all modes.As shown below, these two final states dominate the observed cross section.Also, the ISR photon is simulated to wider angles than the EMC acceptance, reducing the nominal efficiency.For each mode we have 200,000 simulated events from the primary generator.
The selection procedure applied to the data is also applied to the MC-simulated events.Figures 3 and 4 show (a) the m(γγ) distribution and (b) the distribution of m(γγ) vs m(2(π + π − )2π 0 γγ) for the simulated ηπ + π − π 0 π 0 and ωπ 0 η events, respectively.The π 0 peak is not Gaussian in either reaction.Background photons are included in the simulation.Therefore, the simulation accounts for the combinatorial background that arises when background photons are combined with photons from the signal reactions.
The combinatorial background is subtracted using the data from the χ 2 control region.We do not know how large the combinatorial background is in the signal region, and we use a scale factor varying from 1.0 to 1.5 for the subtraction to estimate the uncertainty in the number of signal events.The method is illustrated using simulation in Fig. 5, which shows the m(γγ) distribution with a bin width of 0.02 GeV/c 2 .The solid histograms show the simulated results from the signal region after subtraction of the simulated combinatorial background with the scale factor 1.5.The sum of three Gaussian functions is used to describe the π 0 signal shape.A third-order polynomial function is used to describe the shape of the remaining combinatorial background.The fitted function is shown by the smooth solid curve, while the dashed curve is for the contribution of the remaining combinatorial background.The remaining combinatorial background contribution is almost negligible for the scale factor value 1.5.We obtaine 1122±46 and 1161±55 simulated signal events for each mode, respectively.If the scale factor 1.0 is used, the remaining background is well described by the polynomial function and the signal yield does not change by more than 3%.
Alternatively, for the ηπ + π − π 0 π 0 events, we determine the number of events by fitting the η signal from the η → π + π − π 0 decay: the simulated distribution is shown in Fig. 6(a) (twelve entries per event).The fit functions are again the sum of three Gaussian functions FIG.6: (a) The π + π − π 0 invariant mass distribution for the MC-simulated e + e − → ηπ + π − π 0 π 0 events.The dashed curve is for the combinatorial background.(b) The π + π − π 0 invariant masses for the MC-simulated e + e − → ωπ 0 η events.The histogram shows the π + π − π 0 combination closest to the η mass, while the two remaining combinations (dots) exhibit the ω meson.The curves show the fit functions used to obtain the number of signal (solid) events, and the combinatorial background (dashed) contribution.and a polynomial for the combinatorial background.In total we obtain 1183±49 events.A similar fit of the η signal is performed for the ωπ 0 η final state simulation with 1110±54 selected events.
Similarly, as an alternative for the ωπ 0 η events, the ω mass peak can be used.To reduce the number of combinatorial entries, we require one π + π − π 0 combination to have invariant mass close to the η mass, and fit the remaining two combinations to extract the numbers of signal events with an ω, as shown in Fig. 6(b).In total 1104±71 signal events are found.A Breit-Wigner (BW) function, convolved with a Gaussian distribution to account for the detector resolution, is used to describe the ω signal.A second-order polynomial is used to describe the background.

B. Efficiency evaluation
The mass-dependent detection efficiency is obtained by dividing the number of fitted MC events in each 0.05 GeV/c 2 mass interval of the hadronic system by the number generated in the same interval.The number of signal events in the simulation, obtained by fitting the π 0 , η, or ω signals, is consistent within uncertainties not only in total, but also in every mass interval.We do not see any significant difference in mass-dependent efficiency between the different methods.The uncertainty in the value of the efficiency in each mass bin is dominated by the fluctuation of the combinatorial background.We average the five efficiencies in each 0.05 GeV/c 2 mass interval and fit the result with a third-order polynomial function, shown in Fig. 7.The result of this fit is used for the cross section calculation.
Although the signal simulation accounts for all η decay modes, the efficiency calculation considers the signal η → π + π − π 0 decay mode only.This efficiency estimate takes into account the geometrical acceptance of the detector for the final-state photons and the charged pions, the inefficiency of the detector subsystems, and the event loss due to additional soft-photon emission from the initial and final states.Corrections that account for data-MC differences are discussed below.
From Fig. 7 it is seen that the reconstruction efficiency is about 2.7%, roughly independent of mass.By comparing the results of the five different methods used to evaluate the efficiency, we conclude that the relative overall efficiency does not change by more than 5% because of variations of the functions used to extract the number of events or the use of different models.This value is taken as an estimate of the systematic uncertainty in the acceptance associated with the simulation model used and with the fit procedure.
We do not simulate the 2(π + π − )η and 2(π + π − )2π 0 η intermediate states, which are observed in data (see below) in the η → 3π 0 and η → γγ decays.But our previous studies [13] have demonstrated that, for these and similar decays, the variations in efficiency due to model dependence do not exceed 5%.In combination with the selections above, we assign 7% as a sistematic uncertainty to the detection efficiency.

VI. THE 2(π
The solid histogram in Fig. 8 (a) shows the same m(γγ) distribution of Fig. 1 (b) binned in mass intervals of 0.02 GeV/c 2 .The dashed histogram corresponds instead to the distribution of data from the χ 2 control region, and the dotted histogram is the estimated remaining background from e + e − → 2(π + π − )π 0 π 0 events produced via ISR.No evidence for a peaking background is seen in either of the two background distributions.We subtract the background evaluated using the χ 2 control region with the scale factor 1.0, and vary it to 1.5 to check the stability of the result.The resulting m(γγ) distribution is shown in Fig. 8 (b).
We fit the data of Fig. 8 (b) with a combination of a signal function, taken from simulation, and a background function, taken to be a third-order polynomial.The fit is performed in the m(γγ) mass range from 0.0 to 0.5 GeV/c 2 .The result of the fit is shown by the solid and dashed curves.A total of 12 559 ± 174 events is obtained.Note that this number includes a relatively small peaking background component, due to qq events, which is discussed in Sect.VI B. The same fit is applied to the corresponding m(γγ) distribution in each 0.05 GeV/c 2 interval in the 2(π + π − )2π 0 γγ invariant mass.The resulting number of 2(π + π − )3π 0 event candidates as a function of m(2(π + π − )2π 0 γγ), including the peaking qq background, is reported in Fig. 9 (a).
To normalize the uds simulation, we calculate the diphoton invariant mass distribution of the ISR candidate with all the remaining candidate photon in the event.A π 0 peak is observed, with approximately the same number of events in data and simulation, leading to a normalization factor of 1.0 ± 0.1.The resulting uds background is shown by the squares in Fig. 9 (b), the uds background is negligible below 2 GeV/c 2 , but accounts for more than half of the total event yield around 4 GeV/c 2 and above.

C. Cross section for e
where E c.m. is the invariant mass of the seven-pion system, dN 7πγ is the background-subtracted number of selected events in the interval dE c.m. , and MC 7π (E c.m. ) is the corresponding detection efficiency from simulation.The factor corr 7π accounts for the difference between data and simulation: the MC efficiency is larger by (1.0±1.0)%/percharged track [7] and by (3.0±1.0)% per π 0 [12].The ISR differential luminosity [10], dL, is calculated using the total integrated BABAR luminosity of 469 fb −1 [19].The initial-and final-state soft-photon emission is accounted for by the radiative correction factor (1 + δ R ), which is close to unity for our selection criteria.The cross section results contain the effect of vacuum polarization because this effect is not accounted for in the luminosity calculation.
Our results for the e + e − → 2(π + π − )π 0 π 0 π 0 cross section are shown in Fig. 11.The cross section does not exhibit any clear structures except signals from the J/ψ and ψ(2S) resonances.Because we present our data in bins of width 0.050 GeV/c 2 , compatible with the exper-  imental resolution, we do not apply an unfolding procedure to the data.Numerical values for the cross section are presented in Table I.The J/ψ region is discussed below.

D. Summary of systematic uncertaintes
The systematic uncertainties, presented in the previous sections, are summarized in Table II, along with the corrections that are applied to the measurements.
The three corrections applied to the cross sections sum up to 14.5%.The systematic uncertainties are considered to be uncorrelated and are added in quadrature, summing to 10%.The largest systematic uncertainty arises from the fitting and background subtraction procedures.It is estimated by varying the background levels and the parameters of the functions used.
The distribution of the π + π − π 0 invariant mass (twelve entries per event) is shown in 13(a).Prominent η and ω peaks are seen.The scatter plot in Fig. 13(b) shows one π + π − π 0 vs another π + π − π 0 invariant mass for the same event.Correlated η and ω production from e + e − → ωπ 0 η is seen.A scatter plot of the π + π − π 0 vs the seven-pion mass is shown in Fig. 13(c).A clear signal for a J/ψ peak is also observed.
1.The η2(π + π − ) intermediate state To determine the contribution of the η2(π + π − ) intermediate state, we fit the events of Fig. 12(a) using a triple-Gaussian function to describe the signal peak, as in Fig. 6(a), and a polynomial to describe the background.The result of the fit is shown in Fig. 15(a).We obtain 1410 ± 58 η2(π + π − ) events.The number of η2(π + π − ) events as a function of the seven-pion invariant mass is determined by performing an analogous fit in each 0.05 GeV/c 2 interval of m(2(π + π − )3π 0 ).The resulting distribution is shown in Fig. 16.
The very rich intermediate structures in the η2(π + π − ) mode were carefully studied in our previous paper [11] with significantly larger statistical precision.
Using Eq. ( 1), we determine the cross section for the e + e − → η2(π + π − ) process.Our simulation takes into account all η decays, so the cross section results, shown in Fig. 17 and listed in Table III, correspond to all η decays.Systematic uncertainties in this measurement are  The mass distribution of the 2(π + π − )3π 0 events in the ω peak (circles) correlated with η production in comparison with all 2(π + π − )3π 0 events(squares).
the same as those listed in Table II.Figure 17 shows our measurement in comparison to our previous result [11] and to those from the CMD-3 experiment [32].These previous results are based on different η decay modes than those considered here.The different results are seen to agree within the uncertainties.Including the results of the present study, we have thus now measured the e + e − → η2(π + π − ) cross section in three different η decay BABAR in η → γγ (squares) [13] and from SND (triangles) [29] in η → 3π 0 . modes.
The solid histogram in Fig. 18(a) shows the mass distribution of the π + π − π 0 combination closest to the nominal η mass, while the dotted histogram reports the invariant mass of the remaining two combinations of three pions after selecting the first combination within a window of ±80 MeV/c 2 from the nominal η mass.A fit to the dotted distribution with a sum of a BW for the ω signal and a combinatorial background, as shown in Sect.V, allows the extraction of the ωηπ 0 intermediate state signal, which amounts to 739±51 events.The contribution of the ωπ 0 η intermediate state to all 2(π + π − )3π 0 events is shown in Fig. 18(b).
Using Eq. ( 1), we determine the cross section for the e + e − → ωπ 0 η process.The energy dependence of the  FIG.20: (a) The π + π − π 0 invariant mass for data with the fit function for the ω signal (solid) plus the combinatorial background (dashed curve).The solid histogram shows peaking background from the simulated e + e − → ωη ISR events.(b) The mass distribution of the 2(π + π − )3π 0 events in the ω peak (circles) and estimated contribution from the ωη background (triangles), from ωπ 0 η (up-down triangles), and from uds (squares).
cross section is shown in Fig. 19 by the dots: we are in agreement with our previous measurement [13] and still slightly below the SND result [29].The numerical values of the cross section are listed in Table IV.Again, we have the measurements of this reaction in three different decay modes of η.FIG.21: The energy dependent e + e − → ωπ + π − π 0 π 0 cross section in the 2(π + π − )3π 0 mode (the J/ψ signal is off-scale).

The ωπ
To determine the contribution of the ωπ + π − π 0 π 0 intermediate state, we fit the events of Fig. 13(a) using a BW function to model the signal and a polynomial to model the background.The BW function is convolved with a Gaussian distribution that accounts for the detector resolution, as described for the fit of Fig. 6(b).The result of the fit is shown in Fig. 20(a).We obtain 7808±176 ωπ + π − π 0 π 0 events.The number of ωπ + π − π 0 π 0 events as a function of the seven-pion invariant mass is determined by performing an analogous fit in each 0.05 GeV/c 2 interval of m(2(π + π − )3π 0 ).The resulting distribution is shown by the circle symbols in Fig. 20(b).
For the e + e − → ωπ + π − π 0 π 0 channel, there is a peaking background from e + e − → ωη when ω and η decay to π + π − π 0 .A simulation of this reaction with proper normalization leads to the peaking-background estimation shown by the histogram in Fig. 20(a) and by the triangle symbols in Fig. 20(b).We also have peaking background from the general uds reactions (also shown in Fig. 20(b)).
These background contributions, as well as the events  from the correlated ω and η production from the ωπ 0 η final state, are subtracted from the ωπ + π − π 0 π 0 signal candidate distribution.
The resulting e + e − → ωπ + π − π 0 π 0 cross section, corrected for the ω → π + π − π 0 branching fraction, is shown in Fig. 21 and tabulated in Table V.The uncertainties are statistical only.The systematic uncertainties are about 10%.No previous measurement exists for this process.The cross section exhibits a rise at threshold, a decrease at large E c.m. , and a possibly resonant activity at around 2.3-2.5 GeV.The result by CMD-3 for the significantly lower e + e − → ωπ + π − π + π − cross section [32] is shown by squares.
4. The ηπ + π − π 0 π 0 intermediate state A similar approach is used to determine the contribution of the ηπ + π − π 0 π 0 intermediate state.We fit the events of Fig. 13(a) using the three-Gaussian function for the signal and a polynomial to model the background.The result of the fit is shown in Fig. 22(a).The fitted ηπ + π − π 0 π 0 yield corresponds to 2522 ± 91 events.The signal distribution as a function of the seven-pion invariant mass is determined by performing an analogous fit in each 0.05 GeV/c 2 interval of m(2(π + π − )3π 0 ), and is shown by the circle symbols in Fig. 22(b).
Also in this case a peaking background arises from the process e + e − → ωη when ω and η decay to π + π − π 0 .Its contribution, estimated with MC simulation, is shown by the histogram in Fig. 22(a) and by the triangle symbols in Fig. 22(b).
We also have peaking background from the general uds reactions, shown by squares in Fig. 22(b).And finally, we remove events from the ωπ 0 η final state (up-down triangles).
The e + e − → ηπ + π − π 0 π 0 cross section, corrected for the η → π + π − π 0 branching fraction, is shown in Fig. 23 and tabulated in 0.1 GeV bins in Table VI.The uncertainties are statistical only.The systematic uncertainties are about 10%.We are in good agreement with a recent measurement of this cross section [13] in the η → γγ decay mode.
5. The ρ(770) ± π ∓ π + π − π 0 π 0 intermediate state A similar approach is followed to study events with a ρ ± meson in the intermediate state.Because the ρ meson is broad, a BW function is used to describe the signal shape.There are twelve candidate ρ ± entries per event, leading to a large combinatorial background.To extract the contribution of the ρ ± π ∓ π + π − π 0 π 0 intermediate state we fit the events in Fig. 14(a) with a BW function to describe the signal and a polynomial to describe the background.The parameters of the ρ resonance are taken from Ref. [25].The result of the fit is shown in Fig. 24(a).We obtain 9138±371 ρ ± π ∓ π + π − π 0 π 0 events.The distribution of these events vs the seven-pion invariant mass is shown by the circle symbols in Fig. 24(c), while a similar fit for the uds simulation is shown by squares.The uds background is dominant in all energy regions except for J/ψ and ψ(2S).
In these events more than one ρ ± per event can be expected, indicating a significant production of J/ψ → ρ + ρ − π + π − π 0 .To determine the rate of ρ + ρ − π + π − π 0 events in the J/ψ decays, we perform a fit to determine the number of ρ + in intervals of 0.02 GeV/c 2 in the π − π 0 distribution of Fig. 14(b) for events within ±0.1 GeV/c 2 of the J/ψ mass.The result is shown in Fig. 24(b).Indeed, a small ρ + peak with 415 ± 340 events is observed, compared to 2844 events in the J/ψ peak region, corresponding to about 20% of all decays with one or two ρ ± .However, the uncertainty in this estimate is almost at the same level.
The charmonium region for all intermediate states is discussed below.

F. The sum of intermediate states
We consider whether the 2(π + π − )3π 0 channel contains other intermediate state contributions.The circle symbols in Fig. 25 show the total number of 2(π + π − )3π 0 events, repeated from Fig. 9.We perform a sum of the number of η2(π + π − ), ωπ 0 η, ηπ + π − π 0 π 0 , ωπ + π − π 0 π 0 , and ρ ± π ∓ π + π − π 0 π 0 intermediate state events, found as described in the previous sections, and show this sum by the square symbols in Fig. 25.Based on the results of our study of correlated ρ + ρ − production, we scale the num-  ber of events found from the fit to the ρ peak so that it corresponds to the number of events with either a single ρ ± or with a ρ + ρ − pair.This summed curve is seen to be in agreement with the total number of 2(π + π − )3π 0 events; we conclude there is no significant contribution from other (unobserved) intermediate states.The analogous approach to that described above for e + e − → 2(π + π − )π 0 π 0 π 0 events is used to study e + e − → 2(π + π − )π 0 π 0 η events.We fit the η signal in the thirdphoton-pair invariant-mass distribution (cf., Fig. 1) with the sum of two Gaussians with a common mean, while the relatively smooth background is described by a secondorder polynomial function, as shown in Fig. 26(a).We obtain 1651 ± 50 events.Figure 26(b) shows the mass distribution of these events.The invariant-mass distribution for the 2(π + π − )2π 0 η events obtained from the η signal fit.The contribution of the uds background events is shown by the squares.

B. Peaking background
The major background producing an η peak is the non-ISR background, in particular e + e − → 2(π + π − )π 0 π 0 π 0 η when one of the neutral pions decays asymmetrically, producing a photon interpreted as ISR.The η peak from the uds simulation is visible in Fig. 10.We fit the η peak in the uds simulation in intervals of 0.05 GeV/c 2 in m(2(π + π − )π 0 π 0 γγ).The results are shown by the squares in Fig. 26 (b).
To normalize the uds simulation, we form the diphoton invariant-mass distribution of the ISR candidate with all the remaining photons in the event.Comparing the number of events in the π 0 peaks in data and uds simulation, we assign a scale factor of 1.5 ± 0.2 to the simulation and subtract these events from the data distribution.
C. Cross section for e + e − → 2(π + π − )π 0 π 0 η The cross section for e + e − → 2(π + π − )π 0 π 0 η is determined using Eq. ( 1).The results are shown in Fig. 27 and listed in Table VII.These are the first results for this process.The systematic uncertainties and corrections are the same as those presented in Table II except that the uncertainty in the detection efficiency increases to 13%.
The cross section is approximately zero until well above 2 GeV and so is not useful in the vacuum polarization calculations; we have not yet performed a study of intermediate states for this process.

VIII. THE J/ψ REGION
A. The 2(π + π − )3π 0 final state Figure 28(a) shows an expanded view of the charmonium region from Fig. 9, which has large contributions from the J/ψ and ψ(2S) decays to seven pions.The nonresonant background distribution is flat in this region.
The observed peak shapes are not purely Gaussian because of radiation effects and resolution, as is also seen in the simulated signal distributions shown in Fig. 28(b).The sum of two Gaussians with a common mean is used to describe each peak.We obtain 3391 ± 101 J/ψ events and 290 ± 40 ψ(2S) events.Using these results for the number of events, the detection efficiency, and the ISR luminosity, we determine the product: = (345 ± 10 ± 50) eV , where Γ J/ψ ee is the electronic width, dL/dE = 180 nb −1 / MeV is the ISR luminosity at the J/ψ mass m J/ψ , MC = 0.027 ± 0.002 is the detection efficiency from simulation, corr = 0.85 is the correction, discussed in Sec.VI D, and C = 3.894 × 10 11 nb MeV 2 is a conversion constant [25].We estimate the systematic uncertainty for this region to be 15%.The subscript "7π" for the branching fraction refers to the 2(π + π − )3π 0 final state exclusively.
Note that the J/ψ decay to the ωπ 0 η mode is almost ten times smaller, (0.34 ± 0.17) × 10 −3 [13], and cannot be extracted from our data.where the third uncertainty is associated with the procedure used to determine the correlated ρ + ρ − rate.No other measurements for these processes exist.
There are no previous results for these final states.

C. Summary of the charmonium region study
The rates of J/ψ and ψ(2S) decays to 2(π + π − )3π 0 , 2(π + π − )2π 0 η and several intermediate final states have been measured.The measured products and calculated branching fractions are summarized in Table VIII together with the available PDG values for comparison.

IX. SUMMARY
The excellent performance of the BABAR detector for photon energy and charged-particle resolution, together with its strong particle identification capabilities, allow the reconstruction of the 2(π + π − )3π 0 and 2(π + π − )2π 0 η final states from threshold up to 4.5 GeV via the ISR process.
The analysis shows that the effective luminosity and efficiency have been understood with 10-13% accuracy.The cross section measurements for the e + e − → 2(π + π − )3π 0 and the e + e − → 2(π + π − )2π 0 η reactions has been measured for the first time.
The selected multi-hadronic final states in the broad range of accessible energies provide new information on hadron spectroscopy.
The initial-state radiation events also allow a study of J/ψ and ψ(2S) production and a measurement of the corresponding products of the decay branching fractions and e + e − width for most of the studied channels, the majority of them for the first time.

FIG. 1 :
FIG. 1: (a) Distribution of the invariant mass m(γγ) of the third photon pair vs χ 2 4π2π 0 γγ .The lines define the boundaries of the signal and control regions.(b) Distribution of m(γγ) in the signal region χ 2 4π2π 0 γγ < 50 with the additional selection criteria described in the text.

Figure 1 (
Figure1 (b)  shows the m(γγ) distribution for events in the signal region after the above requirements have been applied.The dip in this distribution at the π 0 mass value is a consequence of the kinematic fit constraint of the best two photon pairs to the π 0 mass, so the third photon pair is always formed from photon candidates that are less well measured.Figure2shows the m(γγ) distribution vs the invariant mass m(2(π + π − )2π 0 γγ) for events in the signal (a) and control (b) region.Events from the e + e − → 2(π + π − )3π 0 and 2(π + π − )2π 0 η processes are clearly seen in the signal region, as well as J/ψ decays to these final states.In the control region no significant structures are seen; we use these events to evaluate background.Our strategy to extract the signals for the e + e − → 2(π + π − )3π 0 and 2(π + π − )2π 0 η processes is to perform a fit for the π 0 and η yields in intervals of 0.05 GeV/c 2 in the distribution of the 2(π + π − )2π 0 γγ invariant mass.This mass interval is about three times wider than the experimental resolution.

FIG. 7 :
FIG. 7:The energy-dependent reconstruction efficiency for e + e − → 2(π + π − )3π 0 events, determined using five different methods: see text.The curve shows the results of a fit to the average values, which is used in the cross section calculation.

FIG. 8 :
FIG. 8: (a) Invariant mass m(γγ) for data in the χ 2 signal (solid) and control (dashed) regions.The dotted histogram shows the estimated background from e + e − → 2(π + π − )π 0 π 0 .(b) The m(γγ) invariant mass for data after the background subtraction.The curves are the fit results as described in the text.

FIG. 9 :
FIG. 9: (a) The invariant-mass distribution of π + π − 3π 0 events, obtained from the fit to the π 0 mass peak.(b) Expanded view of (a) to show the contribution from non-ISR uds background, shown by squares.

FIG. 15 :
FIG.15:(a)The 3π 0 invariant mass for data.The curves show the fit functions.The solid curve shows the η peak (based on MC simulation) plus the non-η continuum background (dashed).

2 )FIG. 18 :
FIG.18: Mass plots for the ωπ 0 η intermediate state: (a) The solid histogram is for the π + π − π 0 invariant mass closest to the η mass, while the dots are for the two remaining combinations of π + π − π 0 .The solid curve shows the fit function for the ω signal plus the combinatorial background (dashed curve).(b) The mass distribution of the 2(π + π − )3π 0 events in the ω peak (circles) correlated with η production in comparison with all 2(π + π − )3π 0 events(squares).

FIG. 24 :
FIG.24: Mass distributions for the ρ ± π ∓ π + π − π 0 π 0 intermediate state: (a) The π ± π 0 invariant mass for data.The dashed curve shows the fit to the combinatorial background.The solid curve is the sum of the background curve and the BW function for the ρ ± .(b) The result of the ρ + fit in bins of 0.02 GeV/c 2 in the ρ − mass.(c) Number of events in bins of Ec.m. from the ρ ± → π ± π 0 (circles) intermediate states.The squares show the event numbers obtained from uds production.

2 )FIG. 25 :
FIG.25:The 2(π + π − )2π 0 γγ mass distribution summed over the intermediate states.The circles show the number of events, determined from the π 0 fit.The squares show the sum of events with η, ω, and ρ production, the latter corrected for the ρ + ρ − production.

FIG. 26 :
FIG. 26:Mass distributions for the 2(π + π − )2π 0 η final state.(a) The third-photon-pair invariant mass for data.The dashed curve shows the fitted background.The solid curve shows the sum of background and the two-Gaussian fit function used to obtain the number of events with an η.(b) The invariant-mass distribution for the 2(π + π − )2π 0 η events obtained from the η signal fit.The contribution of the uds background events is shown by the squares.

FIG. 30 :
FIG.30:The seven-pion invariant mass for events with a three-pion invariant mass in the ω(792) (a) or η (b) mass regions.The curves show the fit functions described in the text.

TABLE III :
Summary of the e + e − → η2(π + π − ) cross section measurement.The uncertainties are statistical only.

TABLE IV :
Summary of the e + e − → ωπ 0 η cross section measurement.The uncertainties are statistical only.