Thermal neutron physics drastic change from 10.5 to 10.6

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?

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

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

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.

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

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

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.

Can anyone confirm whether this fix was added to 10.7 or its patches? The missing phi rotation does not appear to be in the 10.7 source code nor in the 10.7.p1 patch version I downloaded.

It should be fixed, the location in the code of the phi rotation has changed. What happened was that the phi rotation was copied from the general hadron elastic code to each and every specific hadron elastic process (some processes in the future will do something different), but unfortunately thermal scattering was forgotten.

1 Like

I have a similar problem with the geant4-10-07-patch-01 version.
Setup: orb (R=150cm) of G4_C, starting with neutrons of 1eV in the center.
Compared with FLUKA (pointwise + groupwise and Mcnp)
FLUKA+Mcnp agree well, while G4 is lower, not only for the total fluence
but also on the peak at kT=10meV instead of kT=25meV.

Have you tested the behavior with 10.5 and 10.6? I assume that you have followed the instructions for low-energy neutron simulations:
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/physicsProcess.html#high-precision-neutron-interactions-neutronhp

10.6 gives exactly the same result. 10.5 I didn’t try since on my setup I am using the external navigator and was not present in the 10.5.
Just as a note: I am not using the S(a,b) on C, I am using it as amorphous for all setups

If 10.6 and 10.7 give the same result, I think that you are having a different problem. But did you have all the low-energy neutron physics set up correctly?

Hi. I do not think there is a specific example on BNCT (I will write down to implement one). What you need to implement is the proton beam incident on the target and record the secondary radiation field. You should use QGSP_BIC_HP with thermal neutron cross sections activated. What I suggest is to study and start from the extended examples Hadr01 and 02 and change them appropriately.

I hope this helps, let me know how it goes,
cheers
Susanna

I am working on a project that involves studying thermal neutrons generated by exposing protons to beryllium and subsequently analyzing the thermal neutrons for the purpose of application BNCT. I intend to employ the physics list ‘QGSP_BERT_HP’ for this purpose. However, based on my research, it seems necessary to incorporate ‘ENDF/B’ into my project to facilitate the examination of the cross-section of these thermal neutrons. Is this inclusion truly necessary?
Thank you

Hello Reem,

The “HP” part of QGSP_BERT_HP and QGSP_BIC_HP physics describes the incident neutron interactions (from G4NDL). You do not need to add any additional datasets for incident neutron interactions.

As Susanna suggested, you should include the thermal neutron cross section model (to describe elastic interactions in thermal energy range, less than 4 eV). See extended example Hadr04 of how to activate the model (Hadr04 NeutronHPphysics.cc) and how to define a material to use the thermal neutron treatment (Hadr04 DetectorConstruction.cc).

You can also try QGSP_BIC_AllHP, which uses the G4TENDL database to describe the interactions of incident charged particle (including protons) based on evaluated datasets. You can compare the neutron energy yield of the proton+Be interaction to QGSP_BERT_HP/QGSP_BIC_HP.

Regards,
James