Problem 905

Summary: G4UrbanMscModel::ComputeGeomPathLength produce NAN
Product: Geant4 Reporter: vincenzo.innocente
Component: processes/electromagneticAssignee: Vladimir.Ivantchenko
Status: RESOLVED FIXED    
Severity: major CC: f.uhlig
Priority: P4    
Version: 8.1   
Hardware: PC   
OS: Linux   

Description vincenzo.innocente 2006-11-06 09:47:52 CET
I encountered this problem running N03
for 4.8.1p1 on SLC4 compiling geant4 with
g++4 -O2 -funsafe-math-optimizations -fno-math-errno -msse3 -mfpmath=sse,387

0x006f929c in G4UrbanMscModel::ComputeGeomPathLength (this=0x9072850) at src/
G4UrbanMscModel.cc:582
582         zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
(gdb) print tPathLength/currentRange
$2 = 1
(gdb) print tPathLength
$3 = 0.0095235382324220224
(gdb) print currentRange
$4 = 0.0095235382324220224

The problem is solved changing the code as follows
diff /afs/cern.ch/sw/geant4/releases/share/geant4.8.1.p01/source/processes/electromagnetic/
standard/src/G4UrbanMscModel.cc ../../../source/processes/electromagnetic/standard/src/
G4UrbanMscModel.cc
582c582,586
<     zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
---
>     G4double u = 1.-tPathLength/currentRange;
>     if (float(u)>0.0)
>       zmean = (1.-exp(par3*log(u)))/(par1*par3) ;
>     else
>       zmean =  1./(par1*par3) ;
Comment 1 vincenzo.innocente 2006-11-06 10:14:59 CET
more nan
Program received signal SIGFPE, Arithmetic exception.
0x00348a1e in G4UrbanMscModel::SampleCosineTheta (this=0x8867020, trueStepLength=0,
    KineticEnergy=1.4207528465018906) at src/G4UrbanMscModel.cc:763
763       lambdaeff = trueStepLength/currentTau;
(gdb) print trueStepLength
$1 = 0
(gdb) print currentTau
$2 = 0

I protected it with
lambdaeff = (trueStepLength>0) ? trueStepLength/currentTau : 0;
and
at the same place I've also a division by zero
in G4UrbanMscModel::SampleCosineTheta (this=0x8655020,
trueStepLength=0.015376148707796944,
    KineticEnergy=0.47799742059990835) at src/G4UrbanMscModel.cc:763
763       lambdaeff = (trueStepLength>0) ? trueStepLength/currentTau : 0;
(gdb) print trueStepLength>0
$1 = true
(gdb) print trueStepLength
$2 = 0.015376148707796944
(gdb) print currentTau
$3 = -0

I decided Not to trap them...
Comment 2 Vladimir.Ivantchenko 2006-11-27 07:04:59 CET
THe problem was identified and fixed by author; the results were not effected
Comment 3 Florian Uhlig 2007-08-10 11:43:28 CEST
*** Problem 962 has been marked as a duplicate of this problem. ***