_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:
- Initial Kinetic Energy of Secondaries - This shows the energy of secondary particles at the moment of their creation.
- Escaped Secondaries - This represents the kinetic energy of secondaries that leave the detector volume.
- 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;
}
`
