Study of charmonia in four-meson final states produced in two-photon collisions

We report measurements of charmonia produced in two-photon collisions and decaying to four-meson final states, where the meson is either a charged pion or a charged kaon. The analysis is based on a 395fb^{-1} data sample accumulated with the Belle detector at the KEKB electron-positron collider. We observe signals for the three C-even charmonia eta_c(1S), chi_{c0}(1P) and chi_{c2}(1P) in the pi^+pi^-pi^+pi^-, K^+K^-pi^+pi^- and K^+K^-K^+K^- decay modes. No clear signals for eta_c(2S) production are found in these decay modes. We have also studied resonant structures in charmonium decays to two-body intermediate meson resonances. We report the products of the two-photon decay width and the branching fractions, Gamma_{gamma gamma}B, for each of the charmonium decay modes.


Abstract
We report measurements of charmonia produced in two-photon collisions and decaying to fourmeson final states, where the meson is either a charged pion or a charged kaon. The analysis is based on a 395 fb −1 data sample accumulated with the Belle detector at the KEKB electron-positron collider. We observe signals for the three C-even charmonia η c (1S), χ c0 (1P ) and χ c2 (1P ) in the π + π − π + π − , K + K − π + π − and K + K − K + K − decay modes. No clear signals for η c (2S) production are found in these decay modes. We have also studied resonant structures in charmonium decays to two-body intermediate meson resonances. We report the products of the two-photon decay width and the branching fractions, Γ γγ B, for each of the charmonium decay modes.

I. INTRODUCTION
The two-photon decay widths (Γ γγ ) of charge-conjugation (C)-even charmonium states are important observables that are sensitive to the properties of cc quarks inside charmonium bound states. Precise measurements of widths can give valuable constraints on the models that describe the nature of heavy quarkonia [1][2][3].
High luminosity electron-positron colliders are well suited to measurements of the twophoton production of charmonium states, since they provide a large flux of quasi-real photons colliding at two-photon center-of-mass energies covering a wide range. Various charmonium states, including the η c (1S), χ c0 , χ c2 , η c (2S) and χ c2 (2P ), have been observed so far in two-photon production processes.
In resonance formation, the production cross section from quasi-real two-photon collisions at an e + e − collider is proportional to the product of the two-photon decay width of the resonance and its branching fraction (B) to the observed final state: where R is the resonance. For the above charmonium states, the two-photon decay widths are mainly determined from two-photon collision measurements performed in a small number of specific decay modes only since the production rate of the corresponding states is limited. The results have rather large statistical errors and in some cases even larger systematic errors originating from uncertainties in the corresponding charmonium branching fraction to the observed final state. It is necessary to reduce these errors with high-statistics measurements in various decay channels. Through these measurements, we can also obtain ratios of branching fractions for various charmonium decay modes.
In this paper, we report the measurements of production of the charmonium states η c (1S), χ c0 and χ c2 in two-photon collisions from their decays to four-meson final states π + π − π + π − , K + K − π + π − and K + K − K + K − , with the Belle detector. We have studied intermediate meson-resonant structures in these decay modes. We have also searched for production of the η c (2S) in these final states. Previous measurements of the same processes are found in Refs. [4][5][6][7][8][9][10]. Belle previously reported measurements of two-photon production of these charmonia in several decay modes: χ c2 → γJ/ψ [11], χ c0 and χ c2 decays to the two-meson final states π + π − and K + K − [12] as well as K 0 S K 0 S [13], and η c (1S) → pp [14]. Production of the η c (2S) in a two-photon process was observed in the K 0 S K ∓ π ± mode by the CLEO and BaBar experiments [15,16]. This is the only decay mode so far directly observed for this charmonium state.

A. Data and the Belle detector
We use data that corresponds to an integrated luminosity of 395 fb −1 recorded with the Belle detector at the KEKB asymmetric-energy e + e − collider [17]. Since the beam energy dependence of two-photon processes is very small, we combine the on-resonance and off-resonance data samples; the off-resonance data were taken 60 MeV below the Υ(4S) ( √ s = 10.58 GeV). The analysis is made in the "zero-tag" mode, where neither the recoil electron nor positron is detected. We restrict the virtuality of the incident photons to be small by imposing a strict transverse-momentum balance along the beam axis for the finalstate hadronic system.
A comprehensive description of the Belle detector is given elsewhere [18]. We mention here only the detector components essential to the present measurement. Charged tracks are reconstructed from hit information in a central drift chamber (CDC) located in a uniform 1.5 T solenoidal magnetic field. The z axis of the detector and the solenoid are along the positron beam direction, with the positrons moving in the −z direction. The CDC measures the longitudinal and transverse momentum components (along the z axis and in the rϕ plane, respectively). Track trajectory coordinates near the collision point are provided by a silicon vertex detector (SVD). Photon detection and energy measurements are performed with a CsI(Tl) electromagnetic calorimeter (ECL). Kaons are identified using information from the time-of-flight counters (TOF) and silica-aerogel Cherenkov counters (ACC). The ACC provides good separation between kaons and pions or muons at momenta above 1.2 GeV/c. The TOF system consists of a barrel of 128 plastic scintillation counters, and is effective in K/π separation mainly for tracks with momentum below 1.2 GeV/c. Lower energy kaons are also identified using specific ionization (dE/dx) measurements in the CDC. The magnet return yoke is instrumented to form the K L and muon detector (KLM), which detects muon tracks and provides trigger signals.
Signal events are efficiently triggered by several kinds of track-triggers that require two or more CDC tracks with combinations of TOF hits, ECL clusters or summed energies of ECL clusters. No additional cuts on the trigger information are applied. The trigger conditions are complementary for the detection of four-prong events; we obtain an overall trigger efficiency of ∼ 95%.
K/π separation uses a likelihood ratio formed from ACC, TOF and CDC information. No explicit lepton identification requirement is applied, since leptons are not a large background source in exclusive four-prong events.

B. Event selection
Candidate events are selected as follows; all variables in selection criteria (1)-(10) are measured in the laboratory frame. (1) There are exactly four charged particles and each satisfies the following criteria: p t > 0.1 GeV/c, dr < 5 cm, |dz| < 5 cm, where p t is the transverse momentum of a track with respect to the z axis, and dr and dz are the radial and axial distances, respectively, of the closest approach (as seen in the rϕ plane) to the nominal collision point; (2) the net charge of the four tracks is zero; (3) at least two of the four tracks satisfy the following additional criteria: p t > 0.4 GeV/c, dr < 1 cm, and −0.8660 < cos θ < +0.9563, where θ is the polar angle; in addition, the following two criteria are applied for the four-track systems: (4) the scalar sum of the momenta ( i |p i |) of the four tracks must be smaller than 6 GeV/c; (5) the average of the dz values for the four tracks is within 3 cm of the nominal collision point.
In order to reject backgrounds from single-photon annihilation or radiative events, the following criteria are required: (6) the invariant mass of track combinations that satisfy criterion (3) should be smaller than 4.5 GeV/c 2 (here, zero mass is assigned to each of the charged tracks), and (7) the missing mass squared of the recoil against the track combination is larger than 2 (GeV/c 2 ) 2 .
Calorimetry requirements are also applied: (8) the total calorimeter energy in an event is smaller than 6 GeV; (9) there is no electromagnetic cluster unassociated with a track (i.e., a neutral cluster) whose energy exceeds 0.4 GeV; (10) there is no π 0 candidate whose transverse momentum is larger than 0.1 GeV/c.
After all these selection criteria are applied, the four tracks are transformed into the e + e − center-of-mass (c.m.) frame, and the vector sum of their transverse momenta with respect to the beam direction in the c.m. frame, Σp * t , is calculated. The variable Σp * t approximates the transverse momentum of the two-photon-collision system. We require (11) |Σp * t | < 0.1 GeV/c in order to enhance the number of events from quasi-real twophoton collisions. Figure 1 shows the transverse momentum (|Σp * t |) distribution near the selection region.
C. K/π separation We apply particle identification requirements for each track. A likelihood ratio, R = L K /(L K +L π ), is calculated, where the kaon and pion likelihoods, L K and L π , are determined from information provided by the ACC, TOF and the CDC dE/dx system. All tracks with R > 0.8 are assigned to be kaons. This requirement gives a typical identification efficiency of 90% with a probability of 3% for a pion to be misidentified as a kaon. All other tracks are treated as pions.
When a pair of pions (kaons) in the four prong decay of a certain charmonium are misidentified as a pair of kaons (pions), they could produce a broad enhancement in the spectrum of the misidentified channel at a mass very different from that of the original charmonium. However, these backgrounds do not make any pronounced peaks near the three charmonia masses.

D. Rejection of ISR events
Two-photon events with relatively high γγ c.m. energies (above 3.2 GeV in the case of Belle) are contaminated by background processes where a real or virtual photon is emitted in the direction of the incident electron. They are called ISR events or pseudo-Compton events, in case the positron collides with the electron (e + e − → (e + e − * )γ → hadrons + γ) or with the virtual photon (e + e − → (e + γ * )e − → hadrons + e + e − ), respectively. Such events have a kinematical correlation between the p z -component in the laboratory frame (Σp z ) and the invariant mass of the final-state hadron system (W ). We reject them with the following additional requirement: (13) Σp z > (W 2 − 49 GeV 2 /c 4 )/(14 GeV/c 3 ) + 0.6 GeV/c. These backgrounds are not harmful in the measurements of C-even charmonia, since at lowest order they produce only C-odd hadron systems. We show a two-dimensional plot for Σp z vs invariant mass of the hadronic system for the four-pion final state candidate in Fig. 2.

A. Yields of the charmonia
We measure the charmonium yields from the invariant mass distributions in each of the three final-state processes. The distributions are shown in Fig. 3. Peaks at the Ceven charmonium masses are assumed to be from two-photon production processes if the transverse-momentum distribution peaks in the vicinity of |Σp * t | = 0. We discuss the |Σp * t | distributions in Sect. V.C.
We find clear enhancements for the η c (1S) (denoted hereafter as "η c ") at ∼ 2.98 GeV/c 2 , χ c0 at ∼ 3.41 GeV/c 2 , and χ c2 at ∼ 3.555 GeV/c 2 in all the final states. We do not find any clear η c (2S) signals, which would appear at around 3.63-3.65 GeV/c 2 in the invariant mass distributions. A small peak near 3.69 GeV/c 2 in the 4π final state is from ψ(2S) → J/ψπ + π − , J/ψ → l + l − , where the leptons are not identified and are treated as pions. These events are attributed primarily to the double ISR e + e − annihilation process, e + e − → ψ(2S)γγ. Upper limits for η c (2S) production will be derived in Sect. VII.
The invariant-mass distribution in the vicinity of each charmonium peak is fitted to the sum of charmonium and background components using χ 2 minimization. We use a secondorder polynomial for the background component. The fitted regions are 2.8-3.2 GeV/c 2 , 3.3-3.5 GeV/c 2 and 3.5-3.6 GeV/c 2 for the analyses of the η c , χ c0 and χ c2 , respectively. The energy dependence of the efficiency is small within each resonance region and is neglected. We adopt the relativistic Breit-Wigner function for the η c and χ c0 signals, which are smeared with a Gaussian function whose width is fixed to the invariant-mass resolution (ranges in 7-10 MeV/c 2 depending on the decay process) estimated by the signal Monte Carlo (MC) simulation. In the fit to the χ c2 , we use a simple Gaussian function with a floating width, because the natural width is considerably smaller than the mass resolution. The mass resolutions from the fits, 8 -9 MeV/c 2 , are consistent with the MC expectations at the χ c2 mass.
The fit results are shown in Fig. 4 and are summarized in Table I. In the table, the masses are corrected for the effects of a systematic shift seen in the signal MC mainly due to the tails of kaon energy loss in the innermost detector region. The correction size is typically n K × (1 -2 MeV/c 2 ), where n K is the number of kaons. The systematic errors for the masses include the uncertainty due to this effect (a half of the correction size) as well as the uncertainty of the mass scale (2.0 MeV/c 2 ) originating from the uncertainty of the momentum scale, while those for the widths are obtained by changing the mass resolution by 1 MeV/c 2 . We confirm that the mass and width of the η c are robust for the change of the fit region as described in Sect. V.B.
We find a hint of a possible J/ψ contribution from the double ISR background in the 2K2π distribution. This could affect the yield and resonance parameters of the η c determined from the fit. We have tried a fit that includes a peak near the J/ψ mass. The changes of the η c yield (7%) and the total width (1.9 MeV) are taken into account in the systematic errors. There is no detectable shift in the mass. The fit including a J/ψ peak gives the J/ψ mass, 3093.6 ± 2.0 ± 2.4 MeV/c 2 , where we apply the same correction and systematic error as applied for η c → 2K2π. This provides a confirmation of the mass scale in the present measurement.
Possible interference of the charmonia with the continuum component is not taken into account in the present analysis. It is difficult to separate the continuum-component amplitudes that interfere with the charmonium amplitude in each process, because different partial-wave components of the continuum appear within the background. In general, the interference effects could give different mass/width fit results for each decay mode. However, for the three considered processes the values obtained for these parameters from the no-interference fits to the separate decay modes are consistent with each other. The mass values for the η c are systematically higher than the world average value, 2979.8 ± 1.2 MeV/c 2 [19]. The η c width and the parameters of χ c0 and χ c2 are consistent with the previous measurements. Table I also shows the combined results from the measurements of the three decay processes, where the central values are derived from the weighted means of the three fits and the systematic errors are evaluated considering the full correlations between the error sources.

B. Yields of two-body decays
We study the sub-structure of the charmonium signals by searching for quasi-two-body components, in which a charmonium meson decays into two resonances, and each resonance decays to two final-state mesons. We determine the yields of these two-body decays using the following procedure: we first make two-dimensional plots of the two invariant masses constructed from the available two-meson combinations of π + π − , K ± π ∓ or K + K − , for each event. The 4π and 4K samples have two entries from each event. Such plots are made in each charmonium signal region and sideband region. The charmonium sideband distribution is then subtracted from the signal distribution bin-by-bin in two dimensions. We note that the non-charmonium background components are considerable even in the vicinities of the charmonium peaks.
The signal regions are defined as the ranges within ±50 MeV/c 2 (±30 MeV/c 2 for the χ c2 ) of 2.98 GeV/c 2 , 3.41 GeV/c 2 , and 3.56 GeV/c 2 , for the η c , χ c0 and χ c2 , respectively. Next, we define the signal and sideband regions for the ρ 0 → π + π − , f 2 (1270) → π + π − , K * (892) 0 → K ± π ∓ and φ → K + K − mesons, based on their known masses and widths. (Hereafter, we refer to the f 2 (1270), K * (892) 0 and f ′ 2 (1525) as f 2 , K * 0 and f ′ 2 , respectively.) We apply the sideband subtraction at the f 2 side to extract the yield of the f 2 f ′ 2 decay. The signal regions are 0.64-0.88 GeV/c 2 , 1.08-1.40 GeV/c 2 , 0.80-0.96 GeV/c 2 and 1.00-1.04 GeV/c 2 for the ρ 0 , f 2 , K * 0 and φ, respectively. We take sideband regions on both sides adjacent to the signal region. In the φ → K + K − mode the sidebands cover the ranges from the mass threshold to 1.00 GeV/c 2 and 1.04-1.08 GeV/c 2 . For each of the other resonances the sideband has a width that is half the size of the signal region on each side. The sideband TABLE I: Results for the resonance parameters obtained from the fits to the invariant-mass distributions. The third column shows the Breit-Wigner width for the η c and χ c0 , and the Gaussian standard deviation for the χ c2 . The combined results from the three decay processes are also shown. The first and second errors for the masses and widths are statistical and systematic, respectively. The other errors are statistical only. subtraction is applied assuming the continuum yield has a linear dependence on invariant mass. Thus, we obtain a one-dimensional histogram corresponding to a resonance signal component tagged by another resonance. We search for ρ 0 , f 2 , K * 0 , φ and f ′ 2 → K + K − signals in the histograms. In the decays to pairs of same-kind resonances (such as f 2 f 2 or K * 0K * 0 ), we symmetrized the two-dimensional distribution in two directions, by adding it to the same distribution where the two-masses are interchanged. In this case we must divide the peak yield by two to obtain the correct signal yield and take into account the statistical correlation.
We fit the signal to a Breit-Wigner taking into account the angular momentum between the final state mesons and a linear continuum component in the one-dimensional distribution near the resonance masses. As the only exception, however, we count in φφ decays the number of events in the K + K − invariant mass range, 1.00-1.04 GeV/c 2 and subtract backgrounds from a sideband in the 1.04-1.24 GeV/c 2 region, instead of performing a fit. The one-dimensional distributions and the fits for the ρ 0 , f 2 , K * 0 and f ′ 2 resonances in different charmonium decay modes are shown in Figs. 5 -7. Table II summarizes the signal yields. We assume that there are no effects from interference in the sideband subtractions as well as the fits of the invariant mass distributions. We have applied corrections for inefficiencies arising from sideband subtractions in the twomeson resonance and charmonium regions; the targeted resonance component partially leaks outside the signal region and enters the sideband regions. The inefficiencies are calculated by assuming Breit-Wigner forms for the charmonia and light-quark resonances using the known masses and widths [19]. The 90%-confidence-level (C.L.) upper limits correspond to the upper edge of a 90% confidence region derived from Feldman and Cousins's table (Table X in Ref. [20]) assuming statistical fluctuations follow Gaussian distributions (there is a large number of events in both signal and sideband regions in the present case), with the expected number of signal events constrained to be non-negative.
We do not observe a significant ρ 0 ρ 0 component in η c decays in spite of the signal reported for η c → ρρ from measurements of radiative J/ψ decays [19]. Meanwhile, we find f 2 f 2 and f 2 f ′ 2 signals with greater than 4σ statistical significances. It can be noted that the observed yield for η c → f 2 f 2 → 4π is comparable to the total yield of η c → 4π. This indicates that the f 2 f 2 component dominates η c → 4π decays. We confirm the decays to f 2 f 2 in the distribution of the transversity angle ( Fig. 8(a)) defined as the angle between the decay planes of the two f 2 candidates in the η c rest frame. Here, the f 2 candidate is a π + π − combination whose invariant mass lies in the range 1.15-1.39 GeV/c 2 . The distribution shows a characteristic feature of tensor-meson decays. Similarly, η c → 2K2π is dominated by the sum of η c → K * 0K * 0 → 2K2π and η c → f 2 f ′ 2 → 2K2π. The transversity-angle distribution of the f 2 f ′ 2 candidates is shown in Fig. 8(b) (here, the f ′ 2 candidate is a K + K − combination whose invariant mass falls in 1.475 -1.575 GeV/c 2 ). The asymmetry with respect to 90 • can be interpreted as the effect of interference between the above two decay modes; an asymmetry in directional correlations between K + π + and K + π − should not arise in a pure f 2 f ′ 2 process. The (M(K − π + ), M(K + π − )) distributions expected from the f 2 f ′ 2 MC events overlap with the (K * 0 ,K * 0 ) mass region near its lowest mass-combination edge. The histograms are the distributions of the corresponding signal MC events generated assuming pure P-waves between the two tensor mesons. Figures 9 (a) and (b) show the invariant-mass distributions of 4π events that are consistent with ρ 0 ρ 0 and f 2 f 2 , respectively. We find no prominent peaks at the η c mass in the former, but find a clear signal for the η c in the latter. An enhancement at the η c mass is also clearly visible in Fig. 9(c) which shows the invariant-mass distribution for f 2 f ′ 2 → 2K2π candidates.

C. Yields of three-body decays
The distributions corresponding to the processes χ c2 → ρ 0 π + π − in Fig. 5(c) and χ c0 → K * 0 K − π + orK * 0 K + π − in Fig. 6(b) have significant net yields although there are no prominent structures from the two-body decay components, χ c2 → ρ 0 ρ 0 and χ c0 → K * 0K * 0 . We obtain the yields of each three-body decay using the following procedure. For the two processes, we first obtain the numbers of ρ 0 's in χ c2 → 4π and K * 0 's in χ c0 → 2K2π by fitting the invariant-mass distributions of π + π − ( Fig. 10(a), 4 entries/event) and K ± π ∓ ( Fig. 10(b), 2 entries/event) to functional forms corresponding to the sum of ρ 0 and f 2 resonances and a K * 0 resonance, respectively. Here, we also take into account the continuum component parameterized by a second-order polynomial. The results obtained are converted into the yields of three-body decay events after subtracting possible resonance yields from χ c2 → ρ 0 π + π − via ρ 0 ρ 0 and χ c0 → K * 0 K + π − (or c.c.) via K * 0K * 0 using the best estimates of the corresponding two-body yields after doubling the number of events to take account of the multiple entries. The results after these subtractions are summarized in Table II.

IV. DERIVATION OF PRODUCTS OF TWO-PHOTON DECAY WIDTH AND BRANCHING FRACTION
The product of the two-photon decay width and branching fraction is obtained from the following formula: where N is the number of observed events, ε the efficiency, and Ldt the integrated luminosity. F (M R , J) is a factor that is calculated from the two-photon luminosity function (L γγ (M R )), the charmonium mass (M R ) and the charmonium spin (J), using the relation that is valid when the resonance width is small compared to the available kinetic energy in the decay. The F parameter is calculated to be 2.11 fb/eV, 1.15 fb/eV, 4.74 fb/eV and 0.861 fb/eV for the η c , χ c0 , χ c2 and η c (2S), respectively, using the TREPS code [21]. We take the signal yields from Tables I and II. We subtract the small contributions from K 0 S K 0 S decays included in the χ c0 → 4π and χ c2 → 4π samples, 142±6 and 35±5 events respectively, which were determined from the two-dimensional invariant-mass plots mentioned in Sect. III.A. (χ c0 , χ c2 → K 0 S K 0 S results from our experiment are reported separately [13]. η c → K 0 S K 0 S is prohibited from parity conservation.) The efficiency in the detection of each process is determined from the signal MC events generated by the TREPS code [21] and processed by a full detector simulation. Threeand four-body decays were generated according to phase space with a resonance mass distributed according to the Breit-Wigner function, and the K * 0 (ρ 0 ) decaying to K ± π ∓ (π + π − ) isotropically.
is selected on the opposite side. We rescale the number of entries described in Sect. III.B by a factor of 1/2 to normalize the vertical axis to the number of events.
We rescale the number of entries described in Sect. III.B by a factor of 1/2 to normalize the vertical axis to the number of events, except for (d). For two-body decays involving resonances with spin, the decay amplitude takes into account correlations caused by spin, and the full matrix element is symmetrized with respect to identical final-state particles. We assume a pure state with the lowest possible orbital angular momentum that satisfies conservation laws in each decay channel. The helicity of the χ c2 along the incident axis is assumed to be 2, following the previous experiment sidebands, that is, contributions from ρ 0 π + π − etc. may also be included. results [11] and theoretical expectations [22].
The trigger efficiency is taken into account; it is estimated to be around 95% within the acceptance for all the processes. The detection efficiency including that of particle identification is calculated using MC events passed through the detector simulation code. Typical efficiencies are 16%, 10%, 5%, and 9% in the 4π final-state processes, 2K2π finalstate processes, 4K with a phase-space distribution, and the φφ process, respectively. Their dependence on the two-photon c.m. energy, or charmonium mass, is rather small.
We summarize the results on G ≡ Γ γγ B in Table III. The systematic errors and comparison to the previous experimental results also given there are described in detail in the following two sections. In some decay channels, we assume isospin invariance and use necessary branching fractions from Ref. [19]. The "isospin ·BF " factors in the table are the products of the isospin factor and the branching fraction(s) used to obtain the results.

A. Uncertainties in efficiencies
The trigger efficiency is estimated using the Monte-Carlo trigger simulator to be (95 ± 4)%, where the error is estimated from comparison of experimental and MC results in lowmultiplicity events and the efficiency variation under different beam background conditions.
We obtain a 5.5% systematic error from the uncertainty in track reconstruction efficiency for the four-track events, as well as 2% per kaon track from the kaon identification efficiency. The latter is determined from a study of kinematically identified kaons in a sample of D * + → D 0 π + , D 0 → K − π + decays. We check the kaon identification efficiency by loosening the 4K-event selection requirement; we temporarily require only three kaons in an event. We obtain (65±16)% more yield for the η c , χ c0 and χ c2 , on average with the looser condition. The charmonium peaks are less prominent with these selections, due to a larger contamination from non-4K events. The above ratio is consistent with expectations from the MC, ∼ 53%, for the average of the three charmonia decaying to the 4K state.
We use the total widths of the charmonia for the determination of the efficiency in the sideband subtractions. Taking systematic errors in the η c and χ c0 widths into account, we varied the widths of the η c and χ c0 from their nominal values by ±10 MeV and ±7 MeV, respectively, and assigned the obtained efficiency changes, ±12% and ±8%, respectively, as a systematic error. We take a 3% systematic error in the χ c2 case for the uncertainty in the treatment of its long Breit-Wigner tails.
Inefficiencies caused by the sideband subtractions of the isobar resonances decaying to two mesons are calculated using the resonance parameters. We estimated the systematic error from this effect by changing the integration range of the resonance function from ±8Γ to ±5Γ. It is estimated to be 5% for each two-body decay mode.

B. Systematic errors in signal yields
In the fits to the four-meson invariant-mass spectra described in Sect. III.A, we set the fit region for the analyses of the η c to be rather wide, considering the relatively large natural width and the relatively large size of the continuum components. As a result, we have a larger second-order polynomial term contribution from the continuum components, which correlates significantly with the size of the charmonium component and could bias the signal yield. We tentatively narrowed the fit region to 2.85 -3.15 GeV/c 2 and took the variations in the obtained yields as a systematic error. They amount to 7%, 13% and 24% for 4π, 2K2π and 4K channels, respectively. We neglect this error source for the other charmonia, where the quadratic components in the continuum are small.
In the fits to the two-meson invariant-mass distributions used to obtain the yields of the two-body decays, we tried a second-order polynomial constrained to vanish at the two-body threshold for each continuum component. We assign the difference between the yields from this method and the standard fit (to a linear function) to the systematic error. Thus, the systematic errors are 20% for η c → f 2 f ′ 2 and less than 3% for the other two-body channels where we have significant signals.
In the case of the three-body modes, we tested an alternative functional form for the continuum component and used the shape of three-body phase space calculated by the MC generator and the detector simulator. We take the difference in the signal yields from this method as the systematic error, that is 4% and 7% for χ c2 → ρ 0 π + π − and χ c0 → K * 0 K ∓ π ± , respectively.

C. Study of non-exclusive backgrounds
Backgrounds can come from general multi-hadron production in two-photon collisions and qq production in single-photon annihilation. Yields from the cc process in the latter category are estimated from the MC and found to be negligibly small (∼ 1%) compared to the non-charmonium backgrounds from two-photon processes. However, the charmonium production rates from these processes cannot be reliably estimated.
A more reliable signature is the p t -balance nature in the signal events. We expect that non-2γ background processes cannot produce a p t -balanced distribution, except for double-ISR processes where exclusive production of C = + charmonia is prohibited. Even in the decays ψ(2S) → γχ cJ where a ψ(2S) is produced in an ISR or double ISR process, the p t distribution of χ cJ will not peak inside the selected region, |Σp * t | < 0.1 GeV/c. We estimate the yield of this process to be about 1% of the χ c0 and χ c2 signals, based on the observed yields of ψ(2S) decaying to π + π − J/ψ.
We evaluate the non-exclusive backgrounds by examining the observed |Σp * t | distributions. Figure 11 shows the |Σp * t | distributions of the charmonium signal component after the sideband subtractions. The distributions are compared with the signal-MC expectations normalized to the sum of the two leftmost bins where the background contamination is expected to be very small. We expect that the background component has a linear shape in |Σp * The MC is normalized to the sum of the two leftmost bins.

D. Total systematic errors
We consider other sources of systematic errors. The statistics of the signal MC events gives errors of 2-4% depending on the process. The calculation of the luminosity function has a systematic uncertainty of 5%. We take the differences among the efficiency estimates for different spin/angular-momentum assumptions for the χ c2 decays to K * 0K * 0 and φφ as systematic errors. The nominal S-wave hypothesis gives the most conservative upper limit for the χ c2 → ρ 0 ρ 0 case. The total systematic errors are listed in Table III. The size of the systematic error relative to the central value is process dependent and ranges from 11% to 35% in the processes where finite G results are presented. The upper limits in the table are shifted upwards by the 1σ systematic errors, in order to include the effect of systematic uncertainties.

VI. COMPARISON WITH PREVIOUS MEASUREMENTS
Some previous measurements give G results for the processes measured here [4,7,8].
We make direct comparisons with these measurements, which are listed in Table III. We quote the world average value from Ref. [19] when two or more previous measurements are available. We can make such direct comparison only for a limited number of processes, since our results include many processes measured for the first time.
It is possible to make some indirect comparisons by converting the G results measured for a different decay mode from the present measured mode and multiplying by the ratio of the branching fractions; where R → A is a "normalization" process for which previous measurement(s) are available. We adopt η c → KKπ, χ c0 → π + π − /K + K − (the average of the two modes) and χ c2 → γJ/ψ as the normalization processes. Note that the two-photon decay widths for these charmonium states are so far measured only in a few decay modes; there are not enough consistency checks among the different decay modes available so far. We use the PDG's fit when the average is unavailable [19], and treat the errors shown there as independent from each other. We find large differences between the present and previous measurements for η c in the direct comparisons, although they are not inconsistent considering the large errors assigned to the previous measurements. However, in the indirect measurements, for many η c decay modes our G values are consistently smaller by factors of two to four compared to the previous determinations. We do not observe any significant signal in the ρ 0 ρ 0 mode in contrast with previous results [23,24]. We conclude that these discrepancies are not likely to be due to statistical fluctuations, and there could be systematic deviations in previous measurements of the relevant partial decay widths or branching fractions.
In contrast, all the present χ c0 and χ c2 results are in agreement with previous values, and confirm the values of Γ γγ (χ c0 ) and Γ γγ (χ c2 ) results derived in previous experiments [19].

VII. UPPER LIMITS FOR η c (2S) PRODUCTION
We search for the η c (2S) in the invariant-mass distributions of the four-meson final states. Since the π + π − π + π − final-state sample has a prominent ψ(2S) peak just above the η c (2S) mass region, we veto events when the invariant mass of any π + π − combination falls within 0.1 GeV/c 2 of the nominal J/ψ mass. (The two tracks from the J/ψ decay are misidentified leptons. ) We obtain upper limits on η c (2S) yields from the four-meson invariant-mass distributions (Fig. 12). The mass and width of the η c (2S) are not very precisely determined so far. We consider a wide range for the mass and width, 3.62 GeV/c 2 < M(η c (2S)) < 3.67 GeV/c 2 and 10 MeV< Γ(η c (2S)) < 40 MeV, according to previous measurements [15,16,19,25,26]. Fits similar to those made in Sect. III.A are applied for the η c (2S) using a Breit-Wigner function with the mass resolution fixed to 9 MeV and a second-order polynomial for the background component, as well as a Gaussian function for the χ c2 peak. The upper limits for the yields are 340, 164 and 55 events at 90% C.L., for the 4π, 2K2π and 4K final-states, respectively. The curves corresponding to these upper limits are shown by the solid curves in Fig. 12. We calculate the upper limits using χ 2 values from fits. We assume various non-negative values of the η c (2S) yield, and, for each case, we obtain the χ 2 from the best fit while floating all the other parameters. We define the 90% C.L. upper limit of the yield as one whose χ 2 is larger by (1.64) 2 than the minimum χ 2 derived in the different assumptions of the yield. An application of the Feldman-Cousins method (Table X in Ref. [20]) to the yields obtained from the fits gives the values very close to the yield upper limits, 338, 163 and 55 events for the 4π, 2K2π and 4K modes, respectively.
We calculate the efficiency assuming a phase space distribution. The upper limit for G ≡ Γ γγ B for each decay mode of the η c (2S) is shown in Table III. We find that the ra- The results in this paper for G ≡ Γ γγ B and comparisons with previous measurements (shown as "direct" and "indirect" below, see text) [7,8,19]. The results with inequalities correspond to 90%-C.L. upper limits, where we include the systematic errors by shifting the upper limits upwards by 1σ(sys). The "isospin ·BF " factors are used to obtain the G results from the experimental measurements.

VIII. CONCLUSION
The production of the η c (1S), χ c0 and χ c2 charmonium states in two-photon collisions has been observed in all of the four-meson final states, π + π − π + π − , K + K − π + π − and K + K − K + K − . We used data samples with one or two orders of magnitude larger statistics than previous measurements. No clear signature for the η c (2S) is found in any of decay processes, and we obtain the upper limits for the products of its two-photon decay width and the branching fractions. We have studied resonant substructures in these four-meson final states.
For the first time χ cJ signals produced in two-photon collisions are observed in the K + K − π + π − or K + K − K + K − final states.
We also find a new decay mode, η c → f 2 (1270)f ′ 2 (1525). We have obtained products of the two-photon decay width and branching fractions for various decays of charmonium states. The present results for η c are systematically smaller than the derived values from the world averages of previous measurements.