|
| 1 | +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. |
| 2 | +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. |
| 3 | +// All rights not expressly granted are reserved. |
| 4 | +// |
| 5 | +// This software is distributed under the terms of the GNU General Public |
| 6 | +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". |
| 7 | +// |
| 8 | +// In applying this license CERN does not waive the privileges and immunities |
| 9 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 10 | +// or submit itself to any jurisdiction. |
| 11 | + |
| 12 | +#include "VecGeomG4PropagatingNavigator.h" |
| 13 | + |
| 14 | +#include "G4TouchableHistory.hh" |
| 15 | +#include "G4VPhysicalVolume.hh" |
| 16 | + |
| 17 | +#include <VecGeom/base/Transformation3D.h> |
| 18 | +#include <VecGeom/management/GeoManager.h> |
| 19 | +#include <VecGeom/navigation/GlobalLocator.h> |
| 20 | +#include <VecGeom/navigation/VNavigator.h> |
| 21 | +#include <VecGeom/navigation/VSafetyEstimator.h> |
| 22 | +#include <VecGeom/volumes/LogicalVolume.h> |
| 23 | +#include <VecGeom/volumes/PlacedVolume.h> |
| 24 | + |
| 25 | +#include <fairlogger/Logger.h> |
| 26 | + |
| 27 | +#include <cmath> |
| 28 | + |
| 29 | +namespace |
| 30 | +{ |
| 31 | +/// Geant4 abandons a track after ~50 zero steps, so a stalled step is nudged forward by this much |
| 32 | +/// (in Geant4 units) rather than returned as zero. |
| 33 | +const double kNudge = 1.e-3; |
| 34 | +} // namespace |
| 35 | + |
| 36 | +namespace o2::simsetup |
| 37 | +{ |
| 38 | + |
| 39 | +VecGeomG4PropagatingNavigator::VecGeomG4PropagatingNavigator(VecGeomG4Map const& map, bool zeroSafety) |
| 40 | + : VecGeomG4NavigatorBase(map), mZeroSafety(zeroSafety) |
| 41 | +{ |
| 42 | +} |
| 43 | + |
| 44 | +VecGeomG4PropagatingNavigator::~VecGeomG4PropagatingNavigator() |
| 45 | +{ |
| 46 | + LOG(info) << "VecGeom navigation: " << mNudgedSteps << " stalled steps nudged forward, " << mGlobalRelocates |
| 47 | + << " relocations restarted from the world"; |
| 48 | +} |
| 49 | + |
| 50 | +G4double VecGeomG4PropagatingNavigator::ComputeStep(const G4ThreeVector& globalPoint, const G4ThreeVector& direction, |
| 51 | + const G4double proposedStepLength, G4double& newSafety) |
| 52 | +{ |
| 53 | + newSafety = 0.; |
| 54 | + mLastDirection = direction; |
| 55 | + mHaveNextState = false; |
| 56 | + |
| 57 | + auto const* top = mCurState.Top(); |
| 58 | + if (top == nullptr) { // the track is outside the world |
| 59 | + mWouldEnter = mWouldExit = false; |
| 60 | + return kInfinity; |
| 61 | + } |
| 62 | + auto const* navigator = top->GetLogicalVolume()->GetNavigator(); |
| 63 | + |
| 64 | + double limit = proposedStepLength * kG4ToVG; |
| 65 | + if (!(limit < vecgeom::kInfLength)) { |
| 66 | + limit = vecgeom::kInfLength; |
| 67 | + } |
| 68 | + |
| 69 | + // VecGeom's own combined entry point, rather than a step followed by a relocation of our own: it |
| 70 | + // is the one that knows how to descend through an assembly, whose placed volume can never be the |
| 71 | + // answer because it has no DistanceToOut. |
| 72 | + double vgSafety = 0.; |
| 73 | + const double vgStep = navigator->ComputeStepAndSafetyAndPropagatedState( |
| 74 | + toVG(globalPoint), toDir(direction), limit, mCurState, mNextState, !mOnBoundary && !mZeroSafety, vgSafety); |
| 75 | + mOnBoundary = false; |
| 76 | + newSafety = (vgSafety > 0. && !mZeroSafety) ? vgSafety * kVGToG4 : 0.; |
| 77 | + |
| 78 | + G4double step = vgStep * kVGToG4; |
| 79 | + if (mNextState.IsOnBoundary()) { |
| 80 | + // Entering a daughter deepens the state, possibly by more than one level when an assembly |
| 81 | + // stands in between; leaving the current volume does not. |
| 82 | + mWouldEnter = mNextState.GetCurrentLevel() > mCurState.GetCurrentLevel(); |
| 83 | + mWouldExit = !mWouldEnter; |
| 84 | + mNextPoint = globalPoint + step * direction; |
| 85 | + mHaveNextState = true; |
| 86 | + } else { |
| 87 | + mWouldEnter = mWouldExit = false; |
| 88 | + step = kInfinity; |
| 89 | + } |
| 90 | + |
| 91 | + if (vgStep < 0.) { |
| 92 | + // A negative distance means the state and the point disagree. Nudge forward and relocate from |
| 93 | + // the world on the next call rather than propagating the inconsistency. |
| 94 | + mForceReInit = true; |
| 95 | + mHaveNextState = false; |
| 96 | + ++mNudgedSteps; |
| 97 | + mZeroSteps = 0; |
| 98 | + return kNudge; |
| 99 | + } |
| 100 | + if (step < 1.e-10) { |
| 101 | + if (++mZeroSteps > 4) { |
| 102 | + mForceReInit = true; |
| 103 | + mHaveNextState = false; |
| 104 | + ++mNudgedSteps; |
| 105 | + return kNudge; |
| 106 | + } |
| 107 | + } else { |
| 108 | + mForceReInit = false; |
| 109 | + mZeroSteps = 0; |
| 110 | + } |
| 111 | + return step; |
| 112 | +} |
| 113 | + |
| 114 | +G4VPhysicalVolume* VecGeomG4PropagatingNavigator::ResetHierarchyAndLocate(const G4ThreeVector&, const G4ThreeVector&, |
| 115 | + const G4TouchableHistory& history) |
| 116 | +{ |
| 117 | + // Geant4 hands back a touchable it saved earlier, e.g. when resuming a track whose secondaries |
| 118 | + // were followed first. The VecGeom state is rebuilt from it. |
| 119 | + fEnteredDaughter = false; |
| 120 | + fExitedMother = false; |
| 121 | + mWouldEnter = false; |
| 122 | + mWouldExit = false; |
| 123 | + mOnBoundary = false; |
| 124 | + mHaveNextState = false; |
| 125 | + fHistory = *history.GetHistory(); |
| 126 | + if (!stateFromHistory(mCurState) && fHistory.GetVolume(0) != nullptr) { |
| 127 | + LOG(fatal) << "Geant4 handed back a touchable that matches no VecGeom path"; |
| 128 | + } |
| 129 | + mPrevState = mCurState; |
| 130 | + return fHistory.GetTopVolume(); |
| 131 | +} |
| 132 | + |
| 133 | +G4VPhysicalVolume* VecGeomG4PropagatingNavigator::LocateGlobalPointAndSetup(const G4ThreeVector& point, |
| 134 | + const G4ThreeVector*, |
| 135 | + const G4bool relativeSearch, const G4bool) |
| 136 | +{ |
| 137 | + bool onBoundary = fWasLimitedByGeometry; |
| 138 | + if (mHaveNextState && point.diff2(mNextPoint) < 1.e-16) { |
| 139 | + onBoundary = true; |
| 140 | + } |
| 141 | + |
| 142 | + mPrevState = mCurState; |
| 143 | + |
| 144 | + if (!mForceReInit && relativeSearch && onBoundary && mHaveNextState) { |
| 145 | + // The state on the far side of the boundary was already worked out, and relocated, by the step |
| 146 | + // that found it. Adopting it is cheaper than locating again. |
| 147 | + mCurState = mNextState; |
| 148 | + } else if (mForceReInit || !relativeSearch || onBoundary) { |
| 149 | + mCurState.Clear(); |
| 150 | + vecgeom::GlobalLocator::LocateGlobalPoint(vecgeom::GeoManager::Instance().GetWorld(), toVG(point), mCurState, true); |
| 151 | + mForceReInit = false; |
| 152 | + ++mGlobalRelocates; |
| 153 | + } |
| 154 | + // Otherwise the point only moved inside the volume the state already names. |
| 155 | + |
| 156 | + auto* target = historyFromState(mCurState); |
| 157 | + mCrossed = onBoundary; |
| 158 | + if (onBoundary) { |
| 159 | + fExitedMother = mWouldExit; |
| 160 | + fEnteredDaughter = mWouldEnter; |
| 161 | + mOnBoundary = true; |
| 162 | + } |
| 163 | + mHaveNextState = false; |
| 164 | + return target; |
| 165 | +} |
| 166 | + |
| 167 | +void VecGeomG4PropagatingNavigator::LocateGlobalPointWithinVolume(const G4ThreeVector&) |
| 168 | +{ |
| 169 | + // The track moved inside the volume it is already in, so only the boundary flags change. |
| 170 | + mWouldEnter = false; |
| 171 | + mWouldExit = false; |
| 172 | + mOnBoundary = false; |
| 173 | + mCrossed = false; |
| 174 | + mHaveNextState = false; |
| 175 | + fEnteredDaughter = false; |
| 176 | + fExitedMother = false; |
| 177 | +} |
| 178 | + |
| 179 | +G4double VecGeomG4PropagatingNavigator::ComputeSafety(const G4ThreeVector& globalPoint, const G4double, const G4bool) |
| 180 | +{ |
| 181 | + if (mZeroSafety || mOnBoundary || mCrossed || fEnteredDaughter || fExitedMother || mWouldEnter || mWouldExit) { |
| 182 | + return 0.; |
| 183 | + } |
| 184 | + auto const* top = mCurState.Top(); |
| 185 | + if (top == nullptr) { |
| 186 | + return 0.; |
| 187 | + } |
| 188 | + const double safety = top->GetLogicalVolume()->GetSafetyEstimator()->ComputeSafety(toVG(globalPoint), mCurState); |
| 189 | + return (safety > 0.) ? safety * kVGToG4 : 0.; |
| 190 | +} |
| 191 | + |
| 192 | +G4ThreeVector VecGeomG4PropagatingNavigator::GetGlobalExitNormal(const G4ThreeVector& point, G4bool* valid) |
| 193 | +{ |
| 194 | + // The surface just crossed belongs to the volume that was left when the step exited a mother, and |
| 195 | + // to the volume that was entered when it entered a daughter. |
| 196 | + auto const& state = mWouldExit ? mPrevState : mCurState; |
| 197 | + auto const* volume = state.Top(); |
| 198 | + if (volume == nullptr) { |
| 199 | + *valid = false; |
| 200 | + return G4ThreeVector(0., 0., 1.); |
| 201 | + } |
| 202 | + |
| 203 | + vecgeom::Transformation3D m; |
| 204 | + state.TopMatrix(m); |
| 205 | + V3 normal; |
| 206 | + volume->Normal(m.Transform(toVG(point)), normal); |
| 207 | + V3 global = m.InverseTransformDirection(normal); |
| 208 | + |
| 209 | + // Oriented along the direction of motion, as TGeo's FindNormalFast does. |
| 210 | + const V3 dir = toDir(mLastDirection); |
| 211 | + if (global.Dot(dir) < 0.) { |
| 212 | + global = -global; |
| 213 | + } |
| 214 | + // VecGeom's Normal() also answers whether the point was on the surface; Geant4 hands back points a |
| 215 | + // few nanometres off the face, so only a degenerate vector is refused. |
| 216 | + const double mag2 = global.Mag2(); |
| 217 | + *valid = std::isfinite(mag2) && mag2 > 0.25; |
| 218 | + return G4ThreeVector(global[0], global[1], global[2]); |
| 219 | +} |
| 220 | + |
| 221 | +G4ThreeVector VecGeomG4PropagatingNavigator::GetLocalExitNormal(G4bool* valid) |
| 222 | +{ |
| 223 | + // By convention the local normal is expressed in the frame of the final volume. |
| 224 | + const G4ThreeVector global = GetGlobalExitNormal(mNextPoint, valid); |
| 225 | + vecgeom::Transformation3D m; |
| 226 | + mCurState.TopMatrix(m); |
| 227 | + const V3 local = m.TransformDirection(toDir(global)); |
| 228 | + return G4ThreeVector(local[0], local[1], local[2]); |
| 229 | +} |
| 230 | + |
| 231 | +} // namespace o2::simsetup |
0 commit comments