A statistical model to calculate inclusive hadronic cross sections

Hadronic cross sections are important ingredients in many of the ongoing research methods in high energy nuclear physics, and it is always important to measure and/or calculate the probabilities of different types of reactions. In heavy-ion transport simulations at a few GeV energies, these hadronic cross sections are essential and so far mostly the exclusive processes are used, however, if one interested in total production rates the inclusive cross sections are also necessary to know. In this paper, we introduce a statistical-based method, which is able to give good estimates to exclusive and inclusive cross sections as well in the energy range of a few GeV. The method and its estimates for not well-known cross sections, will be used in a Boltzmann-Uehling-Uhlenbeck (BUU) type off-shell transport code to explain charmonium and bottomonium mass shifts in heavy-ion collisions.


Introduction
Hadronic cross sections at a few GeV are usually used as input parameters in heavy-ion transport codes, where up to 10 − 20 GeV the main degrees of freedom are the known baryons, mesons and the well-established resonances [1][2][3][4]. In these simulations mostly 2 → 2 exclusive reactions are used as input, however there were some attempts to include N → M reactions as well [5,6]. It is not straightforward how to include these types of collisions into the codes, so these are mostly omitted in the simulations. At higher energies, partonic degrees of freedoms are also becoming important [7]. Recent calculations show that the charmonium states e.g. J/Ψ, Ψ (3686), Ψ (3770) will acquire mass shifts during their evolution in dense matter, which could indicate a non-zero value for the gluon condensate at specific nuclear densities [8][9][10]. This information, if measured could be really important to non-perturbative quantum chromodynama e-mail: balassa.gabor88@gmail.com (corresponding author) ics. For this purpose the inclusive charmonium production cross sections should be known, however, they are not well (or not at all) measured in the energy range we are interested in (2)(3)(4)(5)(6)(7)(8)(9)(10), so it has to be estimated from some theory. The method presented in this paper is able to predict inclusive cross sections, which will be ultimately used to estimate charmonium production cross sections. The paper is organized as follows. In sections one and two a basic formulation of the method and its capabilities to estimate exclusive cross sections are described. This can be found more rigorously in [11], where the main idea and some examples were also given. After the basic formulation, section three dedicated to the normalization and model error estimation. The fourth section describes a method, which is used to calculate inelastic cross sections without summing over all the possible 2 → N reactions and using only 1-, or 2 fireball event ratios. The model error is also addressed in this section. The fifth section divided into two subsections. In the first subsection six examples are given for the inclusive cross sections compared to their measured values, while the second subsection corresponds to charmonium production and the estimation of the charm quark creational probability. Finally the last section briefly summarizes the paper.

Basic formulation and normalization
The method is based on the so-called Statistical Bootstrap [12][13][14] approach, which idea was widely used to estimate particle multiplicities at high energies [15][16][17]. In our model, we assume that at the initial state of a collision a strongly interacting, compound system is formed (a so-called fireball), which will ultimately decay into hadrons giving some specific final state of the collision. In the reaction probability, the initial fireball formation stage and the hadronization stage can be separated as it is shown in Eq. 1, where we assume that the first stage is described by the total or inelastic cross sections 1 of the reaction and the second stage is described by some mixed statistical and dynamical factors.
where E is the CM energy of the collision (E = √ s), σ n→k is the generalized n-body cross section, p 1 , ... p n are the momenta of the incoming particles, q 1 , ...q k are the momenta of the outgoing particles, R(E, ...) is the function describing the initial state dynamics, and w(E, ...) describes the probability of a specific k-body final state. At the first stage of the hadronization process the fireball with invariant mass M = √ s directly giving hadrons or decay into two subsequent fireballs with smaller masses m 1 , m 2 where m 1 + m 2 = M should be fulfilled. The resulting fireball will hadronize or decay into further fireballs. At the end of the chain, all of the fireballs have to hadronize leaving only hadrons at the end. The probability of the one-, two, three-, etc... fireball scheme can be calculated with an appropriate model described in our previous paper on this topic [11]. The possible number of hadrons coming from one fireball is two or three 2 with the probabilities of P d 2 = 0.69, P d 3 = 0.24 which is described in the original formulation of the statistical bootstrap by Frautschi [18].
An extra ingredient to the full probability (not mentioned in the original formulation, because the normalization was cancelled in the calculated ratios), which is needed to make a proper normalization is to separate stable and nonstable/resonant particles and use a weighting factor to each of the channels corresponding to the full width of the particles. In the model, the stable particles are the ones, which could not decay into other hadrons via the strong interaction. Examples are the N * nucleon resonances, which also could decay strongly to protons/neutrons, so they are not stable. These unstable particles are weighted by a relativistic Breit-Wigner factor [19]: where m i is the i'th particle mass, Γ i is its total width, E i is the invariant mass of the resonance, and the energy integral The appearance of the Breit-Wigner factor can be explained by the propagator of a particle with finite but non negligible width, when calculating Feynman diagrams e.g. as it was described in [20]. To calculate it's contribution to the full probbility it has to be integrated out in the physical region, defined by the masses of the possible decay products of each resonance.
Putting together everything the full probability of a onefireball scheme, which hadronizes into two or three particles can be calculated as it is shown in Eq. 3: where P f b 1 is the probability of the one fireball decay scheme, C Q i is the quark combinatorial probability for final state i, which is proportional to the number of quark-antiquark pairs created at energy E, and to the u, d, s, c, ... quark creational probabilities fitted from measured data. Furthermore, n i stands for the number of hadrons coming from the i'th fireball (2 or 3) and P H,1 n i is the hadronization probability factor given by Eq.4 and Eq.5 for two-and three-body states respectively: where P d 2 = 0.69 and P d 3 = 0.24 are the probabilities of a two or three particle final state decays, s l is the spin of the l'th particle, N I is the number of same particles in the final state, m i are the masses of particles, Φ n (E) is the n-body phase space integral, ans ρ(E) is the density of states (DOS) given by the statistical bootstrap [21,22] : where a,E 0 and T 0 are free parameters, however a always factorized out due to the normalization and E 0 = 500 MeV is previously fitted to the experimentally measured DOS [23]. In our previous paper we also set T 0 = 160 MeV using calculated and measured cross section ratios. All of the factors in Eq. 3 are described in detail in [11]. To simplify the following sections we have introduced the notation: The two-and three-body phase-spaces are weighted by the product of the particles Breit-Wigner factors (if it is a resonance) and integrated out to the corresponding energy region as it is seen in Eq. 8: where q j is the three momenta of the j'th particle, E i is the energy of the i'th particle and F B R r (E r , m r ) is the Breit-Wigner probability factor of the resonance r , and V is the interaction volume. If all of the particles are stable then the integrals in the first bracket and the Breit-Wigner factors are not needed. Otherwise the r index goes through all of the resonances. The Breit-Wigner factor will decrease the phase-space according to the resonance width e.g. for a threebody channel with one resonance and two stable particles (m 1 = m 2 = m 3 = 1 GeV and Γ 1 = 0.2 GeV) the corresponding phase-space ratio compared to the all stable particle configuration can be seen in Fig. 1.
The phase-space with a resonant particle is decreased accordingly to the Breit-Wigner factors inside the phasespace integrals.
The normalization sum in Eq. 3 goes through all the possible one-fireball decays which are allowed by the initial-state quantum numbers. This implies that the normalization sum is different for each reaction e.g. it is not the same for pp and π − p reactions due to their different quantum numbers. The extension to more fireballs is not straightforward, due to the different possibilities of how to include quantum number conservation into the model. In our normalization scheme, only the full final state has to respect the conservation laws and the quantum numbers of the different fireballs could be anything, there are no restrictions to it. The normalization procedure for more than one fireball is straightforward to see from the two fireball probability in Eq. 9: where the normalization sum Z 2 (E) is described by Eq. 10: and i j means a specific hadronic channel for each fireball e.g. i = π + π 0 for the first fireball, and j = pn for the second fireball, so the expression gives us the X → π + π 0 pn probability. P f b 2 is the probability of two fireball formation, x min = 2m π 0 is the threshold energy, i, j, k, l means a specific final state e.g. i = ppπ 0 so n i = 3, N i, j = 1 if the two fireballs give different final states, and N i, j = 2 if the two fireballs are the same i = j. In the normalization sum in the denominator kl ∈ S means it goes through all the possible final states, whose quantum numbers equal the initial state quantum numbers. The normalization sum for W 1 goes through all the possible two-, and three-body final states regardless of the quantum numbers of the initial state. A short example is the pp collision, where we have S = (0, 0, 0, 0) (baryon number, charge, strangeness, charmness). For the two fireball normalization sum, we only consider the channels which at the end give back the initial quantum numbers e.g. Q 1 = (1, 0, −1, 0) and Q 2 = (−1, 0, 1, 0) is a good choice as S = Q 1 + Q 2 . So separately the fireballs do not have to respect the conservation laws, but the whole system has to and give back the initial quantum numbers. It is important to note that the N-fireball probability can be built up by the 1-fireball factors and their corresponding integrals as it can be seen in Eq. 11 for the general case with N-fireballs.
where Z k (E) is the energy-dependent k-fireball normalization factor given by Eq. 12.
The complicated normalization factors can be dropped if the ratio of the two process probabilities are taken with strictly the same number of fireballs. If one channel could come from more fireball decays e.g. 2π + 2π − , which could come from one fireball X → ρ 0 ρ 0 and two fireball decay as well (X 1 , X 2 ) → (2π + , 2π − ) the normalization factors do not drop out, however, if we have two channels, which exclusively come from 1 fireball decay e.g. X → π + π − and X → nn, then taking the ratio of the two probabilities the normalization factor is dropped out. If we take the measured cross section for one of the channels and multiply it with the calculated ratio from the model, we get back the unknown cross section of the other channel. This practically only works for low multiplicities, as from one fireball the maximum number of hadrons, we can get is very limited.
To conclude this section the model is shown to be able to give reasonably good estimates to cross sections with calculating the probability ratios as it is shown in our previous paper. The only free parameters in the model are the quark creational probabilities, the interaction volume V , and the hadronization temperature T 0 . The density of states contains two more previously fitted parameters a, E 0 , where a can always be factored out by normalization, and E 0 is fixed by fitting to the experimentally observed DOS. In our previous paper, we fixed the temperature to T 0 = 160 MeV by calculating cross section ratios. The interaction volume is set to a corresponding interaction radius of r = 0.5 fm, which describes the calculated inclusive cross sections better at higher energies than a larger value. The quark creational probabilities are hidden in the energy-dependent quark combinatorial probability factor. Their fitted values in the energy range of E = 2 − 8 GeV for the three lightest quarks are P u = P d = 0.425 and P s = 0.15. These fits were made without including charmed particles, however due to the smallness of the charm creational probabilities, these values are approximately still valid in the energy range, we are interested in. For the heavier charm quarks and the corresponding charmonium particles, the fit P c is not straightforward as the experimental cross sections are very limited and the cross sections are too small to use the same method which is used to fit P u , P d , P s . The solution to this problem is described in Sect. 5.2.
There are some limitations to the model, however, which has to be addressed. It is not suited to describe the elastic part of the two-body reactions, so that channel is simply omitted from the calculations and the inelastic ratio is used instead of the total cross section. Another limitation is that the model cannot describe the A + B → C reactions, where A, B, C are some particles. This is the consequence of the Frautschi picture of the Bootstrap, where the minimum number of final state hadrons are 2. This problem is solved by introducing an extra contribution to the cross section described by a Breit-Wigner cross section formula σ B R (E): where A, B are the colliding particles, C is the created resonance, which could decay into some specific final state described by it's total decay width Γ tot . The s A,B,C factors are the spins of the corresponding particles, and p i is the center of mass momentum of the initial state. The relevant decay widths can be gathered from the Particle Data Book (PDG) [24]. With this in mind the probability of a specific channel corresponding to the inelastic cross section is expressed in Eq. 14.
where R i (E) is the cross section ratio from our statistical model, where all the possible fireball decay probabilities are summed over (Eq. 15).
where the sum of all the possible k = 1..N fireball schemes for a specific final state are taken. An example is the final state 3π + 3π − where one has to sum the 1-,2-, and 3 fireball contributions as well. If one wants to calculate inelastic cross sections the full normalization sum has to be calculated and possibly an uncertainty analysis has to be done, which is the main focus of the next section.

Normalization and uncertainty estimation
The normalized probability in Eq. 4 should give back the σ I /σ inelastic cross section ratio by definition. If we have a specific channel "i which only could come from a one fireball decay scheme, then only the one fireball normalization sum should be calculated. The following four processes were calculated to test the method: As the measured cross sections and the model also could contain errors it is possible to estimate an error distribution of the model from the distribution of the relative errors. To this purpose, the relative error is calculated as it is shown in Eq. 16.
where r 0 = σ i /σ inel and W 1 is the one fireball model probability. The model calculations does not contain A + B → C → D + E reactions, due to the lack of knowledge of the branching ratios of heavy mesons to proton and antiproton. Using the Gaussian error propagation formula, considering that Δσ i = 0, Δσ inel = 0, and ΔW 1 = 0 the absolute error of the measured relative error Δk is expressed in Eq.17.
The relative error of the model could be expressed from Eq. 17 as follows (Eq. 18): If the Δk uncertainty is approximated by the standard deviation of the relative error distribution, and assuming an energydependent model error, the relative model error distribution could be expressed from Eq. 18 and its histogram can be seen in Fig. 2: The estimated relative error of the model is calculated from this distribution by taking its mean value, which is approximately ΔW 1 /W 1 ≈ 0.4. This value will be the approximated relative error of the model, which will be used in the following sections, where the inelastic cross sections are calculated. The calculated model probability with the measured ratios can be seen in Fig. 3, where the uncertainty bound is also calculated from the previously estimated model error. It is again worth to mention, that in this calculations the possible p p → X → π + π − processes are neglected, where X is a heavy meson above the two proton threshold. This is mainly due to the lack of knowledge of the branching fractions at the PDG for heavy meson decays to proton and antiproton.  Fig. 3 The process pp → π + π − ratio to the inclusive cross section σ inel pp with the full normalization for one fireball. Data taken from [25,26] Fig. 4 The relative error dependence of the two fireball normalization ratios on the number of inelastic-, and inclusive channels (N,M).

Inclusive cross sections
In heavy-ion simulations, if the interesting quantities are particle multiplicities, then the desirable cross sections are the inclusive ones e.g. if one interested in the charmonium spectra which is measured by their dilepton pair decays an interesting background process could be X → DD, where each D meson could decay into an electron(positron) and a Kaon, giving an overall dilepton pair at the end. In this example, the (X → DD+ anything) process should be calculated as every possible channel is giving an extra contribution to the full background. The naive way to calculate inclusive cross sections to sum over all the possible reactions, which respects the quantum number conservation laws, however, this method is not very efficient and there is a much easier way to do this.
Let us take the probability ratio of a channel with two different normalization, which only means the normalization sum contains different channels. The inelastic sum contains all the possible final states allowed by quantum number conservation, while the inclusive sum contains every possibility, but with one specific particle fixed in every possible final state. To make things easy the reference channel should be strictly 1, 2, 3, ...N fireball channel without mixing so that the normalization is clear. Let us take a reference channel I , which is coming from strictly one fireball channel and calculate the ratio like in Eq. 19: where in the ratio the common factors e.g. P f b 1 are dropped out and only the normalization terms in the denominator remained. The Breit-Wigner factors are also considered in the phase-space integrals. In this way the inclusive cross section can be expressed as it is shown in Eq. 20: A problem with using the one-fireball ratio is the smallness of T i at higher energies, which means that the probabilities at higher energies for all of the one fireball possibilities are negligible or practically zero. A solution to this problem is to use the two-fireball ratio due to their non-negligible normalization sums even at the energies up to 15 GeV. The inclusive cross section from strictly two-fireball channels can be expressed in the same way with the corresponding normalization sums as in Eq. 21: where i j∈S inclusive means every possible two-fireball final state from the inclusive set, which respects the conservation laws of the initial state, and for simplicity, we defined H i (x) as: Typically inclusive cross sections like pp → I + X are needed where I is a fixed particle. The inclusive sum then includes all of the possible 2 fireball final states, which contains at least once particle I in its final state, however, it is also possible to calculate inclusive reactions like pp → I +J +X , where I and J are two fixed particles and X is everything else. The generalization is straightforward and the only task is to separate the possible 2 fireball reactions where we have a specific number of fixed particles. For one fireball scheme, the calculation is simple and it is possible to include every possible two and three-body final states into the sum with little effort, however, if one wants to calculate the ratio from more than one fireball decays the possible number of channels are too huge and the integration would take many hours on a standard laptop. Due to the numerical complexity, except for one case, only the one fireball calculations are shown in this paper, however, it can be shown that if one wants to lower the uncertainty of the model, then taking the two fireball normalization ratio is a good option to do so. To see this let us assume that each H i (x) is a uniformly distributed random variable with some mean μ > 0 and a relative width 0.5μ. Taking the two-fireball normalization ratio the following distribution has to be calculated (Eq. 23): where g is a uniformly distributed random variable, and M < N because M corresponds to the inclusive sum, and N corresponds to the inelastic sum. For simplicity let us assume that each g i is energy independent so that f (E) = f will also be energy independent. To further simplify the problem the ratio in Eq. 24 is calculated: where g i is now constant, so the integral can be factored out, and instead of taking every possible i, j combinations, we simply take the square of the generated g i random variables. This will not change the main conclusion, so for the first estimation, it is sufficient to see the relative error distribution of the calculated ratio in Eq. 24. The relative error distribution is calculated in the following steps: 1. Generate N number of uniformly distributed U [0, 1] → g 0 i samples. These will be the noiseless data. 2. Add a random noise to each of the g 0 i 's so that will be the data with a relative error of 0.5. 3. Take M number of samples from the previously generated samples. These will correspond to the inclusive sum. 4. Calculate the ratios with the noiseless samples f 0 , and with the noisy ones f . 5. Calculate the relative error r = | f 0 − f |/ f 0 Do the previous steps many times so at the end a relative error distribution from the calculated r 's is obtained. As has been done previously, the estimated relative error will be the mean of the obtained relative error distribution. The mean relative error dependence on the number of inelastic channels (N) i = 1..100, and the number of inclusive channels (M) j = i..100 can be seen in Fig. 4. The number of inelastic channels has to be larger than the number of inclusive ones e.g. N > M. The relative error at low multiplicity is around 0.5, which means the accumulated relative error is the same as the individual relative errors. If the multiplicity in the numerator and in the denominator increasing the accumulated relative error will be smaller. The main conclusion from this simple simulation is that if one wants to make the estimated relative error smaller it is sufficient to calculate e.g. the two fireball normalization sums instead of the one fireball sums, however, its numerical complexity is much greater. Also, the two fireball cases will cover a different energy range due to the different thresholds and different processes considered.

Non-charmed hadrons
To validate the theory six channels were calculated with the formalism showed in Sect. 3. In three process the one fireball probability ratios were calculated and compared to the measured cross section ratios, while in one case the two fireball ratio is calculated. The errors of the measured-and the calculated cross section ratios are also shown. The first reaction is the pπ − → ρ 0 X , where X is any 1, 2, ...N number of particles, which allowed by the quantum number conservation of the initial pπ − collision, and is allowed by the specific number of fireballs we used in the calculations. To calculate the ratio, first, the inelastic sum has to be done. The second step is to find all the resonances which could decay into ρ 0 e.g. N 1700 → ρ N and select only those channels into the sum (with the direct ρ 0 X channels as well). These are most of the Δ-s and nucleon resonances and are taken from the PDG, where only particles that have at least 3 stars are included. The results can be seen in Fig. 5.
In the second example, the process pπ − → K 0 X is calculated with one-and also with two fireball ratios in order to compare the results. In the inclusive sum, the particles which could decay into K 0 -s have to be included. These are some of the N , Λ, Σ, and K * resonances. The inelastic sum has to include all processes which have one baryon number, zero strangeness, zero charmness, and zero charge. The ratios from the two-and one fireball processes are shown in Fig. 6, where in both cases a really good agreement with the measured data is achieved. The error in the two fireball case is  The process pπ − → K 0 X ratio to the inclusive cross section σ inel pπ − from two-and from one fireball processes. Note the smaller uncertainty bound in the two fireball case. Data taken from [25,26] calculated from the simulation described in Sect. 4, using the number of inelastic, and the number of inclusive channels.
The third example is the channel pp → ρ 0 X , where the inelastic sum will be different than in the previous two reactions because in pp collisions all the important conserved quantum numbers are zero. The results are shown in Fig. 7, where again a really good agreement is achieved even at higher energies (15 GeV).
The fourth calculated process is the pp → ρ 0 X inclusive channel, where the inelastic sum includes all the possible final states with quantum numbers (baryon number, charge, strangeness, charmeness)=(2,2,0,0). The inclusive sum contains all the processes where at least one ρ 0 meson is present and also all the relevant resonances which could decay to ρ 0 . The results are shown in Fig. 8.  Fig. 7 The process pp → ρ 0 X ratio to the inclusive cross section σ inel pp . Data taken from [25,26] The final non-charmed examples are the strange vector meson K * (892) + and K * (892) − production in π − p collisions. The strange vector meson production is an important tool to study the dense hadronic matter in heavy ion collisions [27][28][29] due to their spectral function dependence on the nuclear density and temperature. We expect the model to give back the low energy suppression of the K * (892) − state to the K * (892) + meson. The results can be seen in Fig. 9, where the upper side shows the K * (892) + ratio to the π − p inelastic cross section and the lower side shows the K * (892) − one. Also some measured values are shown from [25,26]. It can be deduced that the main feature (K * (892) − suppression) can be described quite well with the model. It should be worth noted however that for the strange quarks, we used a constant value for the P s quark creational probability, fitted from pp → K + K − process at the energy range E ≈ [1.8, 2.6] GeV, which is a simplification, as the strange quark suppression should be energy dependent, so does P s . Nevertheless, even with this assumption the model gave back the main qualities of the cross sections.
In every aforementioned process, a really good match is achieved with only one fireball process considered. Using two fireball ratios the match is even better, which corresponds to the smaller model error estimated in the previous section.

Charmonium production
One of the models aim is to give estimates to charmonium creation cross sections below 10 GeV, which will be very useful inputs for upcoming transport simulations, regarding charmonium mass shifts in heavy ion collisions. To make reliable estimates the model should describe well the available low energy cross sections, which are very rare at the moment. A good collection of data are used for fitting in [30]. In our  The processes pπ − → K * (892) + X and pπ − → K * (892) − X ratio to the inclusive cross section σ inel pπ − . The upper half of the figure shows the K * (892) + production, while the the lower half shows the K * (892) − case. Data taken from [25,26] model, the necessary parameter, which has to be fitted from experiments is the charm quark creational probability (P c ), which is used to calculate the quark combinatorial factors. Due to the much larger mass of the charm quark than the u,d,s quarks, this value are expected to be much smaller then the others, thus to create more than one charm-anticharm pair is negligible in the calculations. Furthermore we expect it to increase with energy, rather than to be a constant value. In contrast with our previous attempt in [11], where we tried to fix a constant P c value at one point using a fit from [30] near threshold, with the calculation of one exclusive channel, now we have the method to fit for measured inclusive data at higher energies. This was one of the main reasons we extended the model to describe inclusive cross sections.
For the quark creation probability, we expect a linear rise with energy with the simplest assumption possible P c = a E, where only one parameter, the slope a has to be determined. This parameter enters into the quark number probability mass function described by a multinomial distribution shown in Eq. 25 where is the total number of quark-antiquark pairs [31], with hadronization temperature T 0 , and P i is the quark creational probability for the i = (u, d, s, c) type quarks. The expected number of quark-antiquark pairs of type i, will be n i = P i N , which corresponds to the maximum of the probability mass function. Note that the bottom quark is still missing, and we intend to fit the corresponding P b value as well in the future. The expected number of different quarks and antiquarks then build up the hadrons in all possible ways described by simple combinatorics, multiplied by the corresponding probability from the probability mass function, and at the end are normalized with all the possible final states, giving a final probability to create one specific hadronic final state. As the charm quark creation is expected to be highly suppressed only n c = 1 is interesting at the moment, which will corresponds to a constrained maximum value of F(N , n i ) in contrast to the non-constrained maximum, where n c = 0. The suppression ratio is then proportional to the charm-anticharm creational probability, so it can be easily fitted using the suppression to the non-charmed channels. To see this, let us assume that the one charm-anticharm quark pair is created in expense to one strange-antistrange quark pair, so n c = 1, and n s = n s − 1, while n u = n u and n d = n d . The suppression ratio calculated from the probability mass function is then: The derived formula for γ c expresses the fact, that even with a constant P c value, the suppression will be energy dependent if γ c < 1, as the total number of quark-antiquark pairs N are also energy dependent. After γ c reaches 1 the suppression from the probability mass function is gone, so no need to constrain the distribution anymore. The expected number of charm-anticharm quarks then will be n c = P c N , so the transition between the suppressed and non-suppressed probability is continous. The vanishing suppression ratio here only means, that there will be no more suppression from the constrained probability mass function, however if P c is still smaller than the other quark creational probabilities (which will be the case up to a few 100 GeV), then the charm production is still suppressed compared to the up, down and strange quarks.
As we assumed a linear relationship between P c and the energy, the γ c suppression ratio can be expressed as: where a is the slope parameter, which has to be determined from experiments. Using experimental charmonium production cross section data in pN collisions, this parameter is determined to be a = 0.0006 GeV −1 (see below), so the charm-anticharm quark creational probability is given by P c = 0.0006 · E. With the inclusion of another type of quark the previous P u , P d , P s fits were also have to be modified. In the simplest case, we can assume that P i = P i − P c /3 for i = (u, d, s), so that P u + P d + P s + P c = 1 holds, however due to the small value of P c in the energy range below a 100 GeV, we get P u ≈ P u , P d ≈ P d , P s ≈ P s , so the previously used values for the u, d, s quarks can still be used in this energy range. To calculate the cross sections the two fireball probabilities are used, as it gives smaller uncertainties, with the price of taking much longer time to calculate the normalization factors. We intend to describe the available low energy pN and π N data with the model for which, we selected some experimental points below 10 GeV from the collection in [30]. At this stage the uncertainties in the experiments are not really important as the main goal is to describe the 2-3 times larger cross sections of the π N collisions in contrast to the pN collisions. To this end, and to check the model predictability, only one measurement point is used to the fit from the pN data at E ≈ 10 GeV, which is enough, as we only have one free parameter in P c , namely a. The tuning was done by hand to approximately give the correct value at the one point used for the fit. For the cross sections the decays from Ψ (3686), χ c1 and χ c2 mesons are also included, however as there are no bottom quarks in the current model the decays from hadrons containing b quarks are not included. After the tuning of the slope parameter, the full cross sections were calculated for the pp and for the π − p case. The results can be seen in Fig. 10 and in Fig. 11. The calculations show a good match with the data for both cases even with the simple form of the P c and with the restricted one point fit, thus it can be concluded that the main feature of the π N enhancement is described well in the model. To test the model capabilities further the proton-antiproton to proton-proton inclusive J/Ψ cross section ratio is calculated at E = 24.3 GeV CM energy, which is measured in [32] with a good accuracy. The results with the uncertainties are shown in Table 1, where a remarkably good agreement with the measured value is achieved. Inclusive J/Ψ production cross section estimates in π − p collisions. Data taken from [30] The method thus proved to be able to give satisfactory results to low/medium energy charmonium cross sections, however putting more particles (e.g. B mesons) into the model could change the higher energy parts of the cross sections, meaning the P c quark creational probabilities might have to be further fine tuned. Also for P c a very simple linear assumption was made, and it could be interesting to give a theoretical background to its nature. This also applies to the strange quark creational probability P s , where we simply assumed a constant value throught the whole energy range we are interested in.

Conclusion
In this paper, a statistical model is introduced based on the statistical Bootstrap approach, which is able to give reasonable good estimates to hadronic cross section ratios up to a few GeV energy range. The full normalization sum is calculated for one and two fireball decay schemes. The cross section ratios for eight different processes were compared to the model calculations. From the calculated relative errors, the model uncertainty is also estimated giving an overall energy independent relative error for the one fireball decay probability. The calculated normalized probabilities are in good agreement with the measured values. The model is extended to give inclusive cross sections from the ratios of the normalization sums, which is validated through six distinct inclusive non-charmed processes, namely, pπ − → ρ 0 X , pπ − → K 0 X , pπ − → K * + X , pπ − → K * − X , pp → ρ 0 X and pp → ρ 0 X , where in each case a really good match with measured data is achieved. The charm quark creational probability is fitted with the help of the pN → J/Ψ X inclusive data, allowing us to give estimates to charmed final states as well. The model is further validated through two charmed processes, namely pp → J/Ψ X and π − p → J/Ψ X , both giving good results compared to the experimental data. For a last application the ( pp → J/Ψ X ) / (p p → J/Ψ X ) ratio is calculated at √ s = 24.3 GeV giving a really good match with the experimental data. The method will be highly useful in describing previously unknown or not well measured inclusive cross sections e.g. charmonium creational probabilities and will be used in a BUU type [1,33,34] off-shell [35] transport code to study the propagation of these states in dense nuclear matter. included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.