Problem 1481

Created attachment 221 [details]
test.mac is a script for TestEm7 for Geant4 9.6.p02.

The attached file, test.mac, fails when run with TestEm7 using Geant4 9.6.p02.

The specific cause of the failure is a 0./0. (zero divided by zero) condition at line 577 of G4IonParametrisedLossModel.cc. Here is the line:

G4double factor =  (1.0 + (dEdxLimitParam/dEdxLimitBetheBloch-1.0)*(lowEnergyLimit/scaledKineticEnergy));

The failure happens because for certain parameters (ions, energies, small values of the range cuts) dEdxLimitParam == dEdxLimitBetheBloch == 0. It appears that the proper behavior is that the 0./0. should be replaced by 1. Therefore, as an emergency measure until the Geant4 developers produce an official fix, in my own code I have replaced the above line with this:

G4double factor =  ((dEdxLimitParam==dEdxLimitBetheBloch) ? 1. : 
        	(1.0 + (dEdxLimitParam/dEdxLimitBetheBloch-1.0)*(lowEnergyLimit/scaledKineticEnergy)));

This small change is still not defensive against the possibility that dEdxLimitParam is finite while dEdxLimitBetheBloch==0, but it does correct the catastrophic code failure illustrated by the attached test.mac script.

When this failure happens, the range table is populated by NaNs for energies above a certain point. Then the corresponding range table is populated in every bin by the value of the zero bin, which happens to be the same for both opt2 and opt4 when the above temporary code fix is applied.

After applying the change above to my code, opt4 produced tables similar to opt2, and both TestEm7 and my own "mred" code run to completion in all tests that I have done, including for alpha particles. I have had a moderately extensive email discussion with Vladimir Ivanchenko about this already. I have not done quantitative tests on the physics, nor am I interested at this specific moment in code accuracy for small cut values, although it is obviously a critical issue. Here, I am speaking specifically to the mechanical code failure described above. 

I am concerned with this issue and expended the effort to find this code bug because getting the physics right in small volumes and at relatively low energies is a high priority of the space-users community. The need is probably not yet quite as urgent as that of the DNA group, but it is certainly moving in that direction.

Modern FinFet transistors have dimensions in the few nm range, and data that our group will publish shortly confirm that radiation simulations in these volumes are already and important component of understanding their behavior when irradiated. We hope that eventually Geant4 applications can be constructed that will approach the simulation fidelity of Fortran Penelope in small volumes and at low energies.

Final note: My characterization of this as a critical issue is based on our application. I understand that for the larger community it may not be so.