% I touched two places, find them by searching a string %taian

\documentstyle[12pt,epsf,psfig]{article}
\textwidth 432pt
\flushbottom
\textheight 594pt
\topmargin 72pt
\headheight 0pt
\headsep 0pt
\footskip 54pt
\oddsidemargin 0pt
\parindent 0in
\parskip 3ex
\renewcommand{\textfraction}{0.2}    %float (figures) parameters
\renewcommand{\topfraction}{0.8}
\renewcommand{\bottomfraction}{0.4}
\renewcommand{\floatpagefraction}{0.8}
\newcommand{\lsim}{\:\raise -4pt\hbox{$\stackrel{\textstyle <} {\sim}$}\:}
\begin{document}
\begin{titlepage}
\begin{flushright}
LU TP 94-22\\
Nov. 1994
\end{flushright}

\begin{center}
\large
{\bf AN EVENT GENERATOR FOR THE FIRECRACKER MODEL AND THE RESCATTERING
IN HIGH ENERGY 
pA and AA COLLISIONS\\LUCIAE VERSION 2.0}\\
\vspace{0.3in}
\large
Ben-Hao Sa \footnote{Permanent address: China Institute of Atomic Energy, P.O.Box 275, 
Beijing,102413 China.E-mail address: CIAEDNP\%BEPC2@SCS.SLAC.STANFORD.EDU} and An Tai
\footnote{taian@thep.lu.se} \\
\vspace{0.3in}
Department of Theoretical Physics, University of Lund, \\
S\"{o}lvegatan 14 A, S-223 62 Lund, Sweden
\vspace{0.3in}


\vspace{0.3in}
%Submitted to {\bf Computer Physics Communications}\\
\vspace{0.2in}

\end{center}


\vspace{0.4in}
\large
Abstract: 
LUCIAE  is a Monte Carlo program that, connected to FRITIOF implements
both the Firecracker Model (FCM), a possible mechanism for collective
multi-gluon emission from the colour fields of interacting strings, as
well as the reinteraction of the final state hadrons in relativistic pA
and AA collisions. This paper includes a brief presentation of the
dynamics of the model as well as a description of the program.
\end{titlepage}
%\large
\section*{PROGRAM SUMMARY}

\normalsize

{\it Title of Program:} LUCIAE version 2 \\


{\it Program obtainable from:} taian@thep.lu.se (internet)\\

{\it Computer for which the programme is designed:} DECstation,
VAX, IBM and others with a FORTRAN 77 compiler \\

{\it Computer:} DECstation 3100; installation: Department of
Theoretical Physics, University of Lund, Sweden \\

{\it Operating system:} ULTRIX RISC 4.2 \\

{\it Program language used:} FORTRAN 77 \\

{\it High speed storage required:} $\approx 90$k words\\

{\it No. of bits in a word:} 32 \\

{\it Peripherals used:} terminal for input, terminal or printer
for output \\

{\it No. of lines in combined program and test deck:} about 7000 \\

{\it Keywords:} Monte Carlo, hadrons, nuclei, clusters, collective effects, gluon radiation,
                formation time, rescattering. \\

{\it Nature of the physical problem:} The experiments of  relativistic pA and AA collisions reveal that high energy heavy-ion collisions have some
features which could not be understood in a framework where pA and AA collisions correspond to  simple 
superpositions of independent nucleon-nucleon collisions. The low $P_{t}$ enhancement , high$P_{t}$ enhancement,
 strangeness enhancement etc. clearly show that some collective effects are active in the relativistic 
pA and AA collisions. Formation of a QGP is often suggested to be a candidate for these collective effects.
Could we understand those new features in  the pA and AA collisions by conventional physics? The Monte-Carlo
generator, LUCIAE is built in an attempt to  answer this question.\\


{\it Method of solution:} 
When many strings or colour dipoles are formed in relativistic pA and AA collisions, it is natural to ask if there is some 
interaction between those dipoles so that they could emit gluons collectively using their common total
energy density. The Firecracker Model is developed to study such a collective effect. Moreover, many hadrons are produced 
through hadronisation, which implies that they will interact with each other and with the surrounding 
cold spectator matter. The rescattering effect in the final state is also included in  LUCIAE program.

{\it Restriction of complexity of the problem:}
 At very high energies ($\sqrt{s}$ in the TeV range), especially
 for collisions of massive nuclei, certain arrays need to be expanded to
 accommodate the large number of particles produced.   

{\it Typical running time:} Depends on the type of collision and energy.  Two examples running on Alpha Work Station:
\begin{tabbing}
 $^{16}$O $+$ $^{197}$ Au \qquad \= $p_{\rm lab} = 200$ A GeV/$c$:\qquad \=
 2 events/min (central collision)      \kill
$^{32}$S $+$ $^{32}S$ \> $p_{\rm lab} = 200$ A GeV/$c$: \> 
$\sim$ 15 events/min \\
$^{32}$S $+$ $^{207}Pb$ \> $p_{\rm lab} = 200$ A GeV/$c$: \> 
$\sim$ 4 events/min     
 \end{tabbing}
\vspace{0.3in}


%%%%\large
%\normalsize

\newpage
%%%%\baselineskip 0.8cm

%%%%\input frintro.tex
\large
\section*{LONG WRITE-UP}

\normalsize

\section{Introduction} 
One major reason  to study heavy ion collisions at high energies from
both a theoretical and experimental point of view 
is to search for traces of the predicted Quark-Gluon Plasma (QGP) and
to investigate the properties of this new state of matter. The study
is very important, not only because it provides a test of QCD on a new natural 
scale (an interacting system with high energy density), 
but also because it will improve our understanding of other fundamental
concepts such as the nature of confinement and chiral symmetry breaking. 

Many signals have been theoretically suggested  as possible evidence
for the formation of a QGP. So far, the energy density expected to be
needed for the QGP formation ($\sim$ 1-2 GeV/fm$^{3}$ from lattice QCD
predictions) should have been reached in the CERN 200 AGeV/c S-Au
central collisions. There is, however, no unambiguous (``smoking gun'')
signal in the experimental data at present. All the present predictions
based on the formation of the QGP, although sometimes found in the
data, can be explained also by non-QGP models. This means that it is
not only necessary to consider the properties of the QGP itself but
also to develop ``conventional'' models for possible and consistent
dynamics, including different non-QGP collective effects. Such models
should be used in order to provide possible background for the signals
of QGP formation.

FRITIOF is a model for inelastic hadron-hadron scattering and it has
been successful in describing many experimental data from low energies
via the ISR-regime and all the way  to the top SPS energies, \cite{B.N}\cite{H.P1}. This has been achieved by the introduction of a
particular longitudinal momentum transfer scenarium, gluon
bremsstrahlung radiation (The Dipole Cascade Model, DCM, and the Soft
Radiation Model, SRM, discussed in Section \ref{DCMSRM} --- this is
implemented by the use of ARIADNE, \cite{ari}) as well as hard parton
scattering (Rutherford Parton Scattering, RPS --- this is implemented
by the PYTHIA routines, \cite{py}). In FRITIOF, during the collision
two hadrons are excited due to longitudinal momentum transfers and/or a
RPS. The highly excited states will emit gluon bremsstrahlung according
to the SRM and they are afterwards treated as Lund Strings and are
allowed to decay into final state hadrons according to the Lund
prescription as implemented by JETSET, \cite{jet} .

The FRITIOF model has been extended to also describe hadron-nucleus and 
nucleus-nucleus collisions by assuming that the reactions
are superpositions of binary hadron-hadron collisions (the timescales
involved in the nuclear states are much longer then anyone related to
hadronic interactions --- besides the reinteraction of the emerging
final state hadrons, cf below) in which the geometry of the nucleus
plays an important role because the nuclei should then behave as a
``frozen'' bag of nucleons. This 
picture is supported by the fact that the global features of the heavy ion 
collisions are satisfactorily explained by the collision geometry 
together with the independent hadron-hadron collisions.

On the other hand, in the experiments one finds that when both the
baryon density and the energy density are large enough in high energy
heavy ion collisions, some new features emerge. There is e.g. high
$P_{t}$ and low $P_{t}$ enhancement in the particle spectra and there
is also larger strangeness and anti-baryon production than the one
obtained by a straightforward extrapolation from single hadron
collisions. This means that there are some collective effects.

In the FRITIOF string language we find that there are generally many
excited strings formed close by each other during a relativistic
heavy-ion collisions. In case these strings would behave like vortex
lines in a color superconducting QCD vacuum then it is conceivable that
they will interact (``ordinary'' vortex lines in a type II
superconductor are known to exhibit a repulsive interaction). In the
FCM approach we will assume that there are groups of neighbouring
strings which may form interacting quantum states so that both the
emission of gluonic bremsstrahlung as well as the fragmentation
properties can be affected by the large common energy density, \cite{fire}.
 
In Ref \cite{B.A1} a model for interacting-string (a ``Color Rope'')
fragmentation has been developed based on the fragmentation scheme of a
single string in the Lund Model. Besides the ``individual'' string
energy density (the string tension at rest $\kappa \simeq$ 1 GeV/fm),
there is also some interaction energy in the region of the overlapping
strings (and this is further developed in Ref \cite{TABA}). String
interaction in the DPM model has also discussed in, \cite{C.P}.

 Many hadrons are produced 
through the fragmentation, which implies that they would interact with each other and with the surrounding 
cold spectator matter. The rescattering  has been 
considered important in the relativistic heavy-ion collisions, \cite{sa}.


In this paper we will consider a Monte Carlo implementation of the
Firecracker Model (FCM). This model contains a simple collective
mechanism to emit gluons from a system containing several interacting
strings and the FCM is described in Section \ref{FCM}. We will also
include the rescattering of the many hadrons which are formed in high
energy pA and AA collisions. The rescattering dynamics is described in
Section \ref{Rescatt}--\ref{Res2}. The detailed description of the MC
routines will finally be found in Section \ref{descript}.


\section{The Formation of Color Dipoles and their Radiation in the FRITIOF Model}
\label{DCMSRM}

In the FRITIOF model a hadron is assumed to behave like a massless
relativistic string (MRS) corresponding to a confined color force field
of a vortex line character embedded in a type II color superconducting
vacuum. During the collision between two hadrons it is assumed that
there are many momentum transfers and most of these are expected to be
rather small and uncorrelated. It is further assumed that there is no
net color exchange between the hadrons. The result of the many 
momentum transfers, $q_{j}$, is, based upon the properties of the
classical motion of the MRS, assumed to be essentially additive,
i.e.the string reacts as if it had been worked upon by a single (large)
momentum transfer $Q=\sum q_{j}$ in one of its endpoints. In case one 
of the $q_{j}$'s have a large transverse momentum component, it will
correspond to a Rutherford parton-parton scattering (RPS) and the
probability is treated according to the PYTHIA routines. 

In this way there will be color separation in the model, ie there will
for each hadron be a color $\bar{3}$ (a ``diquark'') continuing forward
along the beam direction and a valence quark, a color $3$, moving in
the opposite direction due to the longitudinal momentum transfer (for
the case of an occurring RPS there are further prescriptions,
\cite{H.P1}). This will lead to bremsstrahlung of a dipole character.
In an $e^+ e^-$-annihilation event there is also a dipole formed
between the originally produced $q$ (color $3$) and $\bar{q}$ (color
$\bar{3}$). The bremsstrahlung gluons are emitted according to the
formula, \cite{L.L}
%
\begin{equation}
d\sigma=\frac{C\alpha_{s}}{\pi}\frac{dk_{\perp}^2}{k_{\perp}^2}dy,
\label{f1}
\end{equation}
%
where $k_{\perp}$ and $y$ are the transverse momentum and rapidity of
the emitted gluon, $C$ is a color factor ($C = 4/3$ for a $(q
\bar{q})$-dipole) and $\alpha_s$ the strong coupling constant. Energy
momentum conservation provides the constraint
%
\begin{equation}
k_{\perp}\cosh(y)\leq \frac{M}{2}                                                                                                                 
\label{f2}
\end{equation}
%
with $M$ the invariant mass of the emitting dipole. This relation is
conveniently approximated by
%
\begin{equation}
k_{\perp} \lsim e^{|y|}M,
\label{f3}
\end{equation}
which means that the kinematically allowed phase space for the
bremsstrahlung gluons corresponds to the interior of a triangular
region in the $y-\log(k_{\perp}^2)$ plane as shown in Fig.1a. 

But contrary to an $e^+ e^-$-annihilation reaction the energy momentum
in a hadronic reaction is not localised in a pointlike color charge but
spread over some extended region. The dynamics of FRITIOF will produce
(at least) two transversely extended fronts from each hadron containing
the energy momentum as well as the color $3$ $q$- and the $\bar{3}$
``diquark'' charge. Any bremsstrahlung emission from an extended source
will necessarily be influenced by this extension (``the form factor'').
In the SRM it is advocated that in case the extension of the emitter is
$l$ and the wave length of the radiation is $\lambda$ then only the
part $(\lambda /l)^d$ of the energy momentum in the front (assuming
that the front has $d$ dimensions) is involved. From comparison to DIS
data it has been found (note that in this case there is one presumably
pointlike ``kicked-out'' quark and one extended ``remainder'' front)
that $d \simeq 1$ (ie the ``spectator'' front behaves in a 
``string-like'' scenarium) and $ l \sim 1/\mu_0$ with $\mu_0 \simeq
0.6$ GeV in case we use $\lambda \simeq 1/k_{\perp}$ for the gluonic
radiation. Therefore the  available phase space represented by
Eq.(\ref{f3}) in the rest frame of the string for the gluon emission is
further restricted into the region  
%
\begin{equation}
k_{\perp} \lsim e^{|y|} \frac{\mu}{k_{\perp}} M,
\label{f5}
\end{equation}
see Fig.1b. The parameter $\mu$ is not necessarily equal to the $\mu_0$ in DIS but
we have found that in case we have two fronts (as in FRITIOF) then a
good description of high energy hadron-hadron data is provided by $l
\rightarrow l/2$, ie $\mu \rightarrow 2 \mu_0$. In the DCM (including
the softening according to the SRM) the gluon emissions will occur so
that the original dipole is divided into two new dipoles $\rightarrow$
three dipoles $\rightarrow$ four dipoles etc. Every time a new gluon is
emitted and a dipole is divided. The building of dipoles corresponds in
the single parton cascades (eg Webber-Marchesini) to ``strong angular
order'', ie a coherence condition for the radiation.

For the hadron-nucleus and nucleus-nucleus collisions, the process has
in the FRITIOF model been regarded as incoherent collisions on the
nucleons. Thus a nucleon from the projectile interacts independently with the
encountered target nucleons as it passes through the nucleus. The probability 
distribution for the number of inelastic collisions $\nu$ is taken from
geometric calculations. Each of the sub-collisions is treated in the
same way as an ordinary hadron-hadron collision, although the momentum
transfers will again be additive and every encounter will make the
projectile more excited. If it interacts with $\nu$ nucleons in the target, $\nu+1$ 
excited string states will be formed as a result. These string states 
will in then independently emit associated bremsstrahlung radiation
based on the Soft Radiation Model described above.

\section{The Physics behind the FCM}
\label{FCM}
\subsection{Gluon Emission from The Colour Field of a Multi-String
State in accordance with the FCM}

In accordance with the FRITIOF scenarium many strings are formed inside
a small space-time volume and it is conceivable that these objects
cannot be treated as independent. 

We have already mentioned the possible relationship between the vortex
lines in a type II superconductor and the string-like color force
fields described by the MRS in the Lund Model. In an electro-magnetic
superconductor type II there are three state regions  when  the
magnetic field $\cal H$ is used as
the order parameter. When the field is below a certain number, $\cal H$$_{1}$, 
the Meissner effect works, but in a region  $\cal H$$_{1}$$\leq$
$\cal H$$\leq$$\cal H$$_{2}$ the quantized vortex lines start to penetrate the 
superconducting state and exhibit a repulsive interaction. Above
$\cal H$$_{2}$ the field is so strong that the whole region ``melts''
into a state of the same character as the break-down of a superconductor 
type I, ie a (space-)extended ``bag'' or ``flux-tube'' region. This is
the usual intuitive picture of the QGP.

The FCM is assumed to work inside this intermediate region (with the
energy density as the order parameter). This means that the strings at
this stage interact but still keep some individuality. We note that the
Eq (\ref{f5}) contains the factor $\mu M$ which effectively corresponds
to an energy density over the region 1/$\mu$. The basic assumption in
the FCM is that {\em  in the very early stage of the collisions when
several ($n$) excited strings (with energy momenta $(P_{+j},P_{-j})$
are close by it may be possible for them to emit gluons using their
common total energy density $\mu \sum P_{\pm j}$ }.

In order to estimate the effect we may for simplicity assume that the
($n$) strings are at rest relative to each other with the same masses
$W_{j}=W$. This means that {\em each string may be able to take energy
from the common pool nW.} Thus one may in principle be able to emit $n$
gluons, each with a transverse momentum limited by $\sqrt{\mu nW}$!

\subsection{The partitioning into Interacting Quantum States}

We will next introduce a way to determine the interacting multi-string
states. The partitioning can be done in different ways but it is
evident that only strings which are ``close by'' each other (in some
way to be defined) should be able to interact and form a common energy pool.

In a hadron nucleus collission we assume that all the participant
nucleons from the target form a single interacting system after the
collission. The impinging hadron is assumed to move along a straight
line path and to interact with the nucleons along the path in
accordance with an overlap function mechanism in the impact parameter.
Therefore we assume that all the participating nucleons (of the target)
are close to each other transversely. In a nucleus-nucleus collission,
however, the projectile nucleons have different impact parameters with
respect to the target nucleus. The strings excited in the projectile
and in the target have in general a large relative longitudinal motion
and we will assume that they form independent interacting sets, to be
called ``clusters'', see Fig.2. 

We will then partition the final state strings from each nucleus into
groups with a small transverse motion with respect to each other. In
order to do that consistently we apply the algorithm which was
introduced in Ref \cite{algoyuri}. There are three main steps in this
momentum space clustering scheme:

\begin{enumerate}

\item  For every pair of projectile or nucleons, $k$ and $l$, the
distance variable $d_{kl}$ is defined by means of the pseudorapidities
$\eta$, azimuthal angles $\phi$ and the transverse momenta $p_t$:
\begin{eqnarray}
d_{kl} =min(p_{t k}, p_{t l}) \sqrt{(\eta_k -\eta_l)^2 + (\phi_k - \phi_l)^2}
\end{eqnarray}
Find the pair with the minimum distance $d_{kl}$ !

\item If this minimum $d_{kl}$ is smaller than a value $d_{cut}$ we
combine $k$ and $l$ into a single ``pseudoparticle'' with a
four-momentum $p_{kl}$ given by
\begin{eqnarray}
p_{kl} = p_k + p_l
\end{eqnarray}

\item Repeat the procedure for all nucleons and pseudoparticles until
all remaining pairs have a distance larger than $d_{cut}$. The
remaining objects now define the clusters.

\end{enumerate}

It is possible to introduce a clustering procedure also in the
corresponding (coordinate) impact parameter space. The resulting
clusters will (although not event by event) be very similar to the ones
obtained above. As it is more in the spirit of the model to consider
energy momentum space (and also because for simulation reasons it is
easier to implement it) we chose to use the momentum space clustering.

In case there is a Rutherford parton interaction on one of the strings
we assume that the coherence between that and the rest of the cluster
is lost. Such a string is not taken into account in connection with the
``collective'' gluon emission in the FCM (although we test if the
Rutherford $k_{\perp R}$ is larger than the possible bremsstrahlung
$k_{\perp}$ in order to fulfill the requirements of perturbation theory,
ie a relevant perturbative contribution is always governed by the
largest momentum transfer). 

We now assume that there is inside a small (Lorentz contracted) space
time region a set of $n$ excited FRITIOF strings formed during the high
energy nuclear interaction in accordance with the clustering
assumptions.  We neglect their transverse motion and we parametrise
their longitudinal motion in terms of lightcone coordinates
$(P_{+j},P_{-j})$ along the beam direction.  The total energy momentum
of the cluster is 
\begin{eqnarray}
P_{\pm tot} = \sum_{j=1}^n P_{\pm j}
\end{eqnarray}
We next assume that the set forms a quantum state of $n$ interacting
objects. The repulsive interaction may be small but if it is there it
means that the strings will be aligned in such a way that they all of
them may use the common total energy density $\mu P_{\pm tot}$ in the
emission of bremsstrahlung gluons. We further assume that
energy-momentum in connection with the emissions is ``locally''
conserved. Therefore a gluon emitted according to the FCM can be
related to a particular string only if (in the corresponding string
cms) the Eq (\ref{f2}) is fulfilled.

In case we produce a gluon in accordance with this procedure, using the
total energy density, then there are several possibilities to relate
the gluon to anyone of the strings. We have considered three such schemes
%taian
\begin{itemize}

\item[A] We emit $n$ gluons from a cluster with $n$ strings. Then we stochastically pick  a gluon and a string out
of the set and try to match them. In case the gluon fulfils Eq (\ref{f2}) 
it is ``accepted'' and the emission allowed, otherwise we drop the gluon and consider this string not to have a firecracker
gluon. Then the procedure is repeated for the next gluon and the next string.\\

\item[B]  As A, but we order the $n$ gluons by their invariant $k_{\perp}$ and also the 
original string set  from the largest mass downwards and then try to
match the $j$-th gluon to the $j$-th string ($j$=1, 2, $\cdots, n$).

\item[C] In this case we emit a gluon without any ordering of the
original set and try to find a string in the set for which Eq
(\ref{f2}) is fulfilled. If there are several then we take
the largest mass string and emit the gluon and remove the string from
the set. If there is no string available that ``fits'' to the gluon
then we remove the largest mass string. Next we emit a second gluon and
repeat the procedure.

\end{itemize}

There is in our experience no large differences among the three
procedures and we have therefore only kept the possibility A as an
option in the Monte Carlo. The main idea is to assume that just as for
a Rutherford scattering a hard FCM gluon breaks the coherence in the
interaction set and therefore the corresponding string state no longer
contributes to the FCM mechanism.  

The result is that the original set of strings is partitioned into
those endowed with a collective FCM gluon and those for which the test
failed and therefore they are still ``longitudinal strings''. These
latter strings can then also emit gluon bremsstrahlung but only
according to the ``individual'' application of the SRM. In all cases
the dipole cascades are continued downwards to a suitable cutoff and
afterwards the strings are fragmented according to the Lund
prescriptions. In that case there is no sensitivity to the gluon
emission cutoff, ie {\em the procedure is infrared stable}. 

\section{Initiation of The Rescattering Transport Process}
\label{Rescatt}
\subsection{The Formation Time Concept}
It is evidently not reasonable to expect that there will be no
(re)interaction when there is a large density of final state hadrons
produced in a dense string surrounding  as conceived in the FRITIOF
scenarium. In this and the next section we will consider these
reinteractions and we will start with the notion of ``formation time''.

In the LUND fragmentation model the probability distribution in invariant time, 
$\tau$, of the $q$$\bar{q}$ breakup vertices (which is only partly what
one would like to call ``the formation time of a hadron'') has been
given  in terms of $\Gamma=(\kappa \tau)^2$ as \cite{B.A2}
%
\begin{equation}
P(\Gamma)d\Gamma \propto \Gamma^a exp(-b\Gamma)d\Gamma
\label{f8}
\end{equation}               
%
with its mean value ,
%
\begin{equation}
<\Gamma> = \frac{a+1}{b}.
\label{f8}
\end{equation}               
%
Here $a$ and $b$ are the fragmentation parameters and $\kappa$ is the energy 
density of the string. Using the parametrical values obtained from
experimental fits the mean formation time $<\tau>$ should be around $1$
fm. For simplicity, we will use this value of $<\tau>$ as a parameter
called {\em the formation time of all the hadrons produced from a
string}. In a moving frame the formation time of a hadron with mass $m$
and energy $E$ is $t=<\tau>E/m$. 

\subsection{The Initiation in Space-Time}

We will assume that from the formation time and onwards a particular
final state hadron will be moving ``outwards'' and therefore be
available for rescattering and so will the the ``spectators'' (ie those
nucleons which did not participate in the high energy encounter). We
also make use of classical trajectories for each one of the produced
hadrons assuming that the produced particles are initially randomly
placed at the formation time in space inside the (classical) overlap
region of the participant and the nucleus. The ``spectator'' nucleons
are likewise randomly distributed inside the non-interacting region and
provided with a thermal motion in accordance with a Boltzman
distribution. We require as always that that there is a minimum
distance between such spectators (``the hard intranuclear core
potential'') of $0.5$ fm.

\subsection{The Initiation in Momentum Space}

The produced particles from the fragmentation  have already had their
four-momentum calculated according to the LUND Fragmentation Model. 
For the projectile spectator nucleons, it is reasonable to assume that they 
share the remaining longitudinal momentum of the projectile, ie the one
not used up to produce hadrons. The transverse momentum of the
spectators are sampled according to the two dimensional Gaussian distribution
\begin{equation}
f(p_t)dp_t = Cp_t\exp[-\frac{p_t^2}{<p_t^2>}]dp_t .                      
\end{equation}   
Here $<p_t^2>$ is assumed to be a few tens of MeV. The transverse momentum of 
the last projectile spectator nucleon should be 
decided by transverse momentum conservation so that the energy-momentum 
conservation is guaranteed.

\subsection{The Classical Motion of The Particles and The Collision Times}

The hadrons are treated  as classical, point-like particles, which move
along straight-line trajectories between two consecutive collisions
\begin{equation}
\vec{r}_i = \vec{R}_i + \vec{\beta}_i(t-T),
\end{equation}
where $\vec{R}_i$ and $T$ is the four-coordinate of the $i$-th particle at the 
time of the last collision.

In order to perform the collisions in a consistent way it is necessary
to introduce a time ordering procedure. We have chosen to calculate the
time between the last and the next collision for all pairs of particles
and then chosen the minimal time. 
Two particles will collide if their distance approaches $d_{min}$ with 
\begin{equation}
\label{dmin}
d_{min} \leq \sqrt{\frac{\sigma_{tot}}{\pi}},
\end{equation}
where $\sigma_{tot}$ is the total cross section in fm$^2$. Since the
minimum distance depends on the  reference frame it is necessary to
view this problem in an inertia frame, where the distance of the
closest approach  will be perpendicular to the momentum of both
particles, \cite{inert}. We have chosen the cms frame of the two
colliding particles. If the particles $i$ and $j$ are chosen as an
example, their equations of motion are
\begin{eqnarray}
\vec{r}^{\ast}_i=\vec{R}^{\ast}_i + 
\vec{\beta}^{\ast}_i(t^{\ast}-T^{\ast}_i) & \mbox{and} &
\vec{r}^{\ast}_j=\vec{R}^{\ast}_j + \vec{\beta}^{\ast}_j(t^{\ast}-T^{\ast}_j).
\end{eqnarray}
Thus one may calculate the relative distance squared between the particle $i$ 
and $j$, $d^2(t^{\ast})$. The minimum distance, $d_{min 1}^2$, in the
orbits and also the collision time is obtained from 
\begin{eqnarray}
\frac{\partial (d^2(t^{\ast}))}{\partial t^{\ast}} = 0 & \mbox{ie} &
 t^{\ast}_{ij} =  \frac{[\Delta \vec{R}^{\ast} \cdot \Delta \vec{\beta}^{\ast} + 
     \Delta \vec{\beta}^{\ast} \cdot (\vec{\beta}^{\ast}_jT^{\ast}_j -
     \vec{\beta}^{\ast}_iT^{\ast}_i)]}{(\Delta \vec{\beta}^{\ast})^2}\\
\Delta \vec{R}^{\ast}= \vec{R}^{\ast}_j - \vec{R}^{\ast}_i & \mbox{and} &
\Delta \vec{\beta}^{\ast} = \vec{\beta}^{\ast}_j - \vec{\beta}^{\ast}_i .
\end{eqnarray}
The $d^2_{min 1}$ must satisfy the Eq (\ref{dmin}) and as the
collisions must be causal the collision time must fulfil
\begin{equation}
t^{\ast}_{ij} \geq Max(T^{\ast}_i,T^{\ast}_j) .
\end{equation}
After the collision times of all pairs have been calculated in their
own cms one has to transfer them back to the lab system. 

\section{The Treatment of the Particle-Particle Collisions}
\label{Res2}
\subsection{The Elastic Scatterings}

A particle-particle collision is always performed in the joint cms, for
both elastic and inelastic scatterings. If the total and the elastic cross section 
satisfies 
\begin{equation}
\frac{\sigma_{el}}{\sigma_{tot}} \geq \xi
\end{equation}
where $\xi$ is a random number, then the particles will be elastically
scattered or else the collision will be considered as inelastic, cf
the next subsection. The distribution of the momentum transfer, $t$,
will for elastic scattering be taken as, \cite{elas}
\begin{eqnarray}
\frac{d\sigma}{dt} &\sim &\exp(Bt) ,\\
B &= &\frac{\epsilon^6}{1+\epsilon^6}A ,\\
\epsilon &= &3.65(\sqrt{s}-m_i-m_j) ,\\
A &= &Min[10.3,(1.12<p_t>)^{-2}]  ,
\end{eqnarray}
where $<p_t>$ is the average transverse momentum of the rescattered
particles. Thus $t$ is calculated by means of the random number $\xi$ as
\begin{eqnarray}
t  =  \frac{1}{B}\ln[\xi+(1-\xi)\exp{(Bt_0)}] & \mbox{with} & t_0 = -4(p^\ast_i)^2 .
\end{eqnarray} 
The polar angle of the scattering is then related to $t$ by
\begin{equation}
t = -2(p^\ast_i)^2(1-\cos(\theta^\ast)).
\end{equation}
and the azimuthal angle will be isotropically distributed.

\subsection{The Inelastic Collisions}

As it is impossible to include all inelastic channels in the program we
have chosen to consider only those reactions which we have found to
give noticeable effects in the final state spectra. The following
reactions are included in LUCIAE2.0:

\begin{itemize}

\item 1) $\pi$ + $\pi$$\rightarrow$ $\overline{K}$+$K $

\item 2) $\pi$ + $N$$\rightarrow$ $\Delta$+$\pi$

\item 3) $\pi$ + $N$$\rightarrow$ $\rho$+$N$

\item 4) $N$ + $N$$\rightarrow$ $\Delta$+$N$

\item 5) $\pi + N$ $\rightarrow$$ K+Y $ 

\item 6) $\pi$ + $\overline{N}$ $\rightarrow$$\overline{K}$+$\overline{Y} $ 

\item 7) $\overline{K}$ + $N$ $\rightarrow$ $\pi+Y $  

\item 8) $K$ + $\overline{N}$ $\rightarrow$ $\pi+\overline{Y} $ 

\item 9) $K$ + $\overline{N}$ $\rightarrow$ $\overline{K}$ +$\overline{\Xi}$

\item 10) $\overline{K}$ + $N$ $\rightarrow$ $K$ + $\Xi$

\end{itemize}
where the hyperons are $Y$=$\Lambda$ or $\Sigma$. Another possible
reaction reaction, $N+N $$\rightarrow$ $N+Y+K$, plays almost no role
due to the smallness of the available phase space and the ensuing rate, \cite{stra}.

The relative probabilities for the different channels in eg $(\pi
N)$-scattering, (the reactions 2,3 and 5 above) is used to determine
the outcome of the inelastic encounter. As the reactions introduced
above do not make up the full inelastic cross section, the remainder is
again treated as elastic encounters.

We will for simplicity assume that the total, elastic and inelastic 
cross sections for $NN$, $\pi N$, $KN$ and $\pi\pi$ interactions are
constants and we will use the following values for them:
$\sigma^{NN}_{tot}$ = 40 mb, $\sigma^{\pi{N}}_{tot}$ = 25 mb
$\sigma^{KN}_{tot}$ = 35 mb and $\sigma^{\pi\pi}_{tot}$ = 10mb. The
ratio of the inelastic 
to the total cross section will be taken as 0.75 $\sim$ 0.85. Further the cross 
section for $\pi\pi\rightarrow{KK}$ is conventionally taken to be  3.0 mb.  

An isospin averaged formula parametrised in Ref \cite{stra} is used for
the cross section of the reaction $\pi + N $$\rightarrow{K+Y}$. The
cross section for the $\overline{K}$$N$ induced strangeness exchange
reactions is assumed to be equal to ten times the corresponding cross
section for the $\pi$$N$ induced 
strangeness production reactions, \cite{stra}. Finally the cross
section for the reactions $\pi\overline{N}$ and $K\overline{N}$  are
assumed to be equal to the corresponding cross sections for the
reactions $\pi N$ and 
$\overline{K} N$ interactions, respectively. 

%taian
The cross sections of inelastic reactions given by the isospin averaged parameterization formulas of 
Ref \cite{stra} for $\pi N$ $\rightarrow{K+Y}$ decrease exponentially with cms energy of the two colliding particles. 
On the other hand we know that the total inelastic cross sections of  $\pi N$ are approximately energy-independent, 
which means that more  inelastic reactions
would take place like $\pi^{-} p$ $\rightarrow K^{*0}+\Lambda$ $\rightarrow K+ \pi+\Lambda$ and   $\pi^{-} p$ 
$\rightarrow K^{*0}+\Sigma^{0}$ $\rightarrow K+ \pi+\Sigma^{0}$ etc., when cms energy of the two colliding particles is
increasing. In some sense the reactions we listed above should be looked upon as ``representative channels''
of certain type of reaction. Therefore, instead of applying  the energy-dependent inelastic cross section, 
we could alternatively give a constant ``effective cross section'' to the reactions $\pi N$ $\rightarrow{K+Y}$.\\

We always assume that the highest isospin state contributes to the
inelastic cross sections  in order to calculate the relative cross
sections for the different isospin channels. As an example, for the
reaction $\pi$ + $N$$\rightarrow$ $\Delta$+$\pi$ we assume that the
isospin $I = 3/2$ state is dominant. By giving a relative cross section
to each inelastic channel the isospin conservation is guaranteed in the
rescattering process.

All the inelastic channels included in LUCIAE2.0 can be found in Appendix A 
together with their relative cross sections.

As usual, two body final state, i.e. the energy and momentum of the scattered 
particles of the inelastic reactions above can be calculated from
two-body kinematics, \cite{inel}, in the cms of the colliding particles, i.e.
\begin{eqnarray}
E^{\ast}_i &= &\frac{s-m^2_j+m^2_i}{2\sqrt{s}} ,\\
|\vec{p}^{\ast}_i| &= &\frac{\sqrt{(s-(m_i+m_j)^2)(s-(m_i-m_j)^2)}}{2\sqrt{s}} 
. 
\end{eqnarray} 
The orientation of the three momentum is isotropically sampled in 4$\pi$ for 
$\pi \pi$ and $NN$ collisions as 
suggested in Ref \cite{elas}, but for the rescattering between a meson and a baryon 
we assume that the final state baryon moves along the direction of the
incoming baryon in the cms of the colliding particles. 

We will finally assume that all the unstable particles produced through the 
rescattering processes will survive until the end of the heavy-ion
collision process.

\section{The Description of the Program}
\label{descript}
LUCIAE2.0 is a subroutine package for simulating collective gluon
emission in the Firecracker Model as well as the rescattering of the
produced particles from FRITIOF in a nuclear environment. The program
is written in FORTRAN 77. It should be used together with
FRITIOF7.02, \cite{H.p12}, JETSET7.3, \cite{jet}, PYTHIA 5.5, \cite{py} and
ARIADNE4.02, \cite{ari}. A user should read 
paper \cite{H.p12} first before using LUCIAE2.0 if he is not familiar
with using FRITIOF7.02. To run with LUCIAE2.0, FRITIOF7.02 and
ARIADNE4.02 have been somewhat modified and are called FRITIOF7.02R and
ARIADNE4.02R, respectively.

One important thing that a user has to keep in mind is that the dimension 
size of the common block COMMON/LUJETS/ (KSZJ) and COMMON/FRPARAI (KSZ1) have 
been extended. Now KSZJ=20000 and KSZ1=30. Those who has used
FRITIOF7.02 before will find it is very easy to use FRITIOF7.02R which
is linked to LUCIAE2.0 since the input and output structure of the two
programs are completely the same. One simply starts with the subroutine FREVENT and 
ends with final state particles (after the firecracker gluon emission and 
rescattering) recorded in event record LUJETS. That is all ! You need ,
of course, to write your own analyses program.

The rescattering part of LUCIAE2.0 is quite time consuming. The time
consumption increases rapidly with the mass of the nuclei and the
collision energy. S+Pb collision at a lab frame energy of 200 AGeV/c, for example, 
the time consumption in this part of the program increases compared to
``plain'' FRITIOF by a factor of 3. If a user is not interested in the
distributions highly affected by the rescattering process, like the
distributions related to strangeness production etc he might like to
run LUCIAE without the rescattering. This is done by setting KFR(21)=0.

\subsection{User interface}
This section contains the most essential information a user should know about 
the program, namely the common block FRPARA1 and LUCIDAT2 for input
parameters and switch controls. The input parameters are provided with
sensible default 
values. Several options are also included to make the program more
accessible to user.

\begin{itemize}
\item{PARAMETER (KSZ1=30) \\
      COMMON/FRPARA1/KFR(KSZ1), VFR(KSZ1)}\\
\begin{description}
\item[Note:] KFR(1-14) and VFR(1-16) are used in FRITIOF7.02. We still list them 
here since some small change has been 
made.
\item[KFR(1)] (D=1)
       Fragmentation
      \begin{description}
       \item[=0] Off.
       \item[=1] On.   
      \end{description}
\item[KFR(2)] (D=1) 
       Multiple gluon emission (dipole radiation) 
      \begin{description}
       \item[=0] Off.
       \item[=1] On.
      \end{description}
\item[KFR(3)] (D=0) Event selection for collisions with a nucleus 
      \begin{description}
       \item[=0] Generate minimum bias events (all interactions recorded).
       \item[=1] Generate only events with all projectile nucleons
                  participated.
       \item[=2] Generate only events with impact parameter between
                  $b_{min}$=VFR(1) and $b_{max}$=VFR(2).
       \item[=3] Apply both requirements in 1 and 2.
      \end{description}
\item[KFR(4)] (D=1)
       Fermi motion in nuclei 
      \begin{description}
       \item[=0] Neglected.
       \item[=1] Included.
      \end{description}
\item[KFR(5)] (D=0) 
       Nucleon-nucleon overlap function 
      \begin{description}
       \item[=0] Eikonal.
       \item[=1] Gaussian.
       \item[=2] Gray disc.
      \end{description}
\item[KFR(6)] (D=2)
       Target Nucleus deformation
      \begin{description}
       \item[=0] No deformation.
       \item[=1] Deformed target nucleus.
       \item[=2] Apply deformation only if the target atomic number $A\geq 80$.  
      \end{description}
\item[KFR(7)] (D=1) 
       Rutherford parton scattering processes
      \begin{description}
       \item[=0] Off.
       \item[=1] On. Here only the hardest RPS is used in FRITIOF.    
       \item[=2] On. The full multiple hard scattering scenario of PYTHIA
               is used.
      \end{description}
\item[KFR(8)] (D=1) 
      Hard gluons cause a corner (soft gluon kink) on the string
      \begin{description}
       \item[=0] No kink is formed. 
       \item[=1] Gluon kink is formed.
      \end{description}
\item[KFR(9)] (D=1)
        `Drowning' of Rutherford parton scattering 
      \begin{description}
       \item[=0] Off. Accept all RPS events.    
       \item[=-1] On. Throw away the drowned RPS event completely 
                  and replace it by a purely soft event.
       \item[=1] As in -1, but the 
            transverse momentum transfer of the soft collision is superimposed
            by the $q_T$ of the drowned RPS.  
      \end{description}
\item[KFR(10)] (D=1) 
        SRM parameters in RPS events: $\mu_1=\mu_0/r$, $\mu_2=\mu_0/(1-r)$
      \begin{description}
      \item[=0] $\mu$ remains the same as in a soft event: $\mu_1=\mu_2=\mu_0$.  
      \item[=1] $r$\,=\,VFR(16). 
      \item[=2] $r$ takes a uniform distribution in (0,1).
      \end{description}
\item[KFR(11)] (D=4)
      Write out of a message when the arguments in FREVENT is changed.  
      \begin{description}
      \item[=-1] Write it out every time the change occurs.    
      \item[=$n$ ($n\geq 0$)] The write out is limited to $n$ times.    
      \end{description}
\item[KFR(12)] (D=2) 
        Set up of the dipole cascade and string fragmentation parameters.
      \begin{description}
      \item[=0] No set up.  The default values are used. 
      \item[=1] Set to the values optimised by 
       OPAL collaboration, \cite{opal}:
       PARA(1)=0.20, PARA(3)=1.0, PARJ(21)=0.37, PARJ(41)=0.18, PARJ(42)=0.34.  
      \item[=2] Set to the values optimised by DELPHI collaboration:
       PARA(1)=0.22, PARA(3)=0.6, PARJ(21)=0.405, PARJ(41)=0.23, PARJ(42)=0.34.  
       \end{description}
\item[KFR(13)] (D=0) 
       Compresses the event record to save space in LUJETS.  This switch
       is particularly needed for heavy ion collisions at high energy
       where LUJETS must be compressed before it gets overfilled.  
      \begin{description}
      \item[=0] Do not compress LUJETS.  
      \item[=1-3] LUEDIT(KFR(13)) is called and LUJETS is compressed. 
        Specifically, for KFR(13)=1 fragmented jets and 
        decayed particles are removed, for KFR(13)=2 neutrinos and
        unknown particles are also removed, and for KFR(13)=3
        neutral particles are further excluded.
      \item[=4] A dummy subroutine FREDITD() is provided as an interface
        in which a user may write his own special purpose codes to edit
        and compress LUJETS.    
      \end{description}
\item[KFR(14)] (D=0) 
       Checked for charge and energy-momentum  conservation
      \begin{description}
      \item[=1] The outcome of each event will be checked.
      \item[=0] No checked. 
       \end{description}       
\item[KFR(15)] (D=1) 
       FIRECRACKER model
      \begin{description}
      \item[=1] On.
      \item[=0] Off.
       \end{description}       
\item[KFR(16)] (D=1) 
       Firecracker gluon
      \begin{description}
      \item[=1] Accept firecracker gluons which are kinematically allowed.
      \item[=0] Reject all firecracker gluons. It is equivalent to KFR(15)=0.
       \end{description}     
\item[KFR(17)] (D=1) 
       The 'drowning' of firecracker gluons by bremsstrahlung radiation.
      \begin{description}
      \item[=1] No drowning.
      \item[=0] Drowning.
       \end{description}     
  \item[KFR(18)] (D=1) 
         Size of firecracker cluster 
      \begin{description}
      \item[=1] The size of firecracker cluster = VFR(17)$<P_{t}>$,
         with $<P_{t}>$ being average transverse momentum of nucleons in the 
cluster.
      \item[=0] size of firecracker cluster =VFR(18).
       \end{description}   
\item[KFR(19)] (D=0) 
       Form of output ,only used when KFR(21)=0
      \begin{description}
      \item[=1] Output in event record includes decayed strings.
      \item[=0] Decayed strings have been taken away in output.
       \end{description}       
\item[KFR(20)] (D=1) 
       The energy dependence of the cross section of an inelastic scattering
      \begin{description}
      \item[=1] Constant cross section.
      \item[=0] Energy-dependent cross section.
       \end{description}       
\item[KFR(21)] (D=1) 
       Rescattering
      \begin{description}
      \item[=1] On.
      \item[=0] Off.
       \end{description}     
\item[KFR(22)] (D=1) 
       Rescattering channel
      \begin{description}
      \item[=1] All channels are included.
      \item[=0] Only the reactions with the yield of the resonants $\Delta$ and 
$\rho$ are included.
       \end{description}     
\item[KFR(23)] (D=0) 
       Checked for charge and energy-momentum  conservation
        after rescattering
      \begin{description}
      \item[=1] The outcome of each event will be checked.
      \item[=0] No checked. 
       \end{description}     

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vskip 12pt

\item[VFR(1)] (D=0.0 fm)
       Minimum impact parameter for options KFR(3)=2 or 3.
\item[VFR(2)] (D=0.2 fm)
       Maximum impact parameter for options KFR(3)=2 or 3.
\item[VFR(3)] (D=0.8 fm)
       The minimum allowable distance $R_{min}$ between nucleons in a nucleus.  
\item[VFR(4-5)] (D=0.2, 0.1)
       Dipole and quadrupole deformation coefficients for deformed target
       nucleus.
\item[VFR(6)] (D=0.01 GeV$^2/c^2$)
       The $<Q_T^2>$ for the Gaussian distribution of soft transverse 
       momentum transfer.
\item[VFR(7)] (D=0.20 GeV$^2/c^2$)
       The $<Q^2_{2T}>$ for the Gaussian distribution of primordial transverse
       momenta on the string ends.       
\item[VFR(8)] (D=0.75 GeV)
       Soft radiation coherence parameter $\mu_0$ for projectile hadron or
       nucleon. 
\item[VFR(9)] (D=0.75 GeV)
       Soft radiation coherence parameter $\mu_0$ for target hadron or nucleon. 
\item[VFR(10-11)] (D=0.0, 0.0 mb)
       Projectile-target nucleon total and elastic cross 
       sections, respectively.
       By default, they are taken from the parameterization 
       of Block and Cahn, \cite{block} (MSTP(31)=5 in PYTHIA).  
       The meson-nucleon cross sections are obtained
       simply by scaling down the Block-Cahn fit. The scale factor is
       $(2/3-a/\sqrt s)$, where $a=1.13$ GeV for pions and $a=3.27$ GeV for 
kaons
       are chosen to reproduce the low energy experimental data. For
       all the other baryons, it is treated as a pion if it is a meson
       and it is treated as a proton if it is a baryon.
       User may override the default by setting VFR(17-18) to positive values.
       However, the user assigned cross sections will only affect the
       N-N interaction probability in nucleus collisions. The
       probability for Rutherford parton scattering is not affected.  
\item[VFR(12)] (D=1.0 GeV/$c$)
       The $q_{Tmin}$ for Rutherford parton scattering.  
\item[VFR(13-15)] (D=1/6, 1/3, 1/2)
       The probabilities for assigning various spins and flavours to the diquark
       end of the string.  For example in a proton, VFR(13-15) are the
       probabilities of finding a $ud$ diquark of spin 1,
       a $uu$ diquark of spin 1, and a $ud$ diquark of spin 0, respectively. 
\item[VFR(16)] (D=0.5)
       The fraction $r$ in option KFR(10)=1. 
\item[VFR(17)] (D=0.75)
       The size of firecracker cluster with respect to the average transverse 
momentum of nucleons 
in the cluster.
       See also KFR(18).
\item[VFR(18)] (D=0.35GeV/c)
       The size of firecracker cluster when KFR(18)=0.
\item[VFR(19)] (D=6.5GeV)
       Minimal CMS energy for starting FIRECRACKER model.
\item[VFR(20)] (D=1.0mb)
      Cross section for  $\pi N$ $\rightarrow{K+Y}$ when KFR(20)=1.
\item[VFR(21)] (D=10.0mb)
      Cross section for $\pi$ + $N$$\rightarrow$ $\Delta$+$\pi$ when KFR(20)=1.
\item[VFR(22)] (D=3.5mb)
       Cross section for $\pi$ + $N$$\rightarrow$ $\rho$+$N$ when KFR(20)=1.
\item[VFR(23)] (D=17.0mb)
       Cross section for $N$ + $N$$\rightarrow$ $\Delta$+$N$ when KFR(20)=1.
  
\end{description}
\end{itemize}
\begin{itemize}
\item COMMON/LUCIDAT2/KFMAXT,PARAM(20),WEIGH(150)\\
\begin{description}
\item[KFMAXT] (D=22)
       The maximum number of particle species included in rescattering.
\item[PARAM(1)] (D=40.0mb)
       The total cross-section of reaction $N-N$ .
\item[PARAM(2)] (D=25.0mb)
       The total cross-section of reaction $\pi-N$  .   
\item[PARAM(3)] (D=35.0mb) 
       The total cross-section of reaction $K-N$.    
\item[PARAM(4)] (D=10.0mb) 
       The total cross-section of  reaction $\pi-\pi$.    
\item[PARAM(5)] (D=3.0mb) 
       The cross section of $\pi+\pi \rightarrow K K$.
\item[PARAM(6)] (D=0.85) 
       The ratio  of inelastic cross-section to  total cross-section.
\item[PARAM(7)] (D=1.5fm) 
       The formation time of a hadron at its rest-frame.
\item[PARAM(8)] (D=0.02fm) 
       The minimal step of time elapsing.
\item[PARAM(9)] (D=0.1) 
       The accuracy of four-momentum conservation when sampling spectator 
nucleons.
\item[PARAM(10)] (D=2.2) 
       The size of effective rescattering region is PARAM(10)$R_{t}$ with 
$R_{t}$ being the
       radius of target.
\item[PARAM(11)] (D=0.16/fm$^3$) 
       The nucleon density of nucleus.
\item[PARAM(12)] (D=0.04 GeV$^2$$/c^2$) 
        $<P_{t}^2>$ for the Gaussian distribution of spectator nucleons.
\item[WEIGH(1)-WEIGH(150)] The relative cross sections of the inelastic reactions, 
see Appendix~A.\
\end{description}
\end{itemize} 
\subsection{Other common blocks}
The following common blocks are used for transmitting data internally. While 
they
can be accessed to read out information, they should never be used to input 
data.
\vskip 12pt
\begin{itemize}
\item COMMON/FCRSOUT/MCL(10),INEL(40),NST(2,KSZ4),ACL(2,KSZ4) \\
\begin{description}
\item[MCL(1): ]  
       The number of cluster from projectile.
\item[MCL(2): ]  
       The number of cluster from target.
\item[MCL(3): ]  
       The number of accepted firecracker gluons.
\item[MCL(4): ]  
       The number of the initial collision pairs of rescattering.
\item[INEL(I): ]  
       The total number that I-th inelastic channel takes place in aggregate. 
Order corresponds
                                 to the list in Appendix A.
\item[NST(L,J): ] 
       The number of strings in J-th cluster. L=1,2 - index to label projectile 
(L=1) and target (L=2).
\item[ACL(L,J): ]  
       The invariant mass of J-th cluster. Index L has the same meaning as 
      in NST(L,J) . 
\end{description}
\end{itemize}
\begin{itemize}
\item PARAMETER (mcludi=20000) \\
   common/wz/c17(mcludi,3),ishp(mcludi),tp(mcludi)\\
\begin{description}
\item[c17(i,j):] 
       The coordinates of i-th particle, j=1,2,3 stands for x,y,z, respectively. 
The 
       origin of the coordinates are at the center of target with z along beam 
direction.
\item[tp(i):]     
       The time of i-th particle.
\item[ishp(i):] 
       The status of i-th particle. =1 inside the rescattering region and =0 
outside.
          \end{description}
          \end{itemize} 

Except information which is stored in the particle list (COMMON/JETSETS/), other 
important information
for proceeding rescattering is stored in collision list. The variables relevant 
to the collision list are:

\begin{description}
\item[lc(i,j):]  i indicates the order number in the collision list,\\

j=1 and 2 give the order number of the colliding particles in the particle list, 
respectively,\\
j=3 and 4 give the flavor code of the scattered particles, respectively, \\ 
j=5 indicates which kind of inelastic scattering is selected. \\
\item[tc(i):]   collision time.\\ 
\item[tw(i):]   The ratio of the reaction cross section of the selected inelastic 
process to the total cross section.
\end{description}
\subsection{The main subroutine}
Following subroutines build the main framework of the program. It might be 
worthwhile to know.
\begin{itemize}
\item SUBROUTINE  FRCLUST\\
to cluster the wounded nucleons by their relative transverse momenta.
\item SUBROUTINE  FRARCLU \\
to emit gluons from the FIRECRACKER Model.
\item SUBROUTINE  FRCRGLU \\
 to assign the firecracker gluons to strings.
\item SUBROUTINE FRINGES(KT,PT,VT,KT0,PT0,NB,N0) \\
 to emit bremsstrahlung gluons from the strings without RPS.
\item SUBROUTINE FRINGEH \\
 to emit bremsstrahlung gluons from the strings with RPS.
\item SUBROUTINE FRRESEXE \\
 to administrate the rescattering of final particles.
\item SUBROUTINE even$_{-}$init(lc,tc,tw,i) \\
 to proceed the initiation of the rescattering transport process.
\item SUBROUTINE rescat(i,ij) \\
 to proceed rescattering.
\item SUBROUTINE ctlcre(lc,tc,tw,n)\\
 to create initial collision list.
\item SUBROUTINE his(time,lc,tc,tw,istop) \\
 to perform classical Newton motion.
\item SUBROUTINE coinel(l,l1,ss,b,pi,pj,icp,pii,pjj,lc,tc,tw,winel,ik3,ik4) \\
 to check if an inelastic reaction could take place. Go head if one channel
       is chosen.
\item SUBROUTINE coelas(l,l1,ss,pi,pj) \\
 to perform elastic scattering.
\item SUBROUTINE updatl(l,l1,time,lc,tc,tw,winel,n) \\
 to update the collision list after the rescattering of a pair of particles. 
\end{itemize}
\subsection{Logical flow of the program}
\begin{enumerate}
\item Cluster formation: After the wounded nucleons are produced from FRITIOF 
those participant nucleons from the
projectile and target are clustered, separately.
\item Emission of firecracker gluons: Gluons are emitted from those clusters and 
matched to strings.
\item Bremsstrahlung radiation:\\
\ \ \ \ \ a) The strings without the RPS gluon are sent into ARIADNE for 
bremsstrahlung radiation.\\
\ \ \ \ \ b) A string with the RPS gluon will first compare the hardness of its RPS 
gluon with that of 
         its firecracker 
      gluon if it happens to have one and the RPS gluon is drowned and replaced 
by the
     firecracker gluon if the firecracker gluon is harder. Then the string is 
sent into  ARIADNE for 
     bremsstrahlung radiation.
\item Fragmentation: The strings after bremsstrahlung radiation will go into JETSET 
for fragmentation.
\item Formation of the initial collision list: The rescattering transport process starts 
by the creation of the initial collision list.
\item Find next collision: From the collision list the collision with the 
earliest collision time is selected to perform collision.
\item Move along Newton trajectory: Each particle moves along its own 
trajectory from the moment of the last collision 
to the current one.
\item Adjudge and perform collision: Adjudge what kind of scattering (elastic or 
inelastic) will take place and perform the 
corresponding scattering.
\item Update collision list: All collisions, in the collision list, involving the
scattering particles i or/and j should be removed from the list. The 
collision list is then updated.
\item Loop: Loop over until no collision pair is left in the collision list.
\end{enumerate}
\begin{itemize} 
\begin{description}
\item [Remark:] If a user wishes to include more rescattering channels besides 
the 90 channels in present program, he could
       change the criteria in subroutine `ctlcre' and the corresponding one in 
subroutine `updatl'. At the same time he should also add
       relevant cross sections and reactions in subroutine `prod'.
\end{description}
\end{itemize} 


\section{Acknowledgments}
One of authors, A.Tai would like to thank B. Andersson for initiating the idea of the FIRECRACKER 
Model and the constant instruction of the work. Ben-Hao Sa would like to thank B. Andersson for kind 
invitation. Both of us would like to thank G. Gustafson for helpful discussions with respect to 
the rescattering. Special thanks also go to Hong Pi, from whom we have benefited a lot through 
discussions with respect to FRITIOF program. This research is supported in part by the EU Programme ``Human Capital
and Mobility", Network ``Physics at High Energy Colliders", contract
CHRX-CT93-0357 (DG 12 COMA).


\newpage
\Large
\begin{center}
\bf {Appendix A}
\end{center}
\normalsize
The following inelastic channels are included in the program(numbers behind are 
relative cross sections):\\

a)$\pi$ + $\pi$$\rightarrow$ $\overline{K}$+$K $.
\begin{enumerate}
\item $\pi^{+}\pi^{-}\rightarrow{K^+K^-}$ ,\hspace{3ex} 1/2
\item $\pi^{+}\pi^{-}\rightarrow{K^0\overline{K^0}}$,\hspace{3ex} 1/2
\item $\pi^{+}\pi^{0}\rightarrow{K^+\overline{K^0}}$,\hspace{3ex} 1
\item $\pi^{0}\pi^{0}\rightarrow{K^+K^-}$ ,\hspace{3ex} 0
\item $\pi^{0}\pi^{0}\rightarrow{K^0\overline{K^0}}$,\hspace{3ex} 0
\item $\pi^{-}\pi^{0}\rightarrow{K^-K^0}$,\hspace{3ex} 1 
\end{enumerate}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b) $\pi + N \rightarrow K + \Lambda$.
\begin{enumerate}
\item $\pi^{-}p\rightarrow{K^0\Lambda}$ , \hspace{3ex} 1
\item $\pi^{+}n\rightarrow{K^+\Lambda}$ ,\hspace{3ex} 1
\item $\pi^{0}p\rightarrow{K^+\Lambda}$ ,\hspace{3ex} 1/2
\item $\pi^{0}n\rightarrow{K^0\Lambda}$ ,\hspace{3ex} 1/2\\
\end{enumerate}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c) $\pi + \overline{N} \rightarrow \overline{K} + \overline{\Lambda}$.
\begin{enumerate}
\item $\pi^{+}\bar{p}\rightarrow{\overline{K^0}\hspace{1mm} \overline{\Lambda}}$ 
, \hspace{3ex} 1
\item $\pi^{-}\bar{n}\rightarrow{K^-\overline{\Lambda}}$ , \hspace{3ex} 1
\item $\pi^{0}\bar{p}\rightarrow{K^-\overline{\Lambda}}$ ,  \hspace{3ex} 1/2
\item $\pi^{0}\bar{n}\rightarrow{\overline{K^0} \hspace{1mm} \overline{\Lambda}}$ ,  \hspace{3ex} 1/2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
e) $\pi + N \rightarrow K + \Sigma$.
\begin{enumerate}
\item $\pi^{-}p\rightarrow{K^+\Sigma^-}$ , \hspace{3ex} 1/9
\item $\pi^{-}p\rightarrow{K^0\Sigma^0}$ ,    \hspace{3ex} 2/9
\item $\pi^{+}p\rightarrow{K^+\Sigma^+}$ , \hspace{3ex} 1
\item $\pi^{+}n\rightarrow{K^+\Sigma^0}$ , \hspace{3ex} 2/9
\item $\pi^{+}n\rightarrow{K^0\Sigma^+}$ , \hspace{3ex} 1/9
\item $\pi^{0}p\rightarrow{K^+\Sigma^0}$ , \hspace{3ex} 4/9
\item $\pi^{0}p\rightarrow{K^0\Sigma^+}$ , \hspace{3ex} 2/9
\item $\pi^{0}n\rightarrow{K^+\Sigma^-}$ , \hspace{3ex} 2/9
\item $\pi^{0}n\rightarrow{K^0\Sigma^0}$ , \hspace{3ex} 4/9
\item $\pi^{-}n\rightarrow{K^0\Sigma^-}$ , \hspace{3ex} 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
f) $\overline{K} + N \rightarrow \pi + \Lambda$.
\begin{enumerate}
\item $K^-p\rightarrow{\pi^{0}\Lambda}$ ,\hspace{3ex} 1/2
\item $K^-n\rightarrow{\pi^{-}\Lambda}$ ,\hspace{3ex} 1
\item $\overline{K^0}p\rightarrow{\pi^{+}\Lambda}$ ,\hspace{3ex} 1
\item $\overline{K^0}n\rightarrow{\pi^{0}\Lambda}$ ,\hspace{3ex} 1/2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
g) $\overline{K} + N \rightarrow \pi + \Sigma$.
\begin{enumerate}
\item $K^-p\rightarrow{\pi^{0}\Sigma^0}$ , \hspace{3ex} 0
\item $K^-p\rightarrow{\pi^{-}\Sigma^+}$ ,\hspace{3ex} 1/2
\item $K^-p\rightarrow{\pi^{+}\Sigma^-}$ ,\hspace{3ex} 1/2
\item $\overline{K^0}p\rightarrow{\pi^{+}\Sigma^0}$ ,\hspace{3ex}1
\item $\overline{K^0}p\rightarrow{\pi^{0}\Sigma^+}$ ,\hspace{3ex} 1
\item $\overline{K^0}n\rightarrow{\pi^{0}\Sigma^0}$ ,\hspace{3ex} 0
\item $\overline{K^0}n\rightarrow{\pi^{-}\Sigma^+}$ ,\hspace{3ex} 1/2
\item $\overline{K^0}n\rightarrow{\pi^{+}\Sigma^-}$ ,\hspace{3ex} 1/2
\item $K^-n\rightarrow{\pi^{0}\Sigma^-}$ ,\hspace{3ex} 1
\item $K^-n\rightarrow{\pi^{-}\Sigma^0}$ ,\hspace{3ex} 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
h) $\overline{K} + N \rightarrow K + \Xi$.
\begin{enumerate}
\item $K^-p\rightarrow{K^{0} \Xi^0}$ , \hspace{3ex} 1/4
\item $K^-p\rightarrow{K^{+} \Xi^-}$ , \hspace{3ex} 1/4
\item $K^-n\rightarrow{K^{0} \Xi^-}$ , \hspace{3ex} 1
\item $\overline{K^0}p\rightarrow{K^{+}\Xi^0}$ , \hspace{3ex} 1
\item $\overline{K^0}n\rightarrow{K^{0}\Xi^0}$ , \hspace{3ex} 1/4
\item $\overline{K^0}n\rightarrow{K^{+}\Xi^-}$ , \hspace{3ex} 1/4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
i) $N + N \rightarrow N + \Delta$.
\begin{enumerate}
\item $pp\rightarrow{\Delta^+p}$ , \hspace{3ex} 1/3
\item $pp\rightarrow{\Delta^{++}n}$ , \hspace{3ex} 1
\item $pn\rightarrow{\Delta^{+}n}$ , \hspace{3ex}1/3
\item $pn\rightarrow{\Delta^{0}p}$ , \hspace{3ex} 1/3
\item $nn\rightarrow{\Delta^{0}n}$ , \hspace{3ex} 1/3
\item $nn\rightarrow{\Delta^{-}p}$ , \hspace{3ex} 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
j) $\pi + N \rightarrow \pi + \Delta$. In the reaction those channels, 
$\pi^0p\rightarrow \Delta^{+}\pi^0 $, 
$\pi^-p\rightarrow \Delta^{0}\pi^0 $, $\pi^+n\rightarrow \Delta^{+}\pi^0$ and 
$\pi^0n\rightarrow \Delta^{0}\pi^0$ are
not included because they have a small relative cross section 1/18.
\begin{enumerate}

\item $\pi^+p\rightarrow{\Delta^{++}\pi^0}$ , \hspace{3ex} 1
\item $\pi^+p\rightarrow{\Delta^{+}\pi^+}$ , \hspace{3ex} 2/3
\item $\pi^0p\rightarrow{\Delta^{0}\pi^+}$ , \hspace{3ex} 16/27
\item $\pi^0p\rightarrow{\Delta^{++}\pi^-}$ , \hspace{3ex} 4/9
\item $\pi^+n\rightarrow{\Delta^{++}\pi^-}$ , \hspace{3ex} 2/9
\item $\pi^+n\rightarrow{\Delta^{0}\pi^+}$ , \hspace{3ex} 8/9
\item $\pi^-p\rightarrow{\Delta^{-}\pi^+}$ , \hspace{3ex} 2/9
\item $\pi^-p\rightarrow{\Delta^{+}\pi^-}$ , \hspace{3ex} 8/27
\item $\pi^0n\rightarrow{\Delta^{+}\pi^-}$ , \hspace{3ex} 16/27
\item $\pi^0n\rightarrow{\Delta^{-}\pi^+}$ , \hspace{3ex} 4/9
\item $\pi^-n\rightarrow{\Delta^{-}\pi^0}$ , \hspace{3ex} 1
\item $\pi^-n\rightarrow{\Delta^{0}\pi^-}$ , \hspace{3ex} 2/3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
k) $\pi + N \rightarrow \rho + N$. 
\begin{enumerate}
\item $\pi^-p\rightarrow{\rho^{0}n}$ ,\hspace{3ex}2/9
\item $\pi^-p\rightarrow{\rho^{-}p}$ ,\hspace{3ex}1/9
\item $\pi^-n\rightarrow{\rho^{-}n}$ ,\hspace{3ex}1
\item $\pi^+n\rightarrow{\rho^{0}p}$ ,\hspace{3ex}2/9
\item $\pi^+n\rightarrow{\rho^{+}n}$ ,\hspace{3ex} 1/9
\item $\pi^+p\rightarrow{\rho^{+}p}$ ,\hspace{3ex} 1
\item $\pi^0p\rightarrow{\rho^{+}n}$ ,\hspace{3ex} 2/9
\item $\pi^0n\rightarrow{\rho^{0}n}$ ,\hspace{3ex} 4/9
\item $\pi^0p\rightarrow{\rho^{0}p}$ ,\hspace{3ex} 4/9
\item $\pi^0n\rightarrow{\rho^{-}p}$ ,\hspace{3ex} 2/9
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
l) $K + \overline{N} \rightarrow \pi + \overline{\Lambda}$.
\begin{enumerate}
\item $K^+\bar{p}\rightarrow{\pi^{0}\overline{\Lambda}}$ , \hspace{3ex}1/2
\item $K^+\bar{n}\rightarrow{\pi^{+}\overline{\Lambda}}$ , \hspace{3ex}1
\item $K^0\bar{p}\rightarrow{\pi^{-}\overline{\Lambda}}$ , \hspace{3ex}1
\item $K^0\bar{n}\rightarrow{\pi^{0}\overline{\Lambda}}$ , \hspace{3ex}1/2
\end{enumerate}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m) $K + \overline{N} \rightarrow \overline{K} + \overline{\Xi}$.
\begin{enumerate}
\item $K^+\bar{p}\rightarrow{\overline{\Xi^0}\ \overline{K^0}}$ ,  \hspace{3ex} 
1/4
\item $K^+\bar{p}\rightarrow{\overline{\Xi^-} K^-}$ ,  \hspace{3ex} 1/4
\item $K^+\bar{n}\rightarrow{\overline{\Xi^-}\ \overline{K^0}}$ ,  \hspace{3ex} 
1
\item $K^0\bar{p}\rightarrow{\overline{\Xi^0}K^-}$ ,  \hspace{3ex}1
\item $K^0\bar{n}\rightarrow{\overline{\Xi^0}\ \overline{K^0}}$ ,  
\hspace{3ex}1/4
\item $K^0\bar{n}\rightarrow{\overline{\Xi^-}K^-}$ ,  \hspace{3ex} 1/4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
n) $\pi + \overline{N} \rightarrow \overline{K} + \overline{\Sigma^0}$. The 
following two reactions have a low reaction rate 
because of the poor aboundance of $\overline{N}$, so we only include the 
production of $\overline{\Sigma^0}$ 
which decays into $\overline{\Lambda}$.
\begin{enumerate}
\item $\pi^{+}\bar{p}\rightarrow{\overline{K^0}\ \overline{\Sigma^0}}$ ,  
\hspace{3ex} 2/9
\item $\pi^{-}\bar{n}\rightarrow{K^-\overline{\Sigma^0}}$ ,  \hspace{3ex} 2/9
\item $\pi^{0}\bar{p}\rightarrow{K^-\overline{\Sigma^0}}$ ,  \hspace{3ex}4/9
\item $\pi^{0}\bar{n}\rightarrow{\overline{K^0}\ \overline{\Sigma^0}}$ ,  
\hspace{3ex}4/9
\end{enumerate}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
o) $K + \overline{N} \rightarrow \pi + \overline{\Sigma^0}$.
\begin{enumerate}
\item $K^{+}\bar{p}\rightarrow{\pi^0\overline{\Sigma^0}}$ ,  \hspace{3ex} 0
\item $K^{0}\bar{p}\rightarrow{\pi^-\overline{\Sigma^0}}$ ,  \hspace{3ex} 1
\item $K^{0}\bar{n}\rightarrow{\pi^0\overline{\Sigma^0}}$ ,  \hspace{3ex} 0
\item $K^{+}\bar{n}\rightarrow{\pi^+\overline{\Sigma^0}}$ .  \hspace{3ex} 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{enumerate}
\newpage
\Large
\begin{center}
\bf {Appendix B}
\end{center}
\normalsize
The following is a sample main program.\\
\renewcommand{\thesection}{Appendix \Alph{section}}
\renewcommand{\theequation}{\mbox{A}\arabic{equation}}
\setcounter{equation}{0}
\begin{verbatim}
C..This program generates a few sample  LUCIAE events, and then does
C..histogram for negatively charged particle multiplicity distribution
C..in O+Au collision at 200 GeV/nucleon lab energy.
C..This routine, loaded together with (FRITIOF_7.02R, ARIADNE_4.02R,
C..PYTHIA5.5 and JETSET7.3) can be used to test the installation of programs.  
C---------------------------------------------------------------------

      PARAMETER (KSZJ=20000,KSZ1=30)
C...    **** Be sure to check that all the KSZJ's in MAIN, Fritiof,
C...         Jetset and Pythia are identically set      *****
      COMMON/FRPARA1/KFR(KSZ1),VFR(KSZ1)
      COMMON/FRINTN0/PLI0(2,4),AOP(KSZ1),IOP(KSZ1),NFR(KSZ1)
      COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/LUJETS/N,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
      COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
      DIMENSION MP(0:300)

C...Open a file to take the write out of the program:
      MSTU(11) = 20
      OPEN(MSTU(11),FILE='test.out',STATUS='unknown')

C:::::::Multiplicity distribution for O+AU collision at 200 GeV ::::::::

C...Forbid the decays of Lambda and K_S0: 
      MDCY(LUCOMP(3122),1) = 0
      MDCY(LUCOMP(310),1) = 0

C...Book spaces for the histogram (or use a histogram package):
      DO 50 J=0,300
50    MP(J) = 0

C...Test 10 events (of course a lot more events are needed realistically):
      NEVENT=10
      NTRIG = 0
      DO 100 I=1, NEVENT

      CALL FREVENT('FIXT','O','AU',200.)

C...Output the event using JETSET routine LULIST:
      IF(I.LE.3) CALL LULIST(1)

C...Edit the event record, remove partons or decayed particles:
      CALL LUEDIT(1)

C...Assume a trigger requiring that the energy in the forward cone
C...(theta < 0.3 degree) must be less than 60% of the total beam energy.
C...Also find out the number of negatively charged particles:
      IQTRIG = 0
      EFWD = 0.
      N_=0
      DO 70 J=1, N
      THETA = PLU(J,14)
C...  (PLU is a JETSET function.  Please refer to the JETSET manual.)
      IF(THETA.LT.0.3) EFWD = EFWD+PLU(J,4)

C...Count the negative particles.  Spectator nuclei, which have codes
C...ABS(K(J,2))=10000+N_proton, must be excluded:
      IF(ABS(K(J,2)).LT.10000) THEN
        IF(PLU(J,6).LT.0.) N_ = N_+1
      ENDIF
70    CONTINUE

      EBEAM = 200.*IOP(3)
      IF(EFWD.LT.0.6*EBEAM) IQTRIG = 1

C...Do histogram:
      IF(IQTRIG.EQ.1) THEN
      NTRIG = NTRIG+1
      MP(N_) = MP(N_)+1
      ENDIF

100   CONTINUE

C...Output the histogram data:  
      WRITE(MSTU(11),500) NEVENT, NTRIG
      DO 200 J=0,300
200   WRITE(MSTU(11),*) J, FLOAT(MP(J))/FLOAT(NTRIG)
500   FORMAT(X,'Number of events:',I4,2x,'Triggered events:',I4)

C...Write out the values of the parameters and some statistics:
      CALL FRVALUE(0)
 		
      CLOSE (MSTU(11))
      END
\end{verbatim}


\newpage
\begin{thebibliography}{15}
\bibitem{B.N} B. Andersson, G. Gustafson and B. Nilsson-Almqvist, ``A High 
Energy String Dynamics Model for Hadronic 
Interactions'', $Lund preprint$ LU$-$TP$-87-6$.
\bibitem{H.P1}  B. Andersson, G. Gustafson and H. Pi, $Z. Phys$. {\bf C57}(1993) 
485.
\bibitem{ari} L.\ L\"onnblad, 
  ``{\tt ARIADNE version 4}, A Program for Simulation of QCD Cascades 
implementing the Colour Dipole
     Model'', $DESY$ 92-046.
\bibitem{py} H. -U. Bengtsson and T. Sj\"ostrand, 
   {\it Comput. Phys. Commun.} {\bf 46} (1987) 43 \\
   An updated version of this paper, ``A Manual to The Lund Monte Carlo
   for Hadronic Processes:  PYTHIA version 5.5'', is available upon request
   to the authors.
\bibitem{jet} T. Sj\"ostrand, 
   ``A Manual to The Lund Monte Carlo for Jet Fragmentation and
     $e^+e^-$ Physics: JETSET version 7.3'', available upon request
     to the author.
\bibitem{fire} Bo Andersson, $Phys. Lett. ${\bf B256} (1991)337
\bibitem{B.A1} B. Anderssion and Peter A. Henning. $Nucl. Phys.$ {\bf 
B355}(1991) 82-105.
\bibitem{TABA} An Tai and Bo Andersson, to be published.
\bibitem{C.P} C. Pajares, ``String Production Mechanisms'' in Proceedings of the 
International Workshop on Quark Gluon 
Plasma Signatures in Strasbourg(1990)21.
\bibitem{sa}Sa Ben-Hao, Wang Zhong-Qi, Zhang Xiao-Ze, Song Guang, Lu 
Zhong-Dao \\ and           Zhang Yu-Ming, $Phys. Rev.$ {\bf C 48}(1993)2995
\bibitem{L.L} B. Andersson, G. Gustafson, L. L\"{o}nnblad and U. Petterson, 
$Z.Phys$. {\bf C43}(1989)625.
\bibitem{algoyuri} S. Catani, Yu. L. Dokshitzer, M. H. Seymour and B. R. Webber, 
``Longitudinally$-$Invariant 
$K_{\perp}$$-$Clustering Algorithms fo6r Hadron$-$Hadron Collision'' $Lund 
preprint$ LU$-$TP$-$93-2.
\bibitem{B.A2} B. Andersson, G. Gustafson and B. S\"{o}derberg, $Z. Phys.$ {\bf 
C20}(1983) 317.
\bibitem{inert} S. P. Sorensen, private communication.
\bibitem{elas} J. Cugnon, T. Mizutani, and J. Vandermeulen, $Nucl. Phys$., {\bf 
A352},
      505(1981).
\bibitem{stra} P. Koch, B. M$\ddot{u}$ller and J. Rafeiski, $Phys. Rep.$ {\bf 
142},167(1986).
\bibitem{inel}  Particle data group, $Phys. Lett.$, {\bf B239}, (1990).
\bibitem{H.p12} H. Pi, $Comput. Phys. Commun$. {\bf 71} (1992) 173-192
\bibitem{opal} M. Z. Akrawy {\it et al}, 
   {\it Z. Phys.} {\bf C47} (1990) 505.
\bibitem{block} M. M. Block and R. N. Cahn, 
   in {\it Physics Simulations at High Energy}, edited by V. Barger, 
   T. Gottschalk and F. Halzen (World Scientific, Singapore, 1987), p. 89.


\end{thebibliography}

 \newpage
\parskip 6ex




\begin {center}

\psfig{file=plotfr1.eps,height=7cm,width=8cm}
\mbox{ }

Fig.1a: The phase space limit for the emission of a gluon from a colour dipole \\
\hspace{-5.2cm}with no transverse extension.



\psfig{file=plotfr2.eps,height=7cm,width=8cm} 
\mbox{ }

Fig.1b: The available phase space for gluon emission from an extended colour 
dipole.

\newpage
\parskip 6ex
\psfig{file=plotfr3.eps,height=3cm,width=13cm}

\mbox{ }
Fig.2: A cluster of nucleons after collisions.



%\epsfig{file=plotfr4.eps,height=8cm,width=9cm} 

\mbox{ }
%Fig.3:The available phase space for firecracker gluon emission. 

\end{center}                       


\end{document}

