Feature renormalized photoelectric cross sections and ICRU90 recommended values for key data #322

mainegra · 2017-08-08T20:15:46Z

This is a major PR introducing several important features:

Possibility to simulate the photoelectric effect (PE) using Sabbatucci and Salvat's re-normalized cross sections. These are numerically similar to Scofield's 1973 renormalized values (within 0.05%) except that this compilation uses additional energy points around the ionization edges for accurate linear log-log interpolation.
Ability to generate stopping powers using ICRU90 recommended I values for water and graphite.
The choice of atomic relaxation model is now requested via input file rather than via the macro $EADL_RELAX which is now deprecated.
For consistency with the incoherent scattering default of including binding effects, coherent scattering is now ON by default. Although these effects are negligible in the MV energy range, they have a small impact in simulation time. For instance for a ion chamber simulation in a water phantom irradiated with a Co-60 photon beam the time penalty is about 6% while for a 6 MV beam it is 4%. If this is significant (very long production runs) then the user may turn OFF Rayleigh scattering and use the Klein-Nishina model for incoherent scattering.
EGSnrc aborts if the user requests bound Compton scattering or electron impact ionization without atomic relaxations. In the past, atomic relaxations were turned ON automatically and a warning issued. However relaxations were only turned ON for the specific interaction and not for the photoelectric effect. Now users are forced to turn ON relaxation in their input files.

To use the features in this PR before they are merged with the develop branch users can get the branch feature-shellwise-pe-xsections and create a user configuration for it. See bottom of this PR for more details on this.

To use the re-normalized cross sections one will have to set the input key Photon cross sections in the MC transport parameter input block to one of the following values:

:start MC transport parameter:
  ...
  Photon cross sections = mcdf-xcom # or mcdf-epdl
  Atomic relxations = On # Off, EADL, simple
  ...
:stop MC transport parameter:

These options allow the use of the re-normalized PE cross sections while using the rest of the photon cross sections from either the XCOM or EPDL libraries. Note above that EADL relaxation is turned ON by using the values On or EADL. If the input is simple, then the original relaxation implementation is used. However, when mcdf-xcom or mcdf-epdl is used, no simple atomic relaxation is allowed.

The re-normalized cross sections are available for all atomic sub-shells up to N7 for elements from Z=1 to Z=99 for energies from the ionization threshold up to 1 GeV. However only sub-shells above 1 keV are included in EGSnrc. Taking advantage of the availability of sub-shell PE cross sections a new routine is implemented that allows sampling the photoelectric interaction from all energetically possible sub-shells. Since an adaptive energy grid is used near absorption edges, to ensure accurate interpolation, a binary search algorithm is used in the new photoelectric routine to find the energy interval the photon energy falls in. However, the total PE cross section used for sampling the interaction uses the same equidistant grid as the rest of the photon cross sections for fast interpolation. The time penalty of the new algorithm is about 6% for a BEAMnrc simulation of a 30 kV x-ray tube. However for an egs_fac simulation of HVL attenuation measurements with a free-air chamber (FAC) for 4 attenuation filters including the BEAMnrc x-ray tube simulation no time difference is observed. This can be attributed to the fact that more very low-energy particles are produced in the BEAMnrc simulation that don't make it to the FAC.

To create material data for graphite or water according to ICRU90 recommendations one will have to make use of the following density correction files:

$HEN_HOUSE/pegs4/density_corrections/elements/carbon_graphite_ICRU90_1.7g_cm3.density 
$HEN_HOUSE/pegs4/density_corrections/elements/carbon_graphite_ICRU90_2.0g_cm3.density
$HEN_HOUSE/pegs4/density_corrections/compounds/water_ICRU90.density

The easiest way is to take advantage of the pegsless option which generates the data on-the-fly by adding the following input block to the input file:

:start media definition:
    ### energy transport cutoffs
    ae = 0.512
    ue = 2.511
    ap = 0.001
    up = 2.00
    :start water_icru90:
        density correction file     = water_icru90
        bremsstrahlung correction   = NRC
    :stop water_icru90:
    :start c170_icru90:
      density correction file     =  carbon_graphite_ICRU90_1.7g_cm3
        bremsstrahlung correction   = NRC
    :stop c170_icru90:
    :start c200_icru90:
      density correction file     =  carbon_graphite_ICRU90_2.0g_cm3
        bremsstrahlung correction   = NRC
    :stop c200_icru90:
:stop media definition:

One can of course create a pegs4 data file for these media making sure to use the proper density correction files.

Several checks have been made to ensure that the implementation produces correct values. The graph below shows the ratio of air mass-energy absorption coefficients calculated with the EGSnrc app g.mortran either with the original XCOM PE cross sections or with the re-normalized cross sections of Sabbatucci and Salvat to the Penelope values taken from the ICRU90 report. The first shows similar differences as reported in ICRU90 report while the second demonstrates excellent agreement between both codes except at 3.263 keV (Ar K shell binding energy in Penelope) because of the small differences in binding energies between the codes (EGSnrc uses the XCOM value of 3.2029 keV).

To make use of these features one can clone the git repository for EGSnrc and then all branches are available to the user. To switch to the branch feature-shellwise-pe-xsections one needs to execute

git checkout feature-shellwise-pe-xsections

on the command line in the folder containing the git repository. Once on that branch the user can create a configuration either for Linux or Windows.

Alternatively to using git, one can download this branch directly as a zip archive and then configure it as a separate install using one of the methods mentioned in the previous paragraph. Go to the EGSnrc github page of this branch and download it by clicking the green button on the right Clone or download and click on the Download ZIP link.

Details to take into consideration when using this branch and dosxyznrc:

If one only chooses to create a new configuration for this branch, the user area $EGS_HOME remains untouched.
As a consequence dosxyznrc.mortran is not updated, hence its compilation will fail during the execution of EGSnrc-configure-linux as this file was modified for this branch.
The solution is to update dosxyznrc.mortran from $HEN_HOUSE/user_codes/dosxyznrc/
Setup the environment variable $EGS_CONFIG in .bashrc to point to the new configuration and start a new shell console.

rtownson

Great work Ernesto! We should make obvious note somewhere of all the changes to default parameters.

rtownson · 2017-08-16T14:21:00Z

HEN_HOUSE/omega/progs/gui/beamnrc/beamnrc_params.tcl

 are relaxed via emission of fluorescent X-Rays,\
 Auger and Koster-Cronig electrons.\
 Make sure to turn this option on for low energy\
 applications.\

+Selecting On or EADl is equivalent. A full atomic relaxation cascade\
+is simmulated using EADL transition probabilities.


Typo in simulated

rtownson · 2017-08-16T14:22:27Z

HEN_HOUSE/omega/progs/gui/dosxyznrc/xyznrc_parameters.tcl

+ applications.
+
+Selecting On or EADl is equivalent. A full atomic relaxation cascade\
+is simmulated using EADL transition probabilities.


rtownson · 2017-08-16T14:23:56Z

HEN_HOUSE/src/egs_utilities.mortran

+ of the current simlation. This is due to the fact that energies 
+ below 1 keV are taken from the relaxation database only for 
+ elements requested in the input when EPDL is used. This is not
+ and issue when using XCOM.


Typo: and should be an

rtownson · 2017-08-16T14:30:42Z

HEN_HOUSE/src/get_inputs.mortran

@@ -1690,6 +1713,28 @@ write(i_errors,*)
 ' ******************** end input transport parameter *********************** ';
 write(i_errors,*);

+"Check if EADL relaxation requested. Note that original relaxation"
+"algorithm using <M> and <N> is only turned ON for all regions."
+"Moved passed the $TURN-ON/OFF-IN-REGIONS statement to catch the"


Typo: passed should be past

Thanks a bunch @rtownson for reviewing! Somehow I thought you had made the corrections. Will fix them right away! Updating the documentation is the next step.

blakewalters

Very nice, Ernesto!

ftessier · 2017-08-16T16:07:54Z

@mainegra , do you want to fix the noted typos? You can add a commit at the end; or else if you prefer you can use git rebase -i af4a48a and change pick to edit for the commits where the typos appeared. Then git replays the commits and stop for the ones you wanted to modify. Edit the files, git add, and keep going with git rebase --continue. Loads of fun!

mainegra · 2017-08-16T16:08:45Z

Thanks @blakewalters for reviewing!

mainegra · 2017-08-16T16:10:32Z

@ftessier I will correct the typos found by @rtownson ! Seems easier! :-)

mainegra · 2017-08-16T16:16:29Z

Done!

ftessier · 2017-08-16T16:22:22Z

@mainegra can you rebase origin/develop? There is a conflict which I am not sure how to resolve, viz. should the explanation in the comment remain? I presume not, just want to make sure. Can we remove the commented out IF block?

<<<<<<< f0c54d1ebe26c5a1acd464440846f458990bc101
    "Assign this energy deposition to an electron."
    "Note that this should NOT be treated as though it came from a photon,"
    "even if a photon initiated the relaxations. Rather, energy remaining"
    "in vacancies should be given to an electron or absorbed locally."
    $AUSCALL($SELECTRONA); "call ausgab with iarg=34"
||||||| merged common ancestors
    IF( relax_initiator_q = 0 ) [
        "particle initiating relaxations was photon"
        $AUSCALL($SPHOTONA);  "call ausgab with iarg=33"
    ] ELSE [
        "particle initiating relaxations was electron"
        $AUSCALL($SELECTRONA); "call ausgab with iarg=34"
    ]
=======
    /*IF( relax_initiator_q = 0 ) [
        "particle initiating relaxations was photon"
        $AUSCALL($SPHOTONA);  "call ausgab with iarg=33"
    ] ELSE [
        "particle initiating relaxations was electron"
        $AUSCALL($SELECTRONA); "call ausgab with iarg=34"
    ]*/
    $AUSCALL($SELECTRONA); "call ausgab with iarg=34"
>>>>>>> Implement shellwise photoelectric interaction

ftessier · 2017-08-16T19:07:53Z

Rebased on tip of develop to fix the merge conflict.

rtownson · 2017-08-16T19:56:03Z

Now that you've rebased you can remove the commented out block Fred referenced. relax_initiator_q is no longer used, just the $AUSCALL($SELECTRONA); line should remain.

I also see another instance in egsnrc.mortran where the relax_initiator_q line should be removed:

IF( do_relax ) [
    relax_initiator_q = 0; "a photon is initiating relaxations"
    call egs_eadl_relax(iZ,k);
]

ftessier · 2017-08-16T20:13:37Z

Oh, I forgot to commit the merge conflict resolution.

ftessier · 2017-08-17T19:26:57Z

Updated commit history (reorder, reword, edit, squash, etc.), adjusted contributors in egs_functions.cpp and egs_functions.h, removed EOL white space and ran astyle in egs++ directory. The few code changes can be reviewed with git diff a5aa2b1.

ftessier

The variable relax_initiator_q was used on this branch, but in the mean time it had been deleted on the develop branch.

Add photoelectric (PE) cross sections by Sabbatucci and Salvat based on Dirac–Hartree–Fock–Slater (DHFS) calculations, with a normalization screening correction defined as the ratio of electron densities obtained with DHFS and multi-configuration Dirac-Fock (MCDF) methods. These cross sections are expected to closely reproduce MCDF cross sections. The data base is a 1-byte record length binary file with total and shell-wise photoelectric cross sections. The data is arranged sequentially so as to allow direct access to any given element: NZ fpos(1)...fpos(NZ) NE(1),Nshell(1) E(1,1) SigTot Sig(K) Sig(L1)...Sig(Nshell(1)) ... E(1,NE(1)) SigTot Sig(K) Sig(L1)...Sig(Nshell(1)) ... E(NZ,1) SigTot Sig(K) Sig(L1)...Sig(Nshell(1)) ... E(NZ,NE(1)) SigTot Sig(K) Sig(L1)...Sig(Nshell(1)) ...

- Turn on the use of renormalized photoelectric (PE) cross sections by setting the photon cross sections key to either 'mcdf-xcom' or 'mcdf-epdl': a) Sets variable mcdf_pe_xsections to true. b) Sets photon cross section data base to xcom or epdl, respectively. c) If relaxation is requested, only EADL relaxation is allowed and shell-wise cross sections are used to sample the interacting shell. d) Subroutine egs_shellwise_photo() is called from PHOTO. MCDF stands for "multi-configuration Dirac-Fock" because the normalization screening correction used by Sabbatucci and Salvat brings their DHFS cross sections closer to MCDF results. - If mcdf_pe_xsections is true: a) egs_read_shellwise_pe reads in MCDF shell-wise PE cross sections by Sabbatucci and Salvat only for the elements in the problem, sorted by increasing Z value. b) egs_read_shellwise_pe adjusts cross sections to the EGSnrc binding energies via an energy scale shift. c) egs_read_shellwise_pe normalizes elemental shell cross sections to the total cross section of an element preparing a log(E) log(sigma(E,Z,shell)/sigma(E,Z,total)) table for interpolation when sampling the interacting shell. d) egsi_get_shell_data prepares a database of medium PE cross sections for sampling the interaction using the same energy grid as the cross sections for the other interactions, and normalizes elemental PE cross sections to the medium PE cross section for sampling the element with which the photon interacts, maintaining Sabbatucci and Salvat's energy grid. pe_elem_prob(j,k,medium) = log(pz*sigma(E,Z)/sigma(E,medium)) - egs_shellwise_photo is the new photoelectric routine: a) Determines the energy interval using binary search to preserve structure near edges (one search per element). b) Samples the interacting element by looping through the medium's elements, starting with the highest Z element using pre-calculated pe_elem_prob(j,k,medium) in egsi_get_shell_data. c) The interacting shell is sampled using normalized shell-wise cross sections pe_xsection(j,i,l) for shells with binding energy above a minimum energy, set currently to 1 keV. d) In the case of no atomic relaxations, a photoelectron is simply generated. - Replace macro $EADL_RELAX with integer variable eadl_relax. Acceptable values are: Off -> no atomic relaxation, just generate photoelectron On or eadl -> Use detailed relaxation with EADL transition prob. simple -> Use original relaxation algorithm with <M> and <N> - Move $RELAX-CUTOFF to egsnrc.macros; currently set to 1 keV. $RELAX-CUTOFF defines the lowest energy for which to simulate the photoelectric interaction and create initial vacancies. Note that relaxation transitions can produce sub-threshold particles and vacancies. - Set minimum binding energy requiring relaxation data to $RELAX-CUTOFF rather than the minimum of AE or AP. - In egs_eadl_relax: initial vacancies below L3 are only allowed if using MCDF shell-wise PE cross sections. Note that in this case COMPT and PHOTO will create such vacancies, but not the electron impact ionization (EII) routine eii_sample. This should be expanded when implementing the EII cross sections from Bote and Salvat.

Add an option to select a detailed EADL relaxation and photoelectric (PE) interaction using renormalized cross sections via the variables eadl_relax and mcdf_pe_xsections respectively. The key value pair is: Photon cross section = mcdf-xcom # or mcdf-epdl Note that the only cross sections taken from either the xcom or epdl series are Rayleigh, pair and triplet cross sections. When renormalized PE cross sections are used, only EADL relaxation is allowed, with a warning issued when changing from simple to detailed relaxation. Turning OFF relaxation is also possible.

Force the simulation of atomic relaxations for all interactions with atomic electrons. Before, this was enforced automatically for electron impact ionization (EII) and bound Compton scattering, but not for the photoelectric effect. For consistency, atomic relaxations should be enabled (or disabled) for all atomic interactions.

Use the C++ string copy function instead of strcpy. The latter adds a terminating null character to the character array which interferes with contiguous character arrays in COMMON blocks or C data structures, such as the_media->pxsec and the_media->eiixsec.

Incoherent scattering by atomic electrons is ON by default, hence coherent (Rayleigh) scattering should also be ON by default, for a physically consistent photon scatter model. If atomic relaxations are important, then electron impact ionization should also be turned ON. Currently EII is OFF by default. In cases where scattering by atomic electrons has a negligible effect, such as in MV beam calculations, one can turn OFF coherent scattering and use Klein-Nishina for incoherent scattering. For an ion chamber dose calculation at the reference depth in a water phantom, the time penalty of turning ON photon interactions with atomic electrons is about 6% for a Co-60 beam and 4% for a 6 MV beam.

DOSXYZnrc overrides some of the egsnrc.macros defaults. It explicitly turns OFF Compton binding corrections and atomic relaxations. For consistency, coherent (Rayleigh) scattering is now also turned OFF by default since it doesn't make sense physically to have incoherent scattering off a free electron while including coherent scattering off atomic electrons.

Rename the binary search function in DOSXYZnrc to avoid conflict with the new binary search function in egs_utilities.mortran (needed by the new photoelectric routine). The function name is changed from ibsearch to ibsearch4, reflecting the fact that it takes a real*4 array. Remove the binary search function ibsearch from the EDKnrc application source code. This function is now defined in egs_utilities.mortran, and in fact it is not even used in EDKnrc.

Add newline characters in the message warning the user that the relaxation model is changed to the full EADL algorithm (triggered when the user selects the renormalized photoelectric cross sections but the default relaxation algorithm).

Add gui options to use renormalized photoelectric (PE) cross sections and to select between EADL and the original relaxation algorithm. Also set Rayleigh scattering ON by default. Update help text accordingly.

ICRU Report 90 recommends that density-effect corrections for graphite be calculated for a mean excitation energy I = 81 eV and the crystalline density rho = 2.265 g/cm3. Using the Fortran version of the ESTAR program two files are created, for graphite bulk mass densities of 1.7 g/cm3 and 2.0 g/cm3 (both files contain the same density-effect correction data, calculated for 2.265 g/cm3).

Add a density-effect correction file for water, for a mean excitation energy I = 78 eV, following the recommendation of ICRU Report 90. The data is obtained from the Fortran version of the ESTAR program.

Move the default eadl_relax block below the call to the macro $TURN-ON/OFF-IN-REGIONS which sets the default for atomic relaxations. Otherwise there is a mismatch between eadl_relax and iedgfl. This problem was revealed when testing dosxyznrc.mortran, which overrides egsnrc.macros defaults.

Monte Carlo (MC) transport parameters allowing inputs beyond ON or OFF by region must shift the output_strings array index by 2 to get the proper output when the default values are used, i.e., no input is provided.

Inform the user about the requirement to turn on bound Compton when selecting EADL relaxation, since binding energies below 1 keV are taken from the incoh.data file.

ftessier · 2017-08-18T13:50:29Z

Just squashed the EDKnrc update for ibsearch into commit 26156f5.

mainegra requested review from ftessier, rtownson and blakewalters August 8, 2017 20:16

mainegra added the improvement label Aug 9, 2017

ftessier assigned mainegra Aug 10, 2017

mainegra mentioned this pull request Aug 11, 2017

Atomic relaxation method should be an input option #310

Closed

rtownson approved these changes Aug 16, 2017

View reviewed changes

blakewalters approved these changes Aug 16, 2017

View reviewed changes

ftessier force-pushed the feature-shellwise-pe-xsections branch from 4668d13 to 72f6e9b Compare August 16, 2017 19:07

ftessier force-pushed the feature-shellwise-pe-xsections branch 2 times, most recently from a5aa2b1 to 7e4ba64 Compare August 17, 2017 19:19

ftessier reviewed Aug 17, 2017

View reviewed changes

mainegra added 9 commits August 18, 2017 09:44

Add missing newline characters in warning message

debd967

Add newline characters in the message warning the user that the relaxation model is changed to the full EADL algorithm (triggered when the user selects the renormalized photoelectric cross sections but the default relaxation algorithm).

mainegra and others added 9 commits August 18, 2017 09:48

Update PE and atomic relaxation options in guis

6c7821a

Add gui options to use renormalized photoelectric (PE) cross sections and to select between EADL and the original relaxation algorithm. Also set Rayleigh scattering ON by default. Update help text accordingly.

Add ICRU90 water density-effect correction data

0d0e06d

Add a density-effect correction file for water, for a mean excitation energy I = 78 eV, following the recommendation of ICRU Report 90. The data is obtained from the Fortran version of the ESTAR program.

Set default MC transport parameter output string

3348366

Monte Carlo (MC) transport parameters allowing inputs beyond ON or OFF by region must shift the output_strings array index by 2 to get the proper output when the default values are used, i.e., no input is provided.

Format atomic relaxation error messages

bd11f03

Inform the user about the requirement to turn on bound Compton when selecting EADL relaxation, since binding energies below 1 keV are taken from the incoh.data file.

Fix typos in help text and comments

dd2156f

Run astyle to apply EGSnrc coding style standard

3e0c0cb

Remove legacy variable relax_initiator_q

82e6611

ftessier force-pushed the feature-shellwise-pe-xsections branch from 08cd89b to 82e6611 Compare August 18, 2017 13:49

ftessier approved these changes Aug 18, 2017

View reviewed changes

ftessier merged commit fca7801 into develop Aug 18, 2017

ftessier deleted the feature-shellwise-pe-xsections branch August 18, 2017 15:27

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Feature renormalized photoelectric cross sections and ICRU90 recommended values for key data #322

Feature renormalized photoelectric cross sections and ICRU90 recommended values for key data #322

mainegra commented Aug 8, 2017 •

edited

Loading

rtownson left a comment

rtownson Aug 16, 2017

rtownson Aug 16, 2017

rtownson Aug 16, 2017

rtownson Aug 16, 2017

mainegra Aug 16, 2017 •

edited

Loading

blakewalters left a comment

ftessier commented Aug 16, 2017 •

edited

Loading

mainegra commented Aug 16, 2017

mainegra commented Aug 16, 2017

mainegra commented Aug 16, 2017

ftessier commented Aug 16, 2017 •

edited

Loading

ftessier commented Aug 16, 2017

rtownson commented Aug 16, 2017 •

edited

Loading

ftessier commented Aug 16, 2017 •

edited

Loading

ftessier commented Aug 17, 2017

ftessier left a comment

ftessier commented Aug 18, 2017

Feature renormalized photoelectric cross sections and ICRU90 recommended values for key data #322

Feature renormalized photoelectric cross sections and ICRU90 recommended values for key data #322

Conversation

mainegra commented Aug 8, 2017 • edited Loading

rtownson left a comment

Choose a reason for hiding this comment

rtownson Aug 16, 2017

Choose a reason for hiding this comment

rtownson Aug 16, 2017

Choose a reason for hiding this comment

rtownson Aug 16, 2017

Choose a reason for hiding this comment

rtownson Aug 16, 2017

Choose a reason for hiding this comment

mainegra Aug 16, 2017 • edited Loading

Choose a reason for hiding this comment

blakewalters left a comment

Choose a reason for hiding this comment

ftessier commented Aug 16, 2017 • edited Loading

mainegra commented Aug 16, 2017

mainegra commented Aug 16, 2017

mainegra commented Aug 16, 2017

ftessier commented Aug 16, 2017 • edited Loading

ftessier commented Aug 16, 2017

rtownson commented Aug 16, 2017 • edited Loading

ftessier commented Aug 16, 2017 • edited Loading

ftessier commented Aug 17, 2017

ftessier left a comment

Choose a reason for hiding this comment

ftessier commented Aug 18, 2017

mainegra commented Aug 8, 2017 •

edited

Loading

mainegra Aug 16, 2017 •

edited

Loading

ftessier commented Aug 16, 2017 •

edited

Loading

ftessier commented Aug 16, 2017 •

edited

Loading

rtownson commented Aug 16, 2017 •

edited

Loading

ftessier commented Aug 16, 2017 •

edited

Loading