Monitor all particle tracks at last time step

Particles, Track, Event, Run and Biasing
1

_Geant4 Version:11.2.2
_Operating System:Ubuntu 24
_Compiler/Version:gcc 12.4
_CMake Version:cmake 3.28

#include "EventAction.hh"
#include "G4RunManager.hh"
#include "RunAction.hh"
#include "G4Event.hh"
#include "G4PrimaryVertex.hh"
#include "G4PrimaryParticle.hh"
#include "G4ParticleDefinition.hh"
#include "G4ThreeVector.hh"


namespace B4d
{
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

EventAction::EventAction()
: G4UserEventAction(), fRunAction(nullptr), fInitialPhaseSpaceRecorded(false)
{
  // initialization per event
  fRunAction = const_cast<RunAction*>(static_cast<const RunAction*>(G4RunManager::GetRunManager()->GetUserRunAction()));
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

EventAction::~EventAction()
{
}

void EventAction::BeginOfEventAction(const G4Event* event)
{
 
    
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void EventAction::EndOfEventAction(const G4Event* event)
{
  for (G4int i = 0; i < event->GetNumberOfPrimaryVertex(); ++i) {
        G4PrimaryVertex* vertex = event->GetPrimaryVertex(i);
        
        for (G4int j = 0; j < vertex->GetNumberOfParticle(); ++j) {
            G4PrimaryParticle* particle = vertex->GetPrimary(j);
            
            // Retrieve particle information
            G4double x = vertex->GetX0();
            G4double y = vertex->GetY0();
            G4double z = vertex->GetZ0();
            G4double px = particle->GetMomentum().x();
            G4double py = particle->GetMomentum().y();
            G4double pz = particle->GetMomentum().z();
            G4int particleID = particle->GetPDGcode(); // or another method to get an ID
            G4double timestep = 0. /* Your logic to get the timestep */;

            // Fill the ntuple with the retrieved data
            fRunAction->FillPhaseSpace(x, y, z, px, py, pz, particleID, timestep);
          }
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
void EventAction::RecordStepData(const G4ThreeVector& position, const G4ThreeVector& momentum, G4int particleID, G4double timestep)
{
  // Assuming fRunAction has a method to process step data
  fRunAction->FillPhaseSpace(position.x(), position.y(), position.z(), momentum.x(), momentum.y(), momentum.z(), particleID, timestep);
}

}

#include "RunAction.hh"
#include "G4AnalysisManager.hh"
#include "G4Run.hh"
#include "G4RunManager.hh"
#include "G4UnitsTable.hh"
#include "G4SystemOfUnits.hh"

namespace B4d
{

RunAction::RunAction()
: G4UserRunAction()
{
    G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
    analysisManager->SetVerboseLevel(1);
    analysisManager->SetFileName("phasesapce.root");
}

void RunAction::BeginOfRunAction(const G4Run*)
{
    
    G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
    analysisManager->SetNtupleMerging(true);
    analysisManager->OpenFile();
    // Open the output file directly from the macro
    G4cout << "Opening output file from macro." << G4endl;

    // Create a tree and define branches
    analysisManager->CreateNtuple("phaseSpace", "Phase Space Data");
    analysisManager->CreateNtupleDColumn("x");
    analysisManager->CreateNtupleDColumn("y");
    analysisManager->CreateNtupleDColumn("z");
    analysisManager->CreateNtupleDColumn("px");
    analysisManager->CreateNtupleDColumn("py");
    analysisManager->CreateNtupleDColumn("pz");
    analysisManager->CreateNtupleIColumn("particleID");
    analysisManager->CreateNtupleDColumn("timestep");
    analysisManager->CreateNtupleIColumn("timestepID");

    analysisManager->FinishNtuple();
    
    // Reset timestep ID at the start of each run
    fTimestepID = 0;
}

void RunAction::EndOfRunAction(const G4Run*)
{
    G4AnalysisManager::Instance()->Write();
    G4AnalysisManager::Instance()->CloseFile();

}

void RunAction::FillPhaseSpace(G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz, G4int particleID, G4double timestep)
{
    G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
    
    analysisManager->FillNtupleDColumn(0, x);
    analysisManager->FillNtupleDColumn(1, y);
    analysisManager->FillNtupleDColumn(2, z);
    analysisManager->FillNtupleDColumn(3, px);
    analysisManager->FillNtupleDColumn(4, py);
    analysisManager->FillNtupleDColumn(5, pz);
    analysisManager->FillNtupleIColumn(6, particleID);
    analysisManager->FillNtupleDColumn(7, timestep);
    analysisManager->FillNtupleIColumn(8, fTimestepID);
    
    analysisManager->AddNtupleRow();
    
    fTimestepID++;  // Increment timestep ID for each call
}
}

#include "EventAction.hh"
#include "G4RunManager.hh"
#include "RunAction.hh"
#include "G4Event.hh"
#include "G4PrimaryVertex.hh"
#include "G4PrimaryParticle.hh"
#include "G4ParticleDefinition.hh"
#include "G4ThreeVector.hh"


namespace B4d
{
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

EventAction::EventAction()
: G4UserEventAction(), fRunAction(nullptr), fInitialPhaseSpaceRecorded(false)
{
  // initialization per event
  fRunAction = const_cast<RunAction*>(static_cast<const RunAction*>(G4RunManager::GetRunManager()->GetUserRunAction()));
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

EventAction::~EventAction()
{
}

void EventAction::BeginOfEventAction(const G4Event* event)
{
 
    
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void EventAction::EndOfEventAction(const G4Event* event)
{
  for (G4int i = 0; i < event->GetNumberOfPrimaryVertex(); ++i) {
        G4PrimaryVertex* vertex = event->GetPrimaryVertex(i);
        
        for (G4int j = 0; j < vertex->GetNumberOfParticle(); ++j) {
            G4PrimaryParticle* particle = vertex->GetPrimary(j);
            
            // Retrieve particle information
            G4double x = vertex->GetX0();
            G4double y = vertex->GetY0();
            G4double z = vertex->GetZ0();
            G4double px = particle->GetMomentum().x();
            G4double py = particle->GetMomentum().y();
            G4double pz = particle->GetMomentum().z();
            G4int particleID = particle->GetPDGcode(); // or another method to get an ID
            G4double timestep = 0. /* Your logic to get the timestep */;

            // Fill the ntuple with the retrieved data
            fRunAction->FillPhaseSpace(x, y, z, px, py, pz, particleID, timestep);
          }
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
void EventAction::RecordStepData(const G4ThreeVector& position, const G4ThreeVector& momentum, G4int particleID, G4double timestep)
{
  // Assuming fRunAction has a method to process step data
  fRunAction->FillPhaseSpace(position.x(), position.y(), position.z(), momentum.x(), momentum.y(), momentum.z(), particleID, timestep);
}

}

Above are 3 classes I modified to track particle momentums at various timesteps. I am basically getting the initial distributions and any track that is updated based on some pre-defined geant criteria. It seems that the only time I access the step data is when the particles momentum changes {ie crosses a boundary}

How can I modify this to retrieve every particles info at the very last time step?

I assume I need some sort of maximumTimestep variable that I keep track off. How do I guarantee that I retrieve every particle’s track though? Another thought I had was maybe to just add a post processing step to propagate all particles to their respective positions in the final recorded timestep since the particle information is just a straight line trajectory from the final track information.

2

Hello,

Probably the “action” that you want is the SteppingAction… this is execute for every step which is the one updating the Track… that being said a couple of things come to mind to achieve what you want.

  • You could keep track (pun unintended) of ALL steps overwriting the previous one to simply keep the latest and therefore potentially the last one, OR
  • MAYBE you could check for the TrackStatus after the Step. If the step leaves the track fStopButAlive, assuming you have no AtRest processes (like Radioactive Decay), this could be the last one. However it is early, I haven´t had coffee… this may not makes sense at all.

/Pico

3

Hi Pico,

Thank you for your reply. I actually solved the issue myself. I have a stepping action class. The default output from geant is to output the particles when they have a change in trajectory that will be significant enough to change their next collision boundary (from what I have read). So in my setup the particles get tracked at various boundaries as I want them to. I decided the best course of action was to just progress the particles based on their measured momentum to the maximum timestep observed. This indeed did what I needed to do. I also used the GetTrackID and GetTrackStatus methods to assure that the particles where the ones I was expecting. Thanks for your reply!