Effects of a phase transition on HBT correlations in an integrated Boltzmann+Hydrodynamics approach
Abstract
A systematic study of HBT radii of pions, produced in heavy ion collisions in the intermediate energy regime (SPS), from an integrated (3+1)d Boltzmann+hydrodynamics approach is presented. The calculations in this hybrid approach, incorporating an hydrodynamic stage into the Ultrarelativistic Quantum Molecular Dynamics transport model, allow for a comparison of different equations of state retaining the same initial conditions and final freezeout. The results are also compared to the pure cascade transport model calculations in the context of the available data. Furthermore, the effect of different treatments of the hydrodynamic freezeout procedure on the HBT radii are investigated. It is found that the HBT radii are essentially insensitive to the details of the freezeout prescription as long as the final hadronic interactions in the cascade are taken into account. The HBT radii and and the ratio are sensitive to the EoS that is employed during the hydrodynamic evolution. We conclude that the increased lifetime in case of a phase transition to a QGP (via a Bag Model equation of state) is not supported by the available data.
pacs:
25.75.Gz,25.75.q,24.10.LxI Introduction
One of the main purposes of the research in heavy ion collisions (HICs) at high beam energies is to explore the existence of the quark gluon plasma (QGP) as well as its properties. The equation of state (EoS) of nuclear matter is one of the key points to gain further understanding since the EoS directly provides the relationship between the pressure and the energy at a given netbaryon density. Phase transitions (PT), e.g., from the hadron resonance gas phase (HG) to the colordeconfined QGP (see e.g., (1); (2); (3)), constitute themselves in changes of the underlying EoS.
Although, on the low temperature side (and for low baryochemical potentials ), investigations of the EoS of nuclear matter have been pursued for many years and uncertainties have been largely reduced, on the high temperature side, the EoS of hot and dense QCD matter is still not precisely known. For systems created in the RHIC energy region with high temperatures and low baryochemical potential, lattice quantum chromodynamics (lQCD) (see, e.g., Ref. (4)) calculations predict a crossover transition between the hadron gas and the QGP phase. The additional structures of the phase diagram are still under heavy debate, especially regarding the existence or nonexistence of a critical endpoint (5).
The intermediate SPS energy regime still raises a lot of interest because the onset of deconfinement is expected to occur at those energies and the possibility of a critical endpoint and a firstorder phase transition is not yet excluded. Several beamenergy dependent observables such as the particle ratios (6); (7), the flow (8); (9); (10), the HBT parameters (11); (12); (13) show a nonmonotonic behaviour around GeV and the interpretation remains still unclear. Therefore, future energy scan programs at RHIC, SPS and FAIR are planned to explore the high region of the phase diagram in more detail.
To learn something about the hot and dense stage of the collision from the final state particle distributions, a dynamical modeling of the whole process is necessary. Some of the important ingredients which have to be considered in a consistent manner are

the initial conditions and the initial nonequilibrium dynamics,

the treatment of the phase transition and hadronization, as well as the right degrees of freedom,

viscosity effects in the initial partonic as well as in the hadronic stage of the evolution,

hadronic rescatterings and freezeout dynamics.
We notice that part of these have been pointed out to be of importance especially for the understanding of the HBT results (14); (15); (16).
Combined microscopic+macroscopic approaches are among the most successful ideas for the modeling of the bulk properties of HICs (17); (18); (19). Recently, a transport approach that embeds a full (3+1) dimensional ideal relativistic one fluid evolution for the hot and dense stage of the reaction has been developed and first results are convincing (20); (21). This hybrid model inherits the advantages of the Ultrarelativistic Quantum Molecular Dynamics (UrQMD) model for the dynamic treatment of the initial and the final state by taking into account eventbyevent fluctuations. Furthermore, the hybrid model allows for a dynamical coupling between hydrodynamics and transport calculation in such a way that one can compare calculations with various EoS during the hydrodynamic evolution and with the pure cascade calculations within the same framework.
It is wellknown that by using HBT interferometry techniques one can get detailed information about the spacetime configuration of the system at freezeout. We concentrate here on the two (identical) pion interferometry and test the sensitivity of the HBT results on different stages of the evolution. In our previous investigations on the HBT correlation of various identical particle pairs from HICs at AGS, SPS, and RHIC energies (22); (23); (24); (13); (25); (26), we adopted the UrQMD model but further considered the mean field potentials for both confined and “preformed” particles in the model (13); (25); (26). We found that also initial stage interaction may contribute to a better description of the HBT timerelated puzzle throughout the energies from AGS, SPS, up to RHIC.
In this paper we perform a systematic investigation of the sensitivity of HBT correlation of negatively charged pions to the EoS by applying the newly developed hybrid approach. Similar more focused studies were frequently discussed with microscopic transport or hydrodynamic models before (11); (27); (28); (29). It is also interesting to study if the current set of EoS employed in the hydrodynamic phase support the conclusion about the origin of the HBT timerelated puzzle. In addition, the effects of the hadronic rescattering and of resonance decays (dubbed as “HR”) after the hydrodynamic freezeout on the HBT radii and the ratio deserve more investigation. We have also noticed that some recent progresses of this topic both from an improved hydrodynamic calculation (16) and from the pionopticalpotential point of view (30) have been published which provides additional new insights.
The paper is arranged as follows. In the next section, the UrQMD+hydrodynamics hybrid model is introduced briefly. The set of different EoS that are employed in the hydrodynamic phase are explained. Two different treatments for the transition process from the hydrodynamic evolution to the final state hadronic cascade are discussed. In Section 3, the analyzing program CRAB for constructing the HBT correlator and the corresponding three dimensional (3D) Gaussian fitting process are introduced. In Section 4, the HBT radii , , and , and the ratio of the negatively charged pion source from central Pb+Pb collisions at SPS energies are shown and discussed in the context of the experimental data. Finally, in Section 5, a summary and an outlook are given.
Ii UrQMD+hydrodynamic model
An integrated Boltzmann+hydrodynamics transport approach is applied to simulate the dynamics of the heavy ion collision. To mimic experimental conditions as realistic as possible the initial conditions and the final hadronic freezeout are calculated using the UrQMD approach. Especially for an observable like HBT radii it is important to take care of the complexity of the different effects (14); (15). The nonequilibrium dynamics, e.g. fluctuations of the local baryon and energy density (31), in the very early stage of the collision and the final state hadronic interactions are properly taken into account on an eventbyeventbasis.
UrQMD is a microscopic transport approach based on the covariant propagation of constituent (anti)quarks and diquarks accompanied by mesons and baryons, as well as the corresponding antiparticles, i.e., full baryonantibaryon symmetry is included. It simulates multiple interactions of ingoing and newly produced particles, the excitation and fragmentation of color strings (32); (33); (34) and the formation and decay of hadronic resonances (35); (36). In principle it is also possible to incorporate mean field interactions in the transport calculation, In the present calculation they are neglected in order to test the hydrophase and in the following we will refer to the pure cascade calculation as UrQMD2.3. Studies on the thermodynamic properties of UrQMD can be found in (37); (38); (39).
The coupling between the UrQMD initial state and the hydrodynamical evolution proceeds when the two Lorentzcontracted nuclei have passed through each other, (20). After that, a full (3+1) dimensional ideal hydrodynamic evolution is performed using the SHASTA algorithm (40); (1). The hydrodynamic evolution is stopped, if the energy density of all cells drops below five (default value) times the ground state energy density (i.e. ). This criterion corresponds to a Tconfiguration where the phase transition is expected  approximately MeV at . The hydrodynamic fields are then mapped to particle degrees of freedom via the CooperFrye equation on an isochronous (in the computational frame) hypersurface. The particle vector information is then transferred back to the UrQMD model, where rescatterings and final decays are calculated using the hadronic cascade. We will further refer to this kind of freezeout procedure as the isochronous freezeout (IF). This procedure is explained in detail in (21).
In this paper we introduce another freezeout procedure to account for the large time dilatation that occurs for fluid elements at large rapidities. Faster fluid elements need a longer time to cool down to the same temperatures than the cells at midrapidity since the hydrodynamic calculation is performed in the centerofmass frame of the collision. At higher energies the isochronous hypersurface increasingly differs from an iso hypersurface ( is the proper time). To mimic an iso hypersurface we therefore freeze out transverse slices, of thickness fm, whenever all cells of that slice fulfill our freezeout criterion. For each slice we apply the isochronous procedure described above separately. By doing this we obtain a rapidity independent freezeout temperature even for the highest beam energies. For lower energies ( GeV) the two procedures yield very similar results for the temperature distributions. The hydrodynamic fields are then again mapped to particle degrees of freedom via the CooperFrye equation on this new hypersurface. In the following we will refer to this procedure as “gradual freezeout”(GF). A more detailed description of the hybrid model including parameter tests and results for multiplicities and spectra can be found in (21).
Serving as an input for the hydrodynamical calculation the EoS strongly influences the dynamics of an expanding system. In this work we use three different EoS to investigate their effect on the extracted HBT radii. The first EoS, named the hadron gas (HG), describes a noninteracting gas of free hadrons (41). Included here are all reliably known hadrons with masses up to GeV, which is equivalent to the active degrees of freedom of the UrQMD model (note that this EoS does not contain any form of phase transition). This purely hadronic calculation serves as a baseline calculation to explore the effects of the change in the underlying dynamics  pure transport vs. hydrodynamic calculation. The second EoS, named the Bag Model EoS (BM), follows from coupling a bag model of massless quarks and gluons to a Walecka type of hadron gas including only SU(2) flavours (for details the reader is referred to (1)). This EoS exhibits a strong firstorder phase transition (with large latent heat) for all baryonic chemical potentials . The third EoS, named the chiral+HG (CH), follows from a chiral hadronic Lagrangian and incorporates the complete set of baryons from the lowest flavour octet, as well as the entire multiplets of scalar, pseudoscalar, vector and axialvector mesons (42). Additional baryonic degrees of freedom are included to produce a firstorder phase transition in certain regimes of the  plane, depending on the couplings (43); (28); (44). Using this EoS, a phase structure including a firstorder phase transition and a critical endpoint at finite is obtained (45). This EoS has already been successfully applied to a hydrodynamic calculation (20).
To visualize the differences of these EoS, Fig. 1 shows the average pressure of the expanding system, from central Pb+Pb collisions at GeV (left plot) and GeV (right plot), as a function of time (in the center of mass frame). The vertical line in each plot indicates the starting time of the hydro evolution. The mean value of the pressure has been obtained by weighting the pressure, , in every cell by its energy density, , and integrating over the hydrodynamic grid
(1) 
All curves in Fig. 1 are plotted until the point in time when the isochronous freezeout criterion is fulfilled. Compared to the HG the BMEoS leads to a delayed freezeout time (i.e. a much longer expansion). While in the first few of the evolution, the system obeying the BMEoS expands most violently (due to the high pressure gradient in the QGP phase), once the system enters the mixed phase, its expansion is slowed down considerably. This can be observed as the “kink” in Fig. 1. At the higher beam energy, GeV, this softening of the EoS is even more pronounced. Since the HGEoS does not contain any phase transition, no softening can be observed, resulting in the shortest expansion time. The chiral CHEoS lies in between both extreme cases. Although a small kink can be observed, it is not as pronounced as in the BMEoS. The effect of changes in the EoS on HBT results has been studied before (28), but the great advantages of our approach are the full (3+1) dimensions and the same initial conditions and freezeout for all three cases and all beam energies without adjusting additional parameters.
Iii CRAB analyzing program and the fitting process
To calculate the twoparticle correlator, the CRAB program is adopted (46), which is based on the formula:
(2) 
Here is an effective probability for emitting a particle with 4momentum from the spacetime point . is the relative wave function with being the relative position in the pair’s rest frame. and are the relative momentum and the average momentum of the two particles and .
In this work, we select central ( of the total cross section ) Pb+Pb collisions at SPS energies: , , , and GeV, with a pair rapidity cut ( is the pair rapidity with pion energies and and longitudinal momenta and in the center of mass system). For each EoS about events are calculated. All particles with their phase space coordinates at freezeout are then given into the CRAB analyzing program. Only the negatively charged pions are considered during the analyzing process (for each analysis, one hundred million pion pairs are considered). For the cascade calculations, we take the results from our previous publications as reference (23); (13). We found that the residual Coulomb effect after the hadron freezeout on the HBT radii of the pion source is small (25), therefore we omit it in the present analysis. Finally, we choose the longitudinal comoving system (LCMS) frame of the pair (also called the “OutSideLong” system, in which the longitudinal component of the pair velocity vanishes), which is frequently adopted in recent years, and fit the correlator by a 3D Gaussian distribution
(3) 
Here is the overall normalization factor, the and are the components of the pair relative momentum and homogeneity length (HBT radius) in the direction, respectively. The parameter is called as the incoherence factor or, more correctly, the intercept parameter and lies for BoseEinstein statistics between 0 and 1 for twoboson correlations in realistic HICs. Because the parameter might be influenced by many additional factors, such as contamination, longlived resonances, or the details of the residual Coulomb modification, we regard it as a free parameter and do not show it in this letter. However, In (25) we have found that the calculated factor with UrQMD can be comparable with experimental data at RHIC energies (although somewhat larger than data). At SPS energies, , is also compatible with experimental data which . The represents the crossterm and plays a role at large rapidity. To fit the correlator with Eq. (3), we use ROOT (47) software and minimize .
Iv HBT results
Fig. 2 shows the transverse momentum () dependence of the HBT radii , , and (at midrapidity) of source from central Pb+Pb collisions at SPS energies. The data (solid stars) are from the NA49 Collaboration (48). The pure cascade calculation is depicted by lines while the hybrid model calculations with different EoS (HG, BM and CH) are depicted by dashed lines with open symbols. As was shown before, the cascade calculation gives a fairly good result of the dependence of and values except at quite small , while for , it is slightly larger than data at large . In contrast, the hybrid model calculations show large HBT (for all employed EoS, but to a varying degree) in all directions, especially in the longitudinal direction. The HG and CH are moderately increased and lead to very similar results for all three directions. The large latent heat in the bag model leads to a further strong increase in the longitudinal direction and in the transverse direction at large . This increase in the BM mode becomes more pronounced at higher beam energies. At first glance, this result might be surprising because at least in the transverse direction one would expect a faster expansion including a hydrodynamic evolution. On the other hand, one knows that the system spends a longer time without emitting any particles in the hybrid model calculation.
Fig. 3 exhibits the freezeout time dependence of the
emission in central Pb+Pb at GeV (left
plot) and GeV (right plot). It is clearly seen
that there are almost no pions emitted before fm in the
hybrid model calculations. This is easy to understand because even
in the gradual hydrofreezeout which is applied here, it takes a
while until the first slices have cooled down and are frozen out
from the hydrodynamic evolution. There is no particle emission from
earlier times in contrast to the pure cascade calculation. For the
BMEoS, this effect is present even for a longer time since the
expansion lasts longer
Fig. 4 illustrates the dependence of the HBT radii under various freezeout conditions, which may be divided into two parts: 1) without HR and 2) with HR after the hydrodynamic phase. Here, without HR (lines) means that the evolution is stopped immediately after the CooperFrye freezeout from the hydrodynamic phase (GF and are adopted as default hydrofreezeout criteria), with instantaneous resonance decays. The observed size of the pion source is small at this hydrodynamic freezeout. In previous investigations (21) it has been found that binary baryonmeson and mesonmeson collisions still frequently happen after the hydrodynamic freezeout. In baryonmeson reactions, the most abundant interactions are the excitation and the decay of the resonance (i.e. ), while in mesonmeson collisions, the process is dominant. A large number of these final hadron interactions in which pions are involved contribute significantly to the final HBT radii of pions in this model.
Let us therefore explore if the finally observed HBT radii do depend on the transition criterion from hydrodynamics to the transport model. The full hybrid model calculations (dashed lines with open symbols) are shown with two different cuts of the energy density ( as default and ) for the GF and with the energy density cut for the IF. It is found that the final state hadronic interactions are sufficient that the effects of different treatments of the hydrodynamic freezeout on the final HBT radii are almost totally washed out in all directions and at all investigated energies.
Let us finally explore the dependence of the ratio on the different EoS and freezeout prescriptions. This ratio was expected to be sensitive to the duration time of the homogeneity region. In Fig. 5 the excitation function of the ratio with the different EoS (lines with solid symbols) and freezeout prescriptions (dashed lines with open symbols) are shown. The bin MeV is chosen. The result for the pure cascade calculation is also shown as a baseline (dotted line). It is seen clearly that the ratio is sensitive to the EoS, but not to the various hydrodynamic freezeout prescriptions when including HR (shown as open triangles and open inverted triangles) as it has already been implied from the results of the HBT radii shown in Figs. 2 and 4. With increasing latent heat which corresponds to the softness of EoS implied from Fig. 1, the ratio is increased. The “excessively” large latent heat in BMEoS results in a long duration time of the pion source and hence a large ratio. Although the overall height is largely overpredicted by the BMEoS, the qualitative behaviour of the data (with a maximal lifetime at beam energies around GeV) is well reproduced. In addition, the “peak structure” is less pronounced than in previous predictions (11), due to the different initial state and seems to provide a more reasonable estimate of the magnitude of the lifetime enhancement. The chiral EoS CH exhibits a lower ratio because the firstorder phase transition is less pronounced. The calculation with HG mode (line with solid squares) leads to the smallest ratio due to the most stiffest EoS among the three ones. The result of the cascade calculation lies in between the CH and the BM modes, which implies a relatively soft EoS. It can be understood since in the pure UrQMD model the new particle production is treated either as a resonance decay or a fragmentation of the string, which introduces a finite lifetime and hence leads to a softer EoS. After considering the mean field potentials for both confined and “preformed” particles (13); (25), which gives a strong repulsion at the early stage, the ratio was seen to decrease in line with results obtained here.
For the full hybrid model calculation the different freezeout prescriptions do not affect the final results (when comparing the results by dashed lines with open triangles with that by the line with solid squares). The calculation with HGEoS but without HR (dashed line with open diamonds) seems to provide better description of the data, but, seems clearly unphysical to the authors as a solution to the duration time problem.
V Summary and outlook
In summary, the HBT correlation of the negatively charged pion source created in central Pb+Pb collisions at SPS energies was investigated with a hybrid model that incorporates a (3+1)d hydrodynamic evolution in the UrQMD transport approach. We explored different settings, one where the EoS was varied without changing the initial conditions and the freezeout prescription and another where the EoS was fixed and the treatment of the freezeout was changed. We presented a systematic investigation of these effects on the HBT radii. It was found that the latent heat influences the emission of particles visibly and hence the HBT radii of the pion source. While the final rescatterings result in an independence of the calculated HBTparameters from the transition criterion, they do not improve the quantitative agreement with the experimental data. The details of the hydrofreezeout prescription do not affect the HBT radii as well as the ratio as long as the HR in the subsequent hadronic transport model were taken into account. Overall, the HBT data seem to favor a stiff EoS (51), but one should also keep in mind that viscosity effects are neglected during the hydrodynamic stage and that the particle emission from the early stages should be handled more carefully.
In the future, bulk and shear viscosity will be further considered for the hydrodynamic phase, and the nonequilibrated particle emission should be treated more precisely.
Acknowledgments
We thank S. Pratt for providing the CRAB program and thank D.H. Rischke for providing the hydrodynamics code. We acknowledge support by the Frankfurt Center for Scientific Computing (CSC). This work was supported by the Hessian LOEWE initiative through the Helmholtz International Center for FAIR (HIC for FAIR) and GSI and BMBF. H.P. acknowledges financial support from the Deutsche Telekom Stiftung and support from the Helmholtz Research School on Quark Matter Studies.
Footnotes
References
 D. H. Rischke, Y. Pursun and J. A. Maruhn, Nucl. Phys. A 595 (1995) 383 [Erratumibid. A 596 (1996) 717].
 C. Spieles, H. Stöcker and C. Greiner, Phys. Rev. C 57 (1998) 908.
 M. Bluhm, B. Kampfer, R. Schulze, D. Seipt and U. Heinz, Phys. Rev. C 76 (2007) 034901.
 Z. Fodor and S. D. Katz, JHEP 0203 (2002) 014.
 P. de Forcrand and O. Philipsen, JHEP 0701 (2007) 077.
 S. V. Afanasiev et al. [The NA49 Collaboration], Phys. Rev. C 66 (2002) 054902.
 C. Alt et al. [NA49 Collaboration], Phys. Rev. C 77 (2008) 024903.
 P. F. Kolb, P. Huovinen, U. W. Heinz and H. Heiselberg, Phys. Lett. B 500 (2001) 232.
 M. Bleicher and H. Stöcker, Phys. Lett. B 526 (2002) 309.
 H. Petersen, Q. Li, X. Zhu and M. Bleicher, Phys. Rev. C 74 (2006) 064908.
 D. H. Rischke and M. Gyulassy, Nucl. Phys. A 608 (1996) 479.
 D. Adamova et al. [CERES Collaboration], Phys. Rev. Lett. 90 (2003) 022301.
 Q. Li, M. Bleicher and H. Stöcker, Phys. Lett. B 659 (2008) 525.
 S. Pratt and J. Vredevoogd, Phys. Rev. C 78 (2008) 054906.
 M. A. Lisa and S. Pratt, arXiv:0811.1352 [nuclex].
 W. Broniowski, M. Chojnacki, W. Florkowski and A. Kisiel, Phys. Rev. Lett. 101 (2008) 022301.
 R. Andrade, F. Grassi, Y. Hama, T. Kodama and O. J. Socolowski, Phys. Rev. Lett. 97 (2006) 202302.
 T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey and Y. Nara, Phys. Lett. B 636 (2006) 299.
 C. Nonaka and S. A. Bass, Phys. Rev. C 75 (2007) 014902.
 J. Steinheimer, M. Bleicher, H. Petersen, S. Schramm, H. Stöcker and D. Zschiesche, Phys. Rev. C 77 (2008) 034901.
 H. Petersen, J. Steinheimer, G. Burau, M. Bleicher and H. Stöcker, Phys. Rev. C 78 (2008) 044901.
 Q. Li, M. Bleicher and H. Stöcker, Phys. Rev. C 73 (2006) 064908.
 Q. Li, M. Bleicher, X. Zhu and H. Stöcker, J. Phys. G 33 (2007) 537.
 Q. Li, M. Bleicher and H. Stöcker, J. Phys. G 34 (2007) 2037.
 Q. Li, M. Bleicher and H. Stöcker, Phys. Lett. B 663 (2008) 395.
 Q. Li and M. Bleicher, J. Phys. G 36 (2009) 015111.
 S. Soff, S. A. Bass and A. Dumitru, Phys. Rev. Lett. 86 (2001) 3981.
 D. Zschiesche, S. Schramm, H. Stöcker and W. Greiner, Phys. Rev. C 65 (2002) 064902.
 see, e.g., the Figs. 15 and 16 in the review article: M. A. Lisa, S. Pratt, R. Soltz and U. Wiedemann, Ann. Rev. Nucl. Part. Sci. 55 (2005) 357.
 M. Luzum, J. G. Cramer and G. A. Miller, Phys. Rev. C 78 (2008) 054905.
 M. Bleicher et al., Nucl. Phys. A 638 (1998) 391.
 B. NilssonAlmqvist and E. Stenlund, Comput. Phys. Commun. 43 (1987) 387.
 B. Andersson, G. Gustafson and B. NilssonAlmqvist, Nucl. Phys. B 281 (1987) 289.
 T. Sjostrand, Comput. Phys. Commun. 82 (1994) 74.
 S. A. Bass et al., Prog. Part. Nucl. Phys. 41 (1998) 255.
 M. Bleicher et al., J. Phys. G 25 (1999) 1859.
 S. A. Bass et al., Phys. Rev. Lett. 81 (1998) 4092.
 L. V. Bravina et al., J. Phys. G 25 (1999) 351. [arXiv:nuclth/9810036].
 L. V. Bravina et al., Phys. Rev. C 78 (2008) 014907. [arXiv:0804.1484 [hepph]].
 D. H. Rischke, S. Bernard and J. A. Maruhn, Nucl. Phys. A 595 (1995) 346.
 D. Zschiesche, S. Schramm, J. SchaffnerBielich, H. Stöcker and W. Greiner, Phys. Lett. B 547 (2002) 7.
 P. Papazoglou, D. Zschiesche, S. Schramm, J. SchaffnerBielich, H. Stöcker and W. Greiner, Phys. Rev. C 59 (1999) 411.
 J. Theis, G. Graebner, G. Buchwald, J. A. Maruhn, W. Greiner, H. Stöcker and J. Polonyi, Phys. Rev. D 28 (1983) 2286.
 D. Zschiesche, G. Zeeb, S. Schramm and H. Stöcker, J. Phys. G 31 (2005) 935.
 D. Zschiesche, G. Zeeb and S. Schramm, J. Phys. G 34 (2007) 1665.
 S. Pratt, CRAB version 3, http://www.nscl.msu.edu /pratt/freecodes/crab/home.html.
 http://root.cern.ch/
 C. Alt et al. [NA49 Collaboration], Phys. Rev. C 77 (2008) 064908.
 F. Grassi, Y. Hama and T. Kodama, Phys. Lett. B 355 (1995) 9.
 J. Knoll, arXiv:0803.2343 [nuclth].
 S. Bekele et al., arXiv:0706.0537 [nuclex].