+-------------+ +----------##
| Geant 3.15 | GEANT User's Guide | PHYS010 ##
+-------------+ +----------##
Author(s) : M. Maire Submitted: 20.12.84
Origin : GEANT2 Revised: 23.08.91
The simulation of the processes which affect the propagation of a particle through the matter of the detector (e.g. bremsstrahlung, delta-rays production, Compton scattering and so on) is performed by GEANT in the following steps:
The mean free path of a particle for a given process, lambda, depends on the medium and cannot be used directly to sample the probability of interaction in a heterogeneous detector. The number of mean free paths which a particle travels
N =int ((dx)/ (lambda(x)))
lambda
is, on the other hand, independent of the material traversed. If N is a
R
random variable denoting the number of mean free paths from a given point
until the point of interaction, it can be shown that N has the
R
distribution function
-N
P(N
The total number of mean free paths to travel before the interaction
point, N , is sampled at the beginning of the trajectory as:
lambda
N =-log (eta)
lambda
where eta is a random number uniformly distributed in the range (0,1).
N is updated after each step Deltax according the formula:
lambda
N' =N -((Deltax)/(lambda(x)))
lambda lambda
until the step originating from s(x)= N lambda(x) is the shortest
lambda
and this triggers the specific process.
Common GCPHYS
The variables described above are stored in the common
/GCPHYS/:
COMMON/GCPHYS/IPAIR,SPAIR,SLPAIR,ZINTPA,STEPPA
+ ,ICOMP,SCOMP,SLCOMP,ZINTCO,STEPCO
+ ,IPHOT,SPHOT,SLPHOT,ZINTPH,STEPPH
+ ,IPFIS,SPFIS,SLPFIS,ZINTPF,STEPPF
+ ,IDRAY,SDRAY,SLDRAY,ZINTDR,STEPDR
+ ,IANNI,SANNI,SLANNI,ZINTAN,STEPAN
+ ,IBREM,SBREM,SLBREM,ZINTBR,STEPBR
+ ,IHADR,SHADR,SLHADR,ZINTHA,STEPHA
+ ,IMUNU,SMUNU,SLMUNU,ZINTMU,STEPMU
+ ,IDCAY,SDCAY,SLIFE ,SUMLIF,DPHYS1
+ ,ILOSS,SLOSS,SOLOSS,STLOSS,DPHYS2
+ ,IMULS,SMULS,SOMULS,STMULS,DPHYS3
+ ,IRAYL,SRAYL,SLRAYL,ZINTRA,STEPRA
*
The first 9 processes (from PAIR production up to MUon NUclear
interaction) and the RAYLeigh scattering have the same scheme. Let's take
as an example the first one. For a complete description of the common see
[BASE 030].
- SLPAIR
- last point along the trajectory where the
interaction happened for the current particle. In GEANT now only SLPAIR
(pair production for mu), SLDRAY (delta-ray production), SLBREM
(bremsstrahlung for mu), SLHADR (hadronic interactions), SLMUNU (muon
nucleus interactions) and SLRAYL (Rayleigh effect) are used in the current
version of GEANT.
- ZINTPA
- N =remaining number of mean free paths
lambda
[4] evaluated at the last point where the mechanism was active, i.e.
IPAIR!=0.
- STEPPA
- lambda(x)= value of the mean free path at
the last point where the mechanism was active, i.e. IPAIR!=0;
- SPAIR
- N lambda(x) = remaining step-length
lambda
before interaction, evaluated at the last point where the mechanism was
active, i.e. IPAIR!=0.
- IPAIR
- flag for secondaries:
- 0=
- the process is turned off.
- 1=
- generation of secondaries enabled.
- 2=
- no generation of secondaries.
The evaluation and update of the quantities like STEPPA, SPAIR and ZINTPA
is turned off in the media where the mechanism is not active (IPAIR!=0).
Turning off a mechanism in one tracking medium may give incorrect physics
results because not only the mechanism will not be active, but the
interaction probabilities will not be updated, as if that medium would not
have been traversed at all. This feature of the tracking routines is used
mainly in the vacuum, (defined as a medium with atomic number Z<1), where
all the mechanisms but DECAy are automatically switched off by GEANT
itself.
The DECAy in flight is simpler since the mean life time of the particle,
tau, is not material dependant and can be sampled directly.
- SLIFE
- not used.
- SUMLIF
- proper time left before the decay. At the
beginning of the track SUMLIF= -ctau log (eta).
- SDCAY
- distance left to decay point evaluated at the
last point where the mechanism was active, i.e. IDECA!=0.
The cross-section and energy loss tables
Cross sections and energy loss dE/dx are tabulated for all materials
referenced by a tracking medium under control of the routine GPHYSI. The
values of the energy for which the tabulated quantities are calculated are
stored in the common
/GCMULO/ (see [BASE030]).
To evaluate one of the tabulated quantities for a particle of kinetic
energy E , a linear interpolation is used. Let i be such that:
0
E
Integer IEKBIN in common
/GCTRACK/ is equal to i
during tracking and its value is recomputed by the routine
GEKBIN when the energy of the
particle changes. If the quantity Y has been tabulated so that Y = Y(E )
i i
then the value Y = Y(E ) is calculated as:
0 0
Y =Y +((E -E )/(E -E ))(Y -Y ) = Y (1-((E -E )/(E -E )))+Y ((E -E )/(E - E ))(1)
0 i 0 i i+1 i i+1 i i 0 i i+1 i i+1 0 i i+1 i
Inside the code the following quantities are used:
GEKRAT = ((E - E )/(E - E ))
0 i i+1 i
GEKRT1 = (1- ((E -E )/(E -E )))
0 i i+1 i
where GEKRAT is in common
/GCTRAK/ and GEKRT1 is a
local variable recomputed when needed.
The energy loss tables
Unlike previous versions of GEANT, energy loss and multiple scattering
continous processes are applied at every step for charged particles in
matter (Z>1). The variables SOLOSS, STLOSS, SOMULS, STMULS in common
/GCPHYS/ are obsolete.
They have been kept for backward compatibility, but their use may give now
unpredictable results. The common will be rearranged in a future version
of GEANT.
As explained in [PHYS330] and [PHYS430], energy loss tables are calculated
at initialization time (GPHYSI)
for all the materials which are referenced by a tracking medium. These
tables contain NEK1 values of dE/dx calculated for the corresponding
values of energy in ELOW.
Limitations on the step size
The routine GMULOF called by
GPHYSI creates and fills a table
of NEK1 values corresponding to the ELOW values containing the smaller
upper limit for the step among those imposed by three continuous
processes: energy loss, multiple scattering and magnetic field bending of
the track.
Continuous energy loss can introduce an upper limit on the step via the
variable DEEMAX, an argument to the
GSTMED routine. During tracking
the value of DEEMAX for the current medium is stored in the common
/GCMATE/. DEEMAX is the
maximum fraction of its kinetic energy which a particle can lose in a step
due to continuous ionization (0
step< =((DEEMAX)/(dE/dx))
Multiple scattering as well can limit the step size, see [PHYS320] and
[PHYS325]. The limitation is given as:
step< =min(T , 10X )
Bethe 0
where X is the radiation length.
0
Another upper limit on the step size comes from the magnetic field. The
bending of the particle in magnetic field may be limited by the TMAXFD
argument to the GSTMED routine.
During tracking the value of TMAXFD for the current medium is stored in
the common /GCMATE/.
A lower limitation on the tracking step is not generally imposed. There is
however a protection against the step being reduced to a very small value
by continuous processes. In particular multiple scattering at low energies
(<1 MeV) can impose a very small tracking step with serious consequences
on the tracking time. To avoid this, a lower limit on the step imposed by
continuous processes is introduced: STMIN. The meaning of STMIN is the
following: below 1 Mev the stopping range is usually small. If the
stopping range becomes smaller than STMIN, the constraint imposed by the
multiple scattering is ignored and the minimum is taken between the
reduced stopping range (the distance the particle has to travel to reach
its threshold energy) and STMIN itself. In this sense STMIN is no more
than an accelerator of the tracking for stopping particles.
Another limitation on the step size which is imposed by the GT...
routines at tracking time is the STEMAX parameter of
GSTMED which sets an absolute
upper limit to the size of a traking step for each tracking medium.
Automatic calculation of parameters
The definition of a tracking medium requires the specification by the user
of a set of parameters (see [CONS200]) which can deeply affect the
tracking and hence the physic results of the GEANT Monte-Carlo. To help
the user to find the optimal set of parameters, by default GEANT overrides
the values of STMIN, DEEMAX, STEMAX and TMAXFD. This behaviour is
controlled by the AUTO data record and interactive command. By default
AUTO=1 and automatic evaluation of the parameters is enforced (seen
further down for the partly anomalous behaviour of TMAXFD. When AUTO=0
then only the parameters whose value is <0 are recalculated by GEANT,
while for the others the user input is accepted with minimal checking.
When the automatic calculation of parameters is active, the following
applies:
- DEEMAX
- The formula used by GEANT is the following:
ISVOL =0 andx <2cm 0.25
ISVOL =0 andx0> =2cm 0.25- ((0.2)/(sqrt(x )))
DEEMAX = { 0 0 .
ISVOL! =0 ((0.2)/ (sqrt(x )))
0
ISVOL>0 defines a sensitive detector, while ISVOL<=0 is a non-sensitive
detector (see [HITS]).
- STMIN
- The default value corresponds to a stopping
range of 200 KeV above CUTELE.
o
- TMAXFD
- The default value corresponds to 20 if the
input value to GSTMED is <=0 or
o
if it is >20 .
10
- STEMAX
- The default value is 10 (BIG variable in
common /GCONST/).
These values have been tuned empirically on a variety of setups and users
are invited to start with automatic computation. The values of the
parameters can be checked with a call to GPRINT('TMED',0) after
GPHYSI has been called. Use of
AUTO mode is strongly recommended because it makes GEANT a predictive
tool. In this way all parameters are automatically computed by the system
as opposed to tuning data and Monte-Carlo via the tracking parameters.
Energy loss
Energy loss tables
In previous versions of GEANT (up to 3.13) the energy lost by a charged
particle travelling through matter was calculated by the formula:
DeltaE =((dE)/ (dx)) xstep
# #
The dE/dx quantity was calculated and tabulated for e , mu and protons
and the value for a particle of a given energy E was calculated
0
interpolating linearly the tables:
E
. ((dE)/(dx))| = .((dE)/ (dx))| +((E -E )/(E -E ))(. ((dE)/(dx))| - .((dE)/(dx))| (2)
E=E E=E 0 i i+1 i E=E E=E
0 i i+1 i
= GEKRT1x.((dE)/ (dx))| +GEKRATx.((dE)/ (dx))|
E=E E=E
i i+1
These formula works of course in the assumption that the error is
2
negligible. But the error on the linear approssimation is O(d E/dxdE)
which is far from being negligible at low energy. In this way the energy
loss at low energies was constantly underestimated and the particles could
travel too much before stopping.
Stopping range tables
To correct for this problem, a different approach was introduced in GEANT
version 3.14. The stopping range of a particle is defined as the distance
that the particle will travel before stopping. By definition the stopping
range for a particle of energy E is given by:
0
E
0 0
R =int ((dx)/(dE))dE = int -((dx)/ (dE))dE
E 0
0
Note that in the tables the positive quantity -dE/dx is stored. The
method used was to build a table of stopping ranges based on ELOW by
integrating the inverse of the dE/dx tables in
GRANGI. At tracking time the
algorithm was the following:
- Evaluate the stopping range for the threshold energy (STOPC). This
only once at the beginning of each new track.
- From the energy of the particle derive the stopping range by a linear
interpolation of the range table:
R =GEKRT1 xR +GEKRAT xR
0 i i+1
where E Evaluate the stopping range for the particle after the step:
R' = R -step. If this is less than the stopping range of a particle with
0 0
threshold energy, the particle is terminated as a stopping particle below
the energy cut. Otherwise the following quantities are evaluated:
j such that R
and the final energy is computed as:
E' =GEK1 xELOW(j) +GEK xELOW(j+1)
0
- the energy loss is computed as:
DeltaE =E - E'
0 0
This value is then corrected to take into account the Landau fluctuations
(see [PHYS332]).
This method has two main disadvantages. The first is due to the finite
precision of computers. As the energy loss in a step is calculated as the
difference of two numbers, it is subject to large relative errors. The
effect can be particularly large in case of light materials, particles
near the minimum ionization point or very short steps, where
DeltaE= DESTEP can even result in a negative quantity. As the relative
-6
precision of a 32-bit computer is around 10 , the error on the energy
loss of a 100GeV track can be as large as 100KeV.
The second problem connected with this method can be easily shown if we
compute dE/dx:
((dE)/ (dx)) = ((d)/(dx))(E -E' ) = ((d)/ (dx))(E -GEK1 xELOW(j)-GEK xELOW(j+1))
0 0 0
= ((d)/(dx))(E -E - ((R' -R )/ (R -R ))(E - E )) =-((d)/ (dx))((DeltaE )/(DeltaR ))(R' - R )
0 j 0 j j+1 j j+1 j j j 0 j
= -((DeltaE )/(DeltaR ))((dR' )/(dx)) = ((DeltaE )/(DeltaR ))((step)/ (dx)) = ((DeltaE )/(DeltaR ))
j j 0 j j j j
As we can see the reconstructed dE/dx curve due to continuous energy loss
is a step function constant in each energy bin. This, although the results
obtained with GEANT 3.14 were very satisfactory, this was felt to be an
undesirable feature.
Energy loss in GEANT 3.15
The two problems mentioned above have been solved in GEANT 3.15. As far as
the precision is concerned, the solution was to revert to the algorithm of
GEANT 3.13 every time the relative energy loss in the step (DESTEP/GEKIN)
is smaller than five times the machine precision. This has given good
results without loosing the substantial improvements obtained in GEANT
3.14 with the introduction of the stopping range algorithm. As a matter of
fact the above condition happens only in case of very small step or when
the dE/dx curve is very flat, and in both cases the algorithm based on
the dE/dx tables is a good approximation.
The second problem has been solved changing the interpolation algorithm,
going from a linear to a quadratic interpolation. Instead of supposing a
linear relation between energy and stopping range in every energy bin, we
suppose a quadratic relation of the sort:
2
E = f(R) = aR +bR +c
The only problem is now the determination of the coefficients a, b and c.
To do this we recall that the general formula of the parabola which passes
through the points (x ,y ),(x , y ),(x , y ) is the following:
1 1 2 2 3 3
y = (((x- x )(x-x ))/((x -x )(x -x )))y +(((x- x )(x-x ))/ ((x -x )(x - x )))y + (((x-x )(x-x ))/((x -x )(x -x )))y
2 3 1 2 1 3 1 1 3 2 1 2 3 2 1 2 3 1 3 2 3
2
= x (((y )/((x - x )(x -x ))) +((y )/ ((x -x )(x - x )))+ ((y )/((x - x )(x -x ))))-
1 1 2 1 3 2 2 1 2 3 3 3 1 3 2
x(((y (x +x ))/((x -x )(x -x )))+((y (x +x ))/((x -x )(x -x )))+ ((y (x + x ))/((x - x )(x -x ))))+
1 2 3 1 2 1 3 2 1 3 2 1 2 3 3 1 2 3 1 3 2
(((y x x )/ ((x -x )(x - x )))+ ((y x x )/((x - x )(x -x ))) +((y x x )/ ((x -x )(x - x ))))
1 2 3 1 2 1 3 2 1 3 2 1 2 3 3 1 2 3 1 3 2
This allows us to calculate and tabulate the coefficients a, b and c just
by replacing y by ELOW(I) and x by the correspondent stopping range
i i
ELOW(i). In the routine GRANGI a
table of length 3x(NEKBIN- 1) is created:
Q(JINTRP+3*(I-1)+1) = a(R ,R , R ,E ,E ,E ) = A
I I+1 I+2 I I+1 I+2 I
+2) = 0.5((b(R , R ,R ,E ,E ,E ))/ (a(R ,R , R ,E ,E ,E ))) = B
I I+1 I+2 I I+1 I+2 I I+1 I+2 I I+1 I+2 I
+3) = ((c(R , R ,R ,E ,E ,E ))/ (a(R ,R , R ,E ,E ,E ))) = C
I I+1 I+2 I I+1 I+2 I I+1 I+2 I I+1 I+2 I
where E = ELOW(I) and R = R(ELOW(I)). The calculation of the energy loss
I I
is now done in the following steps:
- Evaluate the stopping range for the threshold energy (STOPC). This
only once at the beginning of each new track.
- From the energy of the particle derive the stopping range by a
quadratic interpolation of the range table:
2
R =-B +((A )/ (|A |))sqrt(B - (C -((E )/ (A ))))
0 I I I I I 0 I
The value of I is chosen according to the following. Let i be such that
ELOW(i)
max(i- 1,1) if((E - ELOW(i))/(ELOW(i+1)-ELOW(i)))<0. 7
0
I = { .
min(i, NEKBIN-1) if((E - ELOW(i))/(ELOW(i+1)-ELOW(i)))> =0.7
0
The value of 0.7 is an empirical number which minimises the
discontinuities of the reconstructed dE/dx curve.
- Evaluate the stopping range for the particle after the step:
R' = R -step. If this is less than the stopping range of a particle with
0 0
threshold energy, the particle is terminated as a stopping particle below
the energy cut. Otherwise, the final energy is computed as:
E' =A (C + R' (2B + R' ))
0 I I 0 I 0
- the energy loss is computed as:
DeltaE =E - E'
0 0
This value is then corrected to take into account the Landau fluctuations
(see [PHYS332]).
This method results in a dE/dx curve which is a set of connected straight
lines and not a step function.