Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix g calc type 1 #1228

Open
wants to merge 4 commits into
base: develop
Choose a base branch
from
Open

Fix g calc type 1 #1228

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
353c994
Harmonize calc types in g
mainegra Oct 14, 2024
a922f4a
Fix bug in command line arguments parser
mainegra Nov 28, 2024
a555c70
Update documentation
mainegra Nov 28, 2024
2203c6d
Move NIST values section up
mainegra Nov 28, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 28 additions & 6 deletions HEN_HOUSE/doc/src/pirs3100-g-app/g.doxy
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,8 +67,9 @@ manual created.

In the current version (1.6) the original algorithm of a `calculation type` 1 has been slightly adjusted to ensure @muen_rho is calculated to a user-requested relative precision @accu and it has been made independent of the number of histories. The target statistical precision now defaults to 1%. The possibility to further optimize `calculation type` 1 in the megavoltage energy range is introduced. The calculation progress messages during a type 1 calculation were made clearer and the documentation has been updated to reflect these changes.

The refinements introduced in 2017 to correctly classify subthreshold events were initially implemented only for `calculation type` 0. As of October 2024, these refinements have also been incorporated into `calculation type` 1 for consistency (see git `commit b9627ee3`). Notably, the NIST \muen_rho values align only when this bug fix is excluded, suggesting that the NIST database may require an update to reflect these corrections.

\section g_method Calculation method
\section g_method Methods

A description of the physical quantities and the equations used for their calculation is presented in this section. Most of these definitions can be found in basic radiation physics textbooks, hence users more interested with practical aspects of using the code can skip to the next section \ref g_calc_types.

Expand Down Expand Up @@ -251,11 +252,15 @@ is scored as part of \Etr.

\section g_calc_types Calculation types

\subsection g_calc_type_0 Calculation type 0

In its default implementation, the \c g code loops over a fixed number of histories
for a source of \a photons or <em>charged particles</em> sampling the initial particle parameters
and invoking the electromagnetic shower routine. This calculation can be explicitly requested by setting the `calculation type` input to 0. Collision kerma is calculated using Eq.(\reflh{eq_kerma_0,10})
and the energy fluence averaged mass energy absorption coefficient \muen_rho_ave is calculated using Eq.(\reflh{eq_muen_0,15})

\subsection g_calc_type_1 Calculation type 1

A `calculation type` 1, for use with photon sources only, is tailored to the efficient calculation of the mass energy absorption coefficient \muen_rho_ave. In this mode the `g` application takes advantage of the faster convergence of \mutr and 1-\g (when \g is small) to a desired statistical precision and collision kerma is calculated with Eq.(\reflh{eq_kerma_1,11})
and the energy fluence averaged mass energy absorption coefficient \muen_rho_ave is calculated using Eq.(\reflh{eq_muen_1,16}).
The calculation finishes as soon as a user-requested relative uncertainty @accu is reached, which by default is set to 0.01 (1%). For this calculation type the `number of histories` input entry is ignored.
Expand All @@ -281,7 +286,7 @@ be placed between delimiters as is customary in most EGSnrc applications. An exa
\latexonly
\\
\endlatexonly
\code
```
####################################################################
# calculation types
#
Expand All @@ -297,6 +302,23 @@ be placed between delimiters as is customary in most EGSnrc applications. An exa
######################
precision = 0.001 # relative precision
#############
```

\subsection g_nist Reproducing NIST values

The refinements introduced in 2017, which account for both radiative and non-radiative subthreshold events, result in discrepancies with the current NIST XCOM \muen_rho values. This suggests that the NIST data may have been derived without explicitly distinguishing between radiative and non-radiative subthreshold events, despite their documentation implying otherwise. Supporting this interpretation is the observation that values calculated without these modifications to the g application align closely with the NIST data.

To request a \muen_rho calculation matching NIST values the user should use the 'CLASSIFY SUB-THRESHOLD RELAXATIONS' input key (case-insensitive):
\latexonly
\\
\endlatexonly
\code
#####
# Option to reproduce NIST values
# Defaults to 'yes'
#####
CLASSIFY SUB-THRESHOLD RELAXATIONS = NO
...
\endcode


Expand All @@ -307,7 +329,7 @@ delimiters. For instance, to define a monoenergetic source emitting particles of
and charge \f$q\f$, the input keys required are:

\verbatim
:Start Source Input:
:Start Source Input:
Incident Charge = q # 0 for photons -1 or 1 for e- or e+
INCIDENT SPECTRUM = mono-energy
INCIDENT KINETIC ENERGY = E # in MeV
Expand All @@ -318,7 +340,7 @@ Several individual energies can be requested by adding comma separated
energy values:

\verbatim
:Start Source Input:
:Start Source Input:
Incident Charge = q # 0 for photons -1 or 1 for e- or e+
INCIDENT SPECTRUM = mono-energy
INCIDENT KINETIC ENERGY = E1,...,EN # in MeV
Expand All @@ -331,7 +353,7 @@ To generate a database of \mutr or \muen values as function of energy, one can d
a logarithmic energy grid. A linear (equidistant) energy grid is defined using the input keys

\verbatim
:Start Source Input:
:Start Source Input:
Incident Charge = q # 0 for photons -1 or 1 for e- or e+
INCIDENT SPECTRUM = mono-energy-lin-range
INCIDENT KINETIC ENERGY = Emin, Emax, deltaE
Expand All @@ -344,7 +366,7 @@ where the input key `deltaE` defines the bin width \f$\Delta\mathrm{E}\f$ defini
Similarly, a logarithmic energy grid is defined using the input keys

\verbatim
:Start Source Input:
:Start Source Input:

INCIDENT SPECTRUM= mono-energy-log-range
INCIDENT KINETIC ENERGY= Emin, Emax, N
Expand Down
2 changes: 1 addition & 1 deletion HEN_HOUSE/doc/src/pirs3100-g-app/makedoc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ while test $# -gt 0; do
dir=*) docdir=`echo $1 | sed 's/dir=//'` ;;
eh=*) eh=`echo $1 | sed 's/eh=//'` ;;
eq=*) update_eq=`echo $1 | sed 's/eq=//'` ;;
latex=*) latex_build=`echo $1 | sed 's/eq=//'` ;;
latex=*) latex_build=`echo $1 | sed 's/latex=//'` ;;
config=*) my_machine=`echo $1 | sed 's/config=//'` ;;
esac
shift
Expand Down
82 changes: 66 additions & 16 deletions HEN_HOUSE/user_codes/g/g.mortran
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -183,13 +183,10 @@ C> @param[out] e_brem/e_tot the average fraction of the kinetic energy
C> subsequently lost in bremsstrahlung (photon-emitting) events
C> @param[out] e_rad/e_tot same as g above
C> @param anorm
C> @param npgi
C> @param npei
C> @param E_ave Average spectrum energy
C> @param factor Converts e_mutr and e_muen scored as
C> \f$\mathrm{MeV}\,\mathrm{cm}^2/\mathrm{g}\f$
C> to \f$\mathrm{Gy}\,\mathrm{cm}^2\f$
C> @param de_pulsei
C> @param mutr The correct definition for a spectrum is that given by Attix
C> of an energy fluence weighted average mass transfer coefficient
C> \f[
Expand Down Expand Up @@ -255,14 +252,16 @@ REPLACE {$CALL-HOWFAR-IN-PHOTON;} WITH {;
REPLACE {;COMIN/SCORE/;} WITH {;
common/score/ e_tot,e_tot2,e_brem,e_brem2,e_bremc,e_rad,e_rad2,e_radc,
etot_tmp,ebrem_tmp,erad_tmp,my_gle,accu,ncase,calc_type,
during_pe_compt,during_eii,verbose,m_balance;
during_pe_compt,during_eii,verbose,m_balance,
classify_sub_threshold_relaxations;
"It is always a good idea to score in double precision!"
real*8 e_tot,e_tot2,e_brem,e_brem2,e_rad,e_rad2,
etot_tmp,ebrem_tmp,erad_tmp,e_radc,e_bremc;
$LONG_INT ncase;
$REAL my_gle,accu;
$INTEGER calc_type,during_pe_compt,during_eii,m_balance;
$LOGICAL verbose;
$LOGICAL classify_sub_threshold_relaxations;
};

REPLACE {$MXENE} WITH {2000}; "maximum number of energies to be looped over"
Expand Down Expand Up @@ -323,16 +322,16 @@ real*8 ert, ert2, ett, ett2, g, g_rel_err;
real*8 accu0, m;
$REAL weight,weight2,cov,covt,rel_e;
$LONG_INT nperbatch,ncasei,icase,nmutr,nmuen;
$REAL de_pulsei, Eave, err_Eave,factor,gmfp;
$REAL Eave, err_Eave,factor,gmfp;
$REAL t_mutr, t_muen;
"For interacting with the source routine and shower"
$REAL ein,uin,vin,win,wtin,xin,yin,zin,ecut_ask,pcut_ask;
$REAL gbr1,gbr2,rnno,err_frac;
$INTEGER iqin,irin,ip;
$REAL gbr1, gbr2, rnno, err_frac;
$INTEGER iqin, irin, ip, iarg;
$INTEGER itimes, ntimes;"loop through different energies"
real*4 cpu,etime,time_array(2);
$INTEGER datcount;
$INTEGER npgi,npei,nspliti,nbini,lmy_gle;
$INTEGER lmy_gle;
$INTEGER nbatch, "$INTEGER and $REAL are defined in egsnrc.macros"
ibatch, "they can be used (if employed consistently throughout"
i,j,i_log; "the user code) e.g. to switch to double precision"
Expand Down Expand Up @@ -540,15 +539,33 @@ IF( calc_type = 1 ) [ "will seek a given precision -- default does not do this"
$EVALUATE gbr1 USING gbr1(my_gle);
$EVALUATE gbr2 USING gbr2(my_gle);
np = 1; e(np) = ein; u(np)=0; v(np)=0; w(np)=1; iq(np) = 0;
$RANDOMSET rnno; edep = 0;
$RANDOMSET rnno; edep = 0; etot_tmp = 0.0;
IF( rnno < gbr1 & ein > rmt2 )[
call pair;
IF (classify_sub_threshold_relaxations)[
$AUSCALL($PAIRAUSB);
]
call pair;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reduce indentation for call pair;

IF (classify_sub_threshold_relaxations)[
$AUSCALL($PAIRAUSA);
]
]
ELSE IF( rnno < gbr2 ) [
IF (classify_sub_threshold_relaxations)[
$AUSCALL($COMPAUSB);
]
call compt;
IF (classify_sub_threshold_relaxations)[
$AUSCALL($COMPAUSA);
]
]
ELSE [
IF (classify_sub_threshold_relaxations)[
$AUSCALL($PHOTOAUSB);
]
call photo;
IF (classify_sub_threshold_relaxations)[
$AUSCALL($PHOTOAUSA);
]
]
/****************************************************
Collects photons' kinetic energy transferred to e-
Expand All @@ -558,12 +575,19 @@ IF( calc_type = 1 ) [ "will seek a given precision -- default does not do this"
transfer is of interest in this block, there is no need
for further classification.
*******************************************************/
DO ip=1,np [
IF( iq(ip) ~= 0 ) edep = edep + e(ip)-prm;
IF (~ classify_sub_threshold_relaxations)[
DO ip=1,np [
IF( iq(ip) ~= 0 ) edep = edep + e(ip)-prm;
]
]
nmutr = nmutr + 1;
IF(gmfp > 0) [
aux = edep/gmfp/rho(1);
IF (classify_sub_threshold_relaxations)[
aux = etot_tmp/gmfp/rho(1);
]
ELSE[
aux = (edep)/gmfp/rho(1);
]
] ELSE [
aux = 0;
]
Expand Down Expand Up @@ -887,7 +911,7 @@ DO icase=1,ncase [
erad_tmp = 0.0; "loses just to annihilation"
"initialize etot_tmp differently for photons and electrons"
IF( iqin .ne. 0 ) [ etot_tmp = (ein - rm); ]
ELSE [ etot_tmp = 0; ]
ELSE [ etot_tmp = 0; ]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is best to avoid if-else one-liner statements.

"============================================================================"

call shower(iqin,ein,xin,yin,zin,uin,vin,win,irin,wtin);
Expand Down Expand Up @@ -1096,6 +1120,7 @@ accu = 0.01;
m_balance = 1;
calc_type = 0;"Defaults to allow photon/electron calculations"
verbose = .false.;
classify_sub_threshold_relaxations = .true.;
"print a summary on unit 6"
ival = 0;
delimeter = 'NONE';
Expand Down Expand Up @@ -1163,7 +1188,7 @@ IF( error_flag = 0 ) [
ELSE [ error_flag = 0; ]

"----------------------------------------------------------------------"
" precision sought and calculation type "
" precision sought and calculation type "
"----------------------------------------------------------------------"
ival = ival + 1;
values_sought(ival) = 'PRECISION';
Expand Down Expand Up @@ -1203,6 +1228,31 @@ $GET_INPUT(IVAL);
IF( error_flag = 0 ) [ calc_type = value(ival,1); ]
ELSE [ error_flag = 0; ]

"-------------------------------------------------"
" Ignore sub-threshold relaxations? "
"-------------------------------------------------"
ival = ival + 1;
values_sought(ival) = 'CLASSIFY SUB-THRESHOLD RELAXATIONS';
nvalue(ival) = 1;
type(ival) = 3;
allowed_inputs(ival,0) = 'NO';
allowed_inputs(ival,1) = 'YES';
$GET_INPUT(IVAL);
IF( error_flag = 0 ) [
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ideally, we would harmonize the IF bracketing style. My perspective is that we ought to settle on conditional branch style with spaces around (...) (no spaces inside the brackets, and one space before the opening [

IF (...) [
  ...
]
ELSE (...) [
  ...
]

IF( value(ival,1) = 0 )[
classify_sub_threshold_relaxations = .false.;
write(*,*);
write(*,*) '==> Sub-threshold relaxations during Compton';
write(*,*) ' and photo-electric events NOT classified!';
write(*,*) ' (Matches NIST values)';
write(*,*);
]
ELSE[
classify_sub_threshold_relaxations = .true.;
]
]
ELSE [ error_flag = 0; ]

"-------------------------------------------------"
" verbose output? "
"-------------------------------------------------"
Expand All @@ -1218,7 +1268,7 @@ IF( error_flag = 0 ) [
verbose = .true.;
]
]
ELSE [ error_flag = 0; ]
ELSE [ error_flag = 0; ]

"---------------------------------------------------------------------"
"STEP 6 DETERMINATION-OF-INICIDENT-PARTICLE-PARAMETERS "
Expand Down