Geant4 Version: 10.7.3 & 11.4.0
Operating System: Ubuntu 24.04
Compiler/Version: C++17
CMake Version: 3.28.03
In the case of my simulation, I want to use an external DEDX table to compute energy losses of heavy ions of ~1MeV/A in Mylar foils. When using this external table in v11.4.0, the computed energy losses were not matching the one found in the table. Therefore I investigated this mismatch for a simple case of a 12C passing through a mylar foil to compute back the DEDX table. The code ion_stopping_repro.cc joined to this post is a minimal portable version of my code.
ion_stopping_repro.cc (14.9 KB)
I designed a simple simulation where a 12C ion is shot for a grid in log scale between 0.01 MeV/A and 100 MeV/A. To be sure that, for each energy point, the mylar thickness is small enough to let all ions pass through, it is changed at each iteration, ranging from 10 nm to 1mm. I benchmarked it in v10.7.3 where I can properly reproduce both icru73 table (look up by G4IonParametrisedLossModel) and my external table (DPASS, by changing the path in the geant4 source code).
However, it shows really surprising and puzzling results when running in v11.4.0. When running using Lindhard-Sorensen parametrization
RegisterPhysics(new G4EmStandardPhysics_option4(0)); // Lindhard-Sorensen model
I get a discontinuity, which sharpness depends on the production cut and insensitive to step limitations.
The most puzzling part is when using G4IonParametrisedLossModel. This issue requires a small patch to G4EmStandardPhysics_option4 (shown below) to expose the useExternalDEDX switch, since stock 11.4.0 hardcodes G4LindhardSorensenIonModel for ions in option4. (this has been implemented in G4 v11.4.1)
////////////////////////////////////////
// PhysicsList.hh
////////////////////////////////////////
class PhysicsList : public G4VModularPhysicsList
{
public:
PhysicsList() : G4VModularPhysicsList()
{
SetVerboseLevel(0);
// Best EM physics for ions (includes nuclear stopping)
// RegisterPhysics(new G4EmStandardPhysics_option4(0));
/RegisterPhysics(new G4EmStandardPhysics_option4(1,"useExternalDEDX")); // allow the use of G4IonParametrisedLossModel through a patch described below
// Allows external step-limit commands
RegisterPhysics(new G4StepLimiterPhysics());
}
~PhysicsList() override = default;
void AddDPASSTable();
};
////////////////////////////////////////
// G4EmStandardPhysics_option4.hh
////////////////////////////////////////
#ifndef G4EmStandardPhysics_option4_h
#define G4EmStandardPhysics_option4_h 1
#include "G4VPhysicsConstructor.hh"
#include "globals.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
class G4EmStandardPhysics_option4 : public G4VPhysicsConstructor
{
public:
explicit G4EmStandardPhysics_option4(G4int ver=1, const G4String& name="");
~G4EmStandardPhysics_option4() override;
void ConstructParticle() override;
void ConstructProcess() override;
private:
G4bool fUseExternalDEDX{false}; // PATCHED PART
};
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#endif
////////////////////////////////////////
// G4EmStandardPhysics_option4.cc
////////////////////////////////////////
// factory
#include "G4PhysicsConstructorFactory.hh"
//
G4_DECLARE_PHYSCONSTR_FACTORY(G4EmStandardPhysics_option4);
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4EmStandardPhysics_option4::G4EmStandardPhysics_option4(G4int ver,
const G4String& name)
: G4VPhysicsConstructor("G4EmStandard_opt4")
{
SetVerboseLevel(ver);
G4EmParameters* param = G4EmParameters::Instance();
param->SetDefaults();
param->SetVerbose(ver);
param->SetGeneralProcessActive(true);
param->SetMinEnergy(100*CLHEP::eV);
param->SetLowestElectronEnergy(100*CLHEP::eV);
param->SetNumberOfBinsPerDecade(20);
param->ActivateAngularGeneratorForIonisation(true);
param->SetStepFunction(0.2, 10*CLHEP::um);
param->SetStepFunctionMuHad(0.1, 50*CLHEP::um);
param->SetStepFunctionLightIons(0.1, 20*CLHEP::um);
param->SetStepFunctionIons(0.1, 1*CLHEP::um);
param->SetUseMottCorrection(true); // use Mott-correction for e-/e+ msc gs
param->SetMscStepLimitType(fUseSafetyPlus); // for e-/e+ msc gs
param->SetMscSkin(3); // error-free stepping for e-/e+ msc gs
param->SetMscRangeFactor(0.08); // error-free stepping for e-/e+ msc gs
param->SetMuHadLateralDisplacement(true);
param->SetFluo(true);
param->SetUseICRU90Data(true);
param->Set3GammaAnnihilationOnFly(true);
param->SetFluctuationType(fUrbanFluctuation);
param->SetMaxNIELEnergy(1*CLHEP::MeV);
param->SetPositronAtRestModelType(fAllisonPositronium);
SetPhysicsType(bElectromagnetic);
if (!name.empty()) { fUseExternalDEDX = true; } // PATCHED PART
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4EmStandardPhysics_option4::~G4EmStandardPhysics_option4() = default;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void G4EmStandardPhysics_option4::ConstructParticle()
{
// minimal set of particles for EM physics
G4EmBuilder::ConstructMinimalEmSet();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void G4EmStandardPhysics_option4::ConstructProcess()
{
...
...
if (fUseExternalDEDX)
ionIoni->SetEmModel(new G4IonParametrisedLossModel());
else
ionIoni->SetEmModel(new G4LindhardSorensenIonModel());
}
After Initialize() and one event, dumping the active EM models on the 12C ion shows a single model: [0] ParamICRU73 0.0001 – 1e+08 MeV covering the full range. So the ICRU 73 model is selected.
In principle I should reproduce exactly the icru table, but here, I can reproduce it but only above a threshold, threshold that depends on the production cut. With production cut = 1 µm, agreement starts above E/A ≈ 5 MeV/u. With cut = 0.1 mm, agreement starts above E/A ≈ 0.45 MeV/u. The threshold scales roughly with the cut.
dedx_icru.pdf (13.9 KB)
Here are some checks that I already discarded:
- Step size:
SetStepFunctionIons(0.1, 1 nm)and(0.01, 1 µm)give identical results. - Energy-loss fluctuations: disabling via
SetLossFluctuations(false)does not change the curve. - Multiple scattering: disabling
mscdoes not change the curve. - Secondary tracking: integrating the kinetic energy of escaping secondaries does not change the curve below E/A ≈ 23 MeV/u (where T_max < cut, no δ-rays produced).
An incidental observation: calling RemoveDEDXTable("ICRU73") followed by AddDEDXTable("DPASS", …) after Initialize() produces a scan that agrees with ICRU 73 across the full energy range — even though the resulting peak height matches ICRU 73’s 8.2 MeV·cm²/mg rather than DPASS’s 9.9 MeV·cm²/mg, suggesting DPASS itself is not actually installed. It appears the table-replacement call triggers a re-initialisation that fixes the underlying issue. It is independent of the path given: passing the icru path “ion_stopping_data/icru“ gives the same result.
I can provide dedx outputs from my simulation if needed.
To me it looks like something in G4 v11.4.0 is not working properly when using G4IonParametrizedLossModel, unless there are some options/parameters that I have to specify in my main?
Any help is welcome!