|
EGSnrc C++ class library
Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
|
The C++ application egs_gammaspec is an advanced EGSnrc application designed for the purpose of calculating gamma spectrometry related quantities including total efficiency, full energy peak efficiency, individual peak efficiencies and summing corrections.
Many of the results are divided into two categories: a 'perfect', and 'non-perfect' style of detector modelling. In this context, 'perfect' means that the scoring regions record energy depositions of every interaction as independently ("perfectly") resolved events. There is no dead time. This is physically unrealistic, because all real detectors are susceptible to detecting multiple different particles that interact at nearly the same time, and summing those contributions as though they came from a single particle. This is the case for the 'non-perfect' model of energy depositions. In 'non-perfect' mode, all energy depositions from a single decay event (all the secondaries that interact in the sensitive region) are combined into a single deposition event. This means that the produced 'non-perfect' spectrum will contain summed peaks, where the energy recorded in the spectrum for an event is the sum of the deposited energy from multiple particles originating from the same decay shower. Note that each decay is still perfectly resolved - there is no cross-contamination between decays. Therefore the amount of summing may still be underestimated for high activity radionuclides. This is left for future work.
This application generally assumes that you are using egs_radionuclide_source to produce decays from a radionuclide. While it is possible to use a different source, the 'perfect' and 'non-perfect' results will be the same unless the source generates secondary particles. For example: when using a point source, all source particles are primaries; with a radionuclide source, one primary decay event may result in many secondary particles (e.g. a series of internal transitions that generate particles).
The spectra are calculated as the energy deposited in each bin, normalized to the number of decays (i.e. disintegration events). The total efficiency is the sum of the spectrum bins. The full energy peak efficiency is the spectrum maximum value.
Individual peak efficiencies are calculated using the bins nearest the peak. If the peak is between two bins, then those bins are added together. The values of the bins from ind1 to ind2+1 are combined, and the background is subtracted. Consider the peak i with energy gammaEnergies[i], and an energy bin width of binWidth. Then we define
ind1 = ceil(gammaEnergies[i]/binWidth-1); ind2 = ceil(gammaEnergies[i]/binWidth-0.5);
where ceil is the ceiling function. Background is subtracted from peak efficiency calculations by subtracting the average of the bins nearby the gamma peak energies (as provided by the user, or the automatic gamma analysis energies). For the spectrum array spectr,
background = (spectr[ind1-1] + spectr[ind2+1])/2;
The summing correction for each peak is just the ratio of the 'perfect' and 'non-perfect' peak efficiency. For each peak i,
summingCorrection = peakEfficiency_perf[i] / peakEfficiency[i];
Here are some future ideas that are not yet implemented:
All of the special inputs for egs_gammaspec are in the scoring options. The scoring options and output spectrum input blocks are required for the simulation to run. Only one output spectrum input block will be used - multiple are not supported. Output spectra can be produced and written to a file, and the regions selected for scoring are also used for calculating efficiencies and summing corrections.
:start scoring options:
:start output spectrum:
scoring regions = crystal_no_dead_label # The region numbers or region label(s) denoting the sensitive regions
minimum spectrum energy = 0.0 # Optional, MeV, default=0, minimum energy in output spectrum
maximum spectrum energy = 0.4 # Optional, MeV, defaults=maximum in source, maximum energy in output spectrum
number of bins = 2000 # Optional, default=1000, number of bins in output spectrum
# Optional, yes or no, default=yes, get the analysis energies automatically, using all of the gamma energies from the radionuclide decay scheme. They will be combined with the manually entered 'gamma analysis energies', so make sure they don't overlap. Only works for egs_radionuclide_source.
automatic analysis energies = no
# Example: these are peaks of interest for Ba-133, in MeV
# Optional, coincidence summing corrections will be calculated for each of these
gamma analysis energies = .0309 .035 .0531 .0796 .0810 .1606 .2232 .2764 .3028 .356 .3838
# Optional, MeV, default=1e-6, the minimum energy from a single shower that can be detected in the sensitive regions. If the total energy deposited by the shower is less than this value, it is not added to the recorded spectrum.
minimum detectable energy = 1e-6
:stop output spectrum:
:stop scoring options:
The output files are automatically named ${basename}-spec.txt and ${basement}-spec-perf.txt, for the 'non-perfect' and 'perfect' spectra, respectively. ${basename} corresponds to the name of the input file, without the .egsinp extension. The format of the output spectrum is simply three columns of floating point numbers. The first column contains the energy (MeV) of the middle of the energy bin (all bins are equal width). The second column contains the energy deposited per decay (MeV per decay). The third is the absolute uncertainty on the energy deposited (MeV per decay). For example:
0.000000e+00 1.200985e-04 8.491990e-05
2.000000e-04 0.000000e+00 0.000000e+00
4.000000e-04 6.004924e-05 6.004924e-05
... etc...
###############################################################################
#
# EGSnrc egs++ ICRM HPGE detector
# Copyright (C) 2019 National Research Council Canada
#
# EGSnrc is free software: you can redistribute it and/or modify it under
# the terms of the GNU Affero General Public License as published by the
# Free Software Foundation, either version 3 of the License, or (at your
# option) any later version.
#
# EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
# more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
#
###############################################################################
#
# Authors: Reid Townson
#
# Contributors:
#
###############################################################################
#
# Input file for simulation of an ICRM HPGE detector
# Used for the GSWG Coincidence Summing Comparison
#
# This input file is compatible with the egs_gammaspec application.
# To use, place this input file in your egs_home/egs_gammaspec directory and
# run with egs_gammaspec -i inputfile
#
# Many options are commented out, that were also used in the comparison.
# Currently, the uncommented geometries and sources correspond to a simulation
# of a Ba-133 source with "filter" material in the "A" configuration (with
# a dead layer in the detector).
#
###############################################################################
##############################################################################
### RUN CONTROL
##############################################################################
:start run control:
ncase = 1e9
:stop run control:
:start rng definition:
initial seeds = 7643 1293
:stop rng definition:
##############################################################################
### GEOMETRY
##############################################################################
:start geometry definition:
# Lead shielding
:start geometry:
name = shielding
library = egs_cones
type = EGS_ConeStack
# Detector is positioned so that the top of the window
# is at the origin, and "down" is in the -z direction
axis = 0 0 15 0 0 -1
# Top of the lead shielding
:start layer:
thickness = 5
top radii = 20
bottom radii = 20
media = lead
:stop layer:
:start layer:
thickness = 30
top radii = 15 20
bottom radii = 15 20
media = air lead
:stop layer:
# Bottom of the lead shielding
:start layer:
thickness = 5
top radii = 20
bottom radii = 20
media = lead
:stop layer:
:stop geometry:
##################################
# Detector A
# The detector WITH dead layer
##################################
:start geometry:
name = detector
library = egs_cones
type = EGS_ConeStack
# Detector is positioned so that the top of the window
# is at the origin, and "down" is in the -z direction
axis = 0 0 0 0 0 -1
# Detector window
:start layer:
thickness = 0.1
top radii = 4
bottom radii = 4
media = aluminum
:stop layer:
# Gap above dead layer
:start layer:
thickness = 0.5
top radii = 3.9 4
bottom radii = 3.9 4
media = vacuum aluminum
:stop layer:
# Top dead layer
:start layer:
thickness = 0.1
top radii = 3 3.9 4
bottom radii = 3 3.9 4
media = germanium vacuum aluminum
:stop layer:
# Crystal above the void
:start layer:
thickness = 1.9
top radii = 2.9 3 3.9 4
bottom radii = 2.9 3 3.9 4
media = germanium germanium vacuum aluminum
:stop layer:
# Crystal with the void
:start layer:
thickness = 4
top radii = 0.5 2.9 3 3.9 4
bottom radii = 0.5 2.9 3 3.9 4
media = vacuum germanium germanium vacuum aluminum
:stop layer:
# Gap below crystal
:start layer:
thickness = 1.3
top radii = 3.9 4
bottom radii = 3.9 4
media = vacuum aluminum
:stop layer:
# Bottom housing
:start layer:
thickness = 0.1
top radii = 4
bottom radii = 4
media = aluminum
:stop layer:
set label = crystal_with_dead_label 15 21
:stop geometry:
##################################
# Detector B
# The detector WITHOUT dead layer
##################################
# :start geometry:
name = detector
library = egs_cones
type = EGS_ConeStack
# Detector is positioned so that the top of the window
# is at the origin, and "down" is in the -z direction
axis = 0 0 0 0 0 -1
# Detector window
:start layer:
thickness = 0.1
top radii = 4
bottom radii = 4
media = aluminum
:stop layer:
# Gap above crystal
:start layer:
thickness = 0.5
top radii = 3.9 4
bottom radii = 3.9 4
media = vacuum aluminum
:stop layer:
# Crystal above the void
:start layer:
thickness = 2
# thickness = 2.1 # This is was mistakenly used previously
top radii = 3 3.9 4
bottom radii = 3 3.9 4
media = germanium vacuum aluminum
:stop layer:
# Crystal with the void
:start layer:
thickness = 4
top radii = 0.5 3 3.9 4
bottom radii = 0.5 3 3.9 4
media = vacuum germanium vacuum aluminum
:stop layer:
# Gap below crystal
:start layer:
thickness = 1.3
top radii = 3.9 4
bottom radii = 3.9 4
media = vacuum aluminum
:stop layer:
# Bottom housing
:start layer:
thickness = 0.1
top radii = 4
bottom radii = 4
media = aluminum
:stop layer:
set label = crystal_no_dead_label 8 13
:stop geometry:
############################
# Source geometries
############################
# The water source geometry
:start geometry:
name = water_source_geom
library = egs_cones
type = EGS_ConeStack
# Source is positioned 1mm above window
axis = 0 0 4.1 0 0 -1
# Top of the lead shielding
:start layer:
thickness = 4
top radii = 4.5
bottom radii = 4.5
media = water
:stop layer:
:stop geometry:
# The soil source geometry
:start geometry:
name = soil_source_geom
library = egs_cones
type = EGS_ConeStack
# Source is positioned 1mm above window
axis = 0 0 2.1 0 0 -1
# Top of the lead shielding
:start layer:
thickness = 2
top radii = 3
bottom radii = 3
media = dirt
:stop layer:
:stop geometry:
# The filter source geometry
:start geometry:
name = filter_source_geom
library = egs_cones
type = EGS_ConeStack
# Source is positioned 1mm above window
axis = 0 0 0.4 0 0 -1
# Top of the lead shielding
:start layer:
thickness = 0.3
top radii = 4
bottom radii = 4
media = cellulose
:stop layer:
:stop geometry:
##########################
# Combining geometries
##########################
# Put the detector in the shielding (for point source)
:start geometry:
name = detector_pointsource_shielding
library = egs_genvelope
base geometry = shielding
inscribed geometries = detector
:stop geometry:
# Put the water source and detector in the shielding
:start geometry:
name = detector_watersource_shielding
library = egs_genvelope
base geometry = shielding
inscribed geometries = detector water_source_geom
:stop geometry:
# Put the soil source and detector in the shielding
:start geometry:
name = detector_soilsource_shielding
library = egs_genvelope
base geometry = shielding
inscribed geometries = detector soil_source_geom
:stop geometry:
# Put the filter source and detector in the shielding
:start geometry:
name = detector_filtersource_shielding
library = egs_genvelope
base geometry = shielding
inscribed geometries = detector filter_source_geom
:stop geometry:
simulation geometry = detector_filtersource_shielding
:stop geometry definition:
##############################################################################
### MEDIA
##############################################################################
:start media definition:
# Cutoff energies
ae = 0.512 # lowest energy for electrons (kinetic+0.511)
ap = 0.001 # lowest energy for photons (kinetic)
# NOTE: You might want to adjust these if you change the radionuclide
ue = 3.512 # maximum energy for electrons (kinetic+0.511)
up = 3.001 # maximum energy for photons (kinetic)
:start air:
density correction file = air_dry_nearsealevel
:stop air:
:start aluminum:
density correction file = aluminum
:stop aluminum:
:start lead:
density correction file = lead
:stop lead:
:start germanium:
density correction file = germanium
:stop germanium:
:start water:
density correction file = water_icru90
:stop water:
:start cellulose:
density correction file = icrm_cellulose
:stop cellulose:
:start dirt:
density correction file = icrm_dirt
:stop dirt:
:stop media definition:
##############################################################################
### SOURCE
##############################################################################
:start source definition:
# NOTE: When you change source energy also change scoring energy max
# Point source
# :start source:
name = point_source_iso
library = egs_point_source
# Positioned 1mm above detector window
position = 0 0 0.1
charge = 0
:start spectrum: # not used when part of radionuclide source
type = monoenergetic
energy = 1
:stop spectrum:
:stop source:
# :start source:
name = point_source
library = egs_radionuclide_source
base source = point_source_iso
:start spectrum:
type = radionuclide
nuclide = Ba-133
atomic relaxations = eadl
extra transition approximation = On
:stop spectrum:
:stop source:
# Water source
# :start source:
name = water_source_iso
library = egs_isotropic_source
charge = 0
geometry = water_source_geom
:start shape:
type = cylinder
radius = 4.5
height = 4
axis = 0 0 1
midpoint = 0 0 2.1
:stop shape:
:start spectrum: # not used when part of radionuclide source
type = monoenergetic
energy = 1
:stop spectrum:
:stop source:
# :start source:
name = water_source
library = egs_radionuclide_source
base source = water_source_iso
:start spectrum:
type = radionuclide
nuclide = Ba-133
atomic relaxations = eadl
extra transition approximation = On
:stop spectrum:
:stop source:
# Soil source
# :start source:
name = soil_source_iso
library = egs_isotropic_source
charge = 0
geometry = soil_source_geom
:start shape:
type = cylinder
radius = 3
height = 2
axis = 0 0 1
midpoint = 0 0 1.1
:stop shape:
:start spectrum: # not used when part of radionuclide source
type = monoenergetic
energy = 1
:stop spectrum:
:stop source:
# :start source:
name = soil_source
library = egs_radionuclide_source
base source = soil_source_iso
:start spectrum:
type = radionuclide
nuclide = Ba-133
atomic relaxations = eadl
extra transition approximation = On
:stop spectrum:
:stop source:
# Filter source
:start source:
name = filter_source_iso
library = egs_isotropic_source
charge = 0
geometry = filter_source_geom
:start shape:
type = cylinder
radius = 4
height = 0.3
axis = 0 0 1
midpoint = 0 0 0.25
:stop shape:
:start spectrum: # not used when part of radionuclide source
type = monoenergetic
energy = 1
:stop spectrum:
:stop source:
:start source:
name = filter_source
library = egs_radionuclide_source
base source = filter_source_iso
:start spectrum:
type = radionuclide
nuclide = Ba-133
atomic relaxations = eadl
extra transition approximation = On
:stop spectrum:
:stop source:
simulation source = filter_source
:stop source definition:
##############################################################################
### AUSGAB
##############################################################################
:start ausgab object definition:
# Score tracks
#:start ausgab object:
name = tracks
library = egs_track_scoring
stop scoring = 2000
:stop ausgab object:
:stop ausgab object definition:
##############################################################################
### SCORING OPTIONS
##############################################################################
# For egs_gammaspec
:start scoring options:
:start output spectrum:
# scoring regions = crystal_no_dead_label
scoring regions = crystal_with_dead_label
number of bins = 2000
# These are peaks of interest for Ba-133, in MeV
# Efficiency and coincidence summing corrections will be calculated for each of these
gamma analysis energies = .0309 .035 .0531 .0796 .0810 .1606 .2232 .2764 .3028 .356 .3838
automatic analysis energies = no
maximum spectrum energy = 0.4
:stop output spectrum:
:stop scoring options:
##############################################################################
### TRANSPORT PARAMETERS
##############################################################################
:start MC transport parameter:
Global ECUT = 0.512 # Global electron transport cutoff
Global PCUT = 0.001 # Global photon transport cutoff
Spin effects = On # On (default),Off
Brems angular sampling = KM # Simple,KM (default)
Brems cross sections = NRC # BH (default),NIST
Pair angular sampling = KM # Off, Simple (default),KM
Bound Compton scattering = On # On (default) ,Off, norej
Radiative Compton corrections = On # On, Off (default)
Photoelectron angular sampling = On # On (default),Off
Atomic relaxations = On # On (default),Off
Photon cross sections = xcom # xcom (default), epdl, si
Rayleigh scattering = On
:stop MC transport parameter:The results expected for this example are as follows:
====================================================== Results output for egs_gammaspec: ====================================================== => last case = 166722516 fluence = 1.66722e+08 Total energy emitted from source = 7.54213e+07 MeV Total energy fraction recorded in raw non-perfect detector spectrum = 0.0738755 Total energy fraction recorded in raw perfect detector spectrum = 0.0757989 === Efficiency for 'non-perfect' detector === Total efficiency = 16.315441 % +- 0.019051 % Full energy peak efficiency = 3.738300 % +- 0.039300 % Gamma energy [MeV] | Peak efficiency (background subtracted) [%] | Uncertainty [%] 0.030900 0.001341 5.603749 0.035000 0.004274 2.532217 0.053100 0.031454 0.487754 0.079600 0.156215 0.213381 0.081000 2.124703 0.052834 0.160600 0.065348 0.349972 0.223200 0.037366 0.454619 0.276400 0.502761 0.109455 0.302800 1.257765 0.068715 0.356000 3.738013 0.039309 0.383800 0.567427 0.102553 === Efficiency for 'perfect' detector === Total efficiency = 17.101547 % +- 0.018706 % Full energy peak efficiency = 3.892380 % +- 0.039127 % Gamma energy [MeV] | Peak efficiency (background subtracted) [%] | Uncertainty [%] 0.030900 0.001631 4.802029 0.035000 0.005278 2.139837 0.053100 0.039016 0.431393 0.079600 0.195182 0.189358 0.081000 2.542644 0.048689 0.160600 0.070338 0.335010 0.223200 0.041473 0.425152 0.276400 0.554652 0.104325 0.302800 1.319757 0.067423 0.356000 3.892176 0.039134 0.383800 0.525332 0.106814 === Coincidence summing correction === Gamma energy [MeV] | Summing correction (perfect/non-perfect) | Uncertainty [%] 0.030900 1.216779 7.379802 0.035000 1.234915 3.315272 0.053100 1.240425 0.651156 0.079600 1.249446 0.285285 0.081000 1.196706 0.071847 0.160600 1.076356 0.484471 0.223200 1.109925 0.622441 0.276400 1.103211 0.151209 0.302800 1.049287 0.096268 0.356000 1.041242 0.055467 0.383800 0.925815 0.148075 Spectrum scoring regions: 21 27 ====================================================== Start of source emissions statistics: Sampled Ba-133 emissions: ======================== Energy | Intensity per 100 decays (adjusted by -0.900500) Beta records: 0.000000 0.000000 0.000000 0.125928 0.000000 0.050653 0.000000 4.412178 0.000000 26.206155 Gamma records (E,Igamma,Ice,Ipp): 0.080998 10.275583 1.7475e+01 0.0000e+00 0.079614 0.827862 1.4188e+00 0.0000e+00 0.160612 0.200889 5.3947e-02 0.0000e+00 0.223237 0.138003 1.3330e-02 0.0000e+00 0.302851 5.620332 2.4794e-01 0.0000e+00 0.383848 2.739988 5.0967e-02 0.0000e+00 0.053162 0.643439 3.6930e+00 0.0000e+00 0.276399 2.201306 1.2311e-01 0.0000e+00 0.356013 19.149180 4.6890e-01 0.0000e+00 Average gamma energy: 0.000000 End of source emissions statistics ======================================================
1.9.1