How to kill particles below a given energy threshold for faster performance

Hi everyone, I want to increase simulation performance by stopping tracking gamma particles, electrons, and positrons (kill them) when their energies drop below some defined energy E (MeV). This doesn’t seem to work with the SetCut() method as I find it just sets a different threshold for secondary particle formation than I would like. So how can I do this?

In SteppingAction :
G4Track* track = (G4Track*)(step->GetTrack());
G4double Ekin = track->GetKineticEnergy();
if (Ekin < Emin) track->SetTrackStatus(fStopAndKill);

1 Like

Thank you for your help. I will try it. Have a good day.

Dear Sir,
I wrote the code below, can you check for me to see if it is correct or not?

maybe you could gain additional performance by replacing the string comparisons with

if (step->GetTrack()->GetParticleDefinition() == G4Gamma::Gamma())

as the stepping action is called for literally every step of the simulation.

Can I do the same with electrons and positrons?

if (step->GetTrack()->GetParticleDefinition() == G4Electron::Electron)
if (step->GetTrack()->GetParticleDefinition() == G4Positron::G4Positron)

sure, that should work perfectly. alternatively:

    G4Track* track = (G4Track*)(step->GetTrack());
    G4double Ekin = track->GetKineticEnergy();
    switch (track->GetParticleDefinition() )
    {
      case G4Gamma::Gamma(): 
        if (Ekin < Emin_gamma) track->SetTrackStatus(fStopAndKill);
        break;
      case G4Electron::Electron():
        /*...*/
        break;
      case G4Positron::G4Positron():
        /*...*/
        break;
    }

In complement to @weller remarks :
1- I realized that cast is not necessary here. It is enough to write :
G4Track* track = step->GetTrack();

2- better to write :
const G4double Emin_gamma = 10*keV; etc …
with const, Emin_gamma is initialized only once and never modified.

Dear @maire and @weller
I have tried compiling the code below. But it shows some errors that I don’t understand.


there is two issues, the easy one:

#include "G4Gamma.hh"
#include "G4Electron.hh"
#include "G4Positron.hh"

the hard one: apparently the switch-case construction does not accept track->GetParticleDefinition(). Sorry I brought this up, thought this would work. Apparently, switch-case is more than just replacing a series of if-elseif-elseif stuff :frowning:

Easiest is to revert to what you suggested (but you still have to include the classes):

    const G4ParticleDefinition* particleDefinition = track->GetParticleDefinition();
    if (particleDefinition == G4Electron::Electron())
     {/* ... */}
    else if (particleDefinition == G4Positron::Positron()) 
     {/* ... */}

Thank you very much, I understood the problem and fixed the code according to your instructions. I wish you a good day.

Another point : as the name of the function tell, Ekin is the kinetic energy of the particle, not its total energy.
Do you really want to kill electrons of kinetic energy = 521 keV ?

Dear Sir, I am simulating Geant4 and comparing the results with MCNP. In MCNP I use the CUT:p and CUT:e cards to set an energy threshold for tracking individual particles. So I just want to put the same condition for two simulation cases.