Pulsed Height Tally in mixed neutron-gamma fields

[kosbor-personal]

Intro

Hi everyone!

Long time listener, first time caller. This is my first Pull Request, so I would appreciate any feedback on how I can improve this contribution and also potential future ones.

All the best,
Bor Kos

Problem Description

We’ve been working on verifying and validating OpenMC for well-logging applications. In this pull request, I wanted to address an issue we encountered when we started analysing so‑called pulsed‑neutron applications, where we are interested in gamma spectra from capture reactions from 14 MeV neutrons thermalized within a rock formation.

Long story short, these have been our observations regarding pulse height tallies (PHT) in gamma-only and coupled neutron–gamma fields:

• Correct behaviour for gamma‑only PHT cases (i.e. perfect match with MCNP and experimental results)
• Excellent agreement between OpenMC and MCNP for gamma spectra in coupled neutron–gamma simulations
• Significant deviations in pulse height tallies for coupled neutron–gamma simulations

This issue was shared on the forum:
https://openmc.discourse.group/t/pulse-height-tally-in-mixed-neutron-gamma-fields/6229

This led me to believe that there is an issue in how neutron‑generated photons are treated within the pulse height accounting logic.

Looking at the code, I believe photons produced by neutrons were incorrectly having their energy subtracted from pht_storage, leading to systematic undercounting in PHT results.

In pht_secondary_particles(), OpenMC subtracts energy from pht_storage for all secondary photons, assuming their parent’s energy had previously been added to the PHT. This assumption is valid for photon/electron/positron parents but invalid for neutron reactions, because:

• Neutron energy is never added to pht_storage
• Neutron‑generated photons therefore had energy subtracted without ever being added

This results in a net loss of recorded pulse height.

Fix

The fix ensures that energy is only subtracted for secondary particles whose parent particle is either a photon, electron, or positron. This is achieved by tracking the parent particle type for each secondary particle and only applying the energy subtraction when the parent is a photon/electron/positron.

Code Changes

1. include/openmc/particle_data.h

   • Added parent_particle to SourceSite
   • Added corresponding member variable and accessors to ParticleData

2. src/particle.cpp

   • Set parent_particle in create_secondary()
   • Propagate via from_source()
   • Apply whitelist logic in pht_secondary_particles()

3. openmc/lib/core.py

   • Added parent_particle to the Python _SourceSite ctypes structure

Testing / Validation

• Pytest suite

  • Failed tests matched failures in unmodified openmc-dev.
  • Suspected cause is compiler settings.
  • None of the failed tests are related to pulse height tallies.

• Gamma‑only PHT

  • Results remain unchanged.
  • Plots attached for two different sources:
    • Potassium (left)
    • Thorium (right)

• Coupled neutron–gamma simulations
  • Pulsed‑neutron tools tested in two environments:
    • Dolomite
    • Brine
  • Results match MCNP.
  • Lesser effect observed in brine, most likely due to high chlorine absorption.

[GuySten]

I think you can add a guard in this line (skip it for neutrons) instead of adding another bank attribute.
This should significantly simplify this PR.

[kosbor-personal]

That's a great idea - so maybe just do this:
if (this->type() == ParticleType::photon || this->type() == ParticleType::electron || this->type() == ParticleType::positron) { bank_second_E() += bank.E; }
whitelist particles instead of blacklisting neutrons for potential future proofing?

I see that the heating tallies also use bank_second_E() - will implementing this break something with them unintentionally?

[GuySten]

I didn't notice that the heating tallies are also using bank_second_E.
So it is likely that my idea does not work as is.

[kosbor-personal]

I am worried about another situation which might arise if we don’t track parent_paricles and we have a particle cascade situation. As far as I can tell, we use the last in, first out principle when banking particles in OpenMC, right? So, imagine this situation (G=gamma, E=electron):

1. Neutron creates G1, G2: bank = [G1, G2]
2. Neutron dies, pop G2: Current particle type = neutron. Correct parent for G2
3. G2 creates E2: bank = [G1, E2]
4. G2 dies, pop E2: Current particle type = photon. Correct parent for E2
5. E2 dies, pop G1: Current particle type = electron. Wrong – G1’s parent was a neutron. And its energy shouldn’t be subtracted from pht_storage.

As the bank can contain secondaries from different parents, I think we have to store the actual parent types.

[GuySten]

I think we can add the photon energy to pht_storage when sampling secondary photons. That way when we remove their energy afterwards it will be okay.
So we won't have to store more data.

[kosbor-personal]

Great suggestion. I’ve been looking at the code and I wanted to confirm what you meant before implementing.
So as far as I can understand the neutron created photons are sampled exclusively in sample_secondary_photons() in physics.cpp, while fluorescence, bremsstrahlung, and annihilation photons are created elsewhere (photon.cpp, bremsstrahlung.cpp). We have two different paths of creation so we can separate the two instances and how they contribute to the pht at the source.

So we could:
Pre-add the photon energy of neutron created photons to pht_storage specifically in sample_secondary_photons() (right after create_secondary) and thus keep the pht_secondary_particles() always subtracting on revival.
Two cases for secondary photons:

1. Neutron photon: pre-add E -> subtract E (revival) -> add E (absorption) = correct energy deposited
2. EM photon: no pre-add -> subtract E (revival) -> add E (absorption) = subtact/add cancel out resulting in the correct energy deposit

Is this the approach you had in mind? So something like this added to sample_secondary_photons around line 1217 is basically all that’s needed (following the same style of pramgming as in pht_collison_energy() and pht_secondary_particles()):

bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon);

// Pre-add energy to pht_storage for neutron-created photons
if (created_photon && !model::active_pulse_height_tallies.empty()) {
  auto it = std::find(model::pulse_height_cells.begin(),
    model::pulse_height_cells.end(), p.lowest_coord().cell());
  if (it != model::pulse_height_cells.end()) {
    int index = std::distance(model::pulse_height_cells.begin(), it);
    p.pht_storage()[index] += E;
  }
}

[GuySten]

You understood me correctly.

[kosbor-personal]

Implemented the suggested fix after testing (Mod. Src. Code. in Figure). Test show same results as previous implementation of fix (Mod. Src. Code. Old. in Figure).

[GuySten]

Looks good to me.
I will merge next week if no one objects.

[paulromano]

Thanks @kosbor-personal for the contribution and @GuySten for the simplification! Looks good to me as well. I also added a quick regression test to cover this particular case so that we don't break it in the future :)