_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 trying to record neutron count and information for when underground muon-induced neutrons pass from a rock boundary into a “lab” space. The issue is that when a neutron enters the lab space, the boundary crossing doesn’t get triggered, despite clear visual proof a neutron entered the lab from the rock. I have tried writing out some debugging to show that the neutrons are indeed first created in the rock and also tagged in the rock volume, and the output is good. I have put below my stepping code. I have a 13.6*m radius world volume which I am propagating muon downward from. My lab space is 3m x 3m x 3.5m, which 5m of rock on the sides and top. The lab is placed such that the bottom is on the equator of the world volume. Any help would be greatly appreciated.
#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 "CLHEP/Units/PhysicalConstants.h"
#include <map>
#include <fstream>
#include <iomanip>
SteppingAction::SteppingAction()
: G4UserSteppingAction(),
fCurrentTrackID(-1) {
//G4cout << "Opening data file..." << G4endl;
fEscapedNeutronData.open("C:/Users/User/GEANT4/share/Geant4/examples/Geant4 mu-induced neutrons/secondaries/escaped_neutron_data.csv", std::ios::out | std::ios::app);
if (!fEscapedNeutronData.is_open()) {
G4cerr << "ERROR: Could not open neutron data file!" << G4endl;
return;
}
if (fEscapedNeutronData.tellp() == 0) {
fEscapedNeutronData << "TrackID,EventID,CreationEnergy (MeV),ExitEnergy (MeV),"
<< "ExitX (mm),ExitY (mm),ExitZ (mm),"
<< "CreationX (mm),CreationY (mm),CreationZ (mm),"
<< "Depth (mm),"
<< "TrackLength (mm),MuonTheta (deg),MuonPhi (deg),MuonInitialEnergy (GeV),creationProcess" << std::endl;
}
}
SteppingAction::~SteppingAction() {
if (fEscapedNeutronData.is_open()) {
fEscapedNeutronData.close();
}
}
void SteppingAction::UserSteppingAction(const G4Step* step) {
G4Track* track = step->GetTrack();
if (!track) return;
G4int trackID = track->GetTrackID();
G4int eventID = G4RunManager::GetRunManager()->GetCurrentEvent()->GetEventID();
// SECTION 1: Primary muon tracking
if (trackID == 1) {
const G4ParticleDefinition* particle = track->GetParticleDefinition();
if (particle && particle->GetParticleName() == "mu-") {
if (track->GetCurrentStepNumber() == 1) {
G4ThreeVector momentum = track->GetVertexMomentumDirection();
G4double cosTheta = std::abs(momentum.z());
G4double phi = std::atan2(momentum.y(), momentum.x());
if (phi < 0) phi += 2 * CLHEP::pi;
MuonAngles angles;
angles.theta = std::acos(cosTheta) * 180.0 / CLHEP::pi;
angles.phi = phi * 180.0 / CLHEP::pi;
angles.initialEnergy = track->GetKineticEnergy();
fMuonAngles[eventID] = angles;
//const G4VProcess* currentProcess = track->GetStep()->GetPostStepPoint()->GetProcessDefinedStep();
//G4cout << "Muon Capture Details:" << G4endl;
//G4cout << " Process Name: "
// << (currentProcess ? currentProcess->GetProcessName() : "Unknown") << G4endl;
//G4cout << " Track ID: " << track->GetTrackID() << G4endl;
//G4cout << " Capture Energy: " << track->GetKineticEnergy() / MeV << " MeV" << G4endl;
//G4cout << "Stored muon info for event " << eventID
// << " theta=" << angles.theta
// << " phi=" << angles.phi
// << " E=" << angles.initialEnergy / GeV << " GeV" << G4endl;
}
}
return;
}
// SECTION 2: Neutron tracking
if (track->GetDefinition()->GetParticleName() != "neutron") return;
// Get volumes with safety checks
G4VPhysicalVolume* preVol = step->GetPreStepPoint()->GetPhysicalVolume();
G4VPhysicalVolume* postVol = step->GetPostStepPoint()->GetPhysicalVolume();
if (!preVol || !postVol) return;
G4LogicalVolume* preVolume = preVol->GetLogicalVolume();
G4LogicalVolume* postVolume = postVol->GetLogicalVolume();
if (!preVolume || !postVolume) return;
// Record new neutrons with boundary conditions
if (track->GetCurrentStepNumber() == 1 && preVolume->GetName() == "rectShell") {
G4ThreeVector pos = track->GetPosition();
//G4double radius = std::sqrt(pos.x() * pos.x() + pos.y() * pos.y() + pos.z() * pos.z());
if (((pos.x() > 1.5 * m || pos.x() < -1.5 * m) &&
(-6.5 * m < pos.y() && pos.y() < 6.5 * m) &&
(pos.z() >= 0 * m)) ||
((pos.y() > 1.5 * m || pos.y() < -1.5 * m) &&
(-6.5 * m < pos.x() && pos.x() < 6.5 * m) &&
(pos.z() >= 0 * m))) {
NeutronTrackInfo info;
info.initialPos = pos;
info.creationEnergy = track->GetVertexKineticEnergy();
fNeutronTracks[trackID] = info;
}
return;
}
/*if (track->GetDefinition()->GetParticleName() == "neutron") {
G4cout << "Neutron at position: "
<< track->GetPosition().x() / m << ", "
<< track->GetPosition().y() / m << ", "
<< track->GetPosition().z() / m
<< " PreVolume: " << preVolume->GetName()
<< " PostVolume: " << postVolume->GetName()
<< G4endl;
}*/
// Process existing neutron tracks
auto trackIter = fNeutronTracks.find(trackID);
if (trackIter == fNeutronTracks.end()) return; // Only process tracked neutrons
// Check for volume boundary crossing
if (preVolume->GetName() == "rectShell" && postVolume->GetName() != "rectShell") {
G4ThreeVector prePos = step->GetPreStepPoint()->GetPosition();
G4ThreeVector postPos = step->GetPostStepPoint()->GetPosition();
//G4double preRadius = std::sqrt(prePos.x() * prePos.x() + prePos.y() * prePos.y() + prePos.z() * prePos.z());
//G4double postRadius = std::sqrt(postPos.x() * postPos.x() + postPos.y() * postPos.y() + postPos.z() * postPos.z());
// Inner boundary crossing (15m)
if ((((prePos.x() >= 1.5 * m || prePos.x() <= -1.5 * m) &&
(-6.5 * m <= prePos.y() && prePos.y() <= 6.5 * m) &&
(prePos.z() >= 0 * m)) ||
((prePos.y() >= 1.5 * m || prePos.y() <= -1.5 * m) &&
(-6.5 * m <= prePos.x() && prePos.x() <= 6.5 * m) &&
(prePos.z() >= 0 * m))) &&
(-1.5 * m < postPos.x() && postPos.x() < 1.5 * m) &&
(-1.5 * m < postPos.y() && postPos.y() < 1.5 * m) &&
(0 * m <= postPos.z() && postPos.z() < 3.5 * m)) {
//G4cout << "Neutron crossing inner boundary at pre-radius: " << preRadius / m
// << " m, post-radius: " << postRadius / m << " m" << G4endl;
// Calculate exact crossing point
//G4double fraction = (15.0 * m - preRadius) / (postRadius - preRadius);
//G4ThreeVector crossingPoint = prePos + fraction * (postPos - prePos);
G4ThreeVector crossingPoint = postPos;
G4double preStepEnergy = step->GetPreStepPoint()->GetKineticEnergy();
G4double postStepEnergy = step->GetPostStepPoint()->GetKineticEnergy();
//G4double crossingEnergy = preStepEnergy + fraction * (postStepEnergy - preStepEnergy);
G4double exitEnergy = step->GetPostStepPoint()->GetKineticEnergy();
//G4cout << "Neutron created with energy " << trackIter->second.creationEnergy / MeV << " MeV "
// << "and exiting with energy " << crossingEnergy / MeV << G4endl;
const G4VProcess* creatorProcess = track->GetCreatorProcess();
//G4cout << "Neutron Creation Details:" << G4endl;
//G4cout << " TrackID: " << track->GetTrackID() << G4endl;
//G4cout << " Creator Process: "
// << (creatorProcess ? creatorProcess->GetProcessName() : "Unknown") << G4endl;
//G4cout << " Parent Track ID: " << track->GetParentID() << G4endl;
//G4cout << " Vertex Energy: " << track->GetVertexKineticEnergy() / MeV << " MeV" << G4endl;
auto trackIter = fNeutronTracks.find(trackID);
if (trackIter == fNeutronTracks.end()) {
//G4cout << "No creation info for neutron " << trackID << G4endl;
return;
}
try {
const G4ThreeVector& initialPos = trackIter->second.initialPos;
G4double trackLength = std::sqrt(
std::pow(crossingPoint.x() - initialPos.x(), 2) +
std::pow(crossingPoint.y() - initialPos.y(), 2) +
std::pow(crossingPoint.z() - initialPos.z(), 2));
// Validate track length
//if (trackLength > 15.0 * m) {
// G4cout << "Warning: Invalid track length " << trackLength / m
// << " m for neutron " << trackID
// << " from position " << initialPos / m << " to " << crossingPoint / m << G4endl;
//return;
//}
//G4double depth = std::sqrt(
// std::pow(initialPos.x(), 2) +
// std::pow(initialPos.y(), 2) +
// std::pow(initialPos.z(), 2)) - 15.0 * m;
// distance from lab calculations
// for z <= 3.5m …. Below the top of the lab
// 1D distance from x faces
if ((1.5 * m <= initialPos.x() || initialPos.x() <= -1.5 * m) &&
(-1.5 * m <= initialPos.y() && initialPos.y() <= 1.5 * m) &&
(initialPos.z() <= 3.5 * m)) {
depth = abs(initialPos.x() - 1.5);
}
// 1D distance from y faces
else if ((1.5 * m <= initialPos.y() || initialPos.y() <= -1.5 * m) &&
(-1.5 * m <= initialPos.x() && initialPos.x() <= 1.5 * m) &&
(initialPos.z() <= 3.5 * m)) {
depth = abs(initialPos.y() - 1.5);
}
// 2D distance from x > 1.5 and y > 1.5
else if ((1.5 * m < initialPos.x()) && (1.5 * m < initialPos.y()) &&
(initialPos.z() <= 3.5 * m)) {
depth = std::sqrt(pow(initialPos.x() - 1.5, 2) + pow(initialPos.y() - 1.5, 2));
}
// 2D distance from x > 1.5 and y < -1.5
else if ((1.5 * m < initialPos.x()) && (-1.5 * m > initialPos.y()) &&
(initialPos.z() <= 3.5 * m)) {
depth = std::sqrt(pow(initialPos.x() - 1.5, 2) + pow(initialPos.y() - 1.5, 2));
}
// 2D distance from x < -1.5 and y < -1.5
else if ((-1.5 * m > initialPos.x()) && (-1.5 * m > initialPos.y()) &&
(initialPos.z() <= 3.5 * m)) {
depth = std::sqrt(pow(initialPos.x() - 1.5, 2) + pow(initialPos.y() - 1.5, 2));
}
// 2D distance from x < -1.5 and y > 1.5
else if ((-1.5 * m > initialPos.x()) && (1.5 * m < initialPos.y()) &&
(initialPos.z() <= 3.5 * m)) {
depth = std::sqrt(pow(initialPos.x() - 1.5, 2) + pow(initialPos.y() - 1.5, 2));
}
// for z > 3.5m ... above the top of the lab
// 2D distance from top of x faces
else if ((1.5 * m <= initialPos.x() || initialPos.x() <= -1.5 * m) &&
(-1.5 * m <= initialPos.y() && initialPos.y() <= 1.5 * m) &&
(initialPos.z() > 3.5 * m)) {
depth = std::sqrt(pow(initialPos.x() - 1.5, 2) + pow(initialPos.z() - 3.5, 2));
}
// 2D distance from top of y faces
else if ((1.5 * m <= initialPos.y() || initialPos.y() <= -1.5 * m) &&
(-1.5 * m <= initialPos.x() && initialPos.x() <= 1.5 * m) &&
(initialPos.z() > 3.5 * m)) {
depth = std::sqrt(pow(initialPos.y() - 1.5, 2) + pow(initialPos.z() - 3.5, 2));
}
// 3D distance from x > 1.5 and y > 1.5
else if ((1.5 * m < initialPos.x()) && (1.5 * m < initialPos.y()) &&
(initialPos.z() > 3.5 * m)) {
depth = std::sqrt(pow(initialPos.x() - 1.5, 2) +
pow(initialPos.y() - 1.5, 2) +
pow(initialPos.z() - 3.5, 2));
}
// 3D distance from x > 1.5 and y < -1.5
else if ((1.5 * m < initialPos.x()) && (-1.5 * m > initialPos.y()) &&
(initialPos.z() > 3.5 * m)) {
depth = std::sqrt(pow(initialPos.x() - 1.5, 2) +
pow(initialPos.y() - 1.5, 2) +
pow(initialPos.z() - 3.5, 2));
}
// 3D distance from x < -1.5 and y < -1.5
else if ((-1.5 * m > initialPos.x()) && (-1.5 * m > initialPos.y()) &&
(initialPos.z() > 3.5 * m)) {
depth = std::sqrt(pow(initialPos.x() - 1.5, 2) +
pow(initialPos.y() - 1.5, 2) +
pow(initialPos.z() - 3.5, 2));
}
// 3D distance from x < -1.5 and y > 1.5
else if ((-1.5 * m > initialPos.x()) && (1.5 * m < initialPos.y()) &&
(initialPos.z() > 3.5 * m)) {
depth = std::sqrt(pow(initialPos.x() - 1.5, 2) +
pow(initialPos.y() - 1.5, 2) +
pow(initialPos.z() - 3.5, 2));
}
else {
G4cerr << "Warning: No valid depth calculation case for position ("
<< initialPos.x() / m << ", " << initialPos.y() / m << ", "
<< initialPos.z() / m << ")" << G4endl;
}
// Get muon info
G4double muonTheta = -1;
G4double muonPhi = -1;
G4double muonInitialEnergy = -1;
auto muonIter = fMuonAngles.find(eventID);
if (muonIter != fMuonAngles.end()) {
muonTheta = muonIter->second.theta;
muonPhi = muonIter->second.phi;
muonInitialEnergy = muonIter->second.initialEnergy;
}
// Write data
if (fEscapedNeutronData.is_open()) {
fEscapedNeutronData << std::fixed << std::setprecision(15)
<< trackID << ","
<< eventID << ","
<< trackIter->second.creationEnergy / MeV << ","
<< exitEnergy / MeV << ","
<< crossingPoint.x() / mm << ","
<< crossingPoint.y() / mm << ","
<< crossingPoint.z() / mm << ","
<< initialPos.x() / mm << ","
<< initialPos.y() / mm << ","
<< initialPos.z() / mm << ","
<< depth / mm << ","
<< trackLength / mm << ","
<< muonTheta << ","
<< muonPhi << ","
<< muonInitialEnergy / GeV << ","
<< (creatorProcess ? creatorProcess->GetProcessName() : "Unknown");
fEscapedNeutronData << std::endl;
}
// Update count and cleanup
if (RunAction* runAction = dynamic_cast<RunAction*>
(const_cast<G4UserRunAction*>(G4RunManager::GetRunManager()->GetUserRunAction()))) {
runAction->AddNeutronEscaped();
}
fNeutronTracks.erase(trackID); // Remove after successful processing
}
catch (const std::exception& e) {
G4cerr << "Error processing neutron escape: " << e.what() << G4endl;
}
track->SetTrackStatus(fStopAndKill);
}
// Outer boundary crossing (20m)
else if ((((prePos.x() >= 1.5 * m || prePos.x() <= -1.5 * m) &&
(-6.5 * m <= prePos.y() && prePos.y() <= 6.5 * m) &&
(prePos.z() >= 0 * m)) ||
((prePos.y() >= 1.5 * m || prePos.y() <= -1.5 * m) &&
(-6.5 * m <= prePos.x() && prePos.x() <= 6.5 * m) &&
(prePos.z() >= 0 * m))) &&
((postPos.x() > 6.5 * m || postPos.x() < -6.5 * m) ||
(postPos.y() > 6.5 * m || postPos.y() < -6.5 * m)) &&
(postPos.z() >= 0 * m)) {
fNeutronTracks.erase(trackID); // Clean up neutrons that exit at outer boundary
}
}
}