Inclusive charged and neutral particle multiplicity distributions in χ cJ and J = ψ decays

Using a sample of 106 million ψ ð Þ decays, ψ ð 3686 Þ → γχ cJ ð J ¼ 0 ; 1 ; 2 Þ and ψ ð 3686 Þ → γχ cJ ; χ cJ → γ J= ψ ð J ¼ 1 ; 2 Þ events are utilized to study inclusive χ cJ → anything, χ cJ → hadrons, and J= ψ → anything distributions, including distributions of the number of charged tracks, electromagnetic calorimeter showers, and π 0 s, and to compare them with distributions obtained from the BESIII Monte Carlo simulation. Information from each Monte Carlo simulated decay event is used to construct matrices connecting the detected distributions to the input predetection “ produced ” distributions. Assuming these matrices also apply to data, they are used to predict the analogous produced distributions of the decay events. Using these, the charged particle multiplicities are compared with results from MARK I. Further, comparison of the distributions of the number of photons in data with those in Monte Carlo simulation indicates that G-parity conservation should be taken into consideration in the simulation.


I. INTRODUCTION
The multiplicity distributions of charged hadrons, which can be characterized by their means and dispersions, are an important observable in high energy collisions and an input to models of multihadron production. Charged particle means from below 2 GeV to LEP energies have been fit as a function of energy with a variety of models in Ref. [1], and a review of theoretical understanding can be found in Ref. [2].
The study of χ cJ ðJ ¼ 0; 1; 2Þ decays is important since they are expected to be an important source of glueballs, and future studies require both more data and better simulation of generic χ cJ decays. Also since χ cJ decays make up approximately 30% of ψð3686Þ decays, a better understanding of χ cJ decays improves that of ψð3686Þ decays.
The branching fractions of ψð3686Þ → γχ cJ and χ cJ → γJ=ψ were measured previously by BESIII using a sample of 106 million ψð3686Þ decays [3]. The accuracy of these measurements depends critically on the ability of the Monte Carlo (MC) simulation to model data well. Since a large fraction of χ cJ hadronic decay modes are still unmeasured [4], it is particularly important to verify the modeling of their inclusive decays, where we rely heavily on the LUNDCHARM model [5] to simulate these events.
In this paper, which is based on the analysis performed in Ref. [3], we report on the "detected" distributions: the efficiency-corrected charged particle multiplicity distributions, as well as the efficiency-corrected distributions of the number of electromagnetic calorimeter showers and π 0 s for χ cJ and J=ψ decays. Our detected distributions are compared with MC simulation, and the results can be used to improve the LUNDCHARM model simulation, in particular for χ cJ hadronic decays.
Information from each MC simulation decay event is used to construct matrices connecting the detected charged particle and photon multiplicity distributions to the input predetection distributions. Assuming the matrices also apply to data, they are used to predict the analogous "produced" distributions of the decay events. Produced charged particles and photons correspond to those coming directly from the χ cJ or J=ψ decays or the decays of their daughter particles. The means of the charged particle a Also at Ankara University, 06100 Tandogan Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP 3 . multiplicity distribution are compared with those of MARK I, which measured the mean charged particle multiplicity for e þ e − → hadrons as a function of center-of-mass energy from 2.6 to 7.8 GeV [6].
In Ref. [3], an electromagnetic calorimeter (EMC) shower (EMCSH) was labeled a "photon", but as described in Sec. IVA, showers include hadronic interactions in the EMC crystals and electronic noise, so here we will explicitly refer to them as EMCSHs. The comparison of data and inclusive ψð3686Þ MC simulation showed good agreement for charged track distributions and most EMCSH energy (E sh ) distributions; however, there was some difference in the distribution of the number of π 0 s [3]. Here, we explore the agreement for χ cJ → anything and χ cJ → hadrons via ψð3686Þ → γχ cJ and J=ψ → anything via χ cJ → γJ=ψ. Recently BESIII observed electromagnetic Dalitz decays χ cJ → l þ l − J=ψðl ¼ e or μÞ [7], so our χ cJ → hadron distributions also include χ cJ → l þ l − J=ψ. However, the branching fractions for these decays are very small, on the order of 10 −4 , which are negligible compared with those of χ cJ → hadrons. Below we will continue to refer to these distributions as χ cJ → hadrons. "Hadrons" is used very loosely and includes all processes except χ cJ → γJ=ψ, such as other χ cJ radiative decays and χ cJ → γγ.
This analysis is based on the 106 million ψð3686Þ event sample gathered in 2009, the corresponding continuum sample with integrated luminosity of 44 pb −1 at ffiffi ffi s p ¼ 3.65 GeV [8] and a 106 million ψð3686Þ event inclusive MC sample.
The paper is organized as follows: In Sec. II, the LUNDCHARM model is described. In Secs. III-V, the distributions of the number of detected charged tracks, EMCSHs, and π 0 s, respectively, are determined and compared with MC simulation. Section VI presents the produced distributions. Section VII discusses systematic uncertainties, while Sec. VIII provides a summary. Additional EMCSH and π 0 tables are included in an appendix.

II. LUNDCHARM MODEL
The LUNDCHARM model is an event generator to produce events for charmonium decaying inclusively to anything [5]. This model, which was inspired by QCD theory, was developed at BESII and migrated to the BESIII experiment. In this model, J=ψ or ψð3686Þ decaying into light hadrons is described as cc quark annihilation into one photon, three gluons or one photon plus two gluons, followed by the photon and gluons transforming into light quarks and further materializing into final light hadron states. To leading order accuracy, the cc quark annihilations are modeled by perturbative QCD [9], while the hadronization of light quark fragmentation is described with the Lund model [10] using a set of parameters to describe the baryon/meson ratio, strangeness and fη; η 0 g suppression, and the distribution of orbital angular momentum, etc.
The LUNDCHARM model is used to generate the unmeasured charmonium decays, while the established decays are exclusively generated with their appropriate BesEvtGen models [11] using branching fractions from the Particle Data Group [4]. The fraction of unmeasured decays for each charmonium state is given in Table I [4]. Since the fractions are quite large for χ cJ decays, the LUNDCHARM model is very important for the simulation of these decays. The parameters of the LUNDCHARM model are optimized using 20 million J=ψ decays accumulated at the BESIII   [12]. Figure 1 shows the comparison between data and MC simulation of the multiplicity of detected charged tracks for J=ψ and ψð3686Þ decays. More comparisons of data and MC simulation for J=ψ decays can be found in Ref. [12] and for ψð3686Þ decays in Refs. [3,12].

A. Method
The basic approach is the same as in Ref. [3]. Charged tracks must be in the active region of the drift chamber and have their points of closest approach consistent with the run-by-run interaction point. Neutral tracks must be in the active regions of the barrel EMC or end cap EMC, satisfy minimum and maximum energy requirements and a time requirement. The basic ψð3686Þ event selection requires at least one charged track (except for the study of the events with no charged tracks, where this requirement is dropped), at least one neutral track, and a minimum event energy. A background filter removes non-ψð3686Þ events, and events consistent with being a ψð3686Þ → ππJ=ψ decay are removed [3]. Following this, the E sh distribution is constructed for the remaining events, where the EMCSH must be in the barrel EMC, not originate from a charged track (δ > 14°, where δ is the angle between the shower and the nearest charged track), and not be a photon from a π 0 decay. Fitting the peaks in the E sh distribution due to ψð3686Þ → γχ cJ and χ cJ → γJ=ψ, as shown in Fig. 2, allows the determination of the number of the inclusive decays and the final branching fractions. Please refer to Ref. [3] for many important details.
To determine the distribution of the number of charged tracks, N ch , ten E sh distributions are constructed for N ch ranging from 0 to 9. These distributions are then fitted to determine the numbers of χ cJ → anything and χ cJ → γJ=ψ; J=ψ → anything events, and these numbers determine the N ch distributions for χ cJ → anything and J=ψ → anything.
In Ref. [3], simultaneous fitting of inclusive and exclusive E sh distributions was performed, but this is not done here, except for the N ch ¼ 0 case, because there are no exclusive E sh distributions versus N ch to be used in such a fit. Another change is that events with N ch ¼ 0 have additional requirements in order to reduce the background in the E sh distributions.
B. N ch = 0 event selection and fit of E sh distributions Events with N ch ¼ 0 were selected in our previous analysis only to determine the systematic uncertainty associated with the N ch > 0 requirement. The photon time requirement was removed since without charged tracks, the event time is not well determined. Although other selection requirements were tightened, the events still had much background [3].
For the current analysis, events with jðP x Þ neu j > 1.0 GeV=c and jðP y Þ neu j > 1.0 GeV=c are removed, since these regions contain much background according to a MC simulation. ðP x Þ neu and ðP y Þ neu are the sum of the momenta of all neutrals in the x and y directions, respectively, where x and y are orthogonal axes perpendicular to the axis of the detector. The E sh distribution with the additional requirements is much cleaner and easily fitted, as shown in Fig. 2. A simultaneous fit with inclusive and exclusive events was used for the previous N ch ¼ 0 systematic uncertainty study since the signal to background ratio was so low, and the same fitting method is used here, as shown in Fig. 2. The χ 2 =ndf for the fit to data is 1.3, where ndf is the number of degrees of freedom.
C. N ch > 0 selection and fitting Figure 3 shows the E sh distributions for all N ch and for individual values of N ch > 0 for data. E sh distributions for different values of N ch for MC simulation and continuum background are constructed similarly. FIG. 2. Simultaneous fits to the E sh distributions of data for N ch ¼ 0. Top set: Inclusive N ch ¼ 0 distribution fit and corresponding pulls. Bottom set: exclusive distribution fit and pull distribution. The five peaks from left to right in the top figure correspond to ψð3686Þ → γχ c2 , γχ c1 , γχ c0 , χ c1 → γJ=ψ and the small χ c2 → γJ=ψ contribution (see arrow). The exclusive modes include ψð3668Þ → γχ cJ ; χ cJ → 2 and 4 charged track events, selected with requirements on the invariant mass of the charged tracks and the angle between the direction of the radiative photon and the recoil momentum from the charged tracks. Here the wide χ cJ → γJ=ψ shapes are described by the inclusive MC shapes, while the narrow ψð3686Þ → γχ cJ shapes are inclusive MC shapes convolved with bifurcated Gaussians. The smooth curves in the two plots are the fit results. The dash-dotted and dotted curves in the top plot are the background distribution from the inclusive ψð3686Þ MC with radiative photons removed and the total background, respectively, where the total is the sum of the MC background and a second order polynomial.
Signal shapes and background shapes used in the fit depend on the value of N ch . In fitting the distributions for N ch > 7, because of the small sample sizes, the signal shapes and background shapes for N ch ¼ 7 are used. The fit result of data for N ch ¼ 5 is shown in Fig. 4, and the χ 2 =ndf is 1.4. Fit results for other values of N ch result in similar χ 2 =ndf values.
The MC simulated sample is fitted as a function of N ch in a similar fashion, but ψð3686Þ → γχ cJ MC signal shapes are fitted without convolution. As described in Ref. [3], the MC events are weighted by wt π 0 × wt trans , where wt π 0 accounts for the difference between data and MC simulation on the number of π 0 s and wt trans accounts for the E 3 γ energy dependence of the radiative photon in the electric dipole transitions for ψð3686Þ → γχ cJ and χ cJ → γJ=ψ.

D. Results
The MC simulated sample is analyzed by counting the number of events versus N ch before applying any selection criteria. The efficiency is then the number of events passing all selection criteria divided by the number of events without imposing any selection versus N ch . Note that N ch here is the "detected" number of charged tracks.
Using the number of detected data events, D, and the MC determined efficiencies, ϵ, which are dependent on N ch , we determine the distribution of the efficiency-corrected number of events in data for χ cJ → anything and χ cJ → γJ=ψ; J=ψ → anything. Results are listed in Table II for χ cJ → anything and Table III for χ c1=2 → γJ=ψ; J=ψ → anything.
For comparison, MC simulation numbers, N MC , are also listed in the tables. N MC corresponds to the N ch distribution before imposing selection requirements. Since the branching fractions of MC simulation are not the same as the measured branching fractions of Ref. [3], the MC numbers are scaled by B BESIII =B MC , where B BESIII and B MC are the  BESIII branching fractions [3] and those used by the MC, respectively, and the N MC in Tables II and III are the scaled MC numbers.
The efficiency corrected N ch distributions for χ cJ → anything contain the χ cJ → γJ=ψ; J=ψ → anything events, as well as the χ cJ → hadrons events. A more interesting comparison between data and the simulated MC sample is with the N ch distributions for χ cJ → hadrons directly. These are obtained by subtracting N ch distributions for χ cJ → γJ=ψ; J=ψ → anything from those of χ cJ → anything. Since we do not have the distribution from data for χ c0 → γJ=ψ; J=ψ → anything, we use the MC distribution TABLE II. Detected data events, D, efficiencies, ϵ, efficiency corrected events, N, and number of scaled simulated events N MC for χ cJ → anything.  Here and below, the first uncertainties are the uncertainties from the fits to the inclusive E sh distributions and the second ones are systematic, described in Sec. VII.  Comparison of fraction of events in percent with N ch for data and the scaled MC simulated sample for χ c1;2 → γJ=ψ; J=ψ → anything. These two sets of measurements describe the same distribution.   Fig. 8 below, the uncertainties shown for MC are the uncertainties from the fits to the inclusive E sh distributions, and the uncertainties for data are those combined in quadrature with the systematic uncertainties, described in Sec. VII. for this process. The branching fraction is small, 1.4%, so the change for χ c0 → anything is small.
The N ch fractions, F, where the fraction is the number of efficiency corrected events with N ch ¼ j (j takes on values from 0 to 9) divided by the sum of all N ch events, are determined and are listed in Table IV for χ cJ → hadrons and Table V for χ c1=2 → γJ=ψ; J=ψ → anything. For comparison, MC simulation numbers, F MC , are also listed in the tables. F MC is calculated in an analogous way as was F using the scaled MC simulation numbers. In Figs. 5(a), 5(c), and 5(e) comparisons of the N ch fractions between data and scaled MC simulated sample are shown, while Figs. 5(b), 5(d), and 5(f) are the corresponding plots in logarithmic scale. Figure 5 shows good agreement between the three χ cJ → anything decay distributions. Data are above MC simulation for N ch ¼ 0 and N ch > 5 and below for N ch ¼ 3 for these distributions. The agreement between data and MC simulation is good for J=ψ → anything (χ c1 and χ c2 → γJ=ψ). Better agreement is expected for those distributions, since MC tuning was performed on the J=ψ → anything events.

A. MC study of EMC energy deposits
The situation for neutral showers is more complicated than for charged tracks. Energy deposits in the EMC from ψð3686Þ → γχ cJ and χ cJ → γJ=ψ events are caused by their radiative photons, photons from the decays of π 0 s from χ cJ and J=ψ hadronic decays and their daughter particles, bremsstrahlung from charged tracks, as well as interactions of hadrons in the EMC crystals and noise. The inclusive MC needs to model all these sources. We are interested in the number of photons, N γ , from the hadronic decays of χ cJ and J=ψ. We can use the MC simulation to determine what fraction of the EMCSHs are due to radiative photons and the photons from the primary and secondary decays. We signify the number of EMCSHs by N sh .
The MC "truth" information tags the radiative photons in the generator model and photons from the generator final particle decays in GEANT4 [13], e.g., π þ → μ þ ν μ γ, as well as final-state radiation photons. The MC truth does not tag the photons produced from the scattering and/or ionization of generator final state particles with the detector materials, simulated by GEANT4. The angles of tagged photons can be compared with the angles of EMCSHs to identify the fraction of showers that are caused by these photons. Figure 6 shows for a small subsample of ψð3686Þ → γχ cJ events the angle D θ , which is the minimum of the difference in angle between an EMCSH and all the MC tagged photons. There is a sharp peak at small D θ corresponding to good shower matches between the MC predictions and the EMCSHs. We define showers with D θ < 0.1 radians as a good shower match. The efficiency of matching photons in the correct angular range (jcos θj < 0.8) and energy range (0.25 GeV < E sh < 2 GeV) is 91.2%.
The fraction of good matches varies from 60% at the lowest energy to 89% at the highest. Figures 7(a) and 7(b) show the number distributions of all and good showers, respectively. In the following, we will compare the N sh distributions of data and MC simulation.

B. N sh distribution
The analysis for the distribution of N sh is similar to that for N ch . N sh is the number of showers satisfying requirements on the energy, polar angle, and time, but no requirement on the angle between the shower and the closest charged track in the event. Here 15 energy distributions are constructed for N sh ranging from 1 to ≥15, where N sh ¼ 1 is because at least one radiative photon must be detected. For more direct comparison of data with MC simulation, MC events are weighted only by wt trans .
As above, using the number of detected data events, D, and the MC determined efficiencies, ϵ, versus N sh , we determine the efficiency corrected N distributions of data for χ cJ → anything and χ cJ → γJ=ψ; J=ψ → anything. Results are listed in the Appendix in Table XX for χ cJ → anything and Table XXI for χ c1=2 → γJ=ψ; J=ψ → anything. The N sh fractions, F, are also determined and are listed in Table VI for χ cJ → hadrons and Table VII Fig. 8 are similar for the three χ cJ decays, and data are above the MC simulation for N sh ¼ 1 and N sh > 7 and below for N sh ¼ 3 and 6. For J=ψ → anything (χ c1 and χ c2 → γJ=ψ), there is only minor disagreement between data and MC simulation for the N sh distributions.
V. Multiplicity distribution of the number of π 0 s An even more complicated case is the distribution of the number of π 0 s, N π 0 . Here, as for the N ch ¼ 0 case, the N π 0 distribution is considered in more detail. The γγ invariant mass, M γγ , distribution of the π 0 candidates is shown in Fig. 9, where there are a large number of γγ miscombinations in the plot. A somewhat better estimate of N π 0 is made with the restrictive requirement 0.120 < M γγ < 0.145 GeV=c 2 , which was the requirement used when vetoing EMCSHs that might be part of a π 0 combination from the E sh distribution used in the fitting for the number of ψð3686Þ → γχ cJ and χ cJ → γJ=ψ events [3]. However, even with this requirement there are still many γγ miscombinations. TABLE VII. Comparison of the fraction of events in percent with N sh between data and the scaled MC simulated sample for χ c1;2 → γJ=ψ → γ anything. These two sets of measurements measure the same distribution and are in agreement within uncertainties.  VI. Comparison of fraction of events in percent with N sh between data and the scaled MC simulated sample for ψð3686Þ → γχ cJ → γ hadrons. To determine the fraction, R, of the π 0 candidates that are valid π 0 s, we fit the M γγ distributions for 0.120 < M γγ < 0.145 GeV=c 2 for each N π 0 for both data and the MC simulated sample to a signal shape and first order Chebychev polynomial background. The basic signal shape was determined using the MC truth information to identify correct γγ combinations in simulated data. For data, the basic signal shape is convolved with a bifurcated Gaussian function to account for the difference in resolution between data and the MC simulated sample. R is the fraction of signal events in the region 0.120 < M γγ < 0.145 GeV=c 2 . The values of R versus N π 0 are listed in Table VIII.
Note that N π 0 may not fully determine the number of valid π 0 s. For instance, N π 0 ¼ 3 may include the cases of three valid π 0 s, two valid π 0 s and one miscombination, one valid π 0 and two miscombinations, and three miscombinations.
The analysis for the detected N π 0 distributions is similar to those for N ch and N sh . Here ten E sh distributions of data are constructed for N π 0 ranging from 0 to ≥9. For more direct comparison of data with MC simulation, MC events are weighted only by wt trans .
Using the number of detected data events, D, the MC determined efficiencies, ϵ, and RðdataÞ versus N π 0 , we determine the efficiency corrected N distributions of data for χ cJ → anything and χ cJ → γJ=ψ; J=ψ → anything, where N ¼ R · D=ϵ, which gives a better representation of the N π 0 distribution. Results are listed in the Appendix in Table XXII for χ cJ → anything and Table XXIII   In Figs. 10(a), 10(c), and 10(e) comparisons of the N π 0 fractions between data and scaled MC simulated samples are shown, and Figs. 10(b), 10(d), and 10(f) provide logarithmic versions. For χ cJ → hadrons, the N π 0 distribution, data are above MC simulation for N π 0 > 2. For J=ψ → anything (χ c1 and χ c2 → γJ=ψ), data are above MC simulation for N π 0 > 5, but the uncertainties are bigger for these decays.

VI. PRODUCED DISTRIBUTIONS
So far, we have only dealt with the distributions of the efficiency-corrected number of detected charged tracks, EMCSHs, or pions. These depend on the geometry and performance of the BESIII detector. Of more interest are the actual physics distributions in the decays of the χ cJ and J=ψ.
To determine these distributions from data, we construct detection matrices using the χ cJ → hadrons and  Comparison of fraction of events in percent with N π 0 for data and scaled MC simulated sample for χ cJ → hadrons. Both F and F MC are based on numbers of events multiplied by R, the fraction of valid π 0 s. The first uncertainties are the uncertainties from the fits to the inclusive E sh distributions and R, and the second are the systematic uncertainties, described in Sec. VII.   9. The M γγ distribution of π 0 candidates reconstructed without the tight π 0 mass selection requirement. Data are represented by dots, and the MC sample by the red and shaded histograms for the MC events weighted by wt π 0 and unweighted events, respectively. χ cJ → γJ=ψ; J=ψ → anything events in the inclusive ψð3686Þ MC events. The matrix (M) times the produced vector (P) determines the detected vector (D), where (P i ) is the number of events with i charged tracks, photons, or π 0 s, etc. 0 The elements of M are determined using the MC "truth" information by tallying the detected versus the produced track information for each event. The detection matrix M is then assumed to apply to data, as well as to MC simulation. Detected histograms are constructed corresponding to each element in the P vector using the matrix equation (1). These are used to give a set of probability density functions (PDFs) with which to perform a χ 2 fit of the detected distributions of data to determine the values for P 0 ; …; P N .
A. P N P ch distributions The results of the fits to the detected charged track distributions of data to determine the produced charged track distributions P N P ch are shown in Fig. 11 for χ cJ → hadrons. Here N P ch refers to the number of produced through 15 below, the produced uncertainties are the uncertainties from the fits for P N P ch combined in quadrature with the systematic errors, described in Sec. VII. The data uncertainties and the fitted fraction uncertainties in (d)-(f) are the uncertainties from the fits to the inclusive E sh distributions and the uncertainties from the fits for P N P ch , respectively. Also shown in these plots are the PDFs used in the fits. The distribution is fitted over bins N ch ¼ 0-8. In Figs. 11(d)-11(f) and 12(c)-12(d), nine bins of detected data are fitted with six parameters (P 0 through P 10 ) and with P 12 fixed to the MC values. Fractions F P of χ cJ → hadrons and χ c1=2 → J=ψ; J=ψ → anything are listed in Tables XI and XII, respectively. The χ 2 =ndf values for the five cases are 65, 52, 85, 18, and 28. Alternative fits with P 12 free give the same results as shown in Tables XI and XII. Comparing the fits and the PDFs in Figs. 11 and 12 suggests that the MC PDFs do not describe data well, which contributes to the large χ 2 =ndf. However, corrections to the PDFs to improve the fits to the detected charged track distributions, as described in Sec. VII, result in small changes to the P N P ch values compared with the systematic uncertainties shown in Table XIX and are neglected.

Mean charged multiplicity and dispersion
We determine values of the mean multiplicity hN P ch i, dispersion D ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi h½N P ch 2 i − hN P ch i 2 p , and hN P ch i=D. Such measurements have been performed for e þ e − → hadrons at LEP [1] and also at lower energies with the MARK I experiment [6]. The results of these measurements from our data are listed in Table XIII. Although we measure J=ψ → anything via χ cJ → γJ=ψ; J=ψ → anything, we can calculate the J=ψ → hadron distribution using the branching fractions of J=ψ → e þ e − and μ þ μ − [4] and assuming that these events populate N P ch ¼ 2 only. The calculated values are also listed in Table XIII. Our values for hN P ch i can be compared with those of MARK I for e þ e − → hadrons [6]. The MARK I values TABLE XI. P N P ch event fractions in percent for data F P χ cJ and MC simulated sample F MC χ cJ for χ cJ → hadrons. In the fit, the value of the fraction for N P ch ¼ 12 is fixed to the MC result. Here and below, the first uncertainties shown are the uncertainties from the fits for P N P ch , and the second are the systematic uncertainties described in Sec. VII.  Fig. 13 along with our values for both J=ψ and χ cJ to hadrons. While our results include statistical and systematic uncertainties, those of MARK I do not include systematic uncertainties, which range from 25% at 2.6 GeV to 15% at 6 GeV and above. Still, the agreement between the results is very good.
B. P N P γ distributions P N P γ distributions are studied in an analogous way. Here the P N P γ distributions correspond to the MC-tagged photons, described in Sec. IV, and the detected distributions are the EMC shower distributions, which include both good and bad shower matches. The results of the fits for the P N P γ distributions are shown in Fig. 14 for χ cJ → hadrons. Shown in Figs. 14(a)-14(c) are the MC fractions and the results from the fits to the detected distributions of data. The radiative photons from ψð3686Þ → γχ cJ are not counted, so the lowest bin is N P γ ¼ 0. Even bins are much larger than odd ones since most photons are from π 0 → γγ decays. Photons from final-state radiation (FSR) and radiative photons from intermediate-state decays are counted and contribute to odd bins. While fit results for bins N P γ ¼ 2, 6, and 10 are smaller than MC, those for N P γ ¼ 0, 4, 8, and 12, which correspond to an even number of π 0 s, are much larger than MC. The detected data fractions as a function of N sh and the fractions determined from the fit results are shown in Figs. 14(d)-14(f).
Results for χ c1=2 → γJ=ψ; J=ψ → anything are shown in Fig. 15. The MC fractions and the results from the fits to the detected distributions of data are shown in Figs. 15(a)-15(b). Since radiative photons from χð3686Þ → γχ cJ and χ cJ → γJ=ψ are not counted, the lowest bin is also N P γ ¼ 0. Here, for bins with N P γ ¼ 2, 6, and 10, which correspond to a preference for an odd number of π 0 s, fit Since at least one EMCSH must be detected, the data fraction for N sh ¼ 0 is empty. results are slightly larger than MC, but uncertainties are large. The detected data fractions versus N sh , and the fit results are shown in Figs. 15(c)-15(d). The G-parity for χ cJ s is positive, suggesting that decays to pions should favor an even number of πs, while G-parity for the J=ψ is negative, implying that decays to pions favor an odd number of πs. These preferences in the distributions of the number of photons are observed above for data, but MC simulation does not reflect this.
In Figs. 14(d)-14(f) and 15(c)-15(d), 14 bins of detected data are being fit with seven parameters (P 0 , P 2 , P 4 , P 6 , P 8 , P 10 , and P 12 ) with in between bins fixed to MC values. The number of degrees of freedom is ndf ¼ 7. Fractions F ndf¼7 of χ cJ → hadrons are listed in Table XIV and χ c1=2 → γJ=ψ; J=ψ → anything are listed in Table XV. The χ 2 =ndf values for the five cases are 17, 7.8, 4.3, 5.7, and 2.9. C. P N P π 0 distributions P N P π 0 distributions are studied in a similar fashion. Here, the situation is complicated because events for a given N π 0  , and MC simulated sample, F MC χ cJ , for χ cJ → hadrons. Odd bins are fixed to MC result values, so only the systematic uncertainties are shown. Here and in Table XV, it is the number of events that is fixed, so the fractions may differ slightly. contain miscombinations as well as real π 0 s. We assume that the MC simulation correctly describes the miscombinations in data and do not multiply by R. Unlike the cases above, the alternate bins are not suppressed, so adjacent PDFs are similar in shape, which results in larger fit uncertainties for the values of P N P π 0 . The low sensitivity that the fit has to P N P π 0 has other consequences. For most of the variations used in Sec. VII to determine N P π 0 systematic uncertainties, the fits of the detected distributions of data fail, with only three successful fits out of a total of nine. See Sec. VII for details. In conclusion, we are not able to determine the systematic uncertainties for the P N P π 0 distributions corresponding to the detected π 0 distributions and therefore the event fractions themselves.

D. Input-output check
The procedures above have been repeated using MC detected distributions as input. The output produced Results from the input-output test. The output (F P ) and input MC (F MC ) fractions of events in percent with N P γ for χ cJ → hadrons (F χ cJ ) and χ c1=2 → γJ=ψ (F J=ψ 1=2 ). The P N P γ values for N P γ odd are fixed to those of the MC in the fitting. for χ cJ → hadrons (F χ cJ ) and χ c1=2 → γJ=ψðF J=ψ 1=2 Þ. The P N P ch values for N P ch ≥ 12 are fixed to those of the MC in the fitting and are not listed. distributions determined by the analyses should then agree closely with the MC truth distributions. We divide the MC data in half and use the first half to construct the detection matrices and use them in fitting the detected distributions of the second half. We compare the fitting results with the MC truth fractions of the second half. For this check, the uncertainties on the detected distributions are taken as the statistical uncertainties on the number of detected events combined in quadrature with the statistical uncertainties on the number of MC events. The output fitted fractions F P and input MC fractions F MC versus N P ch are given in Table XVI, where the agreement is very good. The χ 2 =ndf values for the fits are 1.2, 0.9, 0.8, 0.5, and 0.7.
The output fitted fractions F P and input MC fractions F MC versus N P γ are given in Table XVII. The agreement for these cases is not as good as for the N P ch cases. The χ 2 =ndf values for N P γ are 1.2, 0.5, 0.2, 1.5, and 2.2. In all cases, the differences between input and output are small compared to the systematic uncertainties detailed in Table XIX and are neglected.

VII. SYSTEMATIC UNCERTAINTIES
Extensive studies of systematic uncertainties were carried out in Ref. [3]. For the ψð3686Þ → γχ cJ branching fraction, they are under 4% with the largest contribution coming from fitting the E sh distribution. Many of the uncertainties do not apply here. For the distribution of the number of charged tracks, the uncertainty from N ch > 0 does not apply, since we include events with no charged tracks. The requirement N sh < 17 essentially includes all events, as does E vis > 0.22E cm , which has a small systematic uncertainty. The uncertainty from N ψð3686Þ does not apply since we calculate event fractions. Those for MC signal shape, multipole correction, and j cos θj < 0.80 affect the selection of the radiative photon candidate and the overall number of events, but should not affect the various distributions.
Systematic uncertainties are determined here for detected event fractions in Secs. III-V and for produced event fractions in Secs. VI A-VI B using samples selected with alternate selection criteria and with modified fitting procedures. Systematic uncertainties are the differences from the standard procedure added in quadrature.
For all distributions, the fitting uncertainties are determined by changing the background polynomial, changing the range, and fixing small signals. Background polynomials are changed from second order to first, and the fit ranges are changed from 0.08-0.5 GeV to 0.08-0.35 GeV and 0.2-0.54 GeV.
For the detected charged track event fraction systematic uncertainties in Sec. III and the P N P ch fraction uncertainties in Sec. VI A, (a) the fitting uncertainties are considered, and in addition, uncertainties from (b) the ψð3686Þ background veto, (c) the π þ π − J=ψ veto, (d) the δ > 14°requirement, and (e) the continuum energy difference are determined. The uncertainties from (b)-(d) are determined by removing those requirements and comparing with the analyses making them. The uncertainty from (e) is determined by scaling the EMCSH energies of the continuum events by the ratio of the center-of-mass energies of ψð3686Þ data and the continuum data.
For the detected photon event fraction systematic uncertainties in Sec. IV, the P N P γ event fraction systematic uncertainties in Sec. VI B, and the detected pion event fraction uncertainties in Sec. V, the fitting errors are considered. In addition, uncertainties from the ψð3686Þ background veto, the π 0 π 0 J=ψ veto, the δ > 14°requirement, and continuum energy difference are determined.
An important question is what are the systematic uncertainties associated with the determination of the produced distributions by fitting the detected distributions in Sec. VI. We study this in two ways. In the first, we modify the PDFs used in determining the P N P ch values. In Sec. VI A the fits have large χ 2 =ndf, and the PDFs in Figs. 11(d)-11(f) and 12(c)-12(d) do not appear to describe the detected distributions of data well. Data are above the fit in the highest bin and below in the preceding bin for the N P ch ¼ 4 and N P ch ¼ 6 PDFs in Figs. 11(d)-11(f) and 12(c)-12(d). The PDFs determined from the inclusive MC seem to be too broad.
The PDFs have been modified to see the effect on the χ 2 =ndfs and further to determine the differences in the P N P ch results. We assume that part of the PDFs may be described approximately by a binomial distribution. For instance, we assume that the PDF for N P ch ¼ 4 for N ch ¼ 0 through N ch ¼ 4 can be described by a binomial distribution in terms of an efficiency ϵ 4 , which includes the geometric, tracking, and vertexing efficiencies, and the fraction of N ch ¼ 4 in N ch ¼ 0 through N ch ¼ 4 is given according to a binomial distribution by ϵ 4 4 . The PDFs of the MC being too wide is due to the efficiencies being too small. We estimate the corrected efficiency approximately by comparing the data fractions and the fitted fractions for the N ch ¼ 4 bins in Figs. 11(d)-11(f) and 12(c)-12(d), , where D is the data fraction bin content of N ch ¼ 4 and F is the fitted fraction bin content. For Fig. 11(d), ϵ 4 ¼ 0.8630 and ϵ 4 corr ¼ 0.8685; the difference is only 0.64%. We then use the ratio of the binomial distributions in terms of the two efficiencies in each bin to correct the MC PDFs and use the corrected PDFs to fit the detected distributions. The part of the PDF for N ch > 4 is left unchanged. We do analogous calculations for N P ch ¼ 2 and N P ch ¼ 6. The corrected PDFs fit the detected distributions much better now than they did in Figs. 11 and 12, and the χ 2 =ndfs become 15. 5, 12.3, 17.7, 5.4, and 9.0, which are much reduced compared to those in Sec. VI A (65, 52, 85, 18, and 28). The differences with the uncorrected fractions are small compared with the systematic uncertainties shown in Table XIX and are neglected.
In the second, more quantitative study, we modify the selection criteria for both charged tracks and EMCSHs by requiring that they satisfy j cos θj < 0.8, corresponding to the barrel shower counter. The detected distributions are greatly altered by such a requirement. This is easy to understand: the probability of removing a charged track goes way up when there are a large number of charged tracks. The detected distributions are pushed to lower values, while the produced truth distributions of MC are not affected. For example, the means of the detected N ch distributions are 3.82, 3.58, and 3.68 (3.20, 3.14, and 3.25) for the standard (j cos θj < 0.8) selection for χ c0 , χ c1 , and χ c2 → hadrons, while the means of the produced distributions are 4.27, 4.44, and 4.46 (4.28, 4.44, and 4.48), respectively, for the standard (j cos θj < 0.8) selection. The differences with the standard selection for both P N P ch and P N P γ values are taken as the systematic uncertainties associated with the determination of the produced distributions by fitting the detected distributions.
The detected event fraction uncertainties and the P N P ch=γ event fraction uncertainties are listed in Tables XVIII and XIX, respectively. The uncertainties are the individual uncertainties for all cases added in quadrature.

VIII. SUMMARY
The study of χ cJ decays is important since they are expected to be a source for glueballs, and their simulation is TABLE XVIII. N ch , N sh and N π 0 detected event fraction systematic uncertainties in percent. In the table, χ cJ represents χ cJ → hadrons and J=ψ 1=2 represents χ c1=2 → γJ=ψ; J=ψ → anything.
Using 106 million ψð3686Þ decays, we study χ cJ → anything, χ cJ → hadrons, and J=ψ → anything distributions. Distributions of event fractions for data are compared with the MC simulation versus the number of detected charged tracks, EMCSHs and π 0 s in Figs. 5, 8 and 10, respectively. For all comparisons, the agreement is reasonable. However, there are differences.
To start with χ cJ → anything, for the N ch distributions, data are above MC simulation for N ch ¼ 0 and N ch > 5 and below for N ch ¼ 3 and 6. For the N sh distributions, data are above MC simulation for N sh ¼ 1 and N sh > 7 and below for N sh ¼ 3, and for the N π 0 distributions, data are above MC simulation for N π 0 > 2.
For J=ψ → anything (χ c1 and χ c2 → γJ=ψ), the agreement between data and MC simulation is good for the N ch distributions. There is some disagreement for the N sh distributions, and for the N π 0 distributions data are above MC simulation for N π 0 > 5, but the uncertainties are bigger. Better agreement is expected for J=ψ → anything distributions, since MC tuning was performed on the J=ψ → anything events.
For χ cJ → hadron charged track distributions, fit results shown in Fig. 11 for P N P ch are below the MC fractions for N P ch ¼ 4 and above for N P ch ¼ 0, 8, and 10. P N P ch results for χ c1=2 → γJ=ψ; J=ψ → anything charged track distributions are shown in Fig. 12. The distributions are similar, and the fit fractions are in reasonable agreement with the MC fractions. The means of the above N P ch distributions in Figs. 11 and 12 are determined and plotted along with results from MARK I for e þ e − → hadrons in the same energy range in Fig. 13. The charmonium decays to hadrons and e þ e − → hadrons results are consistent.
The results for the P N P γ distributions are shown in Figs. 14(a)-14(c) for χ cJ → hadrons. The content of even bins are much larger than those of odd ones since most photons are from the decay of π 0 s. While fit results for bins N P γ ¼ 2, 6, and 10 are smaller than MC, those for N P γ ¼ 0, 4, 8, and 12, which correspond to an even number of π 0 s, are much larger than MC. Results for χ c1=2 → γJ=ψ; J=ψ → anything for photons are shown in Fig. 15. Here, bins with N P γ ¼ 2, 6, and 10, which correspond to a preference for an odd number of π 0 s, appear to have fit results slightly larger than MC.
The G-parity for χ cJ s is positive, suggesting that decays should favor an even number of πs, while Gparity for the J=ψ is negative, implying that decays favor an odd number of πs. These preferences in the distributions of the number of produced photons are observed for data, but MC simulation does not adequately reflect this.
While the agreement between data and MC simulation is reasonable at present, it should be improved for future studies of χ cJ decays and measurements of the ψð3686Þ → γχ cJ branching fractions with even larger data sets. This can be accomplished with further MC tuning or by weighting the present or future MC simulation to give better agreement with data. TABLE XXIII. Detected data events, D, efficiencies, ϵ, efficiency corrected events, N, and number of scaled simulated events N MC for χ c1=2 → γJ=ψ; J=ψ → anything. Here N has been multiplied by RðdataÞ and N MC by RðMCÞ, the fractions of valid π 0 s.