Coefficient error in subroutine msdist_pII

[jantolak]

In egsnrc.mortran, subroutine msdist_pII at line 4092, it appears that the polynomial in the denominator is incorrect. Here is the current code at those lines.

4091    "e       = e * (1 - epsilonp*epsilonp*((6+tau*(10+5*tau))/(tau+1)/(tau+2))/24);
4092    e       = e * (1 - epsilonp*epsilonp*(6+10*tau+5*tau2)/(24*tau2+48*tau+72));

Assuming that the expression on line 4091 is correct, the denominator in the last term should be
24*(tau+1)*(tau+2) = 24*(tau2+3*tau+2) = (24*tau2+72*tau+48): the 72 and 48 are not in the right place in the code. Therefore, line 4092 should read

4092    e       = e * (1 - epsilonp*epsilonp*(6+10*tau+5*tau2)/(24*tau2+72*tau+48));

It appears that this error makes very little difference, which is likely why it has not been noticed. The error was actually pointed out to me by a trainee who was working on his own MC code for electron scattering and wanted to compare to an existing code as a baseline.

[ftessier]

Certainly an unfortunate typo!! Fixed in pull request #134.

Accrding to equation A6 of the 2000 Kawrakow paper (Kawrakow, Iwan. "Accurate condensed history Monte Carlo simulation of electron transport. I. EGSnrc, the new EGS4 version." Medical physics 27.3 (2000): 485-498), the denominator should indeed read (24*tau2+72*tau+48), from the expansion of 24*(tau+1)*(tau+2), where tau represents the kinetic energy in units of the rest mass.

For tau > 0 the incurred relative error in the polynomial [ (tau^2+2*tau+3)/(tau^2+3*tau+2) ] - 1 varies between a maximum of 50% at tau = 0 and a minimum of about 0.9 at tau ~= 3.449, and approaches 0 for large energies.

Fortunately, this mistake appears in a second-order term in epsilonp*epsilonp which mitigates its ultimate impact. The relative error in e and other dependent quantities in the code is found to be at most about 1%.