Double Counting electrons in SteppingAction

kaychi · August 28, 2023, 3:47am

Hi, so
I have this loop in my stepping action to try and record the electrons in my simulation. Unfortunately It appears to be double counting the number of electrons:

//SteppingAction.hh
struct ElectronTrack {
    G4int trackID;
    G4double totalLength;
    bool counted;

    ElectronTrack(G4int id) : trackID(id), totalLength(0.0), counted(false) {}
};
private:
  std::map<G4int, ElectronTrack> electronTracks;

//SteppingAction.cc




void SteppingAction::UserSteppingAction(const G4Step* aStep)
{

	...<other code>...
      ElectronTrack& eleTrack = electronTracks[trackID];

if (aStep->GetTrack()->GetDynamicParticle()->GetDefinition() ->GetParticleName() == "e-"
			){
	/* G4cout << "ELECTRON"<<G4endl; */
  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
		G4int eventid = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
        G4int trackID = track->GetTrackID();


        G4double stepLength = step->GetStepLength();
        eleTrack.totalLength += stepLength;

        // Check for exit or absorption
        if (track->GetTrackStatus() == fStopAndKill) {
                G4double trackLength = eleTrack.totalLength;
                // Process the track length as needed

				G4ThreeVector position = track->GetPosition();
				G4double energy = track->GetTotalEnergy();
				G4double timing = track->GetGlobalTime();
				G4ThreeVector mom =track->GetMomentum();
				G4double tracklength = track->GetTrackLength();

	G4cout << "ELECTRON, track length = "<< G4BestUnit(trklength,"Length")<< G4endl;

		  analysisManager->FillNtupleIColumn(1,0,eventid);
		  analysisManager->FillNtupleIColumn(1,1,trackID);
		  analysisManager->FillNtupleDColumn(1,2,position[0]);
		  analysisManager->FillNtupleDColumn(1,3,position[1]);
		  analysisManager->FillNtupleDColumn(1,4,position[2]);
		  analysisManager->FillNtupleDColumn(1,5,energy);
		  analysisManager->FillNtupleDColumn(1,6,timing);
		  analysisManager->FillNtupleDColumn(1,7,mom[0]);
		  analysisManager->FillNtupleDColumn(1,8,mom[1]);
		  analysisManager->FillNtupleDColumn(1,9,mom[2]);
		  analysisManager->FillNtupleDColumn(1,10,tracklength);
		  analysisManager->AddNtupleRow(1);
                 }
            }
	}

But it a appears that ~10 x the number of electrons as expected from theory are being stored. I think by counting via the steps without a check for repeating steps might explain this. Any advice is greatly appreciated! Thanks!

_Geant4 Version:_11.0.0
_Operating System:_MacOS
_Compiler/Version:_Apple clang version 14.0.0
_CMake Version:_3.25.1

mkelsey · August 28, 2023, 4:15am

Electrons, like most particles, can take multiple steps before losing all of their energy and being killed by the tracker. You may want to do your counting just for step #1, or keep a buffer of track ID’s and only count each electron once per event (reset the buffer at each new event).

Also, it’s better to test pointers than strings:

if (track->GetDefinition() == G4Electron::Definition()) {
[...]
}

kaychi · August 29, 2023, 12:53am

Thanks so much for the information @mkelsey! If I only count step #1 however, I wouldn’t also be able to count the total track length of the electrons right? how could I guarantee catching the electrons at the end of the event before I record their information? Am I correct if I use if (track->GetTrackStatus() == fStopAndKill) I can retrieve them?

thanks for any advice!

PS After thinnking a little bit you meant something like

 if (track->GetCurrentStepNumber() == 1)

I guess? Thankns!

maire · August 29, 2023, 8:05am

In TrackingAction::PostUserTrackingAction(), you have the particle at end of its tracking

kaychi · August 30, 2023, 8:56am

Thanks @maire for the tip! Essentially though I want to record all the electrons timing, position and energy information, as well as the track length at every step in the simulation, so maybe the tracking action is not the best suited place for this? Or at least, could I split this information returned at each step as well as the final electrons at end of the tracking using the TrackingAction?
Thanks!

maire · August 30, 2023, 1:03pm

If you want informations step per step, you must do every things in SteppingAction …
If you run 1 event with /tracking/verbose 2, you will see what kind of informations you have step per step

kaychi · August 31, 2023, 12:19am

Hi @maire, thanks, but I want to record the data in an n-tuple. In any case I managed to start recording it via the stepping action, but there are some odd features, here is an example record from a step:

//From User Stepping Action
  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
  G4Track* track = aStep->GetTrack();
  G4ParticleDefinition* particleDef = track->GetDefinition();

  /*   // Check if the particle is an electron */
    if (particleDef == G4Electron::Definition()) {
		/* G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); */
        const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
        const G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
        const G4VProcess* process = postStepPoint->GetProcessDefinedStep();
		G4String processName = process->GetProcessName();

G4ParticleDefinition* particle = track->GetDefinition(); // We get the particle definition
G4String name = particle->GetParticleName(); // The particle name


            ElectronKinematicInfo kinematicInfo;
            kinematicInfo.eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
			kinematicInfo.position = track->GetPosition();
            kinematicInfo.momentum = preStepPoint->GetMomentum();
            kinematicInfo.energy = track->GetTotalEnergy();
            kinematicInfo.totalTrackLength = track->GetTrackLength();
            kinematicInfo.stepNumber = track->GetCurrentStepNumber();
            kinematicInfo.time = track->GetGlobalTime();
			kinematicInfo.interaction = processName;
		    kinematicInfo.name = name;


  /*           // Output kinematic information to console */
            G4cout << kinematicInfo.name << G4endl;

			G4cout << "position: " << kinematicInfo.position << G4endl;
            G4cout << "momentum: " << kinematicInfo.momentum << G4endl;
            G4cout << "energy: " << kinematicInfo.energy << G4endl;
            G4cout << "total track length: "
                   << G4BestUnit(kinematicInfo.totalTrackLength, "Length") << G4endl;
            G4cout << "step number: " << kinematicInfo.stepNumber << G4endl;
            G4cout << "time: " << G4BestUnit(kinematicInfo.time, "Time") << G4endl;
			G4cout << "Track ID: " << kinematicInfo.trackID << G4endl;
			G4cout << "Interaction: "<< processName<< G4endl;

e-
position: (-122.324,-279.039,-67.6839)
momentum: (0.104076,0.131221,0.0138439)
energy: 0.537835
total track length: 1.23784 um
step number: 1
time: 3.83792 ns
Track ID: 1563486488
Interaction: msc

So the Track ID and Interaction seem a bit strange? Also, I dumped the data to an N-tuple, and analysed the number of unique Event IDs against their Number of Steps per run.

kaychi · August 31, 2023, 12:22am

After looking a little bit deeper, I found that using the stepping is giving some odd output. After shooting one particle into my simulation, this is what it is recording (this appears for one electron I believe)

    fEvent    fTrackID  StepNumber       fX       fY       fZ  eEnergy   Time  \
0        0  1563486488           1 -181.289 -273.153 -126.052    0.511  5.007   
1        0           1           2 -181.289 -273.153 -126.052    0.511  5.007   
2        0  1563486488           1 -181.291 -273.154 -126.054    0.511  5.007   
3        0  1563486488           2 -181.291 -273.154 -126.054    0.511  5.007   
4        0  1563486488           3 -181.291 -273.154 -126.054    0.511  5.007   
5        0  1563486488           1 -181.288 -273.152 -126.052    0.537  5.007   
6        0  1563486488           2 -181.281 -273.149 -126.052    0.518  5.007   
7        0  1563486488           3 -181.281 -273.149 -126.052    0.511  5.007   
8        0           3           4 -181.281 -273.149 -126.052    0.511  5.007

The trackID changes with some of these steps appearing with interaction type msc?

kaychi · August 31, 2023, 1:03am

OK I fixed the trackID problem (actually it was not storing anything xD so I added)

kinematicInfo.trackID = track->GetTrackID();

and seems to be resolved. Now I just need to confirm the count…

kaychi · September 5, 2023, 7:33am

So I added to my Stepping action the following and it is working well:

void SteppingAction::UserSteppingAction(const G4Step* aStep)
{


	  G4Track* track = aStep->GetTrack();
	  G4ParticleDefinition* particleDef = track->GetDefinition();
	  if (particleDef == G4Gamma::Definition() ){
		  GatherInformation(1,aStep);
	  } else if (particleDef == G4Electron::Definition() ) {
		  GatherInformation(2,aStep);
	  } else if (particleDef == G4Positron::Definition() ) {
		  GatherInformation(3,aStep);
	  } else if (particleDef == G4OpticalPhoton::Definition() )
		  GatherInformation(4,aStep);
}


void SteppingAction::GatherInformation(const G4int histNo, const G4Step* aStep){


	  ElectronKinematicInfo kinematicInfo;
	  G4Track* track = aStep->GetTrack();

    // Check if 'track' is not null before using it.
    if (!track) {
        G4cerr << "Error: Null track in UserSteppingAction." << G4endl;
        return;
    }

		const G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
        const G4StepPoint* postStepPoint = aStep->GetPostStepPoint();
		const G4VProcess* process = track->GetCreatorProcess();

		if (preStepPoint->GetProcessDefinedStep())
			{
				const G4String& processName = preStepPoint->GetProcessDefinedStep()->GetProcessName();
			kinematicInfo.interaction = processName;
			}

		G4ParticleDefinition* particle = track->GetDefinition(); // We get the particle definition
		G4String name = particle->GetParticleName(); // The particle name


            kinematicInfo.eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
			kinematicInfo.positionPre = preStepPoint->GetPosition();
			kinematicInfo.positionPost = postStepPoint->GetPosition();
            kinematicInfo.momentum = preStepPoint->GetMomentum();
            kinematicInfo.energy = track->GetTotalEnergy();
            kinematicInfo.totalTrackLength = track->GetTrackLength();
            kinematicInfo.stepNumber = track->GetCurrentStepNumber();
            kinematicInfo.time = track->GetGlobalTime();
		    kinematicInfo.name = name;
		    kinematicInfo.trackID= track->GetTrackID();

            G4cout << kinematicInfo.name << G4endl;

			G4cout << "position: " << kinematicInfo.positionPre << G4endl;
            G4cout << "momentum: " << kinematicInfo.momentum << G4endl;
            G4cout << "energy: " << kinematicInfo.energy << G4endl;
            G4cout << "total track length: "
                   << G4BestUnit(kinematicInfo.totalTrackLength, "Length") << G4endl;
            G4cout << "step number: " << kinematicInfo.stepNumber << G4endl;
            G4cout << "time: " << G4BestUnit(kinematicInfo.time, "Time") << G4endl;
            G4cout << "interaction: " << kinematicInfo.interaction << G4endl;
			G4cout << "Track ID: " << kinematicInfo.trackID << G4endl;

	  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();

		  analysisManager->FillNtupleIColumn(histNo,0,kinematicInfo.eventID);
		  analysisManager->FillNtupleIColumn(histNo,1,kinematicInfo.trackID);
		  analysisManager->FillNtupleIColumn(histNo,2,kinematicInfo.stepNumber);
		  analysisManager->FillNtupleDColumn(histNo,3,kinematicInfo.positionPre[0]);
		  analysisManager->FillNtupleDColumn(histNo,4,kinematicInfo.positionPre[1]);
		  analysisManager->FillNtupleDColumn(histNo,5,kinematicInfo.positionPre[2]);
		  analysisManager->FillNtupleDColumn(histNo,6,kinematicInfo.energy);
		  analysisManager->FillNtupleDColumn(histNo,7,kinematicInfo.time);
		  analysisManager->FillNtupleDColumn(histNo,8,kinematicInfo.momentum[0]);
		  analysisManager->FillNtupleDColumn(histNo,9,kinematicInfo.momentum[1]);
		  analysisManager->FillNtupleDColumn(histNo,10,kinematicInfo.momentum[2]);
		  analysisManager->FillNtupleDColumn(histNo,11,kinematicInfo.totalTrackLength);
		  analysisManager->FillNtupleSColumn(histNo,12,kinematicInfo.name);
		  analysisManager->FillNtupleSColumn(histNo,13,kinematicInfo.interaction);
		  analysisManager->FillNtupleDColumn(histNo,14,kinematicInfo.positionPost[0]);
		  analysisManager->FillNtupleDColumn(histNo,15,kinematicInfo.positionPost[1]);
		  analysisManager->FillNtupleDColumn(histNo,16,kinematicInfo.positionPost[2]);
		  analysisManager->AddNtupleRow(histNo);
}

However, it is now undercounting electrons according to the Run action aggregate information (I am taking the Hadr04 Run class to count all particles in my simulation).

Why is the number of electrons counted from the stepping action less than the aggregated information from the Run information? Thanks for any information!