Geant4 Version: 4-11-03 Patch 02
Operating System: Mac OS Sequoia 15.7.2
Compiler/Version: Installed via Anaconda 2.4.0 – CMake reports compiler is AppleClang 17.0.0.17000319
CMake Version: 4.2.0
Hi,
I’ve been using ESA’s GRAS for several years to perform trapped radiation analysis for spacecraft design, but am trying to use G4 directly for some work on radioisotope radiation emission characteristics. I’m simulating (α,n) neutron production from Am-241 in AmO₂ , basing my code on a modified version of the rdecay02 example which comes with Geant4.
My simulation produces ~1,200 neutrons per 100 million Am-241 decays, which translates to a neutron yield of ~1.4 × 10⁶ n/s/g of AmO₂. However, the literature value for AmO₂ is 2.7 × 10³ n/s/g (this value is consistently quoted - e.g. Lees & Lindley, Annals of Nuclear Energy, Vol.5, p133, 1978). So my result is ~500× too high.
My geometry is as follows:
- Target is a 1 cm radius solid sphere of AmO2 (density 11.68 g/cm³, mass ~48.9 g)
- Detector is a spherical shell of Ge, inner radius 3cm, outer radius 23cm (20 cm thick).
- A 2cm radial gap (vacuum/G4_Galactic) sits between the Ge shell and the AmO2
- World volume is a 146 x 146 x 146 cm box filled with G4_Galactic
- AmO2 and Ge volumes are concentric spheres centred at the origin.
My materials definition is as follows:
// O-16, O-17, O-18 isotopes
G4Isotope* O16 = new G4Isotope(“O16”, 8, 16, 15.9949*g/mole);
G4Isotope* O17 = new G4Isotope(“O17”, 8, 17, 16.9991*g/mole);
G4Isotope* O18 = new G4Isotope(“O18”, 8, 18, 17.9992*g/mole);
G4Element* elO = new G4Element(“Oxygen_natural”, “O”, 3);
elO->AddIsotope(O16, 99.757*perCent);
elO->AddIsotope(O17, 0.038*perCent);
elO->AddIsotope(O18, 0.205*perCent);
// Am-241 isotope
G4Isotope* Am241 = new G4Isotope(“Am241”, 95, 241, 241.0568*g/mole);
G4Element* elAm = new G4Element(“Americium241”, “Am”, 1);
elAm->AddIsotope(Am241, 100.0*perCent);
// AmO2 material
fTargetMater = new G4Material(“AmO2”, 11.68*g/cm3, 2);
fTargetMater->AddElement(elAm, 1);
fTargetMater->AddElement(elO, 2);
I’m using a General Particle Source generating 100 million Am-241 ions at rest, uniformly distributed throughout the AmO₂ sphere volume, undergoing radioactive decay (and confirmed this is producing ~5.5 MeV alphas):
/gps/particle ion
/gps/ion 95 241 0 0
/gps/energy 0 eV
/gps/pos/type Volume
/gps/pos/shape Sphere
/gps/pos/centre 0 0 0 mm
/gps/pos/radius 1 cm
/gps/pos/confine Target
My nuclear data libraries are as follows - the G4 output at run time confirms these are being read successfully:
G4PARTICLETENDLDATA=/path/to/G4TENDL1.4
G4NEUTRONHPDATA=/path/to/NDL4.7.1
G4PARTICLEXSDATA=/path/to/G4PARTICLEXS3.1
I’m using a custom physics list with G4ParticleHPInelastic to handle alpha inelastic reactions via TENDL data:
void PhysicsList::AddAlphaInelasticHP()
{
G4ParticleDefinition* alpha = G4Alpha::Alpha();
G4ProcessManager* pManager = alpha->GetProcessManager();
// Remove existing alphaInelastic process
G4ProcessVector* processList = pManager->GetProcessList();
for (size_t i = 0; i < processList->size(); ++i) {
if ((*processList)[i]->GetProcessName() == “alphaInelastic”) {
pManager->RemoveProcess((*processList)[i]);
break;
}
}
// Create new inelastic process with ParticleHP
G4HadronInelasticProcess* alphaInelastic =
new G4HadronInelasticProcess(“alphaInelastic”, G4Alpha::Definition());
// ParticleHP model for low energies (uses TENDL data)
G4ParticleHPInelastic* hpModel = new G4ParticleHPInelastic();
hpModel->SetMinEnergy(0.0 * MeV);
hpModel->SetMaxEnergy(20.0 * MeV);
alphaInelastic->RegisterMe(hpModel);
// Binary cascade for higher energies
G4BinaryLightIonReaction* blic = new G4BinaryLightIonReaction();
blic->SetMinEnergy(15.0 * MeV);
blic->SetMaxEnergy(6.0 * GeV);
alphaInelastic->RegisterMe(blic);
// Glauber-Gribov cross-sections
G4ComponentGGNuclNuclXsc* ggNuclNuclXsc = new G4ComponentGGNuclNuclXsc();
G4CrossSectionInelastic* alphaNucleonXsc =
new G4CrossSectionInelastic(ggNuclNuclXsc);
alphaInelastic->AddDataSet(alphaNucleonXsc);
pManager->AddDiscreteProcess(alphaInelastic);
}
The base physics list uses the following models:
- G4EmStandardPhysics
- G4DecayPhysics
- BiasedRDPhysics (custom radioactive decay)
- G4HadronPhysicsQGSP_BIC
- G4HadronElasticPhysics
- G4IonPhysics
Finally my macro is as follows - I’m restricting the decay to AmO2 only - no daughter decays:
/control/verbose 0
/run/verbose 0
/event/verbose 0
/tracking/verbose 0
#
# Set the target/detector geometry
# Note: For sphere, “Length” is interpreted as diameter (radius = length/2)
/rdecay02/det/setTargetRadius 1 cm
/rdecay02/det/setTargetLength 2 cm
/rdecay02/det/setDetectorThickness 20 cm
/rdecay02/det/setDetectorLength 40 cm
#
# Randomise the seed FIRST
/random/setSeeds -1 -1
/run/initialize
#
# Set a very high time threshold to allow all decays to happen
/process/had/rdm/thresholdForVeryLongDecayTime 1.0e+60 year
#
/process/list
#
# rdm is applied only to the target volume
/process/had/rdm/noVolumes
/process/had/rdm/selectVolume Target
#
# Only Am-241 decays (no daughters)
/process/had/rdm/nucleusLimits 241 241 95 95
# ============================================
# GPS Configuration - Volumetric Am-241 source
# ============================================
#
# Define the source particle (Am-241 ion at rest)
/gps/particle ion
/gps/ion 95 241 0 0
/gps/energy 0 eV
#
# Volumetric distribution - sphere
/gps/pos/type Volume
/gps/pos/shape Sphere
/gps/pos/centre 0 0 0 mm
/gps/pos/radius 1 cm
#
# Isotropic angular distribution (for decay products)
/gps/ang/type iso
#
# Confine source to Target logical volume
# This ensures particles only start inside the AmO2 sphere
/gps/pos/confine Target
#
# ============================================
# Analysis and Run Settings
# ============================================
/analysis/h1/set 6 100 0. 15. MeV
#
# Complete chain decays
/analysis/setFileName Am241c_volumetric
/run/printProgress 1000000
/run/beamOn 100000000
My results are as follows:
With TENDL-1.4 +ParticleHP:
- Neutrons generated: ~1,170 (statistical variation across runs)
- Alpha inelastic events: ~3,650
- Calculated yield: 1.4 × 10⁶ n/s/g of AmO₂
With standard G4IonPhysics (no TENDL):
- Neutrons generated: 0-1
- The standard physics completely misses (α,n) reactions on oxygen
With QGSP_BIC_AllHP:
- Neutrons generated: 1
- Also fails to produce (α,n) reactions
During Execution I receive a sequence of warnings like this:
G4IonTable::GetIon() : illegal atomic number/mass Z =-2 A = -3 E = 0 *** G4Exception : hadr01 issued by :G4ParticleHPFinalState:adjust_final_state G4ParticleHPFinalState oneMoreSec_pd Z=95 - 97 A=242 - 245 projectile neutron No adjustment will be done!
I’ve been looking at this for weeks, trying to track down the error or omission, and would appreciate any advice. Some questions:
- Is TENDL-1.4 known to overestimate (α,n) cross-sections for O-17 and O-18?
- Is there a better way to implement alpha inelastic physics for light target nuclei in Geant4 11.3?
- Should I be using a different physics list or data library? (TENDL-2021? Different approach entirely?)
- Could the warnings about illegal Z/A values be causing incorrect neutron counting or double-counting?
- Is there a validated Geant4 example for (α,n) source calculations I should reference?
Any thoughts or advice welcome. Thanks for reading.