Skip to content

Commit

Permalink
Merge pull request #212
Browse files Browse the repository at this point in the history 
Fix kerma calculation in g and DOSRZnrc.
  • Loading branch information
@ftessier
ftessier committed Jan 24, 2017
2 parents 7a9058b + bb055ad commit 5ce91b4
Show file tree
Hide file tree
Showing 4 changed files with 184 additions and 27 deletions.
14 changes: 8 additions & 6 deletions HEN_HOUSE/src/egsnrc.macros
View file
Expand Up @@ -203,10 +203,8 @@ PARAMETER $MXSINC=1002; "ANGLE INTERVALS IN (0,5*PI/2) FOR SINE"
; "---------- BUFFER FLUSH SEMICOLON ----------"

"THE FOLLOWING ARE PARAMETERS FOR AUSGAB CALLS:"
"Ali:photonuc, 2 lines (used to be 29 and 24, now 30, 27"
"EMH: Should actually be 31 and 26"
PARAMETER $MXAUS=31; "CHANGE IF MORE AUSGAB CALLS ARE ADDED"
PARAMETER $MXAUSM5=26; "SET THIS TO $MXAUS VALUE LESS 5"
PARAMETER $MXAUS=33; "CHANGE IF MORE AUSGAB CALLS ARE ADDED"
PARAMETER $MXAUSM5=28; "SET THIS TO $MXAUS VALUE LESS 5"

"FIRST FIVE AUSGAB CALLS BELOW TURNED ON IN BLOCK DATA (DEFAULT)"
PARAMETER $TRANAUSB=0; "BEFORE TRANSPORT"
Expand Down Expand Up @@ -244,6 +242,9 @@ PARAMETER $AUGERTRA=27; "An Auger transition just occured"
" note that 28 is already used for positron annih at rest - see above"
PARAMETER $PHOTONUCAUSB=29; "BEFORE PHOTONUCLEAR EVENT"
PARAMETER $PHOTONUCAUSA=30; "AFTER PHOTONUCLEAR EVENT"
PARAMETER $EIIB=31; "Before EII"
PARAMETER $EIIA=32; "After EII"

; "---------- BUFFER FLUSH SEMICOLON ----------"

REPLACE {$AUSCALL(#);} WITH
Expand Down Expand Up @@ -834,12 +835,13 @@ REPLACE {COMIN/CH-Steps/;} WITH
"------------------------------------------------------------------"
REPLACE {;COMIN/EPCONT/;} WITH
{;
COMMON/EPCONT/EDEP,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,
COMMON/EPCONT/EDEP,EDEP_LOCAL,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,
RHOF,EOLD,ENEW,EKE,ELKE,GLE,E_RANGE,
x_final,y_final,z_final,
u_final,v_final,w_final,
IDISC,IROLD,IRNEW,IAUSFL($MXAUS);
$ENERGY PRECISION EDEP; "energy deposition in MeV"
$ENERGY PRECISION EDEP, "energy deposition in MeV"
EDEP_LOCAL "local energy deposition in MeV"
$REAL TSTEP, "distance to a discrete interaction"
TUSTEP, "intended step length, befor check with geometry"
USTEP, "transport distance calculated from TUSTEP"
Expand Down
57 changes: 50 additions & 7 deletions HEN_HOUSE/src/egsnrc.mortran
View file
Expand Up @@ -27,6 +27,7 @@
" Contributors: Ernesto Mainegra-Hing "
" Blake Walters "
" Frederic Tessier "
" Reid Townson "
" "
"#############################################################################"
" "
Expand Down Expand Up @@ -780,7 +781,7 @@ IF( ibcmpl = 1 | ibcmpl = 3 ) [
call relax(Uj,shn_array(j),iz_array(j));
"relax will put all particles with energies above ecut,pcut on the "
"stack, the remaining energy will be scored in edep and deposited "
"localy (via the call to ausgab below) "
"locally (via the call to ausgab below) "
]
ELSE [ edep = Uj; ]
IF( edep > 0 ) [ $AUSCALL($PHOTXAUS); "generates IARG = 4 call" ]
Expand Down Expand Up @@ -2543,7 +2544,7 @@ $COMIN-MOLLER; "DEFAULT REPLACEMENT PRODUCES THE FOLLOWING:
$DEFINE-LOCAL-VARIABLES-MOLLER;

$REAL sigm,pbrem,rsh,Uj,sig_j;
$INTEGER lelke,iele,ish,nsh,ifirst,i,jj,iZ;
$INTEGER lelke,iele,ish,nsh,ifirst,i,jj,iZ,iarg;

"IRCODE=1; appears to be unused, IK Oct 97"
NPold = NP; "Set the old stack counter"
Expand Down Expand Up @@ -2574,7 +2575,9 @@ IF( eii_flag > 0 & eii_nsh(medium) > 0 ) [
sig_j = sig_j*pz(medium,iele)*eii_cons(medium);
rsh = rsh - sig_j;
IF( rsh < 0 ) [
$AUSCALL($EIIB);
call eii_sample(ish,iZ,Uj);
$AUSCALL($EIIA);
return;
]
]
Expand Down Expand Up @@ -5204,6 +5207,9 @@ IF( energy <= e_cut ) [
edep = edep + energy; "We assume that edep is zeroed "
"(or set to the appropriate value in the routine "
"calling RELAX "

edep_local = energy;
$AUSCALL($AUGERTRA);
return;
]

Expand All @@ -5216,13 +5222,17 @@ e_check = 0; e_array(n_vac) = energy;
shell = vac_array(n_vac); Ei = e_array(n_vac); n_vac = n_vac - 1;
IF( Ei <= e_cut ) [ " Below cut-off -> local absorption "
edep = edep + Ei;

edep_local = Ei;
$AUSCALL($AUGERTRA);

IF( n_vac > 0 ) goto :START: ;
EXIT;
]
"Set the relax_user common block variables, IK March 22 2004"
ish_relax = shell; u_relax = Ei;
IF( shell = $MXSHELL ) [ "This is N-shell vacancy -> just produce Auger"
IF( Ei > e_cut ) [
IF( Ei > ekcut ) [
np = np + 1;
$CHECK-STACK(np,'RELAX');
e(np) = Ei + prm; iq(np) = -1;
Expand All @@ -5237,7 +5247,12 @@ e_check = 0; e_array(n_vac) = energy;
ELSE [ u(np) = 0; v(np) = 0; w(np) = 1; ]
$AUSCALL($AUGERTRA);
]
ELSE [ edep = edep + Ei; ]
ELSE [
edep = edep + Ei;

edep_local = Ei;
$AUSCALL($AUGERTRA);
]
IF( n_vac > 0 ) goto :START: ;
EXIT;
]
Expand Down Expand Up @@ -5276,6 +5291,15 @@ e_check = 0; e_array(n_vac) = energy;
]
IF( Ex <= elcut ) [ "Below cut-off"
edep = edep + Ex;

IF( finala < 10 ) [
edep_local = Ex;
$AUSCALL($FLUORTRA);
]
ELSE [
edep_local = Ex;
$AUSCALL($AUGERTRA);
]
]
ELSE [
np = np + 1;
Expand Down Expand Up @@ -5598,7 +5622,10 @@ rfu_E0 = Evac;
vac = shell; Nvac = 0; np_save = np;
LOOP [
IF( Evac < min_E | relax_ntran(vac) < 1 ) [
edep += Evac; go to :VACANCY:;
edep += Evac;
edep_local = Evac;
$AUSCALL($AUGERTRA);
go to :VACANCY:;
]
new_state = sample_alias_histogram(relax_ntran(vac),
relax_prob(relax_first(vac)),
Expand Down Expand Up @@ -5637,6 +5664,7 @@ LOOP [
]
]
Ex = Evac - Ef;
edep_local = 0;
IF( Ex > Ecc ) [
np += 1; $CHECK-RELAX-STACK(np,'new_relax');
"$egs_info('(a,i9,a,f7.3)','-> np=',np,' iq=',iqf);"
Expand All @@ -5657,7 +5685,17 @@ LOOP [
e(np) = Ex + rm; $AUSCALL($AUGERTRA);
]
"$egs_info('(a,i2,a,f7.3,a)','-> iq=',iqf,' E=',Ex*1000.,' keV');"
] ELSE [ edep += Ex; ]
] ELSE [
edep += Ex;

IF( iqf = 0 ) [
edep_local = Ex;
$AUSCALL($FLUORTRA);
] ELSE [
edep_local = Ex;
$AUSCALL($AUGERTRA);
]
]
:VACANCY:;
IF( Nvac = 0 ) EXIT;
vac = vacs(Nvac); Evac = shell_be(vac); Nvac -= 1;
Expand Down Expand Up @@ -7294,14 +7332,19 @@ sinthe = dsqrt(1-dcosth); costhe = dsqrt(dcosth);
call uphi(2,1);

pese2 = wx - Uj + prm;
edep_local = 0;
IF( pese2 > ae(medium) ) [
$CHECK-STACK(np+1,'eii_sample');
np = np+1; e(np) = pese2;
dcosth = h1*(pese2-prm)/(pese2+prm);
sinthe = -dsqrt(1-dcosth); costhe = dsqrt(dcosth);
iq(np) = -1; call uphi(3,2);
edep = 0;
] ELSE [ edep = wx - Uj; ]
] ELSE [
edep = wx - Uj;
edep_local = edep;
$AUSCALL($AUGERTRA);
]

call relax(Uj,ish,iZ);

Expand Down
66 changes: 58 additions & 8 deletions HEN_HOUSE/user_codes/dosrznrc/dosrznrc.mortran
View file
Expand Up @@ -30,6 +30,7 @@
" Iwan Kawrakow "
" Blake Walters "
" Ernesto Mainegra-Hing "
" Reid Townson "
" "
"#############################################################################"
" "
Expand Down Expand Up @@ -91,6 +92,11 @@
" depth (bins were offset). Added lots of comments as I tried to relearn the "
" code and made many comments in lower case to make them more readable. "
" "
" Reid Townson, 2017: "
" "
" Energy depositions for kerma calculations below the cut-off are now sorted "
" out using AUSGAB to include only Auger from relaxations, not fluorescence. "
" "
"#############################################################################"


Expand Down Expand Up @@ -2078,7 +2084,7 @@ SCKERMA_TMP($MAXZREG,$MAXRADII,$MAXIT),
SCKERMA_TMPOLD($MAXZREG,$MAXRADII,$MAXIT),
SCSTP_TMP,SCDSTP_TMP,
MXNP,IFULL,ISTORE,IKERMA, IWATCH,IOOPTN,IOUTSP,
IPHR($MXREG),MAXBIN,NCOMPT;
IPHR($MXREG),MAXBIN,NCOMPT,BEFORE_PHOTO,BEFORE_EII;
real*8 SCDOSE,SCDOSE2,SCKERMA,SCKERMA2,SCDOSEtoKERMA2,SCPDST,SCPDST2,
SCPCUM,SCPCUM2,SCPTOT,SCPTOT2,SCPHEN,SCPHEN2,
SCDFBK,SCDFBK2,SCDFEP,SCDFEP2,SCDFDIFF,SCDFDIFF2,
Expand All @@ -2089,7 +2095,7 @@ $REAL DFEN, PHENER,WT1OLD,BINTOP,SLOTE,DELTAE,
SCDOSE_TMP,SCKERMA_TMP,SCKERMA_TMPOLD,SCSTP_TMP,SCDSTP_TMP;
$INTEGER MXNP,IFULL,ISTORE,IKERMA,
IWATCH,IOOPTN,IOUTSP,
IPHR,MAXBIN,NCOMPT;
IPHR,MAXBIN,NCOMPT,BEFORE_PHOTO,BEFORE_EII;
}
"
"SCDOSE(2) accumulates energy deposited (energy deposited^2),
Expand Down Expand Up @@ -2547,15 +2553,22 @@ DO IZ=1,NZ[

;
"SET UP AUSGAB CALLS"
DO J=1,5[IAUSFL(J)=1;]DO J=6,25[IAUSFL(J)=0;] "NORMAL EXECUTION"
DO J=1,5[IAUSFL(J)=1;]DO J=6,33[IAUSFL(J)=0;] "NORMAL EXECUTION"

IF(IFULL = 1 | IFULL = 3 | IKERMA = 1) [
"need to call ausgab to set flag after photon interactions"
"IKERMA=1 means scoring KERMA"
/IAUSFL(6),IAUSFL(17),IAUSFL(19),IAUSFL(21)/=1;
"17 => call after pair, 19=>call after compton, 21=>call after photoelectric"
"6 => after each step"
"for KERMA, rayleigh scatter has no effect"
iausfl(6) = 1; "after each step"
iausfl(17) = 1; "after pair"
iausfl(18) = 1; "before compt"
iausfl(19) = 1; "after compt"
iausfl(20) = 1; "before photo"
iausfl(21) = 1; "after photo"
iausfl(10) = 1; "after Moller (to count radiative losses due to EII)"
iausfl(28) = 1; "after Auger"
iausfl(32) = 1; "before eii"
iausfl(33) = 1; "after eii"
]

IF(IFULL = 4) [
Expand Down Expand Up @@ -3774,8 +3787,28 @@ ELSEIF(IFULL = 3)["score scattered dose separately"
IF (IKERMA = 1) ["want to score KERMA"
"we score kerma for scattered component too if ifull=3"
"the kerma is part of scattered kerma if latch of initial photon is not zero"

IF (IARG = 4 & ~BTEST(LATCH(NP),8))["local energy deposition"

"Here we keep track of whether we are currently inside a PE/Compton"
"event, or EII. This lets us use edep_local to account for kerma"
"contributions below threshold energies depending on the interaction type."
IF( iarg = 19 | iarg = 17 ) [
before_photo = 1;
return;
]
IF( iarg = 20 | iarg = 18 ) [
before_photo = 0;
]
IF( iarg = 31 ) [
before_eii = 1;
return;
]
IF( iarg = 32 ) [
before_eii = 0;
return;
]

IF (IARG = 4 & ~BTEST(LATCH(NP),8) & before_eii = 0 & before_photo = 0)[
"local energy deposition"
"include deposited energy as kerma"
" I am not sure why the BTEST was there??? Currently bit 8 is"
"not set, but was there some reason for it to be there?? DR"
Expand All @@ -3784,6 +3817,21 @@ IF (IKERMA = 1) ["want to score KERMA"
$SCOREDK(SCKERMA,(IZD,IXD,2):WT(NP)*EDEP);
]
]

"depositing sub-threshold energy for PE & Compton"
"Auger only"
IF( iarg = 27 ) [
IF( before_photo = 1 ) [
$SCOREDK(SCKERMA,(IZD,IXD,1):WT(NP)*edep_local);
IF(IFULL=3 & (BTEST(LATCH(NP),6) | BTEST(LATCH(NP),7)))[
$SCOREDK(SCKERMA,(IZD,IXD,2):WT(NP)*edep_local);
]
return;
] ELSE [
return;
]
]

IF (IARG = 16)["pair event just occured"
IF(NP>NPold | i_survived_rr > 0)[
"i_survived_rr, 0=> all particles survive RR,"
Expand Down Expand Up @@ -4669,6 +4717,8 @@ TIMMAX=VALUE(NUM_MXTIME,1);
IFULL=VALUE(NUM_IFULL,1);
STATLM=VALUE(NUM_STATLM,1);
IKERMA=VALUE(NUM_SCKERMA,1);
BEFORE_PHOTO=0;
BEFORE_EII=0;

IF(IWATCH=0 & NCASE < $NCASEMIN)[NCASE=$NCASEMIN;]

Expand Down

0 comments on commit 5ce91b4

Please sign in to comment.