Stopping powers dE/dx for ions don't make sense

norach1991 · July 1, 2025, 8:08pm

I have a fairly simple hybrid detector consisting of an ionization chamber filled with isobutane gas and a DSSSD (Double-Sided Silicon Strip Detector) with 300 μm. The ionization chamber has a 0.9 μm Mylar window at its entrance and is segmented into 11 slabs of 1 cm each. The DSSSD has a 0.4 μm aluminum dead layer.

A heavy ion (A < 50) with tens of MeV of energy passes through the system, and I record the energy deposited in each slab. For example, for an 11B ion with 9 MeV of energy, I would expect hundreds of keV to be deposited in each slab and several MeV deposited in the DSSSD. Nevertheless, the ion deposits almost all its energy in the isobutane gas and barely reaches the DSSSD detector.

I wrote a method to print the stopping power values on screen, and when comparing them to LISE++ or SRIM, I noticed significant differences—sometimes up to two orders of magnitude. I am sharing with you the material definitions and PhysicsList implementation. I have also tried changing the PhysicsList (e.g., QGSP_INCLXX_LIV and FTFP_BERT_EMX) and encountered exactly the same issue.

Interestingly, when I simulate an alpha particle, the output of the simulation makes sense.
I’m sharing with you the material and physics definitions, along with a specific portion of the code used to print the dE/dx values on screen.**

Thanks in advance for the help!

PhysicsList.cc (9.0 KB)
DetectorConstruction.cc (40.8 KB)

void SteppingAction::UserSteppingAction(const G4Step* aStep)
{
//get volume of the current step
G4VPhysicalVolume* volume = aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
G4Material* TrackingMaterial = aStep → GetTrack() →
GetVolume() →
GetLogicalVolume() →
GetMaterial();
G4double IonsAlongStep = fElIonPair->SampleNumberOfIonsAlongStep(aStep);

//Get step particle definition
G4ParticleDefinition* stepParticle = aStep-> GetTrack() → GetDefinition();

//collect energy and track length step by step
G4double edep = aStep->GetTotalEnergyDeposit();
G4double edepNonIon = aStep->GetNonIonizingEnergyDeposit();
G4double StepLength = aStep->GetStepLength();

// Kinetic energy of the particle at the beginning of the step
G4double kineticEnergy = aStep->GetPreStepPoint()->GetKineticEnergy();

G4EmCalculator calc;
G4double dedx = calc.ComputeTotalDEDX(
kineticEnergy,
stepParticle,
TrackingMaterial
);

std::cout << "Volume: " << volume->GetName() << " copy " << volume->GetCopyNo() << std::endl;
fDetector->PrintMaterialProperties(TrackingMaterial);
std::cout << "TrackingMaterial: " << TrackingMaterial->GetName() << std::endl;
std::cout << “Kinetic Energy: " << kineticEnergy << " MeV” << std::endl;
std::cout << “edep: " << edep << " MeV” << std::endl;
std::cout << “EdepNonIon : " << edepNonIon << " MeV” << std::endl;
std::cout << “Step length: " << StepLength/cm << " cm” << std::endl;
std::cout << “dE/dx [edep/(StepLength/cm)]: " << (edep/(StepLength/cm)) << " MeV/cm” << std::endl;
std::cout << “dE/dx [G4EMCalculator]: " << dedx*10. << " MeV/cm” << std::endl;
std::cout << “Edep [dedx * StepLength/cm]: " << dedx * StepLength/cm << " MeV” << std::endl;
std::cout << “---------------------------------” << std::endl;

.
.
.
}

Geant4 version: 11.2.2
Linux luispc 6.8.0-60-generic #63~22.04.1-Ubuntu SMP PREEMPT_DYNAMIC Tue Apr 22 19:00:15 UTC 2 x86_64 x86_64 x86_64 GNU/Linux
g++ (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0
cmake version 3.23.0

mcortesg · July 1, 2025, 10:11pm

Dear Luis,

Have you tried a physics list with the _EMZ suffix? Or, equivalently, use G4EmStandardPhysics_option4 from your physics list options? This one incorporates the most accurate combination of models, especially at low energies (~MeV).

Miguel

atolosad · July 3, 2025, 2:52pm

Dear @norach1991

Thank you for your post. To ensure that this is a Geant4 issue, could you please try with the TestEM0 [1], which prints out the cross section and stopping power tables via the EM physics utilities?

Thank you for your time.

Best,
Alvaro
[1] geant4/examples/extended/electromagnetic/TestEm0 at master · Geant4/geant4 · GitHub

norach1991 · July 7, 2025, 7:51am

Hi mcortesg,
Thanks a lot for your help. Please find below an answer to your recommendation.

norach1991 · July 7, 2025, 7:55am

Hi atolosad,
Thanks a lot for your help.

I tried a few things with the TestEM0 example, as you recommended.

In the main() method, choosing a Geant4 PhysicsList package via:

// runManager->SetUserInitialization(new PhysicsList);
G4PhysListFactory factory;
G4VModularPhysicsList* physList = factory.GetReferencePhysList("FTFP_BERT");
physList->ReplacePhysics(new G4EmStandardPhysics_option3());
runManager->SetUserInitialization(physList);

and also with the options:

#/testem/phys/addPhysics local
/testem/phys/addPhysics emstandard_opt0

in the ion.mac macro

always gives me wrong values (off by several orders of magnitude), even if I replace the EM package. For example:

FTFP_BERT
Particle: Li7
Material: G4_Ge
Kinetic Energy: 8 MeV
dE/dx [G4EMCalculator]: 455.119 MeV/mm

FTFP_INCLXX
Particle: Li7
Material: G4_Ge
Kinetic Energy: 8 MeV
dE/dx [G4EMCalculator]: 455.119 MeV/mm

FTFP_BERT_LIV + G4EmStandardPhysics_option3
dE/dx: 245144 MeV/mm

QGSP_BIC + G4EmStandardPhysics_option4
dE/dx: 245144 MeV/mm

FTFP_BERT + G4EmStandardPhysics_option3
dE/dx: 245144 MeV/mm

With dE/dx values calculated via:

G4EmCalculator calc;
G4double dedx_ = calc.ComputeTotalDEDX(energy, particle, material);

std::cout << "Particle: " << partName << std::endl;
std::cout << "Material: " << material->GetName() << std::endl;
std::cout << "Kinetic Energy: " << energy << " MeV" << std::endl;
std::cout << "dE/dx [G4EMCalculator]: " << dedx_ << " MeV/mm" << std::endl;
std::cout << "---------------------------------" << std::endl;

Nevertheless, when I choose:

/testem/phys/addPhysics local
#/testem/phys/addPhysics emstandard_opt0

I get the correct dE/dx values:

Particle: Li7
Material: G4_Ge
Kinetic Energy: 8 MeV
dE/dx [G4EMCalculator]: 0.439829 MeV/mm
Particle: Li7
Material: G4_MYLAR
Kinetic Energy: 8 MeV
dE/dx [G4EMCalculator]: 0.204843 MeV/mm
Particle: Li7
Material: G4_Al
Kinetic Energy: 8 MeV
dE/dx [G4EMCalculator]: 0.301688 MeV/mm
Particle: Li7
Material: G4_Si
Kinetic Energy: 8 MeV
dE/dx [G4EMCalculator]: 0.274226 MeV/mm
Particle: Li7
Material: Isobutane
Kinetic Energy: 8 MeV
dE/dx [G4EMCalculator]: 0.000481997 MeV/mm

So, I guess I have two conclusions from this:

The former prepackaged PhysicsList in Geant4 are not capable of reproducing the dE/dx values correctly by itself.
This is the block of code that guarantees correct dE/dx values in the PhysListEmStandard.cc file:

else if (particleName == "GenericIon") {
  ph->RegisterProcess(new G4hMultipleScattering(), particle);
  G4ionIonisation* ionIoni = new G4ionIonisation();
  ionIoni->SetEmModel(new G4IonParametrisedLossModel());
  ph->RegisterProcess(ionIoni, particle);
  ph->RegisterProcess(new G4NuclearStopping(), particle);
}

atolosad · July 8, 2025, 6:41pm

Hi @norach1991

Thanks a lot for your clear and detailed report.

I tried to crosscheck your findings using the latest Geant4 version (11.4.beta), built with GCC 11 on AlmaLinux 9. I’ll break down my results into two parts.

TestEM0; for 7Li at 8 MeV in G4_Ge
1.a.dEdx ~ 463.5 MeV/mm (EM0), dEdx~ 245.144 GeV/mm (EMZ)
1.b. zero stopping power with local physics lists (I am using the latest version of Geant4 [1]), I cannot reproduce the second part of your findings.

I can reproduce the value (463.5 MeV/mm) found in case a., using the test TestEM5 (with EM0 and EMZ). This test simulates the energy deposited (I setup the thickness to 1 micron). So I guess that TestEM5 and EM0 built-in physics lists give consistent results.

I do not understand why EMZ gives different values in case of TestEM0.

Using builtin physics lists.

In TestEM0, replacing the physics lists by a builtin one, e.g.,

    const G4String plName = "FTFP_BERT_EMZ";
    G4PhysListFactory plFactory;
    G4VModularPhysicsList *pl = plFactory.GetReferencePhysList( plName );
    runManager->SetUserInitialization( pl );

leads to the same numbers as in case a., when the Physics List messenger pass the EM physics lists.

At this point I let the experts on EM physics (@mnovak , @civanch ) to comment.

Best,
Alvaro

civanch · July 9, 2025, 5:07pm

Hello,

Thank you for pointing out to the problem. In TestEm5 for the same parameters all physics configurations provide correct result for local, default, and opt4 physics. In TestEm0 only the default is correct. The reason: G4EmCalculator has problem for both local and opt4 for ions.

This is a bug. Would it be possible to create Bugzilla bug report? Information from this thread is enough - may be copy/paste. More details are not needed - the problem is reproducible.

VI

atolosad · July 10, 2025, 12:07pm

Hi,

I opened a bugzilla ticket,

Please @norach1991 , feel free to add further information if needed

Best,
Alvaro