Skip to content

Commit 853ab5d

Browse files
sawenzelclaude
andcommitted
Block the exited volume only as Geant4 does
This fixes tracks entering the volume they just left without seeing its boundary. - Both VecGeom navigators kept the volume just left blocked for every following step, whatever the direction. - A track turned back into it at the exit point, as by multiple scattering, crossed it unseen in the mother. - In the hollow RB24 copper tubes this let shower electrons pass the wall and gave about 10 % more FDD-A hits. - The volume is now blocked only in the first step after the exit, and only while the direction points away from it, as in G4NormalNavigation. This follows the same fix in G4VecGeomNav. https://gitlab.cern.ch/swenzel/g4vecgeomnav/-/tree/fix/block-and-reflection Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
1 parent 24d46a0 commit 853ab5d

5 files changed

Lines changed: 91 additions & 12 deletions

File tree

‎Detectors/gconfig/src/VecGeomG4Navigator.cxx‎

Lines changed: 23 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -80,7 +80,7 @@ void VecGeomG4Navigator::locateFromWorld(const V3& point)
8080
/// surely as one at normal incidence, and it is kept small: TOF has layers 2.4e-8 cm apart, and a
8181
/// fixed push along the direction steps over them. It never goes below a thousand times the rounding
8282
/// of the largest coordinate, which matters far from the origin.
83-
double VecGeomG4Navigator::boundaryPush(const V3& point, const V3& dir) const
83+
double VecGeomG4Navigator::boundaryPush(const V3& point, const V3& dir)
8484
{
8585
constexpr double kMaxPush = 1.e-4; // cm along the direction
8686
const double big = std::max({std::abs(point.x()), std::abs(point.y()), std::abs(point.z())});
@@ -89,6 +89,7 @@ double VecGeomG4Navigator::boundaryPush(const V3& point, const V3& dir) const
8989
// The face just crossed: the entered daughter's, or the current volume's own.
9090
vecgeom::NavigationState const& st = mWouldEnter ? mNextState : mCurState;
9191
double cosn = 1.;
92+
mPushNormalValid = false;
9293
if (!st.IsOutside() && st.Top() != nullptr) {
9394
vecgeom::Transformation3D m;
9495
st.TopMatrix(m);
@@ -98,6 +99,10 @@ double VecGeomG4Navigator::boundaryPush(const V3& point, const V3& dir) const
9899
if (c > 0. && n.Mag2() > 0.5) {
99100
cosn = c;
100101
}
102+
if (n.Mag2() > 0.5) {
103+
mPushNormal = m.InverseTransformDirection(n);
104+
mPushNormalValid = true;
105+
}
101106
}
102107
return std::max(rounding, std::min(mPushDepth / cosn, kMaxPush));
103108
}
@@ -123,10 +128,20 @@ G4double VecGeomG4Navigator::ComputeStep(const G4ThreeVector& globalPoint, const
123128
newSafety = mLastSafety;
124129
}
125130

126-
// The step is computed on a copy. The state, and the volume blocked by the last exit, stay as
127-
// the locate left them for every call until the next locate, which is what the field propagator
128-
// relies on when it calls this from trial points along the curve.
131+
// The step is computed on a copy. The state stays as the locate left it for every call until the
132+
// next locate, which is what the field propagator relies on when it calls this from trial points
133+
// along the curve. The volume the last crossing left is blocked in the first call only, and only
134+
// while the direction points away from it; a track turning back into it must see its boundary.
129135
mStepState = mCurState;
136+
if (mExitBlockPending) {
137+
const bool block = samePoint(globalPoint, mLastLocatedPoint) &&
138+
(mExitNormalFromPush ? directionLeaves(mPushNormal, toDir(direction))
139+
: directionLeaves(mExitedState, toVG(globalPoint), toDir(direction)));
140+
if (!block) {
141+
mStepState.SetLastExited(mEmptyState.GetLastExitedState());
142+
}
143+
clearLastExited();
144+
}
130145
const double limit = std::min(proposedStepLength * kG4ToVG, static_cast<double>(vecgeom::kInfLength));
131146
double safety = 0.;
132147
double vgStep = navigator->ComputeStepAndSafety(toVG(globalPoint), toDir(direction), limit, mStepState, calcSafety,
@@ -204,6 +219,7 @@ void VecGeomG4Navigator::leaveFlushVolumes(const V3& point, const V3& dir, int m
204219
}
205220
left = mCurState.Top();
206221
mCurState.SetLastExited();
222+
setExited(mCurState, false);
207223
if (mCurState.GetCurrentLevel() <= 1) {
208224
mCurState.Clear(); // nothing to travel in even in the world: the track left it
209225
return;
@@ -307,8 +323,8 @@ G4VPhysicalVolume* VecGeomG4Navigator::LocateGlobalPointAndSetup(const G4ThreeVe
307323
mLocatedOnBoundary = true;
308324
} else if (crossing) {
309325
// Out of the current volume: up until the point is contained, then down, never back into the
310-
// volume just left; that volume stays blocked at zero distance for the next step, as
311-
// G4Navigator's fBlockedPhysicalVolume.
326+
// volume just left; that volume is blocked at zero distance in the next ComputeStep while the
327+
// direction points away from it, as G4Navigator's fBlockedPhysicalVolume.
312328
mReloScratch = mCurState;
313329
if (mCurState.GetCurrentLevel() <= 1) {
314330
mCurState.Clear(); // left the world
@@ -318,6 +334,7 @@ G4VPhysicalVolume* VecGeomG4Navigator::LocateGlobalPointAndSetup(const G4ThreeVe
318334
vecgeom::GlobalLocator::RelocatePointFromPathForceDifferent(m.Transform(q), mCurState);
319335
mReloScratch.SetLastExited();
320336
mCurState.SetLastExited(mReloScratch.GetLastExitedState());
337+
setExited(mReloScratch, true);
321338
}
322339
mLocatedOnBoundary = true;
323340
} else {

‎Detectors/gconfig/src/VecGeomG4Navigator.h‎

Lines changed: 24 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -24,9 +24,10 @@ namespace o2::simsetup
2424
///
2525
/// - ComputeStep leaves the current volume alone. It records whether the step ends on a boundary
2626
/// and, for a daughter hit, the state that enters it.
27-
/// - The locate on a boundary relocates: into the recorded daughter, or out of the current volume,
28-
/// which then stays blocked for the next step. The point is first pushed across the face by a
29-
/// small depth, and afterwards leaves every volume it is flush with and heading out of.
27+
/// - The locate on a boundary relocates: into the recorded daughter, or out of the current volume.
28+
/// The point is first pushed across the face by a small depth, and afterwards leaves every volume
29+
/// it is flush with and heading out of. As in Geant4, the volume left is blocked only in the first
30+
/// ComputeStep after the exit, and only while the direction points away from it.
3031
/// - Safety is zero only at the boundary point itself.
3132
///
3233
/// Unlike TG4VecGeomNavigator, the VecGeom geometry is converted from TGeo, not from Geant4, so one
@@ -66,9 +67,21 @@ class VecGeomG4Navigator : public VecGeomG4NavigatorBase
6667
void leaveFlushVolumes(const V3& point, const V3& dir, int minLevel, vecgeom::VPlacedVolume const* avoid);
6768
/// Sets fEnteredDaughter and fExitedMother from the paths before (mReloScratch) and after a crossing.
6869
void updateCrossingFlags(bool entering);
69-
void clearLastExited() { mCurState.SetLastExited(mEmptyState.GetLastExitedState()); }
70-
/// How far (cm) a boundary point is pushed along the direction before it is located.
71-
double boundaryPush(const V3& point, const V3& dir) const;
70+
void clearLastExited()
71+
{
72+
mCurState.SetLastExited(mEmptyState.GetLastExitedState());
73+
mExitBlockPending = false;
74+
}
75+
/// Records the top of \a st as the volume just left, blocked in the next ComputeStep by Geant4's rule.
76+
void setExited(vecgeom::NavigationState const& st, bool normalFromPush)
77+
{
78+
mExitedState = st;
79+
mExitBlockPending = true;
80+
mExitNormalFromPush = normalFromPush && mPushNormalValid;
81+
}
82+
/// How far (cm) a boundary point is pushed along the direction before it is located. Keeps the
83+
/// global normal of the face in mPushNormal, for the exit block.
84+
double boundaryPush(const V3& point, const V3& dir);
7285
/// The normal of the surface the last geometry-limited ComputeStep ended on, global, unit length.
7386
bool computeExitNormal(const G4ThreeVector& point, V3& globalNormal) const;
7487

@@ -83,6 +96,11 @@ class VecGeomG4Navigator : public VecGeomG4NavigatorBase
8396
vecgeom::NavigationState mHistoryState; ///< the state fHistory was built from
8497
vecgeom::NavigationState mEmptyState; ///< permanently empty; its last-exited entry clears others
8598
vecgeom::NavigationState mNormalState; ///< the volume whose surface the last boundary step ended on
99+
vecgeom::NavigationState mExitedState; ///< the volume the last crossing left
100+
bool mExitBlockPending = false; ///< the next ComputeStep is the first after that exit
101+
bool mExitNormalFromPush = false; ///< mPushNormal is the outward normal of the volume left
102+
V3 mPushNormal; ///< global normal of the face the last push crossed
103+
bool mPushNormalValid = false;
86104
bool mHistoryValid = false;
87105

88106
bool mWouldEnter = false; ///< the last ComputeStep ends by entering a daughter

‎Detectors/gconfig/src/VecGeomG4NavigatorBase.h‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,8 +18,10 @@
1818
#include "G4SystemOfUnits.hh"
1919
#include "G4ThreeVector.hh"
2020

21+
#include <VecGeom/base/Transformation3D.h>
2122
#include <VecGeom/base/Vector3D.h>
2223
#include <VecGeom/navigation/NavigationState.h>
24+
#include <VecGeom/volumes/PlacedVolume.h>
2325

2426
namespace o2::simsetup
2527
{
@@ -48,6 +50,31 @@ class VecGeomG4NavigatorBase : public G4Navigator
4850
/// False if a level matches no VecGeom placement; \a state is then empty.
4951
bool stateFromHistory(vecgeom::NavigationState& state) const;
5052

53+
/// Geant4's rule for the volume a track just left (G4NormalNavigation, G4VoxelNavigation): it is
54+
/// blocked only while the direction points away from it, along its outward normal \a n (global).
55+
[[gnu::always_inline]] static bool directionLeaves(const V3& n, const V3& dir)
56+
{
57+
constexpr double kMinExitingNormalCosine = 1e-3; // as in G4NormalNavigation
58+
return n.Dot(dir) >= kMinExitingNormalCosine;
59+
}
60+
61+
/// The same, with the normal taken from \a exited, which has the volume left on top.
62+
[[gnu::always_inline]] static bool directionLeaves(vecgeom::NavigationState const& exited, const V3& point,
63+
const V3& dir)
64+
{
65+
if (exited.IsOutside() || exited.Top() == nullptr) {
66+
return false;
67+
}
68+
vecgeom::Transformation3D m;
69+
exited.TopMatrix(m);
70+
V3 n;
71+
exited.Top()->GetUnplacedVolume()->Normal(m.Transform(point), n);
72+
if (!(n.Mag2() > 0.5)) {
73+
return false; // no valid exit normal: Geant4 does not block either
74+
}
75+
return directionLeaves(m.InverseTransformDirection(n), dir);
76+
}
77+
5178
VecGeomG4Map const& mMap;
5279
};
5380

‎Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@ namespace o2::simsetup
3939
VecGeomG4PropagatingNavigator::VecGeomG4PropagatingNavigator(VecGeomG4Map const& map, bool zeroSafety)
4040
: VecGeomG4NavigatorBase(map), mZeroSafety(zeroSafety)
4141
{
42+
mEmptyState.Clear();
4243
}
4344

4445
VecGeomG4PropagatingNavigator::~VecGeomG4PropagatingNavigator()
@@ -69,6 +70,16 @@ G4double VecGeomG4PropagatingNavigator::ComputeStep(const G4ThreeVector& globalP
6970
// VecGeom's own combined entry point, rather than a step followed by a relocation of our own: it
7071
// is the one that knows how to descend through an assembly, whose placed volume can never be the
7172
// answer because it has no DistanceToOut.
73+
// The adopted state marks the volume the crossing left. As in Geant4, it is blocked in the first
74+
// step after the exit only, and only while the direction points away from it; a track turning
75+
// back into it must see its boundary.
76+
const bool block = mExitBlockPending && globalPoint.diff2(mLocatedPoint) < 1.e-20 &&
77+
directionLeaves(mPrevState, toVG(globalPoint), toDir(direction));
78+
mExitBlockPending = false;
79+
if (!block) {
80+
mCurState.SetLastExited(mEmptyState.GetLastExitedState());
81+
}
82+
7283
double vgSafety = 0.;
7384
const double vgStep = navigator->ComputeStepAndSafetyAndPropagatedState(
7485
toVG(globalPoint), toDir(direction), limit, mCurState, mNextState, !mOnBoundary && !mZeroSafety, vgSafety);
@@ -140,11 +151,14 @@ G4VPhysicalVolume* VecGeomG4PropagatingNavigator::LocateGlobalPointAndSetup(cons
140151
}
141152

142153
mPrevState = mCurState;
154+
mLocatedPoint = point;
155+
mExitBlockPending = false;
143156

144157
if (!mForceReInit && relativeSearch && onBoundary && mHaveNextState) {
145158
// The state on the far side of the boundary was already worked out, and relocated, by the step
146159
// that found it. Adopting it is cheaper than locating again.
147160
mCurState = mNextState;
161+
mExitBlockPending = mWouldExit;
148162
} else if (mForceReInit || !relativeSearch || onBoundary) {
149163
mCurState.Clear();
150164
vecgeom::GlobalLocator::LocateGlobalPoint(vecgeom::GeoManager::Instance().GetWorld(), toVG(point), mCurState, true);

‎Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,7 @@ class VecGeomG4PropagatingNavigator : public VecGeomG4NavigatorBase
5050
vecgeom::NavigationState mCurState; ///< where the track is now
5151
vecgeom::NavigationState mNextState; ///< where the last computed step would put it
5252
vecgeom::NavigationState mPrevState; ///< where it was before the last boundary crossing
53+
vecgeom::NavigationState mEmptyState; ///< permanently empty; its last-exited entry clears others
5354

5455
bool mZeroSafety = false;
5556
bool mHaveNextState = false; ///< mNextState holds the result of a geometry-limited step
@@ -62,6 +63,8 @@ class VecGeomG4PropagatingNavigator : public VecGeomG4NavigatorBase
6263
bool mOnBoundary = false; ///< the current point sits on a boundary
6364
bool mForceReInit = false; ///< next locate must start from the world, the state is suspect
6465
bool mCrossed = false; ///< the last locate acted on a boundary crossing
66+
bool mExitBlockPending = false; ///< the next step is the first after leaving mPrevState's volume
67+
G4ThreeVector mLocatedPoint{-1e8, -1e8, -1e8}; ///< where the last locate put the track
6568

6669
int mZeroSteps = 0;
6770
long mNudgedSteps = 0;

0 commit comments

Comments
 (0)