In the Geant4 simulation of the D-Be reaction to calculate neutron yields, the simulation results differ from the reference values by two orders of magnitude

_Geant4 Version:_10.5.1
_Operating System:_Debian 12
_Compiler/Version:_gcc 12.2
_CMake Version:_3.25.1

— I am using Geant4 to calculate the neutron yield from the D-Be reaction, with an incident particle number of 6.25×1086.25 \times 10^86.25×108. In the low-energy region (500 keV), no neutrons are produced. At an incident deuteron energy of 2 MeV, the neutron yield is two orders of magnitude smaller than the experimental values. Below is my physics list and geometry/material code. Could you suggest how I should modify it?

// Hadron Elastic scattering
fHadronElastic = new HadronElasticPhysicsHP(verb);
RegisterPhysics(fHadronElastic);

// Hadron Inelastic Physics
fHadronInelastic = new G4HadronPhysicsQGSP_BIC_AllHP(verb);
RegisterPhysics(fHadronInelastic);

// Ion Elastic Physics
fIonElastic = new G4IonElasticPhysics(verb);
RegisterPhysics(fIonElastic);

// Ion Inelastic Physics
fIonInelastic = new G4IonPhysicsPHP( verb );
RegisterPhysics(fIonInelastic);

// stopping Particles
RegisterPhysics( new G4StoppingPhysics(verb));

// Gamma-Nuclear Physics
fGammaNuclear = new GammaNuclearPhysics(“gamma”);
RegisterPhysics(fGammaNuclear);

// EM physics
fElectromagnetic = new ElectromagneticPhysics();
RegisterPhysics(fElectromagnetic);

// Decay
fDecay = new G4DecayPhysics();
RegisterPhysics(fDecay);

// Radioactive decay
fRadioactiveDecay = new RadioactiveDecayPhysics();
RegisterPhysics(fRadioactiveDecay);

G4double density = 1.848 * g / cm3;

G4Isotope* Be = new G4Isotope( "Be9", 4, 9, 9.0122 * g/mole );
G4Element* elBe = new G4Element("Beryllium", "Be", 1); 
elBe->AddIsotope(Be, 100 * perCent);
G4Material* BeMaterial = new G4Material( "BeTarget", density, 1 );
BeMaterial->AddElement( elBe, 100 * perCent );
G4Material* world_mat = nist->FindOrBuildMaterial( "G4_Galactic" );
G4Box* sworld = new G4Box( "World", 0.5*world_x, 0.5*world_y, 0.5*world_z );
G4LogicalVolume* lworld = new G4LogicalVolume( sworld, world_mat, "world" );
G4VPhysicalVolume* pworld = new G4PVPlacement( 0, G4ThreeVector(), lworld, "world", 0, false, 0 );
G4double targetDiameter = 2 * cm;
G4double targetThickness = 0.1 * cm;
G4Tubs* solidTarget = new G4Tubs( "BeTarget", 0., targetDiameter / 2, targetThickness / 2, 0., 360. * deg );
G4LogicalVolume* logicTarget = new G4LogicalVolume( solidTarget, BeMaterial, "Target" );
G4ThreeVector targetPosition( 0., 0., 0. );
new G4PVPlacement( nullptr, targetPosition, logicTarget, "Target", lworld, false, 0, true );

Hello,

Here are some thoughts on your problem:

  • you are already using the NIST-Manager for defining world_mat;
    you can use it to define Beryllium too - your definition of Beryllium looks right, but NIST manager will defintely do it right! This should do the trick:

G4Material* Be = nist->FindOrBuildMaterial( “G4_Be” );

  • Physics list looks okay too - G4HadronPhysicsQGSP_BIC_AllHP should, be the way to go for Be+d. In newer Geant4 versions, this physics list uses the G4TENDL datasets. Your Geant4 version (10.5.1) is 5 years old and there isn’t even a G4TENDL dataset to download on the download page for your version . I strongly recommend updating since these kind of reactions (deuterons on light ion target) seem to be very sensitive to Geant4 version and Dataset version as you can read about in this post.
  • How do you track/count the neutrons? Maybe it’s a counting problem and the right amount of neutrons are produced but you dont detect them.
  • What experimental values are you referring to? Is there an experiment in a paper you are trying to reproduce?

Regards
Max

1 Like

Thank you for your advice.
My current work is to simulate the D-Be reaction using Geant4 to obtain the yield, angular distribution, and energy spectrum data of the neutrons produced. Setting the Be target to 1mm thick
As per your suggestion, I am re-running the simulation under Geant4 version 11.2 and redefining the material in the way you suggested to define the material
After the simulation, I found that when the energy of the incident deuterium ion is 2MeV, the number of neutrons produced matches the actual measurements provided in the paper, but when the energy of the incident deuterium ion is 500KeV or less, there is no neutron production (the simulated number of particles is 6.25*10^8) I presume that there is no current data on the cross section of the reaction of D-Be at low energy conditions in Geant4.
Thanks again for the advice you have provided!

love&peace
Yu

1 Like

Hello again,

I am re-running the simulation under Geant4 version 11.2 (and) I found that when the energy of the incident deuterium ion is 2MeV, the number of neutrons produced matches the actual measurements provided in the paper

Thanks for sharing this information! I am glad to here the version difference was the problem and updating did fix the production rate!

To fix the 500 Kev-Problem you could try to use a comprehensive reference physics list instead of building everything by yourself - maybe you are missing the physics list that does the neutron production at 500 kev. See the physics list guide and the reference physics list documentation. The obvious choice would be to use the QGSP_BIC_AllHP physics list.

I suppose your code where you initialize your physics list looks somehow like this:
runManager->SetUserInitialization(new PhysicsList);

Just change those lines to
G4VModularPhysicsList* physicsList = new QGSP_BIC_AllHP;
runManager->SetUserInitialization(physicsList);

and dont forget to
include "QGSP_BIC_AllHP.hh"
at the start!

Regards
Max

Hello. Sorry to bother you again.
I recompiled the run in the manner you provided and found no change in the results at incident deuterium beam energies below 500 KeV:
List of generated particles (with meanLife ! = 0) : e-: 3 Emean =
e-: 3 Emean = 586.9 keV ( 379.22 keV → 951.79 keV) stable
gamma: 183 Emean = 1.2738 MeV ( 99.553 keV → 3.6501 MeV) stable
And I found that deuterium beam energies below 1.5 MeV have no neutron production when simulating the D-Be reaction, but above 2 MeV, the results are normal. At the moment I don’t know how to modify

love&peace
Yu

Hello again,

List of generated particles (with meanLife ! = 0) : e-: 3 Emean =
e-: 3 Emean = 586.9 keV ( 379.22 keV → 951.79 keV) stable
gamma: 183 Emean = 1.2738 MeV ( 99.553 keV → 3.6501 MeV) stable

How many events did you simulate? Earlier, you said you’ve simulated 6.25*10^8 events but the numbers of generated particles here looks like you simulated just 1000 events or so.
If so, can you run it with 1 million events?

Otherwise i dont really have any idea what else you can try for now.

Regards
Max

Hello
In the previous simulations, in order to reduce the running time and obtain the neutron yield of deuterium beams with different energies, I set the number of simulated particles to 6.25*10^7. The results are as shown in the previous message. Under the energy condition of less than 1.5MeV, only gamma and electrons are produced.

Thanks
Yu

Hello,

if you use QBBC physics would results be the same or different?

VI

particlePHP uses the data you can find in G4DEUTERONHPDATA. As you can see in the README file, Be data are taken from JENDL/DEU-2020. You can edit the 4_9_Berylium.z file. The first step is to convert it to a text file, without compression
zlib-flate -uncompress < Inelastic/F02/4_9_Berylium.z > 4_9_Berylium

The format of this file is explained in the Geant4 Physics Reference Manual, or directy at https://geant4-userdoc.web.cern.ch/ContributionFromUsers/UsefulNotes/particleHPFormat.pdf