Table Of ContentELECTRON-PHOTON SHOWER
DISTRIBUTION FUNCTION
Tables for Lead, Copper and
Air Absorbers
BY
H. MESSEL and D. F. CRAWFORD
Φ
PERGAMON PRESS
OXFORD · LONDON · EDINBURGH · NEW YORK
TORONTO · SYDNEY · PARIS · BRAUNSCHWEIG
Pergamon Press Ltd., Headington Hill Hall, Oxford
4 & 5 Fitzroy Square, London W.1
Pergamon Press (Scotland) Ltd., 2 & 3 Teviot Place, Edinburgh 1
Pergamon Press Inc., Maxwell House, Fairview Park, Elmsford, New York 10523
Pergamon of Canada Ltd., 207 Queen's Quay West, Toronto 1
Pergamon Press (Aust). Pty. Ltd., 19a Boundary Street,
Rushcutters Bay, N.S.W. 2011
Pergamon Press S.A.R.L., 24 rue des Ecoles, Paris 5e
Vieweg & Sohn GmbH, Burgplatz 1, Braunschweig
Copyright © 1970 Pergamon Press Ltd.
First edition 1970
Library of Congress Catalog Card No. 69-16049
Printed in Great Britain by A. Wheaton & Co., Exeter
08 013374 6
CHAPTER 1
I CTION
EVER since Bhabha and Heitler (1937) suggested that With the advent of high-speed digital computers a
multiplicative processes of bremsstrahlung and pair new approach was possible. This was the direct simu
production played a dominant role in the propagation lation of the cascade by the Monte Carlo method.
of high energy electrons and photons through matter, Wilson (1952) had used hand simulation with the aid
much work, both theoretical and experimental, has of a "wheel of chance". However, it requires the
been done in deducing the characteristics of the cascade speed and flexibility of the digital computer to obtain
produced. relatively accurate and comprehensive results. Because
The theoretical work has been roughly divided into it is possible to simulate all aspects of the cascade, we
two branches, one dealing with the one-dimensional can use accurate cross-sections in an exact three-dimen
longitudinal development which considers the number sional model. Furthermore, the type of information we
development of shower particles, the other with the can extract from the simulated cascade is limited only
three-dimensional angular and radial spread of a by our ingenuity in programming the computer. The
cascade shower as it penetrates matter. disadvantage of the computer simulation is that even
In the past, two approximations, usually referred to on the fastest computers many hours are required to
as approximations A and B, have been commonly used get statistically significant results. In this manner the
in dealing with an electromagnetic cascade. In approxi simulation method is very similar to experimental
mation A, full screening cross-sections are used for observations of cascades.
bremsstrahlung and pair production; other processes In 1958 Messel and his collaborators commenced a
are completely neglected. Ionization loss, assumed to full-scale attack on the shower problem using the
be at a constant rate, is introduced for approximation digital computer SILLIAC. Results of some of this
B, but no other change is made. Even in the case of early Monte Carlo work appeared in a series of five
these two drastic approximations the diffusion equa papers, Butcher and Messel (1958, 1960), Messel et al.
tions for the simplest of the number and the angular- (1962) and Crawford and Messel (1962, 1965). The
radial distribution functions are of considerable results here reported are those obtained using SILLIAC
complexity and their solution of such a complicated and the much faster computer KDF9. They represent
nature that it is almost impossible to evaluate the results eight years of hard work, not by just the computers,
numerically. A review of the early work was given by but by a group of highly skilled staff.
Rossi and Griesen (1941) and a series of papers cover It is well known that the basic character of an electro
ing much of the later theoretical work on the number magnetic cascade depends mainly on two types of
distribution were prepared by Messel and his collabor interactions. When a high energy electron passes near a
ators—Janossy and Messel (1950,1951), Messel (1951, nucleus there is a finite probability that a secondary
1956), Messel and Potts (1952). The angular-radial dis photon will be produced. Classically we say that the
tribution problem likewise received much attention; electron has been accelerated in the Coulomb field of
see, for instance, Green and Messel (1952), Green et al. the nucleus with the emission of electromagnetic radia
(1952), Messel and Green (1953), Kamata and Nishi- tion. This process is called bremsstrahlung. Usually
mura (1958),Chartres and Messel (1954,1955,1956), or the photon has a very small energy so that a single
the review articles by Griesen (1956) and Belenki and electron will produce copious quantities of low energy
Ivanenko (1960) for further references. However, in all photons as it traverses an absorber. Though the nucleus
cases, there was a great lack of meaningful numerical is required to receive some energy for conservation of
results, thus making the interpretation of many cosmic- momentum, it receives very little, due to its relatively
ray and machine experiments difficult. large mass. The other major interaction is that of pair
2 ELECTRON-PHOTON SHOWER DISTRIBUTION FUNCTION
production. Here a photon materializes into an elec annihilation. Here a positron annihilates an atomic
tron and positron in the field of a nucleus. As in electron to produce two photons.
bremsstrahlung the nucleus is required to conserve These are not the only electron or photon inter
momentum but receives a negligible amount of actions, but they are the main ones and those which we
energy. have taken into account in our calculations. The others
If we consider only particles above a fixed cut-off such as the nuclear photo-electric effect, photo-meson
energy we see a rapid build up in particle numbers to a production, direct pair production by electrons and
maximum followed by a slower fall-off. From one high single photon positron annihilation have been left out.
energy electron primary the average secondary energy All of these have negligible cross-sections, except for the
decreases uniformly until it is comparable to the ther nuclear photo-electric effect. This may have a cross-
mal electron energy. section as high as 5% (Grodstein 1957) of the total
Of somewhat lesser importance are the inelastic col photon cross-section. However, since it occurs only in
lisions between high energy electrons and the atomic a narrow range of photon energies and does not con
electrons in the absorber. Since these collisions are tribute further electrons or photons to the cascade,
extremely numerous, with very little energy transfer to we have ignored it too. We shall return to a more
each atomic electron, we can consider them as pro detailed discussion of the cross-sections considered, in
ducing a continuous energy loss. Thus we essentially Chapter 3.
ignore the secondary electrons that are produced and In Chapter 2, we will discuss in some detail the
subtract from the high energy electron an energy method of simulating the electron-photon cascade on
proportional to its path length. To be more exact, we a digital computer. It is hoped that the description is
treat only those secondaries with energy less than some sufficiently complete to enable new programs to be
cut-off in this manner, and the continuous energy loss, written.
we call ionization loss. The interactions with energies The results presented in this work and discussed in
higher than the cut-off contribute to the total number Chapter 4 come from two programs. The first written
of electrons in the cascade. These interactions are for the computer SILLIAC has been described in
commonly called inelastic electron scattering. Both the Messel et al. (1962). This program used exact cross-
ionization loss equation and the inelastic scattering sections for bremsstrahlung, pair production and the
cross-sections are different for positrons and electrons. Compton effect. Since the low energy cut-off was 10
Photons also interact with the atomic electrons. MeV, no low energy corrections to these cross-sections
Compton scattering is the collision between a photon were applied. Inelastic electron scattering was included
and a free electron. We can usually consider the absor into the ionization loss expression. Furthermore the
ber electrons to be free because the energy transfer to photo-electric effect and positron annihilation were
them is much higher than their ionization potential. ignored, although exact treatment of multiple Coulomb
A different type of interaction is elastic \Coulomb scattering was included. Within the range of primaries
scattering. Here we consider the elastic scattering of an less than 2000 MeV and secondaries greater than 10
electron in the Coulomb field of a nucleus or of an MeV for which results are given, all of the omissions
electron. Since these scatterings are quite small and cited have a negligible effect with perhaps an exception
very numerous they are treated as a continuous inter at the 10 MeV level.
action using the theory of multiple Coulomb scattering. The second program was written for a much faster
Although there are angular deflections in the other computer, the English-Electric-Leo KDF9. This pro
interactions mentioned, multiple Coulomb scattering gram includes all the low energy corrections and cross-
dominates the spread of the cascade. It may be safely sections together with some very high energy correc
ignored in light absorbers, where usually only one- tions to bremsstrahlung and pair production. As a
dimensional results are required, but cannot be ignored result, it can accurately simulate cascades for any
in heavy absorbers such as lead or where three-dimen primary energy and for secondaries as low as 1 MeV in
sional results are needed. any absorber. Since this new program supersedes the
Of still less importance is the photo-electric effect in SILLIAC program already described in Miessel et al.
which a low energy photon transfers its energy to a K (1962), it is the only one discussed in the subsequent
or L shell electron. Of similar importance is positron chapters.
CHAPTER 2
SIMULATING THE
ELECTRON-PHOTON CASCADE
SINCE the electron-photon cascade is a stochastic pro where (/, m, n) are the direction cosines of the particle.
cess and hence can be described completely in statistical The selection of the distance t depends on the
terms, it readily lends itself to simulation on a digital cross-sections so we now consider in more detail the
computer. The process of simulation involves setting representation of a cross-section as a probability
up a computer model of the cascade by performing distribution and the sampling of random variates from
operations which represent the motions and interac such a probability distribution.
tions of the cascade particles. Since the correspondence The usual differential cross-section d</>(E , E E ) for
0 l9 2
between the model and the actual physical process is a a primary of energy E and secondaries of energy
0
very close one, we shall describe the model in terms of E and E = E —E may be treated as a combination
x 2 Q 1
the physical processes themselves. of two probability distributions. The first is the prob
The basis of the model is that physical cross-sections ability distribution of distance which the primary will
are directly represented by probability distributions travel until it has an interaction. The second distribu
and the use of random numbers enables us to sample tion, which is proportional to αφ(Ε , E E ), is the
0 u 2
from these distributions in order to simulate the probability of producing a secondary of energy be
cross-sections. Techniques which use random numbers tween E and E + dE given that an interaction has
1 x x
in this way to simulate statistical processes are com occurred. Because the probability that an interaction
monly called Monte Carlo methods. will take place within a distance element dt is indepen
The probability distributions mentioned above have dent of Mn a homogeneous absorber, the probability
parameters determined by the energy of the incoming that an interaction takes place between t and t + dt is
particle and the composition of the absorber. The given by
particles are represented inside the computer by a set
/(/) dt = φβ-φί dt
of numbers giving the type, energy, weight, position in
Cartesian coordinates, and orientation in direction where φ is the total cross-section. That is
cosines of the particle. The absorber is regarded as a
semi-infinite flat slab normal to the z-axis. The primary φ = άφ(Ε , E E )
0 l9 2
electron or photon is injected along the z-axis. Since
the control routines can only handle one particle at a
where the integral is taken over all allowed secondary
time any more must be stored temporarily. Whenever
energies. The probability distribution for the secondary
the program finds that it needs to be tracking two
energies is
particles it chooses to store the one with higher energy
in a last-in-first-out stack. This choice ensures that the f(E) dE = (Iff) άφ(Ε , Ε E ).
x x 0 19 2
maximum number of particles in the stack is less than
From this equation it would appear that we need to
or equal to the next integer greater than log (E /E),
2 0 calculate the total cross-section φ before we can sample
where E is the largest and E the smallest particle
0 variates from this distribution. Fortunately this is not
energy which can occur. The moving of the particles is
necessary. Before we explain this paradox we need to
done by updating the coordinates, for example, to
consider the selection of random variates from statisti
move a straight line distance / without changing direc
cal distributions.
tion the new coordinates (χ', y\ z') are given by
The ability to select random variates within a digital
x' = It + x, y1 = mt + y, z' = nt + z computer depends on the availability of a set of
4 ELECTRON-PHOTON SHOWER DISTRIBUTION FUNCTION
random numbers. We are not interested in abstract However, by taking the larger of two random numbers
concepts of randomness, but only require that the we also obtain a variate x from this distribution usually
numbers shall be mutually unbiased and independent in less time than the square root would take. This
of the problem at hand. Most computers now have example suggests that the design of a sampling pro
routines which produce in practice an unlimited num cedure is very dependent on the idiosyncracies of the
ber of random numbers one at a time. Each number required distribution. There is a method which is a
depends uniquely but in a very obscure way upon its combination of "the composition method" and "the
predecessors. By using the same starting values the rejection methods", both described in Butler (1956),
sequence may be repeated—a facility which is of con which is of considerable power.
siderable value in code checking a program. In actual If f(x) can be written
production runs care is taken not to repeat identical
sequences. We used a method of generating random n
numbers due to Coveyou (1960) and tested by Ruehn Ax) = /*ift(x)glx)
t
(1961). This method provides random uniform num
bers, within the range 0 < ξ < 1. The symbol ξ will
always be used to denote such a variate. where a, > 0 / = 1, 2, . . ., n
Let us consider the frequency distribution f(x) dx
0 ^ g (x) < 1
with Cumulative frequency distribution F(x) i.e. t
9
X ft(x) > 0 for all x
F(x) = J7(x) dx.
and fi(x) dx = 1
0
A theorem which makes it possible in principle to
write down a rule for computing variates from any then a random variate x may be chosen by
distribution is that the cumulative distribution F(x) is
itself distributed uniformly. If x is sampled from/(x), (a) select a random integer / (1 < / ^ n) with
the probability that x will be less than xf is F(x'), and probability of selecting / being proportional to
since F(x) is a monotonic increasing function of x we
have (b) select a variate x' from the frequency distribu
tion/,(*),
Prob [F(x) < F(x')] = Prob [x < x']
(c) calculate gi(x') and accept x = x' with prob
= F(*'). ability gi(x')9 if rejected repeat from step (a).
Therefore one method of sampling x is to compute x The value of the method lies in the fact that a de
from composition of the frequency function can usually be
F(x) = ξ 0 ^ ξ < 1 found where n is a small number, the distributions ft(x)
are easily sampled, and gi(x) is never much less than one,
or x = F-1 (0.
this means that a rejection is a rare event. It is of in
For example if terest to note that the probability of accepting any
trial is
F(x) = 1 — exp (—x/x)
n
then x = — x In (1 — ξ)
72-
= x In ξ.
The last step is possible because 1 — ξ is a uniform
Thus the average number of attempts is
variate with the same limits as ξ. This example
gives us the most useful method of sampling from the
exponential distribution. However, this method is not
much help when the cumulative distribution does not
have a simple inverse and even when it has, some other <=1
methods may be better. For instance to sample from
Examples of distributions easily sampled are those
the distribution
which have simple inverses and those which have the
f(x) dx = 2xdx 0 < x < 1 form
the inverse method would have us take
f(x) = (M + a + l)!^ (i _ y o ^^\
x x
ml n\
SIMULATING THE ELECTRON-PHOTON CASCADE 5
where m and n are positive integers. To sample x we we can avoid the calculation of the total cross-section
select m + n + 1 random numbers and order them, φ since we may use
9
i.e.
n αι
2
£l < £2 ^ ?8 · · · ^ im+if+l
then we select for x the number ξ„ .
+1
instead. Then whenever a variate is rejected we
As an example of the use of this sampling method
sample a new distance and move the particle before
and because it is of general interest let us consider the
another variate is sampled. Thus the particle may be
normal distribution
moved several times before a successful interaction is
simulated.
/(*) ■Ae ~x2/2 x>0. We will now prove that this gives the correct cross-
section. Suppose that the correct total cross-section is φ
and that
We need only consider the reduced form and positive
values as it is a trivial matter after we have a sample
variate to perform a scale change and provide a ran
/*t = Φ\
dom sign. A decomposition due to Butcher (Thesis)
which is very efficient is
then clearly the probability of accepting a trial is φ/φ'.
=J\
Thus if k attempts are made before acceptance, the
/w 1 -*2/2
distance t travelled will be distributed as
(ΦΎ
-f- —- 2 e~2ix~1} £-(*-2)2/2 -ΦΊ
V27T " (k - 1)!
where But the probability of having exactly k attempts is
{ΦΙΦ'ΧΙ—ΦΙΦ'Υ-1. Hence the distance travelled will be
- =J! 1
distributed as
> 0
00
Μχ) = 1 Σ (<L·)l·tk-1e-''", 6
W' e Σ. (1 - φ/φγ-
(x) = e-*i*
gl
(k-iy. φ'
and
"Σ
<f>e~
1
a = —
2
V277
/,(*) =2e-2"-« > JC ^ 1. = φ e -4>t
g2(x) = e-«*-«2/* J which is the correct distribution. The method is readily
generalized to the case where there are two or more
In this case the decomposition corresponds to two competing interactions.
3
ranges of x. The value of a + a is ,— ~ 1.2 which In considering a particular cross-section it may well
x 2
be that the decomposition used has values of
is the average number of trials which must be made. It
will be noted that this method of generating random na i
2
normal deviates is exact within the limitations and
accuracies of the computation.
We now come to a point which is of considerable
value in simulating physical cross-sections. If we ex dependent on the primary energy. In this case care
press the differential cross-section in the form must be taken if the primary is continuously degraded
in energy by a process such as ionization loss. Further
Σn care is necessary if one is dealing with an inhomogen-
<**/*(» giM eous absorber. A useful technique is to introduce a
pseudo cross-section which does nothing except reduce
6 ELECTRON-PHOTON SHOWER DISTRIBUTION FUNCTION
the average step length. Whenever this pseudo cross- error is that position of the terminal point may be
section is selected the program immediately returns to slightly wrong. Furthermore, it is in just this low energy
move another step and in so doing recalculates region that limitations on the maximum distance t are
required for multiple Coulomb scatterings. Thus we
have used only the simple form for calculating ioniza
n
tion loss. It should be noted however that this simpli
fication is not necessarily valid in the general case.
i= 1 The other interaction which is best dealt with as a
In the present program because a limit has been placed continuous process is elastic Coulomb scattering of the
on the maximum distance δ, which electrons may electron in the field of the nucleus. It is possible to
move, this technique has not been used. When a dis simulate these scatterings directly as normal cross-
tance t is sampled the electron is moved in steps of δ sections. However, such a scheme would be prohibit
and then the remaining distance (less than δ) to give ively expensive because of the large cross-sections
a total distance t. The reason for introducing a maxi involved, and since adequate theories are available
mum step δ will be returned to later in the discussion which predict the angular distribution of the scattered
on multiple Coulomb scattering. electron after a path length t, we have adopted this
Not all interactions are conveniently represented as approach. Given that the electron with position co
cross-sections. Ionization loss and multiple Coulomb ordinates (x, y, z) and direction cosines (/, m n) is to
9
scattering are best simulated as continuous functions be moved a distance /, we move this distance in its old
of the electron track length. In general the rate of direction and then sample a spatial angle Θ through
energy loss per track length may be expressed as which it is scattered. The distribution of Θ depends
on the electron energy and on the distance t. The
dt angle 0 is measured from the current direction of the
electron. We now have to sample a random azimuthal
where the function/(is) depends only on the energy E.
angle and then calculate the new direction cosines.
Now provided the distance t is small the energy loss
Let μ be the original direction of the electron and
Δ£ incurred by an electron of initial energy E is
0
μ' = (/', m' ri) be the direction after scattering. Then
approximately 9
if a and β are two unit vectors which form an ortho
gonal triad with μ we have
μ' = μ cos© + (Ca + Sß)s'm0,
The error in using this approximation may be esti
where C and S are the cosine and sine of a random
mated from the more accurate expression
angle. These may be found easily by using a method
suggested by von Neumann (1951). Select a pair of
Kdt/Eol 2\dE dtj* \ random uniform variates ξ and ξ such that
0 χ 2
Now as will be seen later the ionization loss can be fi2 + f.1 < 1
expressed by
then
_ ii2 ~ &'
c
fi2 + f.1
where energies are in mc2 units, distances are in radia and
tion lengths and typically b ~ \a and γ ~ J. The
correction term in this case becomes
\ybt E-v-i
0
where S is given an arbitrary sign. The most conven
which is negligible at high energies. At low energies we
ient forms for a and β that have been found are
can show that the correction is less than
« = ( - ^ - - 1 , ^ , 1)
\n + 1 /i+l /
l Im m2 \
This result follows from the observation that the maxi βQ = , — 1Λ, /w 1.
mum distance that the electron can travel is that which
would reduce its energy to zero. Obviously the error is Because of ionization loss we must allow for the
a maximum at low energies where the required accur change in energy of the electron as it traverses the
acy is least. In this region the only result of such distance t. We do this by replacing E2 by E E, β2 by
0
SIMULATING THE ELECTRON-PHOTON CASCADE 7
ßß, etc., where E and β are the initial and £and ß the results. For 5 MeV electrons the correction to the path
0 0 0
final energies and velocities respectively. This substitu length for an average step is only 5% and is much
tion will give us the correct mean square scattering smaller for high energy electrons.
dE There is one further complication to the moving of
angle if the rate of energy loss, — is constant (see
particles and this is because we need to get information
dt
out of the program. As previously mentioned we con
Rossi, 1952, p. 68). A source of error results from the
sider the absorber to be a semi-infinite horizontal flat
assumption that all the scattering takes place at the end
slab with the primary injected vertically to the surface.
of the track. This error is not at all serious provided
We then count the number of secondaries as they pass
the scattering angle is small. Similarly the lateral dis
through a set of horizontal levels at regular depths in
placement due to scattering has been neglected. Both
the absorber. For complete tabulation we note the
these errors are negligible if the maximum step length
exact type, energy, radial distance and angle of each
is kept small enough at low energies where scattering
secondary that traverses each level. The radial distance
is greatest.
is measured from the core of the cascade, defined as the
A related effect has been included in the simulation
extrapolation of the primary track, and its angle is
program. This effect is the increase in path length over
measured relative to the primary direction. There is
the straight line distance due to multiple Coulomb
no need for any other spatial or angular measurements
scattering. It has been shown by Yang (1951) that if t
as theory predicts complete rotational symmetry.
is the total distance travelled, it is related to the straight
To obtain the required data for each secondary the
line distance s by
simulation program must note whenever a particle
would cross a level. Thus every time a particle is moved,
t = + i S\t) dt the distance (along the particle's direction) to the next
s
level is calculated and if the particle would cross the
level it is first moved to the level, recorded and then
where 0\t) is the mean square angle of scattering and moved the remaining distance. It is of interest to note
is a function of the distance t. Since we are only estima that because of the exponential distribution, we could
ting a correction we may use a rough estimate for quite validly forget the remaining distance and re-
02(f). Following Rossi (1952, p. 68) we have sample a new path length. However, it is usually faster
to remember the remaining distance.
W) =— t
9 So far no mention has been made of particles which
ß*E2 are travelling backwards. There is no trouble in simu
where E is the scattering energy, independent of lating these particles provided we do not wish to record
s
absorber if t is measured in radiation lengths and is them. However, once we introduce levels there are
equal to 42 mc2. On integrating we get complications. Since a particle may circle around and
cross a level several times, a decision has to be made
about how to record them. The decision adopted was to
\β*Ε) track all particles in the backwards direction until they
had crossed a level, after being recorded they are
or
tracked no further. This decision it is felt corresponds
most closely to experimental situations. For instance,
L 4 \ß2Ej J in a multiplate cloud chamber, consider one of the
gaps. The program would simulate all forward going
This equation forms the basis for our limitations on the
electrons coming out of the top plate and all backward
maximum step length allowed as mentioned earlier.
going electrons coming out of the bottom plate, but
We require that the bracketed expression shall always
not any electron which emerged from the bottom plate
be greater than \. This gives the requirement
entered the top plate and then re-emerged in the for
ward direction. Not only is the number of such elec
trons very small but in general their numbers would be
altered by the intervening absorber. A further conse
Since the lowest energy to be considered is 1 MeV ~ 2 quence of our treatment is that backwards going parti
mc2, the smallest value of t is 0.001 radiation lengths. cles do not enter a plate from below and emerge from
m
Even this expression for t„, which gives a root mean . its top surface. The results of the simulation will then
square angle of scattering of approximately 80° is too tabulate at each level the number of particles as a
large for all the approximations made. But since these function of Θ, their angle to the primary direction.
very low energy electrons remain in quite a small Forward going particles have 0 < Θ < 90° and back
region until they lose all their energy no serious error ward going particles which come from the absorber
8 ELECTRON-PHOTON SHOWER DISTRIBUTION FUNCTION
below the level have 90° ^ Θ ^ 180°. The radial and factors involved. To prevent large fluctuations which
energy distributions are given only for forward going could arise if one particle achieved a very large weight
particles. Track lengths and energy losses are tabulated we limited the maximum weight to 8. Once a particle
every time a particle moves, and refer to the absorber achieved this weight it was no longer subject to the
between two levels, hence backward going particles are sampling scheme.
included. One necessary requirement of any modification to
One of the problems with direct simulation of the the direct simulation of the cascade is that it does not
electron photon cascade is that while a large number of introduce any bias in the final results. An essential part
particles occur near the cascade maximum and hence of our importance sampling scheme is that the prob
results for this region have high statistical accuracy abilities for duplication and rejection are fixed and that
there are many fewer in the cascade tail with low statis each sampling process is completely independent of all
tical accuracy. We would like to have relatively the others. This means that we only need consider a single
same statistical accuracy for all results derived from sampling process in order to see if there is any bias.
the simulation. Since the number of particles going into Since we are mostly interested in energy and number
the many cells of the radial and angular distributions averages all we need to show for no bias in these cases
are both highly variable and mutually correlated we is that averages taken immediately after a level have the
shall consider only what can be done about the same expectations as those taken before. Because the
one-dimensional distribution. What we would like is scheme does not alter the individual energies, positions
to simulate fewer particles at the cascade maximum and or angles of the particles the averages will still be put
more in the cascade tail. Furthermore we would like into the correct distribution cells. For the high energy
to simulate more high energy particles and fewer low duplication scheme it is apparent that the before and
energy particles. As a rather crude attempt to achieve after averages are in fact identical. Therefore this part
these aims and save computer time a scheme of impor of the scheme does not introduce any bias. For the low
tance sampling has been implemented for the KDF9 energy rejection scheme the two averages also have
portion of our calculations. the same expectation. This was achieved by dividing
This importance sampling scheme gave each particle the weight of surviving particles by their survival
a weight which could only be altered when the particle probability.
crossed a recording level. By correct control of the The effect that the importance sampling scheme has
particle weights we can either increase or decrease the on the statistical accuracy of the results has to be
actual number of particles simulated but still get correct considered. Except for the first few levels where the
estimates of the average energies and numbers. The scheme has little effect its main result is to increase the
scheme makes a distinction between high energy and number of high energy particles simulated and to
low energy particles; only high energy particles may be decrease the number of low energy particles simulated.
duplicated and low energy particles may be retained or But since the low energy particles are the products of
rejected but not duplicated. We have taken the demar- the many high energy collisions it can happen that
kation energy as one tenth of the primary energy. Two there is an overall increase in the number of particles
preset parameters p and p control the operation of simulated compared to cascades in which the sampling
x 2
the scheme. The parameter ρ is the probability that a scheme is not used. We shall return to this point in
Ύ
high energy particle will be duplicated on traversing more detail in Chapter 4 in the discussion of errors.
a level. The duplication is done by making an identical We will also see in Chapter 4 one disadvantage of the
copy of the particle in the current particle stack and importance sampling scheme: we can no longer have
giving to each one, half the weight of the original an exact accounting of the energy distribution. If a
particle. For low energy particles the parameter p is greater than average number of low energy particles
2
the probability that the particle would be allowed to survive it is possible for there to be more energy carried
continue past the level, in which case its weight was by the total of particles than the original primary
divided by p . Particles which failed this selection were energy. Thus at a given level for p = 1/2 it would
2 2
rejected by removing them from the current particle mean that more than one half of the particles with
stack. Because of the cumulative properties of the E < 0.1 E had survived. Since these particles' weights
0
duplication scheme at successive levels, the time taken are all greater than one, they would contribute more
to simulate a cascade is a very sensitive function of the than the original total energy of the low energy particles,
duplication probability p. We chose the values p = even though as noted previously the average energy
x x
0.2 and p = 0.5 as a reasonable compromise of all expectations still remain the same.
2
