Incorrect output energy using LEND

[jamesv]

Hello,

I am doing some simple simulations with mono energetic neutrons comparing hadronic physics lists including QGSP_BICHP, QGSP_BERTHP, QGSP_BICAllHP, QGSP_INCLXXHP, Shielding and ShieldingLEND.

When I am using ShieldingLEND (which is Shielding but with LEND instead of HP data libraries), I have found some issues with boron capture energies.

For a simple test, I used the ‘radioprotection’ example and a 1cm3 cube of natural boron (G4_B from NISTManager) in vacuum. I checked the output using /tracking/verbose 1 and saw that incorrect energies were being produced from the neutronInelastic process with Boron-10, the kinetic energy of the particles are double what they should be, please see ‘screenshot 1’ below. The neutron physics that is built when using ‘ShieldingLEND’ physics list option is shown in ‘screenshot 2’.

For reference, the secondary product energies with the 94% B10(n,a)Li7 capture reaction should be alpha=1.47MeV, Li7=0.84MeV, gamma=0.478MeV. The Q value is 2.78MeV. We see that the gamma is correct (it is stated as a discrete energy in the ENDF/LEND data file), but alpha and Li7 seem to be doing some sort of energy balancing that is failing.

I also tried other materials from NIST such as G4_BORON_CARBIDE and defining the B-10 isotope explicitly using G4Isotope with Z, N, A numbers. I also tried changing the temperature of the material: 0K, 273K, 300K but it was always the same result.

It is occurring for all the monoenergetic neutrons I tested: 0.0025 eV, 0.1 eV, 1eV, 10 eV, 100eV, 1keV, 10keV, 100keV, 1MeV, 2MeV, 5MeV.

Other interactions such as H1(n,G)H2 have the expected results compared to literature, it is just Boron-10 that is wrong. Interestingly, the number of boron capture reactions per incident neutron is the same for Shielding and ShieldingLEND, just the output KE being different for alpha and Li7.

I then thought that there may be a bug with the overlapping of many cross sections seen here for neutronInelastic. So I then used no prebuilt hadronic physics list and added the LEND model to the neutron inelastic process manually to force only this one to be used as shown above in ‘screenshot 3’.

Unfortunately the same result was given for the boron capture reactions after this physics change. So this tells us that the incorrect energies is due to the LEND inelastic model or its data set.

I have tested this using Geant4 versions: 10.4.p02, 10.5.p01, 10.6.p02. I have also tried on different operating systems: Mac OSX 10.15, CentOS 7 and CentOS 8. I have discussed with other users and they produce the same incorrect KE result independently in their simulations.

Is anyone able to provide any information about this?

[jamesv]

Screenshot 1:

[jamesv]

Screenshot 2, neutron physics processes for ShieldingLEND:

[jamesv]

Screenshot 3, single process for neutrons to only use LEND model:

[jamesv]

I have done some more testing and found this affects other reactions that have 2+ bodies in the output, for example, Li6(n,a)H3, which also has incorrect output energy like the B10(n,a)Li7 reaction.

To make it more familiar, I used example Hadr03 for better reproducibility, here is the result using QGSP_BIC_HP and ShieldingLEND with B10 isotope target:

Hadr03 example output using QGSP_BIC_HP:
The run is 100000 neutron of 1e-05 eV through 10 m of B10 (density: 2.46 g/cm3 )
Process calls frequency:
nCapture= 12 neutronInelastic= 99988
.
MeanFreePath: 344.06 nm ± 343.1 nm massic: 0.084638 mg/cm2
CrossSection: 29065 cm^-1 massic: 11815 cm2/g
crossSection per atom: 1.9645e+05 barn
.
Verification: crossSections from G4HadronicProcessStore:
nCapture= 1.5336 cm2/g 25.499 barn
neutronInelastic= 11785 cm2/g 1.9595e+05 barn
total= 11787 cm2/g 1.9598e+05 barn
.
List of nuclear reactions:
neutron + B10 → N gamma or e- + B11: 12 Q = 13.066 MeV
neutron + B10 → N gamma or e- + alpha + Li7: 93670 Q = 2.7904 MeV
neutron + B10 → alpha + Li7: 6318 Q = 2.7901 MeV
number of gamma or e- (ic): N = 1 → 5
.
List of generated particles:
B11: 12 Emean = 3.7641 keV ( 1.3211 keV → 6.4208 keV)
Li7: 99988 Emean = 851 keV ( 831.89 keV → 1.0143 MeV)
alpha: 99988 Emean = 1.4916 MeV ( 1.4716 MeV → 1.7772 MeV)
gamma: 169268 Emean = 265.43 keV ( 1.0001 keV → 11.456 MeV)
Momentum balance: Pmean = 28.836 keV ( 853.5 eV → 82.715 keV)`

Hadr03 example output using ShieldingLEND:
The run is 100000 neutron of 1e-05 eV through 10 m of B10 (density: 2.46 g/cm3 )
Process calls frequency:
nCapture= 18 neutronInelastic= 99982
.
MeanFreePath: 350.66 nm ± 349.85 nm massic: 0.086262 mg/cm2
CrossSection: 28518 cm^-1 massic: 11593 cm2/g
crossSection per atom: 1.9275e+05 barn
.
Verification: crossSections from G4HadronicProcessStore:
nCapture= 1.5126 cm2/g 25.15 barn
neutronInelastic= 11624 cm2/g 1.9328e+05 barn
total= 11626 cm2/g 1.933e+05 barn
.
List of nuclear reactions:
neutron + B10 → N gamma or e- + B11: 18 Q = 13.604 MeV
neutron + B10 → N gamma or e- + alpha + Li7: 93810 Q = 5.3463 MeV
neutron + B10 → alpha + Li7: 6172 Q = 5.3461 MeV
number of gamma or e- (ic): N = 1 → 4
.
List of generated particles:
B11: 18 Emean = 4.7141 keV ( 260.61 eV → 14.003 keV)
Li7: 99982 Emean = 1.7796 MeV ( 1.7688 MeV → 1.9424 MeV)
alpha: 99982 Emean = 3.1186 MeV ( 3.0998 MeV → 3.4038 MeV)
gamma: 93854 Emean = 479.99 keV ( 477.6 keV → 11.456 MeV)
Momentum balance: Pmean = 448.04 keV ( 0.03038 eV → 477.6 keV)

We can see from the cross section summary and number of alpha+Li7 particles produced that the cross section data of G4NDL and LEND are the same. Just the issue of incorrect Q and kinetic energy in the output products.

I have found the cause of the issue by looking at the XML files in the LEND data folders. I found in LEND_GND1.3_ENDF.BVII.1/neutrons/005_B_010/005_B_010.000.xml at approximately line 400 under “Program DICTIN”, it has the definition of the particles with their mass and energy states:

However in the GIDI source files of Geant4, located under source/particles/management/src/MCGIDI_mass.cc,
it has the atomic masses of the nuclei stated, but with different values for alpha and Li7:

We see that the LEND XML data files has alpha=4.0015 amu and Li7=7.0143 amu, while GIDI has them stated again as alpha=4.0026 amu and Li7=7.0160 amu. After I edited the LEND XML file for B10 to have the same particle mass as GIDI, LEND started behaving as expected, with the same result as BIC_HP and literature:

Hadr03 example output using ShieldingLEND AFTER changing LEND XML data file:
The run is 100000 neutron of 1e-05 eV through 10 m of B10 (density: 2.46 g/cm3 )
Process calls frequency:
nCapture= 18 neutronInelastic= 99982
.
MeanFreePath: 350.66 nm ± 349.85 nm massic: 0.086262 mg/cm2
CrossSection: 28518 cm^-1 massic: 11593 cm2/g
crossSection per atom: 1.9275e+05 barn
.
Verification: crossSections from G4HadronicProcessStore:
nCapture= 1.5126 cm2/g 25.15 barn
neutronInelastic= 11624 cm2/g 1.9328e+05 barn
total= 11626 cm2/g 1.933e+05 barn
.
List of nuclear reactions:
neutron + B10 → N gamma or e- + B11: 18 Q = 13.604 MeV
neutron + B10 → N gamma or e- + alpha + Li7: 93810 Q = 2.7896 MeV
neutron + B10 → alpha + Li7: 6172 Q = 2.7897 MeV
number of gamma or e- (ic): N = 1 → 4
.
List of generated particles:
B11: 18 Emean = 4.7141 keV ( 260.61 eV → 14.003 keV)
Li7: 99982 Emean = 850.62 keV ( 839.9 keV → 1.0134 MeV)
alpha: 99982 Emean = 1.4909 MeV ( 1.4721 MeV → 1.7762 MeV)
gamma: 93854 Emean = 479.99 keV ( 477.6 keV → 11.456 MeV)
Momentum balance: Pmean = 448.04 keV ( 0.03036 eV → 477.6 keV)

This fix also corrects other two body output reactions such as Li6(n,a)H3.

Discussion:
The GND1.3 data file used by LEND (the XML file) is the same as one that can be downloaded from the ENDF website for the desired reactions. It also has the different mass definitions for alpha and Li7. If an older evaluation is used, such as ENDF6 the GND file has the same particle masses as GIDI_mass.

I believe that there may have been some slight correction that was made to the GND file to apply a correction, but this hasn’t been realised in the LEND model of Geant4.

My guess is that the GIDI processor for the LEND data files was originally made for ENDF6 and wasn’t updated after some ‘correction’ was made to the GND data file for ENDF7. So when the neutronInelastic process happens, the energy conservation cannot be calculated properly due to different definitions of the particle masses.

What does everyone else think?

This is a major issue that needs to be addressed by the LEND development team. (I have submitted bug report 2300 but haven’t heard back yet).

[siyuan1979]

i always can not use lend data files, can you upload your main.c code and environment variable about lend data files? thanks!

[jamesv]

Hello, I am using the Geant4 example ‘Hadr03’, which you can find in the source code examples:

examples/extended/hadronic/Hadr03

You will then need to edit the PhysicsList.cc to use ‘ShieldingLEND’ elastic and inelastic physics instead of the HP physics.

In PhysicsList.cc, add to the start of file:
#include “G4HadronPhysicsShieldingLEND.hh”
#include “G4HadronElasticPhysicsLEND.hh”
.
In PhysicsList::PhysicsList():G4VModularPhysicsList(), add:
RegisterPhysics( new G4HadronElasticPhysicsLEND(verb,“ENDF/BVII.1”));
RegisterPhysics( new G4HadronPhysicsShieldingLEND(verb));
and comment out:
//RegisterPhysics( new G4HadronElasticPhysicsHP(verb));
//RegisterPhysics( new G4HadronPhysicsQGSP_BIC_HP(verb));

You will need to export the path of where the LEND data is by typing into your terminal window. For example if you downloaded the data and put it with the other datasets (this is one line):

export G4LENDDATA=~/geant4.10.05.p01-install/share/Geant4-10.5.1/data/LEND_GND1.3_ENDF.BVII.1

Then you can run the example using neutron.mac:

./Hadr03 neutron.mac

You will know if the LEND data is working properly by looking at the physics process output for neutron. If you see anything about HP, then check the PhysicsList.cc file again:

Hadronic Processes for neutron
Process: hadElastic
Model: hElasticCHIPS: 19.5 MeV —> 100 TeV
Model: LENDElastic: 0 eV —> 20 MeV
Cr_sctns: LENDElasitcCrossSection: 0 eV —> 20 MeV
Cr_sctns: G4NeutronElasticXS: 0 eV —> 100 TeV
.
Process: neutronInelastic
Model: FTFP: 9.5 GeV —> 100 TeV
Model: BertiniCascade: 19.9 MeV —> 9.9 GeV
Model: LENDInelastic: 0 eV —> 20 MeV
Cr_sctns: LENDInelasitcCrossSection: 0 eV —> 20 MeV
Cr_sctns: JENDLHEInelasticCrossSection: 0 eV —> 100 TeV
Cr_sctns: Barashenkov-Glauber: 0 eV —> 100 TeV
.
Process: nCapture
Model: LENDCapture: 0 eV —> 20 MeV
Model: nRadCapture: 19.9 MeV —> 100 TeV
Cr_sctns: LENDCaptureCrossSection: 0 eV —> 20 MeV
Cr_sctns: G4NeutronCaptureXS: 0 eV —> 100 TeV
.
Process: nFission
Model: LENDFission: 0 eV —> 20 MeV
Model: G4LFission: 19.9 MeV —> 100 TeV
Cr_sctns: LENDFissionCrossSection: 0 eV —> 20 MeV
Cr_sctns: GheishaFissionXS: 0 eV —> 100 TeV

Obviously you will see the incorrect output energies for the BNCT reaction as I described but you can fix it by following the steps I said in the previous message. Depending on your application, you might not see any errors for your reactions if they are simple neutron interactions (not two+ body).

[siyuan1979]

when i use code you provide, the progran can not Initialize kernel,
can you upload your code? thanks!