Geant4 failing to track boundary crossing

Particles, Track, Event, Run and Biasing
1

_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
        }
    }
}
2

You’ll need to post more info than this. Specifically, your stepping action. In general the plan is usually to check that the step encounters a boundary. Similar for before or after the boundary:

step->GetPostStepPoint()->GetStepStatus()==fGeomBoundary

Then usually a check that (in this case) the step after is the lab:

step->GetPostStepPoint()->GetTouchableHandle()->GetVolume()->GetLogicalVolume()!= fLabVolume)

where fLabVolume would be passed from detector constructor which would be the logical volume of the lab. Or if it is everything not the rock is the “lab” than you can do similar logic for checking that you are not in the rock with the presteppoint.

3

I’m a little bit confused by your comment. I have included the entire code from my SteppingAction.cc file in the original post.

If it helps, this is how I have design my shell of rock.

// Dimensions for the rectangular shell
G4double innerLengthX = 3.0 * m;
G4double innerLengthY = 3.0 * m;
G4double innerLengthZ = 3.5 * m;

G4double thicknessTop = 5.0 * m;     // Thickness above the lab space
G4double thicknessSide = 5.0 * m;    // Thickness on sides
G4double thicknessBelow = 0.0 * m;   // No thickness below the lab space

// Outer box dimensions
G4double outerLengthX = innerLengthX + 2.0 * thicknessSide;
G4double outerLengthY = innerLengthY + 2.0 * thicknessSide;
G4double outerLengthZ = innerLengthZ + thicknessTop; // Only add thickness to top

// Create the outer box (positioned so that bottom of inner space is at z=0)
G4Box* outerBox = new G4Box("OuterBox", outerLengthX / 2, outerLengthY / 2,
    (innerLengthZ + thicknessTop) / 2);

// Create the inner box (positioned to center the lab space)
G4Box* innerBox = new G4Box("InnerBox", innerLengthX / 2, innerLengthY / 2,
    innerLengthZ / 2);

// Subtract the inner box from the outer box to create the shell
G4SubtractionSolid* rectangularShell =
    new G4SubtractionSolid("rShell", outerBox, innerBox,
        0, G4ThreeVector(0, 0, -thicknessTop / 2));

// Define logical volume for the shell
G4LogicalVolume* shellLV = new G4LogicalVolume(rectangularShell, standardRock, "rShell");

// place the lab space on the equator
G4VPhysicalVolume* shellPV =
    new G4PVPlacement(0, G4ThreeVector(0, 0, 4.25*m), shellLV, "rectShell", worldLV, false, 0, true);

 
G4VisAttributes* visShell = new G4VisAttributes(G4Colour::Red());
shellLV->SetVisAttributes(visShell);
4

That did not load properly the first time, apologies. As a formality, I’d replace

if (track->GetDefinition()->GetParticleName() != "neutron") return;

with

if(track->GetDefinition() == G4Neutron::Definition())return;

It is faster and less prone to errors from physics list interactions.

I suspect your issue is here and I alluded to it earlier:

 if (preVolume->GetName() == "rectShell" && postVolume->GetName() != "rectShell")

Geant4 tracks will always “stop” at a boundary which I suspect might cause unstable behavior when calling for their name (it might return null or the mother volume which might be the world). Replace with an explicit fGeomBoundary check and see if that works:

if((step->GetPreStepPoint()->GetStepStatus()==fGeomBoundary) && (postVolume->GetName() = "rectShell"))

This will specifically evaluate to True when leaving a boundary and entering the shell (i.e. first step inside).

5

Fantastic, thank you. I will give this a try!