A Monte-Carlo generator for statistical hadronization in high energy e+e- collisions

We present a Monte-Carlo implementation of the Statistical Hadronization Model in e+e- collisions. The physical scheme is based on the statistical hadronization of massive clusters produced by the event generator Herwig within the microcanonical ensemble. We present a preliminary comparison of several observables with measurements in e+e- collisions at the Z peak. Although a fine tuning of the model parameters is not carried out, a general good agreement between its predictions and data is found.


Introduction
The Statistical Hadronization Model (SHM) has proved to be successful in reproducing essential features of the hadronization process, such as the multiplicities and the transverse momentum spectra of many hadronic species in high energy elementary and nuclear collisions (see e.g. [1] and references therein). Recently, it has been shown that this model is also able to satisfactorily reproduce the rates of multi-hadronic exclusive channels in e + e − annihilation at low energy [2]. According to the SHM, hadronization proceeds through the formation of extended massive colorless objects called clusters or fireballs emitting hadrons according to a pure statistical law, i.e. all multihadronic states compatible with conservation laws are equally likely. The a e-mail: christopher.bignamini@pv.infn.it b e-mail: becattini@fi.infn.it c e-mail: fulvio.piccinini@pv.infn.it number and properties of these clusters are determined by the pre-hadronization dynamical process and cannot be predicted within the SHM. Therefore, in e + e − and pp (−) collisions at high energy ( √ s > 10 GeV), the SHM calculations have been carried out with supplementary assumptions, in order to obtain analytical or semi-analytical formulae to be compared directly with the data. Such is the case of the so-called Equivalent Global Cluster picture [1,3], in which a further assumption is invoked concerning the distribution of charges and masses among the hadronizing pre-hadronic clusters; this allows to make calculations in the more manageable canonical ensemble.
While these supplementary assumptions have been very useful to point out the statistical features of hadronization, it would be desirable to perform tests of the SHM independently thereof, at the most fundamental level of formulation. As clusters are generally small, calculations in a basic SHM framework involve the computation of averages in the microcanonical ensemble, i.e. with exact conservation of energy and momentum, besides that of abelian charges (baryon number, strangeness, electric charge, charm and beauty). This kind of calculation has been recently carried out (with the further complication of angular momentum conservation and more) for e + e − annihilation at low energy [2]. Yet, in this case, only one cluster was involved instead of many as it is generally the case for high energy collisions.
In a multi-cluster environment, the only possible way to cope with the problem of statistical hadronization is a Monte-Carlo integration. In this paper a prototype of Monte-Carlo event generator for the SHM is presented and the results obtained by testing it on e + e − collisions at √ s = 91.2 GeV center of mass energy, i.e. at the Z peak, are discussed.

arXiv:1204.2300v2 [hep-ph] 17 Oct 2012
For the present study, the developed Monte-Carlo hadronization code has been interfaced to the Herwig6510 event generator [4], which is used to simulate the e + e − collisions except for the hadronization phase. This step is instead performed using the SHM hadronization module. This particular choice for the external generator is owing to the similarity between the input clusters required by the SHM and those produced by Herwig6510 [5]. In our framework these clusters are used as starting point for the statistical hadronization process. Finally, the unstable hadron decays are performed with the Herwig6510 generator itself.
The paper is organized as follows: in Sec. 2 the general framework of the SHM is summarized and the algorithms for the simulation of single cluster hadronization are presented. In Sec. 3 a comparison between the theoretical formulation of the SHM and those of the Cluster Model implemented in other Monte Carlo generators, such as Herwig6510, is reported. The computational methods for multi-cluster hadronization in the SHM framework are then discussed in Sec. 4. In Sec. 5 the full event simulation setup, involving the external event generator, is described, as well as results obtained for a preliminary adjustment of free parameters. The impact on SHM predictions of variations of some Herwig6510 free parameters involved in the cluster formation procedure is also discussed. In Sec. 6 the results obtained for e + e − collisions at the Z peak with the SHM-based hadronization module are presented and compared with the corresponding Herwig6510 predictions and LEP experimental data. Conclusions and future developments are outlined in Sec. 7.

Statistical hadronization of single clusters
As it has been mentioned in the Introduction, the cornerstones of the SHM are: in a high energy collision, as a result of the prehadronic dynamical process, a set of extended massive colorless objects called clusters or fireballs are formed; all multi-hadronic states confined within the cluster and compatible with conservation laws are equally likely.
The second postulate gives rise to a definite mathematical formula [1,2] involving (pseudo)projector operators associated to conservation laws, for the probability of observing a specific final state from a single cluster decay. In principle, conservation laws to be implemented include energy-momentum, angular momentum, parity and all of the internal symmetries relevant to strong interactions. However, like all other Monte-Carlo generators, we confine ourselves to abelian charges (electric charge, strangeness, baryonic number, charm and beauty), in addition to energy-momentum. In fact, in a multi-cluster environment, relevant to this work, dealing with non-abelian symmetry groups is overwhelmingly difficult if clusters take their origin from a pure quantum state and are entangled with each other at the hadronization stage.
In the simpler scheme of additive abelian charges and energy-momentum conservation, the probability of a multi-hadronic final state for a single cluster can be derived as described in [6,7]. Particularly, neglecting quantum statistics and quantum field effects, the microcanonical weight Ω {Nj } for the multi-hadronic channel {N j } made of N 1 , N 2 ,..., N K hadrons of species 1, . . . , K is given by [6]: is the array collecting its abelian charge values (baryon number, charge, strangeness, etc.), and N is the total number of particles in the channel. P , Q and V are the cluster 4-momentum, abelian charges and volume respectively. Finally, it must be noted that a sum over particle spin states has been performed: the J j quantities in the corresponding contributions to the total weight are the spin modulus of the single hadron species. The microcanonical probability p({N j }) of the considered channel is given by: where Ω is defined as the sum over all possible channels i.e. the microcanonical partition function. Finally, the microcanonical probability density in phase space can be written as: where {p j } = p 1 , p 2 , ..., p N is the set of four-momenta of the particles in the channel {N j }. The analysis of particle multiplicities has shown [1] that particles carrying strange valence quarks need an additional suppression with respect to the pure statistical predictions. This has led to the introduction of the strangeness suppression parameter γ S , which, in order to fit its definition for inclusive multiplicities, has to multiply Ω {Nj } of Eq. 1 as follows: where N j is the number of hadrons of species j in the channel {N j } and s j is the total number of valence strange and antistrange quarks contained in the j-th hadron species. Extra strangeness suppression also applies to light unflavored mesons, like the η meson, which are actually a flavorless superposition of states with a wave function of the general form: with |C u | 2 + |C d | 2 + |C s | 2 = 1. For this kind of hadrons only the ss part of the wave function is supposed to be suppressed, with a suppression factor given by The γ S parameter and the cluster energy density ρ are the only free parameters of the model, in its microcanonical formulation. The cluster proper energy density ρ plays the role of conversion factor between the cluster volume V , present in Eqs. 1-4, and the cluster mass M through the relationship

Single cluster hadronization algorithm
The simulation of the cluster microcanonical hadronization is mainly based on the algorithms described in [8] and here summarized. More in detail, the single cluster hadronization process involves two main steps, namely the sampling of the hadronization channel and of the corresponding kinematical configuration. Due to the large number of available decay channels for a hadronizing cluster, an efficient sampling algorithm is needed to obtain a fast and optimized hadronization simulation: for the present case the sampling function Π has been defined, as discussed in [8], as the multi-species multiplicity distribution of the SHM grandcanonical formulation where ν j is the mean number of particles of the j-th type, K the number of included hadronic species and N j the number of particles of kind j contained in the channel {N j }. The mean multiplicities ν j are free parameters which should be set with the aim of obtaining the most efficient sampling function. The abelian charge conservation is then introduced in the channel sampling procedure based on Eq. 9 in the following way: the whole set of hadrons is divided into eleven sampling groups, namely, (anti)bottomed hadrons, (anti)charmed hadrons, light (anti)baryons, light strange (anti)mesons, light charged (anti)mesons with zero strangeness and neutral light mesons. Moreover, the following feature of the multi-poissonian function is used: considering for simplicity only the baryon, antibaryon and meson hadron groups and the relative charge conservations, the original sampling function with an extra Kronecker delta, added to impose the conservation of abelian charges, can be written as where is the poissonian distribution of the total number of baryons (x = b) or antibaryons (x =b) and where the corresponding mean multiplicity ν x is given by the sum of the single baryon or antibaryon mean multiplicities. The functions P are the conditioned multinomial distributions of the single hadronic species, namely with, again, x = b for baryons and x =b for antibaryons.
With the above decomposition of the original multipoissonian distribution, the following sampling algorithm can be used to perform the hadronization channel generation: 1. extract randomly the number of baryons and antibaryons according to the distributions Π b (N b ) and Πb (Nb); 2. check whether the baryonic number is conserved. If not reject the sampling and go to point 1, otherwise generate the single baryons and antibaryons using the multinomial distribution of Eq. 12; 3. extract the single mesons using the initial multipoisson distributions; 4. check the conservation of the remaining abelian charges, taking into account the already generated baryons and antibaryons. If the check is not passed start again with point 1, otherwise the sampled channel is momentarily accepted and the corresponding phase space availability can be verified.
The adoption of the described sampling procedure, instead of the independent sampling of the single hadronic species, which would follow from Eq. 9, is motivated by the reduction of random number extractions which occur in case of rejection of the sampled channel due to the abelian charge conservation: the independent sampling procedure would require, for each channel sampling attempt, the extraction of K random numbers, where K is the number of hadron species included in the hadronization procedure, which would be lost in case of event rejection. The above algorithm, on the other hand, allows to stop and start again the sampling procedure, in case of charge conservation failure, when a number of random extractions smaller than K has been performed by checking the conservation of the single abelian charges during the sampling.
The only exception to the above sampling procedure is represented by the heavy flavored hadron sampling step, required when a heavy flavored cluster must be hadronized. For these objects, which are supposed to hadronize into channels containing one heavy hadron and a set of light particles, the heavy hadron is sampled with probability P h of the form with a > 0 and where m h is the hadron mass. This particular choice for the heavy flavored hadron sampling function is justified by two reasons: on one side it allows to avoid the strong undersampling which would follow from the usage of the standard sampling method (Eq. 9) also for these hadronic species. On the other side, the above function allows to approximatively take into account the heavy hadron relative production probabilities due to their mass differences. The remaining part of the channel is then randomly chosen using the same algorithm applied for the hadronization of light clusters. This modification in the sampling function for the heavy flavored clusters allows to obtain a more efficient channel generation algorithm with respect to the standard procedure previously described. Given a hadronization channel composition, a kinematical configuration must be generated: in the present work, a sampling algorithm inspired by the multi-particle decay phase space integration described in [9] 1 is used. The starting point to obtain the kinematical configuration and the corresponding weight, for a decay channel of a cluster of mass M , containing N particles with masses m 1 , m 2 , ..., m N , is the phase space integral of Eq. 1: As discussed in [9], the following condition holds for the relativistic N -body phase space element: and The non-relativistic phase space integral of Eq. 14 can be rewritten making use of the relativistic phase space element of Eq. 16 as follows: The integration/sampling procedure adopted here is based on the property of Eq. 15, iteratively applied in order to deal with 2-body decays only as described by with q 0 = (M, 0) and where the integration limits of the q 2 i virtualities are given by By solving the 2-body phase space integrals, Eq. 18 becomes where |p i | is the i-th particle momentum modulus, determined by energy-momentum conservation, in the ith 2-body decay rest frame, namely for i = 1, ..., N − 2, and for i = N − 1 and i = N . Moreover, the energy componentsp i (0) andq i (0) are given bȳ and E i is the i-th particle energy in the cluster rest frame corresponding to the momentump i . It follows from Eq. 20 that the needed N -body phase space configuration can be obtained by sampling N -1 solid angles, corresponding to the particle emission direction in the 2-body decays, and N -2 virtualities. In the present work, the solid angles Ωp i and the q 2 i virtualities are randomly generated within 4π and within the limits defined in Eq. 19 respectively. The particle momenta generated with this procedure, initially defined in the 2-body decay rest frames, are subsequently boosted to the laboratory frame. Finally, the phase space sampling weight w P S is given by Defining the hadronization channel sampling weight w HC as for light clusters and as for heavy flavored clusters, where the hadronization channel sampling functions Π {Nj } and P h are defined in Eqs. 9 and 13 respectively, the total weight w corresponding to a single cluster hadronization event reads

Comparison with the Cluster Model
Although the SHM is based on the general idea of cluster hadronization, it has peculiar differences with respect to other hadronic event generators implementing Cluster Model hadronization such as Herwig++ [11] and Sherpa [12], besides Herwig6510. Let us now discuss the differences between these approaches more in detail starting from the Herwig case.
The most striking difference is concerned with the generation of the channel multiplicity: whilst the Cluster Model forces clusters to decay into hadron pairs, in SHM the final multiplicity is an outcome of the cluster hadronization process. Particularly, the relative production rate of N -body channels is determined by the cluster's finite volume. This happens because the SHM uses the proper measure of the multihadronic phase space [1,13] with a cluster's proper volume V acting as a coupling constant for the cluster decay according to Eq. 1. Conversely, the Cluster Model uses invariant momentum space measure d 3 p/2E without extra factors and therefore the final multiplicity is to be fixed otherwise (see extensive discussion in [13]).
Indeed, in the Herwig event generator, a hadronizing cluster can be transformed in three different ways depending on its mass: splitting into a lighter cluster pair (see discussion in Sec. 5.2); decay into a hadron pair; transformation into the lightest hadron with the same flavor composition of the cluster.
Here we point out a further difference, even for the twobody decays: given a cluster with flavor composition f 1f2 , in the Herwig Cluster Model the decay hadron pair will have the flavor composition f 1f and ff 2 respectively, with the hadrons sharing the components of the ff parton pair extracted from the vacuum to perform the cluster decay in this framework. On the other hand, in the SHM only abelian charges conservation is required, implying the inclusion of a larger set of channels.
With respect to the general aforementioned Herwig hadronization scheme, Sherpa implementation of the Cluster Model features additional differences with respect to SHM. Indeed, relevant changes have been introduced in implementing the Cluster Model in the event generator [14], both for cluster production and decay: inclusion of soft color reconnection during cluster formation and decay; inclusion of the parton/hadron spin degree of freedom in the transition probability calculation; a refinement of the single hadron cluster transformation. Particularly, in the Sherpa case all clusters with mass falling into the so-called hadronic regime are transformed into single hadrons with the same flavor compositions of the hadronizing clusters. Each specific hadron is then randomly chosen among those with a mass smaller than the one of the cluster itself, instead of always choosing the lightest one as in the Herwig case.

Multi-cluster hadronization and event weight
Within cluster hadronization models, high energy e + e − collisions (say √ s ≥ 4 GeV) always involve multiple cluster production. For instance, the Cluster Model [5] implemented in Herwig6510 event generator [4], provides the formation of color connected parton-antiparton pairs at the end of the perturbative QCD shower, giving rise to massive colorless clusters. As discussed in Sec. 3, those clusters are thereafter decayed into pairs of massive hadrons, the probability of the specific decay channel being determined by the two-body phase space availability and by the spin multiplicities of the involved hadrons. Indeed, in the present work, the multi-cluster system formed as a result of the hard scattering process, followed by the perturbative showering, is taken from the Herwig6510 event generator and used as input for the SHM. Specifically, the general simulation setup is as follows: the Herwig6510 event generator is used to simulate the full pre-hadronization process, from initial state radiation emission to cluster production. The Herwig6510 preliminary cluster splitting procedure described in [4], which follows the cluster formation step, is also included in the event simulation (the details of the cluster splitting procedure and its impact on the SHM predictions are discussed in Sec. 5.2). At that stage, the standard hadronization algorithm is replaced by the SHM one, which takes care of hadronizing all produced clusters; finally, the formed hadrons are returned to the external generator (i.e. Herwig6510) which takes care of decaying of the unstable particles.
The above hadronization procedure is applied to each cluster independently, thus, for each simulated collision, the hadronization event weight W reads: where N C is the number of clusters in the event and w i the i-th cluster hadronization weight defined in Eq. 27. Moreover, w i HC and w i P S are the hadronization channel and phase space sampling weights of the i-th cluster (defined in Eqs. 24, 25 and 26) and Ω i its partition function (defined in Eq. 3).
This method in fact involves a complication, related to the calculation of the hadronization event weight of Eq. 28. Specifically, it requires the prior knowledge of each cluster's microcanonical partition function. If we had only single cluster event with fixed mass and charges, its microcanonical partition function would be an irrelevant constant factor cancelling out when calculating mean values. On the other hand, for a multi-cluster environment, the value of the product of microcanonical partition functions in Eq. 28 is not constant, in fact it is a function of the masses and the charges of the particular set of clusters. Therefore, we need to know the specific numerical value of the microcanonical partition functions of the clusters produced in a single event in order to correctly normalize its weight.
Unfortunately, computing the microcanonical partition function of a cluster of given mass and set of charges is a CPU time demanding task that cannot be afforded at the event generation time. To solve this problem, we decided to pre-calculate microcanonical partition functions for a discrete set of masses (and charges) and store their values in a look-up table to be used at event generation time. The number of pre-calculated functions has been determined on the basis of the observation that, in about 10 6 e + e − collision events at the Z peak, around 150 different cluster charge configurations (including baryon number, strangeness, electric charge, charm and beauty) occurred. For each of these charge configurations, a set of microcanonical partition functions has been calculated, each set including different cluster mass and different free parameters (γ S and ρ) values. Specifically, for each charge configuration, a grid in the mass-ρ-γ S space has been defined as follows: and for each grid point an Ω calculated. Then, during event generation, the partition function values of the single clusters, for their specific charge configuration/mass and for the chosen SHM free parameter values, are determined by means of a linear interpolation between the nearest points in the grid.
As it can be seen, a non uniform grid structure along the cluster mass direction has been adopted: this choice was motivated by the non-smooth dependence of the microcanonical partition function on the mass at low mass values. Indeed, 3 GeV as mass value separating the two differently discretized regions proved to be a good trade-off between the need of good accuracy and the computational cost in getting all grid points calculated.
A final remark concerning the calculation of the partition function set is in order: for a given mass, charge configuration and γ S and ρ parameters values, the corresponding microcanonical partition function value is determined by the set of hadrons used in the hadronization process and by their physical properties. It follows that the set of partition functions needs to be rebuilt only when changes are introduced in the hadron set, such as modifications to the hadron physical properties or to the hadron list.
The details of the numerical methods used for the pre-calculation of microcanonical partition function grid can be found in [10]. It is worth mentioning here that we have used the phase space integration optimized algorithm described in [8,15].

Preliminary tests
In this section the results obtained with the SHM code, interfaced to Herwig6510, for e + e − collisions at 91.2 GeV center of mass energy, will be presented and discussed. While a rigorous tuning of the SHM free parameters is needed to assess the goodness of the model, the results presented here, obtained with an approximate parameter adjustment, already show a good accuracy level. Moreover, in view of a global tuning of both SHM and Herwig6510 free parameters, a preliminary analysis of the dependency of SHM predictions on some Herwig6510 parameters involved in the cluster formation step is presented.

Simulation setup
As already discussed, the Herwig6510 event generator has been used for the present work as external code for the collision simulation, replacing its standard hadronization routines with the SHM Monte-Carlo module. As input clusters for the microcanonical hadronization, the clusters normally produced by Herwig6510 have been used, except for two modifications described below. The set of hadrons produced during each microcanonical hadronization process is then processed by Herwig6510, which performs the unstable particle decays.
The modifications introduced in the standard cluster sets produced by Herwig6510 during the event simulation, and above mentioned, are the following: -Cluster merging: the predictions of the SHM show a dependency on the input cluster mass spectrum, which in its turn depends on the Herwig6510 setup.
With the final aim of obtaining a better agreement between the SHM predictions and the experimental data, it is useful to have the possibility to modify the properties of the cluster mass spectrum. In the present work, the above modification is introduced by means of a cluster merging procedure, which, working after the Herwig6510 cluster production step, allows to control their mass spectrum without any changes in the external generator standard setup. Starting from a set of Herwig6510 clusters, the result of the merging procedure is determined by the value of a free parameter (M C ) representing the minimum allowed cluster mass: the set of incoming clusters is analyzed in order to check the presence of objects with mass under the given threshold. If one or more light clusters are found, an iterative fusion process is activated, which merges the cluster pairs containing at least one light object into heavier clusters, repeating this procedure until no light clusters remain. When more than one combination for a light cluster is possible the merged pair is the one with the smallest invariant mass where p i (p j ) is the 4-momentum of the i-th(j-th) cluster. It must be noted that a maximum mass value of 10 GeV is allowed for the new clusters produced during the merging procedure, a condition required to control the broadening of the cluster mass spectrum towards large mass values and the consequent quick increase in the computational cost of the corresponding partition function grid calculation. Because of the above condition, when a cluster pair is chosen for the merging procedure, the corresponding new cluster is created only if the pair invariant mass (Eq. 29) is smaller than 10 GeV: in a limited number of cases this condition can prevent the cluster merging procedure from being fully completed, with the result that clusters with mass below the chosen M C value can still be included in the hadronization procedure. -Baryonic clusters: the second change to the standard Herwig6510 cluster properties is related to the production of baryonic clusters, disabled in a standard run of this generator, and is obtained by means of a proper setting of its QDIQK parameter value. This parameter represents the maximum scale at which gluons can be (non-perturbatively) split into diquarks during the QCD shower process and is set to zero in the default Herwig6510 configuration, thus switching off the possibility of diquarks production. For the present work the value 2m c = 2 × 1.55 GeV has been chosen for the QDIQK parameter, with the result of activating the possibility of gluon splitting into light flavored diquark-antidiquark pairs. During the subsequent Herwig6510 cluster building step, each produced diquark (antidiquark) is coupled to the corresponding color connected quark (antiquark) to form a colorless baryonic cluster. The need for baryonic clusters is strictly related to the correct baryon production during the microcanonical hadronization process: because of the baryonic charge conservation, baryons can be obtained from the standard Herwig6510 non-baryonic clusters only as baryon-antibaryon pairs. However, these configurations are strongly suppressed because of phase space availability reasons. The introduction of baryonic clusters, on the other side, allows to obtain the production of baryonic final states in less restrictive phase space conditions.
In Fig. 1 the effects of the discussed cluster merging procedure on the cluster invariant mass spectrum are shown, in particular for clusters produced in 10 6 e + e − → dd collisions at 91.2 GeV. Primary clusters are the standard Herwig6510 clusters, except for the activation of baryonic cluster production, and their mass spectrum is compared with the one of clusters obtained with the merging procedure for M C = 1.6 GeV: the presence of the cluster minimum mass selection is evident, as well as of a broadening of the mass distribution. At the same time, it can be seen that a residual fraction of clusters, about 1%, have a mass below the M C value: the existence of these objects after the merging procedure is due to the maximum allowed cluster mass value of the merging procedure which, as already discussed, in a limited number of cases prevents the cluster merging from being completed. The mass distribution of the secondary clusters, namely the ones provided in output by the merging procedure, shows also an inflection point around 3.5 GeV. As it can be seen in Fig. 1, this mass value is the maximum allowed for primary clusters, a condition responsible for the observed inflection point: in fact, only new clusters produced by the merging procedure contributes to the mass distribution for values larger than the above limit, while for smaller mass values contributions from the primary cluster mass distribution are also present. These differences in the secondary cluster mass distribution composition, in the two mass ranges separated by the primary cluster maximum mass value, is the origin of the behavior change in the mass distribution.
Finally, we point out that in the hadronization step, among the light flavored hadron states, only those with mass less or equal to 1.8 GeV were included. This choice was motivated by the opportunity of comparing multiplicities with previous results obtained with the SHM where this cutoff was used. It should be noted, anyhow, that the production of light flavored states with mass larger than 1.8 GeV is negligible for essentially all observables.

Study of Herwig cluster splitting
As already stated, the Herwig6510 cluster splitting procedure described in [4] is included in the present event simulation setup. The splitting procedure, which in Her-wig6510 precedes the cluster decay step, operates as follows: a cluster of mass M f1f2 , composed of the partons f 1 andf 2 of mass m f1 and mf 2 respectively, is split if the following condition holds: where CLPOW and CLMAX are phenomenological parameters. For each cluster satisfying the above condition a flavor pair ff is picked from the vacuum and its components used for the construction of the two daughter clusters, whose flavor compositions will be f 1f and ff 2 respectively. The masses M f1f and M ff2 of the newly created clusters are randomly generated according to the equations where m f is the mass of the generated parton pair components, r is an uniformly distributed random number and where PSPLT is a third free parameter of the splitting procedure (in fact, PSPLT is a two component parameter: one component refers to light flavored and charmed clusters and the other to the bottomed ones). While a detailed analysis of the interplay between SHM free parameters and those belonging to Herwig6510 initial collision steps is needed in view of a global tuning, especially concerning the cluster formation and splitting free parameters such as CLPOW, CLMAX and PSPLT, a preliminary evaluation of their impact on the SHM  The obtained results show that only small effects are introduced by the considered variations in the CLMAX and CLPOW parameter values, as can be seen e.g. for the rapidity distribution reported in Fig. 2 and 3 (for the definition of the observables, see Appendix). On the other side a considerable dependency on the PSPLT parameter value is observed in the SHM predictions, as clearly visible in Fig. 4 again for event rapidity. The same situation is present for single particle observables, as reported for example in Figs. 5 -7 for ρ 0 meson scaled energy distribution, and for particle multiplicities (Tabs. 1 -3). Concerning multiplicities, it is worth noting the large variations in the total number of charged particles and in light meson multiplicities due to the modification of the PSPLT parameter, while CLMAX and CLPOW variation effects are essentially limited to the baryonic sector.

Adjustment of the free parameters
The free parameters of the SHM code, namely ρ and γ S , and the cluster minimum mass parameter M C introduced by the merging procedure, need to be rigorously tuned to correctly evaluate the performances of the code in reproducing the experimental data. For the present work only preliminary tests have been performed, using different parameter configurations, in order to obtain a global overview of the SHM predictivity level in a full collision simulation framework. Moreover, an approximate adjustment of the SHM code free parameters has been realized, through a comparison between the SHM predictions and the corresponding experimental data: this comparison has been performed considering a set of inclusive and exclusive observable distributions and mean values. In particular, event shape and single particle momentum/energy distributions, single hadron  10.09 ± 0.08 9.26 ± 0.07 11.6 ± 0.5 π + 8.71 ± 0.16 7.98 ± 0.10 9.9 ± 0.5 η 1.24 ± 0.03 1.14 ± 0.03 1.5 ± 0.2 ∆ ++ 0.081 ± 0.002 0.067 ± 0.001 0.075 ± 0.003 Σ + 0.0277 ± 0.0006 0.0205 ± 0.0004 0.029 ± 0.001 Σ − 0.0243 ± 0.0005 0.0203 ± 0.0004 0.069 ± 0.003 Σ 0 0.0357 ± 0.0008 0.0247 ± 0.0005 0.027 ± 0.001 mean multiplicities and the charged particle number distribution have been included in the analysis. The reference experimental data used for the comparison come from the measurements of LEP experiments for a center of mass energy of 91.2 GeV. The best parameter estimation has been realized by means of an approximate χ 2 minimization procedure: a set of explorative runs with various free parameter configurations has been considered, evaluating for each configuration the discrepancy between theoretical predictions and experimental data in the form of a global χ 2 value. As usual, the parameter best configuration has been identified as the one giving the lowest global χ 2 value. More in detail, for histograms the χ 2 value has been computed as the sum over channels of the discrepancy between the theoretical prediction y i t and the corresponding experimental data y i e , normalized using the theoretical (Monte-Carlo) and experimental errors σ i t and σ i e : where i is the histogram channel index and N the number of channels. The χ 2 value corresponding to the particle multiplicity analysis has been computed in a similar way, summing over the list of considered particles. Finally, the global χ 2 value has been computed as the sum over the single distribution/mean value χ 2 contributions. The performed analysis, whose details are reported in [10], shows that the parameter setup corresponding to the best global agreement between SHM predictions and the corresponding experimental data is:

Numerical results
In this Section, we present the results obtained with the SHM code for the above parameter set. The SHM predictions are compared to the corresponding predictions of Herwig6510 in its release version and to experimental data. As already mentioned, the reported results refer to e + e − collisions at the Z peak. It must be noted that the chosen observables provide a reliable measure for the evaluation of the hadronization model predictivity, because of their sensitivity to the hadronization process. This condition is clear for what concerns exclusive observables involving a single hadron species, such as the π ± moment distribution and mean multiplicity. Nevertheless, the same property holds also for the inclusive event shape observables referring to the event transverse plane, such as transverse momentum: the main contribution to the global hadron distribution in the event phase space comes from the hard scattering, whose quark emission direction approximatively defines the event thrust axis. On the other side, the hadronic phase space configuration projection on the event transverse plane is strongly influenced by the hadronization process. Therefore, the analysis of the transverse plane related observables allows to effectively evaluate the prediction accuracy of the considered hadronization model. The results obtained for a subset of these inclusive and exclusive observables are reported hereafter (for the definition of observables, see Appendix).
-Event shape observables: the comparison between the SHM predictions and LEP experimental data shows a quite good global agreement for the whole set of considered observables (Figs. 8 -11). However the SHM, for the present parameter configuration, seems to fail in reproducing the transverse momentum distribution tail behavior of the experimental data (Figs. 8 -9), which is instead correctly predicted by Herwig6510.
-Single particle observables 2 (Figs. 12 -27): also in this case a good agreement between the SHM predictions and LEP data is present for a large set of the considered single particle observables. Nevertheless, a wrong behavior in the SHM predictions can be seen in the D 0 and D * scaled energy distributions (Figs. 26 -27): in these cases also Herwig6510's predictions show a disagreement with the experimental data, however the Cluster Model seems to be able to reproduce the general behavior of these data better than the SHM. A more detailed analysis regarding these last results is reported hereafter.
-Charmed hadrons: the performed analysis, involving a complete collision simulation, has given the possibility to identify some possible limits of the SHM, not observed in previous studies. In particular, this is the case of the D 0 and D * scaled energy distributions of Figs. 26 and 27 respectively, where the SHM fails in reproducing the experimental distribution shapes. In these distributions two distinguishable regions are present: a low energy region (x E 0.5), whose filling events correspond to the D 0 and D * meson production in b-hadron decay, and a high energy region approximatively corresponding to the scaled energy distribution of the primary D 0 and D * mesons. While the distribution behavior in these two regions is qualitatively well reproduced by the Cluster Model predictions, the SHM in the adopted simulation framework fails in reproducing the high energy region. Even though this anomalous result needs to be further investigated, it is understood that the underestimation of the charmed mesons scaled energy, provided by the SHM, is strictly related to the large mean number of particles produced during the microcanonical hadronization of charmed clusters and to the low energy availability condition which follows. A condition not present in the two body cluster decays performed by the Herwig6510 Cluster Model.
-Charged particle multiplicity: in this case (Fig. 28) a very good agreement between the theoretical predictions of the SHM and the experimental data is present. It is worth noting that for this distribution, in particular for the distribution tails, the predictions of the SHM show a better agreement with the experimental data with respect to Herwig6510's results.
-Particle multiplicities: the comparison on the mean particle multiplicities (Tabs. 4 -6), for this preli-minary analysis, shows a global agreement of the microcanonical predictions with respect to the experimental data within 3σ in the 75.5% of the considered cases. This percentage becomes 100% if only charmed and bottomed hadrons are considered, for which a discrepancy within 1σ is present in the 55.6% of the cases. The corresponding Herwig6510 predictions agree with experimental data within 3σ in the 55.1% of cases for the whole set of particles and in the 72.2% of cases if only heavy flavored hadrons are considered. More in detail, a different behavior between the SHM and the Cluster Model predictions can be observed for the light flavored baryon multiplicities, with the SHM tending to underestimate the baryon yields and the Herwig6510 hadronization model showing the opposite behavior. Finally, a 5σ discrepancy between the SHM prediction of charged particle total number and the data can be noted. A better agreement for this quantity can be obtained with a more refined tuning of the parameters, probably at the cost of a worsening of the agreement between data and predictions of other observables (e.g. event shape distributions and single particle multiplicities). This discrepancy will be further investigated in the next SHM parameter tuning.

Conclusions and outlook
We have presented a Monte-Carlo implementation of Statistical Hadronization Model based on the microcanonical hadronization of massive clusters. The developed hadronization code has been interfaced to the Herwig event generator providing the massive clusters to be hadronized. With this method, a full description of an e + e − collision at high energy can be obtained with statistical model based hadronization. The results obtained in e + e − collisions at the Z peak have been presented and discussed. While a final and rigorous tuning of the hadronization module free parameters has not been performed in this work, the preliminary comparison confirms the good agreement of the SHM predictions with the data, already observed in previous studies. Specifically, the agreement between measured and predicted abundances of hadronic species is confirmed, particularly in the heavy flavor sector. A good agreement, between SHM predictions and experimental data, is also found for a set of hadronization relevant event shape variables and single particle energy/momentum distributions, with the only exception of the charmed meson scaled energy, where a consistent discrepancy is observed. Indeed, the SHM seems to fail in correctly predicting the momentum spectrum at high momentum. While this anomaly needs to be further investigated, it is already clear how its origin lies in the larger mean number of particles produced by statistical microcanonical hadronization of charmed clusters as compared with standard Herwig procedure.
These problems will be further investigated, and a global fine tuning of the free parameters performed. We envisage an extension of this hadronization Monte-Carlo code to high energy pp and pp collisions, with the final goal of achieving a global assessment of the Statistical Hadronization Model in elementary high energy collisions. A first public release of this code with the needed interfaces for its usage with Herwig6510 is forthcoming. 1/N dN/dp out LEP data Herwig6510 SHM Fig. 9 Transverse momentum (out) with respect to thrust axis normalized distribution: comparison among SHM predictions, Herwig6510 predictions and DELPHI data [16].  Fig. 10 Rapidity with respect to thrust axis normalized distribution: comparison among SHM predictions, Herwig6510 predictions and DELPHI data [16].          Fig. 19 π 0 scaled energy distribution: comparison among SHM predictions, Herwig6510 predictions and OPAL data [18].   Appendix: observable definitions -The thrust T is defined as: where n is a unit vector along the thrust axis, N is the number of particles and p i is the i-th particle 3-momentum. Thrust major M and minor m are similarly defined, replacing n with n M (perpendicular to n) and with n m = n M × n respectively.
-The rapidity with respect to the thrust axis y T is defined as: where E and p T are the particle energy and 3-momentum projection along the trust axis respectively. -The in and out components of the transverse momentum with respect to thrust axis, p in T and p out T , are defined as: where p is the particle 3-momentum and n M (n m ) is the thrust major (minor) axis previously defined. -The scaled energy and momentum x E , x p and ξ p are defined as: where E and p are the particle energy and 3-momentum respectively and √ s is the collision center of mass energy. Table 4 Charged particles, photon and hadron multiplicity mean values for e + e − → qq collisions, with q = u, d, s, c, b: comparison among SHM predictions, Herwig6510 predictions and LEP data [9,24] at 91.2 GeV center of mass energy. The last two columns contain the discrepancies, measured in standard deviations, for SHM and Herwig6510 predictions with respect to experimental data.