EGSnrc
EGSnrc copied to clipboard
nrc-cnrc
Bremsstrahlung photon energy sampling: potential error in code?
I found another possible error in the code for bremsstrahlung, this time when doing the energy sampling (specifically if the default ibr_nist=0 option is used, that is the user wants to use the Bethe-Heitler cross sections, Coulomb-corrected above 50 MeV and Born approximation otherwise).
The error turns up when calculating the functions ϕ1 and ϕ2 on lines 435-436 in egsnrc.mortran. The values in the array dl1 are calculated later on in lines 2355-2378 in the subroutine fix_brems as 20.863 plus some some value, all divided by fmax. However, the number 20.863 should instead be 20.867 according to the EGSnrc manual PIRS-701 (equation 2.1.8 in the section on pair and triplet production), and according to the original article from which these approximations to the screening functions came:
This equally applies to fmax1 and fmax2 in lines 2351-2352, and also to the corresponding values for pair production in lines 2381-2404.
Essentially, wherever the value 20.863 is written in the subroutine fix_brems, I believe this value should be replaced with 20.867. That is, unless the value was consciously changed from 20.867 to 20.863 for some other reason.
Good catch. Let's review if there is indeed a reason for using a slightly different numerical value, since this is in the fix_brems routine with is meant to, well, fix the brems calculation! So it is not impossible that this change was intentional. I will dig through the file history and change log for this remote possibility. It is such a small difference though that I doubt it is more than just a typo or rounding off issue.
None of the lines containing 20.863 have been modified since EGSnrc was ported to CVS in 2003, and the only log message at that point was "Blah" (!). I will have to recover the pre-2003 RCS version control history to confirm further...
Interesting! It looks like that portion of the code (fix_brems) was written in January 2000. Given that it hasn't been modified since 2003, I would be surprised if it was changed some time between 2000 and 2003. Although perhaps there is a log when it was first committed indicated why the value of 20.863 was chosen.
Looking back through the log history (in rcs and sccs version control logs!), and I can confirm that the lines containing the 20.863 constant in the fix_brems subroutine first appeared as is in the code on 20 January 2000, and have not been modified since. So in the end I see no documented reason for using 20.863 instead of 20.867. Again, no significant impact on calculations, but it would be best if the code matched the documentation.
The original sccs log reveals:
D 1.15 00/01/20 16:48:07 iwan 15 14 00238/00180/04252
Replaced BREMS
$ sccs sccsdiff -r1.14 -r1.15 egsnrc.mortran
(...)
1747a1717,1803
> subroutine fix_brems;
> "******************************************************************"
> "
> " Calculates the parameter for the rejection function used in
> " the current implementation of bremsstrahlung sampling
> "
> " I Kawrakow, January 2000
> "
> "*******************************************************************"
>
> $IMPLICIT-NONE;
>
> ;COMIN/BREMPR,MEDIA/;
>
> $INTEGER medium,i;
> $REAL Zt,Zb,Zf,Zg,Zv,fmax1,fmax2,Zi,pi,fc,xi,aux,
> XSIF,FCOULC;
>
> DO medium = 1,nmed [
>
> Zt = 0; Zb = 0; Zf = 0;
> DO i=1, NNE(medium) [
> Zi = ZELEM(medium,i); pi = PZ(medium,i);
> fc = FCOULC(Zi); xi = XSIF(Zi);
> aux = pi*Zi*(Zi + xi);
> Zt = Zt + aux;
> Zb = Zb - aux*Log(Zi)/3;
> Zf = Zf + aux*fc;
> ]
> Zv = (Zb - Zf)/Zt; Zg = Zb/Zt;
> fmax1 = 2*(20.863 + 4*Zg) - 2*(20.029 + 4*Zg)/3;
> fmax2 = 2*(20.863 + 4*Zv) - 2*(20.029 + 4*Zv)/3;
> dl1(1,medium) = (20.863 + 4*Zg)/fmax1;
> dl2(1,medium) = -3.242/fmax1;
> dl3(1,medium) = 0.625/fmax1;
> dl4(1,medium) = (21.12+4*Zg)/fmax1;
> dl5(1,medium) = -4.184/fmax1;
> dl6(1,medium) = 0.952;
> dl1(2,medium) = (20.029+4*Zg)/fmax1;
> dl2(2,medium) = -1.93/fmax1;
> dl3(2,medium) = -0.086/fmax1;
> dl4(2,medium) = (21.12+4*Zg)/fmax1;
> dl5(2,medium) = -4.184/fmax1;
> dl6(2,medium) = 0.952;
> dl1(3,medium) = (20.863 + 4*Zv)/fmax2;
> dl2(3,medium) = -3.242/fmax2;
> dl3(3,medium) = 0.625/fmax2;
> dl4(3,medium) = (21.12+4*Zv)/fmax2;
> dl5(3,medium) = -4.184/fmax2;
> dl6(3,medium) = 0.952;
> dl1(4,medium) = (20.029+4*Zv)/fmax2;
> dl2(4,medium) = -1.93/fmax2;
> dl3(4,medium) = -0.086/fmax2;
> dl4(4,medium) = (21.12+4*Zv)/fmax2;
> dl5(4,medium) = -4.184/fmax2;
> dl6(4,medium) = 0.952;
(...)
A similar code block below for pair production, under the "and these in PAIR" comment, was added on 14 February 2000 without further comment:
D 1.26 00/02/14 11:12:07 iwan 26 25 00186/00280/04939
$ sccs sccsdiff -r1.25 -r1.26 egsnrc.mortran
(...)
1936a1943,1971
> "and these in PAIR"
> dl1(5,medium) = (3*(20.863 + 4*Zg) - (20.029 + 4*Zg));
> dl2(5,medium) = (3*(-3.242) - (-1.930));
> dl3(5,medium) = (3*(0.625)-(-0.086));
> dl4(5,medium) = (2*21.12+8*Zg);
> dl5(5,medium) = (2*(-4.184));
> dl6(5,medium) = 0.952;
> dl1(6,medium) = (3*(20.863 + 4*Zg) + (20.029 + 4*Zg));
> dl2(6,medium) = (3*(-3.242) + (-1.930));
> dl3(6,medium) = (3*0.625+(-0.086));
> dl4(6,medium) = (4*21.12+16*Zg);
> dl5(6,medium) = (4*(-4.184));
> dl6(6,medium) = 0.952;
> dl1(7,medium) = (3*(20.863 + 4*Zv) - (20.029 + 4*Zv));
> dl2(7,medium) = (3*(-3.242) - (-1.930));
> dl3(7,medium) = (3*(0.625)-(-0.086));
> dl4(7,medium) = (2*21.12+8*Zv);
> dl5(7,medium) = (2*(-4.184));
> dl6(7,medium) = 0.952;
> dl1(8,medium) = (3*(20.863 + 4*Zv) + (20.029 + 4*Zv));
> dl2(8,medium) = (3*(-3.242) + (-1.930));
> dl3(8,medium) = (3*0.625+(-0.086));
> dl4(8,medium) = (4*21.12+16*Zv);
> dl5(8,medium) = (4*(-4.184));
> dl6(8,medium) = 0.952;
(...)