Energy losses in material depends on when I replace the dE/dx table

Alexis_Francheteau · February 27, 2026, 4:01pm

Hello,

For the sake of my simulation I wish to use a custom dE/dx table, based on the DPASS model, however, depending on when I add the table (using AddDEDXTable()) I get different energy losses.

I have this configuration:

Geant4 Version: v11.4.0
Operating System: Ubuntu 24.04
Compiler/Version:gcc 13.3.0
CMake Version: 3.28.3

I placed my custom table in :

/path/to/geant/install/share/Geant4/data/G4EMLOW8.8/ion_stopping_data/MyTable

Then I create a public function AddDPASSTable() in my PhysicsList, defined as followed:

NFSPhysicsList::NFSPhysicsList() 
    : QGSP_INCLXX_HP(),
      fDpassAdded(false) // flag initialsation
{
    G4cout << "Entered NFSPhysicsList" << G4endl;
    // Replace EM physics with option4
    ReplacePhysics(new G4EmStandardPhysics_option4(1,"useExternalDEDX"));
    
    // Set DPASS directory
    fDpassDir = "ion_stopping_data/dpass21.06";
    fDpassPath = "/path/to/geant/install/share/Geant4/data/G4EMLOW8.8";
}

NFSPhysicsList::~NFSPhysicsList()
{
}

void NFSPhysicsList::ConstructProcess()
{
    // First construct all standard processes
    QGSP_INCLXX_HP::ConstructProcess();
    
}

void PhysicsList::AddDPASSTable()
{
    G4cout << "\n==========================" << G4endl;
    G4cout << "=== ADDING DPASS TABLE ===" << G4endl;
    G4cout << "==========================\n" << G4endl;
    
    // Verify DPASS data exists
    std::ifstream testFile(fDpassPath+ '/' + fDpassDir + "73/z2_14.dat");
    if(!testFile.good())
    {
        G4Exception("PhysicsList::AddDPASSTable()", 
                    "DPASS001", FatalException,
                    ("Cannot find DPASS data at: " fDpassPath+ '/' + fDpassDir).c_str());
        return;
    }
    
    testFile.close();
    G4cout << "✓ DPASS data directory verified: " << fDpassDir << G4endl;
    
    // Get GenericIon
    G4ParticleDefinition* genericIon = G4GenericIon::GenericIon();
    if(!genericIon)
    {
        G4cout << "⚠ GenericIon not found, skipping DPASS" << G4endl;
        return;
    }
    else
        G4cout << "✓ GenericIon Found! " << G4endl;
    
    // Get process manager
    G4ProcessManager* processManager = genericIon->GetProcessManager();
    if(!processManager)
    {
        G4cout << "⚠ No process manager for GenericIon" << G4endl;
        return;
    }
    else
        G4cout << "✓ Process Manager Found! " << G4endl;
    
    // Find ionIoni process
    G4VProcess* ionIoniProcess = processManager->GetProcess("ionIoni");
    if(!ionIoniProcess)
    {
        G4cout << "⚠ ionIoni process not found" << G4endl;

        // DEBUG: List all processes
        G4ProcessVector* pList = processManager->GetProcessList();
        G4cout << "Available processes for GenericIon:" << G4endl;

        for(size_t i = 0; i < pList->size(); i++)
            G4cout << "  - " << (*pList)[i]->GetProcessName() << G4endl;
        
        return;
    }
    else
        G4cout << "✓ ionIoni process Found! " << G4endl;
    
    // Cast to energy loss process
    G4VEnergyLossProcess* energyLossProcess = 
        dynamic_cast<G4VEnergyLossProcess*>(ionIoniProcess);
    if(!energyLossProcess)
    {
        G4cout << "⚠ Failed to cast to G4VEnergyLossProcess" << G4endl;
        return;
    }
    else
        G4cout << "✓ Cast to G4VEnergyLossProcess Successful! " << G4endl;
    
    // Get the model
    G4VEmModel* emModel = energyLossProcess->EmModel(0);
    if(!emModel) {
        G4cout << "⚠ No EM model found" << G4endl;
        return;
    }
    else
        G4cout << "✓ EM model Found!" << G4endl;
    
    G4cout << "Model name: " << emModel->GetName() << G4endl;
    
    // Check for ParamICRU73
    if(emModel->GetName() == "ParamICRU73") {
        
        G4IonParametrisedLossModel* icruModel = 
            dynamic_cast<G4IonParametrisedLossModel*>(emModel);
        
        if(!icruModel) {
            G4cout << "⚠ Failed to cast to G4IonParametrisedLossModel" << G4endl;
            return;
        }
        else
            G4cout << "✓ Cast to G4IonParametrisedLossModel Successful!" << G4endl;

        
        
        // Remove ICRU73 and add DPASS
        bool rem_icru = icruModel->RemoveDEDXTable("ICRU73");
        if (rem_icru)
            G4cout << "✓ ICRU73 table removed" << G4endl;

        else
            G4cout << "ICRU73 table removal returned False !" << G4endl;

        G4cout << "Adding the DPASS table ..." << G4endl;
        bool success = icruModel->AddDEDXTable("DPASS",
                                               new G4IonStoppingData(fDpassDir, false),
                                               new G4IonDEDXScalingICRU73(93, 102));

        if (success)
            G4cout << "✓ AddDEDX returned True value!" << G4endl;
        else
            G4cout << "⚠ WARNING ⚠ : AddDEDX returned False value!" << G4endl;

        
    }
    else
    {
        G4cout << "⚠ ParamICRU73 model not found!" << G4endl;
        G4cout << "  Current model: " << emModel->GetName() << G4endl;
        G4cout << "  This means DPASS cannot be installed." << G4endl;
    }
    
    G4cout << "=== DPASS TABLE SETUP COMPLETE ===\n" << G4endl;
}

which main role is to remove the standard icru table and replace it with my DPASS table.

In order to have this AddDEDXTable() call run properly, I had to patch the G4EmStandardPhysics_option4.* like this: (patch given by V. Ivantchenko)

////////////////////////////////////////
// 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()
...
...

Then I ran the following tests, for an incoming ion (Z=61, A=155, E=69.08 MeV) passing through a Mylar foil (thickness: 0.53 um):

case 1: Using icru73 table

// main.cc

{...}

    // Construct the run manager
    G4RunManager * runManager = new G4RunManager() ;

    // Initialize the geometry

    // Initialize the physics 
    NFSPhysicsList* physicsList = new NFSPhysicsList();
    physicsList->SetVerboseLevel(-1);
    runManager->SetUserInitialization(physicsList);

{...}
    
    // Initialize the run manager
    G4cout << "\nInitialization...\n" << G4endl;
    runManager->Initialize();


    // Call BeamOn(0) to fully initialize tables
    G4cout << "Calling BeamOn(0) to fully initialize the tables..." << G4endl;
    runManager->BeamOn(0);
    

    // Explicitly call for DPASS add
    // physicsList->AddDPASSTable();

DeltaE = 8.22 MeV

case 2: Adding DPASS table after BeamOn(0)

// main.cc

{...}

    // Construct the run manager
    G4RunManager * runManager = new G4RunManager() ;

    // Initialize the geometry

    // Initialize the physics 
    NFSPhysicsList* physicsList = new NFSPhysicsList();
    physicsList->SetVerboseLevel(-1);
    runManager->SetUserInitialization(physicsList);

{...}
    
    // Initialize the run manager
    G4cout << "\nInitialization...\n" << G4endl;
    runManager->Initialize();


    // Call BeamOn(0) to fully initialize tables
    G4cout << "Calling BeamOn(0) to fully initialize the tables..." << G4endl;
    runManager->BeamOn(0);
    

    // Explicitly call for DPASS add
    physicsList->AddDPASSTable();

DeltaE = 6.72 MeV

case 3: Adding DPASS table before BeamOn(0)

// main.cc

{...}

    // Construct the run manager
    G4RunManager * runManager = new G4RunManager() ;

    // Initialize the geometry

    // Initialize the physics 
    NFSPhysicsList* physicsList = new NFSPhysicsList();
    physicsList->SetVerboseLevel(-1);
    runManager->SetUserInitialization(physicsList);

{...}
    
    // Initialize the run manager
    G4cout << "\nInitialization...\n" << G4endl;
    runManager->Initialize();


    // Explicitly call for DPASS add
    physicsList->AddDPASSTable();

    // Call BeamOn(0) to fully initialize tables
    G4cout << "Calling BeamOn(0) to fully initialize the tables..." << G4endl;
    runManager->BeamOn(0);

DeltaE = 8.22 MeV (same as not adding the table)

case 4: Adding DPASS table before BeamOn(0) + explicit G4IonParametrisedLossModel::Initialise() call

// NFSPhysicsList.cc

{...}
        // Remove ICRU73 and add DPASS
        G4DataVector cuts;
        icruModel->Initialise(genericIon,cuts);
        bool rem_icru = icruModel->RemoveDEDXTable("ICRU73");
        if (rem_icru)
            G4cout << "✓ ICRU73 table removed" << G4endl;

        else
            G4cout << "ICRU73 table removal returned False !" << G4endl;

        G4cout << "Adding the DPASS table ..." << G4endl;
        bool success = icruModel->AddDEDXTable("DPASS",
                                               new G4IonStoppingData(fDpassDir, false),
                                               new G4IonDEDXScalingICRU73(93, 102));

{...}


// main.cc

{...}

    // Construct the run manager
    G4RunManager * runManager = new G4RunManager() ;

    // Initialize the geometry

    // Initialize the physics 
    NFSPhysicsList* physicsList = new NFSPhysicsList();
    physicsList->SetVerboseLevel(-1);
    runManager->SetUserInitialization(physicsList);

{...}
    
    // Initialize the run manager
    G4cout << "\nInitialization...\n" << G4endl;
    runManager->Initialize();


    // Explicitly call for DPASS add
    physicsList->AddDPASSTable();

    // Call BeamOn(0) to fully initialize tables
    G4cout << "Calling BeamOn(0) to fully initialize the tables..." << G4endl;
    runManager->BeamOn(0);

DeltaE = 7.45 MeV

case 5: Adding DPASS table after BeamOn(0) + explicit G4IonParametrisedLossModel::Initialise() call

// NFSPhysicsList.cc

{...}
        // Remove ICRU73 and add DPASS
        G4DataVector cuts;
        icruModel->Initialise(genericIon,cuts);
        bool rem_icru = icruModel->RemoveDEDXTable("ICRU73");
        if (rem_icru)
            G4cout << "✓ ICRU73 table removed" << G4endl;

        else
            G4cout << "ICRU73 table removal returned False !" << G4endl;

        G4cout << "Adding the DPASS table ..." << G4endl;
        bool success = icruModel->AddDEDXTable("DPASS",
                                               new G4IonStoppingData(fDpassDir, false),
                                               new G4IonDEDXScalingICRU73(93, 102));

{...}


// main.cc

{...}

    // Construct the run manager
    G4RunManager * runManager = new G4RunManager() ;

    // Initialize the geometry

    // Initialize the physics 
    NFSPhysicsList* physicsList = new NFSPhysicsList();
    physicsList->SetVerboseLevel(-1);
    runManager->SetUserInitialization(physicsList);

{...}
    
    // Initialize the run manager
    G4cout << "\nInitialization...\n" << G4endl;
    runManager->Initialize();

    // Call BeamOn(0) to fully initialize tables
    G4cout << "Calling BeamOn(0) to fully initialize the tables..." << G4endl;
    runManager->BeamOn(0);


    // Explicitly call for DPASS add
    physicsList->AddDPASSTable();

DeltaE = 6.72 MeV (explicit initialise() call doesn’t affect)

I was quite puzzled with these results, and I couldn’t find any similar issue in a previous forum post, so I first tried asking a LLM (Anthropic Claude), which gave me a explanation that I reformulate here:

It seems that Geant4 works with 2 levels:
• Lvl 1: registered tables in memory
• Lvl 2: some mesh of eloss per ion/per material are written in the cache
If you fall on a value in between two cache values, only use them, otherwise, you compute the missing points from the lvl 1 table. This means that if the cache points and the memory table are not from the same table, you have a mixing of eloss tables. According to Anthropic Claude, the cache is supposed to be initialized when the G4IonParametrisedLossModel::Initialise() is called (by default in the first call of G4RunManager::BeamOn()).

Therefore I have three questions:

Is the LLM answer correct? and to what degree?
If it is, how can I be sure that the cache is emptied when I call the AddDEDXTable()?
If not, can someone explain me the difference between the cases?

Alex2671 · March 4, 2026, 11:21am

Could you first of all check if the dedx table has been uploaded correctly? I usually do it by writing something like this in stepping action:


    G4EmCalculator calc;
    G4NistManager* man = G4NistManager::Instance();
    G4Material* material = man->FindOrBuildMaterial("G4_Au");
    G4ParticleDefinition* recoildefinition;
    G4IonTable* ito = G4ParticleTable::GetParticleTable()->GetIonTable();
    int Z = 22;
    int A = 48;
    recoildefinition = ito->GetIon(Z, A);

    size_t N = 20;
    double lowE = 1./100*48;
    double highE = 1.*48;
    std::ofstream file;
    file.open("48TiAu.csv", std::ios_base::app);
    for(unsigned i = 0; i <= N; ++i){
        double KinE = log(lowE) + i*(log(48) - log(48./100))/N;
        auto dEdX = calc.GetDEDX(exp(KinE), recoildefinition, material);
        file << exp(KinE)/A << " " << dEdX << std::endl;
    }
    file.close();
    exit(0);