Error when simulating electron at low energy

tommasorizzi · June 16, 2023, 3:34pm

Hello everyone,

I am using GEANT4 to have some backscattering data on different materials at different energies etc… So I basically shoot electrons against a material to obtain the backscattering yield.

I build my own physics list (attached later) but I have a problem when decreasing energy below 1 keV. In fact, above 1 keV it behaves correctly and the results I obtained are in agreement with literature.

Below 1 keV , instead, I can see no outcome of my simulations, and the electrons do not even interact with the materials. I think the problem is in my physics list, because using the EmStandardSS provided by Geant4 everything works fine.

Here my physics list:

void MyPhysicsList::ConstructEM()
{
G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();

G4ParticleDefinition* ion = G4GenericIon::GenericIon();
G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
G4ParticleDefinition* electron =   G4Electron::ElectronDefinition();
G4ParticleDefinition* positron =  G4Positron::PositronDefinition();

//gamma
G4RayleighScattering *rs = new G4RayleighScattering;
rs ->SetEmModel(new G4PenelopeRayleighModel()); //Penelope

G4ComptonScattering *cs = new G4ComptonScattering;
cs ->SetEmModel(new G4PenelopeComptonModel()); //Penelope

G4PhotoElectricEffect *pe = new G4PhotoElectricEffect;
pe ->SetEmModel(new G4PenelopePhotoElectricModel()); //penelope

G4GammaConversion *gc = new G4GammaConversion;
gc ->SetEmModel(new G4PenelopeGammaConversionModel()); //penelope

ph->RegisterProcess(new G4PhotoElectricEffect(), gamma);
ph->RegisterProcess(new G4ComptonScattering(), gamma);
ph->RegisterProcess(new G4GammaConversion(), gamma);
ph->RegisterProcess(new G4RayleighScattering(), gamma);


//electron

//multiple scattering

G4eMultipleScattering *msc = new G4eMultipleScattering;
msc->SetEmModel(new G4UrbanMscModel()); //G4GoudsmitSaundersonMscModel -

G4UrbanMscModel - G4WentzelVIModel

G4eIonisation *eion = new  G4eIonisation;
eion -> SetEmModel(new  G4PenelopeIonisationModel()); //Penelope

G4eBremsstrahlung *ebr = new G4eBremsstrahlung;
ebr -> SetEmModel ( new G4PenelopeBremsstrahlungModel() ); //penelope

ph->RegisterProcess(new G4eMultipleScattering(), electron);
ph->RegisterProcess(new G4eBremsstrahlung(), electron);
ph->RegisterProcess(new G4eIonisation(), electron);

//positron

G4eMultipleScattering *mscp = new G4eMultipleScattering;
mscp->SetEmModel(new G4UrbanMscModel ());

G4eIonisation *epion = new  G4eIonisation;
epion -> SetEmModel(new  G4PenelopeIonisationModel()); //penelope

G4eBremsstrahlung *epbr = new G4eBremsstrahlung;
epbr -> SetEmModel ( new G4PenelopeBremsstrahlungModel() ); //penelope

G4eplusAnnihilation *eann = new G4eplusAnnihilation;
eann -> SetEmModel (new G4PenelopeAnnihilationModel()); //penelope

ph->RegisterProcess(new G4eMultipleScattering(), positron);
ph->RegisterProcess(new G4eBremsstrahlung(), positron);
ph->RegisterProcess(new G4eIonisation(), positron);
ph->RegisterProcess(new G4eplusAnnihilation(), positron);

}

Does someone know how to help? Thank you

guatelli · June 16, 2023, 9:45pm

Hi, I see you use the PENELOPE EM models and you should be able to model particle interactions down to 100 eV. Can you please add in your Physics List constructor the following:
G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(100eV, 1GeV);
if you do not add this line, the interactions are modelled down to 1 keV. My impression is that this is the problem. Can you please try?

tommasorizzi · June 19, 2023, 5:08pm

Hi!

I think the problem you are mentioning is correct. In fact I can’t see any electrons production below 1 keV.

That’s weird by the way, because I call the specific method SetCuts to establish a cut off of the order of 0.1 nm that should correspond to around 50 eV.
Anyway, I tried to add you command in my constructor, but still not working.

I attach here my complete physics list. Do you know which might be the problem?

MyPhysicsList::MyPhysicsList()
{

G4EmParameters * param = G4EmParameters::Instance();
param->SetApplyCuts(true);
SetDefaultCutValue(0.01*nm);
SetVerboseLevel(2);
//Deexcitation
G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
de->SetFluo(true);
de->SetAuger(true);
de->SetPIXE(true);
de->SetAugerCascade(true);

G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(1*eV, 1*GeV);

}

MyPhysicsList::~MyPhysicsList()
{}

// //define all particles: e-, e+, n, gamma, p
void MyPhysicsList::ConstructParticle()
{
G4Electron::ElectronDefinition();
G4Gamma::GammaDefinition();
G4Positron::PositronDefinition();
G4Proton::ProtonDefinition();
G4ParticleDefinition* particle = G4GenericIon::GenericIon();

}
void MyPhysicsList::ConstructProcess()
{

AddTransportation();

ConstructEM();

G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(1*eV, 1*GeV);

}

void MyPhysicsList::ConstructEM()
{
G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();

G4ParticleDefinition* ion = G4GenericIon::GenericIon();
G4ParticleDefinition* gamma = G4Gamma::GammaDefinition();
G4ParticleDefinition* electron =  G4Electron::ElectronDefinition();
G4ParticleDefinition* positron =  G4Positron::PositronDefinition();

//gamma

G4RayleighScattering *rs = new G4RayleighScattering;
rs ->SetEmModel(new G4PenelopeRayleighModel()); //Penelope
//rs ->SetEmModel(new G4LivermoreRayleighModel()); //Livermore

G4ComptonScattering *cs = new G4ComptonScattering;
cs ->SetEmModel(new G4PenelopeComptonModel()); //Penelope
//cs ->SetEmModel(new G4LivermoreComptonModel()); //livermore

G4PhotoElectricEffect *pe = new G4PhotoElectricEffect;
pe ->SetEmModel(new G4PenelopePhotoElectricModel()); //penelope
//pe ->SetEmModel(new G4LivermorePhotoElectricModel()); //livermore

G4GammaConversion *gc = new G4GammaConversion;
gc ->SetEmModel(new G4PenelopeGammaConversionModel()); //penelope
//gc ->SetEmModel(new G4LivermoreGammaConversionModel()); //livermore


ph->RegisterProcess(new G4PhotoElectricEffect(), gamma);
ph->RegisterProcess(new G4ComptonScattering(), gamma);
ph->RegisterProcess(new G4GammaConversion(), gamma);
ph->RegisterProcess(new G4RayleighScattering(), gamma);

//electron

//multiple scattering

G4eMultipleScattering *msc = new G4eMultipleScattering;
msc->SetEmModel(new G4UrbanMscModel()); //G4GoudsmitSaundersonMscModel - G4UrbanMscModel - G4WentzelVIModel

G4EmParameters* param = G4EmParameters::Instance();


G4eIonisation *eion = new  G4eIonisation;
eion -> SetEmModel(new  G4PenelopeIonisationModel()); //Penelope
//eion -> SetEmModel(new  G4LivermoreIonisationModel); //Livermore

G4eBremsstrahlung *ebr = new G4eBremsstrahlung;
ebr -> SetEmModel ( new G4PenelopeBremsstrahlungModel() ); //penelope
//ebr -> SetEmModel ( new G4LivermoreBremsstrahlungModel ); //Livermore

ph->RegisterProcess(new G4eMultipleScattering(), electron);
ph->RegisterProcess(new G4eBremsstrahlung(), electron);
ph->RegisterProcess(new G4eIonisation(), electron);

//positron

G4eMultipleScattering *mscp = new G4eMultipleScattering;
mscp->SetEmModel(new G4UrbanMscModel());


G4eIonisation *epion = new  G4eIonisation;
epion -> SetEmModel(new  G4PenelopeIonisationModel()); //penelope

G4eBremsstrahlung *epbr = new G4eBremsstrahlung;
epbr -> SetEmModel ( new G4PenelopeBremsstrahlungModel() ); //penelope

G4eplusAnnihilation *eann = new G4eplusAnnihilation;
eann -> SetEmModel (new G4PenelopeAnnihilationModel()); //penelope

ph->RegisterProcess(new G4eMultipleScattering(), positron);
ph->RegisterProcess(new G4eBremsstrahlung(), positron);
ph->RegisterProcess(new G4eIonisation(), positron);
ph->RegisterProcess(new G4eplusAnnihilation(), positron);

G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(1*eV, 1*GeV);

}

void MyPhysicsList::SetCuts()
{
G4VUserPhysicsList::SetCuts();

G4double NewCutValue=0.01*nm;
G4VUserPhysicsList::SetDefaultCutValue(NewCutValue);


G4VUserPhysicsList::SetCutValue(NewCutValue, "gamma");

G4VUserPhysicsList::SetCutValue(NewCutValue, "e+");

G4VUserPhysicsList::SetCutValue(NewCutValue, "e-");

}

Thank you for any help.

guatelli · June 19, 2023, 10:21pm

My suggestion is to use a pre-build EM physics list of Geant4 to check if you observe the same behaviour. For example, look at the physics list of the brachytherapy advanced examples (keep the commands that you included in the constructor). Let me know how it goes,
cheers
Susanna

civanch · June 21, 2023, 8:37am

Hello,

if you need to simulate tracking below 1 keV better to do following:

use EM physics Option4, if you have an application with only EM physics - use G4EmStandardPhysics_option4, if you are using full physics - use predefined QGSP_BIC_EMZ .
for transport of low-energy e- it is needed to reduce tracking cut on electrons, use UI command (there is equivalent C++ interface):

/process/em/lowestElectronEnergy 10 eV

also you may allow to reduce energy threshold for delta electron production as Susanna proposed (there is equivalent UI command): G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(100eV, 1 GeV);

VI

tommasorizzi · July 3, 2023, 2:16pm

Hello,

thank you both for your help and sorry for this late update.

Since I need to benchmark my precise physicsList I would not mind to use the EM physics Option4, which is working properly in my code. However, when I try to use my physics List with Penelope models I still get that no electrons are simulated. I tried to replace Penelope libraries with Livermore, and in that case it works fine.

Is there a specific reason why Penelope processes do not work below 1 KeV? Do I need to add some more specifics?
I also added the tracking cut on electrons as you suggested…

Thanks again for helping,
TR