Energy Tracking Discrepancy in Detector Simulation – Possible Causes

_Geant4 Version:_11.2.1
_Operating System:_Ubuntu
_Compiler/Version:_11.4.0
_CMake Version:_3.22.1
Hello experts,

I hope you’re all doing well.

I’m currently facing an issue with tracking the energy deposited by secondary particles in my detector simulation. I’ve attached three histograms for reference:

  1. Initial Kinetic Energy of Secondaries - This shows the energy of secondary particles at the moment of their creation.
  2. Escaped Secondaries - This represents the kinetic energy of secondaries that leave the detector volume.
  3. Energy Balance for Secondaries - This is the sum of the energy deposited by secondaries within the detector and the energy carried away by those that escape.

Upon comparing the first and third histograms, I’ve noticed a ~12 keV discrepancy in the mean energy. Specifically, the total energy accounted for (deposited + escaped) is lower than the initial kinetic energy, resulting in a noticeable gap.

This discrepancy appears to affect the overall energy deposition spectrum, particularly reducing the number of events with deposited energy above 400 keV. As a consequence, I’m observing a 2% discrepancy in the photoelectric interaction probability compared to expected values.

My question is:
What could be causing this energy gap?



Is there a potential flaw in the logic I’m using to track the energy deposited by secondary particles? Or could there be another explanation for this loss? This is my senstive detector class `

namespace Project
{

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

ScintillatorSD::ScintillatorSD(const G4String& name,
const G4String& hitsCollectionName,
G4int nofCells)
: G4VSensitiveDetector(name),
fNofCells(nofCells),
fBGOEnergy(nofCells, 0.0),
fEJ232Energy(nofCells, 0.0)
{
collectionName.insert(hitsCollectionName);
fBGOPrimarygammaEnergy.resize(nofCells, 0.0);
fEJ232PrimarygammaEnergy.resize(nofCells, 0.0);
fBGOSecondaryEnergy.resize(nofCells, 0.0);
fEJ232SecondaryEnergy.resize(nofCells, 0.0);

fBGOPrimarygammaCompt.resize(nofCells, 0.0);
fEJ232PrimarygammaCompt.resize(nofCells, 0.0);

fBGOSecondariesPhoto.resize(fNofCells, 0.0);
fEJ232SecondariesPhoto.resize(fNofCells, 0.0);
fBGOSecondariesComp.resize(fNofCells, 0.0);
fEJ232SecondariesComp.resize(fNofCells, 0.0);
}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
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();
}
const std::map<G4String, G4double>& ScintillatorSD::GetAllProcessEnergies() const
{
return fProcessEnergies;
}
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
G4double ScintillatorSD::GetPrimaryInitialEnergy() const {
G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
// G4cout << "GetPrimaryInitialEnergy() called at event ID: " << eventID << G4endl;

auto it = fPrimaryGammaEventData.find(eventID);
if (it != fPrimaryGammaEventData.end()) {
 //   G4cout << "Debug: Returning stored primary initial energy = " 
        //   << it->second.initialEnergy / CLHEP::MeV << " MeV" << G4endl;
    return it->second.initialEnergy;
}

// G4cout << "Debug: No data found for event ID " << eventID << “, returning 0.0” << G4endl;
return 0.0;
}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
void ScintillatorSD::Initialize(G4HCofThisEvent* hce)
{
// Create hits collection
fHitsCollection
= new ScintillatorHitsCollection(SensitiveDetectorName, collectionName[0]);

// Add this collection in hce
auto hcID
= G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
hce->AddHitsCollection( hcID, fHitsCollection );

for (G4int i=0; i<fNofCells+1; i++ ) {
fHitsCollection->insert(new ScintillatorHit());
}

// Clear energy deposition vectors
std::fill(fBGOEnergy.begin(), fBGOEnergy.end(), 0.0);
std::fill(fEJ232Energy.begin(), fEJ232Energy.end(), 0.0);
std::fill(fBGOPrimarygammaCompt.begin(), fBGOPrimarygammaCompt.end(), 0.0);
std::fill(fEJ232PrimarygammaCompt.begin(), fEJ232PrimarygammaCompt.end(), 0.0);
std::fill(fBGOSecondariesPhoto.begin(), fBGOSecondariesPhoto.end(), 0.0);
std::fill(fEJ232SecondariesPhoto.begin(), fEJ232SecondariesPhoto.end(), 0.0);
std::fill(fBGOSecondariesComp.begin(), fBGOSecondariesComp.end(), 0.0);
std::fill(fEJ232SecondariesComp.begin(), fEJ232SecondariesComp.end(), 0.0);

std::fill(fBGOPrimarygammaEnergy.begin(), fBGOPrimarygammaEnergy.end(), 0.0);
std::fill(fEJ232PrimarygammaEnergy.begin(), fEJ232PrimarygammaEnergy.end(), 0.0);
std::fill(fBGOSecondaryEnergy.begin(), fBGOSecondaryEnergy.end(), 0.0);
std::fill(fEJ232SecondaryEnergy.begin(), fEJ232SecondaryEnergy.end(), 0.0);

fPrimaryGammaTracker.clear();
ClearProcessEnergies();
fProcessedEventIDs.clear();
fEscapedPrimaryGammaEnergy.clear();
fEscapedSecondaryEnergy.clear();
fSecondaryTracker.clear();

}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

G4bool ScintillatorSD::ProcessHits(G4Step* step, G4TouchableHistory*)
{
// Get track information
G4Track* track = step->GetTrack();
G4String particleName = track->GetDefinition()->GetParticleName();
G4int parentID = track->GetParentID();
G4int trackID = track->GetTrackID();
G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
//////////THis Block Handle hit collection by Processess ///////////

if (particleName == "gamma" && parentID == 0) {

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

// Get the process that created interactions
    const G4VProcess* process = step->GetPostStepPoint()->GetProcessDefinedStep();
    if (!process) return false;
    G4String processName = step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();            
   
    auto edep = step->GetTotalEnergyDeposit();
    G4double stepLength = step->GetStepLength();  
     
// Get geometry info
    auto touchable = step->GetPreStepPoint()->GetTouchable();
    G4int layerID = touchable->GetVolume(1)->GetCopyNo();
        
// Check bounds
if (layerID < 0 || layerID >= fNofCells) {
    G4ExceptionDescription msg;
    msg << "Layer ID out of bounds: " << layerID;
    G4Exception("ScintillatorSD::ProcessHits()", 
               "OutOfBounds", FatalException, msg);
    return false;
}



 /*   // Log detailed information about gamma interaction
    G4cout << "Primary Gamma Interaction: " 
           << "Process = " << processName 
           << ", Layer = " << layerID 
           << ", Energy Deposit = " << edep / CLHEP::keV << " keV" 
           << ", Step Length = " << stepLength << G4endl; */

// Get hit accounting data for this cell
auto hit = (*fHitsCollection)[layerID];
if (!hit) {
    G4ExceptionDescription msg;
    msg << "Cannot access hit " << layerID;
    G4Exception("ScintillatorSD::ProcessHits()",
               "MyCode0004", FatalException, msg);
}

// Get hit for total accounting
auto hitTotal = (*fHitsCollection)[fHitsCollection->entries()-1];
    
// Get material and store energy
G4Material* material = step->GetPreStepPoint()->GetMaterial();
G4String materialName = material->GetName();

if (processName == "phot") {
if (materialName == "BGO") {
    fBGOEnergy[layerID] += edep;
} else if (materialName == "EJ232") {
    fEJ232Energy[layerID] += edep;
}

}

else if (processName == “compt”) {
if (materialName == “BGO”) {
fBGOPrimarygammaCompt[layerID] += edep;
} else if (materialName == “EJ232”) {
fEJ232PrimarygammaCompt[layerID] += edep;
}
}

// Update the hit
hit->Add(edep, stepLength);
hitTotal->Add(edep, stepLength);

// }
}

const G4VProcess* creatorProcess = track->GetCreatorProcess();
if (parentID == 1 && creatorProcess) {    
// Get the process that created this particle

G4String procName = creatorProcess->GetProcessName();

        auto edep = step->GetTotalEnergyDeposit();
        if (edep == 0.) return false;
        
        auto touchable = step->GetPreStepPoint()->GetTouchable();
        G4int layerID = touchable->GetVolume(1)->GetCopyNo();
        if (layerID < 0 || layerID >= fNofCells) return false;
        
        // Store electron energy based on material
        G4Material* material = step->GetPreStepPoint()->GetMaterial();
    if (procName == "phot") {
    if (material->GetName() == "BGO") fBGOSecondariesPhoto[layerID] += edep;
    else if (material->GetName() == "EJ232") fEJ232SecondariesPhoto[layerID] += edep;
} 
else if (procName == "compt") {
    if (material->GetName() == "BGO") fBGOSecondariesComp[layerID] += edep;
    else if (material->GetName() == "EJ232") fEJ232SecondariesComp[layerID] += edep;
} 

}

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…
////// This loop will handle all primaries and secondaries///////

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

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

// 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 preStepEnergy  = step->GetPreStepPoint()->GetTotalEnergy();
    G4double postStepEnergy = step->GetPostStepPoint()->GetTotalEnergy();
G4double stepLength = step->GetStepLength();

G4ThreeVector prestep_momentum = step->GetPreStepPoint()->GetMomentumDirection();
G4ThreeVector prestep_position = step->GetPreStepPoint()->GetPosition();

// 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.hasEnteredDetector = false;
    info.hasExited = false;
  //  info.initialEnergy = track->GetVertexKineticEnergy();
    info.depositedEnergy = 0.0;
    info.transferredToSecondaries = 0.0;
    info.exitEnergy = 0.0;
    info.trackedSecondaries = std::set<G4int>();
info.enteredBGO = false;  
	info.enteredEJ232 = false;  
info.lastMaterial = "";

info.initialMomentum = prestep_momentum;  
info.initialPosition = prestep_position;

    fPrimaryGammaTracker[trackID] = info;
}

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

if (fPrimaryGammaTracker[trackID].isFirstEntry) {
    fPrimaryGammaTracker[trackID].isFirstEntry = false; // Mark as processed
}

// Track which detector material this particle has entered
if (material->GetName() == "BGO") {
    fPrimaryGammaTracker[trackID].hasEnteredDetector = true;
    fPrimaryGammaTracker[trackID].enteredBGO = true;
    fPrimaryGammaTracker[trackID].lastMaterial = "BGO";
    
    if (fPrimaryGammaTracker[trackID].stepLengths.empty()) {
	fPrimaryGammaTracker[trackID].stepLengths.clear();
    }
} else if (material->GetName() == "EJ232") {
    fPrimaryGammaTracker[trackID].hasEnteredDetector = true;
    fPrimaryGammaTracker[trackID].enteredEJ232 = true;
    fPrimaryGammaTracker[trackID].lastMaterial = "EJ232";
    
    if (fPrimaryGammaTracker[trackID].stepLengths.empty()) {
	fPrimaryGammaTracker[trackID].stepLengths.clear();
    }
}

// Only store step lengths for particles that have entered the detector
if (fPrimaryGammaTracker[trackID].hasEnteredDetector) {
    fPrimaryGammaTracker[trackID].stepLengths.push_back(step->GetStepLength());
    fStoredStepLengths[trackID].push_back(step->GetStepLength());
}

// Check if particle is truly exiting the ENTIRE detector system
if ((preStep->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName() == "BGO" || 
     preStep->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName() == "EJ232")) {
    
    // Now check if the next step is outside ALL detector volumes
    bool isExitingAllDetectors = false;
    
    // Case 1: Exiting to the world volume or another non-detector volume
    if (postStep->GetPhysicalVolume() == nullptr || 
	(postStep->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName() != "BGO" &&
	 postStep->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName() != "EJ232")) {
	
	isExitingAllDetectors = true;
    }
    
    // If we're truly exiting all detector volumes
    if (isExitingAllDetectors) {
	// Only update exit info if not already marked AND has entered detector
	if (!fPrimaryGammaTracker[trackID].hasExited && 
	    fPrimaryGammaTracker[trackID].hasEnteredDetector && 
	    preStepEnergy > 0) {  // Only count if it still has energy
	    
	    fPrimaryGammaTracker[trackID].exitEnergy = preStepEnergy;
	    fPrimaryGammaTracker[trackID].hasExited = true;
	    fEscapedPrimaryGammaEnergy[trackID] = preStepEnergy;
	    
	    G4double totalPathLength = 0.0;
	    std::vector<G4double>& steps = fPrimaryGammaTracker[trackID].stepLengths;
	    
	    // Calculate total path length
	    for (G4double step : steps) {
	        totalPathLength += step;
	    }
	    
	    // Optionally print details for debugging
	  /*  G4cout << "Particle truly escaped entire detector: TrackID = " << trackID
	           << ", Exit Energy = " << preStepEnergy / CLHEP::keV << " keV"
	           << ", Path length = " << totalPathLength / CLHEP::mm << " mm" << 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) {
   // G4ThreeVector prestep_Secmomentum = step->GetPreStepPoint()->GetMomentumDirection();
    const G4VProcess* process = postStep->GetProcessDefinedStep();
    G4String processName = "Unknown";
    if (process) processName = process->GetProcessName();
    
    // Skip transportation process
    if (processName == "Transportation") return true;
    
    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();
	//G4double secStepLength = secTrack->GetStepLength();
	 // G4cout << "Secondaries Step Length: " << secStepLength << G4endl;
	  
	// Store the secondary momentum
    	G4ThreeVector secMomentum = secTrack->GetMomentumDirection();
    	fPrimaryGammaTracker[trackID].secondaryMomenta[secTrackID] = secMomentum; 
    	// Store the secondary position (using creation position from step)
       	G4ThreeVector secPosition = secTrack->GetPosition();
        	fPrimaryGammaTracker[trackID].secondaryPositions[secTrackID] = secPosition;
    	//G4cout << "Secondaries Position: " << secPosition << G4endl;
    	     
	processEnergy += secEnergy;
	fPrimaryGammaTracker[trackID].trackedSecondaries.insert(secTrackID);
	/*G4cout << "Process: " << processName 
       << ", Track ID: " << secTrackID << "Storing initial momentum for Secondaries" << secMomentum << ", Kinetic Energy: " << secTrack->GetKineticEnergy()/CLHEP::keV << " keV"
       << G4endl; */
    }
    
   // G4cout << "Process: " << processName << ", Energy: " << processEnergy/CLHEP::keV << " keV" << G4endl;
   // G4cout << "Storing initial momentum for Secondaries" << prestep_Secmomentum << G4endl;
    // Add to process-specific tracking
    fPrimaryGammaTracker[trackID].energyByProcess[processName] += processEnergy;
    fPrimaryGammaTracker[trackID].transferredToSecondaries += processEnergy;
    
    static std::map<G4int, std::set<G4String>> processedProcesses;
    G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
    
    // If this is the first time seeing this process in this event
    if (processedProcesses[eventID].find(processName) == processedProcesses[eventID].end()) {
	AddProcessEnergy(processName, processEnergy);
	processedProcesses[eventID].insert(processName);
        //G4cout << "Process: " << processName << ", Energy: " << processEnergy/CLHEP::keV << " keV" << G4endl;
    }
}
   }

//…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…oooOO0OOooo…

if (parentID > 0) {
G4int primaryParentID = track->GetParentID();
auto edep = step->GetTotalEnergyDeposit();
if (edep == 0.) return false;

if (material->GetName() == "BGO") {
    fBGOSecondaryEnergy[layerID] += edep;
} else if (material->GetName() == "EJ232") {
    fEJ232SecondaryEnergy[layerID] += edep;
}

// Check if this is a secondary we want to track (from our primary gammas)
if (fPrimaryGammaTracker.find(primaryParentID) != fPrimaryGammaTracker.end()) {
    
    // Initialize tracker for new secondaries - only once per track
    if (fSecondaryTracker.find(trackID) == fSecondaryTracker.end()) {
        SecondaryParticleInfo info;
        info.initialEnergy = track->GetKineticEnergy();
        info.particleType = track->GetDefinition()->GetParticleName();
        info.parentID = primaryParentID;
        info.creationProcess = track->GetCreatorProcess() ? 
                              track->GetCreatorProcess()->GetProcessName() : "Unknown";
        info.initialPosition = track->GetPosition();
        info.initialMomentum = track->GetMomentumDirection();
        info.hasExited = false;
        info.hasEnteredDetector = false;
        info.exitEnergy = 0.0;
        info.depositedEnergy = 0.0;
        info.totalPathLength = 0.0;
        info.enteredBGO = false;
        info.enteredEJ232 = false;
        fSecondaryTracker[trackID] = info;
    }
    
    fSecondaryTracker[trackID].depositedEnergy += edep;
    
    // Track which detector materials this particle has entered
    if (material->GetName() == "BGO") {
        fSecondaryTracker[trackID].enteredBGO = true;
        fSecondaryTracker[trackID].hasEnteredDetector = true;
    } else if (material->GetName() == "EJ232") {
        fSecondaryTracker[trackID].enteredEJ232 = true;
        fSecondaryTracker[trackID].hasEnteredDetector = true;
    }
    
    // Only store step lengths for particles in the detector
    if (material->GetName() == "BGO" || material->GetName() == "EJ232") {
        fSecondaryTracker[trackID].stepLengths.push_back(step->GetStepLength());
    }
    
    // Check if the secondary is truly exiting the entire detector
    if (fSecondaryTracker[trackID].hasEnteredDetector && 
        (preStep->GetPhysicalVolume()->GetName() == "BGO" || 
         preStep->GetPhysicalVolume()->GetName() == "EJ232") &&
        (postStep->GetPhysicalVolume() == nullptr || 
         (postStep->GetPhysicalVolume() && 
          postStep->GetPhysicalVolume()->GetName() != "BGO" &&
          postStep->GetPhysicalVolume()->GetName() != "EJ232"))) {
        
        // Only update exit info if not already marked
        if (!fSecondaryTracker[trackID].hasExited) {
            // Only consider as "exited" if it still has energy
            G4double exitEnergy = preStep->GetKineticEnergy();
            bool createdInDetector = (track->GetCreatorProcess() && 
                                     (material->GetName() == "BGO" || material->GetName() == "EJ232"));
            bool passedThroughBoth = (fSecondaryTracker[trackID].enteredBGO && fSecondaryTracker[trackID].enteredEJ232);
            
            // Only mark as truly exited if it's exiting the LAST detector material it was in
            if (exitEnergy > 0 && (passedThroughBoth || createdInDetector)) {
                fSecondaryTracker[trackID].exitEnergy = exitEnergy;
                fSecondaryTracker[trackID].hasExited = true;
                
                // Only store the exit energy ONCE per track
                if (fEscapedSecondaryEnergy.find(trackID) == fEscapedSecondaryEnergy.end()) {
                    fEscapedSecondaryEnergy[trackID] = exitEnergy;
                    
                    // Calculate total path length
                    G4double totalPathLength = 0.0;
                    for (G4double step : fSecondaryTracker[trackID].stepLengths) {
                        totalPathLength += step;
                    }
                    fSecondaryTracker[trackID].totalPathLength = totalPathLength;
                    
                 /*   // Debugging information
                    G4cout << "Secondary Escaped: TrackID = " << trackID
                           << ", Type = " << fSecondaryTracker[trackID].particleType
                           << ", Exit Energy = " << exitEnergy / CLHEP::keV 
                           << " keV, Parent = " << primaryParentID << G4endl; */
                }
            }
        }
    }
}

}

return true;

}

`

here a macro for example TestEm1 (in G4 11.3), an extract of its printout, and the plots.
What is your material, its thickness, and the energy of the gamma beam ?

rehman.mac.txt (561 Bytes)
rehman.out.txt (667 Bytes)

If you include primary particle in your plots, you must get exactly the energy beam.
This will be a test of your code.

Hi Maire,
I am using BGO mit den Dimensionen 3 × 3 × 15 mm³. beam hit 3 × 3 mm² face. Energy of gamma is 511 keV. I have the Plot for primaries but i want to track seprately the secondaries that leave my detector. From verbosity 2 output i can see some secondaries leaving my detector without full energy deposition my code also tracking energy deposition from these secondaries but it fails when trasportation take place and particle enters into world. This is the output


  • G4Track Information: Particle = e-, Track ID = 2, Parent ID = 1

Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName
0 0.389 -1.39 -0.597 0.42 0 0 0 BGO initStep
1 0.384 -1.4 -0.58 0.402 0.0185 0.0215 0.0215 BGO msc
Secondary created in detector: TrackID=2 Material=BGO
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 420.4659999999999 keV
2 0.378 -1.42 -0.572 0.385 0.0172 0.0201 0.0416 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 401.9461509124998 keV
3 0.374 -1.43 -0.58 0.368 0.0167 0.0188 0.0604 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 384.7314472108905 keV
4 0.374 -1.43 -0.595 0.355 0.0135 0.0176 0.078 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 368.021666810087 keV
5 0.373 -1.44 -0.603 0.34 0.014 0.0167 0.0947 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 354.5010939621002 keV
6 0.368 -1.46 -0.601 0.331 0.00901 0.0157 0.11 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 340.4695952756444 keV
7 0.362 -1.47 -0.6 0.316 0.0157 0.0151 0.125 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 331.4622495403154 keV
8 0.365 -1.48 -0.601 0.308 0.00779 0.014 0.139 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 315.728985038857 keV
9 0.363 -1.49 -0.609 0.299 0.00887 0.0135 0.153 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 307.9409117250323 keV
10 0.366 -1.5 -0.611 0.292 0.00665 0.0106 0.164 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 299.0744522953214 keV
11 0.367 -1.5 -0.61 0.291 0.000934 0.00141 0.165 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 292.4269804945865 keV
12 0.367 -1.5 -0.609 0.29 0.00121 0.000616 0.166 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 291.4927479114677 keV
13 0.367 -1.5 -0.609 0.29 2.21e-05 8.22e-05 0.166 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 290.2869441696222 keV
14 0.367 -1.5 -0.609 0.29 0 5.69e-05 0.166 BGO msc
15 0.367 -1.5 -0.609 0.29 5.01e-05 0.000158 0.166 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 290.2648497109943 keV
16 0.367 -1.5 -0.609 0.29 1.29e-05 0.000126 0.166 BGO msc
Secondary Step - TrackID: 2 PV: BGO PreMat: BGO → BGO PostMat: BGO Energy: 290.2147147494251 keV
17 0.367 -1.5 -0.609 0.29 0 7.64e-05 0.166 World Transportation
18 0.386 -1.61 -0.448 0.29 7.28e-05 0.194 0.36 World msc
19 0.397 -1.66 -0.357 0.29 0 0.108 0.468 World msc
20 0.569 -2.57 1 0.29 0.000101 1.64 2.11 World msc
21 0.614 -2.73 1.27 0.29 0 0.321 2.43 World msc
22 0.714 -3.03 1.76 0.29 1.95e-05 0.581 3.01 World msc
23 0.753 -3.48 2.72 0.29 0.000185 1.06 4.07 World msc
24 0.761 -3.54 2.85 0.289 0.000755 0.146 4.21 World msc
25 0.815 -4.37 4.71 0.289 3.45e-05 2.04 6.25 World msc
26 0.737 -5.6 7.3 0.288 0.000745 2.87 9.12 World msc
27 0.533 -7.24 11.1 0.288 0.000469 4.1 13.2 World msc
28 0.589 -9.7 16.2 0.287 0.00105 5.74 19 World msc
29 1.6 -12.3 24.5 0.286 0.00117 8.72 27.7 World msc
30 1.21 -16.7 40.9 0.281 0.00424 17 44.7 World msc
31 -5.38 -27.1 71.7 0.276 0.0058 33.4 78 World msc
32 -9.22 -31 80.8 0.272 0.00141 10.6 88.7 World eIoni
:----- List of 2ndaries - #SpawnInStep= 1(Rest= 0,Along= 0,Post= 1), #SpawnTotal= 1 ---------------
: -9.22 -31 80.8 0.00213 e-
:----------------------------------------------------------------- EndOf2ndaries Info ---------------
33 -35 -66.1 135 0.255 0.0157 69.6 158 World eIoni
:----- List of 2ndaries - #SpawnInStep= 1(Rest= 0,Along= 0,Post= 1), #SpawnTotal= 2 ---------------
: -35 -66.1 135 0.00131 e-
:----------------------------------------------------------------- EndOf2ndaries Info ---------------
34 -42.2 -75.4 149 0.25 0.00366 18.3 177 World eIoni
:----- List of 2ndaries - #SpawnInStep= 1(Rest= 0,Along= 0,Post= 1), #SpawnTotal= 3 ---------------
: -42.2 -75.4 149 0.00133 e-
:----------------------------------------------------------------- EndOf2ndaries Info ---------------
35 -45 -82.2 159 0.245 0.00336 12.5 189 World eIoni
:----- List of 2ndaries - #SpawnInStep= 1(Rest= 0,Along= 0,Post= 1), #SpawnTotal= 4 ---------------
: -45 -82.2 159 0.00144 e-
:----------------------------------------------------------------- EndOf2ndaries Info ---------------
36 -36.7 -101 188 0.18 0.0107 37 226 World eIoni
:----- List of 2ndaries - #SpawnInStep= 1(Rest= 0,Along= 0,Post= 1), #SpawnTotal= 5 ---------------
: -36.7 -101 188 0.0543 e-
:----------------------------------------------------------------- EndOf2ndaries Info ---------------
37 -30.2 -118 200 0.173 0.00755 21.7 248 World msc
38 -27.1 -126 205 0.171 0.00195 10.1 258 World msc
39 -25.6 -130 207 0.17 0.000982 4.79 263 World msc
40 -24.8 -133 208 0.169 0.00125 3.01 266 World msc
41 -24.4 -135 209 0.168 0.00012 1.87 268 World msc
42 -24.3 -135 209 0.168 8.79e-05 0.434 268 World msc
43 -24.3 -135 209 0.168 7.15e-06 0.0307 268 World msc
44 -24.2 -135 209 0.168 0 0.319 268 World msc
45 -24.1 -136 210 0.168 0.00011 0.557 269 World msc
46 -23.8 -137 210 0.168 0.000289 1.1 270 OutOfWorld Transportation "
This is the code Block that i am using to track secondaries escape `// For secondary particles (electrons/positrons)

if (parentID > 0) {

auto edep = step->GetTotalEnergyDeposit();
// Record energy based on material
if (material->GetName() == "BGO") {
	fBGOSecondaryEnergy[layerID] += edep;
    } else if (material->GetName() == "EJ232") {
	fEJ232SecondaryEnergy[layerID] += edep;
    }
    
    // Initialize tracker if new secondary
    if (fSecondaryTracker.find(trackID) == fSecondaryTracker.end()) {
	SecondaryParticleInfo info;
	info.initialEnergy = track->GetKineticEnergy();
	info.particleType = particleName;
	info.parentID = parentID;
	info.creationProcess = track->GetCreatorProcess() ? 
	                      track->GetCreatorProcess()->GetProcessName() : "Unknown";
	info.initialPosition = track->GetPosition();
	info.initialMomentum = track->GetMomentumDirection();
	info.hasExited = false;
	info.exitEnergy = 0.0;
	info.depositedEnergy = 0.0;
	info.totalPathLength = 0.0;
	
	// Set hasEnteredDetector based on creation volume
	const G4LogicalVolume* creatorVolume = track->GetLogicalVolumeAtVertex();
	if (creatorVolume) {
	    G4String creatorMaterial = creatorVolume->GetMaterial()->GetName();
	    info.hasEnteredDetector = (creatorMaterial == "BGO" || creatorMaterial == "EJ232");
	    if (info.hasEnteredDetector) {
	        G4cout << "Secondary created in detector: TrackID=" << trackID 
	               << " Material=" << creatorMaterial << G4endl;
	    }
	}
	
	fSecondaryTracker[trackID] = info;
    }

    // Get current step information
    G4StepPoint* preStep = step->GetPreStepPoint();
    G4StepPoint* postStep = step->GetPostStepPoint();
    
    // Get material information
    G4String preMat = preStep->GetMaterial()->GetName();
    G4String postMat = postStep->GetPhysicalVolume() ? 
	              postStep->GetMaterial()->GetName() : "World";

    // Debug output
    G4cout << "Secondary Step - TrackID: " << trackID 
	   << " PV: " << (preStep->GetPhysicalVolume() ? preStep->GetPhysicalVolume()->GetName() : "nullptr")
	   << " PreMat: " << preMat << " -> "
	   << (postStep->GetPhysicalVolume() ? postStep->GetPhysicalVolume()->GetName() : "World")
	   << " PostMat: " << postMat
	   << " Energy: " << preStep->GetKineticEnergy()/CLHEP::keV << " keV"
	   << G4endl;
	   
    // Store step length if in detector material
    if (preMat == "BGO" || preMat == "EJ232") {
	fSecondaryTracker[trackID].stepLengths.push_back(step->GetStepLength());
    }

    // Escape detection - check transition from detector to World
    if (!fSecondaryTracker[trackID].hasExited && 
	fSecondaryTracker[trackID].hasEnteredDetector &&
	(preMat == "BGO" || preMat == "EJ232") && 
	postMat == "World") {
	
	G4double exitEnergy = preStep->GetKineticEnergy();
	fSecondaryTracker[trackID].exitEnergy = exitEnergy;
	fSecondaryTracker[trackID].hasExited = true;
	fEscapedSecondaryEnergy[trackID] = exitEnergy;
	
	// Calculate total path length in detector
	G4double totalPathLength = 0.0;
	for (G4double step : fSecondaryTracker[trackID].stepLengths) {
	    totalPathLength += step;
	}
	fSecondaryTracker[trackID].totalPathLength = totalPathLength;
	
	G4cout << ">>>> ESCAPE DETECTED: TrackID=" << trackID
	       << " Type=" << particleName
	       << " ExitEnergy=" << exitEnergy/CLHEP::keV << " keV"
	       << " PathLength=" << totalPathLength/CLHEP::mm << " mm"
	       << G4endl;
    }
}`

rehman.mac.txt (530 Bytes)
rehman.out.txt (872 Bytes)

About SteppingVerbose, may I suggest to put in your main.cc :
G4SteppingVerbose::UseBestUnit(precision)

See exampleB1.cc, lines 57-59

I don’t get the logic of how you are getting the energy leakage?
I am tracking the escaped primaries (511 keV gamma), and I am getting the correct values for these escaped particles. The problem is that when I applied similar logic to secondaries, I got few secondaries that escaped from my detector, which means I am not getting the actual number of escaped secondaries. I am interested in escaped tracking of both secondaries and primaries. I am attaching the plot of escaped primaries and secondaries. My main concern is the number of escaped.



Problem is in this logic its not called for all events where secondaries leave the detector.

    // Check if the secondary is truly exiting the entire detector
    if (fSecondaryTracker[trackID].hasEnteredDetector && 
	(preStep->GetPhysicalVolume()->GetName() == "BGO" || 
	 preStep->GetPhysicalVolume()->GetName() == "EJ232") &&
	(postStep->GetPhysicalVolume() == nullptr || 
	 (postStep->GetPhysicalVolume() && 
	  postStep->GetPhysicalVolume()->GetName() != "BGO" &&
	  postStep->GetPhysicalVolume()->GetName() != "EJ232"))) {
	
	// Only update exit info if not already marked
	if (!fSecondaryTracker[trackID].hasExited) {
	    G4double exitEnergy = preStep->GetKineticEnergy();
	  
	    if (exitEnergy > 0 ) {
	        fSecondaryTracker[trackID].exitEnergy = exitEnergy;
	        fSecondaryTracker[trackID].hasExited = true;
	        
	        // Only store the exit energy ONCE per track
	        if (fEscapedSecondaryEnergy.find(trackID) == fEscapedSecondaryEnergy.end()) {
	            fEscapedSecondaryEnergy[trackID] = exitEnergy;

In example TestEm1, the evaluation of leakage is rather recent; it is present in G4 11.3 (December 2024), but not in 11.2 (December 2023).
Today, in my working copy, I have plotted separately leakage from primary and secondaries. (see TrackingAction.cc.txt)
Over 3935 events which have leakage, only 74 of them have leakage from secondaries. As you have already noticed ! …

rehman.mac.txt (592 Bytes)
TrackingAction.cc.txt (2.4 KB)