Problem 1856

Created attachment 397 [details]
A complete geant4.9.2 project

Overview: 

When neutron thermal scattering is simulated with G4NeutronHPThermalScattering, I register G4NeutronHPThermalScattering in the hadron elastic process G4UHadronElasticProcess (note: not G4HadronElasticProcess, which will be discussed) with the following code in the physics list:
{{
        // elastic scattering
        G4UHadronElasticProcess* theNeutronElasticProcess =
          new G4UHadronElasticProcess;
        G4LElastic* theElasticModel1 = new G4LElastic;
        G4NeutronHPElastic * theElasticNeutron = new G4NeutronHPElastic;
        G4NeutronHPThermalScattering* thermal = new G4NeutronHPThermalScattering;//
        theNeutronElasticProcess->RegisterMe(theElasticModel1);
        theElasticModel1->SetMinEnergy(19*MeV);
        theNeutronElasticProcess->RegisterMe(theElasticNeutron);
        theElasticNeutron->SetMinEnergy(4*eV);//
        theNeutronElasticProcess->RegisterMe(thermal);//
        G4NeutronHPElasticData * theNeutronData = new G4NeutronHPElasticData;
        G4NeutronHPThermalScatteringData* thermalData = new G4NeutronHPThermalScatteringData; //
        theNeutronElasticProcess->AddDataSet(theNeutronData);
        theNeutronElasticProcess->AddDataSet(thermalData); //
        pmanager->AddDiscreteProcess(theNeutronElasticProcess);
}}

By printing the neutron vertex and momentum infomation via the stepping action (detailed code in the attachment), a typical output is attached below:
{{
trkID= 1     step# 1     mat= Water_ts process=  energy(eV)= 779711     x= -2000      y= 0          z= 0          px= -3.70824e+07 py= -8.89919e+06 pz= 3.38872e+06
trkID= 1     step# 2     mat= Water_ts process= hElastic energy(eV)= 722146     x= -1991.89   y= 0          z= 0          px= -3.49147e+07 py= -9.64319e+06 pz= -6.74548e+06

...

trkID= 1     step# 57    mat= Water_ts process= hElastic energy(eV)= 0.026832   x= -2087.58   y= 19.9279    z= -38.915    px= 3635.08    py= 6099.76    pz= -5.65986e-13
trkID= 1     step# 58    mat= Water_ts process= hElastic energy(eV)= 0.0216959  x= -2085.54   y= 19.4551    z= -38.915    px= -5563.92   py= 3132.46    pz= 1.25499e-14
trkID= 1     step# 59    mat= Water_ts process= hElastic energy(eV)= 0.055645   x= -2084.36   y= 21.429     z= -38.915    px= 1585.03    py= -10102.1   pz= -1.24557e-14
trkID= 1     step# 60    mat= Water_ts process= hElastic energy(eV)= 0.0633198  x= -2085.5    y= 22.0675    z= -38.915    px= 9631.47    py= 5120.64    pz= 4.34345e-15
trkID= 1     step# 61    mat= Water_ts process= hElastic energy(eV)= 0.0376553  x= -2085.34   y= 21.073     z= -38.915    px= -4905.72   py= 6833.24    pz= -4.47491e-16
trkID= 1     step# 62    mat= Water_ts process= hElastic energy(eV)= 0.0375492  x= -2084.1    y= 21.732     z= -38.915    px= 2067.96    py= -8141.47   pz= 4.15985e-16
trkID= 1     step# 63    mat= Water_ts process= hElastic energy(eV)= 0.0429501  x= -2085.66   y= 23.8957    z= -38.915    px= 8920.6     py= -1063.83   pz= 1.59818e-16
trkID= 1     step# 64    mat= Water_ts process= hElastic energy(eV)= 0.00740901 x= -2085.54   y= 23.4224    z= -38.915    px= -279.164   py= 3720.83    pz= -1.27695e-17
trkID= 1     step# 65    mat= Water_ts process= hElastic energy(eV)= 0.0106628  x= -2082.43   y= 23.0514    z= -38.915    px= 2491.33    py= 2542.32    pz= -2714.17  
trkID= 1     step# 66    mat= Water_ts process= hElastic energy(eV)= 0.00467495 x= -2082.49   y= 23.8685    z= -38.915    px= -2652.71   py= 199.739    pz= 1306.94   
trkID= 1     step# 67    mat= Water_ts process= hElastic energy(eV)= 0          x= -2081.33   y= 25.0499    z= -40.1763   px= -0         py= 0          pz= 0 
}}

As one can see that starting with a sub-MeV kinetic energy, the neutron got scattered and slowed down to sub-eV very soon. The G4NeutronThermalScattering is the major physics in this energy region. Once the neutron got thermlized, the movement should be random in x,y and z-ditection. However, in the above simulation, neutron often stops in z-direciton and the z-momentum is very small. The issue is slightly better in some events, but it still exist.

However, when I repalced G4UHadronElasticProcess with G4HadronElasticProcess. The problem disappear. Amazing! But what's wrong with G4UHadronElasticProcess? I have submitted a tarball of this test project in the attachment for the experts to repeat the problem very quickly.


Steps to Reproduce: 
1) Un-tar the attachment.
This is a project modified from examples/novince/N01. The physics list is borrowed from examples/advanced/underground_physics/src/DMXPhysicsList.cc with some least modification. All the detector material has been changed to water as I need to test G4NeutronHPThermalScattering in water (300kelvin, 1bar).

2) Compile the package and run it.
The information is printed via stepping action, and 10 events will be simulated. One can change the total number of events in the exampleN01.cc.

3) Replace G4UHadronElasticProcess with G4HadronElasticProcess in src/DMXPhysicsList.cc (line 656 and 657).
Compile it again and run it.

Actual Results (from step 2):
See the description in overview. The momentum of neutron in z is almost zero for many steps.

Expected Results (from step 3):
This is a tyipical output by using G4HadronElasticProcess. It looks good.
{{
...

trkID= 1     step# 140   mat= Water_ts process= HadronElastic energy(eV)= 0.0523288  x= -1936.55   y= -22.74     z= 2.37527    px= -9905.26   py= 423.212    pz= -198.393  
trkID= 1     step# 141   mat= Water_ts process= HadronElastic energy(eV)= 0.0395362  x= -1937.63   y= -25.8565   z= 0.799538   px= -5880.33   py= -5624.65   pz= -2842.32  
trkID= 1     step# 142   mat= Water_ts process=  energy(eV)= 0.0382156  x= -1939.3    y= -25.7852   z= 0.766092   px= 4113.36    py= -6180.66   pz= -4085.58  
trkID= 1     step# 143   mat= Water_ts process= HadronElastic energy(eV)= 0.0327191  x= -1940.15   y= -26.6015   z= 0.353599   px= -2006.48   py= -6892.97   pz= 3153.48   
trkID= 1     step# 144   mat= Water_ts process= HadronElastic energy(eV)= 0.0320088  x= -1940.02   y= -26.7957   z= 0.225194   px= 2123.83    py= -6495.21   pz= 3667.47   
trkID= 1     step# 145   mat= Water_ts process= HadronElastic energy(eV)= 0.0307962  x= -1940.16   y= -27.2823   z= 0.447803   px= 2458.02    py= -6821.2    pz= -2302.04  
trkID= 1     step# 146   mat= Water_ts process= HadronElastic energy(eV)= 0.0510382  x= -1939.24   y= -30.108    z= 2.04332    px= 809.259    py= -7246.44   pz= 6537.71   
trkID= 1     step# 147   mat= Water_ts process= HadronElastic energy(eV)= 0.0462862  x= -1939.11   y= -30.4635   z= 1.92335    px= -7965.02   py= -2397.66   pz= -4217.54  
trkID= 1     step# 148   mat= Water_ts process= HadronElastic energy(eV)= 0.0448166  x= -1938.44   y= -36.4526   z= 7.3267     px= -8274.99   py= 2927.28    pz= -2678.03  
trkID= 1     step# 149   mat= Water_ts process= HadronElastic energy(eV)= 0.0450263  x= -1940.06   y= -36.9395   z= 6.47028    px= 3366.29    py= 6233.06    pz= -5867.48  
trkID= 1     step# 150   mat= Water_ts process= HadronElastic energy(eV)= 0.0388195  x= -1940.3    y= -36.8528   z= 6.39098    px= -5254.81   py= 1417.47    pz= -6582.15  
trkID= 1     step# 151   mat= Water_ts process= HadronElastic energy(eV)= 0.0399259  x= -1940.26   y= -36.7811   z= 6.32355    px= -6587      py= 5490.15    pz= 1222.94   
trkID= 1     step# 152   mat= Water_ts process= HadronElastic energy(eV)= 0.035251   x= -1940.37   y= -36.7519   z= 6.18788    px= 4649.5     py= 1648.99    pz= -6473.34  
trkID= 1     step# 153   mat= Water_ts process= HadronElastic energy(eV)= 0.0405012  x= -1940.63   y= -36.5368   z= 6.2358     px= 7135.46    py= 4976.06    pz= -656.569  
trkID= 1     step# 154   mat= Water_ts process= HadronElastic energy(eV)= 0          x= -1939.19   y= -36.0268   z= 4.23378    px= 0          py= 0          pz= -0   
}}

Platform and package versions:
Tested on geant4.9.2 first. Since G4UHadronElasticProcess has been removed from latest geant (>=4.9.5), I only make test up to geant4.9.4.
All my three tests give me consistent result.
1) geant4.9.2.p01, CentOS 6.5, CLHEP 2.0.4.2, G4NDL 3.13
2) geant4.9.2.p01, Ubuntu 14.04, CLHEP 2.0.4.2, G4NDL 3.13
3) geant4.9.4.p02, CentOS 6.2, CLHEP 2.1.0.1, G4NDL 3.14