How to Track Energy transferd by Primaries to Secondaries

Particles, Track, Event, Run and Biasing
1

_Geant4 Version:_11.2.1
_Operating System:_Ubuntu
_Compiler/Version:_11.4.0
_CMake Version:_3.22.1
Hi everyone,
I have a layered detector consisting of two different materials: BGO and EJ232 (Plastic). I am simulating the interaction of 511 keV gamma photons with this detector and tracking the energy deposition. Specifically, I want to track the kinetic energy that a primary photon transfers to secondary particles through both the photoelectric and Compton processes in each event.

However, I am encountering an issue: if a primary photon undergoes both Compton scattering and the photoelectric effect within the same event, my current tracking implementation only records the kinetic energy transferred during the first process (typically Compton scattering). As a result, I am unable to account for the total energy transferred by the primary photon to secondaries, which prevents me from achieving energy balance at the end of the simulation.

I have attached a code snippet of my tracking logic along with terminal output that illustrates the problem. In the output, you can see that the energy transfer is only recorded once, even when multiple interactions occur in the same event. How can I modify my tracking approach to ensure that the total energy transferred by the primary photon to secondaries is correctly recorded for each event?

Screenshot from 2025-02-28 15-23-37
Screenshot from 2025-02-28 15-23-371444×543 107 KB

G4StepPoint* preStep = step->GetPreStepPoint();
	G4StepPoint* postStep = step->GetPostStepPoint();

	// Get kinetic energy and energy deposition
	G4double kinEnergy = preStep->GetKineticEnergy();
	//G4double edep = step->GetTotalEnergyDeposit();

	// Get layer ID and material
	auto touchable = step->GetPreStepPoint()->GetTouchable();
	G4int layerID = touchable->GetVolume(1)->GetCopyNo();
	G4Material* material = step->GetPreStepPoint()->GetMaterial();
	if (layerID < 0 || layerID >= fNofCells) return false;

	// Track primary gamma interactions specifically
	if (particleName == "gamma" && parentID == 0) {

 	auto edep = step->GetTotalEnergyDeposit();
        G4double stepLength = step->GetStepLength();
	
	G4double preStepEnergy  = step->GetPreStepPoint()->GetTotalEnergy();
        G4double postStepEnergy = step->GetPostStepPoint()->GetTotalEnergy();

    // Initialize tracker for new primary gammas - only once per track
    if (fPrimaryGammaTracker.find(trackID) == fPrimaryGammaTracker.end()) {
        PrimaryGammaInfo info;
        info.initialEnergy = preStepEnergy;
        info.isFirstEntry = true;
        info.initialMaterial = material->GetName();
        info.initialLayerID = layerID;
        info.hasExited = false;
        info.depositedEnergy = 0.0;
        info.transferredToSecondaries = 0.0;
        info.exitEnergy = 0.0;
        info.trackedSecondaries = std::set<G4int>();

        fPrimaryGammaTracker[trackID] = info;

        // Print initial energy only once
    /*    G4cout << "Primary gamma entered detector: Track ID = " << trackID
               << ", Initial Energy = " << preStepEnergy / CLHEP::keV << " keV" << G4endl;*/
    }
    
    G4cout << "Step: PreStep Volume = " << preStep->GetPhysicalVolume()->GetName()
       << ", PostStep Volume = " << (postStep->GetPhysicalVolume() ? postStep->GetPhysicalVolume()->GetName() : "nullptr")
       << ", Step Status = " << postStep->GetStepStatus()
       << ", PreStep Energy = " << preStepEnergy / CLHEP::keV << " keV"
       << ", PostStep Energy = " << postStepEnergy / CLHEP::keV << " keV" << G4endl;

    // Add deposited energy to the tracker
    fPrimaryGammaTracker[trackID].depositedEnergy += edep;

     // Check for detector exit (when particle leaves world volume)
    if (postStep->GetStepStatus() == fWorldBoundary || 
        postStep->GetPhysicalVolume() == nullptr ||
        postStep->GetPhysicalVolume()->GetName() == "OutOfWorld") {

        // Only update exit info if it hasn't been marked as exited yet
        if (!fPrimaryGammaTracker[trackID].hasExited) {
            fPrimaryGammaTracker[trackID].exitEnergy = preStepEnergy;
            fPrimaryGammaTracker[trackID].hasExited = true;

            // Debug: Print exit energy
            G4cout << "Primary gamma exited detector: Track ID = " << trackID
                   << ", Exit Energy = " << preStepEnergy / CLHEP::keV << " keV" << G4endl;
        }
    }

    else if (postStep->GetPhysicalVolume()->GetName() != "BGO" && 
             postStep->GetPhysicalVolume()->GetName() != "EJ232") {
        // Gamma is transitioning between detector components
      /*  G4cout << "Primary gamma transitioning: Track ID = " << trackID
               << ", Current Energy = " << postStepEnergy / CLHEP::keV 
               << " keV, Entering volume: " << postStep->GetPhysicalVolume()->GetName() << G4endl; */
    }
    
    // Record energy based on material
    if (material->GetName() == "BGO") {
        fBGOPrimarygammaEnergy[layerID] += edep;
    } else if (material->GetName() == "EJ232") {
        fEJ232PrimarygammaEnergy[layerID] += edep;
    }

    // Track energy transferred to secondaries - avoid double counting
    const G4TrackVector* secondaries = step->GetSecondary();
    if (secondaries && secondaries->size() > 0) {
        G4double newEnergyTransferred = 0.0;
        
        for (size_t i = 0; i < secondaries->size(); i++) {
            G4Track* secTrack = (*secondaries)[i];
            G4int secTrackID = secTrack->GetTrackID();
            
            // Only count secondaries we haven't seen before
            if (fPrimaryGammaTracker[trackID].trackedSecondaries.find(secTrackID) == 
                fPrimaryGammaTracker[trackID].trackedSecondaries.end()) {
                
                newEnergyTransferred += secTrack->GetKineticEnergy();
                fPrimaryGammaTracker[trackID].trackedSecondaries.insert(secTrackID);
            }
        }
        
        
        if (newEnergyTransferred > 0) {
            fPrimaryGammaTracker[trackID].transferredToSecondaries += newEnergyTransferred;
         }   
            // Debug: Print energy transferred to new secondaries
            G4cout << "Primary gamma Track ID = " << trackID
                   << ", New Energy Transferred = " << newEnergyTransferred / CLHEP::keV << " keV"
                   << ", Total Energy Transferred = " << fPrimaryGammaTracker[trackID].transferredToSecondaries / CLHEP::keV << " keV" 
                   << G4endl;
}
}
2

Here is a thread with a similar discussion.
If I am interpreting the logic correctly you could move all of this to PostUserTrackingAction, specifically this part:

// Initialize tracker for new primary gammas - only once per track
    if (fPrimaryGammaTracker.find(trackID) == fPrimaryGammaTracker.end()) {
        PrimaryGammaInfo info;
        info.initialEnergy = preStepEnergy;
        info.isFirstEntry = true;
        info.initialMaterial = material->GetName();
        info.initialLayerID = layerID;
        info.hasExited = false;
        info.depositedEnergy = 0.0;
        info.transferredToSecondaries = 0.0;
        info.exitEnergy = 0.0;
        info.trackedSecondaries = std::set<G4int>();

        fPrimaryGammaTracker[trackID] = info;

You could define these parameters in EventAction and then call it only in the PostUserTrackingAction since 1 and only 1 particle has a trackID of 1. That would simplify your code immensely and lead to a slight improvement since TrackingActions are called far less often than steps.

If you want to differentiate primary compton vs PE events than you will instead put all of that info class structure into your own UserTrackingInfo class. Specifically follow the example of Run and Event 01 and pass that info to all descendants by this cookbook. You could set a field to be “first interaction” that is either “compton” or “photoelectric” and then record all energy deposited in your detector into an ntuple that includes that field so you can filter it however you please afterward.

3

By adjusting my code, I have successfully summed up the energy transferred to secondary particles. However, the issue now is how to access this energy within the EventAction class. The terminal output confirms that the energy summation for secondary particles is correct. The “Total Process Energy” output is accurate when there is only a single process, but when multiple processes (two or three) are involved, the total energy calculation appears to be incorrect. How can I resolve this issue?

Screenshot from 2025-03-02 11-36-02
Screenshot from 2025-03-02 11-36-02916×334 66.2 KB

this is my code `// Track energy transferred to secondaries - avoid double counting

    const G4TrackVector* secondaries = step->GetSecondary();
if (secondaries && secondaries->size() > 0) {
    const G4VProcess* process = postStep->GetProcessDefinedStep();
    G4String processName = "Unknown";
    if (process) processName = process->GetProcessName();
    
    G4double processEnergy = 0.0;
    
    for (size_t i = 0; i < secondaries->size(); i++) {
	G4Track* secTrack = (*secondaries)[i];
	G4int secTrackID = secTrack->GetTrackID();
	G4double secEnergy = secTrack->GetKineticEnergy();		
	processEnergy += secEnergy;
	fPrimaryGammaTracker[trackID].trackedSecondaries.insert(secTrackID);
    }
    
    // Add to process-specific tracking
    fPrimaryGammaTracker[trackID].energyByProcess[processName] += processEnergy;
    fPrimaryGammaTracker[trackID].transferredToSecondaries += processEnergy;
    
    // Add to our process energy tracking
	    AddProcessEnergy(processName, processEnergy);
    G4cout << "Process: " << processName << ", Energy: " << processEnergy/CLHEP::keV << " keV" << G4endl;
}`  `

void ScintillatorSD::AddProcessEnergy(const G4String& processName, G4double energy)
{
fProcessEnergies[processName] += energy;
}

G4double ScintillatorSD::GetTotalPhysicalProcessEnergy() const
{
G4double total = 0.0;
for (const auto& process : fProcessEnergies) {
if (process.first != “Transportation”) {
total += process.second;
}
}
return total;
}

G4double ScintillatorSD::GetEnergyByProcess(const G4String& processName) const
{
auto it = fProcessEnergies.find(processName);
if (it != fProcessEnergies.end()) {
return it->second;
}
return 0.0;
}

void ScintillatorSD::ClearProcessEnergies()
{
fProcessEnergies.clear();
} if (bgoSD && ej232SD) {
// Get process energies from both sensitive detectors
G4double bgoPhysicalEnergy = bgoSD->GetTotalPhysicalProcessEnergy();
G4double ej232PhysicalEnergy = ej232SD->GetTotalPhysicalProcessEnergy();
G4double totalPhysicalEnergy = bgoPhysicalEnergy + ej232PhysicalEnergy;

    // Get process-specific energies
    G4double bgoComptEnergy = bgoSD->GetEnergyByProcess("compt");
    G4double bgoPhotEnergy = bgoSD->GetEnergyByProcess("phot");
    G4double ej232ComptEnergy = ej232SD->GetEnergyByProcess("compt");
    G4double ej232PhotoEnergy = ej232SD->GetEnergyByProcess("phot");
    
    G4double totalComptEnergy = bgoComptEnergy + ej232ComptEnergy;
    G4double totalPhotEnergy = bgoPhotEnergy + ej232PhotoEnergy;
    
    // Debug output
    G4cout << "Event ID: " << event->GetEventID()
           << " Total Physical Process Energy: " << totalPhysicalEnergy / keV << " keV"
           << " (Compton: " << totalComptEnergy / keV << " keV, "
           << "Photo: " << totalPhotEnergy / keV << " keV)" << G4endl;

}`

4

It is unclear what you are trying to measure/separate.
A photo peak event can consist of any number of compton scatters + a final photoeletric absorption. Two options:

  1. Create your own Tracking Info

Create a field that starts at 0 and increments every time you encounter a CS/PE process in your stepping action that is passed to all secondaries. Histogram or record the energy based on the current value of the counter. You could even just make a set inside the TrackingInfo that keeps track of what interaction of the two took place. For photopeak events it should ALWAYS end with a PE (‘CS’, ‘CS’, ‘PE’ and not ‘CS’, ‘CS’).
That is the cleanest, avoids edge cases, and avoids double counting. You could even get more sophisticated and also include the X, Y, Z coordinates of each of these sites.

  1. Just check for PE

Alternatively with your code, you could try just checking if PE occurs and summing the energy for PE. Everything else is lumped into CS. There is no physically real track that will have a compton scatter after PE absorption.