diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 97c3cc4ddb412..57ae4a28da91c 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -36,8 +36,15 @@ 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 */ +}; + +// 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 @@ -49,6 +56,22 @@ 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 + // 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/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) diff --git a/Detectors/gconfig/CMakeLists.txt b/Detectors/gconfig/CMakeLists.txt index a1a2b426f2bb4..d6fffad875401 100644 --- a/Detectors/gconfig/CMakeLists.txt +++ b/Detectors/gconfig/CMakeLists.txt @@ -14,11 +14,40 @@ 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 + src/VecGeomG4PropagatingNavigator.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..8644f9a29e625 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4Navigator.cxx @@ -0,0 +1,513 @@ +// 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) +{ + 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.; + mPushNormalValid = false; + 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; + } + if (n.Mag2() > 0.5) { + mPushNormal = m.InverseTransformDirection(n); + mPushNormalValid = true; + } + } + 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 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, + 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(); + setExited(mCurState, false); + 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 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 + } else { + vecgeom::Transformation3D m; + mCurState.TopMatrix(m); + vecgeom::GlobalLocator::RelocatePointFromPathForceDifferent(m.Transform(q), mCurState); + mReloScratch.SetLastExited(); + mCurState.SetLastExited(mReloScratch.GetLastExitedState()); + setExited(mReloScratch, true); + } + 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..12e281c040016 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4Navigator.h @@ -0,0 +1,134 @@ +// 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. +/// 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 +/// 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. +/// 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 +/// 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()); + 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; + + 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 + 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 + 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..6f86368a2b221 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4NavigatorBase.h @@ -0,0 +1,83 @@ +// 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 +#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; + + /// 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; +}; + +} // namespace o2::simsetup + +#endif diff --git a/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx new file mode 100644 index 0000000000000..4057c9c65dd3e --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.cxx @@ -0,0 +1,245 @@ +// 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) +{ + mEmptyState.Clear(); +} + +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. + // 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); + 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; + 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); + 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..4ccf469966304 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomG4PropagatingNavigator.h @@ -0,0 +1,76 @@ +// 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 + 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 + + 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 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; + long mGlobalRelocates = 0; +}; + +} // namespace o2::simsetup + +#endif diff --git a/Detectors/gconfig/src/VecGeomNavigation.cxx b/Detectors/gconfig/src/VecGeomNavigation.cxx new file mode 100644 index 0000000000000..68dbc96f883c3 --- /dev/null +++ b/Detectors/gconfig/src/VecGeomNavigation.cxx @@ -0,0 +1,144 @@ +// 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 "VecGeomG4PropagatingNavigator.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 && 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"; + } + + 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)); + } + + 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 + // 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