\chapter{Parametrization Driven Models}

\section{Introduction}

Two sets of parameterized models are provided for the simulation of high 
energy hadron-nucleus interactions.  The so-called ``low energy model'' is 
intended for hadronic projectiles with incident energies between 1 GeV and 
25 GeV, while the ``high energy model'' is valid for projectiles between 25 
GeV and 10 TeV.  Both are based on the well-known GHEISHA package of GEANT3.
The physics underlying these models comes from an old-fashioned multi-chain 
model in which the incident particle collides with a nucleon inside the 
nucleus.  The final state of this interaction consists of a recoil nucleon, 
the scattered incident particle, and possibly many hadronic secondaries.  
Hadron production is approximated by the formation zone concept, in which the 
interacting quark-partons require some time and therefore some range to 
hadronize into real particles.  All of these particles are able to re-interact
within the nucleus, thus developing an intra-nuclear cascade. \\ 

\noindent
In these models only the first hadron-nucleon collision is simulated in 
detail.  The remaining interactions within the nucleus are simulated by 
generating additional hadrons and treating them as secondaries from the 
initial collision.  The numbers, types and distributions of the extra 
hadrons are determined by functions which were fitted to experimental data 
or which reproduce general trends in hadron-nucleus collisions.  Numerous 
tunable parameters are used throughout these models to obtain reasonable
physical behavior.  This restricts the use of these models as generators for 
hadron-nucleus interactions because it is not always clear how the parameters
relate to physical quantities.  On the other hand a precise simulation of 
minimum bias events is possible, with significant predictive power for 
calorimetry.  

\section{Low Energy Model}

In the low energy parameterized model the mean number of hadrons produced in 
a hadron-nucleus collision is given by  
\begin{equation} 
\label{pm.eq1}
 N_m = C(s) A^{1/3} N_{ic} 
\end{equation}
where $A$ is the atomic mass, $C(s)$ is a function only of the center of 
mass energy $s$, and $N_{ic}$ is approximately the number of hadrons generated
in the initial collision.  Assuming that the collision occurs at the center of 
the nucleus, each of these hadrons must traverse a distance roughly equal to 
the nuclear radius.  They may therefore potentially interact with a number of 
nucleons proportional to $A^{1/3}$.  If the energy-dependent cross section 
for interaction in the nuclear medium is included in $C$ then Eq. \ref{pm.eq1}
can be interpreted as the number of target nucleons excited by the initial 
collision.  Some of these nucleons are added to the intra-nuclear cascade.
The rest, especially at higher momenta where nucleon production is suppressed,
are replaced by pions and kaons. \\

\noindent
Once the mean number of hadrons, $N_m$ is calculated, the total number of 
hadrons in the intra-nuclear cascade is sampled from a Poisson distribution 
about the mean.  Sampling from additional distribution functions provides
\begin{itemize}
\item the combined multiplicity $w(\vec{a},n_{i})$ for all particles
$i$, $i = \pi^{+}, \pi^{0}, \pi^{-}, p, n, .....$, including the
correlations between them,
\item the additive quantum numbers $E$ (energy), $Q$ (charge), $S$ 
(strangeness) and $B$ (baryon number) in the entire phase space region, and
\item the reaction products from nuclear fission and evaporation.
\end{itemize}
A universal function $f(\vec{b},x/p_{T},m_{T})$ is used for the distribution
of the additive quantum numbers, where $x$ is the Feynman variable, $p_T$ is
the transverse momentum and $m_T$ is the transverse mass.  $\vec{a}$ and 
$\vec{b}$ are parameter vectors, which depend on the particle type of the 
incoming beam and the atomic number $A$ of the target.  Any dependence on 
the beam energy is completely restricted to the multiplicity distribution 
and the available phase space. \\ 

\noindent 
The low energy model can be applied to the $\pi^+$, $\pi^-$, $K^+$, $K^-$,
$K^0$ and $\overline{K^0}$ mesons.  It can also be applied to the baryons $p$,
$n$, $\Lambda$, $\Sigma^+$, $\Sigma^-$, $\Xi^0$, $\Xi^-$, $\Omega^-$, and 
their anti-particles, as well as the light nuclei, $d$, $t$ and $\alpha$.  The 
model can in principal be applied down to zero projectile energy, but the 
assumptions used to develop it begin to break down in the sub-GeV region.


\section{High Energy Model}

The high energy model is valid for incident particle energies from 10-20 GeV 
up to 10-20 TeV.  Individual implementations of the model exist for $\pi^+$, 
$\pi^-$, $K^+$, $K^-$, $K^0_S$ and $K^0_L$ mesons, and for $p$, $n$, 
$\Lambda$, $\Sigma^+$, $\Sigma^-$, $\Xi^0$, $\Xi^-$, and $\Omega^-$ baryons 
and their anti-particles.

\subsection{Initial Interaction}
In a given implementation, the generation of the final state begins with the
selection of a nucleon from the target nucleus.  The pion multiplicities 
resulting from the initial interaction of the incident particle and the target
nucleon are then calculated.  The total pion multiplicity is taken to be a 
function of the log of the available energy in the center of mass of the 
incident particle and target nucleon, and the $\pi^+$, $\pi^-$ and $\pi^0$ 
multiplicities are given by the KNO distribution.

From this initial set of particles, two are chosen at random to be replaced 
with either a kaon-anti-kaon pair, a nucleon-anti-nucleon pair, or a kaon and 
a hyperon.  The relative probabilities of these options are chosen according 
to a logarithmically interpolated table of strange-pair and 
nucleon-anti-nucleon pair cross sections.  The particle types of the 
pair are chosen according to averaged, parameterized cross sections typical at
energies of a few GeV.  If the increased mass of the new pair causes the total
available energy to be exceeded, particles are removed from the initial set as
necessary.

\subsection{Intra-nuclear Cascade}
The cascade of these particles through the nucleus, and the additional
particles generated by the cascade are simulated by several models.  These 
include high energy cascading, high energy cluster production, medium energy 
cascading and medium energy cluster production.  For each event, high energy 
cascading is attempted first.  If the available energy is sufficient, this 
method will likely succeed in producing the final state and the interaction 
will have been completely simulated.  If it fails due to lack of energy or 
other reasons, the remaining models are called in succession until the final 
state is produced.  If none of these methods succeeds, quasi-elastic 
scattering is attempted and finally, as a last resort, elastic scattering is 
performed.  These models are responsible for assigning final state momenta to 
all generated particles, and for checking that, on average, energy and 
momentum are conserved. 

\subsection{High Energy Cascading}

As particles from the initial collision cascade through the nucleus more
particles will be generated.  The number and type of these particles are
parameterized in terms of the CM energy of the initial particle-nucleon
collision.  The number of particles produced from the cascade is given 
roughly by 
\begin{equation}
\label{he.eq1}
 N_m = C(s) [A^{1/3} - 1] N_{ic}
\end{equation}
where $A$ is the atomic mass, $C(s)$ is a function only of $s$, the square of 
the center of mass energy, and $N_{ic}$ is approximately the number of hadrons
generated in the initial collision. This can be understood qualitatively 
by assuming that the collision occurs, on average, at the center of the 
nucleus. Then each of the $N_{ic}$ hadrons must traverse a distance 
roughly equal to the nuclear radius.  They may therefore potentially interact 
with a number of nucleons proportional to $A^{1/3}$.  If the energy-dependent 
cross section for interaction in the nuclear medium is included in $C(s)$ then 
Eq. \ref{he.eq1} can be interpreted as the number of target nucleons excited 
by the initial collision and its secondaries.  

Some of these nucleons are added to the intra-nuclear cascade. The rest, 
especially at higher momenta where nucleon production is suppressed, are 
replaced by pions, kaons and hyperons.  The mean of the total number of
hadrons generated in the cascade is partitioned into the mean number of 
nucleons, $N_n$, pions, $N_\pi$ and strange particles, $N_s$.  Each of these
is used as the mean of a Poisson distribution which produces the randomized
number of each type of particle.  

The momenta of these particles are generated by first dividing the final state
phase space into forward and backward hemispheres, where forward is in the 
direction of the original projectile.  Each particle is assigned to one 
hemisphere or the other according to the particle type and origin: 

\begin{itemize} 

\item the original projectile, or its substitute if charge or strangeness 
exchange occurs, is assigned to the forward hemisphere and the target nucleon 
is assigned to the backward hemisphere;

\item the remainder of the particles from the initial collision are assigned 
at random to either hemisphere;

\item pions and strange particles generated in the intra-nuclear cascade are 
assigned 80\% to the backward hemisphere and 20\% to the forward hemisphere;

\item nucleons generated in the intra-nuclear cascade are all assigned 
to the backward hemisphere.

\end{itemize}

It is assumed that energy is separately conserved for each hemisphere.  If
too many particles have been added to a given hemisphere, randomly chosen 
particles are deleted until the energy budget is met.  The final state 
momenta are then generated according to two different algorithms, a cluster 
model for the backward nucleons from the intra-nuclear cascade, and a
fragmentation model for all other particles.  Several corrections are then 
applied to the final state particles, including momentum re-scaling, effects 
due to Fermi motion, and binding energy subtraction.  Finally the 
de-excitation of the residual nucleus is treated by adding lower energy 
protons, neutrons and light ions to the final state particle list. \\       

\noindent
{\bf Fragmentation Model.} This model simulates the fragmentation of the
highly excited hadrons formed in the initial projectile-nucleon collision.
Particle momenta are generated by first sampling the average transverse 
momentum $p_T$ from an exponential distribution:

\begin{equation}
\label{he.eq2}
 exp [-a {p_T}^b ]
\end{equation}
where 
\begin{eqnarray}
 1.70 \le a \le 4.00 ; \
 1.18 \le b \le 1.67 .
\label{he.eq3} 
\end{eqnarray}

The values of $a$ and $b$ depend on particle type and result from a
parameterization of experimental data.  The value selected for $p_T$ is 
then used to set the scale for the determination of $x$, the fraction of
the projectile's momentum carried by the fragment.  The sampling of $x$
assumes that the invariant cross section for the production of fragments 
can be given by

\begin{equation}
\label{he.eq4}
 E \frac{d^3 \sigma}{dp^3} = \frac{K}{(M^2 x^2 + {p_T}^2)^{3/2}}
\end{equation}  
where $E$ and $p$ are the energy and momentum, respectively, of the produced
fragment, and $K$ is a proportionality constant.  $M$ is the average 
transverse mass which is parameterized from data and varies from 0.75 GeV to 
0.10 GeV, depending on particle type.  Taking $m$ to be the mass of the 
fragment and noting that 
\begin{equation}
   p_z \simeq x E_{proj}
\label{he.eq5}
\end{equation}
in the forward hemisphere and
\begin{equation}
   p_z \simeq x E_{targ}
\label{he.eq6}
\end{equation}
in the backward hemisphere, Eq. \ref{he.eq4} can be re-written to give the 
sampling function for $x$:
\begin{equation}
 \frac{d^3 \sigma}{dp^3} = \frac{K}{(M^2 x^2 + {p_T}^2)^{3/2}} 
           \frac{1}{ \sqrt{m^2 + {p_T}^2 + x^2 E_i^2} } ,
\label{he.eq7}
\end{equation}
where $i = proj$ or $targ$.

$x$-sampling is performed for each fragment in the final-state candidate list.
Once a fragment's momentum is assigned, its total energy is checked to see 
if it exceeds the energy budget in its hemisphere.  If so, the momentum of 
the particle is reduced by 10\%, as is $p_T$ and the integral of the 
$x$-sampling function, and the momentum selection process is repeated.  If 
the offending particle starts out in the forward hemisphere, it is moved to 
the backward hemisphere, provided the budget for the backward hemisphere is 
not exceeded.  If, after six iterations, the particle still does not fit, it
is removed from the candidate list and the kinetic energies of the particles 
selected up to this point are reduced by 5\%.  The entire procedure is 
repeated up to three times for each fragment.

The incident and target particles, or their substitutes in the case of charge-
or strangeness-exchange, are guaranteed to be part of the final state.
They are the last particles to be selected and the remaining energy in their
respective hemispheres is used to set the $p_z$ components of their momenta.
The $p_T$ components selected by $x$-sampling are retained. \\ 

\noindent
{\bf Cluster Model.} This model groups the nucleons produced in the 
intra-nuclear cascade together with the target nucleon or hyperon, and 
treats them as a cluster moving forward in the center of mass frame.  The 
cluster disintegrates in such a way that each of its nucleons is given a 
kinetic energy 
\begin{equation}
40 < T_{nuc} < 600 \rm{MeV} 
\label{he.eq8}
\end{equation}
if the kinetic energy of the original projectile, $T_{inc}$, is 5 GeV or more.
If $T_{inc}$ is less than 5 GeV,
\begin{equation}
 40 ( T_{inc}/5 \rm{GeV} )^2 < T_{nuc} < 600 ( T_{inc}/5 \rm{GeV} )^2 .
\label{he.eq9}
\end{equation}  
In each range the energy is sampled from a distribution which is skewed 
strongly toward the high energy limit.  In addition, the angular distribution 
of the nucleons is skewed forward in order to simulate the forward motion 
of the cluster. \\

\noindent
{\bf Momentum Re-scaling.} Up to this point, all final state momenta have 
been generated in the center of mass of the incident projectile and the 
target nucleon.  However, the interaction involves more than one nucleon 
as evidenced by the intra-nuclear cascade.  A more correct center of mass 
should then be defined by the incident projectile and all of the baryons 
generated by the cascade, and the final state momenta already calculated 
must be re-scaled to reflect the new center of mass.

This is accomplished by correcting the momentum of each particle in the 
final state candidate list by the factor $T_1 / T_2$.  $T_2$ is the total 
kinetic energy in the lab frame of all the final state candidates generated
assuming a projectile-nucleon center of mass.  $T_1$ is the total kinetic 
energy in the lab frame of the same final state candidates, but whose momenta 
have been calculated by the phase space decay of an imaginary particle.  
This particle has the total CM energy of the original projectile and a 
cluster consisting of all the baryons generated from the intra-nuclear 
cascade. \\

\noindent
{\bf Corrections.} Part of the Fermi motion of the target nucleons is taken 
into account by smearing the transverse momentum components of the final 
state particles.  The Fermi momentum is first sampled from an average 
distribution and a random direction for its transverse component is chosen.  
This component, which is proportional to the number of baryons produced in 
the cascade, is then included in the final state momenta.

Each final state particle must escape the nucleus, and in the process reduce
its kinetic energy by the nuclear binding energy.  The binding energy is 
parameterized as a function of $A$:
\begin{equation}
E_B = 25\rm{MeV} \left( \frac{A-1}{120} \right) e^{-(A-1)/120)} .
\label{he.eq10}
\end{equation}

Another correction reduces the kinetic energy of final state $\pi^0$s when 
the incident particle is a $\pi^+$ or $\pi^-$.  This reduction increases 
as the log of the incident pion energy, and is done to reproduce
experimental data.  In order to conserve energy on average, the energy 
removed from the $\pi^0$s is re-distributed among the final state $\pi^+$s,
$\pi^-$s and $\pi^0$s. \\

\noindent
{\bf Nuclear De-excitation.} After the generation of initial interaction
and cascade particles, the target nucleus is left in an excited state.
De-excitation is accomplished by evaporating protons, neutrons, deuterons,
tritons and alphas from the nucleus according to a parameterized model.
The total kinetic energy given to these particles is a function of the
incident particle kinetic energy:
\begin{equation}
T_{evap} = 7.716 \rm{GeV} \left( \frac{A-1}{120} \right) F(T) e^{-F(T) - (A-1)/120} ,
\end{equation}
where
\begin{equation}
F(T) = \rm{max} [ 0.35 + 0.1304 ln(T) , 0.15 ] ,
\end{equation}
and 
\begin{eqnarray}
T = 0.1 \rm{GeV} \quad \rm{for} \quad T_{inc} < 0.1 \rm{GeV} \\
T = T_{inc} \quad \rm{for} \quad 0.1 \rm{GeV} \le T_{inc} \le 4 \rm{GeV} \\
T = 4 \rm{GeV} \quad \rm{for} \quad T_{inc} > 4 \rm{GeV} .
\end{eqnarray}
The mean energy allocated for proton and neutron emission is 
$\overline{T_{pn}}$ and that for deuteron, triton and alpha emission is 
$\overline{T_{dta}}$.  These are determined by partitioning $T_{evap}$ : 
\begin{eqnarray}
\overline{T_{pn}} = T_{evap} R \quad , \quad \overline{T_{dta}} = T_{evap} (1-R) \quad \rm{with} \nonumber
\end{eqnarray}
\begin{eqnarray}
R = \rm{max}[ 1 - (T/4\rm{GeV})^2 , 0.5 ] .
\end{eqnarray} 
The simulated values of $T_{pn}$ and $T_{dta}$ are sampled from normal 
distributions about $\overline{T_{pn}}$ and $\overline{T_{dta}}$ and their 
sum is constrained not to exceed the incident particle's kinetic energy,
$T_{inc}$.

The number of proton and neutrons emitted, $N_{pn}$, is sampled from a
Poisson distribution about a mean which depends on $R$ and the number of
baryons produced in the intranuclear cascade.  The average kinetic energy 
per emitted particle is then $T_{av} = T_{pn}/N_{pn}$.  $T_{av}$ is used
to parameterize an exponential which qualitatively describes the nuclear
level density as a function of energy.  The simulated kinetic energy of 
each evaporated proton or neutron is sampled from this exponential.  Next, 
the nuclear binding energy is subtracted and the final momentum is 
calculated assuming an isotropic angular distribution.  The number of 
protons and neutrons emitted is $(Z/A)N_{pn}$ and $(N/A)N_{pn}$, 
respectively.

A similar procedure is followed for the deuterons, tritons and
alphas.  The number of each species emitted is $0.6 N_{dta}$, $0.3 N_{dta}$
and $0.1 N_{dta}$, respectively. \\

\noindent
{\bf Tuning of the High Energy Cascade} The final stage of the high 
energy cascade method involves adjusting the momenta of the produced 
particles so that they agree better with data.  Currently, five such 
adjustments are performed, the first three of which apply only to 
charged particles incident upon light and medium nuclei at incident 
energies above $\simeq$ 65 GeV.
  
\begin{itemize}

\item If the final state particle is a nucleon or light ion with a 
momentum of less than 1.5 GeV/c, its momentum will be set to zero 
some fraction of the time.  This fraction increases with the logarithm 
of the kinetic energy of the incident particle and decreases with 
$log_{10}(A)$.

\item If the final state particle with the largest momentum happens 
to be a $\pi^0$, its momentum is exchanged with either the $\pi^+$ 
or $\pi^-$ having the largest momentum, depending on whether the 
incident particle charge is positive or negative.

\item If the number of baryons produced in the cascade is a significant 
fraction ($ > 0.3 $) of $A$, about 25\% of the nucleons and light ions 
already produced will be removed from the final particle list, provided 
their momenta are each less than 1.2 GeV/c.

\item The final state of the interaction is of course heavily 
influenced by the quantum numbers of the incident particle, particularly 
in the forward direction.  This influence is enforced by compiling, for 
each forward-going final state particle, the sum
\begin{equation}
S_{forward} = \Delta_M + \Delta_Q + \Delta_S + \Delta_B, 
\end{equation}
where each $\Delta$ corresponds to the absolute value of the 
difference of the quantum number between the incident particle and 
the final state particle.  $M$, $Q$, $S$, and $B$ refer to mass, 
charge, strangeness and baryon number, respectively.  For final 
state particles whose character is significantly different from the 
incident particle ($S$ is large), the momentum component parallel to 
the incident particle momentum is reduced.  The transverse component 
is unchanged.  As a result, large-$S$ particles are driven away from 
the axis of the hadronic shower.  For backward-going particles, a 
similar procedure is followed based on the calculation of $S_{backward}$.

\item Conservation of energy is imposed on the particles in the 
final state list in one of two ways, depending on whether or not 
a leading particle has been chosen from the list.  If all the 
particles differ significantly from the incident particle in 
momentum, mass and other quantum numbers, no leading particle is 
chosen and the kinetic energy of each particle is scaled by the 
same correction factor.  If a leading particle is chosen, its 
kinetic energy is altered to balance the total energy, while all 
the remaining particles are unaltered. 

\end{itemize}

\subsection{High Energy Cluster Production}

As in the high energy cascade model, the high energy cluster model 
randomly assigns particles from the initial collision to either a 
forward- or backward-going cluster.  Instead of performing the 
fragmentation process, however, the two clusters are treated 
kinematically as the two-body final state of the hadron-nucleon 
collision.  Each cluster is assigned a kinetic energy $T$ which is 
sampled from a distribution
\begin{equation}
exp[-aT^{1/b}]
\end{equation}
where both $a$ and $b$ decrease with the number of particles in a 
cluster.  If the combined total energy of the two clusters is 
larger than the center of mass energy, the energy of each cluster
is reduced accordingly.  The center of mass momentum of each 
cluster can then be found by sampling the 4-momentum transfer 
squared, $t$, from the distribution
\begin{equation}
 exp [t (4.0 + 1.6ln(p_{inc}) ) ]
\end{equation}
where $t < 0$ and $p_{inc}$ is the incident particle momentum.
Then,
\begin{equation}
 cos\theta = 1 + \frac{t - (E_c - E_i)^2  + (p_c - p_i)^2}{2p_c p_i},
\end{equation}
where the subscripts $c$ and $i$ refer to the cluster and incident 
particle, respectively.  Once the momentum of each cluster is 
calculated, the cluster is decomposed into its constituents.  The
momenta of the constituents are determined using a phase space decay 
algorithm.

The particles produced in the intra-nuclear cascade are grouped into 
a third cluster.  They are treated almost exactly as in the high 
energy cascade model, where Eq. \ref{he.eq1} is used to estimate the 
number of particles produced.  The main difference is that the cluster
model does not generate strange particles from the cascade.  Nucleon 
suppression is also slightly stronger, leading to relatively higher 
pion production at large incident momenta.  Kinetic energy and 
direction are assigned to the cluster as described in the cluster 
model paragraph in the previous section.  

The remaining steps to produce the final state particle list are 
the same as those in high energy cascading: 
\begin{itemize}
\item re-scaling of the momenta to reflect a center of mass which 
involves the cascade baryons,
\item corrections due to Fermi motion and binding energy,
\item reduction of final state $\pi^0$ energies,
\item nuclear de-excitation and 
\item high energy tuning.
\end{itemize}

\subsection{Medium Energy Cascading}
The medium energy cascade algorithm is very similar to the high 
energy cascade algorithm, but it may be invoked for lower incident 
energies (down to 1 GeV).  The primary difference between the two 
codes is the parameterization of the fragmentation process.  The 
medium energy cascade samples larger transverse momenta for pions 
and smaller transverse momenta for kaons and baryons.

A second difference is in the treatment of the cluster consisting 
of particles generated in the cascade.  Instead of parameterizing 
the kinetic energies and angles of the outgoing particles, the
phase space decay approach is used.  

Another difference is that the high energy tuning of the final 
state distribution is not performed.

\subsection{Medium Energy Cluster Production}
The medium energy cluster algorithm is nearly identical to the
high energy cluster algorithm, but it may be invoked for incident 
energies down to 10 MeV.  There are three significant differences 
at medium energy: less nucleon suppression, fewer particles 
generated in the intra-nuclear cascade, and no high energy tuning 
of the final state particle distributions.    

\subsection{Elastic and Quasi-elastic Scattering}
When no additional particles are produced in the initial 
interaction, either elastic or quasi-elastic scattering is 
performed.  If there is insufficient energy to induce an 
intra-nuclear cascade, but enough to excite the target nucleus, 
quasi-elastic scattering is performed.  The final state is
calculated using two-body scattering of the incident particle 
and the target nucleon, with the scattering angle in the center
of mass sampled from an exponential:
\begin{equation}
exp [ - 2 b p_{in} p_{cm} (1 - cos\theta) ] .
\end{equation} 
Here $p_{in}$ is the incident particle momentum, $p_{cm}$ is the 
momentum in the center of mass, and $b$ is a logarithmic function 
of the incident momentum in the lab frame as parameterized from 
data.  As in the cascade and cluster production models, the 
residual nucleus is then de-excited by evaporating nucleons and
light ions.

If the incident energy is too small to excite the nucleus,
elastic scattering is performed.  The angular distribution of 
the scattered particle is sampled from the sum of two 
exponentials whose parameters depend on $A$.


\section{Status of this document}

7.10.02 re-written by D.H. Wright \\
1.11.04 new section on high energy model by D.H. Wright \\