- Introduction
- Common inputs
- Scoring options
- Variance reduction techniques
- Monte Carlo transport parameters
- Usage
- An input example: NRC MEES FAC corrections.
Introduction
The C++ application egs_fac is an advanced EGSnrc application designed for the purpose of calculating free air chamber (FAC) correction factors. A self-consistent approach for the calculation of such correction factors, recently defined in [Mainegra-Hing, Reynaert and Kawrakow, Med. Phys. 35 (8), 2008, 3650-3660] is implemented which allows the conversion of the chamber reading into the quantity air-kerma at the point of measurement (POM). Correction factors for attenuation 
, scatter 
, electron loss 
, aperture leakage 
, beam geometry 
, lack of CPE 
and for the inconsistency between measured and MC calculated attenuation correction factors 
are directly calculated in egs_fac. Originally implemented in cavity, FAC correction calculations were ported to a stand-alone code due to the increasing complexity of the logic in cavity. As a result, a cleaner, more efficient implementation is obtained. For details on the definition, implementation and calculation of these correction factors the reader is referred to the above mentioned paper.
Scoring options
FAC correction factors can be calculated for several geometries if desired. The input key calculation geometry is used to define as many geometries as needed as shown below.
:start scoring options:
:start calculation geometry:
geometry name = sim_no_tube
cavity regions = list_of_cavity_region_indices
cavity mass = cavity mass in g of a cyl. defined by diaphragm and plates
aperture regions = regions_defining_aperture
front and back regions = f_reg b_reg
POM = 0.45 0.5 # first input is z-position, second is radius.
include scatter = no # yes,no(default)
:stop calculation geometry:
:start calculation geometry:
geometry name = sim_air_tube
...
:stop calculation geometry:
:start calculation geometry:
geometry name = sim_vacuum_tube
...
:stop calculation geometry:
:start calculation geometry:
geometry name = sim_no_fac
...
include scatter = yes # yes,no(default)
:stop calculation geometry:
...
:start calculation geometry:
geometry name = sim_n
...
:stop calculation geometry:
correlated geometries = sim_i sim_j
correlated geometries = sim_no_tube sim_no_fac
Ax calculation = sim_no_tube sim_air_tube sim_vacuum_tube
muen file = E*muen file name
:stop scoring options:
The definition of calculation geometries is very similar to the cavity code, but there are inputs added that define the aperture regions (needed to get the correction), the front/back CV regions (needed for the electron loss calculation) and the point of measurement (POM), needed for KERMA scoring. In the input example above, sim_air_tube,sim_vacuum_tube are geometries to model the evacuated-tube technique for experimentally determining at the NRCC as the ratio of the FAC readings with and without an evacuated tube. To correct for differences between the experimental and the MC calculated attenuation correction , a correction factor, 
, was introduced which can be calculated by using the Ax calculation input key. This input key expects the name of three geometries, the realistic FAC geometry (sim_no_tube), plus the two geometries mentioned above. The names used here for these geometries constitute only examples. Any other experimental procedure used to determine based on the ratio of two cavity doses can be simulated with this application.
The correction for backscatter in air when no FAC is present, 
, needs an extra calculation of the air-kerma at the POM. In the initial implementation this calculation was done separately with the application cavity. Now it can be done in the same calculation with egs_fac . To this end, one must define a geometry without the FAC geometry and allow the scoring of the contribution from scattered photons to the air-kerma at the POM. This is acomplished by setting the input key include scatter to yes as shown in the example above in the calculation geometry definition block for the geometry without the chamber sim_no_fac. The correction factor is then obtained as the ratio of the air-kerma at the POM for both geometries, with and without the FAC. One should allow for enough air behind the POM to avoid underestimating the backscatter.
For the computation of the correction factors, several quantities need to be scored. These quantities are mostly energy deposited in the collecting volume (CV) by electrons originating from different types of photons (primary, scattered or leaked through the aperture). A detailed description of these quantities is given in the output file *.egslog as shown below:
...
The energies listed below have the following definitions:
Eg = Total energy deposited in the collecting volume (CV)
E1 = Total energy deposited in CV excluding energy from particles that have
= visited an aperture region
E2 = Energy deposited in the CV from primary electrons
E3 = E2 corrected for energy loss/gain through the side faces of the CV
E4 = E3 corrected for energy loss/gain through the front/back CV faces
E4a = Same as E4 but computed from a kerma approximation for photons
= passing through the CV
E5 = E4a for an ideal point source
E6 = E5 corrected for attenuation. E6 = collision kerma at POM
...
A correlated scoring of all these energy deposition ratios can be requested by linking any combination of two geometries from the simulation geometries previously defined in the input file. This is done by assigning two geometry names to the correlated geometries key.
Variance reduction techniques
egs_fac implements the same variance reduction techniques used in cavity.
Usage
As any other EGSnrc application, egs_fac can be started from the command line using
egs_fac -i input_file -p pegs_file [-o output_file] [-b] [-s] [-P N -j i]
where the arguments in square brackets are optional. With the -o option one can change the name of the output files (by default input_file.xxx is used, where xxx is .egslog for the log file, .egsdat for the data file, etc.). With -b one specifies a "batch" run, i.e. the output is redirected to output_file.egslog. With -s one can force egs_fac to use a simple RCO instead of a JCF-RCO in parallel runs specified with -P N -j i, where N is the number of parallel jobs and i the job index. Note that one Unix-type systems it is easier to use the exb command to submit parallel jobs
exb egs_fac input_file pegs_file [p=N] [batch=pbs]
where the batch option specifies the queuing system to be used. The EGSnrc GUI can of course be also used, see see PIRS-877 for more details on running EGSnrc applications.
An input example: NRC MEES FAC corrections.
############################################################################
#
# This is an example input file for the egs_fac EGSnrc application.
#
# The incident beam is a full BEAMnrc simulation source that models the
# NRC MXR-320 x-ray tube.
#
# The simulation geometry is a detailed model of the NRC medium energy
# free air chamber and also includes simulation of the vacuum tube used
# to measure the attenuation correction.
#
# The simulation is set to compute all correction factors as described
# in Mainegra-Hing, Reynaert and Kawrakow, Med.Phys. 35(8),2008 3650-3660.
#
# See comments accompanying the various inputs for more details.
#
##############################################################################
:start geometry definition:
#
# The base structure of the simulation geometry will be made using a cone stack.
# The FAC will then be inscribed into one of the regions of the cone stack.
#
#
# We start by defining the lead box surrounding the FAC. The definition
# that follows will only define the side walls of the lead box and the inside
# air, the resulting geometry being infinite in z-direction (the direction of the
# beam). We will get the top/bottom lead box plates, including the aperture opening,
# when this geometry is inscribed into the cone stack mentioned above.
#
#
# Lead box x-planes
#
:start geometry:
library = egs_planes
type = EGS_Xplanes
name = lead_box_xplanes
positions = -28.7 28.7
:stop geometry:
#
# Lead box y-planes
#
:start geometry:
library = egs_planes
type = EGS_Yplanes
name = lead_box_yplanes
positions = -23.7 23.7
:stop geometry:
#
# The lead box, infinite along z
#
:start geometry:
library = egs_ndgeometry
name = lead_box_xy
dimensions = lead_box_xplanes lead_box_yplanes
:start media input:
media = PBICRU512
:stop media input:
:stop geometry:
#
# We now make an air box, also infinite along z, filled with air.
# When we inscribe this air box into the lead box defined above,
# we will get the lead side walls of the FAC.
#
:start geometry:
library = egs_planes
type = EGS_Xplanes
name = air_xplanes
positions = -27.5 27.5
:stop geometry:
:start geometry:
library = egs_planes
type = EGS_Yplanes
name = air_yplanes
positions = -22.5 22.5
:stop geometry:
:start geometry:
library = egs_ndgeometry
name = inside_air
dimensions = air_xplanes air_yplanes
:start media input:
media = AIRICRU512
:stop media input:
:stop geometry:
#
# The air filled lead box, infinite along z.
#
:start geometry:
library = egs_genvelope
name = main_lead_box
base geometry = lead_box_xy
inscribed geometries = inside_air
:stop geometry:
#
# The aluminum guard bars of the FAC are defined in a way similar to the above.
# We first define a box filled with aluminum, then a box filled with air, and then
# inscribe the air box into the aluminum box. We wll add the air spaces
# between the guard bars and the top/bottom openings for the beam later on.
# Note that the box is again infinite in z-direction.
#
# Aluminum box x-planes
#
:start geometry:
library = egs_planes
type = EGS_Xplanes
name = guard_outer_xplanes
positions = -17.5 17.5
:stop geometry:
#
# Aluminum box y-planes
#
:start geometry:
library = egs_planes
type = EGS_Yplanes
name = guard_outer_yplanes
positions = -14.0 14.0
:stop geometry:
#
# Aluminum box
#
:start geometry:
library = egs_ndgeometry
name = guard_outer_xy
dimensions = guard_outer_xplanes guard_outer_yplanes
:start media input:
media = ALICRU512
:stop media input:
:stop geometry:
#
# Air box x-planes
#
:start geometry:
library = egs_planes
type = EGS_Xplanes
name = guard_inner_xplanes
positions = -16.5475 16.5475
:stop geometry:
#
# Air box y-planes
#
:start geometry:
library = egs_planes
type = EGS_Yplanes
name = guard_inner_yplanes
positions = -13.0 13.0
:stop geometry:
#
# Air box
#
:start geometry:
library = egs_ndgeometry
name = guard_inner_xy
dimensions = guard_inner_xplanes guard_inner_yplanes
:start media input:
media = AIRICRU512
:stop media input:
:stop geometry:
#
# The air filled aluminum box
#
:start geometry:
library = egs_genvelope
name = guard_xy
base geometry = guard_outer_xy
inscribed geometries = guard_inner_xy
:stop geometry:
#
# The following will be used to define the air gaps in the aluminum guard bars
#
:start geometry:
library = egs_planes
type = EGS_Yplanes
name = guard_y_divisions
positions = -13.5,
-12.0475 -11.911,
-10.95847826 -10.82195652,
-9.869456522 -9.732934783,
-8.780435 -8.643913,
-7.691413043 -7.5548913043,
-6.60239130434 -6.465869565,
-5.513369565 -5.376847826,
-4.424347826 -4.287826087,
-3.335326087 -3.198804348,
-2.246304348 -2.109782609,
-1.157282609 -1.02076087,
-0.06826086957 0.06826086957,
1.02076087 1.157282609,
2.109782609 2.246304348,
3.198804348 3.335326087,
4.287826087 4.424347826,
5.376847826 5.513369565,
6.465869565 6.60239130434,
7.5548913043 7.691413043,
8.643913 8.780435,
9.732934783 9.869456522,
10.82195652 10.95847826,
11.911 12.0475,
13.5
:start media input:
media = ALICRU512 AIRICRU512
:start input loop:
loop count = 23
loop variable = 0 region 1 2
set medium = $(region) 1
:stop input loop:
#
# The above input loop will expand into a
# a series of 23 inputs of the form
# set medium = 1 1
# set medium = 3 1
# sed medium = 5 1
# ...
# which will set every second region of the guard_y_divisions geometry
# to medium 1 (i.e., to air).
:stop media input:
:stop geometry:
#
# We now use a so called "smart" envelope geometry to inscribe the
# above planes into region 0 of the air filled aluminum box (i.e., into
# the aluminum walls. That way the previously solid aluminum walls will get
# air gaps as defined by the guard_y_divisions geometry. We use a logic type 1,
# which means that the geometry being inscribed does not fit completely inside
# the specified region and forces the smart envelope geometry to do boundary
# checks against the region 0 boundaries (in contrast, a logic type 0 inscription
# means that the geometry being inscribed fits completely in the specified region
# so that when a particle is in a region of the inscribed geometry one does not need
# to check for intersections with the envelope region.
#
:start geometry:
library = egs_smart_envelope
name = guard_xy_divided
base geometry = guard_xy
inscribe geometry = guard_y_divisions 0 1
:stop geometry:
#
# We now need the square top/bottom beam openings of the
# guards. We also need to air the air gaps to the top/bottom parts of the guards.
# We define the openings as 2D geometries made out of x- and y-planes
# (again infinite along z) filled with air. These are inscribed into the
# infinite aluminum box guard_outer_xy defined above. We then use again a smart
# envelope to get the air gaps by inscribing the guard_y_divisions geometry.
# Because the definition is equivalent for the top and bottom
# guard parts (apart from the different opening size), we use an input loop to save some typing.
#
:start input loop:
loop count = 2
loop variable = 2 planes 1 2
loop variable = 2 name top bottom
:start geometry:
library = egs_planes
type = EGS_Xplanes
name = $(name)_opening_xplanes
positions = -$(planes) $(planes)
:stop geometry:
:start geometry:
library = egs_planes
type = EGS_Yplanes
name = $(name)_opening_yplanes
positions = -$(planes) $(planes)
:stop geometry:
:start geometry:
library = egs_ndgeometry
name = $(name)_guard_opening
dimensions = $(name)_opening_xplanes $(name)_opening_yplanes
:start media input:
media = AIRICRU512
:stop media input:
:stop geometry:
:start geometry:
library = egs_genvelope
name = guard_$(name)
base geometry = guard_outer_xy
inscribe in regions = $(name)_guard_opening 0
:stop geometry:
:start geometry:
library = egs_smart_envelope
name = guard_$(name)_divided
base geometry = guard_$(name)
inscribe geometry = guard_y_divisions 0 1
:stop geometry:
:stop input loop:
#
# The above, when expanded by the EGSnrc input reader will give us
# geometries named guard_top_divided and guard_bottom_divided.
#
#
# To define the actual collecting volume (CV), we use a CD geometry
# to inscribe the top/bottom/side guard bars into a set of z-planes.
#
# First the z-planes.
#
:start geometry:
library = egs_planes
type = EGS_Zplanes
name = CV_zplanes
positions = 13.35 14.3025 37.35 47.35 70.3975 71.35
# ^ ^ ^ ^ ^
# top guard bars air not part of CV CV air bottom
:stop geometry:
:start geometry:
library = egs_cdgeometry
name = guard_plus_CV
base geometry = CV_zplanes
new indexing style = 1
set geometry = 0 guard_top_divided
set geometry = 1 guard_xy_divided
set geometry = 2 guard_xy_divided
set geometry = 3 guard_xy_divided
set geometry = 4 guard_bottom_divided
:stop geometry:
#
# We now put the above CV plus guard bars geometry into the air filled lead
# box using a smart envelope. Note inscription logic 0, which is now applicable
# because the geometry being inscribed fits completely into region 1 of the
# air filled lead box.
# Note that the resulting main_fac geometry is still infinite along the z-axis.
# It will belimited and will get its top/bottom walls + aperture when it gets
# inscribed into the cone stack geometry below.
#
:start geometry:
library = egs_smart_envelope
name = main_fac
base geometry = main_lead_box
inscribe geometry = guard_plus_CV 1 0
:stop geometry:
#
# The following input loop defines 3 cone stack geometries named
# fac_no_tube, fac_air_tube and fac_vacuum_tube for the 3 cases
# we want to simulate at once. These cone stacks form the "skeleton"
# of the final simulation geometries where the lead box containing
# the CV and guard bars will be inscribed.
#
:start input loop:
loop count = 3 # i.e., execute the loop 2 times.
# a loop variable of type 2 is simply a list of possible values
loop variable = 2 gname base_no_tube base_air_tube base_vacuum_tube
loop variable = 2 med_be AIRICRU512 BEICRU512 BEICRU512
loop variable = 2 med_tube_wall AIRICRU512 SS3161ICRU512 SS3161ICRU512
loop variable = 2 med_tube AIRICRU512 AIRICRU512 THINAIR512
loop variable = 2 bname fac_no_tube fac_air_tube fac_vacuum_tube
:start geometry:
library = egs_cones
type = EGS_ConeStack
name = $(gname)
axis = 0 0 -51.2 0 0 1
#
# Regions 0-2
# This is used to define the Be window of the vacuum or air filled tube
# when present. When the input loop gets expanded, $(med_be) will be
# replaced by AIRICRU512, BEICRU512 and BEICRU512
#
:start layer:
thickness = 0.1
top radii = 3.2 40.0
bottom radii = 3.2 40.0
media = $(med_be) AIRICRU512
:stop layer:
#
# Regions 3-5
# This is used to define the vacuum/air-filled tube when present.
# We will get for the media
# AIRICRU512 AIRICRU512 AIRICRU512
# AIRICRU512 SS3161ICRU512 AIRICRU512
# THINAIR512 SS3161ICRU512 AIRICRU512
# in the 3 loop cases. We use THINAIR512 (instead of simply vacuum), because
# there is a bug still somewhere that causes a crash when vacuum is used.
# THINAIR512 is a medium with 1e-6 times lower density than air.
#
:start layer:
thickness = 41.9
top radii = 3.0 3.2 40.0
bottom radii = 3.0 3.2 40.0
media = $(med_tube) $(med_tube_wall) AIRICRU512
:stop layer:
#
# Regions 6-8
# Be exit window when present
:start layer:
thickness = 0.1
top radii = 3.0 40.0
bottom radii = 3.0 40.0
media = $(med_be) AIRICRU512
:stop layer:
#
# Regions 9-11
# Air between tube and FAC aperture
:start layer:
thickness = 9.1
top radii = 40.0
bottom radii = 40.0
media = AIRICRU512
:stop layer:
#
# Regions 12-14
# Part of the NRC FAC aperture.
:start layer:
thickness = 0.15
top radii = 0.65 2 40.0
bottom radii = 0.5 2 40.0
media = AIRICRU512 HEVIMETICRU512 AIRICRU512
:stop layer:
#
# Regions 15-17
# Part of the NRC FAC aperture.
:start layer:
thickness = 0.3
bottom radii = 0.5 2 40.0
media = AIRICRU512 HEVIMETICRU512 AIRICRU512
:stop layer:
#
# Regions 18-20
# Part of the NRC FAC aperture.
:start layer:
thickness = 0.8
bottom radii = 0.57 2 40.0
media = AIRICRU512 HEVIMETICRU512 AIRICRU512
:stop layer:
#
# Regions 21-23
# The inner region will remain air for the FAC opening, the outer region
# will turn into the top lead box wall when the lead box is inscribed.
:start layer:
thickness = 2.1
top radii = 1.25 40
bottom radii = 1.25 40
media = AIRICRU512 AIRICRU512
:stop layer:
#
# Regions 24-26
# Here is where the FAC will go.
:start layer:
thickness = 78
top radii = 40
bottom radii = 40
media = AIRICRU512
:stop layer:
#
# Regions 27-29
# The inner region will remain air for the bottom FAC opening, the outer region
# will turn into the bottom lead box wall when the lead box is inscribed.
:start layer:
thickness = 2.1
top radii = 2.25 40
bottom radii = 2.25 40
media = AIRICRU512 AIRICRU512
:stop layer:
:stop geometry:
#
# These are the final simulation geometries.
#
:start geometry:
library = egs_smart_envelope
name = $(bname)
base geometry = $(gname)
inscribe geometry = lead_box_xy 22 1
inscribe geometry = main_fac 24 1
inscribe geometry = lead_box_xy 28 1
:stop geometry:
:stop input loop:
###################################
# The cavity for range rejection
# We make it slightly larger (0.01 cm) to avoid fat
# particles entering the CV due to numerical round off issues.
#
:start geometry:
library = egs_ndgeometry
type = EGS_XYZGeometry
name = cavity
x-planes = -16.55 16.55
y-planes = -13.01 13.01
z-planes = 37.349 47.351
:stop geometry:
simulation geometry = fac_no_tube
:stop geometry definition:
#
# The BEAMnrc simulation source
# The specific input we used was using DBS with a splitting radious of 13 cm
# to avoid the possibilities of having a rare interaction in air of a fat
# particles destroying the statistics. We use the "cutout" option of the
# BEAM simulation source to limit the beam to a 1.6 x 1.6 cm square
# (larger beams give the same results as nothing gets through the lead walls).
# We also don't simulate contaminant electrons from the x-ray tube and
# have a weight window to throw away fat particles outside of the
# DBS splitting field.
#
:start source definition:
:start source:
library = egs_beam_source
name = the_source
beam code = BEAM_mxr-320
pegs file = xcomx
input file = mxr-320-250kVp_tube
particle type = photons
cutout = -1.6 1.6 -1.6 1.6
weight window = 0 0.01
:stop source:
simulation source = the_source
:stop source definition:
:start run control:
ncase = 1600000000
:stop run control:
:start rng definition:
type = ranmar
initial seeds = 1802 4000
:stop rng definition:
:start scoring options:
#
# The definition of calculation geometries is very similar to the
# cavity code, but there are inputs added that define the aperture
# regions (needed to get the Aap correction), the front/back
# CV regions (needed for the electron loss calculation) and the
# point of measurement (POM), needed for KERMA scoring.
# In our simulation source the geometry is defined so that the
# z-position of the particles leaving at the bottom of the geometry is
# 48.35 cm, which coincides with the top of the egs_fac simulation
# geometry. To make these two self consistent, we apply a translation
# of -99.55 cm along z before transporting source particles.
#
#
# The FAC without the air/vacuum filled tube used to measre Aatt.
#
:start calculation geometry:
geometry name = fac_no_tube
cavity regions = 132
aperture regions = 13 16 19 30
front and back regions = 83 181
cavity mass = 0.009462477073 # cylinder defined by diaphragm and plates
:start transformation:
translation = 0 0 -99.55
:stop transformation:
POM = 0.45 0.5 # first input is z-position, second is radius.
:stop calculation geometry:
#
# The FAC wit the tube filled with air.
#
:start calculation geometry:
geometry name = fac_air_tube
cavity regions = 132
aperture regions = 13 16 19 30
front and back regions = 83 181
cavity mass = 0.009462477073 # cylinder defined by diaphragm and plates
:start transformation:
translation = 0 0 -99.55
:stop transformation:
POM = 0.45 0.5
:stop calculation geometry:
#
# The FAC wit the tube filled with vacuum.
#
:start calculation geometry:
geometry name = fac_vacuum_tube
cavity regions = 132
aperture regions = 13 16 19 30
front and back regions = 83 181
cavity mass = 0.009462477073 # cylinder defined by diaphragm and plates
:start transformation:
translation = 0 0 -99.55
:stop transformation:
POM = 0.45 0.5
:stop calculation geometry:
#
# The following will output the ratio of the CV dose for the
# 2 geometries, including an uncertainty estimate that correctly
# takes into account the correlation. This will be the measured
# Aatt correction
#
correlated geometries = fac_air_tube fac_vacuum_tube
#
# To compute the Ax correction factor defined by Mainegra-Hing, Reynaert
# and Kawrakow, one must tell egs_fac the 3 geometries involved (i.e.,
# the actual KERMA measurement geometry, the geometry with a vacuum
# tube filled with air and the geometry with a vacuum tube filled with vacuum.
#
Ax calculation = fac_no_tube fac_air_tube fac_vacuum_tube
#
# One must also provide a tabulation for mu_en/rho of air.
# The format of this file is
# number N of data points
# N lines with energy and mu_en/rho
# Note that although the energies are given, at this point egs_fac
# assumes that the grid points are uniform in log(E) and actually
# ignores them (except for the first and last energy, which define the
# energy interval).
#
muen file = /home/iwan/egsnrc_mp/cavity/muen_air.data
:stop scoring options:
:start variance reduction:
#
# Split each photon 1000 times.
# The optimum splitting number will depend on how much time
# is spent in the source and in the egs_fac simulation. In this
# example the source is a BEAMnrc simulation source, which is quite
# slow for x-ray tube simulations, and therefore a large splitting number
# is appropriate. If you start with a simnle source (e.g. point source
# with a spectrum), splitting numbers in the range 10-100 will give better
# efficiency
#
photon splitting = 1000
:start range rejection:
#
# rejection = 1000 means that all electrons that can not
# reach the CV will be subjected to Russian Roulette with
# survival probability of 1/1000
#
rejection = 1000
Esave = 1
#
# if you don't define a cavity geometry, range rejection/Russian Roulette
# will be on a region-by-region basis only
#
cavity geometry = cavity
#
# The medium to use for range calculation. This must be the medium with
# the lowest stopping power in the simulation geometry
#
rejection range medium = AIRICRU512
:stop range rejection:
:stop variance reduction:
:start MC transport parameter:
Rayleigh scattering = On
#
# Note: using 'Bound Compton Scattering = norej' is essential!!!
# The Bound Compton Scattering = On' option will not work because
# the scoring does not take into account the possibility of rejection
# of Compton interactions when computing the attenuation. This is explained
# in the Mainegra-Kawrakow Med.Phys. paper.
#
Bound Compton Scattering = norej
Photon cross sections = xcom # (could also be epdl)
:stop MC transport parameter:
