Non-Relativistic Quantum Chromodynamics in Parton Showers

Measurements of quarkonia isolation in jets at the Large Hadron Collider (LHC) have been shown to disagree with fixed-order non-relativistic quantum chromodynamics (NRQCD) calculations, even at higher orders. Calculations using the fragmenting jet function formalism are able to better describe data but cannot provide full event-level predictions. In this work we provide an alternative model via NRQCD production of quarkonia in a timelike parton shower. We include this model in the Pythia 8 event generator and validate our parton-shower implementation against analytic forms of the relevant fragmentation functions. Finally, we make inclusive predictions of quarkonia production for the decay of the standard-model Higgs boson.


Introduction
The production of quarkonia, flavourless mesons comprised of a heavy quark and its corresponding antiquark bound by the strong interaction, is still not fully understood in the context of hadron-hadron collisions, despite decades of study .For a complete review, see ref. [22].The reason for this is that our theory of the strong interaction, quantum chromodynamics (QCD), is best understood at large energy scales, whereas the dynamics of quarkonia production spans both large and small scales.
One of the first models of quarkonium production to make quantitative predictions was the coloursinglet model (CSM) [23][24][25][26][27], which factorises these perturbative and non-perturbative regimes.In the CSM, the production of the heavy quark/antiquark pair (QQ) is calculated using perturbative QCD, and these quarks are then bound non-perturbatively.In this second step, the constituents are assumed to be at rest in the quarkonium frame, and the colour and spin state of the QQ pair does not change during the binding.Comparisons of CSM predictions to Tevatron data [7,8] demonstrated that the CSM, when calculated at leading fixed order, significantly underestimated the prompt production of ψ(1S) mesons, a cc 3 S 1 quarkonium state 1 .
The initial fixed-order CSM predictions [27] for the Tevatron data lacked two important sources of ψ(1S) production: one, partonic fragmentation; and, two, feed-down from heavier quarkonium states such as the 3 P J χ c mesons.This first contribution can be significant in hadronic collisions because the rate for inclusive charm production is much larger than that calculated only from fixedorder processes such as qq → cc and gg → cc.Fragmentation functions [28][29][30][31][32] were derived for c → cc 3 S 1 c, considering the simpler Z → cc process, and convolved with Tevatron jet cross-sections.This additional fragmentation contribution brought the CSM prediction into better agreement with data at high p T of the ψ(1S), but the overall agreement with data was still poor.For feed-down contributions from the χ c , predictions were made at fixed-order with the CSM [33][34][35].These calculations contain a kinematic cut-off to regulate an infrared divergence, typically taken as the mass of the bound state.The quantitative CSM results are very sensitive to the exact choice of this cut-off.Even after accounting for feed-down, the CSM prediction still significantly disagrees with the Tevatron data.
With the realization that the CSM was inadequate to describe the data, non-relativistic QCD (NRQCD) theory [36][37][38][39] was developed.While NRQCD is similar to the CSM in that it factorises into a perturbative and non-perturbative component, it does not require the initial partonic state to have the same colour or quantum numbers as the final state, and allows, for example, the production of colour-octet states which then evolve into physical colour-singlet states.
Production rates for quarkonia at both the Tevatron and LHC are, in general, well described by NRQCD.However, NRQCD predicts ψ(1S) polarisation [40][41][42][43][44][45] that is not observed in LHC data [16][17][18].Additionally, fixed-order NRQCD predicts relatively isolated ψ(1S) production, while LHC measurements demonstrate that ψ(1S) mesons are oftentimes produced in association with nearby particles.Recent developments in calculating quarkonium fragmentation functions promise to resolve these issues [46,47].However, a comprehensive approach for calculating these fragmentation functions is lacking.Here, we introduce an alternative approach of using the Pythia 8 [48][49][50][51] parton shower to calculate the production of quarkonia from fragmentation.The idea is inspired by an earlier study using the Ariadne program [52][53][54], but is implemented here with a complete and extensible set of quarkonia processes in the full Pythia 8 framework.
The remainder of this paper is organised as follows.In section 2 we briefly review NRQCD, while in section 3 we describe parton shower formalism.We then include NRQCD quarkonium fragmentation processes into the Pythia 8 parton-shower in section 4. In section 5, we present quantitative predictions using our approach, and finally, we provide conclusions in section 6.

Non-relativistic Quantum Chromodynamics
Non-relativistic QCD is an effective field theory where a perturbative expansion is performed not only in the QCD coupling α s , but also in the relative velocity between the two heavy quarks of the bound state.Including this relative velocity scale results in the expansion of the physical quarkonium state into non-physical Fock states.For a heavy quark state QQ, the general Fock-state expansion is, where v is the relative velocity between the constituent quark/antiquark, (1) in the superscript notates a colour-singlet state, and (8) in the superscript notates a colour-octet state.For a ψ(1S), this can be explicitly written as, and for the χ c states as, where J corresponds to 0, 1, and 2.
The heavy quarkonium cross-sections can then be factorised into short-and long-distance scales, which for initial partons a 1 and a 2 can be written as where the dσ(a 1 a 2 → QQ 2S ′ +1 L ′ J ′ (C ′ ) X) are perturbatively calculated partonic short-distance matrix elements (SDME) and ⟩ are empirically determined long-distance matrix elements (LDME).The sums are over the states of the Fock expansion of eq. ( 1), where the indices correspond to the spin and colour quantum numbers, and an SDME and LDME must be provided for each term.For hadronic beams h 1 and h 2 , the partonic cross-section is then convolved with the parton density functions (PDF) for the hadrons where the first two sums are over all possible partonic combinations a i and a j .The PDFs for the two hadrons are given by f 1 and f 2 , evaluated at the partonic energy scale µ with corresponding momentum fractions x 1 and x 2 .
The LDMEs of the O(1) leading colour-singlet Fock state, i.e. the Fock state with the same quantum numbers as the physical state, can be determined from the wave-function of the physical state at the origin.The LDMEs for the additional Fock states must be determined from fits to data, although velocity scaling rules can be used to relate these LDMEs.Typically, these fits are made to the differential quarkonium production cross-section with respect to the p T of the quarkonium, as the colour-singlet contributions dominate at low p T , while the colour-octet contributions dominate at high p T .Using LDME fits from CDF, the differential cross-section in p T of the ψ(1S) is described well at both the Tevatron [37,55] and the LHC [56][57][58].However, these NRQCD predictions are dominated by colour-octet contributions at high p T , which are typically transversely polarised, unlike colour-singlet contributions which are typically longitudinally polarised.This is in direct conflict with experimental results, which are primarily unpolarised across all p T of the ψ(1S) [15][16][17][18][59][60][61][62].
The relative isolation of quarkonia can also be measured experimentally, where a jet clustering algorithm such as anti-k T can be used to group the quarkonium with surrounding activity, and observables such as p T (quarkonium)/p T (jet) can be constructed.For fixed-order NRQCD predictions, quarkonia are expected to be relatively isolated for colour-singlet contributions.Even for fixed-order colour-octet contributions with a parton shower applied similar to that of Pythia 8, the quarkonia is still expected to be primarily isolated.However, measurements from both LHCb [19] and CMS [20,21] do not show this isolation, and instead match distributions that are more consistent with production of the quarkonia from fragmentation rather than fixed-order calculations.In refs.[47,63], fragmenting jet functions (FJF) calculated via resummation, and similar to the fragmentation probabilities of eq. ( 8) below, were used to model the LHCb data.However, these calculations are incomplete and do not include all sources of ψ(1S) production such as feed-down.
Until this work, Pythia 8 did not include NRQCD-based fragmentation contributions to quarkonium production.This contribution can be written, formally, as, which is analogous to eq. ( 5), but now uses the fixed-order partonic cross-section dσ(a i a j → a k X) for parton a i and a j scattering to produce parton a k inclusively, multiplied by the the fragmentation probability D k→QQ 2S+1 L J (z, µ), which can be for either the sum of all Fock states or a specific Fock state.The precise definition of z varies within the literature, but generally all definitions reduce to the relative energy fraction carried by the emitted quarkonium in the high-energy limit.While this formalism can calculate quarkonium production from fragmentation, it cannot produce realistic event topologies that emulate events observed in experiments.Additionally, this type of calculation typically does not include competition between quarkonium states, i.e. the fragmentation probabilities do not simply factorize.In this work, we provide an alternative method that overcomes both of these shortcomings, by introducing quarkonium splittings into the Pythia 8 parton shower.This technique is formally equivalent to analytic resummation for inclusive observables, but provides more comprehensive results via a flexible Monte Carlo generator interface.To avoid double counting, this NRQCD parton shower should not be used together with fixed-order quarkonium processes without considering appropriate matching2 .

Parton Showers
A complete description of the Pythia 8 parton shower is available in ref. [64].Here, a brief review is given, with a focus on the algorithmic procedure.Parton showers evolve partons produced from fixed-order calculations at high virtuality µ max to low virtuality µ min through a series of probabilistic parton emissions.While various parton showers differ in the choice of the evolution variable µ, it is illustrative to consider it as an energy-like scale.The probability for parton a i to not radiate, i.e. a i → a j a k , between virtualities µ 1 and µ 2 is given by, where z is defined as E j /E i , the fractional energy parton a j takes from a i , and P i→jk is the splitting kernel for the process a i → a j a k .A summation is performed over all available splitting kernels, e.g. for a gluon both the splittings g → gg and g → qq would be summed, where the second splitting is over all kinematically available quark flavours.The parton shower then proceeds through an iterative algorithm.For a given parton dipole with virtuality µ 1 , a no-emission probability is uniformly sampled, ∆ i (µ 1 , µ 2 ) ∈ U(0, 1), yielding an equation that can be solved for the emission scale µ 2 .In practice, the integrands in eq. ( 9) are complicated, and so the emission scale is calculated using an overestimate and rejection method called the veto algorithm.In this case, the summed integrands P i→jk of eq. ( 9) are replaced by O i , is valid for all µ in the range µ 2 to µ 1 , and all z in the range z min to z max .If α s is evolved at first order, the new virtuality µ 2 can be calculated as where b 0 is the one-loop β-function coefficient [65], (33 − 2n f )/6, and Λ is the scale from which α s is being evolved.If, instead, α s is fixed at one value, the new virtuality becomes, After the virtuality µ 2 is determined according to eq. ( 11) or eq.( 12), one of the quarkonium splittings a i → a j a k is chosen randomly with probability O i→jk /O i .A z is then randomly sampled3 between z min and z max .The splitting is accepted if, where the right-hand side is a uniformly sampled random variable.If the splitting is rejected, then: the algorithm returns to eq. ( 11) or eq.( 12); a new µ 1 is set as the old µ 2 , µ 2 → µ 1 , and a new µ 2 is selected; and the procedure continues as before.The parton shower continues until µ 2 drops below the minimum virtuality, µ min .The virtuality variable µ in Pythia 8 [64] is defined as, where √ s is the centre-of-mass energy for the dipole, m r is the mass of the radiating parton in the dipole, and z is a light-cone momentum fraction carried by the emitted parton4 .

Implementation
A full set of splitting kernels to complement the ψ(1S) and χ c fixed-order production implemented in Pythia 8, described in section 2, were derived from fragmentation functions in the literature similar  (54,55,(59)(60)(61) to those used in eq. ( 8).The implemented splitting kernels, as well as their references, are given in table 1.The full forms for all the implemented splitting kernels can be found in the references.To illustrate our procedure for including the NRQCD-based fragmentation functions in the Pythia 8 shower, we give the details for the Q → QQ 3 S 1 (1) Q splitting.The splitting kernel is given as, where the LDME ⟨O QQ 2S+1 L J 3 S 1 (1) ⟩ depends upon the final physical state QQ 2S+1 L J that this splitting kernel is being used to produce, e.g. the ψ(1S) as a cc 3 S 1 state.
This expression differs from that given in ref. [30] for a number of reasons.First, the z convention used there is equivalent to 1 − z in Pythia 8, and so a change of variables is made here.An additional factor of 2π is also included to account for the 1/(2π) in eq. ( 9).There is already a factor of α s (µ) in eq. ( 9), so only one factor of α s is included, evaluated at the scale p T,evol ; in ref. [30], both factors of α s are evaluated at a scale of 3m Q .The choice for this scale is discussed in section 5. Finally, the integration variable of eq. ( 9) is µ 2 ≡ p 2 T,evol , rather than s as provided in ref. [30], so a Jacobian transforming from s to p 2 T,evol has been introduced: Because there is already a factor of 1/p 2 T,evol in eq. ( 9), only the additional factor of (s − m 2 Q ) needs to be included, corresponding to the final line of eq. ( 15).A similar conditioning must be done for all the splitting kernels of table 1.
An overestimate must also be determined for eq.( 15).This can be determined analytically or numerically.For this splitting, we can factorise eq. ( 15) into a z-independent component, and a z-dependent component f (z), given by, A factor of 1/(s−m 2 Q ) 2 is kept in the z-dependent component f (z) to ensure it remains dimensionless.The overestimate can then be written as, where we have used the numerical knowledge that f (z) < 3/2, and we use the parton shower cut-off to evaluate α s .Since (s − m 2 Q ) must always be larger than (m Q in the numerator of the second term.Again, a similar procedure can be performed for the remaining splitting kernels of table 1.Some of the splitting kernels above deviate significantly enough from the example to warrant special mention.The g → QQ 3 S 1 (1) gg kernel is a 1 → 3 splitting, rather than the normal 1 → 2, and so the splitting kernel depends upon another kinematic variable, the digluon mass m 2 gg , which must be sampled when generating the final-state kinematics.Another factor of α s enters into this splitting kernel, which is naturally evaluated at a scale of p T,evol with no choice given to the user.The g → QQ 3 P J (1) g splitting kernel diverges rapidly as p 2 T,evol → 0, and so a p 2 T,evol -dependent overestimate has been introduced to increase the sampling efficiency for this kernel.Finally, it should be noted that table 1 does not include any splittings for the QQ 3 S 1 (8) Fock state.These must be treated specially as described in the following.
There is one QQ 3 S 1 (8) splitting, given by g → QQ 3 S 1 (8) .This is unique because the splitting is 1 → 1, rather than 1 → 2, and arises from a delta function in the splitting kernel.Because the gluon has the same quantum numbers as the QQ 3 S 1 (8) , this splitting can be thought of as an off-shell gluon converting directly to a QQ 3 S 1 (8) state when the virtuality of the gluon matches the corresponding quarkonium mass.The splitting kernel is then given by, where δm is a small mass splitting, which for this work is set as 0.01 GeV.The results of the parton shower are not sensitive to the choice of this mass splitting.

Results
The differential fragmentation probability from which the splitting kernel of eq. ( 15) was derived can be integrated over s to determine the total fragmentation probability.From ref. [30], this is for ψ(1S) production from the c → cc 3 S 1 (1) c splitting at the lowest allowed scale of 3m c .Here, z does not follow the same convention as sections 3 and 4, but rather is defined as 2E QQ 2S+1 L J / √ s where √ s is the centre-of-mass energy for the initial dipole.A comparison can be made directly between the parton-shower implementation of this work, and eq.( 20).This is done by producing a cc dipole with a sufficiently high √ s that mass effects are negligible.The Pythia 8 parton shower, with the inclusion of the individual quarkonium splitting kernels from section 4 with no competition between the quarkonium splittings, is then run on this dipole.All QCD splittings are also switched off.The result can be seen in the left plot of fig. 1, where the Pythia 8 result given by the red histogram matches well with the analytic prediction of eq. ( 20) corresponding to the dashed line.Similar agreement is also seen for the η c , shown in blue.
The DGLAP equations can be applied to eq. ( 20) to numerically evolve the fragmentation probability from the scale of 3m c up to m Z /2.In ref. [30] this is provided numerically5 for the ψ(1S) and η c , as given by the dashed red and blue lines, respectively, in the right plot of fig. 1.The implementation of this work can be compared by showering a cc dipole with √ s = m Z , and keeping all QCD splittings switched on.The results are also given in the right plot of fig. 1.It should be noted that for high z, where mass effects become negligible, the results match well.However, at low z the results differ significantly, because in ref. [30] the masses of the ψ(1S) and η c were not considered.Accounting for this mass can be seen clearly in the Pythia 8 result, where there are no quarkonia produced with z < 2m ψ(1S) /m Z .
Analytic expressions similar to eq. ( 20) are available in the literature for all the splittings of table 1, evaluated at the lowest kinematically allowed scale.For heavy-quark radiators Q, this scale is 3m Q , while for gluon radiators g, the scale is 2m Q .In the left plot of fig. 2 the fragmentation probabilities D g→ψ(1S) (1) (z, 2m c ) and D g→η (1) c (z, 2m c ) for the colour singlet splittings g → cc 3 S 1 (1) g and g → cc 1 S 0 (1) g are given, respectively.There is good agreement between the parton-shower results and the analytic predictions.In the right plot of fig.2, comparisons between the parton-shower results and analytic predictions are given for D g→χ (z, 2m c ), corresponding to the splittings g → cc 3 P 0 (1) g, g → cc 3 P 1 (1) g, and g → cc 3 P 2 (1) g, respectively.All three parton-shower results agree well with their corresponding analytic predictions.
The fragmentation probabilities, evaluated at the energy scale of 3m c , for the production of the χ c from c-initiated splittings, are given in fig. 3. On the left, the fragmentation probabilities D c→χ (1) cJ (z, 3m c ) are given for the colour-singlet splittings c → cc 3 P J (1) c.On the right, the fragmentation probabilities D c→χ (8) cJ (z, 3m c ) are given for the colour-octet splittings c → cc 3 S 1 (8) c.All six parton-shower results are in good agreement with the analytic predictions.
Fig. 3 Production of (left) colour-singlet and (right) colour-octet P -wave states from charm splittings compared between (solid) Pythia 8 and (dashes) analytic expressions at the energy scale of 3mc.Analytic predictions evolved up to the scale of m Z are not in the literature for any of the processes other than c → cc 3 S 1 (1) c and c → cc 1 S 0 (1) c for the ψ(1S) and η c , respectively.Consequently, we are not able to make comparisons similar to the right plot of fig. 1.We instead provide results at the energy scale of m Z /2 in fig. 4 for all the remaining splittings implemented in Pythia 8.The fragmentation probabilities for all the colour-singlet c → cc 2S+1 L J (1) c splittings are given in the upperleft plot, while the equivalent fragmentation probabilities for the colour-octet c → cc 3 S 1 (8)   Fig. 5 Comparison of (left) splitting kernel choices for the colour-octet cc 3 S1 (8) splitting with Pythia 8 at the energy scale of m Z /2, and (right) scale choices for the colour-singlet cc 3 S1 (1) splitting.
of the χ c states are given in the upper-right plot.The fragmentation probabilities for all the coloursinglet g → cc 2S+1 L J (1) g splittings are given in the lower left, and the fragmentation probabilities for the colour-octet g → cc 3 S 1 (8) splittings of eq. ( 19) are given in the lower right.
There is ambiguity in how the colour-octet QQ 3 S 1 (8) states, lower-right plot of fig. 4, are treated in the parton shower after being produced.Because colour-octet states carry colour charge, they must undergo further QCD radiation before transitioning to their physical state, requiring the emission of one or more soft gluons.The default behaviour in Pythia 8, prior to this work, was to treat these states with twice the q → qg splitting kernel, i.e. the two heavy quarks of the quarkonium state would both be treated as radiators.However, it may be more appropriate to consider this colour-octet state as a massive gluon, where the splitting kernel as derived from ref. [68] is where the radiating gluon has mass m g and the emitted gluon is massless.These splittings do not guarantee the gluon emission necessary to transition from an octet to singlet state, so the colour-octet states are required to have slightly larger masses than their corresponding physical states.They are then forced to transition to their physical states through the isotropic emission of a single soft gluon.The difference between the two splittings of eqs.( 21) and ( 22) depends on when the gluon splitting producing the colour-octet state, g → QQ 3 S 1 (8) , occurs in the parton-shower evolution.Prior to this splitting, the gluon can continue to radiate, g → gg.Consequently, the second of the two splittings results in a transition much later in the parton-shower evolution, driven primarily by g → gg splittings.The effects of these two different splitting kernels, as well as the pre-factor used for eq.( 21), can be explored by creating an initial dipole with a ψ(1S) 3 S 1 (8) state at an √ s of m Z , and evolving this downward to the parton shower cut-off scale.The results for this configuration are shown on the left of fig. 5.The difference in the energy radiated away by the colour-octet state depending on the splitting kernel is sizeable.For the previous Pythia 8 default splitting of q → qg with a pre-factor of 2, the colour-octet state remains relatively isolated.If the pre-factor is unphysically doubled again to a pre-factor of 4, the colour-octet state radiates significantly more, as expected.However, for the splitting g → gg, the colour-octet state radiates even more.In this work, this gluon splitting kernel has been used throughout, i.e. to produce the bottom-right plot of fig. 4.
As mentioned in the discussion of eq. ( 15), some choice needs to be made for the scale at which α s is evaluated.Given the diagrams for most quarkonium production, it is natural to choose the first scale as p T,evol , or whatever the virtuality of the parton shower is.However, the second choice of scale is not as clear.Here, we consider three different options for µ: the mass of the quarkonium m, the parton-shower virtuality p T,evol , and the centre-of-mass of the dipole √ s.In the right plot of fig. 5, the impact of this choice is shown for ψ(1S) production from a cc dipole with √ s = m Z , using only the c → cc 3 S 1 (1) c splitting kernel.The default choice of p T,evol is given in blue, and is identical to the corresponding blue curve in the right plot of fig. 1. Switching to m provides a nearly identical result, while using √ s significantly reduces the fragmentation probability, although the shape remains similar, as expected.It should be noted for the g → QQ 3 S 1 (1) gg splitting, only one of the three factors of α s is effected by this scale choice, the other two are always evaluated at p T,evol .
Up to this point, we have only considered parton-shower configurations where no competition is included between the various quarkonium splittings.This competition is critical when generating comprehensive quarkonium predictions.Because of this competition, the different fragmentation contributions do not simply factorize into independent pieces.To demonstrate the full flexibility of our implementation, we calculate the inclusive branching fractions for standard-model Higgs boson decays into quarkonia.When performing this calculation, all possible quarkonium splitting kernels are switched on with full competition, as well as all other possible splittings, i.e.QCD and electroweak.The results are given in table 2. Here, the branching fractions in the first set of columns is separated into quarkonium production from parton-shower feed-down and from non-quarkonium hadron decays, e.g.b hadrons.The second set of columns is production from the parton shower, including feed-down, separated by the splitting type.Here, X (8) is any colour-octet production, g → X (1) is colour-singlet production from a gluon splitting, and Q → X (1) is colour-singlet production from heavy-quark splitting.
For ψ(1S) states, production from non-quarkonium hadron decays dominates the branching fraction, while the parton-shower contribution is nearly two orders of magnitude smaller.This is similar for all the charmonium branching fractions.For bottomonium production there is no production from non-quarkonium hadron decays, as there are not sufficiently high-mass states.For the η b and Υ, the primary production mechanism is colour-singlet production from b splittings, whereas for the χ b states, the primary production mechanism is from colour-octet splittings.

Conclusions
We have implemented a full-featured quarkonium parton shower in the Pythia 8 Monte Carlo event generator, available from version 8.310 onward.This parton shower can be used to model the LHC quarkonium isolation results, as well as explore the production of quarkonia in general.Our implementation includes all the relevant splitting kernels needed for the production of QQ 1 S 0 , QQ 3 S 1 , and QQ 3 P J states with arbitrary radial excitations and for Q ∈ c, b.These splitting kernels have been validated against fragmentation probabilities available in the literature, whenever possible, and found to be in good agreement.We provide first predictions for a number of fragmentation probabilities evolved up to the energy scale of m Z /2.The framework is flexible, so that an arbitrary number of competing splitting kernels can be included.Various choices are available for the radiation of the colour-octet states, as well as a choice of scale for evaluating α s .Finally, we provide comprehensive first predictions for the production of QQ 1 S 0 , QQ 3 S 1 , and QQ 3 P J states from decays of the standard-model Higgs boson.Future work includes developing a method to match these results from the parton shower with fixed-order predictions, including matrix-element corrections, and introducing additional splitting kernels, such as those needed for double-heavy meson and baryon production.

Fig. 1 Fig. 2
Fig. 1 Production of colour-singlet S-wave states from charm splittings compared between (solid) Pythia 8 and (dashes) analytic expressions at the energy scales of (left) 3mc and (right) m Z /2.

Fig. 4
Fig. 4 Production of (left) colour-singlet and (right) colour-octet states from (upper) charm and (lower) gluon splittings with Pythia 8 at the energy scale of m Z /2.

Table 1
Summary of the implemented quarkonium splitting kernels, with their corresponding reference from the literature.The numbers in parentheses indicate the relevant equations from that reference.

Table 2
Branching fractions for Higgs decays into quarkonia calculated with Pythia 8. Here, feed-down is from quarkonium decays produced by the parton shower, while hadrons is production from all non-onia hadron decays, e.g. from b hadrons.The second set of columns is production from the shower separated by splitting type.