QGSP_BIC_AllHP:alphaInelastic on 9Be emits a duplicated half-energy gamma instead of a single discrete line

gitonganjeri · July 31, 2026, 11:42am

Geant4 Version: 11.4.1
Operating System: Rocky Linux 10.1
Compiler/Version: gcc 15.2.0
CMake Version: 4.3.2

Physics list: QGSP_BIC_AllHP

Reaction: ⍺ + ⁹ Be → n + ¹² C (process name alphaInelastic), alpha energies 5.39-5.49 MeV (²⁴¹ Am decay lines

Summary

While simulating an Am-241/Be neutron source from first principles (sampling Am-241’s alpha decay branches directly as primaries and letting the physics list handle the (⍺,n) reaction during transport, with JENDLTENDL data library for alpha cross-sections), I found that whenever the reaction leaves the residual ¹² C nucleus in an excited state, alphaInelastic returns two gamma secondaries with numerically identical energy, rather than one gamma at the correct transition energy. This happens consistently for both of ¹² C’s low-lying excited states, and I have not observed a single case of lone, correctly-valued gamma from this process.

Observed behaviour

For the first excited state (true transition energy 4.4389 MeV):

— alphaInelastic step: 4 secondaries —

neutron  E=5.60587 MeV

C12  E=0.771438 MeV

gamma  E=2.21873 MeV

gamma  E=2.21873 MeV

2.21873 x 2 = 4.43746 - matches the real 4.4389 MeV ¹² C²⁺ line to within ~1 keV.

For the second excited state:

— alphaInelastic step: 4 secondaries —

neutron  E=1.08597 MeV

C12  E=1.49862 MeV

gamma  E=3.82636 MeV

gamma  E=3.82636 MeV

3.82636 x 2 = 7.65272 MeV - closely matching the 7.654 MeV ¹² C ⁰⁺ (Hoyle State) to within a few keV

Across a 10M-event run, a histogram of gamma multiplicity per alphaInelastic step shows only 0 or 2, never 1 - i.e. this is not an occasional artifact, it appears to be systematic whenever a discrete gamma should be emitted at all.

Diagnostic method

Confirmed via track→GetCreatorProcess()→GetProcessName() on the gamma tracks themselves (returns alphaInelastic), track→GetParentID() (matches primary alpha’s track ID), and creation volume (matches the source volume) - so this is not a downstream/unrelated process being misattributed.
Confirmed the duplication is exact: gamma1_E - gamma2_E shows a hard spike at zero across all flagged events, with no spread.
Ruled out double-counting in own analysis code (separate boundary-crossing tally bug, already found and fixed independently) by checking directly at the alphaInelastic steps’s G4Step::GetSecondaryInCurrentStep() list, before any of the downstream code touches the tracks.

Expected behaviour:

The 4.439 MeV level should emit a single gamma at the full transition energy (it decays radiatively essentially 100% of the time).
The Hoyle state (7.654 MeV, 0⁺ ) should predominantly (~99.96%, per published branching ratio measurements, ~4x10^-4 radiative fraction) undergo alpha breakup (¹²C^* → ⍺ + ⁸Be-> 3⍺), producing no discrete gamma at all; the rare radiative branch should cascade via two different energies (≅ 3.214 MeV then 4.439 MeV), not two identical half-energy photons.

Currently, 100% of Hoyle-populating events in the simulation produce a (duplicated, wrong-energy) gamma pair, and none show alpha-breakup secondaries.

Minimal Reproducer

A SteppingAction::UserSteppingAction() hook like the following, attached to a run firing bare alpha primaries (5.4-5.5MeV) into a Be target, reproduces this on our end:

if (track->GetTrackID() == 1 && track->GetDefinition() == G4Alpha::AlphaDefinition()) {
    const G4VProcess* proc = step->GetPostStepPoint()->GetProcessDefinedStep();
    if (proc && proc->GetProcessName() == "alphaInelastic") {
        for (auto sec : *step->GetSecondaryInCurrentStep()) {
            G4cout << sec->GetDefinition()->GetParticleName()
                   << " E=" << sec->GetKineticEnergy()/MeV << " MeV" << G4endl;
        }
    }
}

Happy to share the full generator/geometry setup if useful for reproducing this independently.

Note

This is posted for awareness/potential upstream investigation rather than a request for a workaround - we’ve patched around it on our end (collapsing the duplicate gamma for the 4.439 MeV case, and suppressing it in favour of the dominant alpha-breakup outcome for the Hoyle case). Flagging in case it’s useful for future releases.

maire1 · August 7, 2026, 4:51pm

Here, a run of example Hadr03.

we do not observe those splitted gamma

njeri.mac.txt (381 Bytes)

njeri.out.txt (2.9 KB)

gitonganjeri · August 25, 2026, 2:51pm

Hi @maire1

Thanks for looking into this and testing it with the Hadr03 example.

In my simulation, I am using the standard pre-packaged QGSP_BIC_AllHP physics list, rather than a custom modular list. I noticed that the example registers multiple alternative physics lists, and with your sample macros, tried to see if the same pattern I observe in my simulation appears in the example with the exclusive use of QGSP_BIC_AllHP.

With the only changes I’ve made to the example’s main.cc being:

#include "QGSP_BIC_AllHP.hh"
// #include "PhysicsList.hh"
...
// in main()
G4VModularPhysicsList* phys = new QGSP_BIC_AllHP;
runManager->SetUserInitialization(phys);

I run it with two scenarios:

Scenario 1: The environment variables calls commented out - Best replicates my simulation
Scenario 2: The environment variables as in the original example

comparison.out.txt (5.9 KB) - Gives the extracted data from the ‘.out’ file for both scenarios

And the gamma plots I get are as shown. Scenario2 results in a blank plot. From scenario1, the split of the 4.44 MeV gamma would contaminate the neutron capture peak when introducing a hydrogenous material