Geant4 Version: 11.4.1 (MT build)
Operating System: Debian GNU/Linux “forky” (testing)
Compiler/Version: GCC 15.2.0 (Debian 15.2.0-17)
CMake Version: 4.3.3
Summary
G4EmStandardPhysics_option4 (EM Opt4) produces energy deposits that exceed the physically allowed upper bound in a simple G4Box geometry. The anomaly is direction-dependent, it occurs most when the detector is placed on the positive Y axis but not on the negative Y axis (primaries are shoot towards the detector). The issue is reproducible with electrons, positrons, and photons, with positrons having the most anomalous events. Switching to G4EmStandardPhysics_option3 eliminates the anomaly entirely (this is also true for option3 with fUseSafetyPlus step limit).
Geometry
- World: 2×2×2 m³ box of
G4_AIR - Target: 10×20×10 cm³ box of
G4_CESIUM_IODIDE(CsI) - Target position: (0, 50, 0) cm, the 10×10cm end faces the origin
- Particle gun: at the origin (0, 0, 0), firing toward the detector in direction (0, 1, 0)
- Scoring: total energy deposit per event in the CsI volume
Physics List
G4EmStandardPhysics_option4(EM Opt4)- No other physics processes registered
Primary Particles
Each particle type was simulated with 1,000,000 events at 50 MeV kinetic energy:
| Particle | Kinetic Energy | Direction |
|---|---|---|
| e⁻ | 50 MeV | (0, 1, 0) |
| e⁺ | 50 MeV | (0, 1, 0) |
| γ | 50 MeV | (0, 1, 0) |
Results
Anomalous Event Fraction
I set a energy deposition threshold:
$$E_{\text{threshold}} = (1 + \varepsilon) \times (E_{\text{kin}} + 2 m_e c^2)$$
where $\varepsilon$ is the machine epsilon for float. Any single event with Edep exceeding this threshold is flagged as anomalous.
| Particle | Events Simulated | Anomalous Events | Fraction |
|---|---|---|---|
| e⁻ | 1,000,000 | 15 | 0.0015% |
| e⁺ | 1,000,000 | 380 | 0.0380% |
| γ | 1,000,000 | 153 | 0.0153% |
Energy Deposition Spectra
In the plots below, the red line is EM Opt4 and the blue line is EM Opt3 + fUseSafetyPlus.
Electrons:
Positrons:
Photons:
Observations
-
Directional dependence: When the CsI box is placed at (0, −50, 0) cm (i.e., on the negative Y axis, with the particle still fired toward it from the origin), no anomalous events occur.
-
EM Opt3 is unaffected: Switching the physics list to
G4EmStandardPhysics_option3(optionally withfUseSafetyPlus) eliminates all anomalous events. -
Positron low-energy discrepancy: For positrons, EM Opt3 +
fUseSafetyPlusshows more events in the sub-20 MeV region compared to EM Opt4 (visible as the blue line being above the red line below ~20 MeV). Given that the total event count is the same, this suggests that at least some of the anomalous high-energy events seen with Opt4 may be events whose true energy deposit belongs in the lower-energy region.
Reproduction Code
A reproduction program (main.cxx) is provided below.
Build:
mkdir build && cd build && cmake .. && make
Run:
./opt4bug 1000000 0 50 0 e+ 50
The program writes opt4bug_edep.root with per-event Edep values and prints anomalous events to stdout.
CMakeLists.txt
cmake_minimum_required(VERSION 3.16)
project(Opt4BugRepro)
find_package(Geant4 REQUIRED)
add_executable(opt4bug main.cxx)
target_link_libraries(opt4bug ${Geant4_LIBRARIES})
target_compile_options(opt4bug PRIVATE -Wall -Wextra -pedantic)
main.cxx
// Bug reproduction: G4EmStandardPhysics_option4 + G4Box
// Output: per-event Edep written to ROOT ntuple (opt4bug_edep.root)
//
// Build: mkdir build && cd build && cmake .. && make
// Run: ./opt4bug nEvents x y z particleName energy
// nEvents Number of events
// x y z Box position in cm
// particleName Incident particle (e.g. gamma, e-, proton)
// energy Incident kinetic energy in MeV
//
// Long side of the box always faces the origin.
// Photon always incident on the near small face.
#include "G4Accumulable.hh"
#include "G4AccumulableManager.hh"
#include "G4AnalysisManager.hh"
#include "G4Box.hh"
#include "G4EmStandardPhysics_option3.hh"
#include "G4EmStandardPhysics_option4.hh"
#include "G4Event.hh"
#include "G4LogicalVolume.hh"
#include "G4NistManager.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4PhysicalConstants.hh"
#include "G4PVPlacement.hh"
#include "G4Run.hh"
#include "G4RunManager.hh"
#include "G4RunManagerFactory.hh"
#include "G4Step.hh"
#include "G4SystemOfUnits.hh"
#include "G4Threading.hh"
#include "G4UserEventAction.hh"
#include "G4UserRunAction.hh"
#include "G4UserSteppingAction.hh"
#include "G4VModularPhysicsList.hh"
#include "G4VUserActionInitialization.hh"
#include "G4EmParameters.hh"
#include "G4VUserDetectorConstruction.hh"
#include "G4VUserPrimaryGeneratorAction.hh"
#include <atomic>
#include <cmath>
#include <iomanip>
#include <memory>
#include <string>
namespace {
std::atomic<int> gProgressCounter{0};
constexpr int kProgressInterval{1000};
} // namespace
// ============================================================================
// DetectorConstruction
// ============================================================================
class DetectorConstruction : public G4VUserDetectorConstruction {
public:
explicit DetectorConstruction(const G4ThreeVector& pos) :
fPosition{pos} {}
auto Construct() -> G4VPhysicalVolume* override {
auto nist{G4NistManager::Instance()};
auto checkOverlaps{true};
// --- World: 2x2x2 m^3 box of G4_AIR ---
auto solidWorld(new G4Box{"World", 1.0 * m, 1.0 * m, 1.0 * m});
auto logicWorld(new G4LogicalVolume{solidWorld, nist->FindOrBuildMaterial("G4_AIR"), "World"});
auto physWorld(new G4PVPlacement{nullptr, G4ThreeVector{}, logicWorld, "World", nullptr, false, 0, checkOverlaps});
// --- Target: 10×20×10 cm³ CsI box, long side (Y) faces origin ---
auto solidECAL(new G4Box{"ECAL", 5.0 * cm, 10.0 * cm, 5.0 * cm});
auto logicECAL(new G4LogicalVolume{solidECAL, nist->FindOrBuildMaterial("G4_CESIUM_IODIDE"), "ECAL"});
fTargetVolume = logicECAL;
// Build rotation: local Y axis → dirToOrigin (= -fPosition.unit())
auto colY{-fPosition.unit()};
auto ref((std::abs(colY.x()) < 0.9) ? G4ThreeVector{1, 0, 0} : G4ThreeVector{0, 0, 1});
auto colZ{colY.cross(ref).unit()};
auto colX{colY.cross(colZ).unit()};
auto rot(new G4RotationMatrix{colX, colY, colZ});
new G4PVPlacement{rot, fPosition, logicECAL,
"ECAL", logicWorld, false, 0, checkOverlaps};
return physWorld;
}
auto GetECALVolume() const -> G4LogicalVolume* { return fTargetVolume; }
private:
G4LogicalVolume* fTargetVolume{nullptr};
G4ThreeVector fPosition;
};
// ============================================================================
// PrimaryGeneratorAction
// ============================================================================
class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
public:
PrimaryGeneratorAction(const G4ThreeVector& dir,
const std::string& particleName,
double energy) :
fGun{1} {
fGun.SetParticleDefinition(G4ParticleTable::GetParticleTable()->FindParticle(particleName));
fGun.SetParticleEnergy(energy);
fGun.SetParticlePosition(G4ThreeVector{0.0, 0.0, 0.0});
fGun.SetParticleMomentumDirection(dir.unit());
}
auto GeneratePrimaries(G4Event* event) -> void override {
fGun.GeneratePrimaryVertex(event);
}
auto GetGun() const -> const G4ParticleGun* { return &fGun; }
private:
G4ParticleGun fGun;
};
// ============================================================================
// RunAction
// ============================================================================
class RunAction : public G4UserRunAction {
public:
RunAction() {
G4AccumulableManager::Instance()->Register(fTotalEvents);
G4AccumulableManager::Instance()->Register(fAnomalousEvts);
// --- ROOT analysis setup ---
auto am{G4AnalysisManager::Instance()};
am->SetDefaultFileType("root");
am->SetNtupleMerging(true); // required for MT merging
am->SetVerboseLevel(0); // suppress per-event analysis chatter
am->SetFileName("opt4bug");
// Book ntuple with a single double column: Edep (MeV)
am->CreateNtuple("Edep", "Energy deposit per event");
am->CreateNtupleFColumn("Edep");
am->FinishNtuple();
// Use a separate file for the ntuple (opt4bug_edep.root)
am->SetNtupleFileName(0, "opt4bug_edep");
}
auto BeginOfRunAction(const G4Run*) -> void override {
G4AccumulableManager::Instance()->Reset();
auto am{G4AnalysisManager::Instance()};
am->Reset();
am->OpenFile();
}
auto EndOfRunAction(const G4Run*) -> void override {
G4AccumulableManager::Instance()->Merge();
auto nTotal{fTotalEvents.GetValue()};
auto nAnomal{fAnomalousEvts.GetValue()};
if (IsMaster()) {
auto fraction{nTotal > 0 ? 100.0 * nAnomal / nTotal : 0.0};
G4cout << G4endl
<< "============================================================" << G4endl
<< " Run Summary (Master) — " << nTotal << " events" << G4endl
<< "============================================================" << G4endl
<< " Anomalous events: " << nAnomal
<< " (" << std::fixed << std::setprecision(3) << fraction << "%)" << G4endl
<< " ROOT output: opt4bug_edep.root" << G4endl
<< "============================================================" << G4endl;
}
// Write and close analysis files
auto am{G4AnalysisManager::Instance()};
am->Write();
am->CloseFile(false);
}
auto CountEvent() -> void {
fTotalEvents += 1;
}
auto CountAnomalous() -> void {
fAnomalousEvts += 1;
}
private:
G4Accumulable<int> fTotalEvents{0};
G4Accumulable<int> fAnomalousEvts{0};
};
// ============================================================================
// EventAction
// ============================================================================
class EventAction : public G4UserEventAction {
public:
EventAction(RunAction* runAction, double threshold) :
fRunAction{runAction},
fThreshold{threshold} {}
auto BeginOfEventAction(const G4Event*) -> void override {
fEdep = 0.0;
}
auto EndOfEventAction(const G4Event*) -> void override {
fRunAction->CountEvent();
if (fEdep > fThreshold) {
fRunAction->CountAnomalous();
G4cout << " *** ANOMALOUS! Edep = " << fEdep / MeV
<< " MeV > " << fThreshold / MeV << " MeV ***" << G4endl;
}
// Fill ROOT ntuple
auto am{G4AnalysisManager::Instance()};
am->FillNtupleFColumn(0, fEdep / MeV);
am->AddNtupleRow();
// Print progress every kProgressInterval completed events (across all threads)
auto done{++gProgressCounter};
if (done % kProgressInterval == 0) {
G4cout << ">>> Progress: " << done << " events completed" << G4endl;
}
}
auto AddEdep(double e) -> void { fEdep += e; }
private:
RunAction* fRunAction{nullptr};
double fThreshold;
double fEdep{0.0};
};
// ============================================================================
// SteppingAction
// ============================================================================
class SteppingAction : public G4UserSteppingAction {
public:
SteppingAction(EventAction* evtAction) :
fEventAction{evtAction} {}
auto UserSteppingAction(const G4Step* step) -> void override {
if (!fScoringVolume) {
const auto det{
static_cast<const DetectorConstruction*>(
G4RunManager::GetRunManager()->GetUserDetectorConstruction())};
fScoringVolume = det->GetECALVolume();
}
auto vol{step->GetPreStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume()};
if (vol == fScoringVolume) {
fEventAction->AddEdep(step->GetTotalEnergyDeposit());
}
}
private:
EventAction* fEventAction{nullptr};
G4LogicalVolume* fScoringVolume{nullptr};
};
// ============================================================================
// ActionInitialization
// ============================================================================
class ActionInitialization : public G4VUserActionInitialization {
public:
ActionInitialization(const G4ThreeVector& direction,
const std::string& particleName,
double energy,
double threshold) :
fDirection{direction},
fParticleName{particleName},
fEnergy{energy},
fThreshold{threshold} {}
auto BuildForMaster() const -> void override {
SetUserAction(new RunAction);
}
auto Build() const -> void override {
auto runAction{new RunAction};
auto evtAction(new EventAction{runAction, fThreshold});
SetUserAction(new PrimaryGeneratorAction{fDirection, fParticleName, fEnergy});
SetUserAction(runAction);
SetUserAction(evtAction);
SetUserAction(new SteppingAction{evtAction});
}
private:
G4ThreeVector fDirection;
std::string fParticleName;
double fEnergy;
double fThreshold;
};
// ============================================================================
// main
// ============================================================================
auto main(int argc, char** argv) -> int {
// --- Require exactly 6 positional arguments ---
if (argc != 7) {
G4cout << "Usage: " << argv[0]
<< " nEvents x y z particleName energy(MeV)" << G4endl
<< G4endl
<< " nEvents Number of events" << G4endl
<< " x y z Box position in cm" << G4endl
<< " particleName Incident particle (e.g. gamma, e-, proton)" << G4endl
<< " energy Incident kinetic energy in MeV" << G4endl
<< G4endl
<< "Example: " << argv[0] << " 100000 0 50 0 gamma 60" << G4endl;
return 1;
}
// --- Parse arguments in strict order ---
int nEvents{std::atoi(argv[1])};
G4ThreeVector boxPos{std::atof(argv[2]) * cm,
std::atof(argv[3]) * cm,
std::atof(argv[4]) * cm};
std::string particleName{argv[5]};
double energy{std::atof(argv[6]) * MeV};
if (nEvents <= 0) {
G4cout << "Error: nEvents must be > 0" << G4endl;
return 1;
}
auto direction{boxPos.unit()}; // from origin to box
// --- Run manager (MT) ---
G4MTRunManager runManager;
// --- Detector ---
runManager.SetUserInitialization(new DetectorConstruction{boxPos});
// --- Physics: EM only ---
auto physList{new G4VModularPhysicsList{}};
physList->RegisterPhysics(new G4EmStandardPhysics_option4{});
G4EmParameters::Instance()->SetMscStepLimitType(fUseSafetyPlus);
G4EmParameters::Instance()->SetMscMuHadStepLimitType(fUseSafetyPlus);
runManager.SetUserInitialization(physList);
// --- User actions ---
auto threshold{(1 + FLT_EPSILON) * (energy + 2 * electron_mass_c2)};
runManager.SetUserInitialization(new ActionInitialization{direction, particleName, energy, threshold});
// --- Initialize ---
runManager.SetNumberOfThreads(G4Threading::G4GetNumberOfCores());
runManager.Initialize();
// --- Batch run ---
G4cout << G4endl
<< "============================================================" << G4endl
<< " Starting " << nEvents << " events with G4EmStandardPhysics_option4" << G4endl
<< " Primary: " << particleName
<< ", E_kin = " << energy / MeV << " MeV"
<< ", E_total = " << threshold / MeV << " MeV" << G4endl
<< " Anomaly threshold: Edep > " << threshold / MeV << " MeV" << G4endl
<< " Target: G4_CESIUM_IODIDE box 10*20*10 cm3 (G4Box)" << G4endl
<< " Box position: (" << boxPos.x() / cm << ", " << boxPos.y() / cm << ", " << boxPos.z() / cm << ") cm" << G4endl
<< " Photon direction: (" << direction.x() << ", " << direction.y() << ", " << direction.z() << ")" << G4endl
<< " ROOT file: opt4bug_edep.root" << G4endl
<< "============================================================" << G4endl
<< G4endl;
runManager.BeamOn(nEvents);
return 0;
}