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);

Alexis_Francheteau · March 5, 2026, 9:56am

Thank you for your answer, I implemented what you asked (for Al instead of Au, and putting the output in Mev/(mg/cm2)). Here are the outputs:

// Results without the DPASS table
0.01 0.0825619
0.0125893 2.73545
0.0158489 3.06575
0.0199526 3.43548
0.0251189 3.85079
0.0316228 4.41661
0.0398107 5.08037
0.0501187 5.85899
0.0630957 6.77384
0.0794328 7.84267
0.1 9.0808
0.125893 10.4954
0.158489 12.0792
0.199526 13.7932
0.251189 15.5414
0.316228 17.2618
0.398107 18.9261
0.501187 20.5121
0.630957 21.9438
0.794328 23.0982
1 23.8682

// Results with the DPASS table
0.01 78.6344
0.0125893 86.4868
0.0158489 94.7654
0.0199526 103.126
0.0251189 111.234
0.0316228 118.44
0.0398107 141.519
0.0501187 187.186
0.0630957 245.613
0.0794328 319.3
0.1 412.064
0.125893 529.15
0.158489 677.787
0.199526 865.076
0.251189 1097.72
0.316228 1380.58
0.398107 1710.29
0.501187 2084.68
0.630957 2501.76
0.794328 2952.9
1 3428.19

Obviously something is wrong when I add my DPASS table. Also, I compared the results without my table to the icru73 file:

// icru73/z22_13.dat
0.025 4.12757
 0.03 4.60104
 0.04 5.48804
 0.05  6.3171
 0.06  7.0989
 0.07 7.83895
 0.08 8.54119
 0.09 9.20856
  0.1 9.84352
 0.15 12.5939
  0.2 14.7314
 0.25 16.3811
  0.3 17.6947
  0.4 19.6881
  0.5 21.1346
  0.6 22.1932
  0.7 22.9531
  0.8 23.4845
  0.9 23.8446
    1 24.0777

And it seems that the computation doesn’t give the exact same thing, is it a normal behavior?

Thanks in advance.