Webmasters
Walter R. Nelson
Stanford Linear Accelerator Center
Stanford University
Stanford, CA 94305, U.S.A.
Hideo Hirayama
National Laboratory for High Energy Physics (KEK)
Oho-machi, Tsukuba-gun, Ibaraki, Japan
David W. O. Rogers
National Research Council of Canada
Ottawa K1A 0R6, Canada
31 December 1985
This manual is kept up-to-date on the EGS account at SLAC and is distributed along with other files on the EGS4 distribution tape.
** A. J. Cook, "Mortran3 User's Guide", SLAC Computation Research Group Technical Memorandum CGTM-209 (1983) [see also APPENDIX 4: "EGS4 Users Guide to Mortran3"].
In addition, Mortran is a structured programming language. It is implemented by means of a set of macros which are used by a macro processor to translate the language into FORTRAN. The resultant FORTRAN is then run like any other FORTRAN code. Because the language is structured, the codes that are written tend to be more logical in appearance and much easier to read---especially to an outsider, or the author himself at a later time. Another feature is that one can switch back and forth between Mortran and FORTRAN with ease by means of the %M and %F control statements.
Although there might be some resistance by users of EGS to learn another language, we would like to point out two facts:
We would encourage EGS users not to do the latter, however, for this would truly defeat the real purpose for using Mortran---namely, the macro facility. Admittedly, Mortran macros can be difficult to read and understand, depending on their complexity. In most cases, however, the user should not have to understand the more complicated ones, and should be able to master the few that will be used in each problem.
As described in Chapter 2 of SLAC-265 ("The EGS4 Code System"), the EGS code itself consists of two User-Callable subroutines, HATCH and SHOWER, which in turn call the other subroutines in the EGS code, some of which call two User- Written subroutines, HOWFAR and AUSGAB. This is best illus- trated with the aid of Fig. A2.2.1.
To use EGS the user must write a "User Code." This consists of a MAIN program and the subroutines HOWFAR and AUSGAB, the latter two determining the geometry and output (scoring), respectively. Additional auxiliary subprograms might be included in the User Code to facilitate matters. The user can communicate with EGS by means of various COMMON variables. Usually MAIN will perform any initialization needed for the geometry routine, HOWFAR, and sets the values of certain EGS COMMON variables which specify such things as names of the media to be used, the desired cutoff energies, and the distance unit to be used (e.g., inches, centimeters, radiation lengths, etc.). MAIN then calls the HATCH subroutine which "hatches EGS" by doing necessary once-only initialization and by reading material data for the media from a data set that had been previously created by PEGS. This initialization completed, MAIN may then call SHOWER when desired. Each call to SHOWER results in a gener- ation of one history (often referred to as a "case"). The arguments to SHOWER specify the parameters of the incident particle initiating the cascade.
In addition, macro definitions can be included in MAIN in order to control or over-ride various functions in EGS as well as in the User-Written codes.
+---------+ +-------------+
| User | | Information |
| Control | | Extracted | U
| Data | | From Shower | S
+---------+ +-------------+ E
*** * R
*** ***
*** ***** C
***** *** O
*** *** D
* *** E
+---------+ +--------+ +--------+
| MAIN | | HOWFAR | | AUSGAB |<----- + +---------+ +--------+ +--------+ | | | A A A A | | | | | | | | | | | +-----------+ | |="====" |="=====" |="============" |="==========" |="|=" |="=======" |="=" | | | | | | | | | | +--------+ | | | V V | | | | | +-------+ +--------+ +--------+ +--------+ | | HATCH | | SHOWER |--->| ELECTR | +->| PHOTON |--> + |
+-------+ +--------+ +--------+ | +--------+ | |
* | | | | | E
*** +-------> | --------> + | | G
***** +--+ | | S
*** | +-------+ +-------+ | |
*** +-->| MSCAT | +--| COMPT |<-+ | C *** | +-------+ | +-------+ | | O *** | +--| PAIR |<-+ | D *** | +-------+ | +-------+ | | E *** +-->| ANNIH |--+ | | |
+---------+ | +-------+ | | +-------+ | |
| Media | +-->| BHABHA|--+ | | PHOTO |<-+ | | Data | | +-------+ | | +-------+ | | (PEGS) | +-->| MOLLER|--+ | | |
+---------+ | +-------+ | | | |
+-->| BREMS |--+ | +------> +
+---------+ | +-------+ | | |
| BLOCK | | V V |
| DATA | | +------+ |
|(Default)| +----------> | UPHI |-----------> +
+---------+ +------+
Fig. A2.2.1 Flow Control with User Using EGS
In summary, the user communicates with EGS by means of:
1.) Subroutines HATCH --- to establish media data
SHOWER --- to initiate the cascade
HOWFAR --- to specify the geometry
AUSGAB --- to score and output the
results
2.) COMMON blocks --- by changing values of variables
3.) Macro definitions --- re-definition of pre-defined
features.
To reiterate, we shall refer to the MAIN/HOWFAR/AUSGAB
combination (plus auxiliary subprograms and macros) as the
User Code. The following sections discuss these things in
greater detail.
Listed here are the names of the COMMON blocks relevant to the user (and relevant variables contained in them) with a brief description of their functions. Their useage will be discussed in more detail in subsequent sections. The easiest way to declare any of the COMMON blocks is with the COMIN macro. For example, COMIN/STACK,BOUNDS/ will automatically expand to the correct COMMON/STACK/ and COMMON/BOUNDS/ forms [Note: An asterisk (*) in the left column indicates a change from EGS3].
Table A2.3.1
-----------------------------------------------------------------
COMMON BLOCK VARIABLE FUNCTION
-----------------------------------------------------------------
BOUNDS ECUT Array of regions' charged particle
cutoff energies in MeV.
PCUT Array of regions' photon cutoff
energies in MeV.
VACDST Distance to transport in vacuum
(default=1.E8).
EPCONT EDEP Energy deposited in MeV (Double
Precision).
TSTEP Distance to next interaction (cm)
TUSTEP Total (curved) step length
requested.
USTEP User (straight line) step length
requested and granted.
TVSTEP Actual total (curved) step length
to be transported.
VSTEP Actual (straight line) step
length to be transported.
IDISC User discard request flag (to be set
in HOWFAR). IDISC .GT. 0 means user re-
quests immediate discard, IDISC .LT. 0
means user requests discard after com-
pletion of transport, and IDISC=0 (de-
fault) means no user discard requested.
IROLD Index of previous region.
IRNEW Index of new region.
* RHOF Value of density correction (default=1)
EOLD Charged particle (total) energy
at beginning of step in MeV.
ENEW Charged particle (total) energy
at end of step in MeV.
* BETA2 Beta squared for present particle.
(Note: BETA is no longer included).
IAUSFL Array of flags for turning on
various calls to AUSGAB.
EKE Kinetic energy of charged
particle in MeV.
EPCONT ELKE Natural logarithm of EKE.
GLE Natural logarithm of photon energy.
TSCAT See Eq. 2.14.82 in SLAC-265.
MEDIA NMED Number of media being used
(default=1).
MEDIA Array containing names of media
(default is NaI).
* IRAYLM Array of flags for turning on (=1)
coherent (Rayleigh) scattering in
various media. Set in HATCH based
on values of IRAYLR.
RLC Array containing radiation
lengths of the media in cm.
RLDU Array containing radiation
lengths of the media in distance
units established by DUNIT.
RHO Array containing density of the
media in g/cm**3.
MISC MED Array containing medium index for
each region.
DUNIT The distance unit to be used.
DUNIT=1 (default) establishes all
distances in cm; whereas,
DUNIT=2.54 establishes all
distances in inches.
KMPI FORTRAN unit number (default=12)
from which to read material data.
KMPO FORTRAN unit number (default=8) on
which to "echo" material data (e.g.,
printed output, "dummy" output, etc.).
RHOR Array containing the density for each
region (g/cm**3). If this is differ-
ent than the default density for the
medium for that region, the cross
sections and stopping powers (with
the exception of the density effect)
are scaled appropriately.
* NOSCAT Number of times multiple scattering
has been bypassed in subroutine MSCAT
(initialized to 0 in BLOCK DATA).
* IRAYLR Array of flags for turning on (=1)
coherent (Rayleigh) scattering in
various regions (default=0).
RANDOM IXX Random number generator seed
(default=123456789).
STACK Note: This COMMON block contains the
information about the particles
currently in the shower. All
of the following variables are
arrays except NP.
E Total energy in MeV (Double
Precision).
X,Y,Z Position of particle in units
established by DUNIT.
U,V,W Direction cosines of particle (not
necessarily normalized---see
Section A2.4.1c).
DNEAR A lower bound of distance from
(X,Y,Z) to nearest surface of
current region.
STACK WT Statistical weight of current particle
(default=1.0). To be used in conjunc-
tion with variance reduction tech-
niques as determined by user.
IQ Integer charge of particle
(+1,0,-1).
IR Index of particle's current region.
NP The stack pointer (i.e., the
particle currently being pointed
to). Also, the number of particles
on the stack.
THRESH RMT2 Twice the electron rest mass
energy in MeV.
RMSQ Electron rest mass energy squared
in MeV-squared.
AP Array containing PEGS lower photon
cutoff energy for each medium in
MeV.
UP Array containing PEGS upper photon
cutoff energy for each medium in
MeV.
AE Array containing PEGS lower charged
particle cutoff energy for each
medium in MeV.
UE Array containing PEGS upper charged
particle cutoff energy for each
medium in MeV.
TE Same as AE except kinetic energy
rather than total energy.
THMOLL Array containing the Moller thresh-
hold energy (THMOLL=AE+TE) for
each medium in MeV.
UPHIOT THETA Collision scattering angle (polar).
SINTHE Sine of THETA.
COSTHE Cosine of THETA.
SINPHI Sine of PHI (the azimuthal
scattering angle of the collision).
COSPHI Cosine of PHI.
PI Pi.
TWOPI Twice pi.
USEFUL MEDIUM Index of current medium. If
vacuum, then MEDIUM=0.
MEDOLD Index of previous medium.
RM Electron rest mass energy in MeV.
PRM "Precision" electron rest mass
energy in MeV (Double Precision).
PRMT2 Twice PRM (Double Precision).
* IBLOBE Flag indicating if photon is below
binding energy (EBINDA) after a
photoelectric interaction (yes=1).
The sequence of operations needed for the correct operation of EGS is shown below along with sections in this chapter for additional details.
Step 1. User-Over-Ride-Of-EGS-Macros (A2.4.1).
Step 2. Pre-HATCH-Call-Initialization (A2.4.2).
Step 3. HATCH-Call (A2.4.3).
Step 4. Initialization-For-HOWFAR (A2.4.4).
Step 5. Initialization-For-AUSGAB (A2.4.5).
Step 6. Determination-Of-Incident-Particle-Parameters
(A2.4.6).
Step 7. SHOWER-Call (A2.4.7).
Step 8. Output-Of-Results (A2.4.8).
The following are restrictions on the order of these operations:
Details for the above steps are given in the following sub-sections.
EGS macros which the user might want to over-ride include the following:
$MXMED Maximum number of media (default=10).
$MXREG Maximum number of regions (default=2000).
For example, to extend the number of media to 25, include the statement
REPLACE {$MXMED} WITH {25}
in the User Code.
$RANDOMSET RNUMBR;
For efficiency purposes, EGS currently uses an "in-line" simulation of the SLAC system routine called RAN6. If the user desires to use some other random number generator (e.g., RAN10), this can be most easily accomplished by including the override macro
REPLACE {$RANDOMSET#;} WITH {{P1}=RAN10(0);}
in the User Code at Step 1 (prior to any corres- ponding pattern). Note, however, that not using an in-line random number generator slows EGS down significantly. Of course, the user must make sure that the random number generator is properly initialized. As an example, suppose we wish to "seed" the random number generator in the User Code by the statement
IXX=987654321;
which could be included at any step prior to Step 7. To properly initialize, however, the associated COMMON block, RANDOM, would have to be included in the declaration section of the User Code. This is best done with the statement
COMIN/RANDOM/;If no initialization is provided, EGS automatic- ally seeds with IXX=123456789 in BLOCK DATA (default). SLAC users should note that RAN6 must be re-seeded (with an odd integer number) for any new run of the same problem in order to avoid getting identical statistical results.
REPLACE {$EVALUATE#USING SIN(#);}
WITH {{P1}=SIN({P2});}
REPLACE {$SET INTERVAL#,SINC;} WITH {;}
which are thereby included at Step 1 of the User Code.
[The reader is referred to the Mortran3 User's Guide
as an aid in understanding the macros].
It should be pointed out that due to the precision involved in the table look-up method, the direction cosines can become slightly unnormalized.
Depending on the problem at hand, this can lead to incorrect results---such as when two direction cosines are simultaneously involved in an angular sort of particles. The problem can generally be remedied by renormalizing the direction cosines prior to using them.
REPLACE {$CHARGED-TRANSPORT;} WITH {CALL MYTRAN;}
could be included in Step 1 of the User Code, and an appropriate subroutine MYTRAN would need to be provided by the user.
Another pattern that has been included in ELECTR is $TMXS-OVER-RIDE, which is there in case the user wants to impose an additional constraint on the electron transport step size. It can be used, for example, to revert back to the method used in EGS1 in which step sizes were taken discretely (rather than randomly as is done in subsequent EGS versions). To do this one uses the macro
REPLACE {$TMXS-OVER-RIDE;}
WITH {;IF(E(NP).GT.5.0) [TMXS=0.01*RLDU(MEDIUM);]
ELSE [TMXS=0.001*RLDU(MEDIUM);]}
Note: Always prefix an additional semicolon to IF-
statements that begin replacement strings.
The macro
REPLACE {$TMXS-OVER-RIDE;}
WITH {TP=1.335E-3*EKE**2*RLDU(MEDIUM);
TMXS=AMIN1(TMXS,TP);}
can be used to limit path length corrections (see Chapter 2 of SLAC-265) to values less than 15 per cent (whether or not this is important depends on the problem being solved). On the other hand, to remove the path length correction entirely (as in EGS1 and EGS2), one can use the macro-pair
REPLACE {$SET-USTEP;} WITH {USTEP=TUSTEP;}
REPLACE {$SET-TVSTEP;} WITH {TVSTEP=VSTEP;}
Finally, the macro
REPLACE {$TMXS-OVER-RIDE;}
WITH {TP=200.0*TEFF0(MEDIUM); TMXS=AMIN1(TMXS,TP);}
limits the step size essentially in the manner done in EGS2, and we use this as the default macro for EGS3, and now for EGS4, for reasons discussed and demon- strated at the end of Section 3.5 of SLAC-210. For low energy problems, other step size algorithms should be used (see end of Section 2.14 of SLAC-265).
This step consists of setting EGS COMMON variables that are used by HATCH in its initialization operations. All of these variables are initialized to some reasonable value in the BLOCK DATA subprogram. Therefore, if different values are desired they should be set with executable code (as opposed to another BLOCK DATA). Concurrently, the various COMMON blocks (i.e., BOUNDS, MEDIA, MISC) will have to be included in the declaration section of the MAIN program of the User Code. These variables are:
NMED --- This must be initialized to the number of media to be used in the shower generation (default=1).
MEDIA --- This array contains the names of the media required and is dimensioned MEDIA(24,$MXMED), where $MXMED is an EGS macro that is currently defined to be 10 (default), and whose value is the maximum number of media for which array space has been allocated (see Section A2.4.1a above). The media names are stored in MEDIA in alphameric field specification A1 to ensure transportability. Each medium name is 24 characters long. For the convenience of users compiling with EGS' macros, there is a macro to generate A1 strings. For example,
$S'STRING' expands to 'S','T','R','I','N','G'.
One way of implementing this in the User Code is
demonstrated in the next example, which is for
three media: lead, steel, and air at NTP.
A temporary array is declared and initialized in MAIN by
$TYPE TEMP(24,3)/$S'PB', 22*' ',
$S'STEEL', 19*' ',
$S'AIR AT NTP',14*' '/;
Note: $TYPE is recognized as a macro template and gets expanded to CHARACTER*4 for FORTRAN 77 and to INTEGER for FORTRAN IV.
Then at Step 2 one puts
NMED=3; "NUMBER OF MEDIA USED"
DO J=1,NMED [DO I=1,24 [MEDIA(I,J)=TEMP(I,J);]]
MED --- This array, which is dimensioned MED($MXREG), contains the medium indices for each region (default values are 1 for all $MXREG). A medium index of zero means a region is filled with a vacuum. For instance, if we consider the three media example above along with vacuum to define four regions, we might have
MED(1)=3; "FIRST REGION IS AIR AT NTP"
MED(2)=1; "SECOND REGION IS LEAD"
MED(3)=0; "THIRD REGION IS VACUUM"
MED(4)=2; "FOURTH REGION IS STEEL"
in Step 2 of the User Code.
ECUT and PCUT --- These arrays contain the cutoff energies (in MeV) for charged particles and photons, respectively, for each region. They are dimen- sioned ECUT($MXREG) and PCUT($MXREG) and are given temporary (default) values of 0.0 in BLOCK DATA. At the time that data for each medium are gen- erated in the preprocessing code (PEGS), two para- meters (AE and AP) are set to the lowest energies at which it will be desired to transport electrons and photons. When the EGS subroutine HATCH is called, these AE and AP values are read in and HATCH upgrades the values of ECUT and PCUT such that the maximum of the current (ECUT,AE) (and (PCUT,AP)) is chosen. Therefore, by assigning values of ECUT and PCUT prior to the HATCH call, the user can raise (but not lower) the cutoff energies in this manner. For instance, consider the four region example from above. The statement
DO I=1,3 [ECUT(I)=10.0; PCUT(I)=100.0;]
when put in Step 2 of the User Code results in charged particle histories being terminated at 10.0 MeV (total energy) and photon histories being terminated at 100.0 MeV in the first three regions only. In the fourth region the respective cutoffs are set by AE and AP as established by PEGS. Of course COMMON/BOUNDS/ will have to be declared in the routine that calls HATCH in order to pass ECUT and PCUT to HATCH. Combined with COMMON/MEDIA/ and COMMON/MISC/, the macro declara- tion might look like
COMIN/BOUNDS,MEDIA,MISC/;
IRAYLR --- The elements of this array (dimensioned
IRAYLR($MXREG) and passed in COMMON/MISC/), are to
be set to 1 prior to calling HATCH if coherent
(Rayleigh) scattering is to be done in a particular
region. Execution is terminated if Rayleigh data
is not included in the PEGS data, however.
DUNIT --- This parameter determines the unit of distance to be used in the shower simulation (the default is cm if DUNIT=1.0). On input to HATCH, this parameter will be interpreted as follows:
The distance unit used by PEGS is the radiation length. After HATCH interprets DUNIT, it scales all distance- type data from PEGS in the proper way, so that all subsequent operations in EGS will be correctly performed with all distances in units of DUNIT (default value: 1.0 cm).
This step is very simple---HATCH has no arguments, so all one has to do is:
CALL HATCH;
The following is a typical output message when DUNIT has not been changed (and Rayleigh data is included in the file):
RAYLEIGH DATA AVAILABLE FOR MEDIUM 1 BUT OPTION NOT REQUESTED. EGS SUCCESSFULLY 'HATCHED' FOR ONE MEDIUM.
However, if the user has set DUNIT=2.54 prior to calling HATCH, the message will look like the following (two media, no Rayleigh data):
DUNIT REQUESTED&USED ARE: 2.54000E+00 2.54000E+00(CM.) EGS SUCCESSFULLY 'HATCHED' FOR 2 MEDIA.
Failure to "hatch", on the other hand, will result in the message:
END OF FILE ON UNIT 12
PROGRAM STOPPED IN HATCH BECAUSE THE
FOLLOWING NAMES WERE NOT RECOGNIZED:
(list of names)
followed by a STOP in HATCH. [Note: one cannot ask for the same medium twice].
As stated previously, HOWFAR is the routine that deter- mines the geometry of the regions. Although initialization for items that are taken up in HOWFAR can be done at any step prior to calling SHOWER (Step 7), Step 4 allows a space in MAIN to consider if such initialization need be performed. For example, if regions are defined by semi-infinite planes, data defining each plane (e.g., coordinates and unit normal vectors) can be established here. The data may be referred to in HOWFAR or by user-written subprograms called by HOWFAR. It may be that some of the dimensions of the regions are determined at run-time, or the geometry may be so complex that it is desirable to use executable code to generate tables for use by HOWFAR. In such cases, initialization for HOWFAR will probably consist of filling up some user-written COMMON blocks for HOWFAR.
COMMON/TOTALS/ESUM($MXREG);in both the MAIN code and in AUSGAB, and we could add
DO I=1,$MXREG [ESUM(I)=0.0;]to the MAIN code (at Step 5). Then the statement
ESUM(IR(NP))=ESUM(IR(NP)) + EDEP;in AUSGAB could keep a running total of the energy deposited in each region under consideration (Note: the EGS utility subroutines, ECNSV1 and NTALLY, have been developed to facilitate keeping track of where the energy is deposited and how many events are scored, respectively).
This would be a good time to point out that EDEP is a double precision variable---established as such via a macro. Therefore, one might wish to establish ESUM as double pre- cision as well (in both MAIN and AUSGAB). We have experienced situations whereby energy balancing could not be attained due to round-off error difficulties. This was particularly evident for large shower-history problems involving the addition of small energy values to large numbers in various regions. As a result of this experience, certain key energy variables in the EGS code have been defined as double precision. Users may take advantage of this at their discretion.
This step is really self-explanatory---particularly when looked at in conjunction with Step 7 below. A specfic exam- ple of such coding might be useful and is given as follows:
IQI=-1; "INCIDENT PARTICLE IS AN ELECTRON" EI=1000.0; "TOTAL ENERGY (MEV)" XI=0.0; YI=0.0; ZI=0.0; "PARTICLE COORDINATES" UI=0.0; VI=0.0; WI=1.0; "DIRECTION COSINES" IRI=2; "REGION NUMBER 2 IS THE INCIDENT REGION" WTI=1.0; "WEIGHT FACTOR IN IMPORTANCE SAMPLING" IXX=987654321; "RANDOM NUMBER GENERATOR SEED" NCASES=10; "NUMBER OF HISTORIES TO RUN"
The calling sequence for SHOWER is:
CALL SHOWER(IQI,EI,XI,YI,ZI,UI,VI,WI,IRI,WTI);The types of the arguments are given by their starting letter in accordance with standard FORTRAN convention. These arguments specify the charge, total energy, position, direction, region index, and statistical weight of the incident particle, and are used to fill the corresponding stack variables (see COMMON/STACK/ in Section A2.3). Section A2.4.6 above might be of some aid in understanding the parameter list. The subroutine may be called repeatedly by means of statements like
DO I=1,NCASES [CALL SHOWER(IQI,EI,XI,....,etc.);] .
The statistical weight, WTI, is generally taken as unity unless variance reduction techniques are employed by the user. It should noted, however, that if IQI is assigned the value of 2, subroutine SHOWER recognizes this as a pi-zero meson decay event, and two photons are added to the stack with energies and direction cosines appropriately obtained by sampling.
This step has been added for completeness and is self-explanatory.
USTEP, IDISC, IRNEW, and DNEAR(NP).
Except for the last variable (which is in COMMON/STACK/), these are available to the user via COMMON/EPCONT/. The ways in which these may be changed, and the way EGS will interpret these changes, will now be discussed in detail. [Note, flow diagrams for subroutines ELECTR and PHOTON have been included in Appendix 1 of SLAC-265 for the user who requires a more complete understanding of what actu- ally takes place during particle transport].
If the user decides that the current particle should be discarded, then IDISC must be set nonzero (the usual convention is to set IDISC=1).
A positive value for IDISC will cause the particle to be discarded immediately. A negative value for IDISC will cause EGS to discard the particle when it completes the transport. EGS initializes IDISC to zero, and if left zero no user requested discard will take place. For example, the easiest way to define an infinite, homo- geneous medium is with the HOWFAR routine:
SUBROUTINE HOWFAR;
RETURN;
END;
In this case, particle transport will continue to take place until energy cutoffs are reached. However, a common procedure is to set IDISC=1 whenever the particle reaches a discard region.
If immediate discard has not been requested, then the user should check to see whether transport by distance USTEP will cause a region boundary to be crossed. The presence of the region index for the current particle, IR(NP), should make this task much easier than if only the posi- tion of the particle were known. If no boundary will be crossed, then USTEP and IRNEW may be left as they are. If a boundary will be crossed, then USTEP should be set to the distance to the boundary from the current position along the current direction, and IRNEW should be set to the region index of the region on the other side of the boundary. For sophisticated geometries, this is the most complex part of the User Code.
The setting of DNEAR(NP) by the user is optional. However, in many situations a significant gain in efficiency will result by defining DNEAR in HOWFAR. First of all, it is clear that boundary checking takes time and should be avoided whenever possible. If EGS had no way of knowing how far it was from a boundary, then it would have to ask the user how far to go every time it wanted to transport a particle. This would not be too serious for photons since they usually go relatively far on each transport. However, the transport of a charged particle from one interaction to the next requires the path-length to be split up into smaller lengths in order to simulate properly the multiple scattering process. If, relative to the small steps being taken, the particle is a fairly good distance from the nearest boundary, then checking for boundary crossings on each transport is a waste of time. In order to avoid this inefficiency, each particle has stored on the stack a variable called DNEAR, which is used by EGS to store a lower bound to the distance from the particle's current position to the nearest region boundary. This variable is used by EGS in the following ways:
| |
| |
Region | Region | Region
| |
1 | 2 | 3
| |
| |
| (X,Y,Z) |
| x |
| . |
Vacuum | . | Air at NTP
| . |
O---------.+----------------x Z
| |
| |.
| Iron | .
| | .
| | .
| | .
| | .
| | .
| | x
| | (U,V,W)
| |
| |
| |
--x| ZTHICK |x--
| |
|
|
V
X
(Y into paper)
Fig. A2.5.1 A Three Region Geometry Example for HOWFAR
Consider, as an example of how to write a HOWFAR sub- program, the three region geometry in Fig. A2.5.1. A particle is shown in Region 2 with coordinates (X,Y,Z) and direction cosines (U,V,W). We will assume that the slab of thickness ZTHICK is semi-infinite (x and y-directions), and that particles are immediately discarded whenever they go into Region 1 or Region 3. The following HOWFAR code is then applicable:
SUBROUTINE HOWFAR; COMIN/EPCONT,STACK/; "COMMON BLOCKS NEEDED IN CALCULATIONS" COMMON/PASSIT/ZTHICK; "SLAB THICKNESS DEFINED IN MAIN" IF(IR(NP).NE.2) [IDISC=1; RETURN;] "MIGHT AS WELL SET DNEAR NEXT" DNEAR(NP)=AMIN1(Z(NP),ZTHICK-Z(NP)); IF(W(NP).EQ.0.0) [RETURN; "PARTICLE GOING PARALLEL TO PLANES"] "CHECK FORWARD PLANE FIRST SINCE SHOWER HEADING THAT WAY" "MOST OF THE TIME" IF(W(NP).GT.0.0) [DELTAZ=(ZTHICK-Z(NP))/W(NP); IRNEXT=3;] "OTHERWISE, PARTICLE MUST BE HEADING IN BACKWARDS DIRECTION" ELSE [DELTAZ=-Z(NP)/W(NP); IRNEXT=1;] "NOW CHECK WITH USTEP AND RESET THINGS IF NECESSARY" IF(DELTAZ.LE.USTEP) [USTEP=DELTAZ; IRNEW=IRNEXT;] RETURN; END;
Note: A number of geometry subprograms and their macro equivalents are distributed with the EGS Code System in order to make it easier to write HOWFAR. For example, subroutine PLAN2P, or its equivalent macro $PLAN2P, could have been used in place of several lines above and the program would have been easier to read.
The subroutine AUSGAB is called by EGS with the statement:
CALL AUSGAB(IARG);
The argument IARG indicates the situation under which AUSGAB is being called. IARG can take on 25 values starting from zero (i.e., IARG=0 through IARG=24), although only the first five are called in the default version of EGS. The remaining 20 IARG values must be "switched-on" by means of the array IAUSFL. The value for IARG and the corresponding situations are given in Table A2.6.1.
Table A2.6.1
-----------------------------------------------------------------
IARG Situation
-----------------------------------------------------------------
0 Particle is going to be transported by distance TVSTEP.
1 Particle is going to be discarded because its energy
is below the cutoff ECUT (for charged particles)
or PCUT (for photons)---but its energy is larger
than the corresponding PEGS cutoff AE or AP,
respectively.
2 Particle is going to be discarded because its energy
is below both ECUT and AE (or PCUT and AP).
3 Particle is going to be discarded because the user
requested it (in HOWFAR usually).
4 A photoelectric interaction has occurred and either:
a) the energy of the incident photon was below
the K-edge binding energy and it is going to
be discarded, or
b) a (fluorescent) photon is going to be dis-
carded with the K-edge binding energy.
The above IARG values are the ones generally required in the majority of situations in which EGS is used to simulate electromagnetic cascade shower development. In particular, IARG=0 is useful whenever track lengths are being calculated or when charged particle ionization loss is needed. Also, as a check on energy conservation, EDEP can be summed in AUSGAB for all IARG values less than 5. The extended IARG range allows the user to extract additional information without making changes to the EGS coding. To do this we have created the integer flag array, IAUSFL(J), for J=1 through 25. It takes on values of 1 or 0 depending on whether AUSGAB is called or not, respectively. For J=1 through 5, which corresponds to IARG=0 through 4, IAUSFL(J)=1 (default). In other words, AUSGAB is always called for the situations listed in Table A2.6.1. For the remaining values of J, corresponding to IARG=5 through 24, IAUSFL(J)=0 (default). The value for IARG and the correspond- ing situations for this upper set of IARG values are shown in Table A2.6.2.
Table A2.6.2
-----------------------------------------------------------------
IARG IAUSFL Situation
-----------------------------------------------------------------
5 6 Particle has been transported by distance TVSTEP.
6 7 A bremsstrahlung interaction is to occur and a call
to BREMS is about to be made in ELECTR.
7 8 Returned to ELECTR after a call to BREMS was made.
8 9 A Moller interaction is to occur and a call to
MOLLER is about to be made in ELECTR.
9 10 Returned to ELECTR after a call to MOLLER was made.
10 11 A Bhabha interaction is to occur and a call to
BHABHA is about to be made in ELECTR.
11 12 Returned to ELECTR after a call to BHABHA was made.
12 13 An in-flight annihilation of the positron is to
occur and a call to ANNIH is about to be made
in ELECTR.
13 14 Returned to ELECTR after a call to AHHIH was made.
14 15 A positron has annihilated at rest.
15 16 A pair production interaction is to occur and a
call to PAIR is about to be made in PHOTON.
16 17 Returned to PHOTON after a call to PAIR was made.
17 18 A Compton interaction is to occur and a call to
COMPT is about to be made in PHOTON.
18 19 Returned to PHOTON after a call to COMPT was made.
19 20 A photoelectric interaction is to occur and a call
to PHOTO is about to be made in PHOTON.
20 21 Returned to PHOTON after a call to PHOTO was made
(assuming NP is non-zero).
21 22 Subroutine UPHI was just entered.
22 23 Subroutine UPHI was just exited.
23 24 A coherent (Rayleigh) interaction is about to occur.
24 25 A coherent (Rayleigh) interaction has just occurred.
As an example of how to write an AUSGAB subprogram, consider the previous three region geometry (Fig. A2.5.1). Suppose that we wish to score (i.e., output on the line printer) only photons that emanate from Region 2 into Region 3. The AUSGAB subprogram that will accomplish this is given below. In this example we print out the stack variables plus IARG.
SUBROUTINE AUSGAB(IARG);
COMIN/STACK/;
"ONLY OUTPUT INFORMATION FOR PHOTONS THAT ARE DISCARDED"
"(BY THE USER) IN REGION 3"
IF(IARG.EQ.3.AND.IQ(NP).EQ.0.AND.IR(NP).EQ.3) [
OUTPUT E(NP),X(NP),Y(NP),Z(NP),U(NP),V(NP),W(NP),
IQ(NP),IR(NP),IARG; (7G15.7,3I5);]
RETURN; END;
The following User Code, called UCSAMPL4, simulates electro- magnetic cascade showers initiated by 1 GeV electrons that are incident (normally) on a 3 cm, semi-infinite slab of iron. The upstream region of the slab is vacuum and the downstream region is air at NTP. A particle is discarded whenever it leaves the slab (on either side), or whenever its total energy falls below a preset cutoff energy of 100 MeV [Note: the medium assigned to Region 3 is really not important in this example, and was included solely for purposes of illustration]. Some of the stack variable information E(NP), Z(NP), W(NP), IQ(NP), IR(NP), plus the IARG value, is output on the line printer (first 15 lines only) for photons reaching Region 3.
The maximum number of regions is changed from 2000 (default) to 20 by means of an over-ride macro at Step 1 in the User Code. Furthermore, the distance unit (DUNIT) is changed (at Step 2) so that distances are in radiation lengths of iron rather than in centimeters. A total of 10 cases of incident electrons is run and the total energy fraction for each region is summed and printed out at the end of the run for an energy balance check. Finally, the last random number generator (integer) that was used is printed out for possible use in avoiding statistically identical results in future runs of the same problem.
The UCSAMPL4 User Code (FORTRAN 77 compatible) is given below.
"******************************************************************"
"*************************** STANFORD LINEAR ACCELERATOR CENTER"
"*** U C S A M P L 4 *** "
"*************************** EGS4 USER CODE -- 23 NOV 1985/1830"
"******************************************************************"
" "
"************"
"*** MAIN ***"
"************"
"STEP 1. USER-OVER-RIDE-OF-EGS-MACROS"
" THE FOLOWING RANDOM NUMBER GENERATOR CAN BE USED ON A VAX. IT IS"
" COMMENTED OUT BELOW BECAUSE THE DEFAULT ONE PROVIDED WITH THE EGS4"
" MACROS (IBM COMPATIBLE) WAS THE ONE ACTUALLY USED IN THIS EXAMPLE."
"REPLACE {;COMIN/RANDOM/;} WITH {;COMMON/RANDOM/IXX;}"
"REPLACE {$RANDOMSET#;} WITH"
"{IXX=IXX*663608941;{P1}=0.5 + IXX*0.23283064E-09;}"
REPLACE{$MXREG} WITH {20} "OVER-RIDE MAXIMUM NO. OF REGIONS"
;COMIN/BOUNDS,MEDIA,MISC,USEFUL/; "COMMONS NEEDED"
COMMON/PASSIT/ZTHICK; "SLAB THICKNESS....NEEDED IN HOWFAR"
COMMON/LINES/NLINES,NWRITE; "TO KEEP TRACK OF LINES-PRINTED"
COMMON/TOTALS/ESUM($MXREG); "FOR ENERGY CONSERVATION CHECK"
$ENERGY PRECISION EI,ESUM,EKIN,TOTKE,ETOT; "DOUBLE PRECISION"
"CREATE A TEMPORARY ARRAY AND DEFINE THE MEDIA, NEXT"
$TYPE TEMP(24,2)/$S'FE',22*' ',$S'AIR AT NTP',14*' '/;
COMIN/RANDOM/; "LOCATED HERE TO AVOID FORTRAN 77 DIAGNOSTIC"
"STEP 2. PRE-HATCH-CALL-INITIALIZATION"
NREG=3; "THE NUMBER OF REGIONS---A LOCAL VARIABLE ONLY"
NMED=2; "TWO MEDIA WILL BE USED"
DO J=1,NMED [DO I=1,24 [MEDIA(I,J)=TEMP(I,J);]]
MED(1)=0; "REGION 1 IS VACUUM"
MED(2)=1; "REGION 2 IS IRON"
MED(3)=2; "REGION 3 IS AIR AT NTP"
"SET ENERGY CUTOFFS FOR EACH REGION NEXT"
DO I=1,NREG [ECUT(I)=100.0; PCUT(I)=100.0;]
"STEP 3. HATCH-CALL"
CALL HATCH;
"STEP 4. INITIALIZATION-FOR-HOWFAR"
ZTHICK=3.0; "SLAB THICKNESS IN CENTIMETERS"
"STEP 5. INITIALIZATION-FOR-AUSGAB"
DO I=1,$MXREG [ESUM(I)=0.D0;] "ZERO THE ENERGY BALANCE ARRAY"
NLINES=0; "INITIALIZE THE NLINES-COUNTER"
NWRITE=15; "THE NUMBER OF LINES TO PRINT OUT"
"STEP 6. DETERMINATION-OF-INCIDENT-PARTICLE-PROPERTIES"
IQI=-1; "INCIDENT PARTICLE IS AN ELECTRON"
EI=1000.D0; "INCIDENT ENERGY (TOTAL) IN MEV"
EKIN=EI-PRM; "K.E. OF ELECTRON---PRM IS THE REST MASS"
XI=0.0; YI=0.0; ZI=0.0; "COORDINATES OF INCIDENT PARTICLE"
UI=0.0; VI=0.0; WI=1.0; "DIRECTION COSINES---ALONG Z=AXIS"
IRI=2; "INCIDENT PARTICLE STARTS OUT IN REGION 2---IRON"
WTI=1.0; "WEIGHT FACTOR---NOT USED IN CALCULATION, BUT"
" IS A PARAMETER IN SUBROUTINE SHOWER; HENCE DEFINE"
" AS UNITY"
IXX=987654321; "RANDOM NUMBER GENERATOR SEED"
NCASES=10; "NUMBER OF HISTORIES (CASES) TO RUN"
ICODE=-1; "AN OUTPUTING PARAMETER, INVENTED TO MARK THE"
" INCIDENT PARTICLES"
"STEP 7. SHOWER-CALL"
OUTPUT; (/,' SHOWER RESULTS:',///,7X,'E',14X,
'Z',14X,'W',10X,'IQ',3X,'IR',2X,'IARG',/);
DO I=1,NCASES [
IF (NLINES.LT.NWRITE) [
OUTPUT EI,ZI,WI,IQI,IRI,ICODE;
(3G15.7,3I5);
NLINES=NLINES+1;]
CALL SHOWER(IQI,EI,XI,YI,ZI,UI,VI,WI,IRI,WTI);
"END OF SHOWER-CALL LOOP"]
"STEP 8. OUTPUT-OF-RESULTS"
TOTKE=NCASES*EKIN; "TOTAL K.E. INVOLVED IN RUN"
OUTPUT EI,ZTHICK,NCASES,IXX;
(//,' INCIDENT TOTAL ENERGY OF ELECTRON=',F12.1,' MEV',/,
' IRON SLAB THICKNESS=',F6.3,' CM',/,
' NUMBER OF CASES IN RUN=',I3,/,' LAST RANDOM NUMBER=',
I12,//,' ENERGY DEPOSITION SUMMARY:',//);
"CALCULATE AND PRINT OUT THE FRACTION OF ENERGY"
"DEPOSITED IN EACH REGION"
ETOT=0.D0;
DO I=1,NREG [
ETOT=ETOT+ESUM(I);
ESUM(I)=ESUM(I)/TOTKE; "FRACTION IN EACH REGION"
OUTPUT I, ESUM(I); (' FRACTION IN REGION',I3,'=',F10.7);
]
ETOT=ETOT/TOTKE; "THE TOTAL FRACTION OF ENERGY IN RUN"
OUTPUT ETOT; (//,' TOTAL ENERGY FRACTION IN RUN=',G15.7,/,
' WHICH SHOULD BE CLOSE TO UNITY');
STOP;
END; "LAST STATEMENT OF MAIN"
SUBROUTINE AUSGAB(IARG);
COMIN/EPCONT,STACK/; "COMMONS NEEDED IN AUSGAB"
COMMON/LINES/NLINES,NWRITE; "TO KEEP TRACK OF LINES-PRINTED"
COMMON/TOTALS/ESUM($MXREG); "FOR ENERGY CONSERVATION CHECK"
$ENERGY PRECISION ESUM; "DOUBLE PRECISION"
"KEEP A RUNNING SUM OF THE ENERGY DEPOSITED IN EACH REGION"
ESUM(IR(NP))=ESUM(IR(NP)) + EDEP;
"PRINT OUT THE FIRST NLINES OF STACK INFORMATION, ETC."
"BUT, ONLY FOR PHOTONS THAT ARE DISCARDED IN REGION 3"
IF (NLINES.LT.NWRITE) [
IF (IARG.EQ.3.AND.IQ(NP).EQ.0.AND.IR(NP).EQ.3) [
OUTPUT E(NP),Z(NP),W(NP),
IQ(NP),IR(NP),IARG; (3G15.7,3I5);
NLINES=NLINES+1; ]]
RETURN;
END; "LAST STATEMENT OF SUBROUTINE AUSGAB"
SUBROUTINE HOWFAR;
COMIN/EPCONT,STACK/; "COMMON NEEDED IN HOWFAR"
COMMON/PASSIT/ZTHICK; "SLAB THICKNESS DEFINED IN MAIN"
IF (IR(NP).NE.2) [IDISC=1; RETURN;]
"MIGHT AS WELL SET DNEAR NEXT"
DNEAR(NP)=AMIN1(Z(NP),ZTHICK-Z(NP));
IF (W(NP).EQ.0.0) [RETURN; "PARTICLE GOING PARALLEL TO PLANES"]
"CHECK FORWARD PLANE FIRST SINCE SHOWER HEADING THAT WAY"
"MOST OF THE TIME"
IF (W(NP).GT.0.0) [DELTAZ=(ZTHICK-Z(NP))/W(NP); IRNXT=3;]
"OTHERWISE, PARTICLE MUST BE HEADING IN BACKWARDS DIRECTION"
ELSE [DELTAZ=-Z(NP)/W(NP); IRNXT=1;]
"NOW CHECK WITH USTEP AND RESET THINGS IF NECESSARY"
IF (DELTAZ.LE.USTEP) [USTEP=DELTAZ; IRNEW=IRNXT;]
RETURN;
"**********************************************************"
"******************** END OF UCSAMPL4 *********************"
"**********************************************************"
END; "LAST STATEMENT OF EGS4 USER CODE UCSAMPL4"
The results of running this User Code on the IBM-3081 at SLAC are given below:
EGS SUCCESSFULLY 'HATCHED' FOR 2 MEDIA.
SHOWER RESULTS:
E Z W IQ IR IARG
1000.000 0.0000000E+00 1.000000 -1 2 -1
163.9946 2.999998 0.9999203 0 3 3
236.1791 2.999999 0.9993725 0 3 3
147.1071 2.999999 0.9993142 0 3 3
123.6069 2.999999 0.9997734 0 3 3
194.9122 2.999999 0.9987941 0 3 3
1000.000 0.0000000E+00 1.000000 -1 2 -1
107.8515 2.999999 0.9981139 0 3 3
113.2713 2.999999 0.9998847 0 3 3
1000.000 0.0000000E+00 1.000000 -1 2 -1
133.4546 2.999999 0.9977404 0 3 3
1000.000 0.0000000E+00 1.000000 -1 2 -1
339.5281 2.999999 0.9993687 0 3 3
448.6772 2.999999 0.9997988 0 3 3
1000.000 0.0000000E+00 1.000000 -1 2 -1
INCIDENT TOTAL ENERGY OF ELECTRON= 1000.0 MEV
IRON SLAB THICKNESS= 3.000 CM
NUMBER OF CASES IN RUN= 10
LAST RANDOM NUMBER= 1903435093
ENERGY DEPOSITION SUMMARY:
FRACTION IN REGION 1= 0.0000000
FRACTION IN REGION 2= 0.3319190
FRACTION IN REGION 3= 0.6680810
TOTAL ENERGY FRACTION IN RUN= 1.000000
WHICH SHOULD BE CLOSE TO UNITY
last updated 10/04/01