Inclusive Lambda_c Production in e^+e^- Annihilations at sqrt{s}=10.54 GeV and in Upsilon(4S) Decays

We present measurements of the total production rates and momentum distributions of the charmed baryon $\Lambda_c^+$ in $e^+e^- \to$ hadrons at a center-of-mass energy of 10.54 GeV and in $\Upsilon(4S)$ decays. In hadronic events at 10.54 GeV, charmed hadrons are almost exclusively leading particles in $e^+e^- \to c\bar{c}$ events, allowing direct studies of $c$-quark fragmentation. We measure a momentum distribution for $\Lambda_c^+$ baryons that differs significantly from those measured previously for charmed mesons. Comparing with a number of models, we find none that can describe the distribution completely. We measure an average scaled momentum of $\left= 0.574\pm$0.009 and a total rate of $N_{\Lambda c}^{q\bar{q}} = 0.057\pm$0.002(exp.)$\pm$0.015(BF) $\Lambda_c^+$ per hadronic event, where the experimental error is much smaller than that due to the branching fraction into the reconstructed decay mode, $pK^-\pi^+$. In $\Upsilon (4S)$ decays we measure a total rate of $N_{\Lambda c}^{\Upsilon} = 0.091\pm$0.006(exp.)$\pm$0.024(BF) per $\Upsilon(4S)$ decay, and find a much softer momentum distribution than expected from B decays into a $\Lambda_c^+$ plus an antinucleon and one to three pions.

(Dated: June 4, 2018) We present measurements of the total production rates and momentum distributions of the charmed baryon Λ + c in e + e − → hadrons at a center-of-mass energy of 10.54 GeV and in Υ (4S) decays. In hadronic events at 10.54 GeV, charmed hadrons are almost exclusively leading particles in e + e − → cc events, allowing direct studies of c-quark fragmentation. We measure a momentum distribution for Λ + c baryons that differs significantly from those measured previously for charmed mesons. Comparing with a number of models, we find none that can describe the distribution completely. We measure an average scaled momentum of xp = 0.574±0.009 and a total rate of N qq Λc = 0.057±0.002(exp.)±0.015(BF) Λ + c per hadronic event, where the experimental error is much smaller than that due to the branching fraction into the reconstructed decay mode, pK − π + . In Υ (4S) decays we measure a total rate of N Υ Λc = 0.091±0.006(exp.)±0.024(BF) per Υ (4S) decay, and find a much softer momentum distribution than expected from B decays into a Λ + c plus an antinucleon and one to three pions.

I. INTRODUCTION
The production properties of charmed baryons in e + e − annihilations into cc and in decays of bottom (b) hadrons probe different aspects of strong interaction physics. Experiments running at and below the Υ (4S) resonance are uniquely positioned to explore each of these processes in detail. The CLEO experiment has made precise studies of charmed mesons in this way [1], and the larger data samples available at the B factories have allowed improved studies of charmed mesons [2] and the first precise studies of charmed baryons [2,3].
Heavy hadrons (H) produced in e + e − annihilations provide a laboratory for the study of heavy-quark (Q = c, b) jet fragmentation, in terms of both the relative production rates of hadrons with different quantum numbers and their associated spectra. The latter can be characterized in terms of a scaled energy or momentum, such as x p ≡ p * H /p * max , where p * H is the hadron momentum in the e + e − center-of-mass (c.m.) frame and p * max = s/4 − m 2 H is the maximum momentum available to a particle of mass m H at a c.m. energy of √ s. For √ s ≫ 2m H it has been observed [4] that the x p distributions P (x p ) for heavy hadrons peak at relatively high values, and that very few b hadrons are produced apart from those containing initial b quarks, so that one can probe leading b-hadron production directly. This is also the case for charmed (c) hadrons when √ s < 2m B , where m B is the mass of the lightest b meson, but above BB threshold a large fraction of the c hadrons are b-hadron decay products.
Since the hadronization process is intrinsically nonperturbative, P (x p ) cannot be calculated using perturbative Quantum Chromodynamics (QCD). However, a high quark mass provides a convenient cut-off point and the distribution of the scaled momentum of the heavy * Also with Università di Perugia, Dipartimento di Fisica, Perugia, Italy † Also with Università della Basilicata, Potenza, Italy quark before hadronization, x Q ≡ 2p * Q / √ s, can be calculated [5,6,7,8]. The observable P (x p ) is thought to be related by a simple convolution or hadronization model. Several phenomenological models of heavy-quark fragmentation have been proposed [9,10,11,12]. Predictions depend on the mass of the heavy quark, with P (x p ) being much harder for b hadrons than c hadrons, and in some cases on the mass and quantum numbers of H. Hadrons containing the same heavy quark type are generally predicted to have quite similar P (x p ), although differences between mesons and baryons have been suggested [13,14]. Measurements of P (x p ) serve to constrain perturbative QCD and these model predictions. Furthermore, measurements for a given c or b hadron at different √ s can test QCD evolution, and comparisons of c-and bhadron distributions can test heavy-quark symmetry [15].
The inclusive b-hadron scaled energy distribution and its average value of 0.71 have been measured precisely [16] by experiments at the Z 0 , using partial reconstruction techniques. However, these techniques do not distinguish the different types of b hadrons. The relative production of B − u , B 0 d , B 0 s , excited b mesons and b baryons have been measured [4,17], but with limited precision and no sensitivity to differences in their x p distributions. Several c mesons have been studied at the Z 0 [18], but it is difficult to disentangle the leading charm and b-decay contributions and neither component is measured precisely. Recent measurements below BB threshold [1,2] have good precision over the full x p range and show substantial differences between the pseudoscalar D and vector D * meson states. P (x p ) has been measured below BB threshold for two charmed baryons, Λ + c by CLEO [19] and Belle [2], and Ξ 0 c by BABAR [3], but with limited statistics, especially at low x p . In this article we use the excellent particle identification of the BABAR experiment to isolate Λ + c baryons (the inclusion of charge conjugate states is implied throughout) in a large data sample collected below BB threshold, at √ s = 10.54 GeV. We measure P (x p ) precisely, and compare our results with available predictions and previous measurements of heavy hadrons.
The large b-hadron masses allow many hadronic decay modes, a small fraction of which have been studied in de-tail. The Υ (4S) resonance provides a unique laboratory, in which no b baryons are produced and decays of mesons into baryons can be studied directly. Many c baryons have been observed in inclusive Υ (4S) decays [4] and the low rate of associated leptons and high rate of "wrongsign" Λ + c [20] suggest interesting dynamics. However only a few exclusive decays with c baryons have been observed [21,22,23]. Again, momentum distributions have been measured only for the Λ + c [2,24] and Ξ 0 c [3] with limited precision. Here we use data collected on the Υ (4S) resonance ( √ s = 10.58 GeV) and subtract the e + e − → cc contribution measured at √ s = 10.54 GeV to make a precise measurement of the Λ + c momentum distribution in B meson decays, which we compare with a number of possible models.
In section II, we describe the BABAR detector, in particular the particle identification capabilities essential to these measurements. In sections III and IV, we discuss the selection of Λ + c candidates and the measurement of their x p distributions, respectively. We interpret the results for e + e − → qq events and Υ (4S) decays in sections V and VI, respectively, and summarize in section VII.

II. APPARATUS AND EVENT SELECTION
In this analysis, we use data samples corresponding to 9.5 fb −1 of integrated luminosity at √ s = 10.54 GeV and 81 fb −1 on the Υ (4S) resonance, √ s = 10.58 GeV. The BABAR detector is located at the PEP-II asymmetricenergy e + e − collider at the Stanford Linear Accelerator Center and is described in detail in Ref. [25]. We use charged tracks measured in the five-layer, double sided silicon vertex tracker (SVT) and the 40-layer drift chamber (DCH). In a 1.5 T axial magnetic field, they provide a combined resolution on the momentum p T transverse to the beam axis of [σ(p T )/p T ] 2 = [0.0013p T ] 2 + 0.0045 2 , where p T is measured in GeV/c. Charged particle identification uses a combination of the energy loss (dE/dx) measured in the DCH, and information from the detector of internally reflected Cherenkov light (DIRC). The DCH gas is helium:isobutane 80:20 and a typical cell size is 18 mm. A truncated mean algorithm gives a dE/dx value for each track with an average resolution of 7.5%, from which we calculate a set of five relative likelihoods L DCH i for the particle hypotheses i = e, µ, π, K, and p. Differences between the log-likelihoods l DCH are used as input to the particle identification algorithm.
The DIRC comprises 144 fused silica bars that guide Cherenkov photons to an expansion volume filled with water and equipped with 10752 photomultiplier tubes. The fused silica refractive index is 1.473, corresponding to Cherenkov thresholds of 128, 458 and 867 MeV/c for pions, kaons and protons, respectively. A particle well above threshold yields 20-75 measured photons, each with a Cherenkov angle resolution of about 10 mrad. A global likelihood algorithm considers all the recon-structed charged tracks and detected photons in each event and assigns each track a set of likelihoods L DIRC i . The DCH provides excellent K-π and p-K separation for momenta in the laboratory frame below 0.5 and 0.9 GeV/c, respectively. The DIRC provides very good separation for momenta above 1.0 and 1.5 GeV/c, respectively. In both detectors the separation power is lower for tracks with polar angles near 90 • in the laboratory than for more forward or backward tracks. To minimize the systematic errors in this analysis, the identification efficiencies must not vary rapidly as a function of momentum or polar angle. We have therefore developed an algorithm that uses linear combinations of l DCH ij and l DIRC ij chosen to minimize such variations. It is described in detail in Ref. [26] and its performance for tracks used in this analysis is shown as a function of momentum in Fig. 1. The identification efficiencies are better than 99% at low momenta and above 90% for the majority of Λ + c decay products. They are seen to vary smoothly with momentum, and are almost independent of polar angle except near 0.8 GeV/c (1.2 GeV/c) for pions and kaons (protons), where they are as much as 10% lower for central tracks than for forward/backward tracks. The misidentification rates depend strongly on polar angle in the momentum regions 0.6-0.8, 1.1-1.3 and 2.5-3.5 GeV/c. They are below 5% everywhere except that the rate for kaons (protons) to be misidentified as pions reaches 11% at 0.7 (1.2) GeV/c for the most central tracks. These rates have negligible effects on the results. About 16% of the selected tracks have good DCH information but are outside the DIRC fiducial acceptance. These can be identified with essentially the same efficiencies as in Fig. 1 for pion and kaon (proton) momenta below 0.6 (0.9) GeV/c.
The event selection requires three or more charged tracks in the event, which retains any e + e − → qq event or Υ (4S) decay containing a reconstructable Λ + c → pK − π + decay and suppresses beam-related backgrounds. We evaluate its performance using a number of simulations, each consisting of a generator for a certain type of event combined with a detailed simulation of the BABAR detector [27]. For e + e − → qq events we use the JETSET [28] generator and for Υ (4S) events we use our own generator, EVTGEN [29], in which the Υ (4S) decays into a BB pair, then the B and B decay using a combination of measured exclusive and semi-exclusive modes, and a b→ cW − model tuned to the world's inclusive data. We study large samples of simulated two-photon, τ -pair and radiative e-and µ-pair events, and find their contributions to both signal and background to be negligible.
We construct Λ + c candidates from charged tracks that are consistent with originating at the e + e − interaction point and have good tracking and particle identification information. Each track must have: (i) at least 20 measured coordinates in the DCH; (ii) at least 5 coordinates in the SVT, including at least three in the direction along the e − beam; (iii) a distance of closest approach to the beam axis below 1 mm; and (iv) a z-coordinate at this point within 10 cm of the nominal interaction point. These criteria ensure good quality information from the DCH and a well measured entrance angle into the DIRC. If the extrapolated trajectory intersects a DIRC bar then the track is accepted if it is identified as a pion, kaon or proton by the combined DCH and DIRC algorithm. If not, then it is accepted if it is identified as a proton (pion or kaon) using DCH information only and has a momentum below 1.2 (0.6) GeV/c.
We consider a combination of three charged tracks as a Λ + c candidate if the total charge is +1, one of the positively charged tracks is identified as a proton and the other as a pion, and the negatively charged track is identified as a kaon. With the appropriate particle type assigned to each track, we correct their measured momenta for energy loss and calculate their combined four momentum from their momenta at their points of closest approach to the beam axis. The distributions of invariant mass for the candidates in the on-and off-resonance data are shown in Fig. 2; Λ + c signals of about 137,000 and 13,000 decays, respectively, are visible over nearly uniform backgrounds.
The Λ + c reconstruction efficiency depends primarily on the momenta p and polar angles θ of the daughter tracks in the laboratory frame. To reduce systematic uncertainty, we apply an efficiency correction to each candidate before boosting it into the e + e − c.m. frame. The efficiencies for reconstructing and identifying tracks from pions, kaons and protons are determined from large con-trol samples in the data as two-dimensional functions of (p, θ). We use these efficiencies in dedicated simulations of qq and Υ (4S) events containing a Λ + c baryon that is decayed into pK − π + . From these we calculate the Λ + c selection efficiency ε as a smooth two-dimensional function of (p, θ) of the Λ + c . We check that the efficiency does not depend on other track or event variables, in particular that it is the same in simulated qq and Υ (4S) events for given values of (p, θ). The resolutions on the Λ + c momentum and polar angle are much smaller than the bin sizes used below, so we include resolution effects by defining the efficiency as the number of Λ + c reconstructed within a given (p, θ) range divided by the number generated in that range, using ranges smaller than the relevant bin sizes. We test the efficiency using a number of simulations, and find biases to be below 1%.
The efficiency varies rapidly near the edges of the detector acceptance and at very low momenta in the laboratory. We make the tight fiducial requirement that θ * , the polar angle of the Λ + c candidate in the e + e − c.m. frame, satisfy −0.7< cos θ * <0.2, which reduces model dependence and rejects all candidates in regions with efficiency below about 5%, including those with a total laboratory momentum below about 0.7 GeV/c. A feature of the boosted c.m. system is that true Λ + c baryons with low c.m. momentum p * are boosted forward and often have all three tracks in the detector acceptance, giving efficient access to the full p * range.
We define A k as the fraction of the Λ + c in events of type k produced within our fiducial range −0.7< cos θ * <0.2. In Υ (4S) → BB → Λ + c X decays, the true cos θ * distribution is uniform and A Υ = 0.45. In cc events the angular distribution of the initial c-quark follows 1 + cos 2 θ * , and we use the JETSET simulation to calculate the distribution for Λ + c after QCD radiation and hadronization. Soft Λ + c are produced predominantly in events with hard gluon radiation, which flattens the distribution considerably. The resulting value of A cc is 0.46 at p * = 0, and falls with increasing p * toward an asymptotic value of 0.38.
We bin candidates according to their reconstructed values of x p = p * /p * max , where p * max is calculated for each event from its c.m. energy and the nominal Λ + c mass [4]. Figure 3 shows the average value in each x p bin of the product A cc (p * ) · ε(p, θ) for selected candidates in the off-resonance data. It ranges from 8% at low x p to 19% at high x p . The error bars represent the statistical uncertainty on the efficiency calculation. The corresponding quantity for Λ + c from Υ (4S) decays, A Υ · ε(p, θ) , is slightly higher at low x p due to a small dependence of ε on θ, and rises faster with increasing x p since A Υ is constant, whereas A cc decreases. We give each candidate a weight equal to the inverse of either A cc · ε or A Υ · ε, as specified below. The RMS deviation of the weights in each bin is always much smaller than the average value.

IV. SIGNAL EXTRACTION
To estimate the number of Λ + c → pK − π + decays in each x p bin in the data, we fit the weighted invariant mass distribution with a function comprising signal and background components. Based on the simulated mass distributions, we describe the signal with a sum of two Gaussian functions of common mean value, one of which has 1.5 times the width and one quarter of the area of the other, and correct for a 1.3% residual bias in the fitted area. We check the simulated bias by comparing with a single Gaussian signal function. The change in the yields is 1.2% in both data and simulation. The simulation predicts a nearly uniform background over the pK − π + mass range shown in Fig. 2. We search for reflections in the data by changing the particle mass assignments. We observe signals for D + → K − π + π + (π + misidentified as p) and D + s → K − K + π + (K + misidentified as p) at very low levels consistent with the predictions of our detector simulation, but no unexpected structure. From these studies we calculate that reflections known to give broad structures in the vicinity of the Λ + c peak, such as D * + → K − π + π + (π + misidentified as p) contribute a number of entries in each bin much smaller than the statistical fluctuations. We also study processes such as Σ ++ c → Λ + c π + , with the wrong π + included in the Λ + c candidate, and find their contributions to be negligible.
In each x p bin, a linear function describes the mass distribution in the data over a wide range away from the Λ + c peak region, so we use a linear background function and perform fits over the range 2235-2335 MeV/c 2 .
We first fit the full data sample in each x p bin in order to study mass resolution and bias. These fits yield Λ + c mass values that vary slightly with x p in a manner consistent with the simulation and our recent measurement of the Λ + c mass [30]. The fitted mass resolutions (RMS width of the signal function) are shown as a function of x p in Fig. 4. The simulation is consistent with the data at low x p , and is slightly optimistic at high x p . The effect of this difference on the efficiency estimate is negligible.
Next we fix the mean and width of the signal function in each x p bin to values from linear parametrizations, and perform fits to the on-and off-resonance data separately. Dividing the signal yields by the integrated luminosity, bin width and branching fraction B pK − π + ≡ B(Λ + c → pK − π + ) = 5.0 ± 1.3% [4] gives the differential production cross sections shown in Fig. 5 with statistical errors only. We use the e + e − → cc acceptance factor A cc for the off-resonance data and, for purposes of this comparison, an average value of A in each bin weighted by the relative production rates measured below (see Tables I and V) for the on-resonance data. There are two broad peaks in the on-resonance cross section, corresponding to the contributions from Υ (4S) decays at low x p and from e + e − → cc events at high x p . For x p > 0.47, the kinematic limit for a B decay including a Λ + c and an antiproton, the two cross sections are consistent, indicating no visible contribution from Υ (4S) events.
We extract the cross section for Υ (4S) decays by repeating the analysis using A Υ for both data sets. In each x p bin we then subtract the off-resonance cross section, scaled down by 0.8% to account for the dependence of the cross section on the c.m. energy, from the onresonance cross section. We divide the off-resonance and Υ (4S) cross sections by the e + e − → hadrons and effective e + e − → Υ (4S) cross sections, respectively, to yield the differential production rates per event discussed in the following sections.
We consider several sources of systematic uncertainty. We propagate the uncertainties on the measured track finding and particle identification efficiencies by recalculating ε(p, θ) with all efficiencies of a given type varied simultaneously and repeating the fits. Tracking gives an uncertainty of 2.5% and the particle identification contributions total 1.5-2.1%, depending on x p . We obtain a 0.9% uncertainty due to the resonant substructure of the Λ + c → pK − π + decay similarly. The uncertainty on our integrated luminosity is 1.0%. Simulation statistics contribute 2-4% where the rate is significantly nonzero. We check the fitting procedure by floating the signal mean and/or width, fixing them to nominal or fitted values, using a single Gaussian signal function, using a quadratic background function, and varying the bin size. All changes in the signal yields are less than the corre- sponding statistical errors, and we take the largest change in each x p bin as a systematic uncertainty. Each simulated bias is varied by ±50%; together they contribute 0.7% to the systematic uncertainty. The imperfect x p distribution used in the efficiency calculation can affect the result if the efficiency varies over the width of an x p bin. We recalculate the efficiency with the input distribution shifted by plus and minus our bin width to derive a conservative limit on any such effect of 0.5-1.9%, depending on x p , which we take as a systematic uncertainty.
We also perform several systematic checks of the results. The cross sections measured separately for Λ + c and Λ − c , which have very different efficiencies at low laboratory momentum, are consistent. Cross sections measured in six different regions of cos θ * are consistent with each other. Due to the boosted c.m. system, these would be affected differently by any deficiency in the detector simulation, especially at low x p . They also have very different A cc values, and these studies indicate that the uncertainties due to both the production angle model in cc events and the value of the boost are negligible. The differential Λ + c production rate per hadronic event (1/N qq )(dN qq Λc /dx p ) is tabulated in Table I and compared with previous charmed baryon measurements in Fig. 6. We distinguish between systematic uncertainties that affect the shape of the cross section and those that affect only its normalization. The former include both uncertainties that are uncorrelated between bins and the parts of the correlated uncertainties whose values depend on  6: Differential Λ + c production rate per e + e − → qq event compared with previous measurements. The error bars include statistics and those systematic errors that affect the shape. Each experiment has a normalization uncertainty of a few percent, and there is an overall 26% uncertainty due to the Λ + c → pK − π + branching fraction.
x p . The uncertainty from the fitting procedure has negligible correlation between bins, and those from the particle identification and the shift of the simulated distribution have only very short-range correlations, so we include them in the uncorrelated category. We express the error matrix for the efficiency calculation as the sum of a diagonal "uncorrelated" matrix and a remainder matrix, in which the elements of the former are as large as possible but no correlation coefficient in the latter exceeds unity. The sum in quadrature of the uncorrelated uncertainties is listed in the "independent" column of Table I. It is typically 3% in the peak region, increasing to 10% where the cross section is one-third of its peak value, and becoming relatively large at the ends of the x p range. The square roots of the diagonal elements of the remainder matrix are listed in the "correlated" column of Table I, and included with the independent and statistical components in the error bars in the figures. All other uncertainties are fully correlated between bins and very nearly independent of x p , so are considered experimental normalization uncertainties. They total 2.9%, dominated by the track-finding efficiency. There is a 26% uncertainty on B pK − π + that also affects the normalization. The integral of the differential rate, taking the correlation in the errors into account, gives the total rate, listed at the bottom of Table I along with the normalization uncertainties. The product of the total rate per event and branching fraction of N qq Λc · B pK − π + = 2.84 ± 0.04 (stat.) ± 0.09 (syst.) × 10 −3 is consistent with, and more precise than, previous measurements. The normalization uncertainties are not included in any of the figures, and all rates shown assume the same value of B pK − π + .  Assuming the Λ + c are produced predominantly in e + e − → cc events, the total rate corresponds to a rate of N c Λc = 0.071±0.003 (exp.)±0.018 (BF) Λ + c per c-quark jet. Roughly 10% of the particles in high-energy jets have generally been observed to be baryons [4], and our measurement is consistent with 10% of c jets producing a c baryon, with a large fraction of these decaying via a Λ + c . All known non-strange charmed baryons decay predominantly through a Λ + c , whereas no Ω c states and only the heaviest observed Ξ c states [31] are known to decay through a Λ + c , so that 70-85% of inclusive charmed baryons would be expected to decay through a Λ + c in current hadronization models.
The shape of the differential production rate is consistent with previous results and measured more precisely. It is quite hard, as expected, peaking near x p = 0.6. We normalize the rate to unit area to obtain P (x p ), and compare it with previously measured distributions for Ξ 0 c baryons and D and D * mesons in Fig. 7. The Ξ 0 c distribution is normalized to have the same peak height as the Λ + c distribution, and, since it was measured on the Υ (4S) resonance, is shown only above the kinematic limit for B-meson decays. The two charmed baryons have similar distributions, with that for the heavier baryon shifted up in x p by roughly 0.05. Although qualitatively similar, the D meson distributions show broader peaks than the baryon distributions and differ greatly in the way they fall toward zero at high x p . The charmed baryon and meson distributions are all much softer than the inclusive B-hadron distribution at c.m. energies well above bb threshold, which peaks around x p = 0.75 [16].
The average x p value is often used in comparisons be- tween different heavy hadrons, and the higher moments of the distribution are of theoretical interest. In Table II we list values of the first six moments of the x p distribution, calculated by summing over bins. They are consistent with previous measurements from Belle [2]; all are 1-2 standard deviations lower, but the moments are strongly correlated with each other. The x p value of 0.574±0.009 is consistent with those measured [2] for D 0 and D + mesons, and about 5% lower than those for D * 0 and D * + mesons.

B. Tests of c-Quark Fragmentation Models
Testing models of heavy-quark fragmentation can be problematic since the predictions are usually functions of a variable z that is not accessible experimentally, such as z 1 = (E + p ) H /(E + p ) Q , z 2 = p H /p Q or z 3 = p H /p Hmax (p Q ), where p represents a momentum projection on the flight direction of the heavy quark before it hadronizes. Monte Carlo event generators use similar internal variables, and in some cases can be made to produce events according to a given input function f (z, β), where β represents the set of model parameters.
In this way one can test the large-scale features of any model, although the detailed structure may not be reproduced exactly. We consider the perturbative QCD calculations of Collins & Spiller (CS) [5] and Braaten et al. (BCFY) [7], as well as the phenomenological models of Kartvelishvili et al. for mesons (KLP-M) [9] and baryons (KLP-B) [13], Bowler [10], Peterson et al. [11], the Lund group [12], the UCLA group [14] and the HERWIG group [32]. The latter two include heavy quark fragmentation within their own generators, and the other seven predict the functional forms listed in Table III. We implement each of these functions f (z, β) within the JETSET generator. JETSET uses z 1 as its internal variable, but z 2 and z 3 are very similar at high x p where we are most sensitive to the shape. All distributions are affected by JETSET's simulation of hard and soft gluon radiation.
We test each model against our measured P (x p ) using a binned χ 2 where n is the number of bins, P data i (P MC i ) is the fraction of the measured (modelled) distribution in bin i, and V is the full error matrix formed from the errors on the data (Table I) with the statistical errors on the simulation added in quadrature to the diagonal elements. For each function in Table III we minimize this χ 2 with respect to the set of parameters β by scanning over a wide range of possible β values and generating a large sample of qq events at each of several points. We then generate additional sets near each minimum. All other model parameters are fixed to their default values.
The functions with one free parameter (CS, BCFY, KLP-M, KLP-B and Peterson) all show a single, well behaved minimum. The two parameters in the Bowler and Lund functions are strongly correlated, and the χ 2 shows a single narrow valley. The UCLA model has internal parameters a and b that control production of all particles simultaneously; since it has been suggested that their values should be different for mesons and baryons, we vary them by the same procedure, again finding a strong correlation and a narrow χ 2 valley. HERWIG has no free parameter controlling heavy hadron production, so we consider only the default parameter values.
We compare the fitted distributions with the data in Fig. 8, and list the χ 2 , the fitted parameter values and the average x p of each fitted distribution in Table IV  and Bowler models give the best descriptions of the data, with respective χ 2 confidence levels of 0.15, 0.11 and 0.06. However their fitted distributions are systematically below the data at the lowest x p values and above the data just below the peak region. The UCLA distribution is qualitatively similar to these three models but falls more rapidly at low x p , resulting in poor agreement with the data. The CS, BCFY and KLP-M models predict distributions that are much too broad, and the Peterson distribution is also too broad. The HERWIG distribution is consistent with the data in the peak region, but cuts off too sharply at high x p .
The fitted values of the parameter a for the UCLA, JETSET+Lund and JETSET+Bowler models are larger than those that describe the production of inclusive light hadrons and charmed mesons (a ≈ 1.2 for UCLA and 0.1 ≤ a ≤ 0.6 for the other two models). Differences between baryon and meson distributions have been suggested on the basis of quark counting [13,14]. Crossing of the diagram for leading hadron production in e + e − annihilation gives a deep inelastic scattering diagram, calculations for which depend on the number of spectator quarks, N s ; in the limit z →1, one expects f (z) ∝ (1 − z) 2Ns−1 [13,33], and 2N s − 1 = 3 for baryons. This is the form of the KLP-B function, which provides a much better description of the data than its counterpart for mesons. The UCLA model and the Lund and Bowler functions also contain (1 − z) a terms. For UCLA, the fitted value of a = 2.9 is close to 3, as anticipated [14].
The models predict P (x p ) for primary leading charmed hadrons, whereas the data also contain secondary charmed hadrons from the splitting of hard gluons, and some of the reconstructed Λ + c are decay products of other charmed baryons. Both of these effects are included in the JETSET, HERWIG and UCLA models, but it is important to consider the effects of possible mismodelling. The fraction of qq events containing a gluon splitting into a cc pair has been measured at √ s = 92 GeV to be about 0.01 [34]. At our lower c.m. energy this rate is ex- pected to be reduced, and the fraction of these that produce charmed baryons is expected to be lower than that for primary c and c quarks. The JETSET, UCLA and HERWIG models predict overall contributions of only (0.009±0.004)%, (0.017±0.005)% and (0.034±0.012)%, respectively, concentrated at low x p . The uncertainties in our lowest-x p bins are large enough to accommodate such a contribution. Adjusting the models to remove or double this contribution does not change the results of the fits significantly.
Currently there are three known Σ c states (with masses 2455, 2520 and 2880 MeV/c 2 ) that decay into Λ + c π, and four Λ c states (masses 2593, 2625, 2765 and 2880 MeV/c 2 ) that decay into Λ + c π + π − . In decays of such baryons into a slightly lighter baryon and one or two pions, the daughter baryon carries most of the momentum, so the effect of any such decay is to soften P (x p ) slightly without distorting it substantially. In the JETSET+Lund, JETSET+Bowler and UCLA models, this effect is partially compensated by the fact that the heavier baryons are generated with slightly harder distributions. Collectively, the effect is to broaden P (x p ) slightly and shift its average value down. Combining our Λ + c candidates with additional pions in the same event, we see clear signals for all of these states, and can compare the relative contributions to the detected Λ + c with the simulation. The largest contribution of about 7% from Σ c (2455) is well simulated in all models, but the Σ c (2520) rate is too high by a factor of ∼3, and the excited Λ c states are not in any simulation. Removing two-thirds of the Σ c (2520) narrows the simulated distributions slightly, but does not improve any χ 2 value significantly. Similarly, adding excited Λ c states broadens all distributions slightly, with no change in the conclusions of the model tests. No Ξ c or Ω c states are known to decay to Λ + c [4], except two recently reported by Belle [31]. The latter are observed at very low rates, so should have negligible effect on P (x p ).

A. Charmed Baryon Production
The differential production rate per Υ (4S) decay is shown in Fig. 9 and listed in Table V. The errors are as for the qq results, with an additional 1.5% normalization uncertainty due to the qq subtraction procedure. The kinematic limit for Υ (4S)→ BB→ Λ + c p decays is x p = 0.47, and above this value the rate is consistent with zero. This region is omitted from the table and only partly shown in the figure. We calculate the total rate by integrating the differential rate over the kinematically allowed bins. The resulting product of total rate and branching fraction, N Υ Λc · B pK − π + = 4.56 ± 0.09 (stat.) ± 0.31 (syst.) × 10 −3 , is consistent with previous measurements [2,24]. It corresponds to N Υ Λc = 0.091 ± 0.006 (exp.)±0.024 (BF) Λ + c per Υ (4S) decay, and 4.5±1.2% of B 0 /B − decays including a Λ + c baryon, assuming the Υ (4S) decays predominantly to BB.
Our results on the shape are consistent with, and more precise than, previous results, which are also shown in Fig. 9. The x p distribution is quite soft. In particular, the data drop rapidly above the peak and are consistent with zero above x p ≈ 0.35. This is the range expected for quasitwo-body decays into a Λ + c or Σ c plus an antibaryon such as a p, n or ∆, and includes much of the range expected for decays involving one or two additional pions. The B 0 → Λ + c p decay has been observed [21] at a very low rate consistent with our inclusive data.
A soft x p distribution was also seen in our recent study of Ξ 0 c baryons [3], but the statistics of that study did not allow a meaningful direct subtraction of the contribution from qq events. As an exercise, we assume a smooth distribution from qq events by choosing an empirical function that describes both the Λ + c data in Sec. V and the high-x p Ξ 0 c data (Fig. 7a). We fit this function to the high-x p Ξ 0 c data and subtract the result in all x p bins. The resulting approximate x p distribution for Ξ 0 c in Υ (4S) decays is also shown in Fig. 9, normalized to have roughly the same peak height as the Λ + c data. It is also quite soft, similar in shape to the Λ + c but shifted slightly downward in x p , and consistent with zero above x p ≈ 0.35. Because of the e + e − → cc subtraction procedure, the error bars cannot be compared with those for the Λ + c data, but the noted features do not depend on the details of the subtraction procedure. FIG. 9: Differential Λ + c production rate per Υ (4S) decay compared with previous measurements. Normalization errors are not shown, and the Ξ 0 c rate is normalized to match the peak Λ + c rate.

B. Model Tests
Existing models of B meson decays into c baryons were developed with little data on the x p distribution. In Fig. 10a we compare P (x p ) (normalized over the range x p < 0.475) with the predictions of the models in the JETSET generator [28] and our internal generator [29]. The former has been tuned to measured multiplicities and momentum distributions of stable B decay products, and predicts a distribution that is much too soft. The latter includes many measured exclusive decays involving D mesons and postulates analogous few-body decays involving charmed baryons. Its predicted x p distribution is similar in shape to the data, but shifted to higher x p values. Neither generator produces simple (quasi-)twobody decays such as B→ Λ + c p, Λ + c n, Λ + c ∆ or Σ c p. In Fig. 10b we compare the same data with simulated events of the type B→ Λ + c p(mπ) for selected values of m, the number of pions in the decay in addition to the Λ + c and antiproton. The distributions are insensitive to the charges of the pions or B meson, or to replacing the antiproton with an antineutron. Decays via a ∆ or strange antibaryon are not included; they give Λ + c distributions only slightly different from those shown with m − 1 pions. For m = 1, 2, 4 and 6, the distributions shown are from the JETSET simulation; phase space decays give similar distributions. The spread in the distribution for m = 0 is due to the finite momentum of the B meson in the e + e − c.m. frame. The measured P (x p ) is described adequately by the simulation with m = 4. Adding contributions from m = 3 and m = 5 improves the χ 2 of a comparison with the data, but no further contributions are helpful. That m is restricted to such a narrow range suggests that different types of decay modes are needed.
An intriguing possibility is that there is a large con- have recently been observed, and we have previously measured an unexpectedly high rate of inclusive Λ − c in B − decays [20]. As an exercise, we model two-body two-c-baryon decays based on the simplest internal W diagrams, i.e. those of the forms where the parentheses indicate Cabibbo-suppressed modes. Here, Ξ c represents any of the states Ξ c (2470), Ξ ′ c (2570) or Ξ c (2645), Σ c represents Σ c (2455) or Σ c (2520), Λ c represents Λ c (2285), Λ c (2593) or Λ c (2625), and we consider all kinematically allowed combinations at relative rates determined by phase space. We decay all Σ c baryons into Λ + c π and all excited Λ c baryons into Λ + c ππ, with 73% of the Λ + c (2593) baryons decaying through a Σ c π intermediate state and all others via phase space.
The x p distribution of the Λ + c from this simulation is compared with the data in Fig. 10c. Although it is too narrow to describe the data completely, it appears that such processes could contribute substantially to the overall rate. Combining this simulation with those for the B→ Λ + c p(mπ) modes and assuming a smooth, broad distribution of m, we can describe the data with as much as a 50% contribution from these two-c-baryon decays. For example, the "composite model" in Fig. 10c comprises 35% two-c-baryon decays and (12,25,12,9,7)% of m =(2, 3, 4, 5, 6), and provides an excellent description of the data. Decays with two charmed baryons and additional pions or kaons could also contribute at low x p , shifting the m distribution downward. Measurements of many exclusive baryonic B decays, including both oneand two-c-baryon modes, are needed to understand the dynamics in detail.

VII. SUMMARY AND CONCLUSIONS
We use the excellent tracking and particle identification capabilities of the BABAR detector to reconstruct large, clean samples of Λ + c baryons over the full kinematic range. We measure their total production rates and inclusive scaled momentum distributions in both e + e − → hadrons events at √ s = 10.54 GeV and Υ (4S) decays. Our results are consistent with those published previously and more precise.
In e + e − → qq events we measure a total rate per event times branching fraction into the pK − π + mode of N qq Λc · B pK − π + = 2.84 ± 0.04 (stat.) ± 0.09 (syst.) × 10 −3 , where the first error is statistical and the second systematic. The uncertainty on the total rate per qq event, N qq Λc = 0.057 ± 0.002 (exp.) ± 0.015 (BF), is dominated by the uncertainty on B pK − π + . The corresponding value of N c Λc = 0.071±0.003 (exp.)±0.018 (BF) Λ + c per c-quark jet is consistent with the hypothesis that roughly 10% of c jets produce a c baryon and a large fraction of these decay via a Λ + c . The scaled momentum distribution peaks at x p ≈ 0.6. It is similar in shape to those measured previously for D and D * mesons, but peaks more sharply and drops toward zero more rapidly as x p → 1. We measure an average value for Λ + c of x p = 0.574 ± 0.006 (stat.) ± 0.006 (syst.), which is consistent with values measured for ground state D mesons, but about 5% lower than those for D * mesons. We use this distribution to test several models of heavy quark fragmentation, none of which provides a complete description of the data. The baryon-specific model of Kartvelishvili et al. and the models of Lund and Bowler have acceptable χ 2 values, but all show a steeper slope on the low side of the peak than the data. The UCLA model shows similar qualitative features, but worse agreement with the data. The HERWIG model is far too narrow, and all others are too broad. In previous model tests using specific c mesons [2] and inclusive b hadrons [16] (a mix of roughly 90% mesons and 10% baryons), the Lund, Bowler and Kartvelishvili models generally gave the best description of the data, and UCLA described the b-hadron data, whereas the other models showed discrepancies similar in form to those reported here. The Kartvelishvili and UCLA models postulate different spectra for mesons and baryons. Their strong preference for their respective baryonic forms, combined with the observed differences in shape between our Λ + c spectrum and previously measured D meson spectra, indicate a difference in the underlying dynamics.
In Υ (4S) decays, we measure a total rate per event times branching fraction of N Υ Λc · B pK − π + = 4.56 ± 0.09 (stat.) ± 0.31 (syst.) × 10 −3 , corresponding to N B Λc = 0.045 ± 0.003 (exp.)±0.012 (BF) per B 0 /B − decay. The spectrum is softer than predicted by our B decay model, and much harder than that predicted by JETSET. It can be described by models of B→ Λ + c p(mπ) decays only if the m = 3-5 contributions dominate. Alternatively, a model including a large contribution from decays involving both a charmed and an anti-charmed baryon can describe the data in conjunction with a broad distribution of m. Additional studies of exclusive modes are needed to understand the details of B-meson decays into baryons.