Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!
Geant4 Version: 11.3.2
Operating System: Ubuntu 24.04.4 LTS, run on WSL
Compiler/Version: g++ 13.3.0
CMake Version: 3.28.3
I have a Biopsy sample that I am modelling radioactive decay within, and to mimic the surrounding tissue, I have created a periodic boundary condition usng a custom PostStepDoIt to wrap interacting particles back into to the volume.
The condition is working as expected, in that the particle exits one side, reappears at the other and continues its journey with no energy lost during the jump.
However, when I visualise the particle I can see a track through the sample during the jump. Is there a way of turning off the trajectory just for this part?
Below is my PostStepDoIt section;
G4VParticleChange* PeriodicBoundaryProcess::PostStepDoIt(const G4Track& aTrack,
const G4Step& aStep)
{
// Safety check
if (!fDimensionsFound) {
FindBiopsyDimensions();
// If still not found, can't apply PBC
if (!fDimensionsFound) return &fParticleChange;
}
fParticleChange.Initialize(aTrack);
const G4StepPoint* postStepPoint = aStep.GetPostStepPoint();
//Only act at geometry boundary
if (postStepPoint->GetStepStatus() != fGeomBoundary) return &fParticleChange;
fOldPosition = postStepPoint->GetPosition();
G4double tol = 100.0 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
bool onBoundaryX = fPeriodicX && (std::abs(std::abs(fOldPosition.x()) - fHalfLengths.x()) < tol);
bool onBoundaryY = fPeriodicY && (std::abs(std::abs(fOldPosition.y()) - fHalfLengths.y()) < tol);
bool onBoundaryZ = fPeriodicZ && (std::abs(std::abs(fOldPosition.z()) - fHalfLengths.z()) < tol);
// If internal hit
if (!onBoundaryX && !onBoundaryY && !onBoundaryZ) {
// Silent return
return &fParticleChange;
}
// WRAPPING LOGIC
const G4DynamicParticle* dp = aTrack.GetDynamicParticle();
fOldMomentum = dp->GetMomentumDirection();
fNewPosition = fOldPosition;
fNewMomentum = fOldMomentum;
G4bool valid;
G4ThreeVector normal = G4TransportationManager::GetTransportationManager()
->GetNavigatorForTracking()
->GetGlobalExitNormal(fNewPosition, &valid);
if (!valid) return &fParticleChange;
//Push particle strictly inside volume
G4double epsilon = std::max(1e-5*mm, 100.0*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance());
bool wrapped = false;
G4String wrappedAxes = "";
// --- X ---
if (onBoundaryX && std::abs(normal.x()) > 0.5) {
fNewPosition.setX(-fOldPosition.x() + (normal.x() > 0 ? +epsilon : -epsilon));
wrapped = true;
wrappedAxes += "X";
}
// --- Y ---
if (onBoundaryY && std::abs(normal.y()) > 0.5) {
fNewPosition.setY(-fOldPosition.y() + (normal.y() > 0 ? +epsilon : -epsilon));
wrapped = true;
wrappedAxes += (wrappedAxes.empty() ? "Y" : ",Y");
}
// --- Z ---
if (onBoundaryZ && std::abs(normal.z()) > 0.5) {
fNewPosition.setZ(-fOldPosition.z() + (normal.z() > 0 ? +epsilon : -epsilon));
wrapped = true;
wrappedAxes += (wrappedAxes.empty() ? "Z" : ",Z");
}
if (!wrapped) {
return &fParticleChange;
}
// Update the ParticleChange
fParticleChange.ProposePosition(fNewPosition);
fParticleChange.ProposeMomentumDirection(fNewMomentum);
G4Navigator* nav = G4TransportationManager::GetTransportationManager()
->GetNavigatorForTracking();
nav->LocateGlobalPointAndSetup(fNewPosition, &fNewMomentum, false, false);
aTrack.GetStep()->GetPostStepPoint()->SetStepStatus(fUndefined);
return &fParticleChange;
}