Unexpected high counts in the energy distribution of alpha particles from Am‑241 decay

Hsiangfan_Mu · July 15, 2025, 11:04am

Dear All,
I’m simulating the energy distribution of particles emitted from radioactive decay of an Am-241 volume source and encountering an unexpected issue. At 0.025 MeV, I’m getting 2,000,179 counts, which is significantly higher than my run of 1e6 primary particles. This seems physically unrealistic.

My SteppingAction code is designed to record alpha particles that exit the AmSource_phys volume and enter the scintillator_phys volume:

if (step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName() == "AmSource_phys") {
    if (step->GetTrack()->GetParentID() == 1 && step->GetTrack()->GetParticleDefinition() == G4Alpha::Definition()) {
        if (step->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetName() == "scintillator_phys") {
            if (step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) {
                G4double kin = step->GetPostStepPoint()->GetKineticEnergy();
                analysisManager->FillNtupleDColumn(0, 2, kin / MeV);
                analysisManager->AddNtupleRow(0);
            }
        }
    }
}

My PhysicsList includes:

PhysicsList::PhysicsList() {
    RegisterPhysics(new G4EmStandardPhysics_option4());
    RegisterPhysics(new G4OpticalPhysics());
    RegisterPhysics(new G4DecayPhysics());
    RegisterPhysics(new G4RadioactiveDecayPhysics());
}

Is there an issue with my boundary condition check? Any insights or suggestions would be greatly appreciated!

Geant4 Version:11.2.2

mkelsey · July 15, 2025, 2:01pm

Did you configure the RadioactiveDecay process to only decay the Am-241 nucleus? Generally, Geant4 follows the whole decay chain to the end, so you’ll see not just the initial decay, but all of the subsequent daughters as well.

Hsiangfan_Mu · July 16, 2025, 1:04am

Hi mkelsey,

Yes, I had already restricted the RadioactiveDecay process to only decay the Am-241 nucleus, using

/process/had/rdm/nucleusLimits 241 241 95 95

Here’s my GPS macro:

/gps/particle ion
/gps/ion 95 241
/process/had/rdm/nucleusLimits 241 241 95 95
/gps/energy 0
/gps/pos/type   Volume
/gps/pos/shape  Para
/gps/pos/halfx  0.5 cm
/gps/pos/halfy  0.5 cm
/gps/pos/halfz  0.5 um
/gps/pos/confine AmSource_phys
/gps/pos/centre 0 0 0 cm
/gps/ang/type iso

And the geometry definition is:

G4Box *AmSource_solid =
    new G4Box("AmSource_solid", 0.5 * cm, 0.5 * cm, 0.5 * um);
G4LogicalVolume *AmSource_logical =
    new G4LogicalVolume(AmSource_solid, AmO2, "AmSource_logical");
new G4PVPlacement(nullptr, G4ThreeVector(), AmSource_logical, "AmSource_phys",
                  world_logical, false, 0, checkOverlaps);

G4Box *scintillator_solid =
    new G4Box("scintillator_solid", 0.5 * cm, 0.5 * cm, 0.5 * um);
G4LogicalVolume *scintillator_logical =
    new G4LogicalVolume(scintillator_solid, ZnS_Cu, "scintillator_logical");
new G4PVPlacement(nullptr, G4ThreeVector(0, 0, 1 * um), scintillator_logical,
                  "scintillator_phys", world_logical, false, 0,
                  checkOverlaps);

Since daughter decays should be suppressed, I’m puzzled about where this excess of low-energy alpha particles is coming from.

mkelsey · July 16, 2025, 3:00am

I agree with your puzzlement. That said, I notice that you’ve defined your source volume as actually made of AmO2; is that right? Something I don’t know is what the reaction Am-241 + alpha will do. Will it be inelastic, triggerring emission of another alpha?

… Except that your physics list doesn’t include any of the hadronic processes, so you shouldn’t be seeing secondary interactions.

jrellin · July 16, 2025, 4:30am

Simplify all of your nested-ifs (which makes it hard to debug) with something like:

if (step->GetTrack()->GetParentID() != 1){
  // Your particle gun only shoots alpha, just return
  return
}

if (step->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetName() == "scintillator_phys") && (step->GetPreStepPoint()->GetStepStatus() == fGeomBoundary){
  // This is an alpha crossing a boundary entering into the 'scintillator'
  G4double kin = step->GetPostStepPoint()->GetKineticEnergy();
  analysisManager->FillNtupleDColumn(0, 2, kin / MeV);
  analysisManager->AddNtupleRow(0);
}

What does this look like in the UI? I wonder if you aren’t also getting undefined floating point errors at the boundary. You should make one box of one of the materials that covers both volumes and then make the second a daughter volume of the first.

Hsiangfan_Mu · July 16, 2025, 7:01am

Yes, the source material is AmO2. In principle, Am-241 can undergo α‑induced reactions such as 241Am(α, 2n)243Bk and 241Am(α, 3n)242Bk, α‑particles from natural Am-241 decay do not have enough energy to open those channels.

Hsiangfan_Mu · July 16, 2025, 7:02am

Here’s what the UI looks like:

The blue region—the world volume, defined as air.
The green region—the AmO2 source and the ZnS:Cu scintillator.

Hsiangfan_Mu · July 17, 2025, 1:07am

Here’s an update on the issue:
After switching the output from ntuple fills to a histogram (using FillH1() instead of FillNtupleDColumn()), the spectrum looks perfectly normal.
This strongly suggests the problem lies in the ntuple writing step, but I’m still not sure what’s going wrong there.

if (step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetName() ==
      "AmSource_phys") {
    if (step->GetTrack()->GetParentID() == 1 &&
        step->GetTrack()->GetParticleDefinition() == G4Alpha::Definition()) {
      if (step->GetPostStepPoint()
              ->GetTouchableHandle()
              ->GetVolume()
              ->GetName() == "scintillator_phys") {
        if (step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) {
          G4double kin = step->GetPostStepPoint()->GetKineticEnergy();
          //          analysisManager->FillNtupleDColumn(0, 2, kin / MeV);
          //          analysisManager->AddNtupleRow(0);
          analysisManager->FillH1(10, kin);
        }
      }
    }
  }

jrellin · July 17, 2025, 5:13am

Interesting…

What do you want your ntuple to look like? Right now, without context, this looks like you are filling the second (“2”) column with your energy values. But you are not also filling out the first (or any other?) column of values at the same time before you add the row?

Is this meant to be a single 1D column of values? If so, you need to get rid of the 2. You might be running into a problem where it doesn’t know what to do with the 2, it becomes a null value that is then histogrammed as 0.

If there are other columns somewhere you will run into an issue where you don’t have values for each column. NTuples are rectangular, not “ragged” arrays. That could also cause 0s to be put in a lot of places to “fill” the space. But I am not sure.