Skip to content

Commit

Permalink
Fix atomic relaxation vacancy dependency on ecut
Browse files Browse the repository at this point in the history 
Check only against the RELAX-CUTOFF variable for energy cutoff on
vacancy creation in both relaxation routines, in order to make the
calculation of g independent of the cutoff energy ecut. Also fix a bug
in the output of the g application.

Previously, the energy cutoff for vacancy creation was the maximum of
ecut and pcut, or the RELAX-CUTOFF. Thus changing ecut caused a
different chain of relaxations, hence results related to sub-threshold
contributions (e.g., the value of g) were ecut dependent.
  • Loading branch information
@rtownson @ftessier
rtownson authored and ftessier committed May 25, 2017
1 parent babf4de commit 6d487e7
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 59 deletions.
96 changes: 39 additions & 57 deletions HEN_HOUSE/src/egsnrc.mortran
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5212,7 +5212,7 @@ $REAL e_array($MXVAC), "array with vacancy energies "
Ex, "kinetic energy of emitted particle "
eta, "a random number "
e_check, "energy conservation check "
e_cut,ekcut,pkcut,elcut,max_Ec,max_Pc; "cut-off energies "
min_E,ekcut,pkcut,elcut; "cut-off energies "

$REAL xphi,yphi,xphi2,yphi2,rhophi2,
cphi,sphi; "for azimuthal angle selection"
Expand Down Expand Up @@ -5256,31 +5256,22 @@ IF( n < 1 | n > $MXSHELL ) [ return; ] "unknown vacancy"
iz_relax = iZ;
irl = ir(np);
ekcut = ecut(irl)-rm; pkcut = pcut(irl);
max_Ec = Max($RELAX-CUTOFF,ekcut);
max_Pc = Max($RELAX-CUTOFF,pkcut);
e_cut = Min(ekcut,pkcut);
e_cut = Max($RELAX-CUTOFF,e_cut);

IF( relax_initiator_q = 0 ) [ "particle initiating relaxations was photon"
IF( energy <= max_Pc ) [
edep = edep + energy; "We assume that edep is zeroed "
"(or set to the appropriate value in the routine "
"calling RELAX "

edep_local = energy;
$AUSCALL($SPHOTONA);
return;
]
] ELSE [ "particle initiating relaxations was electron"
IF( energy <= max_Ec ) [
edep = edep + energy; "We assume that edep is zeroed "
"(or set to the appropriate value in the routine "
"calling RELAX "
min_E = $RELAX-CUTOFF;

IF( energy <= min_E ) [
edep = edep + energy; "We assume that edep is zeroed "
"(or set to the appropriate value) in the routine "
"calling RELAX "
edep_local = energy;

edep_local = energy;
"particle initiating relaxations was photon"
IF( relax_initiator_q = 0 ) [
$AUSCALL($SPHOTONA);
] ELSE [ "particle initiating relaxations was electron"
$AUSCALL($SELECTRONA);
return;
]

return;
]

" Set-up the array of vacancies for the relaxation cascade "
Expand All @@ -5291,26 +5282,20 @@ e_check = 0; e_array(n_vac) = energy;

shell = vac_array(n_vac); Ei = e_array(n_vac); n_vac = n_vac - 1;

IF( relax_initiator_q = 0 ) [ "particle initiating relaxations was photon"
IF( Ei <= max_Pc ) [ " Below cut-off -> local absorption "
edep = edep + Ei;

edep_local = Ei;
$AUSCALL($SPHOTONA);
IF( Ei <= min_E ) [ " Below cut-off -> local absorption "
edep = edep + Ei;

IF( n_vac > 0 ) goto :START: ;
EXIT;
]
] ELSE [ "particle initiating relaxations was electron"
IF( Ei <= max_Ec ) [ " Below cut-off -> local absorption "
edep = edep + Ei;
edep_local = Ei;

edep_local = Ei;
"particle initiating relaxations was photon"
IF( relax_initiator_q = 0 ) [
$AUSCALL($SPHOTONA);
] ELSE [ "particle initiating relaxations was electron"
$AUSCALL($SELECTRONA);

IF( n_vac > 0 ) goto :START: ;
EXIT;
]

IF( n_vac > 0 ) goto :START: ;
EXIT;
]

"Set the relax_user common block variables, IK March 22 2004"
Expand Down Expand Up @@ -5670,12 +5655,11 @@ $declare_max_medium;
$COMIN-RELAX-EADL;
;COMIN/RELAX-USER/;

$REAL Ec,Pc,max_Ec,max_Pc,min_E,rnno,Evac,Ef,Ef1,Ef2,Ex,Ecc,
$REAL Ec,Pc,min_E,rnno,Evac,Ef,Ef1,Ef2,Ex,Ecc,
cost,sint,cphi,sphi;
"Ec ecut as k.e. for current region"
"Pc pcut for current region"
"max_Ec maximum of Ec and $RELAX-CUTOFF"
"max_Pc maximum of Pc and $RELAX-CUTOFF"
"min_E the minimum energy for transitions to continue"
"rnno a random number"
"Evac binding energy of current vacancy"
"Ef the sum of the binding energies of the new vacancies"
Expand Down Expand Up @@ -5710,9 +5694,8 @@ IF( shell < 1 | shell > $MAXSHELL ) [ return; ] "unknown vacancy"
irl = ir(np);
Ec = ecut(irl) - rm;
Pc = pcut(irl);
max_Ec = max($RELAX-CUTOFF,Ec);
max_Pc = max($RELAX-CUTOFF,Pc);
min_E = min(Ec,Pc); min_E = max($RELAX-CUTOFF,min_E);
min_E = $RELAX-CUTOFF;

Evac = shell_be(shell); "provides relevant binding energy for this vacancy"
"store some information in comin RELAX-FOR-USER. This just duplicates"
"various pieces of info so the user can access them."
Expand Down Expand Up @@ -5745,21 +5728,20 @@ IF (shell_egs > 4) [
vac = shell; Nvac = 0; np_save = np;
LOOP [ "from here to end of routine over all vacancies created"

"check if energy of vacancy < cutoffs OR no transitions from this shell"
IF( relax_initiator_q = 0 ) [ "particle initiating relaxations was photon"
IF( Evac < max_Pc | relax_ntran(vac) < 1 ) [
edep += Evac; "add energy of vacancy to edep"
edep_local = Evac; "set value of edep_local to energy of vacancy"
"check if energy of vacancy < cutoff OR no transitions from this shell"
IF( Evac < min_E | relax_ntran(vac) < 1 ) [
edep += Evac; "add energy of vacancy to edep"
edep_local = Evac; "set value of edep_local to energy of vacancy"

IF( relax_initiator_q = 0 ) [
"particle initiating relaxations was photon"
$AUSCALL($SPHOTONA); "call ausgab with iarg=33"
go to :VACANCY:; "exit loop and if Nvac still 0, exit routine"
]
] ELSE [ "particle initiating relaxations was electron"
IF( Evac < max_Ec | relax_ntran(vac) < 1 ) [
edep += Evac; "add energy of vacancy to edep"
edep_local = Evac; "set value of edep_local to energy of vacancy"
] ELSE [
"particle initiating relaxations was electron"
$AUSCALL($SELECTRONA); "call ausgab with iarg=34"
go to :VACANCY:; "exit loop and if Nvac still 0, exit routine"
]

go to :VACANCY:; "exit loop and if Nvac still 0, exit routine"
]

"prepare_alias_histogram is called from subroutine egs_init_relax and"
Expand Down
4 changes: 2 additions & 2 deletions HEN_HOUSE/user_codes/g/g.mortran
View file
Original file line number Diff line number Diff line change
Expand Up @@ -524,10 +524,10 @@ IF( e_muen > 0 ) [
OUTPUT e_tot,e_tot2, 100*e_tot2/e_tot;
(' Ave energy released per particle: ',1PE12.4, ' MeV +/-', 1PE10.3,
' =', 0PF7.3,' %');
OUTPUT e_brem,e_brem2, 100*e_brem2/e_tot;
OUTPUT e_brem,e_brem2, 100*e_brem2/e_brem;
(' Ave energy lost to bremsstrahlung: ',1PE12.4, ' MeV +/-', 1PE10.3,
' =', 0PF7.3,' %');
OUTPUT e_rad,e_rad2, 100*e_rad2/e_tot;
OUTPUT e_rad,e_rad2, 100*e_rad2/e_rad;
(' Ave energy lost to all radiation: ',1PE12.4, ' MeV +/-', 1PE10.3,
' =', 0PF7.3,' %');

Expand Down

0 comments on commit 6d487e7

Please sign in to comment.