Difference in particle distribution using fStopAndKill

Particles, Track, Event, Run and Biasing
1

Geant4 Version: geant4-11-02-patch-01
Operating System: Ubuntu 22.04
Compiler/Version: gcc version 11.4.0
CMake Version: cmake version 3.28.4

Dear colleagues,

Currently, I’m facing a problem. I am modelling electromagnetic showers with particle energies above the Cherenkov threshold. To speed up, I applied the fStopAndKill track status for particles whose starting point is below the Cherenkov threshold. As a result, the code began to work faster (which is expected), but the average particle distribution above the threshold along the shower axis has noticeably changed. What could be the reason for this?
Below I attach the code from my SteppingAction.cc. Additionally, I attach a graph of the distributions for 100 GeV electrons I received with and without the fStopAndKill flag.

output
output850×552 23.1 KB

I’d really appreciate any help offered!

#include "SteppingAction.hh"
#include "HDF5Manager.hh"
#include "G4Step.hh"
#include "G4Track.hh"
#include "G4RunManager.hh"
#include "G4SystemOfUnits.hh"

void SteppingAction::UserSteppingAction(const G4Step* step) {
    G4Track* track = step->GetTrack();

    if(!HDF5Manager::Instance()->IsFileOpen()) return;
    
    G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
    
    if(track->GetParticleDefinition()->GetPDGCharge() == 0) return;

    if(step->GetPreStepPoint()->GetVelocity() / CLHEP::c_light < 1. / 1.34) {
        track->SetTrackStatus(fStopAndKill);
        return;
    };

    const G4int currentStep = track->GetCurrentStepNumber();

    if(currentStep == 1) {
        G4StepPoint* prePoint = step->GetPreStepPoint();

        std::vector<double> data_1 = {
        static_cast<double>(prePoint->GetPosition().x() / m),
        static_cast<double>(prePoint->GetPosition().y() / m),
        static_cast<double>(prePoint->GetPosition().z() / m),
        static_cast<double>(prePoint->GetGlobalTime() / ns),
        static_cast<double>(prePoint->GetKineticEnergy() / MeV),
        static_cast<double>(prePoint->GetVelocity() / CLHEP::c_light),
        static_cast<double>(0.0 / m)
    };

    HDF5Manager::Instance()->WriteStepData(
        G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID(),
        track->GetTrackID(), track->GetDefinition()->GetPDGEncoding(),
        data_1
    );

    }
    
    G4StepPoint* postPoint = step->GetPostStepPoint();
    
    std::vector<double> data = {
        static_cast<double>(postPoint->GetPosition().x() / m),
        static_cast<double>(postPoint->GetPosition().y() / m),
        static_cast<double>(postPoint->GetPosition().z() / m),
        static_cast<double>(postPoint->GetGlobalTime() / ns),
        static_cast<double>(track->GetTotalEnergy() / MeV),
        static_cast<double>(track->GetVelocity() / CLHEP::c_light),
        static_cast<double>(step->GetStepLength() / m)
    };

    HDF5Manager::Instance()->WriteStepData(
        G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID(),
        track->GetTrackID(), track->GetDefinition()->GetPDGEncoding(),
        data
    );
}