Thermal neutron physics drastic change from 10.5 to 10.6

Physics Processes, Models and Cross Sections
#1

Quick summary: the behaviour of thermal neutrons in my simulation is very different when using Geant4 10.5 compared to Geant4 10.6.

I am simulating a monoenergetic 10 keV neutron source on a water phantom, and scoring the fluence of photons and thermal, epithermal and fast neutrons on a voxel grid. The following figures show the thermal and epithermal neutron fluence profiles respectively along the beam axis, normalized so that the maximum is 1. Results are provided for two fully independent code bases. A result calculated with MCNP used as reference is also shown.

lthermal_neutron_fluence_compare_geant4_versions

epithermal_neutron_fluence_compare_geant4_versions

Although it might not be obvious from the images, for epithermal neutrons all calculations agree, while for thermal neutrons all calculations using Geant4 10.6.0 or later have a different profile (also different to the MCNP reference).

All results have been calculated using G4ENDL-4.5. G4ENDL-4.6 has shown similar results when tested.

Any ideas? Is this a bug or a feature?

#2

I looked at the neutron cross-sections and they are the same, in case this helps troubleshooting.

#3

Hello,

I would advice to compared initialization printout for hadronic physics in Geant4 10.5 and 10.6 and be sure, that in 10.6 thermal neutrons are enabled. I mean printout at the beginning of simulation.

VI

#4

Thank you Vladimir, good point!

This is from 10.5:

====================================================================
                  HADRONIC PROCESSES SUMMARY (verbose level 1)

---------------------------------------------------
                           Hadronic Processes for gamma

  Process: photonNuclear
        Model:            BertiniCascade: 0 eV  ---> 3.5 GeV
     Cr_sctns:            PhotoNuclearXS: 0 eV  ---> 100 TeV

---------------------------------------------------
                           Hadronic Processes for neutron

  Process: hadElastic
        Model:          NeutronHPElastic: 4 eV  ---> 20 MeV
        Model: NeutronHPThermalScattering: 0 eV  ---> 4 eV 
     Cr_sctns: NeutronHPThermalScatteringData: 0 eV  ---> 4 eV 
     Cr_sctns:        NeutronHPElasticXS: 0 eV  ---> 20 MeV
     Cr_sctns:            GheishaElastic: 0 eV  ---> 100 TeV

  Process: neutronInelastic
        Model:        NeutronHPInelastic: 0 eV  ---> 20 MeV
     Cr_sctns:      NeutronHPInelasticXS: 0 eV  ---> 20 MeV
     Cr_sctns:          GheishaInelastic: 0 eV  ---> 100 TeV

  Process: nCapture
        Model:          NeutronHPCapture: 0 eV  ---> 20 MeV
     Cr_sctns:        NeutronHPCaptureXS: 0 eV  ---> 20 MeV
     Cr_sctns:        G4NeutronCaptureXS: 0 eV  ---> 100 TeV

================================================================

And this is from 10.6:

====================================================================
                  HADRONIC PROCESSES SUMMARY (verbose level 1)

---------------------------------------------------
                           Hadronic Processes for gamma

  Process: photonNuclear
        Model:            BertiniCascade: 0 eV  ---> 3.5 GeV
     Cr_sctns:            PhotoNuclearXS: 0 eV  ---> 100 TeV

---------------------------------------------------
                           Hadronic Processes for neutron

  Process: hadElastic
        Model:          NeutronHPElastic: 4 eV  ---> 20 MeV
        Model: NeutronHPThermalScattering: 0 eV  ---> 4 eV 
     Cr_sctns: NeutronHPThermalScatteringData: 0 eV  ---> 4 eV 
     Cr_sctns:        NeutronHPElasticXS: 0 eV  ---> 20 MeV

  Process: neutronInelastic
        Model:        NeutronHPInelastic: 0 eV  ---> 20 MeV
     Cr_sctns:      NeutronHPInelasticXS: 0 eV  ---> 20 MeV

  Process: nCapture
        Model:          NeutronHPCapture: 0 eV  ---> 20 MeV
     Cr_sctns:        NeutronHPCaptureXS: 0 eV  ---> 20 MeV
     Cr_sctns:        G4NeutronCaptureXS: 0 eV  ---> 100 TeV

================================================================

The difference between the two seems to be that GheishaElastic and GheishaInelastic are gone. I have no idea what these are… Do you think that might be the cause? I did not do anything to remove those though.

#5

Hi,

I think I found the source of the problem.

Going from 10.5 to 10.6, the random rotation around the z-axis was removed in /source/processes/hadronic/processes/src/G4HadronElasticProcess.cc (i.e. “outdir.rotate(phi, it);” and “pdir.rotate(phi, it);” were removed prior to applying “rotateUz(indir)”). This removal seems to cause the neutrons to change direction only in the transverse plane, example (10 keV neutrons in water, thermal neutron recipe in use, G4_10.06.p02, G4NDL4.6):
20201215_G4_10.6.2_G4NDL4.6_longitudinalview

Running the same user code with and older G4 version (G4_10.05.p01, G4NDL4.6) the longitudinal scattering looks normal, i.e. not confined to the transverse plane.

When the azimuthal rotation described above is reintroduced into 10.6, the longitudinal scattering looks again normal, similar to 10.5.

With this type of bugfix, the G4_10.6.p02 and MCNP cell fluxes again agree nicely for both epithermal and thermal neutrons as well as photons as a function of depth along the center of the beam axis.

/LW

1 Like
#6

The visualization with the same code in 10.06.p02, but reintroducing the azimuthal rotation looks like the following:

20201216_G4_10.6.2_bugfix_G4NDL4.6_longitudinalview
20201216_G4_10.6.2_bugfix_G4NDL4.6_longitudinalview783Ă—522 4.98 KB

#7

I reported a bug:

This figure summarizes the problem that @lauriw’s figures revealed:
angle_z_histogram

The x-axis is output angle of neutrons in hadronic ellastic collisions with respect to the Z-axis. The y axis is probability density.