From 09cd96ed6ab37730e0bfd66a9637a3bfd4dee975 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 28 Sep 2026 21:39:21 +0200 Subject: [PATCH 1/5] Set up the VecGeom BVH navigator and safety estimators for every volume This lets GeometryManager set up VecGeom for Geant4 navigation as well as for the material budget. - buildVecGeomGeometry() converts once and takes the assembly flattening as a parameter; the material budget keeps flattening. - With a VecGeom that has BVHNavigatorV, volumes with more than two daughters get the BVH navigator and level locator, and every volume gets an explicit safety estimator. - Without it the VecGeom v2 setup is unchanged. https://gitlab.cern.ch/VecGeom/VecGeom/-/merge_requests/1547 Co-Authored-By: Claude Opus 5.5 --- .../include/DetectorsBase/GeometryManager.h | 6 +++ Detectors/Base/src/GeometryManager.cxx | 39 ++++++++++++++++--- 2 files changed, 39 insertions(+), 6 deletions(-) diff --git a/Detectors/Base/include/DetectorsBase/GeometryManager.h b/Detectors/Base/include/DetectorsBase/GeometryManager.h index 93f3931e203d4..d325a68929dd7 100644 --- a/Detectors/Base/include/DetectorsBase/GeometryManager.h +++ b/Detectors/Base/include/DetectorsBase/GeometryManager.h @@ -135,6 +135,12 @@ class GeometryManager : public TObject /// Mean material budget between two points, using the VecGeom backend. On first call, /// lazily converts the currently loaded TGeo geometry to VecGeom (once per process). static o2::base::MatBudget vecGeomMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1); + /// Converts the currently loaded TGeo geometry to VecGeom and assigns a navigator, a level + /// locator and a safety estimator to every logical volume. Does the work once per process; later + /// calls, whatever they ask for, return the geometry already built, so a caller that needs a + /// particular assembly treatment must come first. \param flattenAssemblies dissolves TGeo + /// assemblies into their content. + static void buildVecGeomGeometry(bool flattenAssemblies); #else static constexpr bool isVecGeomAvailable() { return false; } #endif diff --git a/Detectors/Base/src/GeometryManager.cxx b/Detectors/Base/src/GeometryManager.cxx index 5d6a8def8e7c3..a3b61217aee01 100644 --- a/Detectors/Base/src/GeometryManager.cxx +++ b/Detectors/Base/src/GeometryManager.cxx @@ -49,6 +49,14 @@ #include #include #include +// The BVH navigator of the VNavigator family, which Geant4 navigation needs on every volume. +#if __has_include() +#define O2_VECGEOM_HAS_BVH_VNAVIGATOR +#include +#include +#include +#include +#endif #endif using namespace o2::detectors; @@ -574,15 +582,22 @@ bool usesBvhAcceleration(vecgeom::LogicalVolume const* vol) /// process, the first time the VecGeom backend is requested. Not part of loadGeometry(), /// which every job calls regardless of whether it ever uses the VecGeom backend. void ensureVecGeomWorldBuilt() +{ + GeometryManager::buildVecGeomGeometry(true); +} +} // namespace + +void GeometryManager::buildVecGeomGeometry(bool flattenAssemblies) { static std::once_flag onceFlag; - std::call_once(onceFlag, []() { + std::call_once(onceFlag, [flattenAssemblies]() { if (!gGeoManager) { LOG(fatal) << "Cannot build VecGeom geometry: no TGeo geometry loaded (call GeometryManager::loadGeometry() first)"; } // Translate geometry and material pointers, then build acceleration structures. tgeo2vecgeom::RootGeoManager::Instance().SetMaterialConversionHook([](TGeoMaterial const* m) { return (void*)m; }); - tgeo2vecgeom::RootGeoManager::Instance().SetFlattenAssemblies(true); + LOG(info) << "VecGeom conversion: flattenAssemblies=" << flattenAssemblies; + tgeo2vecgeom::RootGeoManager::Instance().SetFlattenAssemblies(flattenAssemblies); tgeo2vecgeom::RootGeoManager::Instance().LoadRootGeometry(); // Acceleration structures must be built before the navigators/locators reference them. @@ -593,15 +608,28 @@ void ensureVecGeomWorldBuilt() // Builds a BVH per logical volume. vecgeom::BVHManager::Init(); - // For each logical volume, set both a navigator (used for ComputeStep) and a matched - // level locator (used for point relocation after a boundary crossing via GlobalLocator). + // For each logical volume, set a navigator (used for ComputeStep), a matched level locator + // (used for point relocation after a boundary crossing via GlobalLocator) and, where the + // VNavigator family is complete, the safety estimator LogicalVolume::GetSafetyEstimator() + // hands out, which is separate from the one a navigator uses internally. for (auto& lvol : vecgeom::GeoManager::Instance().GetLogicalVolumesMap()) { auto* vol = lvol.second; if (!usesBvhAcceleration(vol)) { vol->SetNavigator(vecgeom::NewSimpleNavigator<>::Instance()); +#ifdef O2_VECGEOM_HAS_BVH_VNAVIGATOR + vol->SetLevelLocator(vol->ContainsAssembly() ? vecgeom::SimpleAssemblyLevelLocator::GetInstance() + : vecgeom::SimpleLevelLocator::GetInstance()); + vol->SetSafetyEstimator(vecgeom::SimpleSafetyEstimator::Instance()); +#else vol->SetLevelLocator(vecgeom::SimpleLevelLocator::GetInstance()); +#endif } else { -#if VECGEOM_VERSION >= 0x020000 +#if defined(O2_VECGEOM_HAS_BVH_VNAVIGATOR) + vol->SetNavigator(vecgeom::BVHNavigatorV<>::Instance()); + vol->SetLevelLocator(vol->ContainsAssembly() ? vecgeom::BVHAssemblyAwareLevelLocator::GetInstance() + : vecgeom::BVHLevelLocator::GetInstance()); + vol->SetSafetyEstimator(vecgeom::BVHSafetyEstimator::Instance()); +#elif VECGEOM_VERSION >= 0x020000 // VecGeom 2 turned BVHNavigator into a plain class with static entry points instead of a // VNavigator singleton, so there is nothing to attach: vecGeomMaterialBudget() calls it // directly. @@ -627,7 +655,6 @@ void ensureVecGeomWorldBuilt() } }); } -} // namespace //_____________________________________________________________________________________ o2::base::MatBudget GeometryManager::vecGeomMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1) From 9ae3442076ad91897a97b688f8320c11b30fd321 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Tue, 29 Sep 2026 08:25:12 +0200 Subject: [PATCH 2/5] Add a VecGeom navigation mode to the Geant4 engine This adds G4.navmode=kVecGeom, in which Geant4 answers every navigation query from VecGeom. - The geometry, materials and touchables stay the ones g4root builds from TGeo. - VecGeomG4Map pairs each VecGeom placement with the chain of g4root volumes it stands for, so volume ids, copy numbers and CurrentVolOffID are unchanged. - VecGeomG4Navigator relocates at the boundary locate, with the volume just left blocked for the next step and the point pushed across the face by a small depth, as G4VecGeomNav's TG4VecGeomNavigator does. - G4.vecgeomCheckRays, vecgeomCheckLocation and vecgeomCheckVolumes compare VecGeom with TGeo before transport. - It is built only when TGeo2VecGeom and a VecGeom with BVHNavigatorV are found. https://gitlab.cern.ch/VecGeom/g4vecgeomnav/-/merge_requests/25 Co-Authored-By: Claude Opus 5.5 --- Common/SimConfig/include/SimConfig/G4Params.h | 19 +- Detectors/gconfig/CMakeLists.txt | 32 +- Detectors/gconfig/g4Config.C | 9 + .../include/SimSetup/VecGeomNavigation.h | 32 ++ Detectors/gconfig/src/VecGeomChecks.cxx | 257 +++++++++ Detectors/gconfig/src/VecGeomChecks.h | 37 ++ Detectors/gconfig/src/VecGeomG4Map.cxx | 187 +++++++ Detectors/gconfig/src/VecGeomG4Map.h | 124 +++++ Detectors/gconfig/src/VecGeomG4Navigator.cxx | 496 ++++++++++++++++++ Detectors/gconfig/src/VecGeomG4Navigator.h | 116 ++++ .../gconfig/src/VecGeomG4NavigatorBase.cxx | 128 +++++ .../gconfig/src/VecGeomG4NavigatorBase.h | 56 ++ Detectors/gconfig/src/VecGeomNavigation.cxx | 137 +++++ 13 files changed, 1626 insertions(+), 4 deletions(-) create mode 100644 Detectors/gconfig/include/SimSetup/VecGeomNavigation.h create mode 100644 Detectors/gconfig/src/VecGeomChecks.cxx create mode 100644 Detectors/gconfig/src/VecGeomChecks.h create mode 100644 Detectors/gconfig/src/VecGeomG4Map.cxx create mode 100644 Detectors/gconfig/src/VecGeomG4Map.h create mode 100644 Detectors/gconfig/src/VecGeomG4Navigator.cxx create mode 100644 Detectors/gconfig/src/VecGeomG4Navigator.h create mode 100644 Detectors/gconfig/src/VecGeomG4NavigatorBase.cxx create mode 100644 Detectors/gconfig/src/VecGeomG4NavigatorBase.h create mode 100644 Detectors/gconfig/src/VecGeomNavigation.cxx diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 97c3cc4ddb412..63ecf63ab06c6 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -36,8 +36,9 @@ enum class EG4Physics { // enumerating possible geometry navigation modes // (understanding that geometry description is always done with TGeo) enum class EG4Nav { - kTGeo = 0, /* navigate with TGeo */ - kG4 = 1 /* navigate with G4 native geometry */ + kTGeo = 0, /* navigate with TGeo */ + kG4 = 1, /* navigate with G4 native geometry */ + kVecGeom = 2 /* navigate with VecGeom, on the G4 geometry built from TGeo */ }; // parameters to influence the G4 engine @@ -49,6 +50,20 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { EG4Nav navmode = EG4Nav::kTGeo; // geometry navigation mode (default TGeo) + // Settings for navmode == kVecGeom; ignored otherwise. + double vecgeomPushDepth = 1.e-9; // cm; how far past a face, measured across it, a boundary + // point is pushed before it is located + bool vecgeomZeroSafety = false; // answer zero to every safety query; conservative, but it + // shortens steps and so changes the random history + bool vecgeomFlattenAssemblies = true; // dissolve TGeo assemblies into their content when converting + // to VecGeom; the Geant4 touchable keeps the assembly levels + int vecgeomCheckRays = 0; // if > 0, step this many rays out of the interaction point with + // TGeo and VecGeom and report the volumes they enter differently + int vecgeomCheckLocation = 0; // if > 0, locate this many random points with both and report + // the volumes they disagree on + std::string vecgeomCheckVolumes = ""; // comma-separated volumes to cross-check by sampling inside + // their placements + std::string fluenceWeightFile = ""; // file containing the scoring weights (pdg, ekin, weight) std::string const& getPhysicsConfigString() const; diff --git a/Detectors/gconfig/CMakeLists.txt b/Detectors/gconfig/CMakeLists.txt index a1a2b426f2bb4..276a55af61617 100644 --- a/Detectors/gconfig/CMakeLists.txt +++ b/Detectors/gconfig/CMakeLists.txt @@ -14,11 +14,39 @@ o2_add_library(G3Setup PUBLIC_LINK_LIBRARIES MC::Geant3 FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup ) +# Optional VecGeom navigation for Geant4 (G4.navmode=kVecGeom). It needs TGeo2VecGeom and a VecGeom +# with the BVH navigator of the VNavigator family (BVHNavigatorV). Linked PRIVATE: VecGeom types +# never appear in G4Setup's public headers. +find_package(TGeo2VecGeom CONFIG QUIET) +set(G4SETUP_WITH_VECGEOM OFF) +if(TGeo2VecGeom_FOUND) + find_path(O2_VECGEOM_BVHNAVIGATORV_INCLUDE VecGeom/navigation/BVHNavigatorV.h + HINTS ${VecGeom_INCLUDE_DIRS} ${VecGeom_DIR}/../../../include $ENV{VECGEOM_ROOT}/include) + if(O2_VECGEOM_BVHNAVIGATORV_INCLUDE) + set(G4SETUP_WITH_VECGEOM ON) + endif() +endif() + +set(G4SETUP_SOURCES src/G4Config.cxx src/G4RunConfiguration.cxx src/G4LocalFieldConstruction.cxx + src/VecGeomNavigation.cxx) +if(G4SETUP_WITH_VECGEOM) + list(APPEND G4SETUP_SOURCES src/VecGeomG4Map.cxx src/VecGeomChecks.cxx src/VecGeomG4NavigatorBase.cxx + src/VecGeomG4Navigator.cxx) +endif() + o2_add_library(G4Setup - SOURCES src/G4Config.cxx src/G4RunConfiguration.cxx src/G4LocalFieldConstruction.cxx - PUBLIC_LINK_LIBRARIES MC::Geant4VMC MC::Geant4 FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup O2::FastSim + TARGETVARNAME targetG4Setup + SOURCES ${G4SETUP_SOURCES} + PUBLIC_LINK_LIBRARIES MC::Geant4VMC MC::Geant4 FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup O2::FastSim O2::DetectorsBase ) +if(G4SETUP_WITH_VECGEOM) + target_compile_definitions(${targetG4Setup} PRIVATE O2_WITH_VECGEOM) + target_link_libraries(${targetG4Setup} PRIVATE TGeo2VecGeom::TGeo2VecGeom) +else() + message(STATUS "G4.navmode=kVecGeom not built: it needs TGeo2VecGeom and a VecGeom with BVHNavigatorV") +endif() + o2_add_library(FLUKASetup SOURCES src/FlukaConfig.cxx PUBLIC_LINK_LIBRARIES FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup diff --git a/Detectors/gconfig/g4Config.C b/Detectors/gconfig/g4Config.C index 83a932e674e5b..3cc2133961ab4 100644 --- a/Detectors/gconfig/g4Config.C +++ b/Detectors/gconfig/g4Config.C @@ -66,6 +66,7 @@ R__LOAD_LIBRARY(libgeant4vmc) #include "G4VScoringMesh.hh" #include #include "SimSetup/G4RunConfiguration.h" +#include "SimSetup/VecGeomNavigation.h" #endif #include "commonConfig.C" @@ -115,6 +116,10 @@ void Config() geomNavStr = "geomRoot"; } else if (g4Params.navmode == o2::conf::EG4Nav::kG4) { geomNavStr = "geomVMC+RootToGeant4"; + } else if (g4Params.navmode == o2::conf::EG4Nav::kVecGeom) { + // The geometry, its materials and the touchable stay the ones g4root builds from TGeo; + // only the navigator is swapped, once the engine below has built that hierarchy. + geomNavStr = "geomRoot"; } else { LOG(fatal) << "Unsupported geometry navigation mode"; } @@ -137,6 +142,10 @@ void Config() TGeant4* geant4 = new TGeant4("TGeant4", "The Geant4 Monte Carlo", runConfiguration); std::cout << "Geant4 has been created." << std::endl; + if (g4Params.navmode == o2::conf::EG4Nav::kVecGeom) { + o2::simsetup::installVecGeomNavigator(); + } + // setup the stack stackSetup(geant4, FairRunSim::Instance()); diff --git a/Detectors/gconfig/include/SimSetup/VecGeomNavigation.h b/Detectors/gconfig/include/SimSetup/VecGeomNavigation.h new file mode 100644 index 0000000000000..a26a2700f8cb9 --- /dev/null +++ b/Detectors/gconfig/include/SimSetup/VecGeomNavigation.h @@ -0,0 +1,32 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_SIMSETUP_VECGEOMNAVIGATION_H_ +#define O2_SIMSETUP_VECGEOMNAVIGATION_H_ + +namespace o2::simsetup +{ + +/// Whether this build of O2 has the VecGeom navigation backend, i.e. whether TGeo2VecGeom and a +/// VecGeom with BVHNavigatorV were found when O2 was configured. +bool isVecGeomNavigationAvailable(); + +/// Replaces Geant4's tracking navigator by one that answers every navigation query from +/// VecGeom. The Geant4 geometry, its materials and the touchable the scoring code reads stay +/// the ones g4root built from TGeo, so only navigation changes. +/// +/// Call after the TGeant4 engine has been constructed: the Geant4 hierarchy this needs to map +/// onto is built while TG4RunManager configures itself. Aborts if the backend is missing. +void installVecGeomNavigator(); + +} // namespace o2::simsetup + +#endif diff --git a/Detectors/gconfig/src/VecGeomChecks.cxx b/Detectors/gconfig/src/VecGeomChecks.cxx new file mode 100644 index 0000000000000..f17fd12b99c77 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomChecks.cxx @@ -0,0 +1,257 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "VecGeomChecks.h" + +#include "TGeoBBox.h" +#include "TGeoManager.h" +#include "TGeoMatrix.h" +#include "TGeoNode.h" +#include "TGeoVolume.h" +#include "TRandom3.h" + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +using V3 = vecgeom::Vector3D; + +namespace o2::simsetup +{ + +std::size_t checkVecGeomLocation(std::size_t samples) +{ + auto* world = vecgeom::GeoManager::Instance().GetWorld(); + auto const* box = dynamic_cast(gGeoManager->GetTopVolume()->GetShape()); + if (world == nullptr || box == nullptr) { + LOG(warning) << "Cannot cross-check the VecGeom location: no world"; + return 0; + } + TRandom3 rnd(12345); + vecgeom::NavigationState state; + std::map byVolume; + std::size_t bad = 0; + for (std::size_t i = 0; i < samples; ++i) { + const double x = box->GetOrigin()[0] + box->GetDX() * (2. * rnd.Rndm() - 1.); + const double y = box->GetOrigin()[1] + box->GetDY() * (2. * rnd.Rndm() - 1.); + const double z = box->GetOrigin()[2] + box->GetDZ() * (2. * rnd.Rndm() - 1.); + auto* node = gGeoManager->FindNode(x, y, z); + const std::string tgeoName = (node != nullptr) ? node->GetVolume()->GetName() : ""; + state.Clear(); + vecgeom::GlobalLocator::LocateGlobalPoint(world, V3(x, y, z), state, true); + auto const* top = state.Top(); + const std::string vgName = (top != nullptr) ? top->GetLogicalVolume()->GetName() : ""; + if (tgeoName != vgName) { + ++bad; + ++byVolume[tgeoName + " -> " + vgName]; + } + } + LOG(info) << "VecGeom location cross-check: " << bad << " of " << samples << " points land in a " + << "different volume than TGeo puts them in"; + std::vector> worst; + for (auto const& e : byVolume) { + worst.emplace_back(e.second, e.first); + } + std::sort(worst.rbegin(), worst.rend()); + for (std::size_t i = 0; i < worst.size() && i < 15; ++i) { + LOG(info) << " " << worst[i].first << " " << worst[i].second; + } + return bad; +} + +std::size_t checkVecGeomVolume(const char* name, std::size_t perPlacement, std::size_t maxPlacements) +{ + auto* world = vecgeom::GeoManager::Instance().GetWorld(); + auto* vol = gGeoManager->GetVolume(name); + if (world == nullptr || vol == nullptr) { + LOG(warning) << "Cannot cross-check volume " << name << ": not in the geometry"; + return 0; + } + auto const* box = dynamic_cast(vol->GetShape()); + if (box == nullptr) { + LOG(warning) << "Cannot cross-check volume " << name << ": its shape has no bounding box"; + return 0; + } + LOG(info) << "VecGeom sizes: " << vecgeom::GeoManager::Instance().GetRegisteredVolumesCount() + << " logical, " << vecgeom::GeoManager::Instance().GetPlacedVolumesCount() << " placed, " + << vecgeom::VPlacedVolume::GetIdCount() << " ids handed out"; + // Divisions and other generated placements are where the two trees are most likely to differ + // in shape rather than in position, so report the daughter counts before sampling anything. + if (auto* lv = vecgeom::GeoManager::Instance().FindLogicalVolume(name)) { + LOG(info) << " " << name << ": TGeo " << vol->GetNdaughters() << " daughters, VecGeom " + << lv->GetDaughters().size(); + } + { + TIter nextVol(gGeoManager->GetListOfVolumes()); + TGeoVolume* mother = nullptr; + while ((mother = static_cast(nextVol())) != nullptr) { + bool found = false; + for (int i = 0; i < mother->GetNdaughters(); ++i) { + if (mother->GetNode(i)->GetVolume() == vol) { + found = true; + break; + } + } + if (!found) { + continue; + } + auto* mlv = vecgeom::GeoManager::Instance().FindLogicalVolume(mother->GetName()); + LOG(info) << " mother " << mother->GetName() << ": TGeo " << mother->GetNdaughters() + << " daughters, VecGeom " << (mlv != nullptr ? (long)mlv->GetDaughters().size() : -1) + << (mother->GetFinder() != nullptr ? " (divided)" : ""); + break; + } + } + TRandom3 rnd(4321); + vecgeom::NavigationState state; + std::map byResult; + std::size_t placements = 0, tested = 0, bad = 0, badTransform = 0; + double worstTransform = 0.; + TGeoIterator it(gGeoManager->GetTopVolume()); + TGeoNode* node = nullptr; + while ((node = it.Next()) != nullptr && placements < maxPlacements) { + if (node->GetVolume() != vol) { + continue; + } + ++placements; + TGeoHMatrix matrix = *it.GetCurrentMatrix(); + for (std::size_t k = 0; k < perPlacement; ++k) { + double local[3] = {box->GetOrigin()[0] + box->GetDX() * (2. * rnd.Rndm() - 1.), + box->GetOrigin()[1] + box->GetDY() * (2. * rnd.Rndm() - 1.), + box->GetOrigin()[2] + box->GetDZ() * (2. * rnd.Rndm() - 1.)}; + if (!vol->GetShape()->Contains(local)) { + continue; + } + double global[3]; + matrix.LocalToMaster(local, global); + ++tested; + auto* found = gGeoManager->FindNode(global[0], global[1], global[2]); + const std::string tgeoName = (found != nullptr) ? found->GetVolume()->GetName() : ""; + state.Clear(); + vecgeom::GlobalLocator::LocateGlobalPoint(world, V3(global[0], global[1], global[2]), state, true); + auto const* top = state.Top(); + const std::string vgName = (top != nullptr) ? top->GetLogicalVolume()->GetName() : ""; + if (tgeoName != vgName) { + ++bad; + ++byResult[tgeoName + " -> " + vgName]; + continue; + } + // The volume is right; check that the state also composes the right transform. That goes + // through the navigation index table, which is built separately from the daughter lists + // the locators use, so it can be wrong where containment looks perfect. + vecgeom::Transformation3D trans; + state.TopMatrix(trans); + const auto vgLocal = trans.Transform(V3(global[0], global[1], global[2])); + const double d = std::sqrt((vgLocal[0] - local[0]) * (vgLocal[0] - local[0]) + + (vgLocal[1] - local[1]) * (vgLocal[1] - local[1]) + + (vgLocal[2] - local[2]) * (vgLocal[2] - local[2])); + worstTransform = std::max(worstTransform, d); + if (d > 1.e-6) { + ++badTransform; + } + } + } + LOG(info) << "VecGeom volume cross-check " << name << ": worst local-point deviation " + << worstTransform << " cm, " << badTransform << " points above 1e-6 cm"; + LOG(info) << "VecGeom volume cross-check " << name << ": " << placements << " placements, " << tested + << " points, " << bad << " located differently than TGeo"; + for (auto const& e : byResult) { + LOG(info) << " " << e.second << " " << e.first; + } + return bad; +} + +void checkVecGeomRays(std::size_t rays) +{ + auto* world = vecgeom::GeoManager::Instance().GetWorld(); + if (world == nullptr) { + return; + } + TRandom3 rnd(97531); + std::map tgeoSeen, vgSeen; + constexpr std::size_t kMaxSteps = 20000; + + for (std::size_t r = 0; r < rays; ++r) { + const double cost = 2. * rnd.Rndm() - 1.; + const double sint = std::sqrt(1. - cost * cost); + const double phi = 2. * M_PI * rnd.Rndm(); + const double dir[3] = {sint * std::cos(phi), sint * std::sin(phi), cost}; + + gGeoManager->InitTrack(0., 0., 0., dir[0], dir[1], dir[2]); + for (std::size_t k = 0; k < kMaxSteps && !gGeoManager->IsOutside(); ++k) { + ++tgeoSeen[gGeoManager->GetCurrentVolume()->GetName()]; + gGeoManager->FindNextBoundaryAndStep(); + } + + vecgeom::NavigationState cur, next; + V3 pos(0., 0., 0.); + const V3 vdir(dir[0], dir[1], dir[2]); + vecgeom::GlobalLocator::LocateGlobalPoint(world, pos, cur, true); + for (std::size_t k = 0; k < kMaxSteps && cur.Top() != nullptr; ++k) { + ++vgSeen[cur.Top()->GetLogicalVolume()->GetName()]; + double safety = 0.; + auto const* nav = cur.Top()->GetLogicalVolume()->GetNavigator(); + const double step = nav->ComputeStepAndSafetyAndPropagatedState(pos, vdir, vecgeom::kInfLength, cur, next, + false, safety); + if (!(step < vecgeom::kInfLength)) { + break; + } + pos = pos + step * vdir; + cur = next; + } + } + + // A volume one engine never enters is the sharpest signal: a hit can only be made in a + // volume a track actually reaches, so these are the ones that lose a detector its hits. + std::vector neverVG, neverTGeo; + std::vector> diff; + for (auto const& e : tgeoSeen) { + const std::size_t vg = vgSeen.count(e.first) ? vgSeen[e.first] : 0; + if (vg == 0) { + neverVG.push_back(e.first); + } + const long d = static_cast(e.second) - static_cast(vg); + if (d != 0) { + diff.emplace_back(std::labs(d), e.first); + } + } + for (auto const& e : vgSeen) { + if (tgeoSeen.count(e.first) == 0) { + neverTGeo.push_back(e.first); + } + } + std::sort(diff.rbegin(), diff.rend()); + LOG(info) << "VecGeom ray cross-check over " << rays << " rays: " << tgeoSeen.size() << " volumes seen by TGeo, " + << vgSeen.size() << " by VecGeom, " << diff.size() << " entered a different number of times"; + LOG(info) << " never entered by VecGeom (" << neverVG.size() << "):"; + for (std::size_t i = 0; i < neverVG.size() && i < 60; ++i) { + LOG(info) << " " << neverVG[i] << " (TGeo " << tgeoSeen[neverVG[i]] << ")"; + } + LOG(info) << " never entered by TGeo (" << neverTGeo.size() << "):"; + for (std::size_t i = 0; i < neverTGeo.size() && i < 30; ++i) { + LOG(info) << " " << neverTGeo[i] << " (VecGeom " << vgSeen[neverTGeo[i]] << ")"; + } +} + +} // namespace o2::simsetup diff --git a/Detectors/gconfig/src/VecGeomChecks.h b/Detectors/gconfig/src/VecGeomChecks.h new file mode 100644 index 0000000000000..2f1efa9922928 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomChecks.h @@ -0,0 +1,37 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_SIMSETUP_VECGEOMCHECKS_H_ +#define O2_SIMSETUP_VECGEOMCHECKS_H_ + +#include + +namespace o2::simsetup +{ + +/// Samples random points in the world and compares the volume TGeo locates them in with the +/// volume VecGeom locates them in, which tests the conversion and the level locators on their +/// own, with no Geant4 and no stepping involved. Returns the number of disagreements. +std::size_t checkVecGeomLocation(std::size_t samples); + +/// The same comparison, but sampling inside the placements of one named volume rather than over +/// the world. A thin sensitive volume is never sampled often enough by a scan over the world, +/// so this is what tells you whether such a volume is located correctly. +std::size_t checkVecGeomVolume(const char* name, std::size_t perPlacement, std::size_t maxPlacements); + +/// Shoots rays from the interaction point and steps them to the world edge with TGeo and with +/// VecGeom, comparing how often each engine reports being in each volume. Stepping is what a +/// containment scan cannot test, and a volume a ray never enters is invisible to a hit count. +void checkVecGeomRays(std::size_t rays); + +} // namespace o2::simsetup + +#endif diff --git a/Detectors/gconfig/src/VecGeomG4Map.cxx b/Detectors/gconfig/src/VecGeomG4Map.cxx new file mode 100644 index 0000000000000..7285e45ce5fec --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4Map.cxx @@ -0,0 +1,187 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "VecGeomG4Map.h" + +#include "TG4RootDetectorConstruction.h" + +#include "G4PhysicalVolumeStore.hh" +#include "G4VPhysicalVolume.hh" + +#include "TGeoManager.h" +#include "TGeoNode.h" +#include "TGeoVolume.h" +#include "TString.h" + +#include "TGeo2VecGeom/RootGeoManager.h" +#include +#include + +#include + +#include +#include +#include + +namespace o2::simsetup +{ + +namespace +{ +/// Reproduces tgeo2vecgeom's assembly flattening for the daughters of one TGeo volume, so that +/// each VecGeom daughter can be paired with the chain of TGeo nodes it stands for. The +/// converter walks a volume's daughters in order and, when flattening, descends into an +/// assembly instead of placing it, appending the leaves it reaches; the placed daughters of +/// the VecGeom logical volume come out in exactly that order. +void flattenDaughter(TGeoNode* node, bool flatten, std::vector& chain, + std::vector>& out) +{ + chain.push_back(node); + auto* assembly = dynamic_cast(node->GetVolume()); + if (flatten && assembly != nullptr) { + for (int i = 0; i < assembly->GetNdaughters(); ++i) { + flattenDaughter(assembly->GetNode(i), flatten, chain, out); + } + } else { + out.push_back(chain); + } + chain.pop_back(); +} + +std::vector> flattenedDaughters(TGeoVolume const* volume, bool flatten) +{ + std::vector> out; + std::vector chain; + for (int i = 0; i < volume->GetNdaughters(); ++i) { + flattenDaughter(volume->GetNode(i), flatten, chain, out); + } + return out; +} +} // namespace + +void VecGeomG4Map::registerPair(vecgeom::VPlacedVolume const* parent, vecgeom::VPlacedVolume const* pv, + std::vector const& nodes, TG4RootDetectorConstruction const& dc) +{ + const auto vgId = static_cast(pv->id()); + if (vgId >= mChainBegin.size()) { + mChainBegin.resize(vgId + 1, 0); + mChainSize.resize(vgId + 1, 0); + } + if (mChainSize[vgId] != 0) { // a logical volume placed more than once shares its daughters + return; + } + mChainBegin[vgId] = static_cast(mChain.size()); + bool first = true; + mChainSize[vgId] = static_cast(nodes.size()); + for (auto* node : nodes) { + auto* g4pv = dc.GetG4VPhysicalVolume(node); + if (g4pv == nullptr) { + LOG(fatal) << "TGeo node " << node->GetName() << " has no Geant4 counterpart; the two conversions " + << "of the geometry do not agree"; + } + mChain.push_back(g4pv); + if (first && parent == nullptr) { + mWorld = pv; + } + first = false; + if (node != nodes.back()) { + continue; // an intermediate assembly level, dissolved on the VecGeom side + } + const auto g4Id = static_cast(g4pv->GetInstanceID()); + if (g4Id >= mG4ToVG.size()) { + mG4ToVG.resize(g4Id + 1, nullptr); + } + if (mG4ToVG[g4Id] == nullptr) { + mG4ToVG[g4Id] = pv; + } else if (mG4ToVG[g4Id] == ambiguous()) { + mAmbiguous[static_cast(g4Id)].push_back(pv); + } else if (mG4ToVG[g4Id] != pv) { + // Reached through more than one flattened chain: the content of an assembly placed in + // more than one place. They are told apart by the chain itself, not by this volume. + auto& list = mAmbiguous[static_cast(g4Id)]; + list.push_back(mG4ToVG[g4Id]); + list.push_back(pv); + mG4ToVG[g4Id] = ambiguous(); + } + } + ++mPairs; +} + +void VecGeomG4Map::build(TG4RootDetectorConstruction const& dc, bool flattenAssemblies) +{ + auto& rootGeoMgr = tgeo2vecgeom::RootGeoManager::Instance(); + auto& vgMgr = vecgeom::GeoManager::Instance(); + // Ids are handed out by a global counter, so the largest one can exceed the number of + // volumes the manager holds; sizing by the count alone reads past the end. + const auto reserve = std::max(vgMgr.GetPlacedVolumesCount(), vecgeom::VPlacedVolume::GetIdCount()) + 1; + mChainBegin.assign(reserve, 0); + mChainSize.assign(reserve, 0); + mChain.reserve(reserve); + mG4ToVG.assign(G4PhysicalVolumeStore::GetInstance()->size() + 1, nullptr); + + auto* topNode = gGeoManager->GetTopNode(); + auto const* topPV = rootGeoMgr.Lookup(topNode); + if (topPV == nullptr) { + LOG(fatal) << "The VecGeom geometry has no counterpart for the TGeo top node"; + } + registerPair(nullptr, topPV, {topNode}, dc); + + // A TGeoVolume placed twice shares one set of daughter nodes, and so do both conversions of + // it, so each volume's daughters are paired up exactly once. + std::unordered_set seen; + std::function walk = + [&](vecgeom::VPlacedVolume const* pv, TGeoNode* node) { + auto* volume = node->GetVolume(); + if (!seen.insert(volume).second) { + return; + } + const auto chains = flattenedDaughters(volume, flattenAssemblies); + auto const& vgDaughters = pv->GetLogicalVolume()->GetDaughters(); + if (vgDaughters.size() != chains.size()) { + LOG(fatal) << "Volume " << volume->GetName() << " has " << chains.size() + << " daughters after flattening but its VecGeom counterpart has " << vgDaughters.size() + << "; the flattening reproduced here does not match the converter's"; + } + for (std::size_t i = 0; i < chains.size(); ++i) { + auto const* daughter = vgDaughters[i]; + // The pairing is by position, so check it against what the converter recorded. The + // converter names a node it synthesised while flattening after the original one, so + // the recorded node must be either that node itself or a flattened copy of it. + auto const* recorded = rootGeoMgr.tgeonode(daughter); + auto const* expected = chains[i].back(); + const bool ok = recorded != nullptr && + (recorded == expected || + (recorded->GetVolume() == expected->GetVolume() && + TString(recorded->GetName()).BeginsWith(TString(expected->GetName()) + "_assemblyinternalcount_"))); + if (!ok) { + ++mMispaired; + if (mMispaired <= 5) { + LOG(warning) << "VecGeom daughter " << i << " of " << volume->GetName() << " is " + << (recorded != nullptr ? recorded->GetName() : "unknown") << ", expected " + << expected->GetName(); + } + } + registerPair(pv, daughter, chains[i], dc); + walk(daughter, chains[i].back()); + } + }; + walk(topPV, topNode); + + if (mMispaired != 0) { + LOG(error) << mMispaired << " VecGeom placements were paired with the wrong TGeo node; the " + << "order the converter places flattened daughters in is not the order assumed here"; + } + LOG(info) << "VecGeom navigation: paired " << mPairs << " placements with " << mChain.size() + << " Geant4 volumes (flattenAssemblies=" << flattenAssemblies << "), " << ambiguousCount() + << " of them shared by more than one placement"; +} + +} // namespace o2::simsetup diff --git a/Detectors/gconfig/src/VecGeomG4Map.h b/Detectors/gconfig/src/VecGeomG4Map.h new file mode 100644 index 0000000000000..298d3c2a223cd --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4Map.h @@ -0,0 +1,124 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_SIMSETUP_VECGEOMG4MAP_H_ +#define O2_SIMSETUP_VECGEOMG4MAP_H_ + +#include + +#include +#include +#include +#include + +class G4VPhysicalVolume; +class TG4RootDetectorConstruction; +class TGeoNode; + +namespace o2::simsetup +{ + +/// The correspondence between the VecGeom tree tgeo2vecgeom built and the Geant4 tree g4root +/// built. Both are conversions of the same TGeoNode hierarchy, but not necessarily level for +/// level: with assembly flattening one VecGeom placement stands for a whole assembly subtree, +/// so a placement maps to a *chain* of Geant4 physical volumes rather than to one. Pushing the +/// whole chain is what keeps the Geant4 touchable at the depth the TGeo navigator produces, +/// and hence keeps every volume id, copy number and CurrentVolOffName offset unchanged. +/// +/// Without flattening every chain has length one. Lookups are vector indexing: VecGeom +/// placed-volume ids and Geant4 physical-volume instance ids are both dense. +class VecGeomG4Map +{ + public: + /// Walks the TGeo hierarchy and pairs up the two conversions of it. Both must exist. + /// \param flattenAssemblies must match what the converter was told, since it decides how + /// many TGeo nodes stand behind one VecGeom placement. + void build(TG4RootDetectorConstruction const& dc, bool flattenAssemblies); + + /// The Geant4 volumes a VecGeom placement stands for, outermost first. A size of zero means + /// the placement was never paired up, which is a bug rather than a legal state. + G4VPhysicalVolume* const* chain(int vgId, unsigned& size) const + { + if (static_cast(vgId) >= mChainSize.size()) { + size = 0; + return nullptr; + } + size = mChainSize[vgId]; + return mChain.data() + mChainBegin[vgId]; + } + + /// The VecGeom placement whose chain ends at the Geant4 volume \a g4Id. Null means that + /// volume is only ever an intermediate assembly level, which flattening dissolved on the + /// VecGeom side, so the walk extends past it. A plain vector index; the few Geant4 volumes + /// reachable through more than one flattened chain answer ambiguous() and go to candidates(). + vecgeom::VPlacedVolume const* toVecGeom(int g4Id) const { return mG4ToVG[g4Id]; } + + /// Marks a Geant4 volume that several VecGeom placements reach. + static vecgeom::VPlacedVolume const* ambiguous() + { + return reinterpret_cast(std::uintptr_t{1}); + } + + /// The placements that share an ambiguous Geant4 volume, to be told apart by their chains. + std::vector const& candidates(int g4Id) const + { + static const std::vector empty; + const auto it = mAmbiguous.find(g4Id); + return it == mAmbiguous.end() ? empty : it->second; + } + + /// Whether a placement's chain is exactly the given run of Geant4 volumes. + bool chainMatches(vecgeom::VPlacedVolume const* pv, G4VPhysicalVolume* const* first, std::size_t n) const + { + unsigned size = 0; + auto* const* own = chain(pv->id(), size); + if (size != n) { + return false; + } + for (std::size_t i = 0; i < n; ++i) { + if (own[i] != first[i]) { + return false; + } + } + return true; + } + + /// How many Geant4 levels a VecGeom placement accounts for. + unsigned chainSize(int vgId) const + { + return static_cast(vgId) < mChainSize.size() ? mChainSize[vgId] : 0; + } + + /// The outermost placement, standing for the Geant4 world volume. + vecgeom::VPlacedVolume const* world() const { return mWorld; } + + /// How many Geant4 volumes share more than one VecGeom placement. + std::size_t ambiguousCount() const { return mAmbiguous.size(); } + + std::size_t size() const { return mPairs; } + + private: + void registerPair(vecgeom::VPlacedVolume const* parent, vecgeom::VPlacedVolume const* pv, + std::vector const& nodes, TG4RootDetectorConstruction const& dc); + + std::vector mChainBegin; + std::vector mChainSize; + std::vector mChain; + std::vector mG4ToVG; + std::unordered_map> mAmbiguous; + vecgeom::VPlacedVolume const* mWorld = nullptr; + std::size_t mPairs = 0; + std::size_t mMispaired = 0; +}; + +} // namespace o2::simsetup + +#endif diff --git a/Detectors/gconfig/src/VecGeomG4Navigator.cxx b/Detectors/gconfig/src/VecGeomG4Navigator.cxx new file mode 100644 index 0000000000000..66418efe7de64 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4Navigator.cxx @@ -0,0 +1,496 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "VecGeomG4Navigator.h" + +#include "G4Exception.hh" +#include "G4TouchableHistory.hh" +#include "G4VPhysicalVolume.hh" + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +namespace +{ +const G4ThreeVector kNoPoint(-1e8, -1e8, -1e8); + +bool samePoint(const G4ThreeVector& a, const G4ThreeVector& b) { return a.diff2(b) < 1e-20; } + +vecgeom::VPlacedVolume const* topOf(vecgeom::NavigationState const& st) +{ + return st.IsOutside() ? nullptr : st.Top(); +} +} // namespace + +namespace o2::simsetup +{ + +VecGeomG4Navigator::VecGeomG4Navigator(VecGeomG4Map const& map, double pushDepth, bool zeroSafety) + : VecGeomG4NavigatorBase(map), mPushDepth(pushDepth), mZeroSafety(zeroSafety) +{ + mEmptyState.Clear(); +} + +VecGeomG4Navigator::~VecGeomG4Navigator() +{ + LOG(info) << "VecGeom navigation: zero steps " << mZeroStepCount << ", stuck pushes " << mStuckPushCount + << ", abandoned " << mAbandonCount << ", negative safeties " << mNegativeSafetyCount + << ", unmappable touchables " << mUnmappableHistoryCount << ", relocated resumes " + << mRelocatedResumeCount << ", missing exit normals " << mNoNormalCount; +} + +G4VPhysicalVolume* VecGeomG4Navigator::updateG4History() +{ + // The history is a function of the VecGeom state alone, so it is rebuilt only when that changed. + if (mHistoryValid && mCurState.HasSamePathAsOther(mHistoryState)) { + return fHistory.GetTopVolume(); + } + mHistoryState = mCurState; + mHistoryValid = true; + return historyFromState(mCurState); +} + +void VecGeomG4Navigator::locateFromWorld(const V3& point) +{ + mCurState.Clear(); + vecgeom::GlobalLocator::LocateGlobalPoint(vecgeom::GeoManager::Instance().GetWorld(), point, mCurState, true); +} + +/// The push is set across the face, as a depth, so that a grazing track is moved off the face as +/// surely as one at normal incidence, and it is kept small: TOF has layers 2.4e-8 cm apart, and a +/// fixed push along the direction steps over them. It never goes below a thousand times the rounding +/// of the largest coordinate, which matters far from the origin. +double VecGeomG4Navigator::boundaryPush(const V3& point, const V3& dir) const +{ + constexpr double kMaxPush = 1.e-4; // cm along the direction + const double big = std::max({std::abs(point.x()), std::abs(point.y()), std::abs(point.z())}); + const double rounding = 1.e3 * big * std::numeric_limits::epsilon(); + + // The face just crossed: the entered daughter's, or the current volume's own. + vecgeom::NavigationState const& st = mWouldEnter ? mNextState : mCurState; + double cosn = 1.; + if (!st.IsOutside() && st.Top() != nullptr) { + vecgeom::Transformation3D m; + st.TopMatrix(m); + V3 n; + st.Top()->GetUnplacedVolume()->Normal(m.Transform(point), n); + const double c = std::abs(n.Dot(m.TransformDirection(dir))); + if (c > 0. && n.Mag2() > 0.5) { + cosn = c; + } + } + return std::max(rounding, std::min(mPushDepth / cosn, kMaxPush)); +} + +G4double VecGeomG4Navigator::ComputeStep(const G4ThreeVector& globalPoint, const G4ThreeVector& direction, + const G4double proposedStepLength, G4double& newSafety) +{ + newSafety = 0.; + mWouldEnter = false; + mWouldExit = false; + + auto const* top = topOf(mCurState); + if (top == nullptr) { // the track is outside the world + return kInfinity; + } + auto const* navigator = top->GetLogicalVolume()->GetNavigator(); + + // On the point a boundary locate left the track on, the safety is zero; a point seen before + // reuses its safety. Otherwise the navigator computes it with the step. + bool calcSafety = !mZeroSafety && !(mLocatedOnBoundary && samePoint(globalPoint, mLastLocatedPoint)); + if (calcSafety && samePoint(globalPoint, mSafetyOrig)) { + calcSafety = false; + newSafety = mLastSafety; + } + + // The step is computed on a copy. The state, and the volume blocked by the last exit, stay as + // the locate left them for every call until the next locate, which is what the field propagator + // relies on when it calls this from trial points along the curve. + mStepState = mCurState; + const double limit = std::min(proposedStepLength * kG4ToVG, static_cast(vecgeom::kInfLength)); + double safety = 0.; + double vgStep = navigator->ComputeStepAndSafety(toVG(globalPoint), toDir(direction), limit, mStepState, calcSafety, + safety, true); + if (calcSafety) { + if (safety < 0.) { + ++mNegativeSafetyCount; + safety = 0.; + } + newSafety = safety * kVGToG4; + mSafetyOrig = globalPoint; + mLastSafety = newSafety; + } + const bool boundaryLimited = vgStep < limit; + G4double step = std::max(vgStep, 0.) * kVGToG4; + const bool entering = mStepState.GetCurrentLevel() > mCurState.GetCurrentLevel(); + + // A track that is not moving, handled as G4Navigator does: after ten zero steps the step is + // lengthened by 100 kCarTolerance, after twenty-five the event is aborted. + if (step < 0.05 * kCarTolerance) { + ++mZeroStepCount; + if (++mNzeroSteps >= kActionThresholdNoZeroSteps) { + ++mStuckPushCount; + step += 100. * kCarTolerance; + if (mNzeroSteps >= kAbandonThresholdNoZeroSteps) { + ++mAbandonCount; + std::ostringstream msg; + msg << "Track stuck or not moving: " << mNzeroSteps << " zero steps in " << top->GetLabel() << " at (" + << globalPoint.x() << ", " << globalPoint.y() << ", " << globalPoint.z() << ") mm. Event aborted, as " + << "G4Navigator does."; + mNzeroSteps = 0; + G4Exception("VecGeomG4Navigator::ComputeStep()", "GeomNav0003", EventMustBeAborted, msg.str().c_str()); + } + } + } else { + mNzeroSteps = 0; + } + + if (boundaryLimited) { + mWouldEnter = entering; + mWouldExit = !entering; + mNextPoint = globalPoint + step * direction; + if (entering) { + mNextState = mStepState; + } + // The surface this step ends on, for the exit normal: the entered daughter's or our own. + mNormalState = entering ? mNextState : mCurState; + mNormalEnter = entering; + mNormalPoint = mNextPoint; + mNormalValid = true; + } else { + step = kInfinity; + mNormalValid = false; + } + return step; +} + +/// Geant4's rule for a point on several coincident faces: go up past every volume the point is on +/// the surface of and heading out of, then look down again from there, never back into the volume +/// just left, which stays blocked for the next step. Bounded, because popping and descending can +/// meet another flush face. It never pops to or above minLevel: a daughter the step entered stays +/// entered, because the step decided that with the chord direction, and the direction a locate gets +/// under a field is the momentum. Nor does it descend back into `avoid`, the volume the crossing +/// exited: on a face a helix touches tangentially chord and momentum disagree about the side. +void VecGeomG4Navigator::leaveFlushVolumes(const V3& point, const V3& dir, int minLevel, + vecgeom::VPlacedVolume const* avoid) +{ + for (int round = 0; round < 4; ++round) { + vecgeom::VPlacedVolume const* left = nullptr; + while (topOf(mCurState) != nullptr && static_cast(mCurState.GetCurrentLevel()) > minLevel) { + vecgeom::Transformation3D m; + mCurState.TopMatrix(m); + if (mCurState.Top()->GetUnplacedVolume()->DistanceToOut(m.Transform(point), m.TransformDirection(dir)) > 0.) { + break; + } + left = mCurState.Top(); + mCurState.SetLastExited(); + if (mCurState.GetCurrentLevel() <= 1) { + mCurState.Clear(); // nothing to travel in even in the world: the track left it + return; + } + mCurState.Pop(); + } + if (left == nullptr) { + return; + } + auto const* mother = mCurState.Top(); + vecgeom::Transformation3D m; + mCurState.TopMatrix(m); + const auto level = mCurState.GetCurrentLevel(); + const auto blocked = mCurState.GetLastExitedState(); + mCurState.Pop(); + vecgeom::GlobalLocator::LocateGlobalPointExclVolume(mother, left, m.Transform(point), mCurState, false); + mCurState.SetLastExited(blocked); + if (mCurState.GetCurrentLevel() == level) { + return; // no daughter holds the point + } + if (avoid != nullptr && mCurState.GetCurrentLevel() > level) { + // Undo a descent into the volume the crossing exited, at whatever depth the search put it. + mPathScratch = mCurState; + while (mPathScratch.GetCurrentLevel() > level && mPathScratch.Top() != avoid) { + mPathScratch.Pop(); + } + if (mPathScratch.GetCurrentLevel() > level) { + while (mCurState.GetCurrentLevel() > level) { + mCurState.Pop(); + } + return; + } + } + } +} + +/// Geant4's two flags describe the transition, not a change of depth: a track leaving one volume +/// straight into a touching sibling has both exited and entered. Everything below the common prefix +/// of the old path (mReloScratch) was left, everything below it on the new path entered. +void VecGeomG4Navigator::updateCrossingFlags(bool entering) +{ + const int preLevel = mReloScratch.IsOutside() ? 0 : static_cast(mReloScratch.GetCurrentLevel()); + const int curLevel = mCurState.IsOutside() ? 0 : static_cast(mCurState.GetCurrentLevel()); + mPathScratch = mReloScratch; + mStepState = mCurState; + int la = preLevel, lb = curLevel; + while (la > lb) { + mPathScratch.Pop(); + --la; + } + while (lb > la) { + mStepState.Pop(); + --lb; + } + while (la > 0 && !mPathScratch.HasSamePathAsOther(mStepState)) { + mPathScratch.Pop(); + mStepState.Pop(); + --la; + } + fExitedMother = !entering || la < preLevel; + fEnteredDaughter = entering || la < curLevel; +} + +G4VPhysicalVolume* VecGeomG4Navigator::LocateGlobalPointAndSetup(const G4ThreeVector& point, + const G4ThreeVector* direction, + const G4bool relativeSearch, const G4bool) +{ + // A boundary is being crossed when Geant4 says the last step was limited by the geometry (or the + // point is where the last ComputeStep put the boundary) and that step found one. + const bool onBoundary = relativeSearch && (fWasLimitedByGeometry || samePoint(point, mNextPoint)); + const bool crossing = onBoundary && (mWouldEnter || mWouldExit) && direction != nullptr; + const bool entering = crossing && mWouldEnter; + const V3 p = toVG(point); + const V3 dir = direction != nullptr ? toDir(*direction) : V3(0., 0., 0.); + const double push = (onBoundary && direction != nullptr) ? boundaryPush(p, dir) : 0.; + const V3 q = p + push * dir; + + fWasLimitedByGeometry = false; + fEnteredDaughter = false; + fExitedMother = false; + mLocatedOnBoundary = false; + mLastLocatedPoint = point; + mSafetyOrig = kNoPoint; + clearLastExited(); + + if (!relativeSearch) { + mNzeroSteps = 0; // a new track: nothing of the previous one applies + locateFromWorld(p); + } else if (topOf(mCurState) == nullptr) { + locateFromWorld(p); + } else if (entering) { + // Into the daughter the step hit, then down inside it. + mReloScratch = mCurState; + mCurState = mNextState; + auto const* daughter = mCurState.Top(); + mCurState.Pop(); + vecgeom::Transformation3D m; + mCurState.TopMatrix(m); + vecgeom::GlobalLocator::LocateGlobalPoint(daughter, daughter->GetTransformation()->Transform(m.Transform(q)), + mCurState, false); + mLocatedOnBoundary = true; + } else if (crossing) { + // Out of the current volume: up until the point is contained, then down, never back into the + // volume just left; that volume stays blocked at zero distance for the next step, as + // G4Navigator's fBlockedPhysicalVolume. + mReloScratch = mCurState; + if (mCurState.GetCurrentLevel() <= 1) { + mCurState.Clear(); // left the world + } else { + vecgeom::Transformation3D m; + mCurState.TopMatrix(m); + vecgeom::GlobalLocator::RelocatePointFromPathForceDifferent(m.Transform(q), mCurState); + mReloScratch.SetLastExited(); + mCurState.SetLastExited(mReloScratch.GetLastExitedState()); + } + mLocatedOnBoundary = true; + } else { + // Anywhere else: from the current path, up until contained, then down. + vecgeom::Transformation3D m; + mCurState.TopMatrix(m); + vecgeom::GlobalLocator::RelocatePointFromPath(m.Transform(q), mCurState); + if (topOf(mCurState) == nullptr) { + locateFromWorld(p); + } + mLocatedOnBoundary = onBoundary; + } + if (crossing) { + leaveFlushVolumes(q, dir, entering ? static_cast(mNextState.GetCurrentLevel()) : 0, + entering ? nullptr : topOf(mReloScratch)); + updateCrossingFlags(entering); + } + mWouldEnter = false; + mWouldExit = false; + return updateG4History(); +} + +G4VPhysicalVolume* VecGeomG4Navigator::ResetHierarchyAndLocate(const G4ThreeVector& point, const G4ThreeVector&, + const G4TouchableHistory& history) +{ + // A track resumes from a stored touchable, usually a secondary starting where its parent's step + // ended. The touchable names its volume; it is kept unless it does not hold the point. + fWasLimitedByGeometry = false; + fEnteredDaughter = false; + fExitedMother = false; + mWouldEnter = false; + mWouldExit = false; + mLocatedOnBoundary = false; + mLastLocatedPoint = point; + mSafetyOrig = kNoPoint; + mNzeroSteps = 0; + fHistory = *history.GetHistory(); + const V3 p = toVG(point); + bool kept = false; + if (stateFromHistory(mCurState)) { + vecgeom::Transformation3D m; + mCurState.TopMatrix(m); + if (mCurState.Top()->GetUnplacedVolume()->Contains(m.Transform(p))) { + kept = true; + } else { + ++mRelocatedResumeCount; + vecgeom::GlobalLocator::RelocatePointFromPath(m.Transform(p), mCurState); + if (topOf(mCurState) == nullptr) { + locateFromWorld(p); + } + } + } else { + if (++mUnmappableHistoryCount <= 10) { + LOG(warning) << "VecGeom navigation: a touchable matches no VecGeom path; locating from the world"; + } + locateFromWorld(p); + } + clearLastExited(); + if (kept) { + // The touchable's history is the path of the state: keep it rather than rebuild it. + mHistoryState = mCurState; + mHistoryValid = true; + return fHistory.GetTopVolume(); + } + mHistoryValid = false; + return updateG4History(); +} + +void VecGeomG4Navigator::LocateGlobalPointWithinVolume(const G4ThreeVector& position) +{ + // The caller guarantees the point is in the current volume, so the state stays; only what + // described the last crossing is dropped, as in G4Navigator. + mLastLocatedPoint = position; + mLocatedOnBoundary = false; + mWouldEnter = false; + mWouldExit = false; + fEnteredDaughter = false; + fExitedMother = false; + clearLastExited(); +} + +G4double VecGeomG4Navigator::ComputeSafety(const G4ThreeVector& globalPoint, const G4double, const G4bool) +{ + if (mZeroSafety) { + return 0.; + } + if (mLocatedOnBoundary && samePoint(globalPoint, mLastLocatedPoint)) { + return 0.; + } + if ((mWouldEnter || mWouldExit) && samePoint(globalPoint, mNextPoint)) { + return 0.; + } + if (samePoint(globalPoint, mSafetyOrig)) { + return mLastSafety; + } + auto const* top = topOf(mCurState); + if (top == nullptr) { + return 0.; + } + auto const* estimator = top->GetLogicalVolume()->GetSafetyEstimator(); + if (estimator == nullptr) { + return 0.; + } + double safety = estimator->ComputeSafety(toVG(globalPoint), mCurState); + if (safety < 0.) { + ++mNegativeSafetyCount; + safety = 0.; + } + mSafetyOrig = globalPoint; + mLastSafety = safety * kVGToG4; + return mLastSafety; +} + +/// The normal of the surface the last ComputeStep ended on, in global coordinates, with Geant4's +/// convention: out of the volume being left, or into the daughter being entered. The field +/// propagator's intersection locator compares it with the momentum to see whether a curved track +/// turns back through the surface, so the sign is the geometric one, not the direction of travel. +/// VecGeom's Normal() also reports whether the point was on the surface; Geant4 hands back points a +/// few nanometres off the face, so only a degenerate vector is refused. +bool VecGeomG4Navigator::computeExitNormal(const G4ThreeVector& point, V3& globalNormal) const +{ + if (!mNormalValid || topOf(mNormalState) == nullptr) { + return false; + } + vecgeom::Transformation3D m; + mNormalState.TopMatrix(m); + V3 localNormal(0., 0., 0.); + mNormalState.Top()->GetUnplacedVolume()->Normal(m.Transform(toVG(point)), localNormal); + globalNormal = m.InverseTransformDirection(localNormal); + const double mag = globalNormal.Mag(); + if (!(mag > 0.5)) { + return false; + } + globalNormal /= mag; + if (mNormalEnter) { + globalNormal = -globalNormal; + } + return true; +} + +G4ThreeVector VecGeomG4Navigator::GetLocalExitNormal(G4bool* valid) +{ + V3 n; + if (!computeExitNormal(mNormalPoint, n) || topOf(mCurState) == nullptr) { + ++mNoNormalCount; + if (valid != nullptr) { + *valid = false; + } + return G4ThreeVector(); + } + vecgeom::Transformation3D m; + mCurState.TopMatrix(m); + const auto local = m.TransformDirection(n); + if (valid != nullptr) { + *valid = true; + } + return G4ThreeVector(local[0], local[1], local[2]); +} + +G4ThreeVector VecGeomG4Navigator::GetGlobalExitNormal(const G4ThreeVector& point, G4bool* valid) +{ + V3 n; + if (!computeExitNormal(point, n)) { + ++mNoNormalCount; + if (valid != nullptr) { + *valid = false; + } + return G4ThreeVector(); + } + if (valid != nullptr) { + *valid = true; + } + return G4ThreeVector(n[0], n[1], n[2]); +} + +} // namespace o2::simsetup diff --git a/Detectors/gconfig/src/VecGeomG4Navigator.h b/Detectors/gconfig/src/VecGeomG4Navigator.h new file mode 100644 index 0000000000000..54a3325672ec8 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4Navigator.h @@ -0,0 +1,116 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_SIMSETUP_VECGEOMG4NAVIGATOR_H_ +#define O2_SIMSETUP_VECGEOMG4NAVIGATOR_H_ + +#include "VecGeomG4NavigatorBase.h" + +namespace o2::simsetup +{ + +/// A Geant4 tracking navigator that answers every navigation query from VecGeom while keeping the +/// Geant4 navigation history, the touchable the scoring code reads, in step with the VecGeom state. +/// It works as G4VecGeomNav's +/// TG4VecGeomNavigator does: +/// +/// - ComputeStep leaves the current volume alone. It records whether the step ends on a boundary +/// and, for a daughter hit, the state that enters it. +/// - The locate on a boundary relocates: into the recorded daughter, or out of the current volume, +/// which then stays blocked for the next step. The point is first pushed across the face by a +/// small depth, and afterwards leaves every volume it is flush with and heading out of. +/// - Safety is zero only at the boundary point itself. +/// +/// Unlike TG4VecGeomNavigator, the VecGeom geometry is converted from TGeo, not from Geant4, so one +/// VecGeom placement can stand for a chain of g4root volumes (VecGeomG4Map). +class VecGeomG4Navigator : public VecGeomG4NavigatorBase +{ + public: + /// \param pushDepth how far past a face (cm, measured across it) a boundary point is pushed before + /// it is located. \param zeroSafety answers zero to every safety query. + VecGeomG4Navigator(VecGeomG4Map const& map, double pushDepth, bool zeroSafety); + ~VecGeomG4Navigator() override; + + G4double ComputeStep(const G4ThreeVector& globalPoint, const G4ThreeVector& direction, + const G4double proposedStepLength, G4double& newSafety) override; + + G4VPhysicalVolume* ResetHierarchyAndLocate(const G4ThreeVector& point, const G4ThreeVector& direction, + const G4TouchableHistory& history) override; + + G4VPhysicalVolume* LocateGlobalPointAndSetup(const G4ThreeVector& point, const G4ThreeVector* direction = nullptr, + const G4bool relativeSearch = true, + const G4bool ignoreDirection = true) override; + + void LocateGlobalPointWithinVolume(const G4ThreeVector& position) override; + + G4double ComputeSafety(const G4ThreeVector& globalPoint, const G4double proposedMaxLength = DBL_MAX, + const G4bool keepState = true) override; + + // Both point out of the volume left; the local one is in the frame of the final volume. + G4ThreeVector GetLocalExitNormal(G4bool* valid) override; + G4ThreeVector GetGlobalExitNormal(const G4ThreeVector& point, G4bool* valid) override; + + private: + /// Rewrites the Geant4 history from mCurState, unless it already stands for it. + G4VPhysicalVolume* updateG4History(); + void locateFromWorld(const V3& point); + /// Go up past every volume the point is on the surface of and heading out of, then down again. + void leaveFlushVolumes(const V3& point, const V3& dir, int minLevel, vecgeom::VPlacedVolume const* avoid); + /// Sets fEnteredDaughter and fExitedMother from the paths before (mReloScratch) and after a crossing. + void updateCrossingFlags(bool entering); + void clearLastExited() { mCurState.SetLastExited(mEmptyState.GetLastExitedState()); } + /// How far (cm) a boundary point is pushed along the direction before it is located. + double boundaryPush(const V3& point, const V3& dir) const; + /// The normal of the surface the last geometry-limited ComputeStep ended on, global, unit length. + bool computeExitNormal(const G4ThreeVector& point, V3& globalNormal) const; + + double mPushDepth = 1.e-9; ///< cm + bool mZeroSafety = false; + + vecgeom::NavigationState mCurState; ///< the volume the track is in; changed by the locates only + vecgeom::NavigationState mNextState; ///< the state entering the daughter the last ComputeStep hit + vecgeom::NavigationState mStepState; ///< scratch: a trial step + vecgeom::NavigationState mReloScratch; ///< scratch: the state before a crossing + vecgeom::NavigationState mPathScratch; ///< scratch: path comparisons + vecgeom::NavigationState mHistoryState; ///< the state fHistory was built from + vecgeom::NavigationState mEmptyState; ///< permanently empty; its last-exited entry clears others + vecgeom::NavigationState mNormalState; ///< the volume whose surface the last boundary step ended on + bool mHistoryValid = false; + + bool mWouldEnter = false; ///< the last ComputeStep ends by entering a daughter + bool mWouldExit = false; ///< the last ComputeStep ends by leaving the current volume + G4ThreeVector mNextPoint{-1e8, -1e8, -1e8}; + G4ThreeVector mLastLocatedPoint{-1e8, -1e8, -1e8}; + bool mLocatedOnBoundary = false; + G4ThreeVector mSafetyOrig{-1e8, -1e8, -1e8}; ///< the last point a safety was computed for + double mLastSafety = 0.; ///< mm + bool mNormalEnter = false; + bool mNormalValid = false; + G4ThreeVector mNormalPoint{-1e8, -1e8, -1e8}; + + /// Geant4's thresholds for a track that is not moving. + static constexpr int kActionThresholdNoZeroSteps = 10; + static constexpr int kAbandonThresholdNoZeroSteps = 25; + int mNzeroSteps = 0; + + // What the navigator had to work around, reported at the end. + long mZeroStepCount = 0; + long mStuckPushCount = 0; + long mAbandonCount = 0; + long mNegativeSafetyCount = 0; + long mUnmappableHistoryCount = 0; + long mRelocatedResumeCount = 0; + long mNoNormalCount = 0; +}; + +} // namespace o2::simsetup + +#endif diff --git a/Detectors/gconfig/src/VecGeomG4NavigatorBase.cxx b/Detectors/gconfig/src/VecGeomG4NavigatorBase.cxx new file mode 100644 index 0000000000000..7c6ebb22e98f1 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4NavigatorBase.cxx @@ -0,0 +1,128 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "VecGeomG4NavigatorBase.h" + +#include "G4VPhysicalVolume.hh" + +#include + +#include + +namespace +{ +/// Longest run of Geant4 levels one flattened VecGeom placement can stand for. +constexpr std::size_t kMaxChain = 16; +/// Deepest Geant4 touchable the ALICE geometry can produce, with headroom. +constexpr int kMaxDepth = 64; +} // namespace + +namespace o2::simsetup +{ + +G4VPhysicalVolume* VecGeomG4NavigatorBase::historyFromState(vecgeom::NavigationState const& state) +{ + // Collect the Geant4 volumes the state stands for, then keep the levels the history already has + // right and rebuild only from the first difference down: NewLevel composes a transform per level. + G4VPhysicalVolume* want[kMaxDepth]; + int n = 0; + if (!state.IsOutside()) { + const int levels = state.GetCurrentLevel(); + for (int l = 0; l < levels && n < kMaxDepth; ++l) { + auto const* placed = state.At(l); + if (placed == nullptr) { + break; + } + unsigned size = 0; + auto* const* chain = mMap.chain(placed->id(), size); + if (size == 0) { + LOG(fatal) << "VecGeom placement " << placed->GetLabel() << " (id " << placed->id() << ") at level " << l + << " has no Geant4 counterpart"; + } + for (unsigned c = 0; c < size && n < kMaxDepth; ++c) { + want[n++] = chain[c]; + } + } + } + + if (n == 0) { + // Outside the world. A null first entry is how G4NavigationHistory says so. + fHistory.Reset(); + fHistory.SetFirstEntry(nullptr); + return nullptr; + } + + const int depth = static_cast(fHistory.GetDepth()); + int common = 0; + while (common <= depth && common < n && fHistory.GetVolume(common) == want[common]) { + ++common; + } + if (common == 0) { + fHistory.Reset(); + fHistory.SetFirstEntry(want[0]); + common = 1; + } else if (depth >= common) { + fHistory.BackLevel(depth - common + 1); + } + for (int k = common; k < n; ++k) { + // The copy number is what TG4StepManager::CurrentVolID and CurrentVolOffID report, so it has to + // be filled exactly as the TGeo navigator fills it. + fHistory.NewLevel(want[k], kNormal, want[k]->GetCopyNo()); + } + return fHistory.GetTopVolume(); +} + +bool VecGeomG4NavigatorBase::stateFromHistory(vecgeom::NavigationState& state) const +{ + // A run of levels that ends in a Geant4 volume several placements share is settled by comparing + // the run against the placements' recorded chains, which is exact. + state.Clear(); + if (mMap.world() == nullptr || fHistory.GetVolume(0) == nullptr) { + return false; + } + state.Push(mMap.world()); + const std::size_t depth = fHistory.GetDepth(); + std::size_t l = mMap.chainSize(mMap.world()->id()); + while (l <= depth) { + G4VPhysicalVolume* run[kMaxChain]; + vecgeom::VPlacedVolume const* found = nullptr; + std::size_t n = 0; + for (std::size_t e = l; e <= depth && n < kMaxChain; ++e, ++n) { + run[n] = fHistory.GetVolume(e); + auto const* candidate = mMap.toVecGeom(run[n]->GetInstanceID()); + if (candidate == nullptr) { + continue; // an assembly level; the chain reaches further down + } + if (candidate == VecGeomG4Map::ambiguous()) { + for (auto const* c : mMap.candidates(run[n]->GetInstanceID())) { + if (mMap.chainMatches(c, run, n + 1)) { + found = c; + break; + } + } + } else if (mMap.chainMatches(candidate, run, n + 1)) { + found = candidate; + } + if (found != nullptr) { + break; + } + } + if (found == nullptr) { + state.Clear(); + return false; + } + state.Push(found); + l += mMap.chainSize(found->id()); + } + return true; +} + +} // namespace o2::simsetup diff --git a/Detectors/gconfig/src/VecGeomG4NavigatorBase.h b/Detectors/gconfig/src/VecGeomG4NavigatorBase.h new file mode 100644 index 0000000000000..8c276f2f55fd6 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4NavigatorBase.h @@ -0,0 +1,56 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_SIMSETUP_VECGEOMG4NAVIGATORBASE_H_ +#define O2_SIMSETUP_VECGEOMG4NAVIGATORBASE_H_ + +#include "VecGeomG4Map.h" + +#include "G4Navigator.hh" +#include "G4SystemOfUnits.hh" +#include "G4ThreeVector.hh" + +#include +#include + +namespace o2::simsetup +{ + +/// What the two VecGeom navigators share: the correspondence between a VecGeom navigation state and +/// the Geant4 navigation history, which is the touchable the scoring code reads, and the unit +/// conversion. Geant4 works in millimetres, the VecGeom geometry is converted from TGeo in +/// centimetres; every point crossing this boundary is scaled, directions are not. +class VecGeomG4NavigatorBase : public G4Navigator +{ + protected: + using V3 = vecgeom::Vector3D; + static constexpr double kG4ToVG = 1. / CLHEP::cm; + static constexpr double kVGToG4 = CLHEP::cm; + + explicit VecGeomG4NavigatorBase(VecGeomG4Map const& map) : mMap(map) {} + + static V3 toVG(const G4ThreeVector& p) { return {p.x() * kG4ToVG, p.y() * kG4ToVG, p.z() * kG4ToVG}; } + static V3 toDir(const G4ThreeVector& d) { return {d.x(), d.y(), d.z()}; } + + /// Rewrites fHistory to stand for \a state, keeping the levels it already has right. Returns the + /// top volume, or null if the state is outside the world. + G4VPhysicalVolume* historyFromState(vecgeom::NavigationState const& state); + + /// Builds \a state from fHistory, extending past the Geant4 levels assembly flattening dissolved. + /// False if a level matches no VecGeom placement; \a state is then empty. + bool stateFromHistory(vecgeom::NavigationState& state) const; + + VecGeomG4Map const& mMap; +}; + +} // namespace o2::simsetup + +#endif diff --git a/Detectors/gconfig/src/VecGeomNavigation.cxx b/Detectors/gconfig/src/VecGeomNavigation.cxx new file mode 100644 index 0000000000000..83ee6d9ef1f3f --- /dev/null +++ b/Detectors/gconfig/src/VecGeomNavigation.cxx @@ -0,0 +1,137 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "SimSetup/VecGeomNavigation.h" + +#include +#include +#include + +#ifdef O2_WITH_VECGEOM + +// The O2 headers come first on purpose: VecGeom's build interface defines a VECGEOM macro, +// which would otherwise eat the VECGEOM enumerator of MatbudGeomBackend. +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/GeometryManagerParam.h" +#include "SimConfig/G4Params.h" + +#include "VecGeomChecks.h" +#include "VecGeomG4Map.h" +#include "VecGeomG4Navigator.h" + +#include "TG4RootDetectorConstruction.h" +#include "TG4RootNavMgr.h" + +#include "G4EventManager.hh" +#include "G4FieldManager.hh" +#include "G4PropagatorInField.hh" +#include "G4SteppingManager.hh" +#include "G4TrackingManager.hh" +#include "G4TransportationManager.hh" + +#include "TMCManager.h" +#include "TStopwatch.h" + +#endif + +namespace o2::simsetup +{ + +#ifdef O2_WITH_VECGEOM + +bool isVecGeomNavigationAvailable() { return true; } + +void installVecGeomNavigator() +{ + auto const& g4Params = o2::conf::G4Params::Instance(); + + if (TMCManager::Instance() != nullptr) { + LOG(fatal) << "G4.navmode=kVecGeom cannot be used with the multi-engine TMCManager: restoring a " + "geometry state across engines goes through the TGeo navigator"; + } + if (o2::GeometryManagerParam::Instance().useParallelWorld) { + LOG(fatal) << "G4.navmode=kVecGeom cannot be used with GeometryManagerParam.useParallelWorld: VecGeom " + "has no equivalent of the TGeo priority world"; + } + auto* navMgr = TG4RootNavMgr::GetInstance(); + if (navMgr == nullptr || navMgr->GetDetConstruction() == nullptr) { + LOG(fatal) << "G4.navmode=kVecGeom needs the Geant4 geometry built from TGeo by g4root, which is what " + "the geomRoot option provides; no TG4RootNavMgr was found"; + } + auto* detConstruction = navMgr->GetDetConstruction(); + if (!detConstruction->IsConstructed()) { + LOG(fatal) << "The Geant4 geometry has not been built yet; installVecGeomNavigator must be called after " + "the TGeant4 engine has been created"; + } + + if (!g4Params.vecgeomFlattenAssemblies) { + LOG(fatal) << "G4.vecgeomNavigator=kStrict needs G4.vecgeomFlattenAssemblies=true: it enters a daughter by " + "locating inside it, which an assembly cannot answer"; + } + + TStopwatch timer; + timer.Start(); + o2::base::GeometryManager::buildVecGeomGeometry(g4Params.vecgeomFlattenAssemblies); + timer.Stop(); + LOG(info) << "VecGeom geometry built in " << timer.RealTime() << " s"; + + timer.Start(); + // Owned here for the lifetime of the process; the navigator keeps a reference to it. + static VecGeomG4Map map; + map.build(*detConstruction, g4Params.vecgeomFlattenAssemblies); + timer.Stop(); + LOG(info) << "VecGeom to Geant4 map built in " << timer.RealTime() << " s"; + + if (!g4Params.vecgeomCheckVolumes.empty()) { + std::stringstream names(g4Params.vecgeomCheckVolumes); + std::string one; + while (std::getline(names, one, ',')) { + checkVecGeomVolume(one.c_str(), 20, 200); + } + } + if (g4Params.vecgeomCheckRays > 0) { + checkVecGeomRays(static_cast(g4Params.vecgeomCheckRays)); + } + if (g4Params.vecgeomCheckLocation > 0) { + checkVecGeomLocation(static_cast(g4Params.vecgeomCheckLocation)); + } + + auto* navigator = new VecGeomG4Navigator(map, g4Params.vecgeomPushDepth, g4Params.vecgeomZeroSafety); + navigator->SetWorldVolume(detConstruction->GetTopPV()); + + // Same sequence TG4RootNavMgr::SetNavigator uses, run here because by the time the engine + // exists the navigator manager considers itself connected and refuses to swap. + auto* trMgr = G4TransportationManager::GetTransportationManager(); + trMgr->SetNavigatorForTracking(navigator); + auto* fieldMgr = trMgr->GetPropagatorInField()->GetCurrentFieldManager(); + delete trMgr->GetPropagatorInField(); + trMgr->SetPropagatorInField(new G4PropagatorInField(navigator, fieldMgr)); + trMgr->ActivateNavigator(navigator); + if (auto* evtMgr = G4EventManager::GetEventManager()) { + evtMgr->GetTrackingManager()->GetSteppingManager()->SetNavigator(navigator); + } + + LOG(info) << "VecGeom navigator registered with the Geant4 transportation manager"; +} + +#else + +bool isVecGeomNavigationAvailable() { return false; } + +void installVecGeomNavigator() +{ + LOG(fatal) << "G4.navmode=kVecGeom needs O2 built against TGeo2VecGeom and a VecGeom with the BVH navigator of " + "the VNavigator family (BVHNavigatorV), which were not found at configure time"; +} + +#endif + +} // namespace o2::simsetup From 24d46a052018cbd6d7fe9954263a8c83f22ce060 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Tue, 29 Sep 2026 08:25:12 +0200 Subject: [PATCH 3/5] Add a lighter VecGeom navigator that adopts the propagated state This adds G4.vecgeomNavigator=kPropagated, a second VecGeom navigator kept for comparison. - It adopts the state VecGeom propagated during the step instead of relocating at the boundary. - It returns zero safety for the step after a crossing and nudges stuck steps forward. - The default navigator becomes G4.vecgeomNavigator=kRelocating. Co-Authored-By: Claude Opus 5.5 --- Common/SimConfig/include/SimConfig/G4Params.h | 8 + Detectors/gconfig/CMakeLists.txt | 3 +- Detectors/gconfig/src/VecGeomG4Navigator.h | 2 +- .../src/VecGeomG4PropagatingNavigator.cxx | 231 ++++++++++++++++++ .../src/VecGeomG4PropagatingNavigator.h | 73 ++++++ Detectors/gconfig/src/VecGeomNavigation.cxx | 13 +- 6 files changed, 325 insertions(+), 5 deletions(-) create mode 100644 Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx create mode 100644 Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 63ecf63ab06c6..57ae4a28da91c 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -41,6 +41,12 @@ enum class EG4Nav { kVecGeom = 2 /* navigate with VecGeom, on the G4 geometry built from TGeo */ }; +// the Geant4 navigator used with navmode kVecGeom +enum class EVecGeomNav { + kRelocating = 0, /* relocates at the boundary locate, blocking the volume just left (default) */ + kPropagated = 1 /* adopts the state VecGeom propagated during the step; less work per crossing */ +}; + // parameters to influence the G4 engine struct G4Params : public o2::conf::ConfigurableParamHelper { EG4Physics physicsmode = EG4Physics::kFTFP_BERT_EMV_optical; // default physics mode with which to configure G4 @@ -51,6 +57,8 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { EG4Nav navmode = EG4Nav::kTGeo; // geometry navigation mode (default TGeo) // Settings for navmode == kVecGeom; ignored otherwise. + // which of the two VecGeom navigators + EVecGeomNav vecgeomNavigator = EVecGeomNav::kRelocating; double vecgeomPushDepth = 1.e-9; // cm; how far past a face, measured across it, a boundary // point is pushed before it is located bool vecgeomZeroSafety = false; // answer zero to every safety query; conservative, but it diff --git a/Detectors/gconfig/CMakeLists.txt b/Detectors/gconfig/CMakeLists.txt index 276a55af61617..d6fffad875401 100644 --- a/Detectors/gconfig/CMakeLists.txt +++ b/Detectors/gconfig/CMakeLists.txt @@ -31,7 +31,8 @@ set(G4SETUP_SOURCES src/G4Config.cxx src/G4RunConfiguration.cxx src/G4LocalField src/VecGeomNavigation.cxx) if(G4SETUP_WITH_VECGEOM) list(APPEND G4SETUP_SOURCES src/VecGeomG4Map.cxx src/VecGeomChecks.cxx src/VecGeomG4NavigatorBase.cxx - src/VecGeomG4Navigator.cxx) + src/VecGeomG4Navigator.cxx + src/VecGeomG4PropagatingNavigator.cxx) endif() o2_add_library(G4Setup diff --git a/Detectors/gconfig/src/VecGeomG4Navigator.h b/Detectors/gconfig/src/VecGeomG4Navigator.h index 54a3325672ec8..a35bf4e1b0542 100644 --- a/Detectors/gconfig/src/VecGeomG4Navigator.h +++ b/Detectors/gconfig/src/VecGeomG4Navigator.h @@ -19,7 +19,7 @@ namespace o2::simsetup /// A Geant4 tracking navigator that answers every navigation query from VecGeom while keeping the /// Geant4 navigation history, the touchable the scoring code reads, in step with the VecGeom state. -/// It works as G4VecGeomNav's +/// Selected with G4.vecgeomNavigator=kRelocating, the default. It works as G4VecGeomNav's /// TG4VecGeomNavigator does: /// /// - ComputeStep leaves the current volume alone. It records whether the step ends on a boundary diff --git a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx new file mode 100644 index 0000000000000..a33edad3d174f --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx @@ -0,0 +1,231 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "VecGeomG4PropagatingNavigator.h" + +#include "G4TouchableHistory.hh" +#include "G4VPhysicalVolume.hh" + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace +{ +/// Geant4 abandons a track after ~50 zero steps, so a stalled step is nudged forward by this much +/// (in Geant4 units) rather than returned as zero. +const double kNudge = 1.e-3; +} // namespace + +namespace o2::simsetup +{ + +VecGeomG4PropagatingNavigator::VecGeomG4PropagatingNavigator(VecGeomG4Map const& map, bool zeroSafety) + : VecGeomG4NavigatorBase(map), mZeroSafety(zeroSafety) +{ +} + +VecGeomG4PropagatingNavigator::~VecGeomG4PropagatingNavigator() +{ + LOG(info) << "VecGeom navigation: " << mNudgedSteps << " stalled steps nudged forward, " << mGlobalRelocates + << " relocations restarted from the world"; +} + +G4double VecGeomG4PropagatingNavigator::ComputeStep(const G4ThreeVector& globalPoint, const G4ThreeVector& direction, + const G4double proposedStepLength, G4double& newSafety) +{ + newSafety = 0.; + mLastDirection = direction; + mHaveNextState = false; + + auto const* top = mCurState.Top(); + if (top == nullptr) { // the track is outside the world + mWouldEnter = mWouldExit = false; + return kInfinity; + } + auto const* navigator = top->GetLogicalVolume()->GetNavigator(); + + double limit = proposedStepLength * kG4ToVG; + if (!(limit < vecgeom::kInfLength)) { + limit = vecgeom::kInfLength; + } + + // VecGeom's own combined entry point, rather than a step followed by a relocation of our own: it + // is the one that knows how to descend through an assembly, whose placed volume can never be the + // answer because it has no DistanceToOut. + double vgSafety = 0.; + const double vgStep = navigator->ComputeStepAndSafetyAndPropagatedState( + toVG(globalPoint), toDir(direction), limit, mCurState, mNextState, !mOnBoundary && !mZeroSafety, vgSafety); + mOnBoundary = false; + newSafety = (vgSafety > 0. && !mZeroSafety) ? vgSafety * kVGToG4 : 0.; + + G4double step = vgStep * kVGToG4; + if (mNextState.IsOnBoundary()) { + // Entering a daughter deepens the state, possibly by more than one level when an assembly + // stands in between; leaving the current volume does not. + mWouldEnter = mNextState.GetCurrentLevel() > mCurState.GetCurrentLevel(); + mWouldExit = !mWouldEnter; + mNextPoint = globalPoint + step * direction; + mHaveNextState = true; + } else { + mWouldEnter = mWouldExit = false; + step = kInfinity; + } + + if (vgStep < 0.) { + // A negative distance means the state and the point disagree. Nudge forward and relocate from + // the world on the next call rather than propagating the inconsistency. + mForceReInit = true; + mHaveNextState = false; + ++mNudgedSteps; + mZeroSteps = 0; + return kNudge; + } + if (step < 1.e-10) { + if (++mZeroSteps > 4) { + mForceReInit = true; + mHaveNextState = false; + ++mNudgedSteps; + return kNudge; + } + } else { + mForceReInit = false; + mZeroSteps = 0; + } + return step; +} + +G4VPhysicalVolume* VecGeomG4PropagatingNavigator::ResetHierarchyAndLocate(const G4ThreeVector&, const G4ThreeVector&, + const G4TouchableHistory& history) +{ + // Geant4 hands back a touchable it saved earlier, e.g. when resuming a track whose secondaries + // were followed first. The VecGeom state is rebuilt from it. + fEnteredDaughter = false; + fExitedMother = false; + mWouldEnter = false; + mWouldExit = false; + mOnBoundary = false; + mHaveNextState = false; + fHistory = *history.GetHistory(); + if (!stateFromHistory(mCurState) && fHistory.GetVolume(0) != nullptr) { + LOG(fatal) << "Geant4 handed back a touchable that matches no VecGeom path"; + } + mPrevState = mCurState; + return fHistory.GetTopVolume(); +} + +G4VPhysicalVolume* VecGeomG4PropagatingNavigator::LocateGlobalPointAndSetup(const G4ThreeVector& point, + const G4ThreeVector*, + const G4bool relativeSearch, const G4bool) +{ + bool onBoundary = fWasLimitedByGeometry; + if (mHaveNextState && point.diff2(mNextPoint) < 1.e-16) { + onBoundary = true; + } + + mPrevState = mCurState; + + if (!mForceReInit && relativeSearch && onBoundary && mHaveNextState) { + // The state on the far side of the boundary was already worked out, and relocated, by the step + // that found it. Adopting it is cheaper than locating again. + mCurState = mNextState; + } else if (mForceReInit || !relativeSearch || onBoundary) { + mCurState.Clear(); + vecgeom::GlobalLocator::LocateGlobalPoint(vecgeom::GeoManager::Instance().GetWorld(), toVG(point), mCurState, true); + mForceReInit = false; + ++mGlobalRelocates; + } + // Otherwise the point only moved inside the volume the state already names. + + auto* target = historyFromState(mCurState); + mCrossed = onBoundary; + if (onBoundary) { + fExitedMother = mWouldExit; + fEnteredDaughter = mWouldEnter; + mOnBoundary = true; + } + mHaveNextState = false; + return target; +} + +void VecGeomG4PropagatingNavigator::LocateGlobalPointWithinVolume(const G4ThreeVector&) +{ + // The track moved inside the volume it is already in, so only the boundary flags change. + mWouldEnter = false; + mWouldExit = false; + mOnBoundary = false; + mCrossed = false; + mHaveNextState = false; + fEnteredDaughter = false; + fExitedMother = false; +} + +G4double VecGeomG4PropagatingNavigator::ComputeSafety(const G4ThreeVector& globalPoint, const G4double, const G4bool) +{ + if (mZeroSafety || mOnBoundary || mCrossed || fEnteredDaughter || fExitedMother || mWouldEnter || mWouldExit) { + return 0.; + } + auto const* top = mCurState.Top(); + if (top == nullptr) { + return 0.; + } + const double safety = top->GetLogicalVolume()->GetSafetyEstimator()->ComputeSafety(toVG(globalPoint), mCurState); + return (safety > 0.) ? safety * kVGToG4 : 0.; +} + +G4ThreeVector VecGeomG4PropagatingNavigator::GetGlobalExitNormal(const G4ThreeVector& point, G4bool* valid) +{ + // The surface just crossed belongs to the volume that was left when the step exited a mother, and + // to the volume that was entered when it entered a daughter. + auto const& state = mWouldExit ? mPrevState : mCurState; + auto const* volume = state.Top(); + if (volume == nullptr) { + *valid = false; + return G4ThreeVector(0., 0., 1.); + } + + vecgeom::Transformation3D m; + state.TopMatrix(m); + V3 normal; + volume->Normal(m.Transform(toVG(point)), normal); + V3 global = m.InverseTransformDirection(normal); + + // Oriented along the direction of motion, as TGeo's FindNormalFast does. + const V3 dir = toDir(mLastDirection); + if (global.Dot(dir) < 0.) { + global = -global; + } + // VecGeom's Normal() also answers whether the point was on the surface; Geant4 hands back points a + // few nanometres off the face, so only a degenerate vector is refused. + const double mag2 = global.Mag2(); + *valid = std::isfinite(mag2) && mag2 > 0.25; + return G4ThreeVector(global[0], global[1], global[2]); +} + +G4ThreeVector VecGeomG4PropagatingNavigator::GetLocalExitNormal(G4bool* valid) +{ + // By convention the local normal is expressed in the frame of the final volume. + const G4ThreeVector global = GetGlobalExitNormal(mNextPoint, valid); + vecgeom::Transformation3D m; + mCurState.TopMatrix(m); + const V3 local = m.TransformDirection(toDir(global)); + return G4ThreeVector(local[0], local[1], local[2]); +} + +} // namespace o2::simsetup diff --git a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h new file mode 100644 index 0000000000000..108730f8af552 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h @@ -0,0 +1,73 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_SIMSETUP_VECGEOMG4PROPAGATINGNAVIGATOR_H_ +#define O2_SIMSETUP_VECGEOMG4PROPAGATINGNAVIGATOR_H_ + +#include "VecGeomG4NavigatorBase.h" + +namespace o2::simsetup +{ + +/// A lighter VecGeom navigator for Geant4, selected with G4.vecgeomNavigator=kPropagated. ComputeStep +/// uses VecGeom's combined step-and-relocate, and a boundary locate adopts the state that step +/// propagated instead of locating again. Safety is zero for the whole step after a crossing, a stuck +/// step is nudged forward, and the exit normal follows the direction of motion. It does less work per +/// step and per crossing than VecGeomG4Navigator. +class VecGeomG4PropagatingNavigator : public VecGeomG4NavigatorBase +{ + public: + VecGeomG4PropagatingNavigator(VecGeomG4Map const& map, bool zeroSafety); + ~VecGeomG4PropagatingNavigator() override; + + G4double ComputeStep(const G4ThreeVector& globalPoint, const G4ThreeVector& direction, + const G4double proposedStepLength, G4double& newSafety) override; + + G4VPhysicalVolume* ResetHierarchyAndLocate(const G4ThreeVector& point, const G4ThreeVector& direction, + const G4TouchableHistory& history) override; + + G4VPhysicalVolume* LocateGlobalPointAndSetup(const G4ThreeVector& point, const G4ThreeVector* direction = nullptr, + const G4bool relativeSearch = true, + const G4bool ignoreDirection = true) override; + + void LocateGlobalPointWithinVolume(const G4ThreeVector& position) override; + + G4double ComputeSafety(const G4ThreeVector& globalPoint, const G4double proposedMaxLength = DBL_MAX, + const G4bool keepState = true) override; + + G4ThreeVector GetLocalExitNormal(G4bool* valid) override; + G4ThreeVector GetGlobalExitNormal(const G4ThreeVector& point, G4bool* valid) override; + + private: + vecgeom::NavigationState mCurState; ///< where the track is now + vecgeom::NavigationState mNextState; ///< where the last computed step would put it + vecgeom::NavigationState mPrevState; ///< where it was before the last boundary crossing + + bool mZeroSafety = false; + bool mHaveNextState = false; ///< mNextState holds the result of a geometry-limited step + + G4ThreeVector mNextPoint{-1e8, -1e8, -1e8}; ///< where the last computed step ends + G4ThreeVector mLastDirection{0, 0, 1}; ///< direction of the last computed step + + bool mWouldEnter = false; ///< the last step ends by entering a daughter + bool mWouldExit = false; ///< the last step ends by leaving the current volume + bool mOnBoundary = false; ///< the current point sits on a boundary + bool mForceReInit = false; ///< next locate must start from the world, the state is suspect + bool mCrossed = false; ///< the last locate acted on a boundary crossing + + int mZeroSteps = 0; + long mNudgedSteps = 0; + long mGlobalRelocates = 0; +}; + +} // namespace o2::simsetup + +#endif diff --git a/Detectors/gconfig/src/VecGeomNavigation.cxx b/Detectors/gconfig/src/VecGeomNavigation.cxx index 83ee6d9ef1f3f..68dbc96f883c3 100644 --- a/Detectors/gconfig/src/VecGeomNavigation.cxx +++ b/Detectors/gconfig/src/VecGeomNavigation.cxx @@ -26,6 +26,7 @@ #include "VecGeomChecks.h" #include "VecGeomG4Map.h" #include "VecGeomG4Navigator.h" +#include "VecGeomG4PropagatingNavigator.h" #include "TG4RootDetectorConstruction.h" #include "TG4RootNavMgr.h" @@ -72,8 +73,8 @@ void installVecGeomNavigator() "the TGeant4 engine has been created"; } - if (!g4Params.vecgeomFlattenAssemblies) { - LOG(fatal) << "G4.vecgeomNavigator=kStrict needs G4.vecgeomFlattenAssemblies=true: it enters a daughter by " + if (!g4Params.vecgeomFlattenAssemblies && g4Params.vecgeomNavigator == o2::conf::EVecGeomNav::kRelocating) { + LOG(fatal) << "G4.vecgeomNavigator=kRelocating needs G4.vecgeomFlattenAssemblies=true: it enters a daughter by " "locating inside it, which an assembly cannot answer"; } @@ -104,7 +105,13 @@ void installVecGeomNavigator() checkVecGeomLocation(static_cast(g4Params.vecgeomCheckLocation)); } - auto* navigator = new VecGeomG4Navigator(map, g4Params.vecgeomPushDepth, g4Params.vecgeomZeroSafety); + G4Navigator* navigator = nullptr; + if (g4Params.vecgeomNavigator == o2::conf::EVecGeomNav::kPropagated) { + LOG(info) << "VecGeom navigation: propagating navigator"; + navigator = new VecGeomG4PropagatingNavigator(map, g4Params.vecgeomZeroSafety); + } else { + navigator = new VecGeomG4Navigator(map, g4Params.vecgeomPushDepth, g4Params.vecgeomZeroSafety); + } navigator->SetWorldVolume(detConstruction->GetTopPV()); // Same sequence TG4RootNavMgr::SetNavigator uses, run here because by the time the engine From 853ab5d95f5b7fa0f15f3fdb4a757f3df01abb20 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Wed, 30 Sep 2026 17:42:01 +0200 Subject: [PATCH 4/5] 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 --- Detectors/gconfig/src/VecGeomG4Navigator.cxx | 29 ++++++++++++++---- Detectors/gconfig/src/VecGeomG4Navigator.h | 30 +++++++++++++++---- .../gconfig/src/VecGeomG4NavigatorBase.h | 27 +++++++++++++++++ .../src/VecGeomG4PropagatingNavigator.cxx | 14 +++++++++ .../src/VecGeomG4PropagatingNavigator.h | 3 ++ 5 files changed, 91 insertions(+), 12 deletions(-) diff --git a/Detectors/gconfig/src/VecGeomG4Navigator.cxx b/Detectors/gconfig/src/VecGeomG4Navigator.cxx index 66418efe7de64..8644f9a29e625 100644 --- a/Detectors/gconfig/src/VecGeomG4Navigator.cxx +++ b/Detectors/gconfig/src/VecGeomG4Navigator.cxx @@ -80,7 +80,7 @@ void VecGeomG4Navigator::locateFromWorld(const V3& point) /// surely as one at normal incidence, and it is kept small: TOF has layers 2.4e-8 cm apart, and a /// fixed push along the direction steps over them. It never goes below a thousand times the rounding /// of the largest coordinate, which matters far from the origin. -double VecGeomG4Navigator::boundaryPush(const V3& point, const V3& dir) const +double VecGeomG4Navigator::boundaryPush(const V3& point, const V3& dir) { constexpr double kMaxPush = 1.e-4; // cm along the direction 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 // The face just crossed: the entered daughter's, or the current volume's own. vecgeom::NavigationState const& st = mWouldEnter ? mNextState : mCurState; double cosn = 1.; + mPushNormalValid = false; if (!st.IsOutside() && st.Top() != nullptr) { vecgeom::Transformation3D m; st.TopMatrix(m); @@ -98,6 +99,10 @@ double VecGeomG4Navigator::boundaryPush(const V3& point, const V3& dir) const if (c > 0. && n.Mag2() > 0.5) { cosn = c; } + if (n.Mag2() > 0.5) { + mPushNormal = m.InverseTransformDirection(n); + mPushNormalValid = true; + } } return std::max(rounding, std::min(mPushDepth / cosn, kMaxPush)); } @@ -123,10 +128,20 @@ G4double VecGeomG4Navigator::ComputeStep(const G4ThreeVector& globalPoint, const newSafety = mLastSafety; } - // The step is computed on a copy. The state, and the volume blocked by the last exit, stay as - // the locate left them for every call until the next locate, which is what the field propagator - // relies on when it calls this from trial points along the curve. + // The step is computed on a copy. The state stays as the locate left it for every call until the + // next locate, which is what the field propagator relies on when it calls this from trial points + // along the curve. The volume the last crossing left is blocked in the first call only, and only + // while the direction points away from it; a track turning back into it must see its boundary. mStepState = mCurState; + if (mExitBlockPending) { + const bool block = samePoint(globalPoint, mLastLocatedPoint) && + (mExitNormalFromPush ? directionLeaves(mPushNormal, toDir(direction)) + : directionLeaves(mExitedState, toVG(globalPoint), toDir(direction))); + if (!block) { + mStepState.SetLastExited(mEmptyState.GetLastExitedState()); + } + clearLastExited(); + } const double limit = std::min(proposedStepLength * kG4ToVG, static_cast(vecgeom::kInfLength)); double safety = 0.; 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 } left = mCurState.Top(); mCurState.SetLastExited(); + setExited(mCurState, false); if (mCurState.GetCurrentLevel() <= 1) { mCurState.Clear(); // nothing to travel in even in the world: the track left it return; @@ -307,8 +323,8 @@ G4VPhysicalVolume* VecGeomG4Navigator::LocateGlobalPointAndSetup(const G4ThreeVe mLocatedOnBoundary = true; } else if (crossing) { // Out of the current volume: up until the point is contained, then down, never back into the - // volume just left; that volume stays blocked at zero distance for the next step, as - // G4Navigator's fBlockedPhysicalVolume. + // volume just left; that volume is blocked at zero distance in the next ComputeStep while the + // direction points away from it, as G4Navigator's fBlockedPhysicalVolume. mReloScratch = mCurState; if (mCurState.GetCurrentLevel() <= 1) { mCurState.Clear(); // left the world @@ -318,6 +334,7 @@ G4VPhysicalVolume* VecGeomG4Navigator::LocateGlobalPointAndSetup(const G4ThreeVe vecgeom::GlobalLocator::RelocatePointFromPathForceDifferent(m.Transform(q), mCurState); mReloScratch.SetLastExited(); mCurState.SetLastExited(mReloScratch.GetLastExitedState()); + setExited(mReloScratch, true); } mLocatedOnBoundary = true; } else { diff --git a/Detectors/gconfig/src/VecGeomG4Navigator.h b/Detectors/gconfig/src/VecGeomG4Navigator.h index a35bf4e1b0542..12e281c040016 100644 --- a/Detectors/gconfig/src/VecGeomG4Navigator.h +++ b/Detectors/gconfig/src/VecGeomG4Navigator.h @@ -24,9 +24,10 @@ namespace o2::simsetup /// /// - ComputeStep leaves the current volume alone. It records whether the step ends on a boundary /// and, for a daughter hit, the state that enters it. -/// - The locate on a boundary relocates: into the recorded daughter, or out of the current volume, -/// which then stays blocked for the next step. The point is first pushed across the face by a -/// small depth, and afterwards leaves every volume it is flush with and heading out of. +/// - The locate on a boundary relocates: into the recorded daughter, or out of the current volume. +/// The point is first pushed across the face by a small depth, and afterwards leaves every volume +/// it is flush with and heading out of. As in Geant4, the volume left is blocked only in the first +/// ComputeStep after the exit, and only while the direction points away from it. /// - Safety is zero only at the boundary point itself. /// /// Unlike TG4VecGeomNavigator, the VecGeom geometry is converted from TGeo, not from Geant4, so one @@ -66,9 +67,21 @@ class VecGeomG4Navigator : public VecGeomG4NavigatorBase void leaveFlushVolumes(const V3& point, const V3& dir, int minLevel, vecgeom::VPlacedVolume const* avoid); /// Sets fEnteredDaughter and fExitedMother from the paths before (mReloScratch) and after a crossing. void updateCrossingFlags(bool entering); - void clearLastExited() { mCurState.SetLastExited(mEmptyState.GetLastExitedState()); } - /// How far (cm) a boundary point is pushed along the direction before it is located. - double boundaryPush(const V3& point, const V3& dir) const; + void clearLastExited() + { + mCurState.SetLastExited(mEmptyState.GetLastExitedState()); + mExitBlockPending = false; + } + /// Records the top of \a st as the volume just left, blocked in the next ComputeStep by Geant4's rule. + void setExited(vecgeom::NavigationState const& st, bool normalFromPush) + { + mExitedState = st; + mExitBlockPending = true; + mExitNormalFromPush = normalFromPush && mPushNormalValid; + } + /// How far (cm) a boundary point is pushed along the direction before it is located. Keeps the + /// global normal of the face in mPushNormal, for the exit block. + double boundaryPush(const V3& point, const V3& dir); /// The normal of the surface the last geometry-limited ComputeStep ended on, global, unit length. bool computeExitNormal(const G4ThreeVector& point, V3& globalNormal) const; @@ -83,6 +96,11 @@ class VecGeomG4Navigator : public VecGeomG4NavigatorBase vecgeom::NavigationState mHistoryState; ///< the state fHistory was built from vecgeom::NavigationState mEmptyState; ///< permanently empty; its last-exited entry clears others vecgeom::NavigationState mNormalState; ///< the volume whose surface the last boundary step ended on + vecgeom::NavigationState mExitedState; ///< the volume the last crossing left + bool mExitBlockPending = false; ///< the next ComputeStep is the first after that exit + bool mExitNormalFromPush = false; ///< mPushNormal is the outward normal of the volume left + V3 mPushNormal; ///< global normal of the face the last push crossed + bool mPushNormalValid = false; bool mHistoryValid = false; bool mWouldEnter = false; ///< the last ComputeStep ends by entering a daughter diff --git a/Detectors/gconfig/src/VecGeomG4NavigatorBase.h b/Detectors/gconfig/src/VecGeomG4NavigatorBase.h index 8c276f2f55fd6..6f86368a2b221 100644 --- a/Detectors/gconfig/src/VecGeomG4NavigatorBase.h +++ b/Detectors/gconfig/src/VecGeomG4NavigatorBase.h @@ -18,8 +18,10 @@ #include "G4SystemOfUnits.hh" #include "G4ThreeVector.hh" +#include #include #include +#include namespace o2::simsetup { @@ -48,6 +50,31 @@ class VecGeomG4NavigatorBase : public G4Navigator /// False if a level matches no VecGeom placement; \a state is then empty. bool stateFromHistory(vecgeom::NavigationState& state) const; + /// Geant4's rule for the volume a track just left (G4NormalNavigation, G4VoxelNavigation): it is + /// blocked only while the direction points away from it, along its outward normal \a n (global). + [[gnu::always_inline]] static bool directionLeaves(const V3& n, const V3& dir) + { + constexpr double kMinExitingNormalCosine = 1e-3; // as in G4NormalNavigation + return n.Dot(dir) >= kMinExitingNormalCosine; + } + + /// The same, with the normal taken from \a exited, which has the volume left on top. + [[gnu::always_inline]] static bool directionLeaves(vecgeom::NavigationState const& exited, const V3& point, + const V3& dir) + { + if (exited.IsOutside() || exited.Top() == nullptr) { + return false; + } + vecgeom::Transformation3D m; + exited.TopMatrix(m); + V3 n; + exited.Top()->GetUnplacedVolume()->Normal(m.Transform(point), n); + if (!(n.Mag2() > 0.5)) { + return false; // no valid exit normal: Geant4 does not block either + } + return directionLeaves(m.InverseTransformDirection(n), dir); + } + VecGeomG4Map const& mMap; }; diff --git a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx index a33edad3d174f..4057c9c65dd3e 100644 --- a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx +++ b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx @@ -39,6 +39,7 @@ namespace o2::simsetup VecGeomG4PropagatingNavigator::VecGeomG4PropagatingNavigator(VecGeomG4Map const& map, bool zeroSafety) : VecGeomG4NavigatorBase(map), mZeroSafety(zeroSafety) { + mEmptyState.Clear(); } VecGeomG4PropagatingNavigator::~VecGeomG4PropagatingNavigator() @@ -69,6 +70,16 @@ G4double VecGeomG4PropagatingNavigator::ComputeStep(const G4ThreeVector& globalP // VecGeom's own combined entry point, rather than a step followed by a relocation of our own: it // is the one that knows how to descend through an assembly, whose placed volume can never be the // answer because it has no DistanceToOut. + // The adopted state marks the volume the crossing left. As in Geant4, it is blocked in the first + // step after the exit only, and only while the direction points away from it; a track turning + // back into it must see its boundary. + const bool block = mExitBlockPending && globalPoint.diff2(mLocatedPoint) < 1.e-20 && + directionLeaves(mPrevState, toVG(globalPoint), toDir(direction)); + mExitBlockPending = false; + if (!block) { + mCurState.SetLastExited(mEmptyState.GetLastExitedState()); + } + double vgSafety = 0.; const double vgStep = navigator->ComputeStepAndSafetyAndPropagatedState( toVG(globalPoint), toDir(direction), limit, mCurState, mNextState, !mOnBoundary && !mZeroSafety, vgSafety); @@ -140,11 +151,14 @@ G4VPhysicalVolume* VecGeomG4PropagatingNavigator::LocateGlobalPointAndSetup(cons } mPrevState = mCurState; + mLocatedPoint = point; + mExitBlockPending = false; if (!mForceReInit && relativeSearch && onBoundary && mHaveNextState) { // The state on the far side of the boundary was already worked out, and relocated, by the step // that found it. Adopting it is cheaper than locating again. mCurState = mNextState; + mExitBlockPending = mWouldExit; } else if (mForceReInit || !relativeSearch || onBoundary) { mCurState.Clear(); vecgeom::GlobalLocator::LocateGlobalPoint(vecgeom::GeoManager::Instance().GetWorld(), toVG(point), mCurState, true); diff --git a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h index 108730f8af552..4ea5407d348ad 100644 --- a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h +++ b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h @@ -50,6 +50,7 @@ class VecGeomG4PropagatingNavigator : public VecGeomG4NavigatorBase vecgeom::NavigationState mCurState; ///< where the track is now vecgeom::NavigationState mNextState; ///< where the last computed step would put it vecgeom::NavigationState mPrevState; ///< where it was before the last boundary crossing + vecgeom::NavigationState mEmptyState; ///< permanently empty; its last-exited entry clears others bool mZeroSafety = false; bool mHaveNextState = false; ///< mNextState holds the result of a geometry-limited step @@ -62,6 +63,8 @@ class VecGeomG4PropagatingNavigator : public VecGeomG4NavigatorBase bool mOnBoundary = false; ///< the current point sits on a boundary bool mForceReInit = false; ///< next locate must start from the world, the state is suspect bool mCrossed = false; ///< the last locate acted on a boundary crossing + bool mExitBlockPending = false; ///< the next step is the first after leaving mPrevState's volume + G4ThreeVector mLocatedPoint{-1e8, -1e8, -1e8}; ///< where the last locate put the track int mZeroSteps = 0; long mNudgedSteps = 0; From 132c49d43192655e41b0aaee292c3ef4c7f3f83d Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 30 Sep 2026 17:04:34 +0000 Subject: [PATCH 5/5] Please consider the following formatting changes --- .../gconfig/src/VecGeomG4PropagatingNavigator.h | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h index 4ea5407d348ad..4ccf469966304 100644 --- a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h +++ b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h @@ -47,9 +47,9 @@ class VecGeomG4PropagatingNavigator : public VecGeomG4NavigatorBase G4ThreeVector GetGlobalExitNormal(const G4ThreeVector& point, G4bool* valid) override; private: - vecgeom::NavigationState mCurState; ///< where the track is now - vecgeom::NavigationState mNextState; ///< where the last computed step would put it - vecgeom::NavigationState mPrevState; ///< where it was before the last boundary crossing + vecgeom::NavigationState mCurState; ///< where the track is now + vecgeom::NavigationState mNextState; ///< where the last computed step would put it + vecgeom::NavigationState mPrevState; ///< where it was before the last boundary crossing vecgeom::NavigationState mEmptyState; ///< permanently empty; its last-exited entry clears others bool mZeroSafety = false; @@ -58,11 +58,11 @@ class VecGeomG4PropagatingNavigator : public VecGeomG4NavigatorBase G4ThreeVector mNextPoint{-1e8, -1e8, -1e8}; ///< where the last computed step ends G4ThreeVector mLastDirection{0, 0, 1}; ///< direction of the last computed step - bool mWouldEnter = false; ///< the last step ends by entering a daughter - bool mWouldExit = false; ///< the last step ends by leaving the current volume - bool mOnBoundary = false; ///< the current point sits on a boundary - bool mForceReInit = false; ///< next locate must start from the world, the state is suspect - bool mCrossed = false; ///< the last locate acted on a boundary crossing + bool mWouldEnter = false; ///< the last step ends by entering a daughter + bool mWouldExit = false; ///< the last step ends by leaving the current volume + bool mOnBoundary = false; ///< the current point sits on a boundary + bool mForceReInit = false; ///< next locate must start from the world, the state is suspect + bool mCrossed = false; ///< the last locate acted on a boundary crossing bool mExitBlockPending = false; ///< the next step is the first after leaving mPrevState's volume G4ThreeVector mLocatedPoint{-1e8, -1e8, -1e8}; ///< where the last locate put the track