FORSCHUNGSZENTRUM KARLSRUHE
T e c h n i k u n d U m w e l t
Wissenschaftliche Berichte
FZKA 6019
CORSIKA: A Monte Carlo Code to
Simulate Extensive Air Showers
D. Heck1
, J. Knapp2
, J.N. Capdevielle3
, G. Schatz1
, T. Thouw1
1
Institut f¨ur Kernphysik, Forschungszentrum Karlsruhe, D 76021 Karlsruhe
2
Institut f¨ur Experimentelle Kernphysik, Universit¨at Karlsruhe, D 76021 Karlsruhe
3
Physique Corpusculaire et Cosmologie, Coll`ege de France, F 75231 Paris Cedex 05,
France
Forschungszentrum Karlsruhe GmbH, Karlsruhe
1998
Copyright Notice
Copyright and any other appropriate legal protection of these computer programs
and associated documentation reserved in all countries of the world.
These programs or documentation may not be reproduced by any method without
prior written consent of FORSCHUNGSZENTRUM KARLSRUHE or its delegate.
FORSCHUNGSZENTRUM KARLSRUHE welcomes comments concerning the
CORSIKA code but undertakes no obligation for maintenance of the programs, nor
responsibility for their correctness, and accepts no liability whatsoever resulting from
the use of its programs.
Trademark notice: All trademarks appearing in this report are acknowledged as
such.
Abstract
CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers
CORSIKA is a program for detailed simulation of extensive air showers initiated
by high energy cosmic ray particles. Protons, light nuclei up to iron, photons, and
many other particles may be treated as primaries. The particles are tracked through
the atmosphere until they undergo reactions with the air nuclei or - in the case of
instable secondaries - decay. The hadronic interactions at high energies may be described
by five reaction models alternatively: The VENUS, QGSJET, and DPMJET
models are based on the Gribov-Regge theory, while SIBYLL is a minijet model.
HDPM is a phenomenological generator and adjusted to experimental data wherever
possible. Hadronic interactions at lower energies are described either by the more
sophisticated GHEISHA interaction routines or the rather simple ISOBAR model. In
particle decays all decay branches down to the 1% level are taken into account. For
electromagnetic interactions the shower program EGS4 or the analytical NKG formulas
may be used. Options for the generation of Cherenkov radiation and neutrinos
exist.
Zusammenfassung
CORSIKA: Ein Monte Carlo Programm zur Luftschauersimulation
CORSIKA ist ein Programm zur detaillierten Simulation von ausgedehnten Luftschauern,
die durch hochenergetische kosmische Strahlung ausgel¨ost werden. Als
Prim¨arteilchen k¨onnen Protonen, leichte Kerne bis Eisen, Photonen und viele andere
Teilchen behandelt werden. Die Teilchen werden durch die Atmosph¨are verfolgt, bis
sie mit den Kernen der Luft reagieren oder - im Falle von instabilen Sekund¨arteilchen zerfallen.
Die hadronischen Wechselwirkungen bei hohen Energien k¨onnen wahlweise
von f¨unf verschiedenen Reaktionsmodellen beschrieben werden: Die Modelle VENUS,
QGSJET und DPMJET basieren auf der Gribov-Regge Theorie, w¨ahrend SIBYLL
ein Mini-Jet Modell ist. HDPM ist ein ph¨anomenologischer Generator und angepaßt
an experimentellen Daten, wo immer das m¨oglich ist. Hadronische Wechselwirkungen
bei niedrigeren Energien werden entweder durch die detaillierten GHEISHA Wechselwirkungsroutinen
oder durch das ziemlich einfache Isobaren-Modell beschrieben. Bei
Teilchenzerf¨allen werden alle Zerfallskan¨ale bis herab zu 1% H¨aufigkeit ber¨ucksichtigt.
Elektromagnetische Prozesse k¨onnen mit dem Schauerprogramm EGS4 oder mit den
analytischen NKG-Formeln behandelt werden. Es gibt Optionen f¨ur die Erzeugung
von Cherenkov-Strahlung und von Neutrinos.
i
ii
Contents
1 Introduction 1
2 Program frame 5
2.1 Control and run . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Coordinate system . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.4 Atmosphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.5 Random number generator . . . . . . . . . . . . . . . . . . . . . . . . 8
3 Particle tracking 9
3.1 Ionization energy loss . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.2 Coulomb multiple scattering . . . . . . . . . . . . . . . . . . . . . . . 10
3.2.1 Moli`ere scattering . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2.2 Plural Coulomb scattering . . . . . . . . . . . . . . . . . . . . 12
3.2.3 Gaussian approximation . . . . . . . . . . . . . . . . . . . . . 12
3.3 Deflection in the Earth’s magnetic field . . . . . . . . . . . . . . . . . 13
3.4 Time of flight . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.5 Longitudinal development . . . . . . . . . . . . . . . . . . . . . . . . 14
3.6 Thin sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4 Mean free path 17
4.1 Muons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.2 Nucleons and nuclei . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.2.1 Nucleon-air cross section at high energies . . . . . . . . . . . . 20
4.2.2 Nucleon-air cross section at low energies . . . . . . . . . . . . 21
4.2.3 Nucleus-nucleus cross section . . . . . . . . . . . . . . . . . . 21
4.3 Pions and kaons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.4 Other hadrons and resonances . . . . . . . . . . . . . . . . . . . . . . 25
4.5 Neutrinos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5 Hadronic interactions 29
5.1 Strong interactions at high energies . . . . . . . . . . . . . . . . . . . 29
5.1.1 VENUS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
iii
5.1.2 QGSJET . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
5.1.3 DPMJET . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.1.4 SIBYLL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.1.5 HDPM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.2 Strong interactions at low energies . . . . . . . . . . . . . . . . . . . 33
5.2.1 GHEISHA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.2.2 ISOBAR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
5.3 Nuclear fragmentation . . . . . . . . . . . . . . . . . . . . . . . . . . 35
6 Particle decays 37
6.1 πo
decays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6.2 π±
decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
6.3 Muon decay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
6.4 Kaon decays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
6.5 η decays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.6 Strange baryon decays . . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.7 Resonance decays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
7 Electromagnetic interactions 45
7.1 Muonic interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
7.1.1 Muonic bremsstrahlung . . . . . . . . . . . . . . . . . . . . . . 45
7.1.2 Muonic e+
e−
-pair production . . . . . . . . . . . . . . . . . . 46
7.2 Electromagnetic subshowers . . . . . . . . . . . . . . . . . . . . . . . 47
7.2.1 Electron gamma shower program EGS4 . . . . . . . . . . . . . 48
7.2.2 Nishimura-Kamata-Greisen option . . . . . . . . . . . . . . . . 51
7.3 Cherenkov radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
8 Outlook 57
A Atmospheric parameters 59
B Muon range for horizontal showers 63
C Default cross sections 65
C.1 Nucleon-nucleon cross sections . . . . . . . . . . . . . . . . . . . . . . 65
C.2 Pion-nucleon and kaon-nucleon cross sections . . . . . . . . . . . . . . 66
D HDPM 67
D.1 Nucleon-nucleon interactions . . . . . . . . . . . . . . . . . . . . . . . 68
D.2 Nucleon-nucleus interactions . . . . . . . . . . . . . . . . . . . . . . . 78
D.3 Pion-nucleus and kaon-nucleus interactions . . . . . . . . . . . . . . . 79
D.4 Nucleus-nucleus interactions . . . . . . . . . . . . . . . . . . . . . . . 79
Acknowledgements 81
iv
Bibliography 82
v
vi
Chapter 1
Introduction
Analyzing experimental data on Extensive Air Showers (EAS) or planning corresponding
experiments requires a detailed theoretical modelling of the cascade which
develops when a high energy primary particle enters the atmosphere. This can only
be achieved by detailed Monte Carlo calculations taking into account all knowledge
of high energy strong and electromagnetic interactions. Therefore, a number of computer
programs has been written to simulate the development of EAS in the atmosphere
and a considerable number of publications exists discussing the results of such
calculations. A common feature of all these publications is that it is difficult, if not
impossible, to ascertain in detail which assumptions have been made in the programs
for the interaction models, which approximations have been employed to reduce computer
time, how experimental data have been converted into unmeasured quantities
required in the calculations (such as nucleus-nucleus cross sections, e.g.) etc. This
is the more embarrassing, since our knowledge of high energy interactions – though
much better today than fifteen years ago – is still incomplete in important features.
This makes results from different groups difficult to compare, to say the least. In
addition, the relevant programs are of a considerable size which – as experience shows
– makes programming errors almost unavoidable, in spite of all undoubted efforts of
the authors. We therefore encourage the groups engaged in this work to make their
programs available to (and, hence, checkable by) other colleagues. This procedure
has been adopted in high energy physics and has proved to be very successful. It
is in the spirit of these remarks that we want to describe in this report the physics
underlying the CORSIKA program.
CORSIKA (COsmic Ray SImulations for KAscade) is a detailed Monte Carlo
program to study the evolution of EAS in the atmosphere initiated by photons, protons,
nuclei, or any other particle. It was originally developed to perform simulations
for the KASCADE experiment [1, 2] at Karlsruhe and has been refined over the past
years. CORSIKA meanwhile has developed into a tool that is used by many groups.
Its applications range from Cherenkov telescope experiments (E0 ≈ 1012
eV ) up to
the highest energies observed (E0 > 1020
eV ). The development of CORSIKA is
guided by the idea to predict not only correct average values of observables with this
1
EAS simulation code, but also to reproduce the correct fluctuations around the average
value. Therefore, where available, we include all known processes which might
have a noticeable influence on the observable quantities of EAS to create a reference
program that treats all processes to the present state of our knowledge. This concerns
the transport of particles through the atmosphere as well as their interactions with
air as a target. All secondary particles are tracked explicitely along their trajectories
and their parameters are stored on tape when reaching an observation level. This
allows a detailed analysis of all features of simulated showers.
The CORSIKA code has originally been developed on the basis of three well
established program systems. The first was written in the 1970s by Grieder [3]. Its
general program structure has been adopted in CORSIKA and its ISOBAR routines,
a simple hadronic interaction model, may be used as a quick alternative for the
hadronic interactions at low energies. The second part, the interaction generator
HDPM, was developed by Capdevielle [4] inspired by the Dual Parton Model [5]. It
describes the hadronic interactions of protons at high energies in good agreement
with the measured collider data. The third program deals with the simulation of
the electromagnetic part of an air shower. We incorporated the code EGS4 [6] used
successfully in the detector simulation of particle physics experiments. It was slightly
modified to the requirements of air shower simulations. These programs were merged
together and formed the first version of CORSIKA in 1989.
Basing on this primordial program many extensions and improvements have been
performed since that time. The most serious problem of EAS simulation programs is
the extrapolation of hadronic interaction to higher energies and into rapidity ranges
which are not covered by experimental data. The extreme forward direction is not
accessible by present collider experiments, rather the particles vanish undetected in
the beam pipe. But just these particles are of highest importance in the development
of EAS, as they are the most energetic secondary particles, which bring the largest
energy fraction of each collision deep into the atmosphere. Also, pp-colliders are, and
will be, limited in their maximum attainable energy to values much lower than those
found in cosmic rays. Therefore one has to rely on extrapolations based on theoretical
models. To study the systematics of such models, we have made five different
hadronic interaction models available in CORSIKA to simulate the hadronic interactions
at high energies: VENUS [7], QGSJET [8], and DPMJET [9] which describe
the inelastic hadronic interactions in the theoretically well founded manner of the
Gribov-Regge formalism, and the minijet model SIBYLL [10, 11]. These four models
give us an alternative to the phenomenological HDPM generator. Also the GHEISHA
routines [12] have been coupled which represent a more sophisticated replacement of
the ISOBAR model for the treatment of low energy hadronic interactions. With the
advent of the announced successor of VENUS [13] we plan to make it also available
within CORSIKA.
The fragmentation of nuclei in a collision may be treated in various ways, including
two options of particle evaporation from the residual nucleus. The photoproduction
of muon pairs and hadrons is incorporated into EGS4 [6]. This allows the calculation
2
of the muon content of photon induced showers. An alternative way of treating the
electromagnetic component is to use the improved and adapted form [14, 15, 16] of the
analytic NKG formula for each electron or photon produced in the hadronic cascade.
It allows to study the longitudinal development of the electromagnetic cascade and
the electron density at particular coordinates at the observation level. This option
enables a fast but less accurate simulation of hadronic showers.
Further options treat the detailed longitudinal development of various particle
species within an air shower, the production of electron and muon neutrinos and antineutrinos,
and the production of Cherenkov radiation. The elemental composition of
the atmosphere is included as well as the density variation with altitude for several
seasonal days. For nearly horizontal showers various zenith angle dependent density
profiles are provided which take into account the influence of the Earth’s curvature
in our flat atmospheric model.
With this program many calculations have been performed with p, α, O, Fe and
γ primaries in an energy range of 1011
eV ≤ E0 ≤ 1016
eV by the KASCADE group.
Also various laboratories around the world use CORSIKA to interpret and understand
their cosmic ray experiments [17, 18]. Particle numbers for electrons, muons,
and hadrons, their lateral and energy distributions, arrival times, and many other
features have been evaluated from the simulations and compared with experimental
data, where available. The agreement gives us confidence to have with CORSIKA a
useful and flexible tool to study cosmic rays and their secondaries at high primary
energies. We invite all colleagues interested in EAS simulation to propose improvements,
point out errors, or bring forward reservations concerning assumptions or
approximations which we have made.
The scope of this report is a description of the physical basis and the parametrizations
actually used in CORSIKA and to show its capabilities and limitations. The
program is in permanent modification by improvements, additions, and further details
[19, 20]. This report refers to the actual version1
and overrides the older description
[21]. Most recent changes, however, were shown to have minor effects on
the global features of simulated showers.
Additionally there exists a user’s guide [22], describing how to install the program
and how to handle input and outputs. This technical report is permanently updated
and available2
together with the complete program package from the anonymous
ftp installed at the server ftp-ik3.fzk.de . This package includes the source codes of
CORSIKA and all interaction programs, the necessary datasets, an input example,
and the user’s guide. It enables the user to setup the program with the desired
options and to run it with suitable parameters. Information on CORSIKA may be
found also in the world wide web at page http://www-ik3.fzk.de/˜heck/corsika/ .
1
Version 5.60 released in Dec. 1997
2
Requests for ftp-access should be directed to <heck@ik3.fzk.de> or <knapp@ik1.fzk.de> by
e-mail.
3
4
Chapter 2
Program frame
2.1 Control and run
At the beginning of the calculation a variety of parameters can be chosen to control
the simulation. The primary particle type has to be defined, and its energy can be
prechosen or selected at random in a particular energy range with a given slope of
the energy spectrum. This allows a realistic simulation of the shoower rate falling
steeply with rising energy. The primary angle of incidence may be defined to a
fixed value or picked at random within an angular range in a manner giving the
experimentally observable intensity of an equal particle flux from all directions of the
sky penetrating through a horizontal detector area [23]. The atmospheric parameters
may be selected to study the influence of the seasons. Up to 10 observation levels
can be defined and data on all particles penetrating these levels are recorded as long
as the energy exceeds a cut-off, specified for hadrons, muons, electrons, or photons
separately. Several flags select and control the hadronic interaction models at high
energies and the respective cross sections, one flag selects the low energy hadronic
interaction model. Various possibilities exist to simulate the fragmentation of the
primary nucleus, and two flags switch on or off the two options for the simulation of
electromagnetic cascades. Using the ‘thin sampling’ option the thinning level may
be specified. All these controls, also on the printing of various lists and tables, are
described extensively in Ref. [22].
2.2 Particles
The CORSIKA program recognizes 50 elementary particles. These are γ, e±
, µ±
,
πo
, π±
, K±
, Ko
S/L, η, the baryons p, n, Λ, Σ±
, Σo
, Ξo
, Ξ−
, Ω−
, the corresponding
anti-baryons, the resonance states ρ±
, ρo
, K∗±
, K∗o
, K
∗o
, ∆++
, ∆+
, ∆o
, ∆−
, and the
corresponding anti-baryonic resonances. Optionally the neutrinos νe and νµ and antineutrinos
νe and νµ resulting from π, K, and µ decay may be generated explicitly. In
addition nuclei up to A = 56 can be treated. Within the program they are identified
5
by their numbers of protons and neutrons. All these particles can be tracked through
the atmosphere, are able to interact, annihilate or decay, and produce secondary
particles. They are fully defined in the program by their particle identification, the
Lorentz factor, the zenith and azimuthal angle of the trajectory, the time since the
first interaction of the primary, and the three spatial coordinates x, y, and z. The
number of inelastic hadronic reactions or decays which the parent particles have
suffered is protocoled as the generation of a particle. The particle masses and charge
states are stored in an array for fast access during the calculations. Particle identifiers
and masses of elementary particles are taken from the GEANT3 detector simulation
code [24] with extensions for 4 neutrino species and resonance states; the resonance
masses correspond with Ref. [25]. The nuclear masses are taken as the sum of the
constituent nucleon masses neglecting binding energies. If nuclear binding effects
should be regarded, we provide as an alternative for nuclei with Z < 15 the isotopic
masses from the mass table of Ref. [26] with corrections for the electron masses and
for other nuclear masses a calculation according to the Bethe-Weizs¨acker formula.
Projectile nuclei are assumed to be completely stripped, i.e. their charge state q is
set equal to their atomic number Z.
2.3 Coordinate system
The coordinates in CORSIKA are defined with respect to a Cartesian coordinate
system with the positive x-axis pointing to the magnetic north, the positive y-axis to
the west, and the z-axis upwards. The origin is located at sea level. This definition
is necessary, because the Earth’s magnetic field is taken into account. By default it
is implemented for the location of Karlsruhe (49o
N, 8o
E) as described in section 3.3.
The zenith angle θ of a particle trajectory is measured between the particle momentum
vector and the negative z-axis, and the azimuthal angle φ between the positive
x-axis and the x-y-component of the particle momentum vector (i.e. with respect to
north) proceeding counterclockwise.
2.4 Atmosphere
The atmosphere adopted consists of N2, O2, and Ar with the volume fractions of
78.1%, 21.0%, and 0.9% [27]. The density variation of the atmosphere with altitude
is modeled by 5 layers. The boundary of the atmosphere in this model is defined at the
height In the lower four of them the density follows an exponential dependence on the
altitude leading to a relation between the mass overburden T(h) of the atmosphere
and the height h of the form
T(h) = ai + bie−h/ci
i = 1, . . . , 4 . (2.1)
6
-25
-20
-15
-10
-5
0
5
10
15
20
0 5 10 15 20 25 30 35
height (km)
pressuredifference(hPa)
AT115
AT223
AT511
AT616
AT822
AT1014
AT1224
Figure 2.1: Pressure difference of the atmosphere above Stuttgart (Germany) at 7
days distributed over the year 1993 relativ to U.S. standard atmosphere.
In the fifth layer the mass overburden decreases linearly with height
T(h) = a5 − b5
h
c5
. (2.2)
where the mass overburden T(h) vanishes, which is at h = 112.8 km. The parameters
ai, bi, and ci are selected in a manner that the function T(h) is continuous at the
layer boundaries and can be differentiated continuously. Various atmospheres are
foreseen: U.S. standard atmosphere parametrized according to J. Linsley [28], 7 typical
atmospheres as measured above Stuttgart (about 60 km away from Karlsruhe)
at various days of 1993, transmitted by Deutscher Wetterdienst Offenbach and parametrized
according Ref. [29]. Out of a large ensemble of measured atmospheres these
7 sets have been selected such that characteristic seasonal differences show up. The
parameter values of the available atmospheres are listed in Tables A.1 to A.8. Fig.
2.1 shows the pressure difference of the 7 atmospheres relative to the U.S. standard
atmosphere. Therefore also the pressure at ground level varies from parameter set
to parameter set, as listed in Table A.9.
In CORSIKA always a flat atmosphere is adopted. In the simulation of nearly
horizontal showers with θ ≥ 75o
the influence of the curvature of the Earth’s surface
7
is no longer negligible. To avoid the lengthy calculations in a spherical coordinate
system, optionally the analytical description of the atmosphere may be replaced by
tabulated mass overburden distributions, calculated for that angle. Thus we simulate
those nearly horizontal showers with θ = 0, but with an atmospheric profile which is
present along the axis of a nearly horizontal shower.
The passage of the primary particle through the atmosphere starts at the upper
border of the atmospheric model. From this starting point the place of the first
interaction is calculated. The height and the target nucleus of this interaction are
selected at random. Optionally both selections may be fixed by input values. The
coordinates of the point of first interaction are set to (0, 0, z0). At each observation
level the x and y coordinates are shifted such that the shower axis retains the
coordinates (0, 0, zobs). This is done to facilitate later analysis.
2.5 Random number generator
The Monte Carlo method is essentially based on random numbers and, hence, a
random number generator that meets the requirements of the today’s increasingly
long and complex calculations is indispensable. CORSIKA is operated with the
random number generator RANMAR [30] in the version as implemented in the CERN
program library [31] which represents the state of the art in computational physics.
It is a pseudo random number generator delivering uniformly distributed numbers. It
offers the possibility to generate simultaneously up to 9 · 108
independent sequences
with a sequence length of ≈ 2 · 1044
each. The generator is written in standard
FORTRAN and, thus, portable to all types of computers where bit-identical results
are obtained. It satisfies very stringent tests on randomness and uniformity and it is
sufficiently fast.
8
Chapter 3
Particle tracking
For propagating particles between two interaction points their space and time coordinates
as well as their energy have to be updated. For electrons and photons this
is done in EGS4 as described in Ref. [6] and section 7.2.1. Charged particles lose
energy by ionization whereas neutral particles proceed without energy loss. Because
of the large penetration depth of µ±
a deflection due to multiple Coulomb scattering
is taken into account. This is neglected for charged hadrons. All charged particle
trajectories are bent in the Earth’s magnetic field. The time update is handled for
all particles in the same straightforward way. If particles cross an observation level
while being tracked to the next interaction point, their space, momentum, and time
coordinates are computed for the observation level and transferred to the particle
output file.
3.1 Ionization energy loss
The energy loss by ionization of a charged particle which traverses matter of thickness
λ is described by the Bethe-Bloch stopping power formula
dEi =
λ z2
β2
κ1 ln(γ2
− 1) − β2
+ κ2 =
λγ2
z2
γ2 − 1
κ1 ln(γ2
− 1) − β2
+ κ2
where β = v/c is the velocity of the particle in the laboratory in units of the velocity
of light, γ is its Lorentz factor, z is the charge of the ionizing particle in units of e.
The two constants κ1 = 0.153287 MeV g−1
cm2
and κ2 = 9.386417 are derived from
the tables [32] for dry air. This expression is used to compute the ionization energy
loss along the particle trajectory. High energy muons with a Lorentz factor γ > 2·104
suffer from an additional energy loss by bremsstrahlung (see section 7.1.1) and direct
e+
e−
pair formation (see section 7.1.2). The energy loss of muons as function of their
energy is given in Figure 3.1. Whenever, after updating the energy, the corresponding
Lorentz factor is below the cut-off value, the particle is dropped from the calculation.
9
1
10
10 2
1 10 10
2
10
3
10
4
10
5
10
6
γ
dE/dx[MeVg
-1
cm
2
]
Figure 3.1: Energy loss of muons in air as function of Lorentz factor. The contributions
from ionization (dashed line) and direct pair production (dotted line) are
indicated.
3.2 Coulomb multiple scattering
Charged particles are scattered predominantly in the electric Coulomb field of the
nuclei in the (traversed medium) air. As these nuclei are generally much more massive
than the scattered particles, the direction of flight might be altered, but not the
energy. In CORSIKA the process of Coulomb multiple scattering is considered only
for muons and only once for each tracking step in the middle of the tracking distance.
The angular distribution of the multiple scattering is described by Moli`ere’s theory
[33] or may alternatively (selectable) may be approximated by a Gaussian function.
Only about 2% of the events with large scattering angle occur more frequently than
predicted by a Gaussian. The procedure to select the scattering angle is taken from
the detector simulation code GEANT3 [24]. For heavy particles at high energies
multiple scattering is not important. The multiple scattering of electrons is treated
very detailed in EGS4 according to Moli`ere’s theory.
10
3.2.1 Moli`ere scattering
The determination of the scattering angle follows Ref. [24]. First, the number Ω0 of
scatterings along the traversed matter of thickness λ is calculated according to
Ω0 = 6702.33
λ
β2
Zs
mair
e(Ze−Zx)/Zs
where β is the velocity of the muon in units of the speed of light and mair = 14.54
is the average atomic weight of air. The quantities Zs, Ze, and Zx depend on the
atomic fractions ni of atoms of type i with charge number Zi in air:
Zs =
3
i=1
niZi(Zi + 1)
Ze =
3
i=1
niZi(Zi + 1) ln Z
−2/3
i
Zx =
3
i=1
niZi(Zi + 1) ln (1 + 3.34(Ziα)2
) .
Here α is the fine structure constant. The index i represents the 3 components of air
(see section 2.4).
If the number of scattering events is low, i.e. Ω0 ≤ 20, the total scattering angle
is taken as geometrical sum of individual scatterings, which are simulated according
to section 3.2.2. Otherwise the polar scattering angle θ is sampled from
f(θ) θ dθ =
sin θ
θ
fr(η) dη .
Here the first three terms of Bethe’s expansion [34] are used for the function fr(η)
fr(η) = f(0)
r (η) + f(1)
r (η)B−1
+ f(2)
r (η)B−2
.
Tabulated values of the three functions f(k)
r are contained in CORSIKA for the range
0 ≤ η ≤ 13 of the reduced angle η which is defined by
η =
θ
χc
√
B
.
The quantity B is evaluated from
B − ln B = ln Ω0 ,
and the critical angle χc is defined by
χc =
0.00039612
√
Zs
β2E
√
mair
√
λ
11
with E the muon energy. The actual value of fr is derived by a four point continued
fraction interpolation between the tabulated values. Scattering angles with θ > π
are rejected.
The radial deviation from the straight trajectory is computed and the azimuthal
angle is selected at random from a uniform distribution. Ther particle path then
approximated by two stright lines following the incident direction to the midpoint of
the tracking step and the new direction thereafter.
3.2.2 Plural Coulomb scattering
If the number of scatters is low, Moli`ere’s theory is not applicable and we assume a
Poisson distributed number j of scattering events around Ω0. Assuming a Rutherford
cross section σ for single elastic scattering in the material with charge Z we have
f(θ) θ dθ =
dσ
θ dθ
1
σ
θ dθ with
dσ
θ dθ
= 2π
2 Z e2
E β2 c
2
1
(θ2 + χ2
α)
with the screening angle χα which is calculated from
χ2
α =
0.00039612
√
Zs mair
1.167E2 β2 6702.33 Zs e(Ze−Zx)/Zs
.
For each individual scattering we sample θj from
f(θj) θj dθj = χ2
α
2θj dθj
(θ2
j + χ2
α)2
,
which leads, with the random number RD, to
θj = χ2
α
1
RD
− 1
according to Ref. [24]. To get the total polar scattering angle we take the azimuth
angles φj at random from a uniform distribution and add up the projections of θj
onto the x and y plane. With this total polar scattering angle the further calculation
proceeds as described in the previous section 3.2.1.
3.2.3 Gaussian approximation
In the Gaussian approximation the mean square value of the polar scattering angle
θ of a given path is calculated according to the expression of Ref. [35]
θ2
= λθ2
s with θ2
s =
1
λs
Es
mγβ2
2
.
12
Here Es = 0.021 GeV is the scattering constant, m, γ, and β represent the mass,
Lorentz factor and velocity in the laboratory, respectively. λ is the amount of matter
traversed by the particle and λs = 37.7 g/cm2
. The value for θ is picked at random
according to the Gaussian distribution
P(θ, λ) =
1
πθ2
s λ
e−θ2/(λ θ2
s)
.
The azimuthal angle is selected at random from a uniform distribution. The computation
of the arrival coordinates x and y is performed analogously to the Moli`ere
case.
3.3 Deflection in the Earth’s magnetic field
The Earth’s magnetic field is characterized by its strength BE
, its declination angle
δ, and its inclination angle ϑ. For the KASCADE location these values are at present
BE
= 47.80 µT δ = −9 and ϑ = 64o
44
corresponding with the components
Bx = 20.40 µT and Bz = −43.23 µT ,
while By = 0 by definition (see section 2.3). Because of the small value of δ the
deviation of the x-direction from the geographic north may be neglected. Values
for other locations may be obtained from the program Geomag, which is available
on-line in the world wide web [36].
A particle with charge Z and momentum p travelling along the path length in
the magnetic field B suffers a deflection which points to the direction normal to the
plane spanned by B and p. The direction is changed by the angle α which, for small
deflection angles, is approximately given by
α ≈ Z
p × B
p2
.
3.4 Time of flight
At the first interaction of the primary in the atmosphere the timing of the shower is
started. The time interval dt which elapses as the particle moves along its path is
computed by dividing the particular path length by the average particle velocity
βave. Thus,
dt =
cβave
where βave is the arithmetic mean of the laboratory velocities of the particle at
beginning and end of the trajectory.
The total time elapsed since the first interaction is the sum of all time intervals
accumulated by the successive particles to the particular observation level.
13
3.5 Longitudinal development
At a selectable step width (in g/cm2
) the longitudinal development of EAS may be
followed by sampling the number of gammas, positrons, electrons, negative and positive
muons, hadrons, all charged particles, nuclei, and Cherenkov photons traversing
each of these sampling altitudes. A function of the type
N(T) = Nmax
T − T0
Tmax − T0
Tmax−T
a+bT +cT 2
may be fitted to the ‘all charged’ distribution to describe the dependence on the
atmospheric depth T. The resulting 6 parameters Nmax, T0, Tmax, a, b, and c and
the χ2
/dof are stored. If EGS4 is not selected, such a fit is tried for the levels which
are used to determine the NKG longitudinal distribution (see section 7.2.2).
3.6 Thin sampling
In Monte Carlo programs for EAS simulation the computing times roughly scale with
the primary energy, thus becoming excessive for showers initiated by particles with
E0 > 1016
eV . A way out of this problem is the ‘thin sampling’ [37] or ‘variance
reduction’ [6]. All secondary particles below an adjustable fraction of the primary
energy (thinning level εth = E/E0) are exposed to this procedure. If the energy sum
of all j secondary particles 1
falls below the thinning energy
εthE0 >
j
Ej
only one of the secondary particles is followed, selected at random according to its
energy Ei with the probability
pi = Ei/
j
Ej
while all other secondaries are discarded. An appropriate weight wi = 1/pi is attributed
to the surviving particle, thus conserving energy. More correctly the initial
weight of the particle is multiplied by wi.
If only in part the energy of secondary particles falls below the thinning level, the
corresponding particles survive with a probability pi given by
pi = Ei/(εthE0)
and, in case of surviving, get the weight factor wi = 1/pi. The latter procedure is
also applied if the energy sum of the corresponding particles exceeds the thinning
level, thus enabling more than one particle to survive.
1
Particles not taken into account like neutrinos or falling below the cut-off are excluded from
the sum.
14
εth none 10−6
10−5
10−4
10−3
Time (min) 98 51 7.2 1.2 0.16
Table 3.1: Computing times for various thinning levels.
By this selection mechanism only a rather constant number of particles must
be followed in the low energy portion of EAS instead of an exponentially growing
bulk. This mechanism is implemented in the hadronic part as well as in the EGS4
routines of Ref. [6]. Table 3.6 shows the dependence of the computing time on εth for
E0 = 1015
eV proton showers at vertical incidence employing the QGSJET hadronic
interaction model and the EGS4 routines at default CORSIKA parameters.
15
16
Chapter 4
Mean free path
The distance a particle travels before it undergoes its next inelastic interaction or
decay is determined by the cross section for a hadronic reaction together with the
atmospheric density distribution along the flight path, and the probability to decay.
Stable particles can only interact, for unstable ones the two processes compete. A
decay length and an interaction length are determined independently at random and
the shorter one is taken as the actual path length. By this procedure it is also decided
whether a particle decays or interacts.
Electrons and photons are treated in the EGS4 routines. A description of the
processes they may suffer can be found in [6]. Additional information is given in
section 7.2.1.
4.1 Muons
As inelastic hadronic interactions of muons are very rare they are omitted in CORSIKA.
Therefore muons may only decay or undergo bremsstrahlung and e+
e−
-pair
production interactions. The interaction cross sections σint for bremsstrahlung and
e+
e−
-pair production are calculated in Ref. [38]. We use the parametrizations from
the detector simulation code GEANT3 [24]. The mean free path λint for these interactions
is given by
λint = mair/σint
where mair = 14.54 is the average atomic weight of air in g/mol and λint is given in
g/cm2
. The probability of the muon to traverse an atmospheric layer of thickness λ
without corresponding interaction is then
Pint(λ) =
1
λint
e−λ/λint
. (4.1)
Relativistic muons may propagate through a large fraction of the atmosphere. As
we use an atmosphere composed from different layers, a transportation step may only
reach at maximum to the next layer boundary or a detection level. Additionally, the
17
transportation step is limited to 10 · λs = 377 g/cm2
to treat the multiple scattering
correctly (see section 3.2.3).
If decay determines the path length , its mean is defined by
D
= cτµγµβµ (4.2)
where c is the vacuum speed of light, τµ is the muon life time at rest, γµ is its Lorentz
factor and βµ its velocity in units of c. The probability of a muon to travel the
distance before it decays is then
PD
( ) =
1
D
e− / D
(4.3)
and the path length of a muon may be chosen at random from this distribution. It
should be noted that in the above formulas has the dimension of cm. The path
lengths are expressed in units of g/cm2
by taking into account the actual density and
its variation along the trajectory. For a given path length in cm one obtains the
path length λ in g/cm2
as
λ = f( , h0, θ) =
T(h) − T(h0)
cos θ
with h = h0 − cos θ. It depends on the altitude of its origin h0 and the direction
of propagation. The mass overburden T(h) is given by Equation 2.1 or 2.2. The
probability distribution for the decay distance λ in g/cm2
is then
PD
(λ) = PD
( )
d
dλ
= P(f−1
(λ))
df−1
(λ)
dλ
where f−1
represents the inverse function of f. This consideration is only valid, if
the kinetic energy of the muons does not change during transport.
In reality in Equation 4.2 γµ and βµ are functions of itself and one obtains a
shorter decay length. In good approximation [39] we determine in a manner that
0
d
β( )γ( )
= −cτµ ln(RD) (4.4)
where the term at the right side determines at what time in the muon rest frame the
decay happens. The Lorentz factor γ( ) results from
γ( ) = γ0 +
dE
dx
(γ0)
T(h) − T(h0)
mµ cos θ
with dE
dx
(γ0) the ionization loss per radiation length, T(h0) the mass overburden of
the atmosphere and γ0 the Lorentz factor of the muon at the starting point. T(h)
gives the thickness in the altitude h = h0 − cos θ, mµ is the muon mass, and θ the
18
zenith angle of the muon. Neglecting the change of dE/dx along a transport step
and assuming β ≈ 1 the muon range can be calculated analytically in an atmosphere
with an exponential density profile1
. With this assumption Equation 4.4 results in
−cτµ ln(RD) =
0
d
βγ
≈
0
d
γ( )
=
ci
di cos θ
ln
γ0(di − γ)
γ(di − γ0)
(4.5)
where di is defined by
di = γ0 −
dE
dx
(γ0)
T(h0) − ai
mµ cos θ
= γ0 + i .
ai and ci are the atmospheric parameters for layer i defined by Equation 2.1. Solving
Equation 4.5 for the Lorentz factor γ delivers
γ =
γ0di
γ0 + i exp(cτµ ln(RD)di cos θ/ci)
. (4.6)
At the bottom boundary hi of the atmospheric layer i the Lorentz factor is reduced
by ionization energy loss (using the Bethe-Bloch stopping power formula Equation
3.1 with the constants κ1 and κ2) to
γi = γ0 −
T(hi) − T(h0)
mµ cos θ
γ2
0
γ2
0 − 1
κ1 ln(γ2
0 − 1) − β2
+ κ2 . (4.7)
If the Lorentz factor of Equation 4.6 exceeds that of Equation 4.7, the muon decays
within the layer i, otherwise for the layer i − 1 below that boundary the quantity
cτµ ln(RD) is replaced by
[cτµ ln(RD)]i−1 = [cτµ ln(RD)]i −
h0 − hi + ci ln(γ0/γi)
di cos θ
(4.8)
and starting with the parameters at the lower boundary of layer i the muon is followed
down through the next layer beneath until it decays, undergoes a pair formation
or bremsstrahlung event, or reaches the observation level. As the interactions are
stochastic processes, they are treated separately in sections 7.1.1 and 7.1.2.
4.2 Nucleons and nuclei
The slowing down of the hadronic projectiles along their path through the atmosphere
is taken into account. Because of the reduced energy this leads to a modified cross
section and hence to a modified free path. In contrast to the muon decay length
this effect is omitted for hadronic interactions. At energies above ≈ 100 GeV
1
For other density profiles see Appendix B.
19
the cross sections decrease with decreasing energy and the ionization energy losses
would result in a slightly increased range of the hadrons. As the cross sections vary
very moderately with energy and the ionization energy losses are small compared
to the kinetic energy of the hadrons, it is justified to neglect the ionization energy
losses. The free paths of nucleons or nuclei as stable particles are determined by their
inelastic nucleon-air or nucleus-air cross section only.
4.2.1 Nucleon-air cross section at high energies
Depending on the employed hadronic interaction model different parametrizations of
the cross sections are used. In all cases the cross section is taken as the weighted
sum over the cross sections σn−Ni
of the components of air
σn−air =
3
i=1
niσn−Ni
(4.9)
250
300
350
400
450
500
550
600
10
2
10
3
10
4
10
5
10
6
10
7
10
8
10
9
VENUS
QGSJET
SIBYLL
HDPM
DPMJET
plab (GeV/c)
σinel
p-air
(mb)
Mielke et al.
Yodh et al.
Gaisser et al.
Baltrusaitis et al.
Honda et al.
Figure 4.1: Inelastic proton-air cross sections for the models used in CORSIKA as
a function of projectile momentum. The experimental data are taken from Refs.
[40, 41, 42, 43, 44]. The shaded band represents the result of a fit of the form
σp−air
inel = a log(p) + b log2
(p) + c to the data with p < 105
GeV/c. For HDPM see
footnote in Appendix C.
20
with ni the atomic fraction of component i. The available nucleon-air cross sections
are shown in Fig. 4.1 together with experimental data. They include the diffractive
contribution as given by each of the different models. If not specified, the default
cross sections (HDPM) as described in Appendix C are used. The interaction mean
free path λint
is obtained from the nucleon-air cross sections by
λint
= mair/σn−air (4.10)
where mair = 14.54 is the average atomic weight of air in g/mol and λint
is given
here in g/cm2
. The probability of the projectile to traverse an atmospheric layer of
the thickness λ without interaction is then
Pint
(λ) =
1
λint
e−λ/λint . (4.11)
From this distribution, the path lengths of nucleons are chosen at random. The
interacting target nucleus is randomly selected according to the contribution of each
air component to the sum of Equation 4.9.
4.2.2 Nucleon-air cross section at low energies
Using the GHEISHA package below Elab ≤ 80 GeV , the interaction cross sections are
interpolated from tables supplied with GHEISHA [12]. They comprise elastic and
inelastic interactions and include slow neutron capture processes in a rather realistic
manner as derived from experimental data. They are plotted in Fig. 4.2.
Employing the ISOBAR model, the inelastic cross sections as given in Appendix
C are used. The measured inelastic nucleon-nucleon cross section drops rapidly with
decreasing energy. Therefore below plab = 10 GeV/c the ISOBAR model allows only
elastic reactions with a constant cross section. For antinucleons an annihilation with
nucleons can occur in addition, leading to a contribution to the inelastic cross section
which rises with decreasing energy. These cross sections are also shown in Fig. 4.2.
4.2.3 Nucleus-nucleus cross section
In EAS nucleons or complex nuclei are reacting with air nuclei. In the energy range of
interest no experimental data exist on the relevant quantities, such as inelastic cross
section and number of target and projectile nucleons participating in the reaction. So
they have to be calculated from the nucleon-nucleon cross section following Glauber
theory [45, 46]. The input nucleon distributions of nuclei are derived from measured
charge distributions [47] unfolding the finite size of the proton with a mean square
charge radius of r2
p
1/2
= 0.862 fm. For nuclei with mass number below 20 the charge
distributions are assumed to be Gaussian and the radius of the nucleon distribution
is
rm
2
= rch
2
− rp
2
.
21
300
400
500
600
700
800
900
1000
2000
10
-1
1 10 10
2
10
3
plab (GeV/c)
σinel
n-air(mb)
p, n
p
–
, n
–
p
n
p
–
n
–
GHEISHA
ISOBAR
HDPM
Figure 4.2: Inelastic nucleon-air cross sections at lower energies. The cross sections
of the GHEISHA model are drawn for protons, anti-protons (dashed lines), neutrons,
and anti-neutrons (dashed-dotted lines) together with the cross sections of the ISOBAR
and HDPM model for nucleons (solid line) and anti-nucleons (dotted line).
The dotted vertical line marks the boundary between high energy and low energy
hadronic interaction models.
For nuclei with A > 20 the charge distributions are approximated by the Fermi
function. Unfolding was done by folding a correspondingly parametrized nucleon
distribution with the proton charge distribution to obtain the measured radius and
diffuseness of the charge distribution of the nucleus. From the Glauber method the
inelastic cross sections for all projectile nuclei in the stability valley with A = 1 . . . 56
colliding with the target nuclei 14
N, 16
O, and 40
Ar were calculated for three different
values of the nucleon-nucleon cross section (30, 45, and 60 mb; they correspond
to nucleon-nucleon collisions at laboratory energies of 120 GeV , 66.5 TeV ,
and 5.87 PeV , respectively, for the HDPM generator). Values for mass numbers for
which no experimental charge distributions were available have been interpolated.
The results are tabulated in CORSIKA and a quadratic interpolation is performed
between the table values with respect to σn−n to obtain intermediate values of the
cross section σN−Ni
of a nucleus with component i of air. Then σN−air is obtained
22
20
30
40
50
60
70
80
90
100
10
2
10
3
10
4
10
5
10
6
10
7
10
8
10
9
10
10
VENUS
QGSJET
SIBYLL
HDPM
DPMJET
plab (GeV/c)
σinel
h-p(mb)
p
π
K
Figure 4.3: Inelastic cross sections of protons, pions, and kaons with protons for
the models used in CORSIKA as function of projectile momentum. For HDPM see
footnote in Appendix C.
from the weighted sum over the three components of air
σN−air =
3
i=1
niσN−Ni
. (4.12)
The SIBYLL model provides its own nucleus-nucleus cross section table including
an interpolation routine, while for the other models we use the model-specific nucleonnucleon
cross section and apply our Glauber formalism to get the desired nucleusnucleus
cross section. The used proton-proton cross sections are shown in Fig. 4.3.
The resulting nucleus-air cross sections are shown in Figure 4.4. For DPMJET and
VENUS it is possible to calculate model specific nucleus-nucleus cross sections in a
time consuming manner. They differ from the cross sections shown in Fig. 4.4 by
< 8 %.
The interaction mean free path λint
is obtained by
λint
= mair/σN−air
where mair = 14.54 is the average atomic weight of air in g/mol and λint
is given
here in g/cm2
. This is done for high energy interaction models. In GHEISHA and
23
300
400
500
600
700
800
900
1000
2000
10 10
2
10
3
10
4
10
5
10
6
10
7
10
8
10
9
10
10
VENUS
QGSJET
SIBYLL
HDPM
DPMJET
plab (GeV/c)
σinel
A-air(mb)
p
He
O
Fe
Figure 4.4: Inelastic nucleus-air cross sections of various projectile nuclei for the
models used in CORSIKA as function of projectile momentum. For HDPM see
footnote in Appendix C.
ISOBAR only nucleon-nucleon interactions are treated. The path lengths of projectile
nuclei are sampled from an appropriate distribution. The interacting target nucleus
is randomly selected according to the contribution of each air component to the sum
of Equation 4.12.
In the treatment of hadron-nucleus and nucleus-nucleus collisions by the various
models, the number of interacting target and projectile nucleons is determined by the
selected model. SIBYLL and HDPM treat nucleus-nucleus collision as a superposition
of nucleon-nucleus collisions.
4.3 Pions and kaons
Charged pions and all four kinds of kaons are particles where decay and nuclear
interaction compete. Their decay lengths are determined in the same way as for
muons just replacing the free muon life time τµ by the pion and kaon life times τπ
and τK, respectively. For pions and charged kaons the ionization energy loss is also
taken into account in the same manner as for muons.
The interaction lengths are treated in analogy to those of nucleons. They are
24
200
300
400
500
600
10 10
2
10
3
10
4
10
5
10
6
10
7
10
8
10
9
VENUS
QGSJET
SIBYLL
HDPM
DPMJET
plab (GeV/c)
σinel
h-air(mb)
p
π
K
Figure 4.5: Inelastic hadron-air cross sections at higher energies of various interaction
models as used in CORSIKA. In SIBYLL and DPMJET the pion cross sections are
also used for kaons. For HDPM see footnote in Appendix C.
determined according to Equations 4.10 and 4.11 using the model specific mesonair
cross sections as shown in Fig. 4.5. Also the selection of the interacting target
nucleus follows section 4.2.1.
Employing the GHEISHA package below Elab ≤ 80 GeV the interaction cross
sections are interpolated from tables supplied with GHEISHA. The resulting pionair
and kaon-air cross sections are shown in Figure 4.6. The actual free path of the
meson is taken as the minimum of a random decay length and a random interaction
length. By this selection also the decision between decay and interaction is made.
Due to their short life time of neutral pions of τ = 8.3·10−17
sec their probability
for interaction is omitted for Elab < 1014
eV and they decay at their origin. Above
they are tracked by analogy with the charged pions. Above Elab > 1018
eV their
decay length becomes comparable to their interaction length in the lower atmosphere.
4.4 Other hadrons and resonances
The η meson is treated in the same manner as neutral π mesons. Because of their
short life time of order of 10−23
sec the ρ, K∗
, and ∆ resonance states decay imme-
25
200
300
400
σinel
π-air(mb)
GHEISHA
ISOBAR
HDPM
π
±
π
+ π
–
10
2
10 3
10
-1
1 10 10
2
10
3
plab (GeV/c)
σinel
K-air(mb)
K
K
–
K
0
L
K
+K
0
S
a)
b)
Figure 4.6: Inelastic meson-air cross sections at lower energies for a) pions and b)
kaons. The solid lines give the cross section as calculated for the ISOBAR and HDPM
model, the dashed and dotted lines lines with signatures give the cross sections as
used in the GHEISHA option. The dotted vertical line marks the boundary between
high energy and low energy hadronic interaction models.
diately after their production without being tracked.
Baryons with strangeness ±1, ±2, and ±3 are produced by most models explicitly,
while HDPM knows only hadrons with one strange quark or anti-quark. In spite of
their short life time of ≈ 10−10
sec their mean free path is too long to be neglected.
Their decay length is determined analogously to that of π±
and kaons, taking into
account the ionization energy loss for the charged strange baryons. Their interaction
cross sections are adopted to be the same as for nucleons. In the low energy range
only the GHEISHA routines are able to treat strange baryons, the corresponding
26
cross sections are interpolated from the GHEISHA cross section tables.
4.5 Neutrinos
As all kinds of neutrinos have extremely small interaction cross sections, they are
assumed to traverse the atmosphere without interaction.
27
28
Chapter 5
Hadronic interactions
Hadronic interactions are simulated within CORSIKA by several models depending
on energy. If the energy is high enough, the interaction is treated alternatively
with one of the models VENUS [7], QGSJET [8], DPMJET [9], SIBYLL [10], or
HDPM. The former are well documented and the interested reader is refered to the
literature. The last one is described in detail in Appendix D. The high energy models
reach their limit, if the cm energy available for generation of secondary particles
drops below a certain value. This value is presently set at 12 GeV corresponding
to a laboratory energy of 80 GeV for the transition to GHEISHA. If GHEISHA is
replaced by the ISOBAR model the transition energy is lowered to Ecm = 10 GeV
rsp. Elab = 50 GeV . In an EAS, however, the bulk of particles interacts at cm
energies far below these values.
Below these transition energies the user may select between the GHEISHA routines
[12] or the ISOBAR model. The GHEISHA routines as implemented in CORSIKA
are taken from the detector simulation code GEANT3 [24]. This detector
simulation code is used by many high energy experimental groups, and therefore
much experience on the validity of the GHEISHA routines exists. The simple ISOBAR
model enables fast calculations and is used for the hadronic interactions as explained
in section 5.2.2. It should be noted that the selection of DPMJET, QGSJET,
SIBYLL, or VENUS needs (and automatically forces) the GHEISHA option, as the
baryons with strangeness ± 2 and ± 3 generated by those models cannot be treated
by the ISOBAR model.
5.1 Strong interactions at high energies
Table 5.1 gives an overview on essential features of the available high energy hadronic
interaction models. A comparison of the interaction models is given in Refs. [48, 49,
50] and the effects on observables of EAS are discussed there in detail.
29
VENUS QGSJET DPMJET SIBYLL HDPM
Gribov-Regge + + + +
Minijets + + + +
Sec. Interactions +
N-N Interaction + + +
Superposition + +
Max. Energy (GeV) 2 × 107
> 1011
> 1011
> 1011
108
Table 5.1: Basic features of the interaction models used.
5.1.1 VENUS
The VENUS code [7] (Very Energetic NUclear Scattering) is mainly designed to
treat nucleon-nucleon, nucleon-nucleus, and nucleus-nucleus scattering at ultrarelativistic
energies. It is based essentially on the Gribov-Regge theory, which considers single
or multiple Pomeron exchange as the basic process in high energy hadron-hadron
scattering. A Pomeron is represented by a cylinder of gluons and quark loops according
to the topological expansion of quantum chromodynamics. Particle production
in inelastic collisions amounts to cutting Pomerons, i.e. cutting cylinders, which in
VENUS is realized by forming colour strings which subsequently fragment into colour
neutral hadrons. At high densities, for collisions of heavy nuclei, when independent
binary interactions become unlikely, massive quark-matter droplets may be formed.
Final state interactions are taken into account. Diffractive and non-diffractive collisions
as well as mesonic projectiles are described with the same formalism. This
leads to a consistent and theoretically well founded treatment of all types of hadronic
reactions involved in an air shower cascade. As jets which become important at extremely
high energies are not contained within this model, we recommend an upper
limit of Elab < 2 · 1016
eV for the projectile.
All mesons and baryons known to the CORSIKA program (see section 2.2) may
be used as projectiles in VENUS. Moreover also nuclei are admitted. Neutral Ko
SL
mesons are regarded as K
o
or K
o
particles with equal probabilities. Photons undergoing
hadronic reactions are treated as πo
or η with equal chance (see section
7.2.1). The parameter wproj which enters into the dermination of the probability of
diffractive interaction of the projectile, is set to 0.32 for baryons, to 0.24 for strange
mesons, and to 0.20 for pions. The latter two values have been adjusted to reproduce
experimental Feynman x distributions of secondaries at projectile energies of
Elab = 250 GeV [51]. The spectator nucleons surviving from primary nucleus projectiles
may be treated by various fragmentation options as described in section 5.3.
30
Structure function integrals
The presently implemented version1
of VENUS uses the Duke-Owens parametrization
[52] of the parton distribution functions for sea and valence quarks with flavour
i ∈ {u, d, s, u, d, s} as
f
sea/val
i (z) =
Aizai
(1 − z)bi
(1 + αiz)
z2 + 4µ2/s
where Ai, ai, bi, and αi are parameters adjusted to experimental values; µ is a cutoff
parameter of the order of 1 GeV , and s is the squared energy in the cm system. The
original VENUS code uses a rather time consuming numerical integration to evaluate
the structure function integrals I(s, x)
I(s, x) =
x
0
fi(z) dz
at the beginning of a run to perform many collisions at a constant energy. In EAS
the energy varies from collision to collision, and hence we need these integrals several
thousand times for a complete shower. To accelerate the calculation we make use
of special properties of the integrals. For the sea quarks the exponents bi attain approximately
integer values and the integrals are (to good approximations) evaluated
analytically. In the valence quark case we keep the integrals as tables for a set of x
values in the range 0 ≤ x ≤ 1, calculated for s → ∞ with high precision. For
each individual case of s = E2
cm we add a correction term ∆I, which depends on the
collision energy like
∆I = − exp (c1 + c2 ln plab)
with fitted constants c1 and c2 and plab the projectile laboratory momentum. This
term is analytically evaluable, and the integrals are sufficiently well approximated
for x ≥ 5xcut = 20µ2
/s. For the few cases with x values below this limit the valence
quark structure function integrals must be evaluated numerically. Mainly by these
means, the computing time for a single collision at Elab = 1015
eV is reduced to
≤ 4 % for proton projectiles and to ≤ 25 % for iron nuclei projectiles.
5.1.2 QGSJET
QGSJET (Quark Gluon String model with JETs) is an extension of the QGS model
[53, 54], which describes hadronic interactions on the basis of exchanging supercritical
Pomerons. Pomerons are cut according the Abramovski˘i-Gribov-Kancheli rule and
form two strings each. These strings are fragmented by a procedure similar to the
Lund algorithm [55, 56], but with deviating treatment of the momenta at the string
1
VENUS version 4.12, released April 9, 1993
31
ends. Additionally QGSJET [8] includes minijets to describe the hard interactions
which are important at the highest energies.
In the case of nucleus-nucleus collisions the participating nucleons are determined
geometrically by Glauber calculations, assuming a Gaussian distribution of the nuclear
density for the light nuclei with A ≤ 10 and a Woods-Saxon distribution for the
heavier nuclei. The collision is treated by application of the percolation-evaporation
fragmentation mechanism [57] of the spectator parts of the involved nuclei. Thus
in peripheral collisions spallation like reactions happen, while in central collisions a
more or less complete fragmentation into many small fragments takes place. To these
fragments the various fragmentation options may be applied as described in section
5.3.
5.1.3 DPMJET
DPMJET (Dual Parton Model with JETs Version II.4) is based on the two component
Dual Parton Model [9, 58, 59] and contains multiple soft chains as well as multiple
minijets. As VENUS and QGSJET, it relies on the Gribov-Regge theory and the
interaction is described by multi-Pomeron exchange. Soft processes are described by
a supercritical Pomeron, while for hard processes additionally hard Pomerons are introduced.
High mass diffractive events are described by triple Pomerons and Pomeron
loop graphs, while low mass diffractive events are modelled outside the Gribov-Regge
formalism. Cutting a Pomeron gives two strings, which are fragmented by the JETSET
routines [60] according to the Lund algorithm [55, 56].
If nuclei are involved in a collision the number of participating nucleons as well as
the number of nucleon interactions is obtained by Glauber theory using the algorithm
of Ref. [61]. The refined treatment of the residual nuclei by the formation zone
intranuclear cascade model [62] takes into account the nuclear excitation energy,
models nuclear evaporation, high energy fission and break-up of light nuclei, and
emission of deexcitation photons for projectile and target nuclei. This might be
of some importance for projectiles in light-ion induced EAS [63]. For comparative
studies also the various fragmentation options are available (see section 5.3). The
effects of the various fragmentation options on the measurable quantities of EAS are
discussed in Ref. [64].
Short living secondaries not known within CORSIKA are to decay within DPMJET.
DPMJET produces also charmed hadrons which cannot be treated by CORSIKA.
Therefore within these charmed hadrons the charm quark is replaced by a
strange quark and the modified strange hadrons are tracked and undergo interactions
or decays within CORSIKA.
5.1.4 SIBYLL
SIBYLL (Version 1.6) [10, 11] is a minijet model essentially designed for usage in
EAS Monte Carlo programs. For hadronic soft collisions both projectile and target
32
particles fragment into a quark-diquark rsp. quark-antiquark system, that forms
a triplet and anti-triplet of colour. The components of opposite colour of the two
hadrons are then combined to form two colour strings that are fragmented according
to a slightly modified version of the Lund algorithm [55, 56]. Additional strings
originating from hard collisions with minijet production with high pT are considered.
The rise of the inelastic cross section with energy is fully attributed to an increasing
number of minijets, while the contribution of the soft component is assumed to be
energy independent. Diffractive events are modelled independently of soft or hard
collisions.
In hadron-nucleus collisions the number ν of interacting target nucleons determines
the number of soft strings: Each participating target nucleon is split into two
components, while the projectile is split into 2 ν components: two valence components
that carry the the quantum numbers of the incident projectile, and (ν − 1)
quark-antiquark component pairs. Again these partons are joined in pairs of strings
with opposite colour coupling projectile and target components, and these strings
are fragmented.
Nucleus-nucleus interactions are treated in a semi-superposition model, where the
number of interacting projectile nucleons is determined by Glauber theory, while the
projectile spectator nucleons fragment into light to medium heavy nuclei according
to a thermal model [11]. Moreover the various fragmentation options as described in
section 5.3 are available.
Short lived secondary particles produced by SIBYLL decay immediately into particles
which are known within CORSIKA (see section 2.2), but only nucleons and
antinucleons, charged pions, and all four species of kaons can be treated as projectiles
by SIBYLL. Other particles like strange baryons are tracked but only decay and
no interaction is admitted. In photonuclear interactions the incident gamma ray is
substituted by a charged π.
5.1.5 HDPM
As an alternative model of the interactions between hadrons and nuclei at high
energies, a phenomenological generator called HDPM may be used. For historical
reasons it is the default. It is developed by Capdevielle [4] and inspired by the Dual
Parton Model [5]. The HDPM routines are usable up to Elab ≈ 1017
eV . A detailed
description of the physical ideas and the parameters adjusted to results of pp-collider
experiments up to Elab ≈ 2 · 1015
eV is given in Appendix D.
5.2 Strong interactions at low energies
5.2.1 GHEISHA
The GHEISHA package is recommended for the treatment of low energy hadronic
interactions. We coupled it to CORSIKA in the same way as it is implemented into
33
the detector simulation code GEANT3 [24]. The GHEISHA routines [12] are designed
for laboratory energies up to some hundred GeV . We use them only at laboratory
energies below 80 GeV to treat the interactions of hadronic projectiles with the nuclei
of the atmosphere. All hadronic projectiles including baryons with strangeness ±1,
±2, and ±3 can be handled, however nuclear fragments emerging from evaporation
processes like d, t, and α particles cannot be treated by GHEISHA. The cross sections
for elastic and inelastic interactions are interpolated and extrapolated from tabulated
values derived from experimental data and stored within the GHEISHA package.
Neutron capture is taken into account as a third process for neutrons with Elab ≤
0.033 GeV . As the atmosphere does not contain fissile materials, we have eliminated
the routines treating nuclear fission. From the cross sections the type of interaction
is chosen at random, and the multiplicity and kinematic parameters of the secondary
particles are sampled with the appropriate GHEISHA reaction routine.
The GHEISHA routines treat low energy neutrons in a very consistent way. This
must be compared with the ISOBAR model, in which the low energy neutrons are
scattered around without energy loss in a rather unrealistic manner, thus resulting in
numerous low energy neutrons. Therefore we recommend the use of the GHEISHA
option despite the longer calculation times.
It should be noted, that in GHEISHA only the elements H, Al, Cu, and Pb
are tabulated as target materials and that the interesting cross section data for the
target elements N, O, and Ar which compose the atmosphere must be detained by
interpolation, with a loss of accuracy. In standard high energy physics experiments
air is not employed as target ior detector material, therefore a check on the validity
of this interpolation is impeded.
5.2.2 ISOBAR
The routines of the ISOBAR model of Grieder [3] work at cm energies between
0.3 GeV and 10 GeV . In this model the hadron-nucleus collisions are approximated
by hadron-nucleon reactions. The non-interacting nucleons of the target are
neglected. The hadron-nucleon reactions are assumed to take place via several intermediate
states which are decaying immediately into up to 5 secondaries. The
intermediate state can be a single particle, a heavy non-strange or strange meson, or
a light or a heavy ISOBAR. Intermediate states produced simultaneously share the
available cm energy according to their masses and move only forward or backward
with respect to the laboratory direction. These intermediate states cannot be identified
with single well established particles or resonances, but are to represent the
manifold of short lived states observed in this energy region which decay mostly into
few secondaries. Details on the mass parameters of these isobars, on their decay into
secondary particles, on the selection of transverse momenta, and on the annihilation
of antinucleons are given in Ref. [21]
34
5.3 Nuclear fragmentation
For the non-interacting nucleons of a projectile nucleus, the so called spectators, various
options exist [65]. First, they can be regarded as free nucleons with their initial
velocity. They are stored on the internal particle stack and are processed further at
a later time. This option assumes the ‘total fragmentation’ of the projectile nucleus
in the first interaction. A second option, ‘no fragmentation’, keeps all spectators together
as one nucleus propagating further through the atmosphere and reacting later
on. These two options, being the limiting cases of what really happens in nature,
allow to estimate the influence of the nuclear fragmentation on the results of the air
shower simulations. Our calculations show that the differences between the two cases
are small and details of fragmentation are smeared out by fluctuations.
Two further options keep also the spectator nucleons together, but calculate the
excitation energy of this remaining nucleus from the number of wounded nucleons in
a way as described in Ref. [66]. Each wounded nucleon is assumed to be removed out
of the Woods-Saxon potential of the original projectile nucleus, which leads to hole
states in the nucleonic energy level system. The depth of the nuclear potential well
amounts typically to 40 MeV . As hole states may occur in any arbitrary nuclear
level, each wounded nucleon contributes with an energy between 0 and 40 MeV to
the total excitation with equal probability. From this total excitation energy the
number of evaporating nucleons or α particles is calculated assuming a mean energy
loss per evaporated nucleon of 20 MeV [67]. The validity of this treatment is testified
[68] in a comparison with data from the CERN EMU07 experiment [69].
The emitted nucleons carry a transverse momentum, which differs in the two
available evaporation options following the parametrizations of Refs. [70, 71]. But
the global effect of the latter two options - a widening of the lateral distribution of
the whole shower compared with the total fragmentation or no-fragmentation options
- is usually small as compared to statistical fluctuations [64].
35
36
Chapter 6
Particle decays
Most of the particles produced in a high energy interaction are unstable and may
decay into other stable or unstable particles. Neutral pions and η mesons as well as
all resonance states have such a short life time that interaction is negligible before
they decay. Muons are prevented from penetrating the complete atmosphere by
decay only. Neutrons are treated as stable particles due to their long life time. For
all the other unstable particles there is a competition between interaction and decay
processes and the decision is taken when calculating the actual free path as described
in the beginning of chapter 4. If several decay modes exist, all known modes with a
branching ratio down to 1% are taken into account. In this section we describe the
treatment of particle decays in CORSIKA.
6.1 πo
decays
Neutral pions decay predominantly (98.8 %) into 2 photons πo
−→ γ +γ. This decay
is isotropic in the cm system of the πo
. There, the photon energy is Eγ cm = mπo /2
and the angle with respect to the direction of motion in the laboratory system is θcm.
In the laboratory system which moves with βπo with respect to the cm system the
energies and angles of the photons can be found by Lorentz transformation
Ei
γ lab =
1
2
γπo mπo (1 ± βπo cos θcm)
cos θi
lab =
βπo ± cos θcm
1 ± βπo cos θcm
i = 1, 2 .
The values of cos θcm and the angle φ around the incident direction are selected at
random to get a uniform distribution over the whole solid angle.
The Dalitz decay πo
−→ e+
+e−
+γ happens in 1.2 % of the cases only. It is a bit
more complicated than the decay into two photons, because the particle energies in
the cm system are not fixed, but vary in the kinematically allowed range. However,
momentum and energy conservation restrict the three secondary particles to lie in one
plane and to share the cm energy. The energies are calculated to fall inside a Dalitz
37
plot which gives a probability density function for a decay depending on the variables
p12, p13 where pik = pi + pk = −pl are the sums of the momenta of two particles
each. This probability density function directly depends on the matrix element |M|2
of the decay. As we do not know this matrix element for the πo
decay, we assume
a constant value. The energies of the secondary particles are taken at random from
this distribution. The cm values of cos θ and angle φ of the first secondary particle
are chosen at random uniformly in the full solid angle. The angle ψ of the reaction
plane around the first particle direction is also chosen at random. As final step, the
particle kinematical parameters are transformed to the laboratory system.
6.2 π±
decay
The decay π±
−→ µ±
+ νµ is a two-body decay, isotropic in the cm system of the
pion. Therefore, cos θcm and φcm of the muon are taken from a uniform distribution
and the energy is shared between muon and neutrino in a way that their cm momenta
add up to zero. This leads to
Eµ cm =
m2
π + m2
µ
2mπ
= mµγµ cm = 1.039mµ
and after Lorentz transformation into the laboratory system
γµ lab = γπ(γµ cm + βπ cos θµ cm γ2
µ cm − 1)
cos θµ lab =
γπγµ lab − γµ cm
γπβπ γ2
µ lab − 1
.
The muon carries a longitudinal polarization
ξ =
1
βµ
Eπ lab
Eµ lab
2r
1 − r
−
1 + r
1 − r
with r = (mµ/mπ)2
as given in [72]. We calculate the muon spin direction relative
to the laboratory frame and assume that this direction is maintained until the muon
decays. The calculation of the neutrino kinematical parameters is optional and gives
Eν lab = mπγπ − mµγµ
cos θν lab =
βπ − cos θµ cm
1 − βπ cos θµ cm
.
6.3 Muon decay
At the end of its track, a muon can only decay via µ±
−→ e±
+νe +νµ. The electron
energy distribution in the cm system is [73]
dNe
dEe cm
∝ 3
m2
µ + m2
e
2mµ
E2
e cm − 2E3
e cm
38
from which the electron cm energy Ee cm is taken at random. The direction correlation
of this three-body decay is governed by the longitudinal polarization of the
muon. The electron emission direction relative to the muon spin is determined with
the uniformly distributed angle cos δ and with ζ = 1 for µ±
to be
cos θe cm = ζ
1 + A(2 cos δ + A) − 1
A
with A =
1 − 2x
2x − 3
where x is the ratio of the electron energy to its maximum value
x =
2mµEe cm
m2
µ + m2
e
.
By angular addition to the muon polarization direction we get the electron emission
angle θ∗
e cm relative to the muon flight direction. The Lorentz transformation into
the laboratory system with the mean velocity βµ leads to the laboratory energy and
direction
Ee lab = meγe lab = γµ(Ee cm + βµpe cm cos θ∗
e cm)
cos θe lab =
γµ
me γ2
e lab − 1
(pe cm cos θ∗
e cm + βµEe cm) .
The optional calculation of the neutrinos parameters follows an ansatz of Jarlskog
[74] and uses the muon spin direction. In the cm system the angle ϕ between the
electron direction and the muonic neutrino direction defines the quantity ˆϕ by
ˆϕ =
1
2
(1 − cos ϕ) .
The angle ψ describes the rotation of the plane, in which this three-body decay
takes place, relative to the plane given by the muon spin and the electron emission
direction. The probability distribution for a decay with the electron energy x and its
emission angle θe cm is given [75] to
dP ∝
ˆϕ
(1 − ˆϕx)4
1 − ˆϕ(2x − x2
) − ζ cos θe cm 1 − ˆϕ(2 − 2x + x2
)
−ζ sin θe cm cos ψ2(x − 1)
√
ˆϕ − ˆϕ2 d ˆϕdψ .
From this distribution the two angles ϕ and ψ are taken at random and, hence, the
kinematical parameters of the neutrinos are calculated. Finally, they are transformed
to the laboratory system in analogy with the electron parameters.
6.4 Kaon decays
Kaon decays produce a variety of final states consisting mostly of two or three particles.
The dominant decays and their branching ratios are listed in Table 6.1.
39
Decay mode Branching Decay mode Branching
ratio (%) ratio (%)
K±
−→ µ±
+ ν 63.5 Ko
S −→ π+
+ π−
68.6
K±
−→ π±
+ πo
21.2 Ko
S −→ 2πo
31.4
K±
−→ π±
+ π±
+ π 5.6 Ko
L −→ π±
+ e + ν 38.7
K±
−→ πo
+ e±
+ ν 4.8 Ko
L −→ π±
+ µ + ν 27.1
K±
−→ πo
+ µ±
+ ν 3.2 Ko
L −→ 3πo
21.8
K±
−→ πo
+ πo
+ π±
1.7 Ko
L −→ π+
+ π−
+ πo
12.4
Table 6.1: Decay modes and branching ratios for kaons.
Two-body decays
The two-body decays are isotropic in the cm system and, hence, can be treated in
analogy to the charged pion decay (see section 6.2). The secondary particles are
emitted back-to-back and their Lorentz factors are
γ1 cm =
m2
K
+ m2
1 − m2
2
2mK
m1
and γ2 cm =
m2
K
− m2
1 + m2
2
2mK
m2
.
After transformation to the laboratory system, the γ factors and angles are
γi lab = γK
γi cm + βK
cos θi cm γ2
i cm − 1
cos θi lab =
γK
γi lab − γi cm
γK
βK
γ2
i lab − 1
i = 1, 2 .
In the leptonic two-body kaon decay, the muon polarization direction is calculated
analogously to the π±
decay (see section 6.2). For the optional parameter calculation
of the neutrino emerging from this decay we use
Eν lab = mK
γK
− mµγµ
cos θν lab =
βK
− cos θµ cm
1 − βK
cos θµ cm
.
Three-body decays
The situation for the three-body decays is treated in analogy with the neutral pion
Dalitz decay (see section 6.1). The probability density function represented by the
matrix element |M|2
is parametrized in Ref. [76] for the decay of kaons into 3 pions
by a series expansion of the form
|M|2
∝ 1 + g
s3 − s0
m2
π+
+ h
s3 − s0
m2
π+
2
+ j
s2 − s1
m2
π+
+ k
s2 − s1
m2
π+
2
+ · · ·
40
Decay mode g h k
K±
−→ π±
+ π±
+ π -0.22 0.01 -0.01
K±
−→ πo
+ πo
+ π±
0.59 0.035 0.0
Ko
L −→ π+
+ π−
+ πo
0.67 0.08 0.01
Ko
L −→ 3πo
0.0 -0.00033 0.0
Table 6.2: Coefficients of the parametrization of K −→ 3π.
Decay mode λ+ λ0
K±
−→ πo
+ e±
+ νe 0.028 0.0
K±
−→ πo
+ µ±
+ νµ 0.033 0.004
Ko
L −→ π±
+ e + νe 0.03 0.0
Ko
L −→ π±
+ µ + νµ 0.034 0.025
Table 6.3: Coefficients of the parametrization of K −→ π + + ν.
where
si = (pK
− pi)2
= (mK
− mi)2
− 2mK
Ti i = 1, 2, 3
s0 =
1
3 i
si =
1
3
(m2
K
+ m2
1 + m2
2 + m2
3) .
If CP invariance holds, j must be zero in good agreement with measurements. The
values of the other coefficients are given in Table 6.2.
The leptonic three-body kaon decays show a probability density function that can
be parametrized as described in Ref. [77]
|M|2
∝ G2
+ mK
(2E cmEν cm − mK
Eπ) + m2 1
4
Eπ − Eν cm +
Hm2
Eν cm −
1
2
Eπ + H2 1
4
m2
Eπ
with m being the mass of the lepton and
H =
m2
K
− m2
π
m2
π
(λ0 − λ+)G−
G± = 1 ± λ+
m2
K
+ m2
π − 2mK
Eπ cm
m2
π
Eπ =
m2
K
+ m2
π − m2
2mK
− Eπ cm .
The parameters λ+ and λ0 that fit the data best are given in Table 6.3.
41
In the muonic three-body kaon decays, the muon polarization direction is calculated
using the formulas given in Refs. [78, 79, 80]. In the leptonic three-body kaon
decays, the kinematical parameter calculation of the emerging neutrino is optional.
6.5 η decays
The η mesons decay to two photons or – by three-body decays – to pions, occasionally
accompanied by γ radiation [81]. The most frequent decay modes used in CORSIKA
are listed in Table 6.4. The decay into two photons proceeds analogously to the πo
decay as described in section 6.1, with exception of the larger decaying rest mass.
The three-body decays are handled in analogy with the πo
Dalitz decay (see section
6.1) using the Dalitz density algorithm. A constant matrix element |M|2
is assumed
for these decays, except for the η −→ π+
+π−
+πo
decay. For the latter the matrix
element |M|2
is expanded to [81, 82]
|M|2
∝ 1 + ay + by2
+ · · · with y =
3T0
mη − mπ+ − mπ− − mπo
− 1
where T0 is the kinetic energy of the neutral pion. The coefficients of the series
expansion are adopted [82] to a = −1.07 and b = 0.
Decay mode Branching
ratio (%)
η −→ γ + γ 39.13
η −→ 3 πo
32.09
η −→ π+
+ π−
+ πo
23.84
η −→ π+
+ π−
+ γ 4.94
Table 6.4: Decay modes and branching ratios for η.
6.6 Strange baryon decays
All strange baryons decays are two-body decays, which are isotropic in the cm system
and, hence, can be treated in analogy with the charged pion decay (see section 6.2).
The decay modes and branching ratios as used in CORSIKA are listed in Table 6.5.
The decay modes of the strange anti-baryons are not listed, as they correspond
to the baryons, just interchanging the particles by their anti-particles.
42
Decay mode Branching Decay mode Branching
ratio (%) ratio (%)
Λ −→ p + π−
64.2 Ξo
−→ Λ + πo
100.0
Λ −→ n + πo
35.8 Ξ−
−→ Λ + π−
100.0
Σ+
−→ p + πo
51.64
Σ+
−→ n + π+
48.36 Ω−
−→ Λ + K−
67.8
Σo
−→ Λ + γ 100.0 Ω−
−→ Ξo
+ π−
23.6
Σ−
−→ n + π−
100.0 Ω−
−→ Ξ−
+ πo
8.6
Table 6.5: Decay modes and branching ratios for strange baryons.
6.7 Resonance decays
All resonances decay without tracking into two daughter particles. These decays are
isotropic in the cm system and, hence, can be treated in analogy with the charged
pion decay (see section 6.2). The decay modes and branching ratios for the resonance
decays are derived from the combination of the quark content of the resonance
with the various qq-pairs from the sea. They are listed in Table 6.6. The decay
modes of the anti-baryon resonances are not listed, as they correspond to the baryon
resonances, just interchanging the particles by their anti-particles.
Decay mode Branching Decay mode Branching
ratio ratio
∆++
−→ p + π+
1 K∗o
−→ K+
+ π−
2/3
∆+
−→ p + πo
2/3 K∗o
−→ Ko
L + πo
1/6
∆+
−→ n + π+
1/3 K∗o
−→ Ko
S + πo
1/6
∆o
−→ n + πo
2/3 K∗±
−→ K±
+ πo
2/3
∆o
−→ p + π−
1/3 K∗±
−→ Ko
L + π±
1/6
∆−
−→ n + π−
1 K∗±
−→ Ko
S + π±
1/6
K
∗o
−→ K−
+ π+
2/3
ρo
−→ π+
+ π−
1 K
∗o
−→ Ko
L + πo
1/6
ρ±
−→ π±
+ πo
1 K
∗o
−→ Ko
S + πo
1/6
Table 6.6: Decay modes and branching ratios for resonances.
43
44
Chapter 7
Electromagnetic interactions
7.1 Muonic interactions
Muons may suffer from bremsstrahlung and e+
e−
-pair production. Both processes
are negligible below 2 TeV , but become important with increasing energy. The
programming of both types of interactions is taken from the detector simulation
code GEANT3 [24]. The production of subthreshold electromagnetic particles is
subsumed into the continuous energy loss by ionization, while only such interactions
are treated explicitly, which generate electromagnetic particles above the threshold.
7.1.1 Muonic bremsstrahlung
The simulation of discrete bremsstrahlung emitted from a muon of energy E is based
on Ref. [38]. The differential cross section for the emission of a bremsstrahlung
photon of energy k is given by
dσ
dv
=
N
v
4
3
−
4
3
v + v2
Φ(δ) (7.1)
with the energy fraction v = k/E. The constant N is a normalization factor and
δ =
m2
µ
2E
v
1 − v
.
Dependent on the charge number Z of the traversed medium the function Φ(δ) is
calculated from
Φ(δ) =



ln
189mµ
meZ1/3
− ln
189
√
e
meZ1/3
δ + 1 : Z ≤ 10
ln
189mµ
meZ1/3
− ln
189
√
e
meZ1/3
δ + 1 + ln
2
3
Z−1/3
: Z > 10
45
with e = 2.718 . . . Euler’s constant. The energy k and hence the fraction v is limited
by the photon threshold kc and by the maximum kinematically possible value vmax
by
vc =
kc
E
≤ v ≤ 1 −
3
4
√
e
mµ
E
Z1/3
= vmax .
Factorizing Equation 7.1 we write
dσ
dv
= f(v)g(v) (7.2)
with
f(v) = [v ln (vmax/vc)]−1
g(v) = 1 − v + 3
4
v2
Φ(δ)/Φ(0) .
The sampling of the photon energy is performed with two independent random numbers
RD1 and RD2 distributed uniformly between 0 and 1 by sampling v from
v = vc (vmax/vc)RD1
.
With this fraction v the function g(v) is calculated and compared with the second
random number RD2. If g(v) ≤ RD2 the fraction v is accepted and the photon
energy k calculated, otherwise a new set of random numbers is drawn. The energy
of the muon is reduced by the photon energy, but the muon direction is unchanged.
The angle between photon and muon momentum is sampled from a universal angular
distribution function, which approximates the real distribution function of Ref. [83].
The photon azimuthal angle is distributed isotropically around the muon direction.
7.1.2 Muonic e+
e−
-pair production
The double differential cross section for radiating off an e+
e−
pair in a medium with
charge number Z by a muon of energy E and mass mµ is given [38] by
d2
σ
dν dρ
= α4 2
3π
(Zλ)2 1 − ν
ν
Φe + (me/mµ)2
Φµ
with the fine structure constant α = 1/137, the electron Compton wave length
λ = 3.8616·10−11
cm, the fraction of energy transferred to the pair ν = (E+
+E−
)/E
with the energies E±
and mass me of the e±
particles. The complete formulas for Φe
and Φµ may be found in the Appendix of Ref. [38]. The quantity ρ gives the energy
asymmetry of the e+
e−
-pair
ρ =
E+
− E−
E+ + E−
.
The kinematic ranges of ν and ρ are defined by
4me/E = νmin ≤ ν ≤ νmax = 1 − 3
4
√
eZ1/3
mµ/E
0 = ρmin ≤ |ρ(ν)| ≤ 1 −
6m2
µ
E2(1−ν)
1 − 4me
νE
46
with e = 2.718 . . . Euler’s constant. Under the assumptions that the shapes of d2σ
dν dρ
and dσ
dν
dρ d2σ
dν dρ
do not depend on Z, the dominant contribution comes from the low
ν region
νmin = 4me/E ≤ ν ≤ 100νmin ,
and that d2σ
dν dρ
is independent of of ρ in this region, the differential cross section may
be approximated [24] by
dσ
dν
= dρ
Kd2
σ
dν dσ
≈
1
νa
1 −
νmin
ν
(7.3)
with a = 2 − 0.1 ln E, where E is given in GeV . Analogously to Equation 7.2 we
factorize Equation 7.3
dσ
dν
≈ f(ν)g(ν) where f(ν) =
(a − 1)
1
νa−1 − 1
νmax
a−1
1
νa
is the normalized distribution between Ec/E = νc ≤ ν ≤ νmax (with the threshold
energy Ec for e±
) and
g(ν) = 1 −
νmin
ν
is the rejection function. With a valid ν the maximum energy asymmetry is calculated
ρmax = 1 −
6m2
µ
E2(1 − ν)
1 −
4me
Eν
and the actual value of the asymmetry parameter ρ is chosen at random uniformly
in the range −ρmax ≤ ρ ≤ + ρmax.
For the polar angle θ of the e±
momentum relative to the muon momentum the
approximate average value θ = mµ/E is taken, while the azimuthal angle φ+
is taken
to be uniformly distributed and φ−
= φ+
+ π. The muon gets the final energy
Eµ = E − E+
− E−
, while its original direction is kept.
7.2 Electromagnetic subshowers
Electron and photon reactions are treated with EGS4 (Electron Gamma Shower
system version 4) or with the analytic NKG (Nishimura Kamata Greisen) formula.
The former delivers detailed information (momentum, space coordinates, propagation
time) of all electromagnetic particles, but needs extended computing times increasing
linearly with the primary energy, while the latter works fast, but gives only electron
densities at selected points in the detection plane.
47
7.2.1 Electron gamma shower program EGS4
The EGS4 option enables a full Monte Carlo simulation of the electromagnetic component
of showers by calling the EGS4 package which for electrons or positrons
treats annihilation, Bhabha scattering, bremsstrahlung, Møller scattering, and multiple
scattering (according to Moli`ere’s theory). Gamma rays may undergo Compton
scattering, e+
e−
-pair production, and photoelectric reaction. The programming of
these standard interactions is well documented in Ref. [6] and therefore not described
here. The direct µ+
µ−
-pair production and the photonuclear reaction with protons
and neutrons of nuclei of the atmosphere have been added. Despite their small cross
sections, these two processes are essential for the muon production in gamma ray
induced showers.
Photoproduction of muons and hadrons
The µ+
µ−
pair formation is treated in full analogy with the e+
e−
pair formation
replacing the electron rest mass by the muon rest mass. In the high energy limit the
cross section for this process approaches
σµ+µ− =
m2
e
m2
µ
σe+e−
and reaches 11.4 µb above Eγ = 1 TeV .
The photonuclear reaction cross section [84] of protons is shown in Fig. 7.1.
Its parametrization comprises three resonances at Eγ = 0.32, 0.72, and 1.03 GeV
superimposed on a continuum which slightly increases with energy
σγp = 73.7s0.073
+ 191.7/s0.602
1 − s0/s .
Here s and s0 are the squared cm energy rsp. pion production threshold energy in the
cm system in GeV 2
and σγp is given in µb. The measurements of H1 [85] and ZEUS
[86] confirm the extrapolation to higher energies. The photonuclear cross section
of air is calculated from the proton cross section by multiplication with the factor
A0.91
= 11.44 [84, 87]. The cross sections and branching ratios for all processes are
provided in a cross section file as usual in EGS4.
In photonuclear reactions the target nucleons are treated as free particles with
the assumption that only one nucleon is involved in the photonuclear process. To
generate secondary particles various possibilities exist which are selected depending
on the energy of the gamma ray. Below 0.4 GeV only one pion is generated, while in
the subsequent range up to 1.4 GeV the chance to generate one pion decreases linearly
in favour of the generation of two pions. The choice between these two possibilities
is made at random. Between 1.4 GeV and 2 GeV always two pions are produced.
Within the range of 2 GeV < Eγ < 3 GeV the selection between two pion generation
and the HDPM option is made at random with a linearly decreasing chance for two
pion generation. Above 3 GeV multi-particle production by the HDPM is always
48
0
0.1
0.2
0.3
0.4
0.5
0.6
10
-1
1 10 10
2
10
3
10
4
10
5
10
6
10
7
Eγ (GeV)
σγp(mb)
H1
ZEUS
Figure 7.1: Photoproduction cross section of protons as function of energy. Data
points are from Refs. [87, 85, 86].
assumed and at energies above 80 GeV the selected high energy hadronic interaction
model is employed.
In the production of one single pion, the recoil nucleon may undergo charge
exchange with 33% probability according to the decay modes of the ∆+
and ∆o
resonances (see Table 6.6). The pion azimuthal angle is chosen at random from a
uniform distribution, while the polar angle is selected by a rejection method which
produces a dipole or a quadrupole radiation characteristics depending on the energy
of the gamma ray. If a charged pion is produced, these characteristics are modified
to meet approximately the experimentally determined angular distributions [88].
When two pions together with a recoil nucleon are produced the particle energies
are chosen to fall inside a Dalitz plot with a constant probability density. The
treatment is analogous to the πo
Dalitz decay into three secondary particles (see
section 6.1). Charge exchange of the recoil nucleon is allowed giving altogether 6 exit
channels respecting charge conservation. The branching ratios are listed in Table 7.1.
The production of more than two secondary particles is treated by the HDPM as
described in Appendix D. In this case a neutral pion is assumed to be the interacting
particle and diffraction is suppressed, but charge exchange to π±
is admitted. If the
49
Reaction mode Branching Reaction mode Branching
ratio (%) ratio (%)
γ + p −→ πo
+ πo
+ p 15 γ + n −→ πo
+ πo
+ n 15
γ + p −→ π+
+ π−
+ p 15 γ + n −→ π+
+ π−
+ n 15
γ + p −→ π+
+ πo
+ n 20 γ + n −→ π+
+ πo
+ p 20
Table 7.1: Branching ratios for photonuclear reactions leading to two pions.
photon energy exceeds 80 GeV the selected high energy interaction model is used to
describe the production of secondary hadrons, and the incoming photon is replaced
by a πo
or an η meson with equal chance (rsp. by a π+
in SIBYLL).
Modifications of the standard EGS4
The essential modifications of the standard EGS4 code [6] are summarized as follows.
The barometric density dependence of air as described in section 2.4 is implemented
into the particle tracking used within the EGS4 routines. Important is the
path length correction of the mean free path to the next interaction according to a
medium with increasing density.
The deflection of electrons and positrons in the Earth’s magnetic field (see section
3.3) is calculated by an approximation only valid for small deflection angles [6].
As low energy particles at high altitude may have considerable path lengths and,
hence, large deflection angles, the step size is limited to keep the deflection angle
below 11.5o
for each step. The propagation time (including fast renormalization of
direction cosines) is calculated for the total curved path length of the particles also
in the case of magnetic field deflection.
The pressure dependence of the Sternheimer correction [89] for ionization losses
in air is modified. The standard EGS4 cross section files contain the continuous
energy loss dE/dx of electrons and positrons for energies above 1MeV by ionization
in gaseous matter depending on the pressure at a fixed density. The deposited energy
per radiation length X0 in air rises linearly with the logarithm of the energy (in MeV )
as
dE/dx = (61.14 + 5.58 ln E)MeV/X0
until it saturates at an energy that depends on pressure. Expressed as a function of
height h in cm, the saturation energy loss is
(dE/dx)sat = (86.65 + 810−6
h)MeV/X0 .
This approximates the pressure dependence of the energy loss by ionization to better
than 5%.
The reduction of computing time becomes important with increasing primary energy,
as this time increases linearly with the energy dissipated in the atmosphere.
50
Therefore, the probability of electrons or gamma rays to produce a charged particle
at the next observation level is estimated [90] as a function of their altitude and energy.
If the probability remains below a preselected lower limit depending on shower
size, this particle is discarded unless it is closer than 3 radiation lengths to the next
observation level, the gamma ray energy exceeds the pion production threshold of
0.152 GeV , or the electron energy twice exceeds this threshold value. The latter
conditions assure, that the production of pions which may decay into muons with
large penetration depth is not suppressed. This discarding mechanism eliminates the
numerous calculations of low-energy subshowers which do not contribute essentially
to observable particles, thus reducing the computation time by a factor of 3. The
number of lost electrons typically amounts to 3 to 5% for showers initiated by a
1015
eV primary and explains the difference in the electron numbers emerging from
simulations with EGS4 and NKG options. However, when calculating the longitudinal
shower profiles, using the Cherenkov option, or using the thinning option this
discarding mechanism is switched off to keep all charged particles, as they contribute
to the profile and to the emitted Cherenkov radiation or – in the case of thinning –
substantially enlarge the number of arriving electrons by their higher weight.
As a second measure for reduction of computing time, the maximum step size
between two multiple scattering events of electrons and positrons has been increased
by a factor of 10 relative to the value of Ref. [6]. Tests [90] showed, that the influence
on the lateral electron distribution is negligible for showers as would be measured by
the KASCADE experiment. To keep the CORSIKA program flexible to problems,
in which a more frequent treatment of the electron multiple scattering is essential –
this is the case for the lateral distribution of the Cherenkov photon density initiated
by gamma ray primaries with energies as low as 20 to 1000 GeV [91] – the step size
increase factor may be set individually to lower values. A detailed discussion on the
step length is also given in Refs. [6, 92].
To enable simulations of EAS with the highest observed energies (≈ 1020
eV ),
the cross sections and branching ratios are extented also to 1020
eV , assuming that
QED is valid. The influence of the Landau-Pomeranchuk-Migdal effect [93] is small
[94] for proton induced showers even at primary energies of 5 · 1020
eV . This would
not hold for γ-induced showers of these energies. We consider the effect in a manner
as it is done in the programs MOCCA and AIRES [95].
7.2.2 Nishimura-Kamata-Greisen option
In the NKG option the electromagnetic component of air showers is calculated by
an analytical approach [96] without a full Monte Carlo simulation. The advantage
of relatively modest computer time requirements for the analytical treatment is paid
for with less accurate information about the electromagnetic particles. Coordinates
with arrival time, location, and momenta of single electromagnetic particles cannot
be obtained, but only total electron numbers at various atmospheric depths together
with some parameters that give information about the general development of the
51
electromagnetic component of a shower. At one or two observation levels lateral
electron densities are computed for a grid of points around the shower axis (see
below).
Longitudinal electromagnetic shower development
The longitudinal development of the electromagnetic part of showers is obtained
by calculating the total number of electrons for ≤ 10 values of atmospheric depth
separated by 100 g/cm2
down to the lowest observation level. For each subshower
initiated by gamma rays (from πo
or η decays) or by electrons of energy E, the age
si of this subshower at each interesting level i in depth Ti (in g/cm2
) is calculated
si =
3Ti/X0
Ti/X0 + 2 ln(E/Ecrit)
(7.4)
with Ecrit = 82 MeV being the critical energy and X0 = 37.1 g/cm2
the radiation
length in air 1
. The electron number of a subshower at level i amounts to [96]
Ne i(Ee > 0) =
0.31
ln(E/Ecrit)
e(1−1.5 ln si)Ti/X0
. (7.5)
To improve Equation 7.5 to account for an energy threshold Ethr for the electrons,
the Ne i are multiplied with a correction factor
Ne i(Ee > Ethr) = KNe i(Ee > 0) .
The correction factor K has the form [16]
K = Ne i(Ee > 0, Ek2)/Ne i(Ee > 0, Ek1)
which is the quotient of electron numbers which are gained with different critical
energies Ek1 = 0.4Ecrit and Ek2 = Ek1 + Ethr. This correction factor amounts to 4,
10, and 30 % for threshold energies of 1, 3, and 10 MeV . At each interesting depth
value, these electron numbers Ne i are summed up for all subshowers.
Longitudinal age parameter
In order to describe the shower development of the overall electromagnetic or hadronic
cascade in the same way as the parameter s does for individual electromagnetic
subshowers, an artificial parameter, the global longitudinal age slong, is introduced.
1
The Particle Data Group [97] suggests a radiation length of X0 = 36.7 g/cm2
which would
be consistent with the value used within EGS. Also a critical energy of Ecrit = 86 MeV may be
derived form this reference. These values will be adopted in the next release.
52
Using the age parameter s determined according to Equation 7.4, the parameters
a(s), b(s), C1(s), and C2(s) of the NK structure functions are calculated [98, 14]
a(s) =
4
s
e0.915(s−1)
b(s) = 0.15 +
1
1 + s
C1(s) =
as/b
2π
Γ
s
b
+
4Γ(s+1
b
)
sa1/b
−1
C2(s) =
a(s+1)/b
2π
Γ
s + 1
b
+
4Γ(s+2
b
)
sa1/b
−1
.
With the ratio of the coefficients C1(s)/C2(s), the parameter R at each depth is
derived by summing over all subshowers j
R = j (NjC1(s)/C2(s))
j Nj
.
with Nj the electron numbers according the improved Equation 7.5. Finally, slong is
defined as
slong =
B2 − 4A(C − R) − B
2A
.
In this relation the coefficients A, B, and C depend on R and are given in Table 7.2.
Range of R A B C
0.0191 – 0.1796 0.3109 0.2146 -0.0055
0.1796 – 0.5364 0.3667 0.1639 0.0060
0.5364 – 1.0332 0.1460 0.6317 -0.2420
1.0332 – 1.4856 -0.3376 2.0903 -1.3438
Table 7.2: Parameters of the longitudinal age formula.
Lateral electron distribution
The lateral distribution of electromagnetic showers in different materials scales well
with the Moli`ere radius rmol = 0.0212 GeV · ξ0/Ecrit. Here ξ0 is the radiation length
in cm. In the atmosphere ξ0 varies with the density, hence, rmol = 9.6 g cm−2
/ρair.
About 90% of the energy of a shower is deposited inside a cylinder around the shower
axis with radius rmol. In CORSIKA the electron distribution is determined for the two
53
lowest observation levels. The electron density ρe at distance r from the subshower
axis is calculated according to Ref. [99]
ρe =
Ne
2πr2
mols2
m
Γ(4.5 − s)
Γ(s)Γ(4.5 − 2s)
r
rmolsm
s−2
1 +
r
rmolsm
s−4.5
(7.6)
with the modulation function of Lagutin et al. [14, 15]
sm = 0.78 − 0.21s .
The electron densities are calculated for 80 reference points centered around the
shower axis on a circular grid extending in 8 directions spaced by 45o
and with 10
radial distances in each direction, covering the range from 1 m to a maximum radius
which has to be specified (see Ref. [22]) in logarithmic steps. We are aware of the fact
that Equation 7.6 is correct only in the Landau approximation [100] which restricts
(besides others) the distance r between the considered grid point and the subshower
axis to
r
rmol
>
Ecrit
E
.
This condition is violated when the subshower axis comes close to one of the grid
points. The densities of all subshowers are summed up at each reference point of the
grid to get the local densities of the total shower.
7.3 Cherenkov radiation
Basing on a program extension written by HEGRA collaborators [101] we have implemented
an option to simulate Cherenkov radiation. This radiation is emitted, if
the velocity v of charged particles – mainly electrons, but also muons and charged
hadrons – exceeds the local speed of light, which is given by the local refractive index
n and the vacuum speed of light c. Therefore we examine each transportation step
of a charged particle for the condition
nv/c = nβ > 1 .
The refractive index n can be approximated [102] by the local density ρ(h) (in g/cm3
)
n = 1 + 0.000283 ρ(h)/ρ(0)
neglecting the wavelength dependence of n. The number NC
of photons which are
emitted per path length s at an angle θC
is calculated to
dNC
ds
= 2πα
sin2
θC
λ2
dλ
54
where the integral extends over the wavelength band, within which the Cherenkov
detector is sensitive; α is the fine structure constant, and the angle θC
relative to the
charged particle direction is given by
θC
= arccos
1
βn
.
To adapt the Cherenkov simulation easily to various detector types, a wavelength
band may be specified by data input.
Each charged particle path element is subdivided into smaller fractions such that,
within each fractional element, the number of emitted photons is less than a predefined
bunchsize. This photon bunch is treated as a whole rather than each single
Cherenkov photon, thus reducing the computational effort considerably. For each
photon bunch the azimuthal emission angle is taken at random, and its arrival coordinates
at the detector plane are calculated.
Atmospheric absorption of the Cherenkov photons is not taken into account. But
by writing the origin height of each photon bunch onto the Cherenkov output, the
absorption may be introduced later when analyzing the output data.
55
56
Chapter 8
Outlook
The CORSIKA program in its present state models EAS initiated by photons, protons
and nuclei up to the highest primary energies. It employs a number of theoretical
models of high energy hadronic interactions which are adjusted to reproduce experimental
data wherever possible. CORSIKA is a useful and flexible tool to study high
energy cosmic ray interactions, to support the interpretation of EAS measurements,
and to optimize the design of future cosmic ray experiments.
However, the CORSIKA program is under continuous evolution and many details
of the shower development are subject to uncertainties and approximations.
Wherever we are aware of such an uncertainty, we try to improve it. Some of the
improvements to be implemented in the near future have already been mentioned in
the text.
Unfortunately for EAS, the collider results have to be extrapolated into energy
and angular regions where the interactions are supposed to change. Gluons instead
of quarks become the most abundant reaction partners, heavy quarks and minijets
are produced, and the collider events might look different in the forward region
compared to the region with high transverse momentum where the collider detectors
are located. Here theoretically founded approaches realized in different computer
codes are employed. By their coupling with CORSIKA we hope to get rid of one
single model, and the variations between the various models give an idea on the error
of the Monte Carlo predictions of measurable quantities. We expect advances in the
theories to describe high energy hadronic interactions, stimulated by the advent of
new experimental data.
57
58
Appendix A
Atmospheric parameters
Layer i Altitude h (km) ai (g/cm2
) bi (g/cm2
) ci (cm)
1 0 . . . 4 −186.5562 1222.6562 994186.38
2 4 . . . 10 −94.919 1144.9069 878153.55
3 10 . . . 40 0.61289 1305.5948 636143.04
4 40 . . . 100 0.0 540.1778 772170.16
5 > 100 0.01128292 1 109
Table A.1: Parameters of the U.S. standard atmosphere.
Layer i Altitude h (km) ai (g/cm2
) bi (g/cm2
) ci (cm)
1 0 . . . 4 −118.1277 1173.9861 919546.
2 4 . . . 10 −154.258 1205.7625 963267.92
3 10 . . . 40 0.4191499 1386.7807 614315.
4 40 . . . 100 5.4094056 · 10−4
555.8935 739059.6
5 > 100 0.01128292 1 109
Table A.2: Parameters of the AT115 atmosphere (January 15, 1993).
59
Layer i Altitude h (km) ai (g/cm2
) bi (g/cm2
) ci (cm)
1 0 . . . 4 −195.837264 1240.48 933697.
2 4 . . . 10 − 50.4128778 1117.85 765229.
3 10 . . . 40 0.345594007 1210.9 636790.
4 40 . . . 100 5.46207 · 10−4
608.2128 733793.8
5 > 100 0.01128292 1 109
Table A.3: Parameters of the AT223 atmosphere (February 23, 1993).
Layer i Altitude h (km) ai (g/cm2
) bi (g/cm2
) ci (cm)
1 0 . . . 4 −253.95047 1285.2782 1088310.
2 4 . . . 10 −128.97714 1173.1616 935485.
3 10 . . . 40 0.353207 1320.4561 635137.
4 40 . . . 100 5.526876 · 10−4
680.6803 727312.6
5 > 100 0.01128292 1 109
Table A.4: Parameters of the AT511 atmosphere (May 11, 1993).
Layer i Altitude h (km) ai (g/cm2
) bi (g/cm2
) ci (cm)
1 0 . . . 4 −208.12899 1251.474 1032310.
2 4 . . . 10 −120.26179 1173.321 925528.
3 10 . . . 40 0.31167036 1307.826 645330.
4 40 . . . 100 5.591489 · 10−4
763.1139 720851.4
5 > 100 0.01128292 1 109
Table A.5: Parameters of the AT616 atmosphere (June 16, 1993).
60
Layer i Altitude h (km) ai (g/cm2
) bi (g/cm2
) ci (cm)
1 0 . . . 4 − 77.875723 1103.3362 932077.
2 4 . . . 10 −214.96818 1226.5761 1109960.
3 10 . . . 40 0.3721868 1382.6933 630217.
4 40 . . . 100 5.5309816 · 10−4
685.6073 726901.3
5 > 100 0.01128292 1 109
Table A.6: Parameters of the AT822 atmosphere (August 22, 1993).
Layer i Altitude h (km) ai (g/cm2
) bi (g/cm2
) ci (cm)
1 0 . . . 4 −242.56651 1262.7013 1059360.
2 4 . . . 10 −103.21398 1139.0249 888814.
3 10 . . . 40 0.3349752 1270.2886 639902.
4 40 . . . 100 5.527485 · 10−4
681.4061 727251.8
5 > 100 0.01128292 1 109
Table A.7: Parameters of the AT1014 atmosphere (October 14, 1993).
Layer i Altitude h (km) ai (g/cm2
) bi (g/cm2
) ci (cm)
1 0 . . . 4 −195.34842 1210.4 970276.
2 4 . . . 10 − 71.997323 1103.8629 820946.
3 10 . . . 40 0.3378142 1215.3545 639074.
4 40 . . . 100 5.48224 · 10−4
629.7611 731776.5
5 > 100 0.01128292 1 109
Table A.8: Parameters of the AT1224 atmosphere (December 24, 1993).
61
Altitude AT115 AT223 AT511 AT616 AT822 AT1014 AT1224 U.S.
(m) stand.
0 1055.9 1044.6 1031.3 1043.3 1025.5 1020.1 1015.1 1036.1
1000 916.8 900.9 900.7 909.9 895.6 888.9 879.2 901.3
2000 810.4 789.9 799.8 807.0 796.7 787.4 774.3 797.6
3000 715.0 690.2 707.7 713.7 707.9 695.0 679.7 703.7
4000 629.3 600.6 623.7 628.9 628.1 611.0 594.4 618.9
6000 483.0 451.0 479.3 483.8 489.8 467.5 450.6 473.9
8000 364.1 335.9 362.7 366.8 374.2 352.9 337.9 358.4
Table A.9: Air pressure values (in hPa) at low altitude for various atmospheres.
62
Appendix B
Muon range for horizontal showers
The disadvantage of the precise formulation of the muon range in section 4.1 is its
dependence form the special atmospheric model with a density increasing exponentially
along the muon path. This is no longer the case for atmospheric profiles along
nearly horizontal showers with θ > 75o
.
Starting from Equation 4.2 we consider the slowing down of the muon along its
path replacing γµβµ in Equation 4.3 rsp. 4.2 by a suitable function. In a first
approximation the function γµ( )βµ( ) is expanded into a power series, omitting
higher order terms:
γµ( )βµ( ) = γ0 µβ0 µ + ∂(γµβµ)
∂
+ · · ·
= γ0 µβ0 µ + ∂γµ
β0 µ∂
+ · · · .
(B.1)
Here γ0 µ is the Lorentz factor and β0 µ the velocity at the starting point. For
the ionization loss of a muon which traverses air of thickness λ we take the simple
expression
dEµ = −
λ
β2
κ = −
λγ2
γ2 − 1
κ
with the energy loss by ionization of κ = 2 MeV g−1
cm2
. Thus we get
∂γµ
∂
=
∂Eµ
mµ∂
= −
κρ(h0)
mµβ2
0 µ
where ρ(h0) gives the atmospheric density at the altitude of the muon starting point,
and mµ is the muon rest mass. Solving Equation 4.3 for the approximation of Equation
B.1 leads to the individual range of the muon of
=
− ln(RD)cτµγ0 µβ0 µ
1 − ln(RD)cτµκρ(h0)/(mµβ3
0 µ)
(B.2)
with RD a random number distributed uniformely between 0 and 1. In contrast to
section 4.1 this approximation is independent of a special atmospheric model, as the
63
-200
0
200
400
600
800
1000
1200
1400
0 2 4 6 8 10 12 14 16
altitude (km)
difference(m)
0.20.25
0.2
0.25
0.15
0.2
0.25
0.3
0.350.4
0.45
0.5
0.55
0.6
0.35
0.4
0.45
0.5
0.55
0.15
0.2
0.25
0.3
0.35
0.4
γ = 50
γ = 100
x10
Figure B.1: Differences in the muon range calculated for vertical muons of γ = 50
(open symbols) and 100 (filled symbols) starting at 16 km altitude. The difference
to the real penetration depth is shown for different approximations: Equation 4.2,
Equation B.2, and • (with x10 enlarged scale) Equation 4.5. The numbers give the
fraction of the mean life time in the muon rest system.
density at the starting point may be obtained e.g. by interpolation from numerical
tables. Therefore this approximation is advantageously applied in the horizontal
shower version (see section 2.4), despite the limited precision.
This less precise treatment leads to an enlarged number of muons at sea level (ca.
10% for Eµ < 15 GeV ) compared with the approximation of Equation 4.8. Figure
B.1 shows the deviations of the different approximations from the real penetration
depth of vertical muons.
64
Appendix C
Default cross sections
C.1 Nucleon-nucleon cross sections
This cross section is experimentally available for nucleon laboratory momenta plab up
to 1000 GeV/c [87] which corresponds to a center of mass (cm) energy of 44.7 GeV .
The measured nucleon-nucleon cross section [87] can well be parametrized as
σn−n(plab) = A + BpN
lab + C ln2
plab + D ln plab (C.1)
where A, B, C, D, and N are free parameters of the fit. Their values are given in
Table C.1. The momentum plab is given in GeV/c and σn−n is in mb.
For larger momenta the cross section is extrapolated by
σn−n(plab) = 22.01(p2
lab + m2
)0.0321
. (C.2)
This represents an empirical fit1
to the proton-antiproton inelastic cross section which
is known up to 1.8 TeV cm energy and which is expected to be equal to the nucleonnucleon
cross section at these energies.
Employing the ISOBAR model, the measured inelastic cross section drops rapidly
at low energies. Below plab = 10 GeV/c the ISOBAR model allows only elastic
reactions with a constant cross section of
σn−n(plab ≤ 10 GeV/c) = 29.9 mb .
For antinucleons an annihilation with nucleons can occur in addition, leading to a
contribution to the inelastic cross section which is parametrized in Ref. [87] by
σan(plab) = 0.532 + 63.4p−0.71
lab . (C.3)
1
The most recent issue of the Review of Particle Physics [97] gives an exponent of 0.0395.
With a reduced normalization constant of 19.87 we get an identical cross section at 1000 GeV .
Corresponding constants in Equation C.4 for π − n and K − n reactions are 13.25 and 11.01. These
values will be adopted in the next release. This more recent approximation shifts up the HDPM
curves at 109
GeV in Fig. 4.1 by ≈ 8%, in Fig. 4.3 by ≈ 18%, and in Fig. 4.5 by 8 − 12%.
65
Param. n − n π − n K − n
A 30.9 24.3 12.3
B -28.9 -12.3 -7.7
C 0.192 0.324 0.0326
D -0.835 -2.44 0.738
N -2.46 -1.91 -2.12
Table C.1: Parameters of the hadron-nucleon cross section parametrization.
From these nucleon-nucleon cross sections the nucleon-air and nucleus-air cross
sections are calculated in the same manner as described in section 4.2.3 using Glauber
theory.
C.2 Pion-nucleon and kaon-nucleon cross sections
Existing measurements of π−n and K−n reactions show a similar dependence on the
laboratory momentum plab as nucleons. Therefore, Equation C.1 can be fitted to the
measured π and K data as well. The results of such fits have been taken from Ref. [87]
and are listed in Table C.1. In the momentum region above 1000 GeV/c the cross
sections for π and K are assumed to rise with the same momentum dependence as
for nucleons. In order to get a continuous transition between the two energy regimes,
only the scaling factors were modified compared to the nucleon case in Equation C.2
σπ−n(plab) = 14.70(p2
lab + m2
)0.0321
σK−n(plab) = 12.17(p2
lab + m2
)0.0321
.
(C.4)
Again in the ISOBAR model the cross sections are taken to be constant below plab =
5 GeV/c for pions and plab = 10 GeV/c for kaons to account for elastic scattering in
this energy region. The used cross section values are
σπ−n(plab ≤ 5 GeV/c) = 20.64 mb
σK−n(plab ≤ 10 GeV/c) = 14.11 mb .
From these meson-nucleon cross sections the meson-air cross sections are calculated
in the same manner as described in section 4.2.3 using Glauber theory.
66
Appendix D
HDPM
The HDPM generator describes the interactions between hadrons and nuclei at high
energies and is worked out by Capdevielle [4] inspired by the Dual Parton Model
(DPM) [5]. It is a phenomenological description of the interaction based on the
picture that two dominant colour strings are formed between the interacting quarks of
two hadrons. For instance, in a nucleon-nucleon collision, two chains (colour strings)
are stretched between the fast valence di-quark of the projectile and one valence
quark of the target and vice versa between the slow valence quark of the projectile
and the di-quark of the target [5, 58]. The strings separate and fragment into many
colour neutral secondaries that are produced around the primary quark directions.
Such particle jets have been observed in many high energy physics experiments.
Recent experiments at pp colliders have improved the understanding of such reactions
up to cm energies of 1.8 TeV . Unfortunately, the collider data contain mainly
particles that are produced under large angles with respect to the direction of incidence
(central rapidity region). But the major part of the energy escapes with
the particles in the beam pipe. For the development of EAS, however, the particles
emitted in forward direction are the most important ones, because they carry the
energy down through the atmosphere.
In the central region, many quantities, such as the number and type of secondaries,
the longitudinal and transverse momentum distributions and the spatial energy flow
have been measured and correlated with each other and with the available energy.
The rich data collection of collider experiments has been used to build an interaction
model that reproduces the collider results as well as possible.
A difficulty arises, because air shower simulations need a description of nucleonnucleus
or even nucleus-nucleus collisions rather than nucleon-nucleon interactions.
In the following sections, it is described how these interactions are modelled on the
basis of the present knowledge about nucleon-nucleon reactions.
67
D.1 Nucleon-nucleon interactions
Number of secondaries
The average charged-particle multiplicity (including the colliding particles) in nucleon-nucleon
collisions has been measured up to cm energies of 1.8 TeV . It can be
parametrized [103] as a function of s = E2
cm (inGeV 2
) by
nch =



0.57 + 0.584 ln s + 0.127 ln2
s :
√
s ≤ 187.5 GeV
6.89s0.131
− 6.55 : 187.5 GeV <
√
s ≤ 945.5 GeV
3.4s0.17
: 945.5 GeV <
√
s .
(D.1)
The actual charged-particle multiplicity nch for each event fluctuates around the
average value nch . The fluctuations follow a negative binomial distribution
P(nch, nch , k) =
nch + k − 1
nch
nch /k
1 + nch /k
nch
1
1 + nch /k
k
where P gives the probability to obtain nch particles for the parameters nch and k.
The dependence of k on the cm energy
√
s (in GeV ) is parametrized [104] by
1/k = −0.104 + 0.058 ln
√
s .
From this distribution the actual number of charged particles nch is picked at random.
At lower energies the particle numbers are dominated by pions, therefore the
average number of neutral particles nneu produced should be around nch /2, and
the average total multiplicity around N = 1.5 nch at low energies. At high
energies, however, a larger fraction of photons than expected by nγ = 2nneu = nch
has been observed, mainly due to η meson production with subsequent decay into
pions and photons.
The average number of gamma quanta nγ adopted to reproduce this excess
follows the parametrization [105]
nγ =
−1.27 + 0.52 ln s + 0.148 ln2
s :
√
s ≤ 103 GeV
−18.7 + 11.55s0.1195
:
√
s > 103 GeV .
(D.2)
where nγ is the average over many collisions, i.e. over all inclusive data. In contrast
to nγ , we distinguish nγ as the average number of photons in collisions with the same
number of charged secondaries nch. The value of nγ is deduced from the correlation
nγ = 2 + anch with a =
0.0456 log s + 0.464 :
√
s < 957 GeV
1.09 :
√
s ≥ 957 GeV .
We prefer this energy dependence of a compared to the constant value a = 1.03 given
by the UA5 Collaboration [106], because it better describes the variations observed
68
in the energy range 200 GeV ≤
√
s ≤ 900 GeV . For low energies, in any case, nneu
is forced to be at least equal to nch /2.
The actual number nγ for a given collision is derived [107] from the observed
relation between nch and nγ in collider data [108]. For energies
√
s > 200 GeV ,
the probability distribution of nγ around nγ = nγ z is described by a truncated
Gaussian distribution, whose mean m and variance σ depend on z = nch/ nch which
is used as a convenient scaling variable
m = nγ(0.982 − 0.376e−z
)0.88
σ = m (0.147 + 2.532e−z
)0.88 .
For
√
s < 200 GeV the fluctuations are taken to be the same as for the charged
particles nγ = z nγ .
The parent particles of the photons are assumed to be mainly neutral pions, but
also ρ and η mesons, kaons, and hyperons can contribute to the γ component.
Particle ratios
The abundances of kaons, nucleons, Λ and Σ particles in nucleon-nucleon collisions
were measured by the UA5 collaboration [109]. We adopt the UA5 parametrization
[106] of the ratio of charged kaon to charged pion numbers
nK± /nπ± = 0.024 + 0.0062 ln s (D.3)
and the ratio of the number of nucleons to the number of all charged particles
nN /nch = −0.008 + 0.00865 ln s . (D.4)
Λ, Σo
, and Σ±
particles are produced with the same probability and their ratio to
the number of all charged particles is
nΛ
nch
=
nΣo
nch
=
nΣ+ + nΣ−
nch
=
1
3
(−0.007 + 0.0028 ln s) .
A noticeable contribution to photon production originates from η mesons [110]. Their
abundance relative to the πo
mesons is assumed [111] to be slightly energy dependent
to
nη/nπo = 0.06 + 0.006 ln s + 0.0011 ln2
s . (D.5)
Taking into account all these particle ratios and the specific decay modes of
the particular particles, all the particle numbers of nucleons, pions, kaons, etas,
and hyperons are determined to meet the previously selected charged and neutral
multiplicity nch and nγ for each single collision. It should be noted that kaons,
nucleons, and hyperons are always produced as particle-antiparticle pairs.
69
Rapidity distribution
In hadronic reactions, jets of secondary particles are generated by the hadronizing
colour strings. The kinematics of particles of a jet are described by their transverse
momenta pT
and rapidities y where the latter are defined by
y =
1
2
ln
E + pL
E − pL
with E being the particle energy and pL
the longitudinal momentum. In the representation
of rapidity and transverse momentum, the rapidities of the particles of a jet
are approximated by a Gaussian distribution as suggested by Klar and H¨ufner [112]
and by inelastic hadron-hadron and lepton-hadron scattering data [113, 114]. The
two principal jets of a collision are back-to-back in the center of mass and therefore
positioned symmetrically around ycm = 0 in rapidity space. The average positions
of the centers of the respective Gaussians on the rapidity axis in the cm system
ymean and the average width σy are parametrized based on experimental data for
nucleon-nucleon collisions [4]
ymean = ±(0.146 ln(s − 2m2
N ) + 0.072)
σy = 0.12 ln(s − 2m2
N ) + 0.18 .
(D.6)
The rapidity of the cm system in the laboratory frame is expressed by
ycm =
1
2
ln
Elab + mN + plab
Elab + mN − plab
with mN being the nucleon mass and plab the nucleon momentum in the laboratory
system. The amplitude of the rapidity distributions is determined by deducing the
central rapidity density from data. Experimentally, however, only the pseudorapidity
η = − ln tan
θ
2
with θ being the cm production angle, is directly observed. The average pseudorapidity
density in the central region dN
dη (η=0)
as a function of the cm energy is obtained
from non single diffractive data [103] as
dN
dη (η=0)
=
0.82s0.107
:
√
s ≤ 680 GeV
0.64s0.126
:
√
s > 680 GeV .
(D.7)
The conversion from the measured mean pseudorapidity density in the central
region to the needed mean central rapidity density is performed by
dN
dy (y=0)
= fη→y
dN
dη (η=0)
70
where fη→y is kept constant at 1.25 for
√
s < 19.4 GeV . Above, its energy dependence
is deduced [115] from calculations with the DPM [58] to
fη→y = 1.28852 − 0.0065 ln(s − 2m2
N ) .
The central rapidity density dN
dy (y=0)
eventually is calculated from the average central
pseudorapidity density dN
dη (η=0)
in dependence from the scaling variable z as [4, 116]
dN
dy (y=0)
=
dN
dη (η=0)
fη→y
(0.487z + 0.557)2
: z ≥ 1.5
(0.702z + 0.244)2
: z < 1.5 .
(D.8)
The amplitude Ay of the Gaussians is deduced from the requirement that all
particles belong to the Gaussians
+∞
−∞
dN
dy
dy = nch leading to Ay = nch/(σy
√
8π) .
The position of the Gaussians on the rapidity axis is now calculated such that the
central rapidity density seen in semi-inclusive data is obtained by adding the two
rapidity distributions in the center
dN
dy
(y = 0) = 2Aye−y2
mean/2σ2
y .
Thus, ymean is computed by
ymean = σy 2 ln 2Ay/
dN
dy (y=0)
.
The quantity σy is taken to be σy from Equation D.6. The advantage of this
procedure is to find immediately the natural position of the set of rapidities.
In case of parent particles of photons, the same procedure as for charged particles
is applied to fix the position of the Gaussians as required to reproduce the theoretical
central rapidity densities. Optionally a slightly modified procedure may be adopted
to achieve a better agreement with experimental results [105, 117] which suggest a
narrower rapidity distribution for photons. The central rapidity densities dN
dy
γ
y=0
are
determined from Equations D.8 by replacing the scaling variable z by zγ = nγ/ nγ
and multiplying it with 0.5 to account for the average ratio of neutral to charged
pions and with an energy dependent factor g(s) given by
g(s) =



1 :
√
s ≤ 50 GeV
1 + 0.18 ln(
√
s/50 GeV ) : 50 GeV <
√
s ≤ 200 GeV
1.25 : 200 GeV <
√
s .
71
As in case of charged particles, the calculations of Attallah et al. [118] suggest a
more refined treatment of the conversion factor from pseudorapidity to rapidity also
for the neutrals
fγ
η→y =
dN
dy
γ
(y=0)
dN
dη
γ
(η=0)
=



1.1 :
√
s ≤ 19.4 GeV
1.33 − 0.0391 ln(s − 2m2
N ) : 19.4 GeV <
√
s ≤ 900 GeV
0.8 : 900 GeV <
√
s
to take into account UA5 results [105].
For each particle its rapidity yi in the cm system is chosen from the appropriate
Gaussian distribution at random.
Transverse momentum of the secondaries
The transverse momentum distribution of secondaries in nucleon-nucleon collisions
is well described [119] by
d2
N
dpxdpy
∝
p0
p0 + pT
n
. (D.9)
With pT
= p2
x + p2
y one obtains the probability density function
dN
dpT
=
(n − 1)(n − 2)
p2
0
p0
p0 + pT
n
pT
(D.10)
where p0 = 1.3 GeV/c for pions and the parameter n depends on the central pseudorapidity
density dN
dη (|η|<D)
= ρ|η|<D of charged secondaries as [120]
n = 7.4 + 3.67/ρ0.435
|η|<D .
The pseudorapidity density ρ|η|<D is calculated from the central region of width 2D
with
D = 0.67(2.95 + 0.0302 ln s) .
For collisions with
√
s > 500 GeV a recent analysis [121] of the UA1 (minimum
bias) experiment indicates a new correlation of pT
with the central pseudorapidity
density ρ|η|<D. For fixed p0 = 1.3 GeV/c this dependence is parametrized [122] to
pT
=
0.0033(ρ|η|<D − 1.56)2
+ 0.406 : ρ|η|<D < 3
0.0109ρ|η|<D + 0.383 : ρ|η|<D ≥ 3 .
To sample pT
from Equation D.9 for Ecm > 500 GeV the inverse integral method is
applied which leads to a transcendent equation of the form
RD
n − 1
pT
p0
+ 1
n−1
+
1
n − 1
−
pT
p0
= 0 . (D.11)
72
RD is a random number uniformly distributed between 0 and 1, the parameter n is
connected with pT
by
n =
2.64
pT
+ 3 .
Within CORSIKA Equation D.11 is solved by a rejection method.
According to these distributions, the transverse momenta are determined for secondary
pions only. As there is experimental evidence [123] on differences between the
transverse momentum distributions of secondary pions, kaons, nucleons, η mesons,
and strange baryons, this is accounted for by applying energy dependent correction
factors pK
T
/ pπ
T
, pN
T
/ pπ
T
, pη
T
/ pπ
T
, and psb
T
/ pπ
T
with
pπ
T
(s) =
0.3 + 0.00627 ln s : Ecm < 132 GeV
(0.442 + 0.0163 ln s)2
: Ecm ≥ 132 GeV
pK
T
(s) =
0.381 + 0.00797 ln s : Ecm < 131 GeV
(0.403 + 0.0281 ln s)2
: Ecm ≥ 131 GeV
pN
T
(s) =
0.417 + 0.00872 ln s : Ecm < 102 GeV
(0.390 + 0.0341 ln s)2
: Ecm ≥ 102 GeV
pη
T
(s) = 0.88 pK
T
+ 0.12 pN
T
psb
T
(s) = 1.45 pN
T
− 0.45 pK
T
.
A slight inconsistency should be noted. Instead of using the energy dependence of
pπ
T
from Equation D.1 the ln s parametrization is used in the correction factors
only. The advantage of this method is to follow closely the correlation of pT
with
the average central rapidity density.
The sum SpT
= N
i pT,i
of the transverse momenta of all secondaries is calculated
and the pT
values of the particles are reduced by SpT
/N to fulfill transverse
momentum conservation.
Energy of secondaries and leading particles
The laboratory energy of particle i is calculated according to
Ei = p2
T,i
+ m2
i cosh(yi + ycm)
for all but two particles. These extra particles are regarded as the remainder of
the collision partners, and their energies are treated in a differing manner. For the
‘anti-leader’ (the remnant of the participating target nucleon) the energy is taken
at random from a Feynman x distribution, which is parametrized [124] within three
regions as
dN
dxF
=



cxF
: 0 ≤ xF
< x1
cx1 : x1 ≤ xF
< x2
cx1 e−α(xF
−x2 )
: x2 ≤ xF
< 1
73
√
s < 13.8 GeV 13.8 GeV ≤
√
s < 5580 GeV 5580 GeV ≤
√
s
x1 0.2 0.71 + 0.00543 ln(s − 2m2
N ) 0.265
x2 0.65 0.8175 − 0.032 ln(s − 2m2
N ) 0.265
α 1.265 1.14 + 0.022 ln(s − 2m2
N ) 1.14 + 0.022 ln(s − 2m2
N )
Table D.1: Energy dependence of x1 , x2 , and α.
where c is a constant to normalize the distribution. The boundaries x1 and x2 and
the parameter α depend on energy as given in Table D.1.
The ‘leader’ (residue of the projectile) gets the remainder of energy after subtraction
of the energies of anti-leader and all secondaries. In case this is not possible,
because there is not enough energy left, the particle generation is repeated with a
new set of rapidities.
In high energy collisions, about 50% of the cm energy is carried away by secondary
particles, which is usually noted as an inelasticity parameter k ≈ 0.5 . This value is
reproduced rather well by the HDPM generator without any additional constraint.
By attributing the remaining energy to the ‘leader’ particle, this will in general be
the most energetic one.
Two further alternatives have been proposed in literature to determine the leading
particle rapidity. Alner et al. [106] attribute the largest of the randomly selected
rapidities of the secondaries to the leading particle. The second alternative picks the
leading particle rapidity from a separate distribution resulting from DPM calculations
[125] for the valence quarks recombined in the final state. There is no decisive
argument yet in favor of one of the three treatments.
To balance energy and longitudinal momenta of leader, anti-leader, and of all
secondaries simultaneously, all rapidities are slightly modified using the algorithm of
Jadach [126], which delivers the corrected rapidity yc i of particle i to
yc i = A + Byi .
The size of A is predominantly determined by longitudinal momentum conservation
and B by energy conservation; both parameters are approximated by an iterative
adjustment procedure for each collision. The transverse momentum of the leading
particles is chosen analogously to the secondary particles, depending on the particle
type.
Charge exchange and resonance formation
The leading particle after the collision is correlated with the primary incoming particle
due to the fact that the fast spectator quarks of the interaction move almost with
their initial velocity and most likely form the fastest secondary particle together with
an additional quark. This picture limits the possible types of the leading particles
74
leader other pion leader other pion
π−
+ πo
−→ πo
+ π−
π+
+ πo
−→ πo
+ π+
πo
+ π−
−→ π−
+ πo
πo
+ π+
−→ π+
+ πo
K−
+ πo
−→ Ko
L/S + π−
K+
+ πo
−→ Ko
L/S + π+
Ko
L/S + π+
−→ K+
+ πo
Ko
L/S + π−
−→ K−
+ πo
p + πo
−→ n + π+
n + π+
−→ p + πo
p + πo
−→ n + π−
n + π−
−→ p + πo
Table D.2: Charge exchange of leading particles. To conserve charge, another pion
has also to change its charge state.
that can appear for a given primary. The leading particle may be of the same type
as the incoming one, undergo charge exchange, or may be excited to a resonance
state. For reasons of charge conservation the charge of a further (secondary) pion is
changed in both processes. The possible charge exchange and resonance formation
combinations which are taken into account are listed in Table D.2 and Table D.3.
If more than one possibility exists for a leading particle to perform charge exchange
or resonance formation, a random selection is made respecting phase space considerations.
After charge exchange or resonance production the number of positive,
negative, and neutral pions is adjusted in a manner, that the total charge involved in
the collision is conserved, and that after decay of resonances the number of charged
and neutral pions approaches the experimentally observed values as close as possible.
The probability for charge exchange varies with energy as
Pex =



0.10 :
√
s ≤ 19.4 GeV
0.10 + 0.0345 ln(Elab/200 GeV ) : 19.4 GeV <
√
s ≤ 105 GeV
0.45 − 0.0537 ln(Elab/200 GeV ) : 105 GeV <
√
s ≤ 969 GeV
0.03 : 969 GeV <
√
s
(D.12)
and similarly the probability for resonance formation is taken to
Prf =



0.35 :
√
s ≤ 105 GeV
0.08819 ln(Elab/200 GeV ) : 105 GeV <
√
s ≤ 969 GeV
0.69 : 969 GeV <
√
s .
(D.13)
75
leader secondary resonance branching
pion ratio
π−
+ πo
−→ ρ−
1/2
π−
+ π+
−→ ρo
1/2
π+
+ π−
−→ ρo
1/2
π+
+ πo
−→ ρ+
1/2
K−
+ πo
−→ K∗−
1/2
K−
+ π+
−→ K
∗o
1/2
K+
+ πo
−→ K∗+
1/2
K+
+ π−
−→ K∗o
1/2
Ko
L/S + π−
−→ K∗−
1/4
Ko
L/S + π+
−→ K∗+
1/4
Ko
L/S + πo
−→ K∗o
1/4
Ko
L/S + πo
−→ K
∗o
1/4
p + π+
−→ ∆++
1/2
p + πo
−→ ∆+
1/3
p + π−
−→ ∆o
1/6
n + π+
−→ ∆+
1/6
n + πo
−→ ∆o
1/3
n + π−
−→ ∆−
1/2
p + π−
−→ ∆
−−
1/2
p + πo
−→ ∆
−
1/3
p + π+
−→ ∆
o
1/6
n + π−
−→ ∆
−
1/6
n + πo
−→ ∆
o
1/3
n + π+
−→ ∆
+
1/2
Table D.3: Resonance formation of leading particles. To preserve the final number
of pions, another pion must be subsumed into the resonance.
Both the ‘leader’ and the ‘anti-leader’ may independently undergo charge exchange or
resonance formation with the same probabilities. The charge exchange and resonance
formation reactions may be suppressed by a control flag.
The resonance formation of the leading particles involves important consequences
for the hadron cascading in EAS [127, 128] due to subsequent decay or modified
penetration depth in case of decay to the electromagnetic channel.
76
Diffractive processes
Most of the particles do not experience completely central collisions. In peripheral
collisions it may happen that the projectile is just excited by a rather small energy
and momentum transfer from a target nucleon. The excited projectile subsequently
decays and forms secondary particles. Such interactions are called diffractive processes.
Their topology is different from the non-diffractive events, mainly due to the
reduced amount of energy ESD
that is available for production of secondary particles.
As suggested by experimental data [129], a fraction of 15% of all interactions
is assumed to be diffractive in the present version of the program. Strictly speaking,
the ratio of diffractive to total cross section is slightly energy dependent following
the parametrization by [58]
σSD
= (1.77 ln0.7
s − 2.38) mb .
In principle, diffractive interactions are treated in the same way as non-diffractive
processes, with the following differences:
First, experimental results [130] and theoretical predictions [131] indicate the
excitation energy ESD
to follow
dσSD
d(E2
SD
/s)
∝
1
E2
SD
/s
.
ESD
must be large enough to produce at least one additional pion, but is limited to
at maximum 5% of the cm energy.
Second, the position and width of the Gaussians in rapidity space are calculated
as indicated in Equation D.6, however, replacing s by sSD
= E2
SD
and shifting the
zero point by
y0 SD
= ycm ± ln(Ecm/ESD
)
with the positive sign for projectile diffraction and negative sign for target diffraction.
With the same substitution for s the excess of photons from decaying secondaries is
described as explained in Equation D.2, and the particle ratios of kaons and nucleons
to pions are calculated following Equations D.3 and D.4.
Third, hyperon production in diffractive interactions is neglected, while η mesons
are produced according to Equation D.5.
The average central pseudorapidity density is taken in analogy to Equation D.7
to
ρη=0 = 0.82E0.214
SD
neglecting the parametrization for energies above 680 GeV . The average number of
charged particles varies with the energy in the same way as for the non-diffractive
case [111, 132]. Therefore, we adopt the parametrization of Equation D.1 and replace
s by sSD
.
With these parametrizations the same procedure is followed to generate secondaries,
their rapidities and transverse momenta. The energy of the diffracted leading
77
particle is taken from the Gaussian shaped rapidity distributions at random as for
the secondaries, while the rapidity of the non-diffracting leading particle is calculated
from the primary energy reduced by the mass of the diffracting system. Charge
exchange and resonance formation is only considered for the diffracting collision partner,
replacing s by sSD
in Equations D.12 and D.13. Due to the smaller amount of
energy available the overall number of secondaries is smaller. Energy and momentum
conservation is accomplished also by the Jadach filtering method, as described above.
D.2 Nucleon-nucleus interactions
In an EAS the incoming particle does not collide with free nucleons but with nuclei
of the air target. Consequently, the interaction model has to be adapted to this
situation. Since the nucleon-nucleon interactions are experimentally well studied, we
try to construct the nucleon-nucleus interaction basically in terms of the nucleonnucleon
interaction.
The treatment of nucleon-nucleus or nucleus-nucleus interactions starts with the
identification of the type of target nucleus. Therefore, the relative contributions of
the various air nuclei to the total inelastic cross section have been calculated and the
choice is made at random according to these contributions.
When a high energy nucleon hits an air nucleus, it does not interact with the
whole nucleus, but with a few target nucleons only. The number nT
of (‘wounded’)
target nucleons hit by the projectile can be determined in two different manners.
Either, a parametrization of nT
depending on the target mass number Atarget and
the square of the cm energy s is used [4]
nT
= (0.56 + 0.0236 ln(s − 1.76)) A0.31
target (D.14)
neglecting the fluctuations of nT
around its mean value. In the second option nT
is
explicitly selected according to its probability distribution [133] which is obtained by
Glauber calculations.
In case of diffractive interactions, we either set nT
= 1 or calculate it according
to Equation D.14 with s being replaced by sSD
.
Thus, the primary particle is assumed to interact with nT
nucleons of the target
successively. Obviously the multiple interactions in one nucleus are not independent
of each other. Our approach accounts for the multiple interactions in the target by
several corrections based on an analysis by Klar and H¨ufner [112]. The main one is
the production of additional secondaries, the target excess. This excess was measured
by observing the extra negative particles ∆n− from nucleon-nucleus collisions and
was parametrized [112] as
∆n− =
0.285(nT
− 1) nch : Ecm ≤ 137 GeV
0.25(nT
− 1) nch : Ecm > 137 GeV .
78
The neutral excess is then ∆nneu = ∆n− and the number of additional charged
particles from the target excess is ∆nch = 2 ∆n− . The additional particles
originate from a third string which is modelled by a third Gaussian distribution in
rapidity space.
When choosing the energies of the particles from the target excess, the rapidities
are taken at random from the third Gaussian. The particle types of the target excess
are determined following the same ratios as are used for the other secondaries.
The parametrization of position ym and width σ of this Gaussian, depending on
the number of reacting target nucleons nT
, is given [4] by
ym = −3 + 2.575e−0.082nT
σ = 1.23 + 0.079 ln nT
.
The final position of the third Gaussian is chosen in full analogy with the procedure
described above such that the particle excess in the center of rapidity equals the
observed values ρ3
y=0 with [134]
ρ3
y=0
= ρy=0
nT
− 1
2
.
The target excess lies at negative rapidity values in all cases.
According to HELIOS results [135] we assume the ratio of diffractive to total
inelastic cross section to be the same for nucleon-nucleon and nucleon-nucleus colli-
sions.
D.3 Pion-nucleus and kaon-nucleus interactions
The interactions of pions and kaons with a nucleus are simulated in strict analogy with
the nucleon-nucleus interaction. Only for the calculation of the available cm energy
and the determination of the number of target nucleons involved in the interaction
process, the different masses and cross sections of pions and kaons are taken into
account. All other features described in the nucleon case are the same for pions and
kaons.
D.4 Nucleus-nucleus interactions
The probability of nP projectile nucleons interacting and the probability of a projectile
nucleon hitting nT target nucleons are evaluated by Glauber calculations [45, 133]
and kept as tables (see section 4.2.3). The actual probabilities are interpolated from
these tables in analogy with the cross section values and used for a random selection
of nP and nT . The further reaction is now regarded as a superposition of nP
nucleon-nucleus reactions which are simulated as described in section D.2. The ratio
79
of interacting protons and neutrons is assumed to be equal to the ratio in the parent
nucleus.
The treatment of the non-interacting nucleons of the projectile (spectators) may
be selected to ‘total fragmentation’, ‘no fragmentation’, or evaporation. For the
evaporation treatment three options are available which differ only in the selection
of the transverse momenta (see section 5.3).
80
Acknowledgements
We would like to thank all colleagues in the various laboratories around the world
who have contributed to the development of the CORSIKA program by reporting
their problems, experience, suggestions, and detected errors.
Special thanks go to the authors of the various hadronic interaction programs R.S.
Fletcher, T.K. Gaisser, N.N. Kalmykov, P. Lipari, S.S. Ostapchenko, J. Ranft, and
K. Werner for the permission to use their programs and for many helpful discussions
and hints in coupling their codes with CORSIKA.
The contribution of the Cherenkov routines by the HEGRA collaborators F. Arqueros,
J. Cortina, S. Martinez (Madrid), and M. Rozanska (Cracow) is acknowledged
as well as the design of some system routines by G. Trinchero (Torino), by D.
Horn and M. Raabe (Hamburg). We are indebted to C. Pryke (Chicago) for many
communications, which helped us to tune CORSIKA to the highest energies.
We especially thank our Karlsruhe colleagues H.J. Gils, J. Oehlschl¨ager, and H.
Rebel, who regularly discussed with us all the problems arising in the simulations
with CORSIKA.
81
Bibliography
[1] P. Doll et al., The Karlsruhe Cosmic Ray Project KASCADE, Report KfK
4686 (1990), Kernforschungszentrum Karlsruhe; Nucl. Phys. B (Proc. Suppl.)
14A (1990) 336
[2] H.O. Klages for the KASCADE Collaboration, Nucl. Phys. B (Proc. Suppl.)
52B (1997) 92; Proc. 25th
Int. Cosmic Ray Conf., Durban, 6 (1997) 141
[3] P.K.F. Grieder, Report INS–J125 (1970), Inst. for Nuclear Studies, Univ. of
Tokyo; P.K.F. Grieder, Proc. 16th
Int. Cosmic Ray Conf., Kyoto, 9 (1979) 161
[4] J.N. Capdevielle, J. Phys. G: Nucl. Part. Phys. 15 (1989) 909
[5] A. Capella and J. Tran Thanh Van, Phys. Lett. B93 (1980) 146
[6] W.R. Nelson, H. Hirayama and D.W.O. Rogers, Report SLAC 265 (1985),
Stanford Linear Accelerator Center
[7] K. Werner, Phys. Rep. 232 (1993) 87
[8] N.N. Kalmykov, S.S. Ostapchenko, Yad. Fiz. 56 (1993) 105; Phys. At. Nucl.
56(3) (1993) 346; N.N. Kalmykov, S.S. Ostapchenko and A.I. Pavlov, Bull.
Russ. Acad. Sci. (Physics) 58 (1994) 1966
[9] J. Ranft, Phys. Rev. D51 (1995) 64
[10] R.S. Fletcher, T.K. Gaisser, P. Lipari and T. Stanev, Phys. Rev. D50 (1994)
5710
[11] J. Engel, T.K. Gaisser, P. Lipari and T. Stanev, Phys. Rev. D46 (1992) 5013
[12] H. Fesefeldt, Report PITHA–85/02 (1985), RWTH Aachen
[13] S.S. Ostapchenko, T. Thouw and K. Werner, Nucl. Phys. B (Proc. Suppl.) 52B
(1997) 3; K. Werner, private communication (1997)
[14] J.N. Capdevielle and J. Gawin, J. Phys. G: Nucl. Part. Phys. 8 (1982) 1317
[15] A.A. Lagutin et al., Proc. 16th
Int. Cosmic Ray Conf., Kyoto, 7 (1979) 18
82
[16] J.N. Capdevielle, KASCADE Collaboration, Proc. 22nd
Int. Cosmic Ray Conf.,
Dublin, 4 (1991) 405
[17] Proc. 24th
Int. Cosmic Ray Conf., Rome, (1995)
[18] Proc. 25th
Int. Cosmic Ray Conf., Durban, (1997)
[19] J. Knapp et al., Proc. 24th
Int. Cosmic Ray Conf., Rome, 1 (1995) 403
[20] D. Heck for the KASCADE Collaboration, Proc. 25th
Int. Cosmic Ray Conf.,
Durban, 6 (1997) 245
[21] J.N. Capdevielle et al., Report KfK 4998 (1992), Kernforschungszentrum
Karlsruhe
[22] J. Knapp and D. Heck, Report KfK 5196B (1993), Kernforschungszentrum
Karlsruhe; for an up-to-date version see http://www-ik3.fzk.de/˜heck/corsika/
[23] O. C. Allkofer and P. K. F. Grieder, Cosmic Rays on Earth, in: Physics Data
25/1, H. Behrens and G. Ebel eds. (Fachinformationszentrum Karlsruhe, Germany,
1983) chpt. 1.1.2
[24] R. Brun et al., GEANT3, Report CERN DD/EE/84–1 (1987), CERN,
Geneva
[25] M. Aguilar-Benitez et al. (Particle Data Group), Phys. Rev. D45 (1992) II, 6;
VIII, 40
[26] P.E. Hauenstein, Atomic Nucl. Data Tables 39 (1988) 290 (The values of A.H.
Wapstra, G. Audi, R. Hoekstra are taken.)
[27] Handbook of Chemistry and Physics, 67th
Edition, R.C. Weast ed. (The Chemical
Rubber Co., Cleveland, 1986) F141
[28] J. Linsley, private communication by M. Hillas (1988)
[29] H. Ulrich, Diploma thesis University Karlsruhe (1997) (unpublished)
[30] G. Marsaglia and A. Zaman, Report FSU–SCRI–87–50 (1987), Florida State
University
[31] F. James, Report CERN DD/88/22 (1988), CERN, Geneva
[32] R.M. Sternheimer, M.J. Berger and S.M. Seltzer, Atomic Nucl. Data Tables 30
(1984) 261
[33] G.Z. Moli`ere, Z. Naturforsch. 2a (1947) 133; Z. Naturforsch. 3a (1948) 78
[34] H.A. Bethe, Phys. Rev. 89 (1953) 1256
83
[35] B. Rossi, High Energy Particles (Prentice Hall, Englewood Cliffs, New Jersey,
1952)
[36] http://www.ngdc.noaa.gov/seg/potfld/geomag.html
[37] M. Hillas, Nucl. Phys. B (Proc. Suppl.) 52B (1997) 29
[38] W. Lohmann, R. Kopp and R. Voss, Report CERN 85-03 (1985), CERN,
Geneva
[39] G. Schatz and D. Heck, internal report,unpublished (1995)
[40] H.H. Mielke et al., J. Phys. G: Nucl. Part. Phys. 20 (1994) 637
[41] G.B. Yodh et al., Phys. Rev. D27 (1983) 1183
[42] T.K. Gaisser et al., Phys. Rev. D36 (1987) 13503
[43] R.M. Baltrusaitis et al., Phys. Rev. Lett. 52 (1993) 1380
[44] M. Honda et al., Phys. Rev. Lett. 70 (1993) 525
[45] R.J. Glauber and G. Matthiae, Nucl. Phys. B21 (1970) 135
[46] K.G. Boreskov and A.B. Kaidalov, Sov. J. Nucl. Phys. 48 (1988) 367
[47] H. de Vries et al., Atomic Nucl. Data Tables 36 (1987) 495
[48] J. Knapp, D. Heck and G. Schatz, Report FZKA 5828 (1996), Forschungszentrum
Karlsruhe
[49] J. Knapp, Report FZKA 5970 (1997), Forschungszentrum Karlsruhe
[50] J. Knapp et al., to be published
[51] M. Adamus et al., EHS–NA22 Collaboration, Z. Phys. C39 (1988) 311
[52] D.W. Duke and J.F. Owens, Phys. Rev. D30 (1984) 49
[53] A.B. Kaidalov and K.A. Ter-Martirosyan, Sov. J. Nucl. Phys. 39 (1984) 979
[54] A.B. Kaidalov, K.A. Ter-Martirosyan and Yu.M. Shabelsky, Sov. J. Nucl. Phys.
43 (1986) 822
[55] B. Andersson et al., Phys. Rep. 97 (1983) 31
[56] T. Sj¨ostrand, Comp. Phys. Comm. 39 (1986) 347
[57] N.N. Kalmykov and S.S. Ostapchenko, Sov. J. Nucl. Phys. 50 (1989) 315
84
[58] P. Aurenche, F.W. Bopp, A. Capella, J. Kwiecinski, M. Maire, J. Ranft and J.
Tran Thanh Van, Phys. Rev. D45 (1992) 92
[59] F.W. Bopp, D. Petermann, R. Engel and J. Ranft, Phys. Rev. D49 (1994) 3236
[60] T. Sj¨ostrand, Comp. Phys. Comm. 82 (1994) 74
[61] S.Y. Shmakov, V.V. Uzhinski and A.M. Zadoroshny, Comp. Phys. Comm. 54
(1989) 125
[62] A. Ferrari, J. Ranft, S. Roesler and P.R. Sala, Z. Phys. C70 (1996) 413; Z.
Phys. C71 (1996) 75
[63] G. Battistoni et al., Nucl. Phys. B (Proc. Suppl.) 52B (1997) 120
[64] D. Heck, internal report, unpublished (1997)
[65] G. Schatz and D. Heck, in: Report KfK 5027 (1992) 35, J. Knapp and H.
Rebel eds., Kernforschungszentrum Karlsruhe
[66] J.J. Gaimard, Th`ese Universit´e Paris 7 (1990); Report GSI–90–27 (1990),
Gesellschaft f¨ur Schwerionenforschung, Darmstadt
[67] X. Campi and J. H¨ufner, Phys. Rev. C24 (1981) 2199
[68] J.N. Capdevielle, Proc. 23rd
Int. Cosmic Ray Conf., Calgary, 4 (1993) 52
[69] A. D¸abrowska et al., Phys. Rev. D47 (1993) 1751
[70] A.S. Goldhaber, Phys. Lett. 53B (1974) 306
[71] T.H. Burnett et al., JACEE Collaboration, Phys. Rev. D35 (1987) 824
[72] G. Barr et al., Phys. Rev. D39 (1989) 3532
[73] T. Kinoshita and A. Sirlin, Phys. Rev. 107 (1957) 593
[74] C. Jarlskog, Nucl. Phys. 75 (1966) 659
[75] D. Heck, internal report, unpublished (1993)
[76] M. Aguilar-Benitez et al. (Particle Data Group), Phys. Rev. D45 (1992) VII, 77
[77] M. Aguilar-Benitez et al. (Particle Data Group), Phys. Rev. D45 (1992) VII, 84
[78] L. Jauneau, in: Methods in Subnuclear Physics, Vol. III, M. Nikoli˘c ed. (Gordon
& Breach, New York, 1969) 123
[79] L.M. Chounet et al., Phys. Rep. 4 (1972) 199, App. 1
85
[80] N. Cabbibo and A. Maksymowicz, Phys. Lett. 9 (1964) 352; Phys. Lett. 11
(1964) 360; Phys. Lett. 14 (1965) 72
[81] M. Aguilar-Benitez et al. (Particle Data Group), Phys. Rev. D45 (1992) VII, 6
[82] J.G. Layter et al., Phys. Rev. Lett. 29 (1972) 316
[83] Y.S. Tsai, Rev. Mod. Phys. 46 (1974) 815 and erratum Rev. Mod. Phys. 49
(1977) 421
[84] T. Stanev et al., Phys. Rev. D32 (1985) 1244
[85] T. Ahmed et al., (H1 Collaboration), Phys. Lett. B299 (1993) 374
[86] M. Derrick et al., (ZEUS Collaboration), Phys. Lett. B293 (1992) 465
[87] A. Baldini et al., in: Landolt-B¨ornstein, New Series I/12 a+b (Springer,
Berlin, 1987)
[88] H. Genzel et al., in: Landolt-B¨ornstein, New Series I/8 (Springer, Berlin, 1973)
[89] R.M. Sternheimer et al., Phys. Rev. B26 (1982) 6067
[90] J. Spitzer, private communication (1988)
[91] R. Attallah, private communication (1995)
[92] R.L. Ford and W.R. Nelson, Report SLAC 210 (1978), Stanford Linear Accelerator
Center
[93] L.D. Landau and I.Ya. Pomeranchuk, Dokl. Akad. Nauk SSSR 92 (1953) 535
& 735; A.B. Migdal, Phys. Rev. 103 (1956) 1811
[94] K. Kasahara, Proc. Int. Symp. Extremely High Energy Cosmic Rays, Tanashi,
Tokyo (Japan) (1996) 221
[95] S. Sciutto, private communication (1997)
[96] K. Greisen, in: Prog. Cosmic Ray Physics Vol. III, J.G. Wilson ed. (North
Holland Publishing Co., Amsterdam 1965) 1
[97] R.M. Barnett et al. (Particle Data Group), Phys. Rev. D54 (1996) 1
[98] K. Kamata and J. Nishimura, Suppl. Progr. Theoret. Phys. 6 (1958) 93
[99] M.F. Bourdeau et al., J. Phys. G: Nucl. Phys. 6 (1980) 901
[100] J. Nishimura, in: Handbuch der Physik, Vol. XLVI/2, S. Fl¨ugge ed. (Springer,
Berlin, 1967) 1
86
[101] F. Arqueros, S. Martinez and M. Rozanska, Proc. 23rd
Int. Cosmic Ray Conf.,
Calgary, 4 (1993) 738; S. Martinez et al., Nucl. Instr. Meth. A357 (1995) 567
[102] Handbook of Chemistry and Physics, 67th
Edition, R.C. Weast ed. (The Chemical
Rubber Co., Cleveland, 1986) E373
[103] J.N. Capdevielle, J. Phys. G: Nucl. Part. Phys. 16 (1990) 1539
[104] G.J. Alner et al., UA5 Collaboration, Phys. Lett. B167 (1986) 476
[105] R.E. Ansorge et al., UA5 Collaboration, Z. Phys. C43 (1989) 75
[106] G.J. Alner et al., UA5 Collaboration, Nucl. Phys. B291 (1987) 445
[107] J.N. Capdevielle and S. Zardan, Proc. 20th
Int. Cosmic Ray Conf., Moscow, 5
(1987) 160
[108] K. Alpg˚ard et al., UA5 Collaboration, Phys. Lett. B115 (1982) 71; G. Arnison
et al., UA1 Collaboration, Phys. Lett. B122 (1983) 189; G.J. Alner et al., UA5
Collaboration, Phys. Lett. B180 (1986) 415
[109] R.E. Ansorge et al., UA5 Collaboration, Nucl. Phys. B328 (1989) 36
[110] T. ˚Akesson et al., AFS Collaboration, Phys. Lett. B178 (1986) 447
[111] C. Geich-Gimbel, Int. J. Mod. Phys. 4 (1989) 1527
[112] A. Klar and J. H¨ufner, Phys. Rev. D31 (1985) 491
[113] C. De Marzo et al., Phys. Rev. D26 (1982) 1019
[114] M. Arneodo et al., Z. Phys. C31 (1986) 1
[115] R. Attallah et al., Proc. 23rd
Int. Cosmic Ray Conf., Calgary, 4 (1993) 48
[116] J.N. Capdevielle, Proc. 5th
Int. Symp. on Very High Energy Cosmic Ray Interactions,
Lodz, (1987) 1
[117] M. Aguilar-Benitez et al., NA27 Collaboration, Z. Phys. C50 (1991) 405
[118] R. Attallah et al., Proc. 22nd
Int. Cosmic Ray Conf., Dublin, 4 (1991) 157
[119] R. Hagedorn, Rev. Nuovo Cim. 6 (1983) 1
[120] G. Arnison et al., UA1 Collaboration, Phys. Lett. B118 (1982) 167
[121] G. Bocquet et al., Phys. Lett. B366 (1996) 434
[122] J.N. Capdevielle, Nucl. Phys. B (Proc. Suppl.) 52B (1997) 146
87
[123] G.J. Alner et al., UA5 Collaboration, Phys. Rep. 154 (1987) 247
[124] R. Attallah, private communication (1990)
[125] J.A. Capella, Proc. 22nd
Int. Cosmic Ray Conf., Dublin, 5 (1991) 15
[126] S. Jadach, Comp. Phys. Comm. 9 (1975) 297
[127] J.N. Capdevielle and T. Thouw, J. Phys. G: Nucl. Part. Phys. 18 (1992) 143
[128] D. Heck et al., Proc. 23rd
Int. Cosmic Ray Conf., Calgary 4 (1993) 60
[129] R.E. Ansorge et al., UA5 Collaboration, Z. Phys. C33 (1986) 175
[130] P. Bernard et al., UA4 Collaboration, Phys. Lett. B166 (1986) 459
[131] V. Innocente et al., Phys. Lett. B169 (1986) 285
[132] J.N. Capdevielle and J. Gawin, J. Phys. G: Nucl. Part. Phys. 12 (1986) 465
[133] G. Schatz, in: Report KfK 5027 (1992) 32, J.Knapp and H.Rebel eds., Kernforschungszentrum
Karlsruhe
[134] A. Capella, Report LPTHE 91/53 (1991), Orsay, France
[135] T. ˚Akesson et al., HELIOS Collaboration, Z. Phys. C49 (1991) 355
88
List of Tables
3.1 Computing times for various thinning levels . . . . . . . . . . . . . . 15
5.1 Basic features of the interaction models used . . . . . . . . . . . . . . 30
6.1 Decay modes and branching ratios for kaons . . . . . . . . . . . . . . 40
6.2 Coefficients of the parametrization of K −→ 3π . . . . . . . . . . . . 41
6.3 Coefficients of the parametrization of K −→ π + + ν . . . . . . . . 41
6.4 Decay modes and branching ratios for η . . . . . . . . . . . . . . . . 42
6.5 Decay modes and branching ratios for strange baryons . . . . . . . . 43
6.6 Decay modes and branching ratios for resonances . . . . . . . . . . . 43
7.1 Branching ratios for photonuclear reactions leading to two pions . . . 50
7.2 Parameters of the longitudinal age formula . . . . . . . . . . . . . . . 53
A.1 Parameters of the U.S. standard atmosphere . . . . . . . . . . . . . . 59
A.2 Parameters of the AT115 atmosphere (January 15, 1993) . . . . . . . 59
A.3 Parameters of the AT223 atmosphere (February 23, 1993) . . . . . . . 60
A.4 Parameters of the AT511 atmosphere (May 11, 1993) . . . . . . . . . 60
A.5 Parameters of the AT616 atmosphere (June 16, 1993) . . . . . . . . . 60
A.6 Parameters of the AT822 atmosphere (August 22, 1993) . . . . . . . . 61
A.7 Parameters of the AT1014 atmosphere (October 14, 1993) . . . . . . 61
A.8 Parameters of the AT1224 atmosphere (December 24, 1993) . . . . . 61
A.9 Air pressure values at low altitude for various atmospheres . . . . . . 62
C.1 Parameters of the hadron-nucleon cross section parametrization . . . 66
D.1 Energy dependence of x1 , x2 , and α . . . . . . . . . . . . . . . . . . . 74
D.2 Charge exchange of leading particles . . . . . . . . . . . . . . . . . . 75
D.3 Resonance formation of leading particles . . . . . . . . . . . . . . . . 76
89
List of Figures
2.1 Pressure difference of atmosphere relativ to U.S. standard atmosphere 7
3.1 Energy loss of muons as function of Lorentz factor . . . . . . . . . . . 10
4.1 Inelastic proton-air cross sections at higher energies . . . . . . . . . . 20
4.2 Inelastic nucleon-air cross sections at lower energies . . . . . . . . . . 22
4.3 Inelastic hadron-proton cross sections . . . . . . . . . . . . . . . . . . 23
4.4 Inelastic nucleus-air cross sections . . . . . . . . . . . . . . . . . . . . 24
4.5 Inelastic hadron-air cross sections at higher energies . . . . . . . . . . 25
4.6 Inelastic meson-air cross sections at lower energies . . . . . . . . . . . 26
7.1 Photoproduction cross section . . . . . . . . . . . . . . . . . . . . . . 49
B.1 Differences in muon range calculated with different approximations . 64
90