Multi threading changes simulation result

Geometry, Fields and Transportation
1

Dear Geant4 forum,
I am experiencing issues with multi threading. It seems that running my simulation with more than 1 thread, like using:

/run/numberOfThreads 4

changes the output.

Let me take one step back. I am validating a “phase space” (PS) that makes use of the GPS. Please note that I am using a series of macro files to define specific distribution for each particle and quantity, it is not like a “traditional” phase space.
While validating it I have noticed an artifacts on the $\theta_P$ angle (momentum along the Z axis). An anomalous peak is present around $\theta_P > 2.55$ rad.
It seems that single threaded runs are not affected and can reproduce the original distribution.
To clarify the issue I have made the following plot comparing all:

  • the original phase space, directly from the macro (green).
  • The $\theta_P$ angle scored single threaded (blue).
  • The $\theta_P$ angle scored multi threaded (orange)
    Comparison
    Comparison1920×1440 92.1 KB

The same seed and the same number of primary events were used.

I have made a git repo capable of recreating the issue here:

It is based on the B1 example with few modifications:

  • GPS as primary generator
  • Dedicated PrimaryHit class to store events at the beginning and end of EventAction:
...
void EventAction::BeginOfEventAction(const G4Event *event) {
  auto etheta =
      event->GetPrimaryVertex()->GetPrimary()->GetMomentum().getTheta() * rad;
  fPrimaryHit = new G4TPrimaryHit(etheta);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void EventAction::EndOfEventAction(const G4Event *) {
  auto analysisManager = G4AnalysisManager::Instance();
  analysisManager->FillH1(1, fPrimaryHit->GetTheta());
}
...
  • A dummy physics list using only the particles transport and definition

Further notes:

  • this seems to affect only the $\theta_P$ angle. Other quantities like kinetic energy, $\phi$ angle, and position are correctly reproduced.
  • The origin of the phase space is on the X-Y plate at Z = -10 mm.
  • Using QBBC_EMY do not change the outcome.
  • I suspect an issue with seeding?

I have tested and experienced the same effect on both on debian and arch (see detail below).

Geant4 Version: 11.4.3
Operating System: Linux SADPC198 7.2.14-arch1-1 and Linux lxir117 6.12.94+deb13-amd64
Compiler/Version: g++ (GCC) 16.2.1 20260810 and g++ (Ubuntu 15.2.0-4ubuntu4) 15.2.0
CMake Version: 4.4.5 and 3.31.6

Is the issue related to the modification I made to the example or is this a bug of geant4?

Let me know if additional details are needed and thank you for your help.
Enrico