_Geant4 Version: 4-11.2
_Operating System: Windows 11
_Compiler/Version: Visual Studio 17.9
_CMake Version: 3.29.2
Hello,
I am working on a simulation in which I am propagating muons through 10m of rock (downward from the +z axis) and sampling and recording the number of secondaries produced every 4cm, as well as recording the distance each secondary travels in the x-y plane. The physics list I am using QGSP_BERT_HP.
I know that for neutrons, when they undergo inelastic scattering, Geant4 essentially ‘kills’ the primary neutron and ‘creates’ a new trackID for it. For example, if a 100 MeV neutron undergoes neutronInelastic scattering and the end product is two neutrons, one with 70 MeV and the other with 30 MeV energy, I am treating the 70 MeV neutron as the continuation of the primary 100 MeV neutron. This is done in order to avoid overcounting neutrons produced and to also make sure the creation information is passed to the correct neutron post interaction.
My question is, does Geant4 also do this for other secondaries? Would I need to handle all secondary cases of inelastic scattering, multiple scattering, abosrption and reemission, etc? The observables I am after are an averaged secondary particle production vs rock depth and averaged transverse radius vs rock depth, so I need to ensure my tracking of secondary information is correct. I have put my SteppingAction.cc code below for reference. Thank you!
#include "G4AntiProton.hh"
#include "SteppingAction.hh"
#include "RunAction.hh"
#include "G4Step.hh"
#include "G4Track.hh"
#include "G4VTouchable.hh"
#include "G4VPhysicalVolume.hh"
#include "G4LogicalVolume.hh"
#include "G4SystemOfUnits.hh"
#include "G4RunManager.hh"
#include "G4ParticleDefinition.hh"
#include "G4VProcess.hh"
#include "G4Neutron.hh"
#include "G4Proton.hh"
#include "G4Gamma.hh"
#include "G4Electron.hh"
#include "G4Positron.hh"
#include "G4PionPlus.hh"
#include "G4PionMinus.hh"
#include "G4KaonPlus.hh"
#include "G4KaonMinus.hh"
#include "G4MuonPlus.hh"
#include "G4MuonMinus.hh"
#include "G4Deuteron.hh"
#include "G4Alpha.hh"
#include "G4Triton.hh"
SteppingAction::SteppingAction()
: G4UserSteppingAction(),
fRunAction(nullptr)
{
}
SteppingAction::~SteppingAction()
{
}
// Calculate thickness based on z-position (distance from top of rock)
G4double SteppingAction::GetThickness(const G4ThreeVector& position) {
// Rock top is at z = 19m
// Thickness = distance from top
return (19.0 * m - position.z()) / m; // Convert to meters
}
G4bool SteppingAction::IsTrackedSecondary(const G4ParticleDefinition* particle) {
return (particle == G4Neutron::Definition() ||
particle == G4Proton::Definition() ||
particle == G4AntiProton::Definition() ||
particle == G4Gamma::Definition() ||
particle == G4Electron::Definition() ||
particle == G4Positron::Definition() ||
particle == G4PionPlus::Definition() ||
particle == G4PionMinus::Definition() ||
particle == G4KaonPlus::Definition() ||
particle == G4KaonMinus::Definition() ||
particle == G4MuonPlus::Definition() ||
particle == G4MuonMinus::Definition() ||
particle == G4Deuteron::Definition() ||
particle == G4Triton::Definition() ||
particle == G4Alpha::Definition());
}
G4bool SteppingAction::IsHighestEnergyNeutron(const G4Step* step) {
const G4TrackVector* secondaries = step->GetSecondary();
if (!secondaries) return true;
G4double thisEnergy = step->GetTrack()->GetKineticEnergy();
G4bool isHighest = true;
for (size_t i = 0; i < secondaries->size(); i++) {
G4Track* secondary = (*secondaries)[i];
if (secondary->GetDefinition() == G4Neutron::Definition()
) {
if (secondary->GetKineticEnergy() > thisEnergy) {
isHighest = false;
break;
}
}
}
return isHighest;
}
void SteppingAction::UserSteppingAction(const G4Step* step) {
if (!fRunAction) {
fRunAction = dynamic_cast<RunAction*>
(const_cast<G4UserRunAction*>(G4RunManager::GetRunManager()->GetUserRunAction()));
}
G4Track* track = step->GetTrack();
const G4ParticleDefinition* particle = track->GetDefinition();
// Only process tracked secondaries and ignore primary particles
if (!IsTrackedSecondary(particle) || track->GetTrackID() == 1) {
return;
}
// Check energy threshold for electrons/positrons
if ((particle == G4Electron::Definition() || particle == G4Positron::Definition())
&& track->GetKineticEnergy() <= 0.0 * MeV) {
return;
}
// Get current volume
G4LogicalVolume* volume = step->GetPreStepPoint()->GetTouchableHandle()
->GetVolume()->GetLogicalVolume();
// Check if we're in the rock volume
if (volume->GetName() == "B3") {
G4int trackID = track->GetTrackID();
G4ThreeVector preStepPos = step->GetPreStepPoint()->GetPosition();
G4ThreeVector postStepPos = step->GetPostStepPoint()->GetPosition();
// Calculate thickness from top of rock
G4double thickness = GetThickness(preStepPos);
// Handle new particle creation
if (track->GetCurrentStepNumber() == 1) {
const G4VProcess* creatorProcess = track->GetCreatorProcess();
G4String processName = creatorProcess ? creatorProcess->GetProcessName() : "";
G4bool shouldTrack = false;
if (particle == G4Neutron::Definition() && processName == "neutronInelastic") {
if (!IsHighestEnergyNeutron(step)) {
shouldTrack = true;
}
}
else if (!processName.empty()) {
shouldTrack = true;
}
if (shouldTrack) {
// Store creation position and register count
fParticleCreationPositions[trackID] = preStepPos;
fParticleLastPosition[trackID] = preStepPos;
fRunAction->AddSecondaryCount(thickness);
// Debug output
//G4cout << "New secondary: "
// << "\n Particle: " << particle->GetParticleName()
// << "\n Process: " << processName
// << "\n Position (x,y,z): (" << preStepPos.x() / m << ","
// << preStepPos.y() / m << "," << preStepPos.z() / m << ") m"
// << "\n Thickness from top: " << thickness << " m"
// << "\n Energy: " << track->GetKineticEnergy() / MeV << " MeV"
// << G4endl;
}
}
// Track end-of-life for registered particles
auto creationPosIt = fParticleCreationPositions.find(trackID);
if (creationPosIt != fParticleCreationPositions.end()) {
if (track->GetTrackStatus() == fStopAndKill ||
step->GetPostStepPoint()->GetStepStatus() == fGeomBoundary) {
// Calculate transverse distance using only x and y components
G4ThreeVector creationPosXY(creationPosIt->second.x(), creationPosIt->second.y(), 0);
G4ThreeVector finalPosXY(postStepPos.x(), postStepPos.y(), 0);
G4double transDistance = (finalPosXY - creationPosXY).mag() / mm;
G4double creationThickness = GetThickness(creationPosIt->second);
// Debug output for end-of-life
//G4cout << "Particle end-of-life:"
// << "\n Particle: " << particle->GetParticleName()
// << "\n Creation thickness: " << creationThickness << " m"
// << "\n Transverse distance: " << transDistance << " mm"
// << "\n Status: " << (track->GetTrackStatus() == fStopAndKill ? "Killed" : "Left Volume")
// << G4endl;
if (transDistance > 0.0001) {
fRunAction->AddTransverseRadius(creationThickness, transDistance);
}
fParticleCreationPositions.erase(trackID);
fParticleLastPosition.erase(trackID);
}
}
}
}