Bug report: rare large TLEDoseActor score spikes in eight-thread OpenGATE 10.1.1 runs

[GitFuchs]

I have encountered a bug in the TLEDose Actor.

For a set of independent photon-energy-component simulations, TLEDoseActor produces rare, narrow upward spikes in energy deposition and much larger spikes in the associated per-event squared sums. The relative statistical uncertainty rises at the same depths. The co-scored classic DoseActor has no comparable spike at the example location.

The described bug appears on almost any high-statistic multi-threaded simulation in (at least) the 100-120 keV energy range.

I have ChatGPT assist me in writing a proper bug report. Please find it attached.

TLEproblem_mod.md

[BishopWolf]

Bug report — #1135

Rare large TLEDoseActor score spikes in multi-threaded runs

Verdict: confirmed as a multi-threaded data race in the TLE actor's photon-attenuation-coefficient lookup path. The exact races exist in 10.1.1, were fixed on master by commit a92a62bf1 "fix thread race issue", and that fix is in no release tag — the same “fixed on master, still broken in the shipped wheel” pattern as #1115. The mechanism reproduces every one of the reported signatures (upward-only spikes in edep, quadratically larger spikes in the squared sums, rising uncertainty, clean classic DoseActor).

DeepSeek v4.1 flash

---

1. Symptom decoded

Reported observation	Explanation in the code
Rare, narrow upward spikes in energy deposition (TLE only)	TLE edep is computed analytically from a tabulated mass-energy-absorption coefficient. A racy lookup returns the coefficient of a different energy (occasionally much larger), which scales edep directly.
Much larger spikes in the per-event squared sums	edep_squared stores (Σ over the event)². A single spiked event enters edep linearly but the squared image quadratically, so the outlier is amplified by its own factor.
Relative statistical uncertainty rises at the same depths	unc = sqrt(<E²>−<E>²)/<E> — one large E-outlier inflates the numerator (<E²>) far more than the mean, exactly at the spiked voxels.
Classic DoseActor has no comparable spike at the same location	The DoseActor scores the actual step->GetTotalEnergyDeposit() and never touches the attenuation tables, so it is structurally immune to this race.
Appears on high-statistic multi-threaded runs, in the 100–120 keV range	More threads ⇒ more interleavings; more events ⇒ more chances. That energy band is where μ_en/ρ varies most steeply (absorption-edge vicinity), so a lookup that returns a neighbouring energy's value is largest there.

---

2. Root cause #1 — shared mutable memoization cache in GateMuTable (the spike mechanism)

In 10.1.1, GateMuTable memoizes the last query in plain member variables (core/opengate_core/opengate_lib/GateMuTables.h, 10.1.1):

double GetMuEn(double energy);
double GetMuEnOverRho(double energy);      // non-const
double GetMu(double energy);
double GetMuOverRho(double energy);
...
double fLastEnergyMu;
double fLastEnergyMuEn;
double fLastMu;
double fLastMuEn;

and the getters are a classic single-consumer cache (GateMuTables.cpp, 10.1.1):

double GateMuTable::GetMuEnOverRho(double energy) {
  if (energy != fLastEnergyMuEn) {          // (1) test shared cache
    fLastEnergyMuEn = energy;               // (2) publish this energy
    energy = log(energy);
    int inf = 0, sup = fSize - 1;
    while (sup - inf > 1) { /* binary search */ }
    double e_inf = fEnergy[inf], e_sup = fEnergy[sup];
    if (energy > e_inf && energy < e_sup)
      fLastMuEn = exp(interpol(e_inf, energy, e_sup, fMuEn[inf], fMuEn[sup]));  // (3)
    else
      fLastMuEn = exp(fMuEn[inf]);                                            // (3)
  }
  return fLastMuEn;                         // (4) read shared cache
}

There is one GateMuTable per material, shared by all worker threads (see §3), and no lock. The interleaving that produces an outlier:

1. Thread A evaluates energy != fLastEnergyMuEn → true, sets fLastEnergyMuEn = E_A, starts the search.
2. Thread B sets fLastEnergyMuEn = E_B and writes fLastMuEn = v_B.
3. Thread A performs return fLastMuEn; → returns v_B, the coefficient of a different energy, for a step at E_A.

The returned value feeds straight into the TLE estimator (GateTLEDoseActor.cpp):

edep = weight * 0.1 * energy * mu_en_over_rho * step_length * density / (CLHEP::g / CLHEP::cm3);

so a wrong mu_en_over_rho is a wrong edep — occasionally by a large factor near an edge, hence a narrow upward spike. The guard at (1) can also be satisfied spuriously (another thread wrote the same value), so a thread may return a cached value computed for an entirely different material/energy.

---

3. Root cause #2 — unprotected singleton, lazy shared initialization, and a second shared “last-used” cache

The handler that owns the tables is a process-wide singleton per (database, emax), and in 10.1.1 its creation is unsynchronized (GateMaterialMuHandler.cpp, 10.1.1):

std::shared_ptr<GateMaterialMuHandler>
GateMaterialMuHandler::GetInstance(std::string database, double energy_max) {
  auto key = std::make_tuple(database, energy_max);
  if (fInstances.find(key) == fInstances.end()) {     // unsynchronized std::map read + mutate
    fInstances[key] = std::make_shared<GateMaterialMuHandler>();
    fInstances[key]->SetDatabaseName(database);
    fInstances[key]->SetEMax(energy_max);
  }
  return fInstances[key];                             // operator[] can insert
}

First use is triggered lazily from whichever worker thread needs it first, and initialization is heavy and mutates the handler (GateMaterialMuHandler.cpp, 10.1.1):

void GateMaterialMuHandler::CheckLastCall(const G4MaterialCutsCouple *couple) {
  if (!fIsInitialized) Initialize();                  // builds fCoupleTable/fElementsTable; may run Geant4 EM processes
  if (couple != fLastCouple) {                        // shared "last used" cache, written by all threads
    fLastCouple = couple;
    fLastMuTable = fCoupleTable[fLastCouple];         // another thread can overwrite between these two lines
  }
}

Consequences in MT:

• Racy first-call: simultaneous fInstances.find/operator[]/make_shared on one std::map → undefined behaviour.
• Concurrent lazy Initialize() from several workers: fCoupleTable, fElementsTable and fIsInitialized are written without any lock, and fIsInitialized = true is set at the end of Initialize() — another thread can observe the flag while construction is still in progress.
• Shared fLastCouple / fLastMuTable: CheckLastCall can hand a thread the GateMuTable of a different material (materials differ by orders of magnitude in μ_en/ρ), and the 10.1.1 getters dereference it unconditionally (fLastMuTable->GetMuEnOverRho(energy)), so a null/partial table is a crash and a wrong-material table is a large erroneous coefficient.

Both #1 and #2 feed the same observable: a wrong, occasionally much-too-large μ_en/ρ for a step → upward edep spike → amplified squared-sum spike.

---

4. The fix — and its release status

A single commit fixes the whole family (a92a62bf1, 2026-09-10, "fix thread race issue"; the only commit touching these files between tag 10.1.1 and current master):

• removes the shared fLastEnergyMu/fLastEnergyMuEn/fLastMu/fLastMuEn caches and makes all GetMu*/GetMuEn* methods const with purely local computation (no shared mutable state on the read path);
• makes the per-couple shortcut thread_local instead of a shared fLastCouple/fLastMuTable, and returns 0.0 instead of dereferencing a possibly-null table;
• adds a mutex to GetInstance;
• calls fMaterialMuHandler->Initialize() explicitly once inside GateTLEDoseActor::BeginOfEventAction, under the existing SetEkinMaxMutex, so table construction happens single-threaded before any worker reads the tables.

// master, GateMuTables.cpp
double GateMuTable::GetMuEnOverRho(double energy) const {
  const double log_energy = log(energy);
  int inf = 0, sup = fSize - 1;
  while (sup - inf > 1) { /* binary search */ }
  ...
  return (log_energy > e_inf && log_energy < e_sup)
             ? exp(interpol(...))
             : exp(fMuEn[inf]);
}

Verification of scope: git tag --contains a92a62bf1 is empty, and it is an ancestor of origin/master. So the races above are present in the 10.1.1 wheel the reporter used, and absent on master.

---

5. Residual thread-safety hazards that remain on master

These are not the reported spike mechanism as far as the reported symptoms indicate (see §6), but they are real and in the same area; they deserve separate follow-up.

(a) One G4EmCalculator shared by all worker threads.
GateTLEDoseActor holds G4EmCalculator *fEmCalc (a single instance, created once in BeginOfEventAction) and calls it from every thread in InitializeCSDAForNewGamma() (fEmCalc->GetCSDARange(...)) and FindEkinMaxForTLE(). G4EmCalculator holds G4LossTableManager* manager, whose Instance() is thread-local in MT — so the object created by one thread references that thread's loss-table manager and is then used by all others. Contrast with GateDoseActor::CalculateSPR, which deliberately keeps a thread-local calculator:

auto &local = fThreadLocalDataEdep.Get();
if (!local.emcalc) local.emcalc = std::make_unique<G4EmCalculator>();
auto &emc = *local.emcalc;

The TLE path should do the same. This can corrupt the CSDA-range decision and, for tle_threshold_type = "average range", the sec_ekin that feeds the μ lookup — i.e. it can produce the same class of error if it triggers.

(b) The same voxel images are protected by two different mutexes.
ScoreTLEDepositStep writes cpp_edep_image / cpp_dose_image under SetPixelTLEMutex (GateTLEDoseActor.cpp:283-287), while the inherited GateDoseActor::SteppingAction — which the TLE actor calls for non-TLE gammas and for non-gamma tracks (GateTLEDoseActor.cpp:235,248) — writes the same images under SetPixelMutex (GateDoseActor.cpp:221-225). Different mutexes ⇒ no mutual exclusion between the two paths ⇒ lost/duplicated += updates on the primary edep/dose images in MT. Note this cannot explain the upward spikes (a lost read-modify-write can only under-count), but it is a genuine correctness defect and should be unified to one mutex. The lock order used (SetPixelTLEMutex → SetPixelMutex via ScoreSquaredValue) has no reverse path, so there is no deadlock.

(c) counts is not scored on the TLE path.
ScoreTLEDepositStep never increments cpp_counts_image, whereas GateDoseActor::SteppingAction does. fCountsFlag is therefore silently ignored for TLE-scored steps.

---

6. What is not the cause (ruled out)

• Not a physics/statistics artefact. The co-scored classic DoseActor at the same voxels is clean; a genuine rare energy-deposition event would appear in both.
• Not the edep/dose image mutex split of §5(b). A racy += on a shared image can lose an increment but can never produce a value larger than the sum of the increments, so it cannot generate upward spikes.
• Not the squared-sum bookkeeping. ScoreSquaredValue/FlushSquaredValues operate on thread-local buffers (fThreadLocalDataEdep/fThreadLocalDataDose) and write the global squared image under SetPixelMutex, so the squared accumulation is correctly serialized; it amplifies a spiked edep rather than creating one.

---

7. Verification status and limitations

• The races are established directly from the code: the shared mutable memoization members and the unprotected singleton/lazy-init are present in the 10.1.1 tag and removed on master by a92a62bf1 (which is in no release tag).
• I did not empirically reproduce the spikes in this session: the defect requires the pre-fix binary, and obtaining it would mean reverting/rebuilding an older tree; moreover the race is interleaving-dependent, so a bounded local run is not a reliable trigger. The attribution of the reported signatures to these races is by mechanism, and it accounts for all five reported observations.
• Confirming closure of the issue requires re-running the reporter's high-statistics 8-thread 100–120 keV case on master (i.e. from a build that includes a92a62bf1) and checking that the spikes vanish; that fix must also reach a release, which it currently has not.

---

8. One-line summary

In 10.1.1 the TLE estimator multiplies the photon energy by a μ_en/ρ coefficient fetched through a shared, unsynchronized memoization cache (GateMuTable::fLastEnergyMuEn/fLastMuEn) on a singleton handler whose creation and lazy Initialize() are also unguarded — so concurrent worker threads can return the coefficient of the wrong energy (or wrong material), producing rare upward edep spikes that the squared-sum image amplifies quadratically, while the classic DoseActor stays clean; commit a92a62bf1 removes these shared caches, adds a lookup mutex, uses thread-local caching and initializes the tables once under lock, and is on master but in no released tag, so 10.1.1 users still hit it.

[BishopWolf]

@nkrah This is solved but not released

[GitFuchs]

thank you for taking a look and fixing it! :-)