From 228199a142cd5f23e09cf9ef0d87ce10c4ef9bd5 Mon Sep 17 00:00:00 2001 From: SeverinDiederichs Date: Wed, 18 Mar 2026 16:35:12 +0100 Subject: [PATCH 1/6] Move HepEm host-data preparation out of transport --- CMakeLists.txt | 1 + include/AdePT/core/AdePTHepEmHostData.hh | 64 +++++++++++++++++++ include/AdePT/core/AsyncAdePTTransport.cuh | 64 +++++++------------ include/AdePT/core/AsyncAdePTTransport.hh | 18 +++--- include/AdePT/core/AsyncAdePTTransport.icc | 53 +++++++-------- .../integration/AdePTGeant4Integration.hh | 6 +- src/AdePTGeant4Integration.cpp | 15 ++--- src/AdePTHepEmHostData.cpp | 54 ++++++++++++++++ src/AdePTTrackingManager.cc | 23 ++++--- 9 files changed, 196 insertions(+), 102 deletions(-) create mode 100644 include/AdePT/core/AdePTHepEmHostData.hh create mode 100644 src/AdePTHepEmHostData.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 00626f47..51297c07 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -444,6 +444,7 @@ endif() # Build Targets #----------------------------------------------------------------------------# set(ADEPT_G4_INTEGRATION_SRCS + src/AdePTHepEmHostData.cpp src/G4HepEmTrackingManagerSpecialized.cc src/AdePTTrackingManager.cc src/G4EmStandardPhysics_AdePT.cc diff --git a/include/AdePT/core/AdePTHepEmHostData.hh b/include/AdePT/core/AdePTHepEmHostData.hh new file mode 100644 index 00000000..1d164763 --- /dev/null +++ b/include/AdePT/core/AdePTHepEmHostData.hh @@ -0,0 +1,64 @@ +// SPDX-FileCopyrightText: 2026 CERN +// SPDX-License-Identifier: Apache-2.0 + +#ifndef ADEPT_HEPEM_HOST_DATA_HH +#define ADEPT_HEPEM_HOST_DATA_HH + +#include + +struct G4HepEmConfig; +struct G4HepEmData; +struct G4HepEmParameters; + +namespace AsyncAdePT { + +/// @brief Owns the host-side HepEm tables rebuilt from a `G4HepEmConfig`. +/// @details +/// The Geant4 integration side prepares one of these objects before the shared +/// transport is created. The helper owns the rebuilt `G4HepEmData`, while the +/// `G4HepEmParameters` pointer is only borrowed from the `G4HepEmConfig`. +/// +/// AdePT also creates the GPU mirror for the borrowed parameters. That GPU-side +/// allocation is released here in the destructor, because this helper owns the +/// corresponding upload lifecycle even though the host parameters themselves are +/// still owned by G4HepEm. +class HepEmHostData { +public: + /// @brief Rebuild all host-side HepEm tables needed by AdePT from the given config. + explicit HepEmHostData(G4HepEmConfig *hepEmConfig); + + /// @brief Release the GPU mirror of the borrowed HepEm parameters. + /// @details This performs a CUDA-side free through `FreeG4HepEmParametersOnGPU`. + ~HepEmHostData(); + + HepEmHostData(const HepEmHostData &) = delete; + HepEmHostData &operator=(const HepEmHostData &) = delete; + + HepEmHostData(HepEmHostData &&) noexcept = default; + HepEmHostData &operator=(HepEmHostData &&) noexcept = default; + + /// @brief Access the owned host-side HepEm data tables. + G4HepEmData *GetData() const { return fData.get(); } + + /// @brief Access the borrowed HepEm parameter block from the Geant4 config. + G4HepEmParameters *GetParameters() const { return fParameters; } + +private: + /// @brief Deletes the outer `G4HepEmData` object after first freeing all tables it owns. + /// @details `FreeG4HepEmData` releases both the host-side tables and any device-side + /// mirrors embedded in the `G4HepEmData` object, but it does not delete the outer + /// `G4HepEmData` allocation itself. This deleter performs both steps. + struct DataDeleter { + void operator()(G4HepEmData *data) const; + }; + + /// Owned host-side HepEm tables rebuilt for AdePT. + std::unique_ptr fData; + + /// Non-owning pointer to the Geant4-owned HepEm parameter block. + G4HepEmParameters *fParameters{nullptr}; +}; + +} // namespace AsyncAdePT + +#endif diff --git a/include/AdePT/core/AsyncAdePTTransport.cuh b/include/AdePT/core/AsyncAdePTTransport.cuh index f832c57c..5c749b23 100644 --- a/include/AdePT/core/AsyncAdePTTransport.cuh +++ b/include/AdePT/core/AsyncAdePTTransport.cuh @@ -45,9 +45,6 @@ using SteppingAction = adept::SteppingAction::Action; #endif #include -#include -#include -#include #include #include #include @@ -536,48 +533,32 @@ void CopySurfaceModelToGPU() #endif } -G4HepEmState *InitG4HepEm(G4HepEmConfig *hepEmConfig) +void UploadG4HepEmToGPU(G4HepEmData *hepEmData, G4HepEmParameters *hepEmParameters) { - // here we call everything from InitG4HepEmState, as we need to provide the parameters from the G4HepEmConfig and do - // not want to initialize to the default values - auto state = new G4HepEmState; - - // Use the config-provided parameters - state->fParameters = hepEmConfig->GetG4HepEmParameters(); - - // Initialize data and fill each subtable using its initialize function - state->fData = new G4HepEmData; - InitG4HepEmData(state->fData); - InitMaterialAndCoupleData(state->fData, state->fParameters); - - // electrons, positrons, gamma - InitElectronData(state->fData, state->fParameters, true); - InitElectronData(state->fData, state->fParameters, false); - InitGammaData(state->fData, state->fParameters); - - G4HepEmMatCutData *cutData = state->fData->fTheMatCutData; - G4cout << "fNumG4MatCuts = " << cutData->fNumG4MatCuts << ", fNumMatCutData = " << cutData->fNumMatCutData << G4endl; + if (hepEmData == nullptr || hepEmParameters == nullptr) { + throw std::runtime_error("UploadG4HepEmToGPU requires non-null HepEm data and parameters."); + } - // Copy to GPU. - CopyG4HepEmDataToGPU(state->fData); - CopyG4HepEmParametersToGPU(state->fParameters); + // Copy the prepared host-side HepEm data to the GPU. + CopyG4HepEmDataToGPU(hepEmData); + CopyG4HepEmParametersToGPU(hepEmParameters); // Create G4HepEmParameters with the device pointer - G4HepEmParameters parametersOnDevice = *state->fParameters; - parametersOnDevice.fParametersPerRegion = state->fParameters->fParametersPerRegion_gpu; + G4HepEmParameters parametersOnDevice = *hepEmParameters; + parametersOnDevice.fParametersPerRegion = hepEmParameters->fParametersPerRegion_gpu; parametersOnDevice.fParametersPerRegion_gpu = nullptr; ADEPT_DEVICE_API_CALL(MemcpyToSymbol(g4HepEmPars, ¶metersOnDevice, sizeof(G4HepEmParameters))); // Create G4HepEmData with the device pointers. G4HepEmData dataOnDevice; - dataOnDevice.fTheMatCutData = state->fData->fTheMatCutData_gpu; - dataOnDevice.fTheMaterialData = state->fData->fTheMaterialData_gpu; - dataOnDevice.fTheElementData = state->fData->fTheElementData_gpu; - dataOnDevice.fTheElectronData = state->fData->fTheElectronData_gpu; - dataOnDevice.fThePositronData = state->fData->fThePositronData_gpu; - dataOnDevice.fTheSBTableData = state->fData->fTheSBTableData_gpu; - dataOnDevice.fTheGammaData = state->fData->fTheGammaData_gpu; + dataOnDevice.fTheMatCutData = hepEmData->fTheMatCutData_gpu; + dataOnDevice.fTheMaterialData = hepEmData->fTheMaterialData_gpu; + dataOnDevice.fTheElementData = hepEmData->fTheElementData_gpu; + dataOnDevice.fTheElectronData = hepEmData->fTheElectronData_gpu; + dataOnDevice.fThePositronData = hepEmData->fThePositronData_gpu; + dataOnDevice.fTheSBTableData = hepEmData->fTheSBTableData_gpu; + dataOnDevice.fTheGammaData = hepEmData->fTheGammaData_gpu; // The other pointers should never be used. dataOnDevice.fTheMatCutData_gpu = nullptr; dataOnDevice.fTheMaterialData_gpu = nullptr; @@ -588,8 +569,6 @@ G4HepEmState *InitG4HepEm(G4HepEmConfig *hepEmConfig) dataOnDevice.fTheGammaData_gpu = nullptr; ADEPT_DEVICE_API_CALL(MemcpyToSymbol(g4HepEmData, &dataOnDevice, sizeof(G4HepEmData))); - - return state; } template @@ -1828,8 +1807,8 @@ std::thread LaunchGPUWorker(int trackCapacity, int leakCapacity, int scoringCapa hasWDTRegions}; } -void FreeGPU(std::unique_ptr &gpuState, G4HepEmState &g4hepem_state, - std::thread &gpuWorker, adeptint::WDTDeviceBuffers &wdtDev) +void FreeGPU(std::unique_ptr &gpuState, std::thread &gpuWorker, + adeptint::WDTDeviceBuffers &wdtDev) { gpuState->runTransport = false; gpuWorker.join(); @@ -1844,9 +1823,10 @@ void FreeGPU(std::unique_ptr // Free resources. gpuState.reset(); - // Free G4HepEm data - FreeG4HepEmData(g4hepem_state.fData); - FreeG4HepEmParametersOnGPU(g4hepem_state.fParameters); + // Note: the GPU mirror of G4HepEmParameters is not released here. + // That cleanup happens in HepEmHostData::~HepEmHostData(), which owns the + // upload lifecycle for the borrowed parameter block and performs the CUDA call + // via FreeG4HepEmParametersOnGPU() when the transport-owned HepEmHostData dies. // Free magnetic field #ifdef ADEPT_USE_EXT_BFIELD diff --git a/include/AdePT/core/AsyncAdePTTransport.hh b/include/AdePT/core/AsyncAdePTTransport.hh index 51ecbc60..9ce2904b 100644 --- a/include/AdePT/core/AsyncAdePTTransport.hh +++ b/include/AdePT/core/AsyncAdePTTransport.hh @@ -10,10 +10,10 @@ #define ASYNC_ADEPT_TRANSPORT_HH #include +#include #include #include #include -#include #include #include // forward declares vecgeom::cxx::VPlacedVolume @@ -27,8 +27,6 @@ class G4Region; class G4VPhysicalVolume; -class G4HepEmConfig; -struct G4HepEmState; namespace AsyncAdePT { struct TrackBuffer; struct GPUstate; @@ -51,7 +49,7 @@ private: unsigned short fMaxWDTIter{5}; ///< Maximum number of Woodcock tracking iterations per step std::unique_ptr fGPUstate{nullptr}; ///< CUDA state placeholder std::unique_ptr fBuffer{nullptr}; ///< Buffers for transferring tracks between host and device - std::unique_ptr fg4hepem_state; ///< The HepEm state singleton + std::unique_ptr fHepEmHostData; ///< Host-side HepEm data and parameter view adeptint::WDTDeviceBuffers fWDTDev{}; ///< device buffers for Woodcock tracking data std::thread fGPUWorker; ///< Thread to manage GPU std::condition_variable fCV_G4Workers; ///< Communicate with G4 workers @@ -74,14 +72,17 @@ private: ///< Needed to stall the GPU, in case the nPartInFlight * fHitBufferSafetyFactor > available HitSlots double fHitBufferSafetyFactor{1.5}; - void Initialize(G4HepEmConfig *hepEmConfig); + void Initialize(adeptint::VolAuxData *auxData, const adeptint::WDTHostPacked &wdtPacked, + const std::vector &uniformFieldValues); void InitBVH(); bool InitializeGeometry(const vecgeom::cxx::VPlacedVolume *world); - bool InitializePhysics(G4HepEmConfig *hepEmConfig); + bool InitializePhysics(); void InitWDTOnDevice(const adeptint::WDTHostPacked &src, adeptint::WDTDeviceBuffers &dev, unsigned short maxIter); public: - AsyncAdePTTransport(AdePTConfiguration &configuration, G4HepEmConfig *hepEmConfig); + AsyncAdePTTransport(AdePTConfiguration &configuration, std::unique_ptr hepEmHostData, + adeptint::VolAuxData *auxData, const adeptint::WDTHostPacked &wdtPacked, + const std::vector &uniformFieldValues); AsyncAdePTTransport(const AsyncAdePTTransport &other) = delete; ~AsyncAdePTTransport(); @@ -93,9 +94,6 @@ public: bool GetCallUserActions() const { return fReturnFirstAndLastStep; } std::vector const *GetGPURegionNames() { return fGPURegionNames; } std::vector const *GetCPURegionNames() { return fCPURegionNames; } - G4HepEmState *GetHepEmState() const { return fg4hepem_state.get(); } - void CompleteInitialization(adeptint::VolAuxData *auxData, const adeptint::WDTHostPacked &wdtPacked, - const std::vector &uniformFieldValues); /// Block until transport of the given event is done. void Flush(int threadId, int eventId, AdePTGeant4Integration &g4Integration); void ProcessGPUSteps(int threadId, int eventId, AdePTGeant4Integration &g4Integration); diff --git a/include/AdePT/core/AsyncAdePTTransport.icc b/include/AdePT/core/AsyncAdePTTransport.icc index 2967c224..c4b0c822 100644 --- a/include/AdePT/core/AsyncAdePTTransport.icc +++ b/include/AdePT/core/AsyncAdePTTransport.icc @@ -21,10 +21,7 @@ #include #include -#include -#include #include -#include #include #include @@ -34,7 +31,7 @@ namespace async_adept_impl { void setDeviceLimits(int stackLimit = 0, int heapLimit = 0); void CopySurfaceModelToGPU(); void InitWDTOnDevice(const adeptint::WDTHostPacked &, adeptint::WDTDeviceBuffers &, unsigned short); -G4HepEmState *InitG4HepEm(G4HepEmConfig *hepEmConfig); +void UploadG4HepEmToGPU(G4HepEmData *hepEmData, G4HepEmParameters *hepEmParameters); std::shared_ptr> GetGPUHits(unsigned int, AsyncAdePT::GPUstate &); std::pair GetGPUHitsFromBuffer(unsigned int, unsigned int, AsyncAdePT::GPUstate &, bool &); void CloseGPUBuffer(unsigned int, AsyncAdePT::GPUstate &, GPUHit *, const bool); @@ -45,7 +42,7 @@ std::unique_ptr InitializeGPU int trackCapacity, int leakCapacity, int scoringCapacity, int numThreads, AsyncAdePT::TrackBuffer &trackBuffer, double CPUCapacityFactor, double CPUCopyFraction, std::string &generalBfieldFile, const std::vector &uniformBfieldValues); -void FreeGPU(std::unique_ptr &, G4HepEmState &, std::thread &, +void FreeGPU(std::unique_ptr &, std::thread &, adeptint::WDTDeviceBuffers &); } // namespace async_adept_impl @@ -72,16 +69,20 @@ std::ostream &operator<<(std::ostream &stream, TrackDataWithIDs const &track) // These definitions live in a header-included .icc file, so they must remain // inline to avoid multiple definitions across translation units. -inline AsyncAdePTTransport::AsyncAdePTTransport(AdePTConfiguration &configuration, G4HepEmConfig *hepEmConfig) +inline AsyncAdePTTransport::AsyncAdePTTransport(AdePTConfiguration &configuration, + std::unique_ptr hepEmHostData, + adeptint::VolAuxData *auxData, const adeptint::WDTHostPacked &wdtPacked, + const std::vector &uniformFieldValues) : fAdePTSeed{configuration.GetAdePTSeed()}, fNThread{(ushort)configuration.GetNumThreads()}, fTrackCapacity{(uint)(1024 * 1024 * configuration.GetMillionsOfTrackSlots())}, fLeakCapacity{(uint)(1024 * 1024 * configuration.GetMillionsOfLeakSlots())}, fScoringCapacity{(uint)(1024 * 1024 * configuration.GetMillionsOfHitSlots())}, fDebugLevel{configuration.GetVerbosity()}, fCUDAStackLimit{configuration.GetCUDAStackLimit()}, fCUDAHeapLimit{configuration.GetCUDAHeapLimit()}, fLastNParticlesOnCPU{configuration.GetLastNParticlesOnCPU()}, - fMaxWDTIter{configuration.GetMaxWDTIter()}, fEventStates(fNThread), fGPUNetEnergy(fNThread, 0.0), - fTrackInAllRegions{configuration.GetTrackInAllRegions()}, fGPURegionNames{configuration.GetGPURegionNames()}, - fCPURegionNames{configuration.GetCPURegionNames()}, fReturnAllSteps{configuration.GetCallUserSteppingAction()}, + fMaxWDTIter{configuration.GetMaxWDTIter()}, fHepEmHostData(std::move(hepEmHostData)), fEventStates(fNThread), + fGPUNetEnergy(fNThread, 0.0), fTrackInAllRegions{configuration.GetTrackInAllRegions()}, + fGPURegionNames{configuration.GetGPURegionNames()}, fCPURegionNames{configuration.GetCPURegionNames()}, + fReturnAllSteps{configuration.GetCallUserSteppingAction()}, fReturnFirstAndLastStep{configuration.GetCallUserTrackingAction() || configuration.GetCallUserSteppingAction()}, fBfieldFile{configuration.GetCovfieBfieldFile()}, fCPUCapacityFactor{configuration.GetCPUCapacityFactor()}, fCPUCopyFraction{configuration.GetHitBufferFlushThreshold()}, @@ -94,12 +95,12 @@ inline AsyncAdePTTransport::AsyncAdePTTransport(AdePTConfiguration &configuratio std::atomic_init(&eventState, EventState::LeakedTracksRetrieved); } - AsyncAdePTTransport::Initialize(hepEmConfig); + AsyncAdePTTransport::Initialize(auxData, wdtPacked, uniformFieldValues); } inline AsyncAdePTTransport::~AsyncAdePTTransport() { - async_adept_impl::FreeGPU(std::ref(fGPUstate), *fg4hepem_state, fGPUWorker, fWDTDev); + async_adept_impl::FreeGPU(std::ref(fGPUstate), fGPUWorker, fWDTDev); } inline void AsyncAdePTTransport::AddTrack(int pdg, uint64_t trackId, uint64_t parentId, double energy, double x, @@ -175,14 +176,19 @@ inline bool AsyncAdePTTransport::InitializeGeometry(const vecgeom::cxx::VPlacedV return success; } -inline bool AsyncAdePTTransport::InitializePhysics(G4HepEmConfig *hepEmConfig) +inline bool AsyncAdePTTransport::InitializePhysics() { - // Initialize shared physics data - fg4hepem_state.reset(async_adept_impl::InitG4HepEm(hepEmConfig)); + if (!fHepEmHostData) { + throw std::runtime_error("AsyncAdePTTransport::InitializePhysics: Missing HepEm host data."); + } + + // Upload the prepared HepEm data and parameters to the device. + async_adept_impl::UploadG4HepEmToGPU(fHepEmHostData->GetData(), fHepEmHostData->GetParameters()); return true; } -inline void AsyncAdePTTransport::Initialize(G4HepEmConfig *hepEmConfig) +inline void AsyncAdePTTransport::Initialize(adeptint::VolAuxData *auxData, const adeptint::WDTHostPacked &wdtPacked, + const std::vector &uniformFieldValues) { if (vecgeom::GeoManager::Instance().GetRegisteredVolumesCount() == 0) throw std::runtime_error("AsyncAdePTTransport::Initialize: Number of geometry volumes is zero."); @@ -196,23 +202,12 @@ inline void AsyncAdePTTransport::Initialize(G4HepEmConfig *hepEmConfig) if (!InitializeGeometry(world)) throw std::runtime_error("AsyncAdePTTransport::Initialize: Cannot initialize geometry on GPU"); - // Initialize G4HepEm - if (!InitializePhysics(hepEmConfig)) + // Upload the prepared HepEm physics data to the device. + if (!InitializePhysics()) throw std::runtime_error("AsyncAdePTTransport::Initialize cannot initialize physics on GPU"); -} - -inline void AsyncAdePTTransport::CompleteInitialization(adeptint::VolAuxData *auxData, - const adeptint::WDTHostPacked &wdtPacked, - const std::vector &uniformFieldValues) -{ - // This is the second half of the split initialization. A non-zero volume count was already - // required in Initialize() before geometry upload, and it remains a hard precondition here - // before uploading any geometry-derived metadata to the device. - const auto numVolumes = vecgeom::GeoManager::Instance().GetRegisteredVolumesCount(); - if (numVolumes == 0) - throw std::runtime_error("AsyncAdePTTransport::CompleteInitialization: Number of geometry volumes is zero."); // Initialize volume auxiliary data on device + const auto numVolumes = vecgeom::GeoManager::Instance().GetRegisteredVolumesCount(); auto &volAuxArray = adeptint::VolAuxArray::GetInstance(); volAuxArray.fNumVolumes = numVolumes; volAuxArray.fAuxData = auxData; diff --git a/include/AdePT/integration/AdePTGeant4Integration.hh b/include/AdePT/integration/AdePTGeant4Integration.hh index 46660cec..0d3466eb 100644 --- a/include/AdePT/integration/AdePTGeant4Integration.hh +++ b/include/AdePT/integration/AdePTGeant4Integration.hh @@ -20,7 +20,7 @@ #include #include -struct G4HepEmState; +struct G4HepEmData; namespace AdePTGeant4Integration_detail { struct ScoringObjects; @@ -53,10 +53,10 @@ public: static void CreateVecGeomWorld(G4VPhysicalVolume const *physvol); /// @brief This function compares G4 and VecGeom geometries and reports any differences - static void CheckGeometry(G4HepEmState *hepEmState); + static void CheckGeometry(G4HepEmData const *hepEmData); /// @brief Fills the auxiliary data needed for AdePT - static void InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmState *hepEmState, + static void InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmData const *hepEmData, G4HepEmTrackingManagerSpecialized *hepEmTM, bool trackInAllRegions, std::vector const *gpuRegionNames, adeptint::WDTHostRaw &wdtRaw); diff --git a/src/AdePTGeant4Integration.cpp b/src/AdePTGeant4Integration.cpp index 8b13fca9..0d15e8e3 100644 --- a/src/AdePTGeant4Integration.cpp +++ b/src/AdePTGeant4Integration.cpp @@ -23,7 +23,6 @@ #include #include -#include #include #include #include @@ -257,7 +256,7 @@ namespace { struct VisitContext { const int *g4tohepmcindex; std::size_t nvolumes; - G4HepEmState const *hepEmState; + G4HepEmData const *hepEmData; }; /// Recursive geometry visitor matching one by one Geant4 and VecGeom logical volumes @@ -329,7 +328,7 @@ void visitGeometry(G4VPhysicalVolume const *g4_pvol, vecgeom::VPlacedVolume cons const int g4mcindex = g4_lvol->GetMaterialCutsCouple()->GetIndex(); const int hepemmcindex = context.g4tohepmcindex[g4mcindex]; // Check consistency with G4HepEm data - if (context.hepEmState->fData->fTheMatCutData->fMatCutData[hepemmcindex].fG4MatCutIndex != g4mcindex) + if (context.hepEmData->fTheMatCutData->fMatCutData[hepemmcindex].fG4MatCutIndex != g4mcindex) throw std::runtime_error("Fatal: CheckGeometry: Mismatch between Geant4 mcindex and corresponding G4HepEm index"); if (vg_lvol->id() >= context.nvolumes) throw std::runtime_error("Fatal: CheckGeometry: Volume id larger than number of volumes"); @@ -344,21 +343,21 @@ void visitGeometry(G4VPhysicalVolume const *g4_pvol, vecgeom::VPlacedVolume cons } } // namespace -void AdePTGeant4Integration::CheckGeometry(G4HepEmState *hepEmState) +void AdePTGeant4Integration::CheckGeometry(G4HepEmData const *hepEmData) { const G4VPhysicalVolume *g4world = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume(); const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); - const int *g4tohepmcindex = hepEmState->fData->fTheMatCutData->fG4MCIndexToHepEmMCIndex; + const int *g4tohepmcindex = hepEmData->fTheMatCutData->fG4MCIndexToHepEmMCIndex; const auto nvolumes = vecgeom::GeoManager::Instance().GetRegisteredVolumesCount(); std::cout << "Visiting geometry ...\n"; - const VisitContext context{g4tohepmcindex, nvolumes, hepEmState}; + const VisitContext context{g4tohepmcindex, nvolumes, hepEmData}; visitGeometry(g4world, vecgeomWorld, context); std::cout << "Visiting geometry done\n"; } -void AdePTGeant4Integration::InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmState *hepEmState, +void AdePTGeant4Integration::InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmData const *hepEmData, G4HepEmTrackingManagerSpecialized *hepEmTM, bool trackInAllRegions, std::vector const *gpuRegionNames, adeptint::WDTHostRaw &wdtRaw) @@ -371,7 +370,7 @@ void AdePTGeant4Integration::InitVolAuxData(adeptint::VolAuxData *volAuxData, G4 const G4VPhysicalVolume *g4world = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume(); const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); - const int *g4tohepmcindex = hepEmState->fData->fTheMatCutData->fG4MCIndexToHepEmMCIndex; + const int *g4tohepmcindex = hepEmData->fTheMatCutData->fG4MCIndexToHepEmMCIndex; // We need to go from region names to G4Region std::vector gpuRegions{}; diff --git a/src/AdePTHepEmHostData.cpp b/src/AdePTHepEmHostData.cpp new file mode 100644 index 00000000..83766822 --- /dev/null +++ b/src/AdePTHepEmHostData.cpp @@ -0,0 +1,54 @@ +// SPDX-FileCopyrightText: 2026 CERN +// SPDX-License-Identifier: Apache-2.0 + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace AsyncAdePT { + +/// @brief Release the tables owned by `G4HepEmData` and then delete the outer object. +void HepEmHostData::DataDeleter::operator()(G4HepEmData *data) const +{ + if (data == nullptr) return; + FreeG4HepEmData(data); + delete data; +} + +/// @brief Rebuild a complete host-side HepEm view from the supplied Geant4 HepEm config. +HepEmHostData::HepEmHostData(G4HepEmConfig *hepEmConfig) + : fData(new G4HepEmData), fParameters(hepEmConfig->GetG4HepEmParameters()) +{ + // Rebuild the HepEm tables from the Geant4-owned config so the host-side + // transport preparation has a complete, self-contained view of the data. + InitG4HepEmData(fData.get()); + InitMaterialAndCoupleData(fData.get(), fParameters); + + // Build all EM species up front so the subsequent host preparation can run + // once before the GPU-side upload/initialization. + InitElectronData(fData.get(), fParameters, true); + InitElectronData(fData.get(), fParameters, false); + InitGammaData(fData.get(), fParameters); + + G4HepEmMatCutData *cutData = fData->fTheMatCutData; + G4cout << "fNumG4MatCuts = " << cutData->fNumG4MatCuts << ", fNumMatCutData = " << cutData->fNumMatCutData << G4endl; +} + +/// @brief Release the GPU mirror of the borrowed HepEm parameters. +/// @details +/// The host-side `G4HepEmParameters` remain owned by the `G4HepEmConfig`. +/// AdePT only owns the device allocation created by `CopyG4HepEmParametersToGPU`, +/// so this destructor releases just that GPU-side mirror. +HepEmHostData::~HepEmHostData() +{ + FreeG4HepEmParametersOnGPU(fParameters); +} + +} // namespace AsyncAdePT diff --git a/src/AdePTTrackingManager.cc b/src/AdePTTrackingManager.cc index b0a17c1e..21951378 100644 --- a/src/AdePTTrackingManager.cc +++ b/src/AdePTTrackingManager.cc @@ -37,7 +37,10 @@ std::weak_ptr &SharedAdePTTransportStorage() } std::shared_ptr CreateSharedAdePTTransport(AdePTConfiguration &conf, - G4HepEmTrackingManagerSpecialized *hepEmTM) + std::unique_ptr hepEmHostData, + adeptint::VolAuxData *auxData, + const adeptint::WDTHostPacked &wdtPacked, + const std::vector &uniformFieldValues) { auto &transport = SharedAdePTTransportStorage(); // weak_ptr::lock() promotes the stored weak reference to a shared_ptr if the @@ -46,8 +49,9 @@ std::shared_ptr CreateSharedAdePTTransport(AdePTConfiguration &c return existing; } - auto created = std::make_shared(conf, hepEmTM->GetConfig()); - transport = created; + auto created = + std::make_shared(conf, std::move(hepEmHostData), auxData, wdtPacked, uniformFieldValues); + transport = created; return created; } @@ -91,23 +95,22 @@ void AdePTTrackingManager::InitializeSharedAdePTTransport() if (fAdePTConfiguration->GetCovfieBfieldFile() == "") std::cout << "No magnetic field file provided!" << std::endl; #endif const auto uniformFieldValues = fGeant4Integration.GetUniformField(); - - // Create the shared AdePT transport engine on the first worker thread. - fAdeptTransport = CreateSharedAdePTTransport(*fAdePTConfiguration, fHepEmTrackingManager.get()); + auto hepEmHostData = std::make_unique(fHepEmTrackingManager->GetConfig()); // Check VecGeom geometry matches Geant4 before deriving any geometry metadata for transport. - fGeant4Integration.CheckGeometry(fAdeptTransport->GetHepEmState()); + fGeant4Integration.CheckGeometry(hepEmHostData->GetData()); // Initialize auxiliary per-LV data and collect the raw WDT metadata on the Geant4 side. auto *auxData = new adeptint::VolAuxData[vecgeom::GeoManager::Instance().GetRegisteredVolumesCount()]; adeptint::WDTHostRaw wdtRaw; - fGeant4Integration.InitVolAuxData(auxData, fAdeptTransport->GetHepEmState(), fHepEmTrackingManager.get(), + fGeant4Integration.InitVolAuxData(auxData, hepEmHostData->GetData(), fHepEmTrackingManager.get(), fAdePTConfiguration->GetTrackInAllRegions(), fAdePTConfiguration->GetGPURegionNames(), wdtRaw); adeptint::WDTHostPacked wdtPacked = adeptint::PackWDT(wdtRaw); - // Finish the shared transport initialization by uploading the prepared metadata to the device. - fAdeptTransport->CompleteInitialization(auxData, wdtPacked, uniformFieldValues); + // Create the shared AdePT transport engine once the full host-side preparation is complete. + fAdeptTransport = CreateSharedAdePTTransport(*fAdePTConfiguration, std::move(hepEmHostData), auxData, wdtPacked, + uniformFieldValues); } void AdePTTrackingManager::InitializeAdePT() From a7f506e5458a0a430cf929c1878e239773537c4f Mon Sep 17 00:00:00 2001 From: SeverinDiederichs Date: Wed, 18 Mar 2026 16:53:43 +0100 Subject: [PATCH 2/6] cleaning and improved doc --- include/AdePT/core/AdePTHepEmHostData.hh | 19 +++++++++++++++---- include/AdePT/core/AsyncAdePTTransport.cuh | 2 +- src/AdePTHepEmHostData.cpp | 14 +++++++++----- 3 files changed, 25 insertions(+), 10 deletions(-) diff --git a/include/AdePT/core/AdePTHepEmHostData.hh b/include/AdePT/core/AdePTHepEmHostData.hh index 1d164763..99ddaa50 100644 --- a/include/AdePT/core/AdePTHepEmHostData.hh +++ b/include/AdePT/core/AdePTHepEmHostData.hh @@ -22,13 +22,22 @@ namespace AsyncAdePT { /// allocation is released here in the destructor, because this helper owns the /// corresponding upload lifecycle even though the host parameters themselves are /// still owned by G4HepEm. +/// +/// Cleanup is intentionally split: +/// - `DataDeleter` performs the deep cleanup of the owned `G4HepEmData` +/// and then deletes the outer `G4HepEmData` allocation. +/// - `~HepEmHostData()` releases only the GPU mirror of the borrowed +/// `G4HepEmParameters`. class HepEmHostData { public: /// @brief Rebuild all host-side HepEm tables needed by AdePT from the given config. explicit HepEmHostData(G4HepEmConfig *hepEmConfig); /// @brief Release the GPU mirror of the borrowed HepEm parameters. - /// @details This performs a CUDA-side free through `FreeG4HepEmParametersOnGPU`. + /// @details + /// This performs a CUDA-side free through `FreeG4HepEmParametersOnGPU`. + /// It does not free the host-side `G4HepEmParameters`, because those remain + /// owned by the original `G4HepEmConfig`. ~HepEmHostData(); HepEmHostData(const HepEmHostData &) = delete; @@ -45,9 +54,11 @@ public: private: /// @brief Deletes the outer `G4HepEmData` object after first freeing all tables it owns. - /// @details `FreeG4HepEmData` releases both the host-side tables and any device-side + /// @details + /// `FreeG4HepEmData` releases both the host-side tables and any device-side /// mirrors embedded in the `G4HepEmData` object, but it does not delete the outer - /// `G4HepEmData` allocation itself. This deleter performs both steps. + /// `G4HepEmData` allocation itself. This deleter performs both steps for the + /// owned `fData` member. It does not touch the separately borrowed parameter block. struct DataDeleter { void operator()(G4HepEmData *data) const; }; @@ -55,7 +66,7 @@ private: /// Owned host-side HepEm tables rebuilt for AdePT. std::unique_ptr fData; - /// Non-owning pointer to the Geant4-owned HepEm parameter block. + /// Non-owning pointer to the G4HepEmParameters that are owned by G4HepEm itself, not AdePT. G4HepEmParameters *fParameters{nullptr}; }; diff --git a/include/AdePT/core/AsyncAdePTTransport.cuh b/include/AdePT/core/AsyncAdePTTransport.cuh index 5c749b23..44ce8276 100644 --- a/include/AdePT/core/AsyncAdePTTransport.cuh +++ b/include/AdePT/core/AsyncAdePTTransport.cuh @@ -536,7 +536,7 @@ void CopySurfaceModelToGPU() void UploadG4HepEmToGPU(G4HepEmData *hepEmData, G4HepEmParameters *hepEmParameters) { if (hepEmData == nullptr || hepEmParameters == nullptr) { - throw std::runtime_error("UploadG4HepEmToGPU requires non-null HepEm data and parameters."); + throw std::runtime_error("UploadG4HepEmToGPU requires non-null G4HepEmData and G4HepEmParameters."); } // Copy the prepared host-side HepEm data to the GPU. diff --git a/src/AdePTHepEmHostData.cpp b/src/AdePTHepEmHostData.cpp index 83766822..c3cf8f00 100644 --- a/src/AdePTHepEmHostData.cpp +++ b/src/AdePTHepEmHostData.cpp @@ -15,6 +15,10 @@ namespace AsyncAdePT { /// @brief Release the tables owned by `G4HepEmData` and then delete the outer object. +/// @details +/// This is the cleanup path for the fully owned `fData` member. It is separate +/// from the class destructor because `std::unique_ptr` needs a deleter for the +/// deep cleanup before the outer `G4HepEmData` allocation itself can be deleted. void HepEmHostData::DataDeleter::operator()(G4HepEmData *data) const { if (data == nullptr) return; @@ -22,17 +26,16 @@ void HepEmHostData::DataDeleter::operator()(G4HepEmData *data) const delete data; } -/// @brief Rebuild a complete host-side HepEm view from the supplied Geant4 HepEm config. +/// @brief Rebuild a complete host-side G4HepEm view from the supplied G4HepEm config. HepEmHostData::HepEmHostData(G4HepEmConfig *hepEmConfig) : fData(new G4HepEmData), fParameters(hepEmConfig->GetG4HepEmParameters()) { - // Rebuild the HepEm tables from the Geant4-owned config so the host-side + // Rebuild the HepEm tables from the G4HepEm-owned config so the host-side // transport preparation has a complete, self-contained view of the data. InitG4HepEmData(fData.get()); InitMaterialAndCoupleData(fData.get(), fParameters); - // Build all EM species up front so the subsequent host preparation can run - // once before the GPU-side upload/initialization. + // Build all EM species InitElectronData(fData.get(), fParameters, true); InitElectronData(fData.get(), fParameters, false); InitGammaData(fData.get(), fParameters); @@ -45,7 +48,8 @@ HepEmHostData::HepEmHostData(G4HepEmConfig *hepEmConfig) /// @details /// The host-side `G4HepEmParameters` remain owned by the `G4HepEmConfig`. /// AdePT only owns the device allocation created by `CopyG4HepEmParametersToGPU`, -/// so this destructor releases just that GPU-side mirror. +/// so this destructor releases just that GPU-side mirror. The owned `G4HepEmData` +/// cleanup is handled independently by `DataDeleter`. HepEmHostData::~HepEmHostData() { FreeG4HepEmParametersOnGPU(fParameters); From b605b08cb4fd35a30a7bd734cb3f0b93babb2e98 Mon Sep 17 00:00:00 2001 From: SeverinDiederichs Date: Wed, 18 Mar 2026 17:26:43 +0100 Subject: [PATCH 3/6] cleanup --- src/AdePTTrackingManager.cc | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/AdePTTrackingManager.cc b/src/AdePTTrackingManager.cc index 21951378..d75e3d53 100644 --- a/src/AdePTTrackingManager.cc +++ b/src/AdePTTrackingManager.cc @@ -36,11 +36,11 @@ std::weak_ptr &SharedAdePTTransportStorage() return transport; } -std::shared_ptr CreateSharedAdePTTransport(AdePTConfiguration &conf, - std::unique_ptr hepEmHostData, - adeptint::VolAuxData *auxData, - const adeptint::WDTHostPacked &wdtPacked, - const std::vector &uniformFieldValues) +std::shared_ptr GetSharedAdePTTransport(AdePTConfiguration &conf, + std::unique_ptr hepEmHostData, + adeptint::VolAuxData *auxData, + const adeptint::WDTHostPacked &wdtPacked, + const std::vector &uniformFieldValues) { auto &transport = SharedAdePTTransportStorage(); // weak_ptr::lock() promotes the stored weak reference to a shared_ptr if the @@ -109,8 +109,8 @@ void AdePTTrackingManager::InitializeSharedAdePTTransport() adeptint::WDTHostPacked wdtPacked = adeptint::PackWDT(wdtRaw); // Create the shared AdePT transport engine once the full host-side preparation is complete. - fAdeptTransport = CreateSharedAdePTTransport(*fAdePTConfiguration, std::move(hepEmHostData), auxData, wdtPacked, - uniformFieldValues); + fAdeptTransport = + GetSharedAdePTTransport(*fAdePTConfiguration, std::move(hepEmHostData), auxData, wdtPacked, uniformFieldValues); } void AdePTTrackingManager::InitializeAdePT() From de275bde7fa4b1789e2a91f937fa2626532a78f4 Mon Sep 17 00:00:00 2001 From: SeverinDiederichs Date: Sat, 21 Mar 2026 13:44:24 +0100 Subject: [PATCH 4/6] fix ownership of G4HepEm parameters, do a deep copy and do not borrow it from worker, as data might be deleted before the AdePTTransport is --- CMakeLists.txt | 2 +- include/AdePT/core/AdePTG4HepEmState.hh | 80 +++++++++++++++++ include/AdePT/core/AdePTHepEmHostData.hh | 75 ---------------- include/AdePT/core/AsyncAdePTTransport.cuh | 8 +- include/AdePT/core/AsyncAdePTTransport.hh | 13 +-- include/AdePT/core/AsyncAdePTTransport.icc | 15 ++-- src/AdePTG4HepEmState.cpp | 100 +++++++++++++++++++++ src/AdePTHepEmHostData.cpp | 58 ------------ src/AdePTTrackingManager.cc | 24 ++--- 9 files changed, 212 insertions(+), 163 deletions(-) create mode 100644 include/AdePT/core/AdePTG4HepEmState.hh delete mode 100644 include/AdePT/core/AdePTHepEmHostData.hh create mode 100644 src/AdePTG4HepEmState.cpp delete mode 100644 src/AdePTHepEmHostData.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 51297c07..f238296a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -444,7 +444,7 @@ endif() # Build Targets #----------------------------------------------------------------------------# set(ADEPT_G4_INTEGRATION_SRCS - src/AdePTHepEmHostData.cpp + src/AdePTG4HepEmState.cpp src/G4HepEmTrackingManagerSpecialized.cc src/AdePTTrackingManager.cc src/G4EmStandardPhysics_AdePT.cc diff --git a/include/AdePT/core/AdePTG4HepEmState.hh b/include/AdePT/core/AdePTG4HepEmState.hh new file mode 100644 index 00000000..a396c67e --- /dev/null +++ b/include/AdePT/core/AdePTG4HepEmState.hh @@ -0,0 +1,80 @@ +// SPDX-FileCopyrightText: 2026 CERN +// SPDX-License-Identifier: Apache-2.0 + +#ifndef ADEPT_G4_HEPEM_STATE_HH +#define ADEPT_G4_HEPEM_STATE_HH + +#include + +struct G4HepEmConfig; +struct G4HepEmData; +struct G4HepEmParameters; + +namespace AsyncAdePT { + +/// @brief Owns the prepared host-side G4HepEm inputs used by transport. +/// @details +/// The Geant4 integration side prepares one of these objects before the shared +/// transport is created. This wrapper owns both: +/// - the rebuilt `G4HepEmData` +/// - a deep copy of the `G4HepEmParameters` taken from the provided config +/// +/// Cleanup is intentionally split: +/// - `DataDeleter` performs the deep cleanup of the owned `G4HepEmData` +/// and then deletes the outer `G4HepEmData` allocation. +/// - `ParametersDeleter` performs the deep cleanup of the owned +/// `G4HepEmParameters`, including the GPU mirror created by AdePT, and then +/// deletes the outer `G4HepEmParameters` allocation. +class AdePTG4HepEmState { +public: + /// @brief Build the AdePT-owned `G4HepEmData` and `G4HepEmParameters` copies from the supplied config. + explicit AdePTG4HepEmState(G4HepEmConfig *hepEmConfig); + + /// @brief Destroy the owned `G4HepEmData` and `G4HepEmParameters` copies. + ~AdePTG4HepEmState(); + + AdePTG4HepEmState(const AdePTG4HepEmState &) = delete; + AdePTG4HepEmState &operator=(const AdePTG4HepEmState &) = delete; + + AdePTG4HepEmState(AdePTG4HepEmState &&) noexcept = default; + AdePTG4HepEmState &operator=(AdePTG4HepEmState &&) noexcept = default; + + /// @brief Access the owned host-side HepEm data tables. + G4HepEmData *GetData() const { return fData.get(); } + + /// @brief Access the owned HepEm parameter copy. + G4HepEmParameters *GetParameters() const { return fParameters.get(); } + +private: + /// @brief Deletes the outer `G4HepEmData` object after first freeing all tables it owns. + /// @details + /// `FreeG4HepEmData` releases both the host-side tables and any device-side + /// mirrors embedded in the `G4HepEmData` object, but it does not delete the outer + /// `G4HepEmData` allocation itself. This deleter performs both steps for the + /// owned `fData` member. It does not touch the separately owned + /// `G4HepEmParameters` copy stored in `fParameters`. + struct DataDeleter { + void operator()(G4HepEmData *data) const; + }; + + /// @brief Deletes the outer `G4HepEmParameters` object after first freeing + /// all host/device allocations it owns. + /// @details + /// The copied parameter block owns its `fParametersPerRegion` host array and + /// the GPU mirror pointed to by `fParametersPerRegion_gpu` after transport + /// upload. `FreeG4HepEmParameters` releases those nested allocations, while + /// this deleter also deletes the outer `G4HepEmParameters` allocation. + struct ParametersDeleter { + void operator()(G4HepEmParameters *parameters) const; + }; + + /// Owned `G4HepEmData` rebuilt for AdePT. + std::unique_ptr fData; + + /// Owned deep copy of `G4HepEmParameters` used to build and upload transport data. + std::unique_ptr fParameters; +}; + +} // namespace AsyncAdePT + +#endif diff --git a/include/AdePT/core/AdePTHepEmHostData.hh b/include/AdePT/core/AdePTHepEmHostData.hh deleted file mode 100644 index 99ddaa50..00000000 --- a/include/AdePT/core/AdePTHepEmHostData.hh +++ /dev/null @@ -1,75 +0,0 @@ -// SPDX-FileCopyrightText: 2026 CERN -// SPDX-License-Identifier: Apache-2.0 - -#ifndef ADEPT_HEPEM_HOST_DATA_HH -#define ADEPT_HEPEM_HOST_DATA_HH - -#include - -struct G4HepEmConfig; -struct G4HepEmData; -struct G4HepEmParameters; - -namespace AsyncAdePT { - -/// @brief Owns the host-side HepEm tables rebuilt from a `G4HepEmConfig`. -/// @details -/// The Geant4 integration side prepares one of these objects before the shared -/// transport is created. The helper owns the rebuilt `G4HepEmData`, while the -/// `G4HepEmParameters` pointer is only borrowed from the `G4HepEmConfig`. -/// -/// AdePT also creates the GPU mirror for the borrowed parameters. That GPU-side -/// allocation is released here in the destructor, because this helper owns the -/// corresponding upload lifecycle even though the host parameters themselves are -/// still owned by G4HepEm. -/// -/// Cleanup is intentionally split: -/// - `DataDeleter` performs the deep cleanup of the owned `G4HepEmData` -/// and then deletes the outer `G4HepEmData` allocation. -/// - `~HepEmHostData()` releases only the GPU mirror of the borrowed -/// `G4HepEmParameters`. -class HepEmHostData { -public: - /// @brief Rebuild all host-side HepEm tables needed by AdePT from the given config. - explicit HepEmHostData(G4HepEmConfig *hepEmConfig); - - /// @brief Release the GPU mirror of the borrowed HepEm parameters. - /// @details - /// This performs a CUDA-side free through `FreeG4HepEmParametersOnGPU`. - /// It does not free the host-side `G4HepEmParameters`, because those remain - /// owned by the original `G4HepEmConfig`. - ~HepEmHostData(); - - HepEmHostData(const HepEmHostData &) = delete; - HepEmHostData &operator=(const HepEmHostData &) = delete; - - HepEmHostData(HepEmHostData &&) noexcept = default; - HepEmHostData &operator=(HepEmHostData &&) noexcept = default; - - /// @brief Access the owned host-side HepEm data tables. - G4HepEmData *GetData() const { return fData.get(); } - - /// @brief Access the borrowed HepEm parameter block from the Geant4 config. - G4HepEmParameters *GetParameters() const { return fParameters; } - -private: - /// @brief Deletes the outer `G4HepEmData` object after first freeing all tables it owns. - /// @details - /// `FreeG4HepEmData` releases both the host-side tables and any device-side - /// mirrors embedded in the `G4HepEmData` object, but it does not delete the outer - /// `G4HepEmData` allocation itself. This deleter performs both steps for the - /// owned `fData` member. It does not touch the separately borrowed parameter block. - struct DataDeleter { - void operator()(G4HepEmData *data) const; - }; - - /// Owned host-side HepEm tables rebuilt for AdePT. - std::unique_ptr fData; - - /// Non-owning pointer to the G4HepEmParameters that are owned by G4HepEm itself, not AdePT. - G4HepEmParameters *fParameters{nullptr}; -}; - -} // namespace AsyncAdePT - -#endif diff --git a/include/AdePT/core/AsyncAdePTTransport.cuh b/include/AdePT/core/AsyncAdePTTransport.cuh index 44ce8276..b27a3d0f 100644 --- a/include/AdePT/core/AsyncAdePTTransport.cuh +++ b/include/AdePT/core/AsyncAdePTTransport.cuh @@ -1823,10 +1823,10 @@ void FreeGPU(std::unique_ptr // Free resources. gpuState.reset(); - // Note: the GPU mirror of G4HepEmParameters is not released here. - // That cleanup happens in HepEmHostData::~HepEmHostData(), which owns the - // upload lifecycle for the borrowed parameter block and performs the CUDA call - // via FreeG4HepEmParametersOnGPU() when the transport-owned HepEmHostData dies. + // Note: the GPU mirror of `G4HepEmParameters` is not released here. + // That cleanup happens when the transport-owned `AdePTG4HepEmState` dies, + // because it owns both the copied `G4HepEmParameters` object and the upload + // lifecycle attached to that copy. // Free magnetic field #ifdef ADEPT_USE_EXT_BFIELD diff --git a/include/AdePT/core/AsyncAdePTTransport.hh b/include/AdePT/core/AsyncAdePTTransport.hh index 9ce2904b..095cf1d3 100644 --- a/include/AdePT/core/AsyncAdePTTransport.hh +++ b/include/AdePT/core/AsyncAdePTTransport.hh @@ -10,7 +10,7 @@ #define ASYNC_ADEPT_TRANSPORT_HH #include -#include +#include #include #include #include @@ -48,10 +48,11 @@ private: unsigned short fLastNParticlesOnCPU{0}; ///< Number N of last N particles that are finished on CPU unsigned short fMaxWDTIter{5}; ///< Maximum number of Woodcock tracking iterations per step std::unique_ptr fGPUstate{nullptr}; ///< CUDA state placeholder - std::unique_ptr fBuffer{nullptr}; ///< Buffers for transferring tracks between host and device - std::unique_ptr fHepEmHostData; ///< Host-side HepEm data and parameter view - adeptint::WDTDeviceBuffers fWDTDev{}; ///< device buffers for Woodcock tracking data - std::thread fGPUWorker; ///< Thread to manage GPU + std::unique_ptr fBuffer{nullptr}; ///< Buffers for transferring tracks between host and device + std::unique_ptr + fAdePTG4HepEmState; ///< Transport-owned wrapper around `G4HepEmData` and copied `G4HepEmParameters` + adeptint::WDTDeviceBuffers fWDTDev{}; ///< device buffers for Woodcock tracking data + std::thread fGPUWorker; ///< Thread to manage GPU std::condition_variable fCV_G4Workers; ///< Communicate with G4 workers std::mutex fMutex_G4Workers; ///< Mutex associated to the condition variable std::vector> fEventStates; ///< State machine for each G4 worker @@ -80,7 +81,7 @@ private: void InitWDTOnDevice(const adeptint::WDTHostPacked &src, adeptint::WDTDeviceBuffers &dev, unsigned short maxIter); public: - AsyncAdePTTransport(AdePTConfiguration &configuration, std::unique_ptr hepEmHostData, + AsyncAdePTTransport(AdePTConfiguration &configuration, std::unique_ptr adeptG4HepEmState, adeptint::VolAuxData *auxData, const adeptint::WDTHostPacked &wdtPacked, const std::vector &uniformFieldValues); AsyncAdePTTransport(const AsyncAdePTTransport &other) = delete; diff --git a/include/AdePT/core/AsyncAdePTTransport.icc b/include/AdePT/core/AsyncAdePTTransport.icc index c4b0c822..53961aa0 100644 --- a/include/AdePT/core/AsyncAdePTTransport.icc +++ b/include/AdePT/core/AsyncAdePTTransport.icc @@ -70,7 +70,7 @@ std::ostream &operator<<(std::ostream &stream, TrackDataWithIDs const &track) // These definitions live in a header-included .icc file, so they must remain // inline to avoid multiple definitions across translation units. inline AsyncAdePTTransport::AsyncAdePTTransport(AdePTConfiguration &configuration, - std::unique_ptr hepEmHostData, + std::unique_ptr adeptG4HepEmState, adeptint::VolAuxData *auxData, const adeptint::WDTHostPacked &wdtPacked, const std::vector &uniformFieldValues) : fAdePTSeed{configuration.GetAdePTSeed()}, fNThread{(ushort)configuration.GetNumThreads()}, @@ -79,8 +79,8 @@ inline AsyncAdePTTransport::AsyncAdePTTransport(AdePTConfiguration &configuratio fScoringCapacity{(uint)(1024 * 1024 * configuration.GetMillionsOfHitSlots())}, fDebugLevel{configuration.GetVerbosity()}, fCUDAStackLimit{configuration.GetCUDAStackLimit()}, fCUDAHeapLimit{configuration.GetCUDAHeapLimit()}, fLastNParticlesOnCPU{configuration.GetLastNParticlesOnCPU()}, - fMaxWDTIter{configuration.GetMaxWDTIter()}, fHepEmHostData(std::move(hepEmHostData)), fEventStates(fNThread), - fGPUNetEnergy(fNThread, 0.0), fTrackInAllRegions{configuration.GetTrackInAllRegions()}, + fMaxWDTIter{configuration.GetMaxWDTIter()}, fAdePTG4HepEmState(std::move(adeptG4HepEmState)), + fEventStates(fNThread), fGPUNetEnergy(fNThread, 0.0), fTrackInAllRegions{configuration.GetTrackInAllRegions()}, fGPURegionNames{configuration.GetGPURegionNames()}, fCPURegionNames{configuration.GetCPURegionNames()}, fReturnAllSteps{configuration.GetCallUserSteppingAction()}, fReturnFirstAndLastStep{configuration.GetCallUserTrackingAction() || configuration.GetCallUserSteppingAction()}, @@ -178,12 +178,13 @@ inline bool AsyncAdePTTransport::InitializeGeometry(const vecgeom::cxx::VPlacedV inline bool AsyncAdePTTransport::InitializePhysics() { - if (!fHepEmHostData) { - throw std::runtime_error("AsyncAdePTTransport::InitializePhysics: Missing HepEm host data."); + if (!fAdePTG4HepEmState) { + throw std::runtime_error("AsyncAdePTTransport::InitializePhysics: Missing AdePT-owned G4HepEm state."); } - // Upload the prepared HepEm data and parameters to the device. - async_adept_impl::UploadG4HepEmToGPU(fHepEmHostData->GetData(), fHepEmHostData->GetParameters()); + // Upload the transport-owned `G4HepEmData` and copied + // `G4HepEmParameters` to the device. + async_adept_impl::UploadG4HepEmToGPU(fAdePTG4HepEmState->GetData(), fAdePTG4HepEmState->GetParameters()); return true; } diff --git a/src/AdePTG4HepEmState.cpp b/src/AdePTG4HepEmState.cpp new file mode 100644 index 00000000..87753306 --- /dev/null +++ b/src/AdePTG4HepEmState.cpp @@ -0,0 +1,100 @@ +// SPDX-FileCopyrightText: 2026 CERN +// SPDX-License-Identifier: Apache-2.0 + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace AsyncAdePT { + +/// @brief Release the tables owned by `G4HepEmData` and then delete the outer object. +/// @details +/// This is the cleanup path for the fully owned `fData` member. It is separate +/// from the class destructor because `std::unique_ptr` needs a deleter for the +/// deep cleanup before the outer `G4HepEmData` allocation itself can be deleted. +void AdePTG4HepEmState::DataDeleter::operator()(G4HepEmData *data) const +{ + if (data == nullptr) return; + FreeG4HepEmData(data); + delete data; +} + +/// @brief Release the copied HepEm parameters and then delete the outer object. +/// @details +/// The copied `G4HepEmParameters` block is fully owned by +/// `AdePTG4HepEmState`. This deep +/// cleanup therefore releases both the host-side per-region array and the +/// device-side mirror created during transport upload before deleting the outer +/// `G4HepEmParameters` allocation itself. +void AdePTG4HepEmState::ParametersDeleter::operator()(G4HepEmParameters *parameters) const +{ + if (parameters == nullptr) return; + FreeG4HepEmParameters(parameters); + delete parameters; +} + +/// @brief Rebuild a complete AdePT-owned set of host-side G4HepEm inputs from the supplied config. +/// @details +/// `AdePTG4HepEmState` owns two different G4HepEm objects: +/// - a deep copy of the `G4HepEmParameters` stored in the supplied `G4HepEmConfig` +/// - a freshly rebuilt `G4HepEmData` derived from that copied parameter block +/// +/// We must copy `G4HepEmParameters` because the original object remains owned by +/// the worker-local `G4HepEmConfig`, while the shared AdePT transport can outlive +/// the worker that first created it. `G4HepEmData` is rebuilt here directly, so +/// it is already fully owned by AdePT and does not need a second copy step. +AdePTG4HepEmState::AdePTG4HepEmState(G4HepEmConfig *hepEmConfig) + : fData(new G4HepEmData), fParameters(new G4HepEmParameters) +{ + if (hepEmConfig == nullptr) { + throw std::runtime_error("AdePTG4HepEmState requires a non-null G4HepEmConfig."); + } + + G4HepEmParameters *sourceParameters = hepEmConfig->GetG4HepEmParameters(); + if (sourceParameters == nullptr) { + throw std::runtime_error("AdePTG4HepEmState requires initialized G4HepEmParameters in the supplied config."); + } + + // Deep-copy the G4HepEmParameters so the shared transport does not keep a + // pointer into a worker-owned G4HepEmConfig. + *fParameters = *sourceParameters; + fParameters->fParametersPerRegion = nullptr; +#ifdef G4HepEm_CUDA_BUILD + fParameters->fParametersPerRegion_gpu = nullptr; +#endif + if (sourceParameters->fNumRegions > 0) { + if (sourceParameters->fParametersPerRegion == nullptr) { + throw std::runtime_error("AdePTG4HepEmState requires initialized per-region G4HepEmParameters."); + } + fParameters->fParametersPerRegion = new G4HepEmRegionParmeters[sourceParameters->fNumRegions]; + std::copy_n(sourceParameters->fParametersPerRegion, sourceParameters->fNumRegions, + fParameters->fParametersPerRegion); + } + + // Rebuild the G4HepEmData tables from the copied G4HepEmParameters so the + // transport owns a complete, self-contained set of host-side inputs. + InitG4HepEmData(fData.get()); + InitMaterialAndCoupleData(fData.get(), fParameters.get()); + + // Build all EM species + InitElectronData(fData.get(), fParameters.get(), true); + InitElectronData(fData.get(), fParameters.get(), false); + InitGammaData(fData.get(), fParameters.get()); + + G4HepEmMatCutData *cutData = fData->fTheMatCutData; + G4cout << "fNumG4MatCuts = " << cutData->fNumG4MatCuts << ", fNumMatCutData = " << cutData->fNumMatCutData << G4endl; +} + +AdePTG4HepEmState::~AdePTG4HepEmState() = default; + +} // namespace AsyncAdePT diff --git a/src/AdePTHepEmHostData.cpp b/src/AdePTHepEmHostData.cpp deleted file mode 100644 index c3cf8f00..00000000 --- a/src/AdePTHepEmHostData.cpp +++ /dev/null @@ -1,58 +0,0 @@ -// SPDX-FileCopyrightText: 2026 CERN -// SPDX-License-Identifier: Apache-2.0 - -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -namespace AsyncAdePT { - -/// @brief Release the tables owned by `G4HepEmData` and then delete the outer object. -/// @details -/// This is the cleanup path for the fully owned `fData` member. It is separate -/// from the class destructor because `std::unique_ptr` needs a deleter for the -/// deep cleanup before the outer `G4HepEmData` allocation itself can be deleted. -void HepEmHostData::DataDeleter::operator()(G4HepEmData *data) const -{ - if (data == nullptr) return; - FreeG4HepEmData(data); - delete data; -} - -/// @brief Rebuild a complete host-side G4HepEm view from the supplied G4HepEm config. -HepEmHostData::HepEmHostData(G4HepEmConfig *hepEmConfig) - : fData(new G4HepEmData), fParameters(hepEmConfig->GetG4HepEmParameters()) -{ - // Rebuild the HepEm tables from the G4HepEm-owned config so the host-side - // transport preparation has a complete, self-contained view of the data. - InitG4HepEmData(fData.get()); - InitMaterialAndCoupleData(fData.get(), fParameters); - - // Build all EM species - InitElectronData(fData.get(), fParameters, true); - InitElectronData(fData.get(), fParameters, false); - InitGammaData(fData.get(), fParameters); - - G4HepEmMatCutData *cutData = fData->fTheMatCutData; - G4cout << "fNumG4MatCuts = " << cutData->fNumG4MatCuts << ", fNumMatCutData = " << cutData->fNumMatCutData << G4endl; -} - -/// @brief Release the GPU mirror of the borrowed HepEm parameters. -/// @details -/// The host-side `G4HepEmParameters` remain owned by the `G4HepEmConfig`. -/// AdePT only owns the device allocation created by `CopyG4HepEmParametersToGPU`, -/// so this destructor releases just that GPU-side mirror. The owned `G4HepEmData` -/// cleanup is handled independently by `DataDeleter`. -HepEmHostData::~HepEmHostData() -{ - FreeG4HepEmParametersOnGPU(fParameters); -} - -} // namespace AsyncAdePT diff --git a/src/AdePTTrackingManager.cc b/src/AdePTTrackingManager.cc index d75e3d53..8b95828d 100644 --- a/src/AdePTTrackingManager.cc +++ b/src/AdePTTrackingManager.cc @@ -36,11 +36,10 @@ std::weak_ptr &SharedAdePTTransportStorage() return transport; } -std::shared_ptr GetSharedAdePTTransport(AdePTConfiguration &conf, - std::unique_ptr hepEmHostData, - adeptint::VolAuxData *auxData, - const adeptint::WDTHostPacked &wdtPacked, - const std::vector &uniformFieldValues) +std::shared_ptr GetSharedAdePTTransport( + AdePTConfiguration &conf, std::unique_ptr adeptG4HepEmState, + adeptint::VolAuxData *auxData, const adeptint::WDTHostPacked &wdtPacked, + const std::vector &uniformFieldValues) { auto &transport = SharedAdePTTransportStorage(); // weak_ptr::lock() promotes the stored weak reference to a shared_ptr if the @@ -50,7 +49,7 @@ std::shared_ptr GetSharedAdePTTransport(AdePTConfiguration &conf } auto created = - std::make_shared(conf, std::move(hepEmHostData), auxData, wdtPacked, uniformFieldValues); + std::make_shared(conf, std::move(adeptG4HepEmState), auxData, wdtPacked, uniformFieldValues); transport = created; return created; } @@ -95,22 +94,23 @@ void AdePTTrackingManager::InitializeSharedAdePTTransport() if (fAdePTConfiguration->GetCovfieBfieldFile() == "") std::cout << "No magnetic field file provided!" << std::endl; #endif const auto uniformFieldValues = fGeant4Integration.GetUniformField(); - auto hepEmHostData = std::make_unique(fHepEmTrackingManager->GetConfig()); + auto adeptG4HepEmState = std::make_unique(fHepEmTrackingManager->GetConfig()); // Check VecGeom geometry matches Geant4 before deriving any geometry metadata for transport. - fGeant4Integration.CheckGeometry(hepEmHostData->GetData()); + fGeant4Integration.CheckGeometry(adeptG4HepEmState->GetData()); // Initialize auxiliary per-LV data and collect the raw WDT metadata on the Geant4 side. auto *auxData = new adeptint::VolAuxData[vecgeom::GeoManager::Instance().GetRegisteredVolumesCount()]; adeptint::WDTHostRaw wdtRaw; - fGeant4Integration.InitVolAuxData(auxData, hepEmHostData->GetData(), fHepEmTrackingManager.get(), + fGeant4Integration.InitVolAuxData(auxData, adeptG4HepEmState->GetData(), fHepEmTrackingManager.get(), fAdePTConfiguration->GetTrackInAllRegions(), fAdePTConfiguration->GetGPURegionNames(), wdtRaw); adeptint::WDTHostPacked wdtPacked = adeptint::PackWDT(wdtRaw); - // Create the shared AdePT transport engine once the full host-side preparation is complete. - fAdeptTransport = - GetSharedAdePTTransport(*fAdePTConfiguration, std::move(hepEmHostData), auxData, wdtPacked, uniformFieldValues); + // Build the AdePT-owned G4HepEm state on the Geant4 side, then move that + // ownership into the shared transport once all host-side preparation is complete. + fAdeptTransport = GetSharedAdePTTransport(*fAdePTConfiguration, std::move(adeptG4HepEmState), auxData, wdtPacked, + uniformFieldValues); } void AdePTTrackingManager::InitializeAdePT() From 361307ecea868c65c6a07de44676796b6cabbba4 Mon Sep 17 00:00:00 2001 From: SeverinDiederichs Date: Sat, 21 Mar 2026 13:49:08 +0100 Subject: [PATCH 5/6] improve comments --- .../AdePT/integration/AdePTTrackingManager.hh | 11 +++++++++++ src/AdePTTrackingManager.cc | 19 ++++++++++++++----- 2 files changed, 25 insertions(+), 5 deletions(-) diff --git a/include/AdePT/integration/AdePTTrackingManager.hh b/include/AdePT/integration/AdePTTrackingManager.hh index f20232ba..93fd49f1 100644 --- a/include/AdePT/integration/AdePTTrackingManager.hh +++ b/include/AdePT/integration/AdePTTrackingManager.hh @@ -62,6 +62,17 @@ private: const G4NavigationHistory *aG4NavigationHistory = nullptr); /// @brief Perform the one-time shared AdePT transport initialization on the first Geant4 worker. + /// @details + /// The first worker prepares all host-side inputs needed by transport: + /// - the uniform magnetic-field values + /// - the AdePT-owned `AdePTG4HepEmState` + /// - geometry consistency checks + /// - `VolAuxData` + /// - packed WDT metadata + /// + /// Once that host-side preparation is complete, the worker creates the + /// shared `AsyncAdePTTransport`. The transport constructor then performs the + /// corresponding one-time device initialization and upload. void InitializeSharedAdePTTransport(); std::unique_ptr fHepEmTrackingManager; diff --git a/src/AdePTTrackingManager.cc b/src/AdePTTrackingManager.cc index 8b95828d..888d2226 100644 --- a/src/AdePTTrackingManager.cc +++ b/src/AdePTTrackingManager.cc @@ -48,6 +48,9 @@ std::shared_ptr GetSharedAdePTTransport( return existing; } + // Create the shared AdePT transport engine on the first worker thread. At + // this point all required host-side inputs have already been prepared, so the + // transport constructor can perform the one-time device initialization. auto created = std::make_shared(conf, std::move(adeptG4HepEmState), auxData, wdtPacked, uniformFieldValues); transport = created; @@ -93,6 +96,9 @@ void AdePTTrackingManager::InitializeSharedAdePTTransport() std::cout << "Reading in covfie file for magnetic field: " << fAdePTConfiguration->GetCovfieBfieldFile() << std::endl; if (fAdePTConfiguration->GetCovfieBfieldFile() == "") std::cout << "No magnetic field file provided!" << std::endl; #endif + // Prepare the complete host-side input package before constructing the + // shared transport. The transport constructor then performs the one-time + // device-side initialization from these prepared inputs. const auto uniformFieldValues = fGeant4Integration.GetUniformField(); auto adeptG4HepEmState = std::make_unique(fHepEmTrackingManager->GetConfig()); @@ -107,8 +113,9 @@ void AdePTTrackingManager::InitializeSharedAdePTTransport() fAdePTConfiguration->GetGPURegionNames(), wdtRaw); adeptint::WDTHostPacked wdtPacked = adeptint::PackWDT(wdtRaw); - // Build the AdePT-owned G4HepEm state on the Geant4 side, then move that - // ownership into the shared transport once all host-side preparation is complete. + // Move the fully prepared host-side package into the shared transport. The + // first worker creates the transport here; later workers only retrieve the + // already-created shared instance. fAdeptTransport = GetSharedAdePTTransport(*fAdePTConfiguration, std::move(adeptG4HepEmState), auxData, wdtPacked, uniformFieldValues); } @@ -130,7 +137,9 @@ void AdePTTrackingManager::InitializeAdePT() static std::condition_variable initCV; static bool commonInitDone = false; - // Global initialization: only done once by the first worker thread + // Global initialization: only done once by the first worker thread. This + // first worker closes the geometry, prepares the shared host-side inputs, + // and creates the shared transport, which performs the one-time device init. std::call_once(onceFlag, [&]() { // get number of threads from config, if available if (fNumThreads <= 0) { @@ -184,8 +193,8 @@ void AdePTTrackingManager::InitializeAdePT() // Now the fNumThreads is known and all workers can initialize fAdePTConfiguration->SetNumThreads(fNumThreads); - // The shared AdePT transport was already created and initialized by the first worker. - // The remaining workers only retrieve the shared pointer here. + // The shared AdePT transport was already created and device-initialized by + // the first worker. The remaining workers only retrieve the shared pointer. fAdeptTransport = GetSharedAdePTTransport(); // Initialize the GPU region list From 347252e0dd7659b3b3de178201b0cba053bc7c06 Mon Sep 17 00:00:00 2001 From: SeverinDiederichs Date: Thu, 19 Mar 2026 09:14:35 +0100 Subject: [PATCH 6/6] factor out geometry part of AdePTGeant4Integration --- CMakeLists.txt | 1 + .../integration/AdePTGeant4Integration.hh | 30 -- .../AdePT/integration/AdePTGeometryBridge.hh | 59 +++ src/AdePTGeant4Integration.cpp | 326 +--------------- src/AdePTGeometryBridge.cpp | 349 ++++++++++++++++++ src/AdePTTrackingManager.cc | 13 +- 6 files changed, 418 insertions(+), 360 deletions(-) create mode 100644 include/AdePT/integration/AdePTGeometryBridge.hh create mode 100644 src/AdePTGeometryBridge.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f238296a..6b46e296 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -445,6 +445,7 @@ endif() #----------------------------------------------------------------------------# set(ADEPT_G4_INTEGRATION_SRCS src/AdePTG4HepEmState.cpp + src/AdePTGeometryBridge.cpp src/G4HepEmTrackingManagerSpecialized.cc src/AdePTTrackingManager.cc src/G4EmStandardPhysics_AdePT.cc diff --git a/include/AdePT/integration/AdePTGeant4Integration.hh b/include/AdePT/integration/AdePTGeant4Integration.hh index 0d3466eb..844da0a2 100644 --- a/include/AdePT/integration/AdePTGeant4Integration.hh +++ b/include/AdePT/integration/AdePTGeant4Integration.hh @@ -17,11 +17,8 @@ #include #include -#include #include -struct G4HepEmData; - namespace AdePTGeant4Integration_detail { struct ScoringObjects; struct Deleter { @@ -40,31 +37,6 @@ public: AdePTGeant4Integration(AdePTGeant4Integration &&) = default; AdePTGeant4Integration &operator=(AdePTGeant4Integration &&) = default; -#ifdef VECGEOM_GDML_SUPPORT - /// @brief Initializes VecGeom geometry - /// @details Currently VecGeom geometry is initialized by loading it from a GDML file, - /// however ideally this function will call the G4 to VecGeom geometry converter - static void CreateVecGeomWorld(/*Temporary parameter*/ std::string filename); -#endif - - /// @brief Construct VecGeom geometry from Geant4 physical volume - /// @details This calls the G4VG converter - /// @throws std::runtime_error if input or output volumes are nullptr - static void CreateVecGeomWorld(G4VPhysicalVolume const *physvol); - - /// @brief This function compares G4 and VecGeom geometries and reports any differences - static void CheckGeometry(G4HepEmData const *hepEmData); - - /// @brief Fills the auxiliary data needed for AdePT - static void InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmData const *hepEmData, - G4HepEmTrackingManagerSpecialized *hepEmTM, bool trackInAllRegions, - std::vector const *gpuRegionNames, adeptint::WDTHostRaw &wdtRaw); - - /// @brief Returns a mapping of VecGeom placed volume IDs to Geant4 physical volumes and a mapping of VecGeom logical - /// volume IDs to Geant4 logical volumes - static void MapVecGeomToG4(std::vector &vecgeomPvToG4Map, - std::vector &vecgeomLvToG4Map); - /// @brief Reconstructs GPU hits on host and calls the user-defined sensitive detector code void ProcessGPUStep(std::span gpuSteps, bool const callUserSteppingAction = false, bool const callUserTrackingaction = false); @@ -127,8 +99,6 @@ private: // helper class to provide the CPU-only data for the returning GPU tracks std::unique_ptr fHostTrackDataMapper; - static std::vector fglobal_vecgeom_pv_to_g4_map; - static std::vector fglobal_vecgeom_lv_to_g4_map; std::unique_ptr fScoringObjects{nullptr}; }; diff --git a/include/AdePT/integration/AdePTGeometryBridge.hh b/include/AdePT/integration/AdePTGeometryBridge.hh new file mode 100644 index 00000000..759bf54f --- /dev/null +++ b/include/AdePT/integration/AdePTGeometryBridge.hh @@ -0,0 +1,59 @@ +// SPDX-FileCopyrightText: 2026 CERN +// SPDX-License-Identifier: Apache-2.0 + +#ifndef ADEPT_GEOMETRY_BRIDGE_HH +#define ADEPT_GEOMETRY_BRIDGE_HH + +#include +#include + +#include + +#include +#include + +class G4LogicalVolume; +class G4VPhysicalVolume; +struct G4HepEmData; + +/// @brief Global bridge between Geant4 host geometry and the VecGeom world used by AdePT. +/// @details +/// The geometry world and the VecGeom-to-Geant4 lookup tables are process-global, +/// so this bridge provides the corresponding global services needed during setup and +/// touchable reconstruction. +class AdePTGeometryBridge { +public: +#ifdef VECGEOM_GDML_SUPPORT + /// @brief Initializes the VecGeom world from a GDML file. + /// @details + /// This is the temporary path used when a standalone VecGeom GDML description is + /// provided instead of converting the Geant4 world through G4VG. + static void CreateVecGeomWorld(std::string filename); +#endif + + /// @brief Converts the Geant4 world volume into the process-global VecGeom world. + /// @throws std::runtime_error if the input Geant4 world or resulting VecGeom world is null. + static void CreateVecGeomWorld(G4VPhysicalVolume const *physvol); + + /// @brief Verifies that the Geant4 and VecGeom geometries match. + static void CheckGeometry(G4HepEmData const *hepEmData); + + /// @brief Fills the auxiliary per-volume data needed by AdePT. + static void InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmData const *hepEmData, + G4HepEmTrackingManagerSpecialized *hepEmTM, bool trackInAllRegions, + std::vector const *gpuRegionNames, adeptint::WDTHostRaw &wdtRaw); + + /// @brief Returns the Geant4 placed volume matching a VecGeom placed volume. + /// @throws std::runtime_error if the VecGeom placed volume is not present in the global lookup table. + static G4VPhysicalVolume const *GetG4PhysicalVolume(vecgeom::VPlacedVolume const *placedVolume); + +private: + /// @brief Builds the lookup tables from VecGeom placed/logical volume ids to Geant4 volumes. + static void MapVecGeomToG4(std::vector &vecgeomPvToG4Map, + std::vector &vecgeomLvToG4Map); + + static std::vector fGlobalVecGeomPvToG4Map; + static std::vector fGlobalVecGeomLvToG4Map; +}; + +#endif diff --git a/src/AdePTGeant4Integration.cpp b/src/AdePTGeant4Integration.cpp index 0d15e8e3..4c3bf3be 100644 --- a/src/AdePTGeant4Integration.cpp +++ b/src/AdePTGeant4Integration.cpp @@ -2,37 +2,25 @@ // SPDX-License-Identifier: Apache-2.0 #include +#include -#include -#ifdef VECGEOM_GDML_SUPPORT -#include -#endif #include #include #include -#include #include -#include #include #include #include -#include #include -#include -#include #include -#include -#include #include #include "G4Electron.hh" #include "G4Positron.hh" #include "G4Gamma.hh" -#include - #include #include @@ -169,318 +157,8 @@ void Deleter::operator()(ScoringObjects *ptr) } // namespace AdePTGeant4Integration_detail -std::vector AdePTGeant4Integration::fglobal_vecgeom_pv_to_g4_map; -std::vector AdePTGeant4Integration::fglobal_vecgeom_lv_to_g4_map; - AdePTGeant4Integration::~AdePTGeant4Integration() {} -void AdePTGeant4Integration::MapVecGeomToG4(std::vector &vecgeomPvToG4Map, - std::vector &vecgeomLvToG4Map) -{ - const G4VPhysicalVolume *g4world = - G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume(); - const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); - - // recursive geometry visitor lambda matching one by one Geant4 and VecGeom logical volumes - typedef std::function func_t; - func_t visitGeometry = [&](G4VPhysicalVolume const *g4_pvol, vecgeom::VPlacedVolume const *vg_pvol) { - const auto g4_lvol = g4_pvol->GetLogicalVolume(); - const auto vg_lvol = vg_pvol->GetLogicalVolume(); - - // Initialize mapping of Vecgeom PlacedVolume IDs to G4 PhysicalVolume* - vecgeomPvToG4Map.resize(std::max(vecgeomPvToG4Map.size(), vg_pvol->id() + 1), nullptr); - vecgeomPvToG4Map[vg_pvol->id()] = g4_pvol; - // Initialize mapping of Vecgeom LogicalVolume IDs to G4 LogicalVolume* - vecgeomLvToG4Map.resize(std::max(vecgeomLvToG4Map.size(), vg_lvol->id() + 1), nullptr); - vecgeomLvToG4Map[vg_lvol->id()] = g4_lvol; - - // Now do the daughters - for (size_t id = 0; id < g4_lvol->GetNoDaughters(); ++id) { - auto g4pvol_d = g4_lvol->GetDaughter(id); - auto pvol_d = vg_lvol->GetDaughters()[id]; - - visitGeometry(g4pvol_d, pvol_d); - } - }; - visitGeometry(g4world, vecgeomWorld); -} - -#ifdef VECGEOM_GDML_SUPPORT -void AdePTGeant4Integration::CreateVecGeomWorld(std::string filename) -{ - // Import the gdml file into VecGeom - vecgeom::GeoManager::Instance().SetTransformationCacheDepth(0); - vgdml::Parser vgdmlParser; - auto middleWare = vgdmlParser.Load(filename, false, mm); - if (middleWare == nullptr) { - std::cerr << "Failed to read geometry from GDML file '" << filename << "'" << G4endl; - return; - } - - // Generate the mapping of VecGeom volume IDs to Geant4 physical volumes - MapVecGeomToG4(fglobal_vecgeom_pv_to_g4_map, fglobal_vecgeom_lv_to_g4_map); - - const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); - if (vecgeomWorld == nullptr) { - std::cerr << "GeoManager vecgeomWorld volume is nullptr" << G4endl; - return; - } -} -#endif - -void AdePTGeant4Integration::CreateVecGeomWorld(G4VPhysicalVolume const *physvol) -{ - // EXPECT: a non-null input volume - if (physvol == nullptr) { - throw std::runtime_error("AdePTGeant4Integration::CreateVecGeomWorld : Input Geant4 Physical Volume is nullptr"); - } - - vecgeom::GeoManager::Instance().SetTransformationCacheDepth(0); - g4vg::Options options; - options.reflection_factory = false; - auto conversion = g4vg::convert(physvol, options); - vecgeom::GeoManager::Instance().SetWorldAndClose(conversion.world); - - // Get the mapping of VecGeom volume IDs to Geant4 physical volumes from g4vg - fglobal_vecgeom_pv_to_g4_map = conversion.physical_volumes; - fglobal_vecgeom_lv_to_g4_map = conversion.logical_volumes; - - // EXPECT: we finish with a non-null VecGeom host geometry - vecgeom::VPlacedVolume const *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); - if (vecgeomWorld == nullptr) { - throw std::runtime_error("AdePTGeant4Integration::CreateVecGeomWorld : Output VecGeom Physical Volume is nullptr"); - } -} - -namespace { -struct VisitContext { - const int *g4tohepmcindex; - std::size_t nvolumes; - G4HepEmData const *hepEmData; -}; - -/// Recursive geometry visitor matching one by one Geant4 and VecGeom logical volumes -void visitGeometry(G4VPhysicalVolume const *g4_pvol, vecgeom::VPlacedVolume const *vg_pvol, const VisitContext &context) -{ - const auto g4_lvol = g4_pvol->GetLogicalVolume(); - const auto vg_lvol = vg_pvol->GetLogicalVolume(); - - // Geant4 Parameterised/Replica volumes are represented with direct placements in VecGeom - // To accurately compare the number of daughters, we need to sum multiplicity on Geant4 side - const size_t nd = g4_lvol->GetNoDaughters(); - size_t nd_converted = 0; - for (size_t daughter_id = 0; daughter_id < nd; ++daughter_id) { - const G4VPhysicalVolume *daughter_pvol = g4_lvol->GetDaughter(daughter_id); - nd_converted += daughter_pvol->GetMultiplicity(); - } - - const auto daughters = vg_lvol->GetDaughters(); - - if (nd_converted != daughters.size()) - throw std::runtime_error("Fatal: CheckGeometry: Mismatch in number of daughters"); - - // Check if transformations are matching - // As above, with Parameterized/Replica volumes, we need to compare the transforms between - // the VG direct placement and that for the Parameterised/Replicated volume given the same copy - // number as that of the VG physical volume. - // NOTE: - // 1. Nasty const_cast as currently known way to get transform is to directly transform the - // Geant4 phys vol before extracting, which is non-const.... - // 2. ...this does modify the physical volume, but this is _probably_ o.k. as actual navigation - // will reset things o.k. - if (G4VPVParameterisation *param = g4_pvol->GetParameterisation()) { - param->ComputeTransformation(vg_pvol->GetCopyNo(), const_cast(g4_pvol)); - } else if (auto *replica = dynamic_cast(const_cast(g4_pvol))) { - G4ReplicaNavigation nav; - nav.ComputeTransformation(vg_pvol->GetCopyNo(), replica); - } - - const auto g4trans = g4_pvol->GetTranslation(); - const G4RotationMatrix *g4rot = g4_pvol->GetRotation(); - G4RotationMatrix idrot; - const auto vgtransformation = vg_pvol->GetTransformation(); - constexpr double epsil = 1.e-8; - - for (int i = 0; i < 3; ++i) { - if (std::abs(g4trans[i] - vgtransformation->Translation(i)) > epsil) - throw std::runtime_error( - std::string("Fatal: CheckGeometry: Mismatch between Geant4 translation for physical volume") + - vg_pvol->GetName()); - } - - // check if VecGeom and Geant4 (local) transformations are matching. Not optimized, this will re-check - // already checked placed volumes when re-visiting the same volumes in different branches - if (!g4rot) g4rot = &idrot; - for (int row = 0; row < 3; ++row) { - for (int col = 0; col < 3; ++col) { - int i = row + 3 * col; - if (std::abs((*g4rot)(row, col) - vgtransformation->Rotation(i)) > epsil) - throw std::runtime_error( - std::string("Fatal: CheckGeometry: Mismatch between Geant4 rotation for physical volume") + - vg_pvol->GetName()); - } - } - - // Check the couples - if (g4_lvol->GetMaterialCutsCouple() == nullptr) - throw std::runtime_error("Fatal: CheckGeometry: G4LogicalVolume " + std::string(g4_lvol->GetName()) + - std::string(" has no material-cuts couple")); - const int g4mcindex = g4_lvol->GetMaterialCutsCouple()->GetIndex(); - const int hepemmcindex = context.g4tohepmcindex[g4mcindex]; - // Check consistency with G4HepEm data - if (context.hepEmData->fTheMatCutData->fMatCutData[hepemmcindex].fG4MatCutIndex != g4mcindex) - throw std::runtime_error("Fatal: CheckGeometry: Mismatch between Geant4 mcindex and corresponding G4HepEm index"); - if (vg_lvol->id() >= context.nvolumes) - throw std::runtime_error("Fatal: CheckGeometry: Volume id larger than number of volumes"); - - // Now do the daughters - for (size_t id = 0; id < g4_lvol->GetNoDaughters(); ++id) { - const auto g4pvol_d = g4_lvol->GetDaughter(id); - const auto pvol_d = vg_lvol->GetDaughters()[id]; - - visitGeometry(g4pvol_d, pvol_d, context); - } -} -} // namespace - -void AdePTGeant4Integration::CheckGeometry(G4HepEmData const *hepEmData) -{ - const G4VPhysicalVolume *g4world = - G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume(); - const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); - const int *g4tohepmcindex = hepEmData->fTheMatCutData->fG4MCIndexToHepEmMCIndex; - const auto nvolumes = vecgeom::GeoManager::Instance().GetRegisteredVolumesCount(); - - std::cout << "Visiting geometry ...\n"; - const VisitContext context{g4tohepmcindex, nvolumes, hepEmData}; - visitGeometry(g4world, vecgeomWorld, context); - std::cout << "Visiting geometry done\n"; -} - -void AdePTGeant4Integration::InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmData const *hepEmData, - G4HepEmTrackingManagerSpecialized *hepEmTM, bool trackInAllRegions, - std::vector const *gpuRegionNames, - adeptint::WDTHostRaw &wdtRaw) -{ - - // Note: the hepEmTM must be passed as an argument despite the member fHepEmTrackingManager, - // as InitVolAuxData is a static function and cannot call member variables - wdtRaw.ekinMin = (float)hepEmTM->GetWDTKineticEnergyLimit(); - - const G4VPhysicalVolume *g4world = - G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume(); - const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); - const int *g4tohepmcindex = hepEmData->fTheMatCutData->fG4MCIndexToHepEmMCIndex; - - // We need to go from region names to G4Region - std::vector gpuRegions{}; - if (!trackInAllRegions) { - for (std::string regionName : *(gpuRegionNames)) { - G4Region *region = G4RegionStore::GetInstance()->GetRegion(regionName); - gpuRegions.push_back(region); - } - } - - // recursive geometry visitor lambda matching one by one Geant4 and VecGeom logical volumes - typedef std::function - func_t; - func_t visitGeometry = [&](G4VPhysicalVolume const *g4_pvol, vecgeom::VPlacedVolume const *vg_pvol, - vecgeom::NavigationState currentNavState) { - const auto g4_lvol = g4_pvol->GetLogicalVolume(); - const auto vg_lvol = vg_pvol->GetLogicalVolume(); - - // Push this placed volume into the running NavState - currentNavState.Push(vg_pvol); - currentNavState.SetBoundaryState(false); - - // Fill the MCC index in the array - int g4mcindex = g4_lvol->GetMaterialCutsCouple()->GetIndex(); - int hepemmcindex = g4tohepmcindex[g4mcindex]; - volAuxData[vg_lvol->id()].fMCIndex = hepemmcindex; - - const int regionId = g4_lvol->GetRegion()->GetInstanceID(); - - // Check if the region is a Woodcock tracking region in G4HepEm - if (hepEmTM->IsWDTRegion(regionId)) { - // check if this logical volume is one of the declared WDT root LVs for this region - const int rootIMC = hepEmTM->GetWDTCoupleHepEmIndex(regionId, g4_lvol->GetInstanceID()); - if (rootIMC >= 0) { - // this placed volume belongs to a WDT root LV -> record a WDTRoot - int idx = (int)wdtRaw.roots.size(); - wdtRaw.roots.push_back(adeptint::WDTRoot{currentNavState, rootIMC}); - wdtRaw.regionToRootIndices[regionId].push_back(idx); - } - } - - // Check if the volume belongs to a GPU region - if (!trackInAllRegions) { - for (G4Region *gpuRegion : gpuRegions) { - if (g4_lvol->GetRegion() == gpuRegion) { - volAuxData[vg_lvol->id()].fGPUregionId = g4_lvol->GetRegion()->GetInstanceID(); - } - } - } else { - volAuxData[vg_lvol->id()].fGPUregionId = g4_lvol->GetRegion()->GetInstanceID(); - } - - if (g4_lvol->GetSensitiveDetector() != nullptr) { - if (volAuxData[vg_lvol->id()].fSensIndex < 0) { - G4cout << "VecGeom: Making " << vg_lvol->GetName() << " sensitive" << G4endl; - } - volAuxData[vg_lvol->id()].fSensIndex = 1; - } - // Now do the daughters - for (size_t id = 0; id < g4_lvol->GetNoDaughters(); ++id) { - auto g4pvol_d = g4_lvol->GetDaughter(id); - auto pvol_d = vg_lvol->GetDaughters()[id]; - - visitGeometry(g4pvol_d, pvol_d, currentNavState); - } - - // Pop NavState before before returning - currentNavState.Pop(); - }; - - // Initialize root NavState - vecgeom::NavigationState rootNavState; - - visitGeometry(g4world, vecgeomWorld, rootNavState); - - auto findRegionName = [](int rid) -> std::string { - for (auto *r : *G4RegionStore::GetInstance()) { - if (r && r->GetInstanceID() == rid) return r->GetName(); - } - return std::string(""); - }; - - std::cout << "\n=== Woodcock tracking summary (host) ===\n"; - std::cout << "KineticEnergyLimit = " << wdtRaw.ekinMin << " [G4 units]\n"; - std::cout << "Total WDT roots found: " << wdtRaw.roots.size() << std::endl; - std::cout << "Regions with WDT: " << wdtRaw.regionToRootIndices.size() << std::endl; - - if (wdtRaw.regionToRootIndices.empty()) { - std::cout << " (none)\n"; - } else { - for (const auto &kv : wdtRaw.regionToRootIndices) { - const int rid = kv.first; - const auto &idxs = kv.second; - std::cout << "\nRegionID " << rid << " (" << findRegionName(rid) << "): " << idxs.size() - << " root placed-volume(s)\n"; - - for (size_t i = 0; i < idxs.size(); ++i) { - const int rootIdx = idxs[i]; - const auto &root = wdtRaw.roots[rootIdx]; - - std::cout << " [" << i << "] hepemIMC=" << root.hepemIMC << "\n"; - std::cout << " NavState (level=" << root.root.GetLevel() << "):\n"; - // vecgeom::NavigationState::Print() prints the full stack - root.root.Print(); - } - } - } - std::cout << "=== End Woodcock tracking summary ===\n\n"; -} - G4Track *AdePTGeant4Integration::ConstructSecondaryTrackInPlace(GPUHit const *secHit) const { auto &so = *fScoringObjects; @@ -696,7 +374,7 @@ void AdePTGeant4Integration::FillG4NavigationHistory(const vecgeom::NavigationSt // While we are in levels shallower than the history depth, it may be that we already // have the correct volume in the history assert(aNavState.At(aLevel)); - pnewvol = const_cast(fglobal_vecgeom_pv_to_g4_map[aNavState.At(aLevel)->id()]); + pnewvol = const_cast(AdePTGeometryBridge::GetG4PhysicalVolume(aNavState.At(aLevel))); assert(pnewvol != nullptr); if (!pnewvol) throw std::runtime_error("VecGeom volume not found in G4 mapping!"); diff --git a/src/AdePTGeometryBridge.cpp b/src/AdePTGeometryBridge.cpp new file mode 100644 index 00000000..1f76cc90 --- /dev/null +++ b/src/AdePTGeometryBridge.cpp @@ -0,0 +1,349 @@ +// SPDX-FileCopyrightText: 2026 CERN +// SPDX-License-Identifier: Apache-2.0 + +#include + +#include +#ifdef VECGEOM_GDML_SUPPORT +#include +#endif +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +std::vector AdePTGeometryBridge::fGlobalVecGeomPvToG4Map; +std::vector AdePTGeometryBridge::fGlobalVecGeomLvToG4Map; + +void AdePTGeometryBridge::MapVecGeomToG4(std::vector &vecgeomPvToG4Map, + std::vector &vecgeomLvToG4Map) +{ + const G4VPhysicalVolume *g4world = + G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume(); + const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); + + // Recursive geometry visitor lambda matching one by one Geant4 and VecGeom logical volumes. + using VisitFn = std::function; + VisitFn visitGeometry = [&](G4VPhysicalVolume const *g4_pvol, vecgeom::VPlacedVolume const *vg_pvol) { + const auto g4_lvol = g4_pvol->GetLogicalVolume(); + const auto vg_lvol = vg_pvol->GetLogicalVolume(); + + // Initialize mapping of VecGeom placed-volume ids to Geant4 physical volumes. + vecgeomPvToG4Map.resize(std::max(vecgeomPvToG4Map.size(), vg_pvol->id() + 1), nullptr); + vecgeomPvToG4Map[vg_pvol->id()] = g4_pvol; + + // Initialize mapping of VecGeom logical-volume ids to Geant4 logical volumes. + vecgeomLvToG4Map.resize(std::max(vecgeomLvToG4Map.size(), vg_lvol->id() + 1), nullptr); + vecgeomLvToG4Map[vg_lvol->id()] = g4_lvol; + + // Now do the daughters. + for (size_t id = 0; id < g4_lvol->GetNoDaughters(); ++id) { + visitGeometry(g4_lvol->GetDaughter(id), vg_lvol->GetDaughters()[id]); + } + }; + + visitGeometry(g4world, vecgeomWorld); +} + +#ifdef VECGEOM_GDML_SUPPORT +/// @brief Initialize the process-global VecGeom world from a GDML file. +void AdePTGeometryBridge::CreateVecGeomWorld(std::string filename) +{ + // Import the GDML file into VecGeom. + vecgeom::GeoManager::Instance().SetTransformationCacheDepth(0); + vgdml::Parser vgdmlParser; + auto middleWare = vgdmlParser.Load(filename, false, mm); + if (middleWare == nullptr) { + std::cerr << "Failed to read geometry from GDML file '" << filename << "'" << G4endl; + return; + } + + // Generate the mapping of VecGeom volume ids to Geant4 physical volumes. + MapVecGeomToG4(fGlobalVecGeomPvToG4Map, fGlobalVecGeomLvToG4Map); + + const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); + if (vecgeomWorld == nullptr) { + std::cerr << "GeoManager vecgeomWorld volume is nullptr" << G4endl; + return; + } +} +#endif + +/// @brief Convert the Geant4 world into the process-global VecGeom world via G4VG. +void AdePTGeometryBridge::CreateVecGeomWorld(G4VPhysicalVolume const *physvol) +{ + // EXPECT: a non-null input volume. + if (physvol == nullptr) { + throw std::runtime_error("AdePTGeometryBridge::CreateVecGeomWorld: Input Geant4 physical volume is nullptr"); + } + + vecgeom::GeoManager::Instance().SetTransformationCacheDepth(0); + g4vg::Options options; + options.reflection_factory = false; + auto conversion = g4vg::convert(physvol, options); + vecgeom::GeoManager::Instance().SetWorldAndClose(conversion.world); + + // Get the mapping of VecGeom volume ids to Geant4 physical volumes from G4VG. + fGlobalVecGeomPvToG4Map = conversion.physical_volumes; + fGlobalVecGeomLvToG4Map = conversion.logical_volumes; + + // EXPECT: we finish with a non-null VecGeom host geometry. + if (vecgeom::GeoManager::Instance().GetWorld() == nullptr) { + throw std::runtime_error("AdePTGeometryBridge::CreateVecGeomWorld: Output VecGeom world is nullptr"); + } +} + +namespace { +struct VisitContext { + const int *g4tohepmcindex; + std::size_t nvolumes; + G4HepEmData const *hepEmData; +}; + +/// @brief Recursively compare the Geant4 and VecGeom geometry trees. +void VisitGeometryForChecks(G4VPhysicalVolume const *g4_pvol, vecgeom::VPlacedVolume const *vg_pvol, + const VisitContext &context) +{ + const auto g4_lvol = g4_pvol->GetLogicalVolume(); + const auto vg_lvol = vg_pvol->GetLogicalVolume(); + + // Geant4 parameterised/replica volumes are represented with direct placements in VecGeom. + // To accurately compare the number of daughters, sum multiplicity on the Geant4 side. + const size_t nd = g4_lvol->GetNoDaughters(); + size_t nd_converted = 0; + for (size_t daughter_id = 0; daughter_id < nd; ++daughter_id) { + nd_converted += g4_lvol->GetDaughter(daughter_id)->GetMultiplicity(); + } + + const auto daughters = vg_lvol->GetDaughters(); + if (nd_converted != daughters.size()) { + throw std::runtime_error("Fatal: CheckGeometry: Mismatch in number of daughters"); + } + + // Check whether transformations are matching. + // As above, with parameterized/replica volumes we need to compare the transforms between + // the VecGeom direct placement and that for the parameterised/replicated volume given the + // same copy number as that of the VecGeom physical volume. + // NOTE: + // 1. This needs a const_cast because the current Geant4 API computes the transform by + // mutating the physical volume. + // 2. This does modify the physical volume, but navigation will reset things afterwards. + if (G4VPVParameterisation *param = g4_pvol->GetParameterisation()) { + param->ComputeTransformation(vg_pvol->GetCopyNo(), const_cast(g4_pvol)); + } else if (auto *replica = dynamic_cast(const_cast(g4_pvol))) { + G4ReplicaNavigation nav; + nav.ComputeTransformation(vg_pvol->GetCopyNo(), replica); + } + + const auto g4trans = g4_pvol->GetTranslation(); + const G4RotationMatrix *g4rot = g4_pvol->GetRotation(); + G4RotationMatrix idrot; + const auto vgtransformation = vg_pvol->GetTransformation(); + constexpr double epsil = 1.e-8; + + for (int i = 0; i < 3; ++i) { + if (std::abs(g4trans[i] - vgtransformation->Translation(i)) > epsil) { + throw std::runtime_error( + std::string("Fatal: CheckGeometry: Mismatch between Geant4 translation for physical volume") + + vg_pvol->GetName()); + } + } + + // Check if VecGeom and Geant4 local transformations are matching. + // This is not optimized and will re-check already-checked placed volumes when + // revisiting the same volumes in different branches. + if (!g4rot) g4rot = &idrot; + for (int row = 0; row < 3; ++row) { + for (int col = 0; col < 3; ++col) { + int i = row + 3 * col; + if (std::abs((*g4rot)(row, col) - vgtransformation->Rotation(i)) > epsil) { + throw std::runtime_error( + std::string("Fatal: CheckGeometry: Mismatch between Geant4 rotation for physical volume") + + vg_pvol->GetName()); + } + } + } + + if (g4_lvol->GetMaterialCutsCouple() == nullptr) { + throw std::runtime_error("Fatal: CheckGeometry: G4LogicalVolume " + std::string(g4_lvol->GetName()) + + std::string(" has no material-cuts couple")); + } + const int g4mcindex = g4_lvol->GetMaterialCutsCouple()->GetIndex(); + const int hepemmcindex = context.g4tohepmcindex[g4mcindex]; + // Check consistency with G4HepEm data. + if (context.hepEmData->fTheMatCutData->fMatCutData[hepemmcindex].fG4MatCutIndex != g4mcindex) { + throw std::runtime_error("Fatal: CheckGeometry: Mismatch between Geant4 mcindex and corresponding G4HepEm index"); + } + if (vg_lvol->id() >= context.nvolumes) { + throw std::runtime_error("Fatal: CheckGeometry: Volume id larger than number of volumes"); + } + + // Now do the daughters. + for (size_t id = 0; id < g4_lvol->GetNoDaughters(); ++id) { + VisitGeometryForChecks(g4_lvol->GetDaughter(id), vg_lvol->GetDaughters()[id], context); + } +} +} // namespace + +/// @brief Compare the Geant4 and VecGeom host geometries for consistency. +void AdePTGeometryBridge::CheckGeometry(G4HepEmData const *hepEmData) +{ + const G4VPhysicalVolume *g4world = + G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume(); + const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); + const int *g4tohepmcindex = hepEmData->fTheMatCutData->fG4MCIndexToHepEmMCIndex; + const auto nvolumes = vecgeom::GeoManager::Instance().GetRegisteredVolumesCount(); + + std::cout << "Visiting geometry ...\n"; + const VisitContext context{g4tohepmcindex, nvolumes, hepEmData}; + VisitGeometryForChecks(g4world, vecgeomWorld, context); + std::cout << "Visiting geometry done\n"; +} + +/// @brief Fill the auxiliary per-volume transport metadata used by AdePT. +void AdePTGeometryBridge::InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmData const *hepEmData, + G4HepEmTrackingManagerSpecialized *hepEmTM, bool trackInAllRegions, + std::vector const *gpuRegionNames, adeptint::WDTHostRaw &wdtRaw) +{ + // Note: the hepEmTM must be passed explicitly here, since this is now a stateless + // global bridge helper and therefore cannot reach any per-thread integration members. + wdtRaw.ekinMin = (float)hepEmTM->GetWDTKineticEnergyLimit(); + + const G4VPhysicalVolume *g4world = + G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume(); + const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld(); + const int *g4tohepmcindex = hepEmData->fTheMatCutData->fG4MCIndexToHepEmMCIndex; + + // We need to go from region names to G4Region objects. + std::vector gpuRegions{}; + if (!trackInAllRegions) { + for (const std::string ®ionName : *gpuRegionNames) { + gpuRegions.push_back(G4RegionStore::GetInstance()->GetRegion(regionName)); + } + } + + // Recursive geometry visitor lambda matching one by one Geant4 and VecGeom logical volumes. + using VisitFn = + std::function; + VisitFn visitGeometry = [&](G4VPhysicalVolume const *g4_pvol, vecgeom::VPlacedVolume const *vg_pvol, + vecgeom::NavigationState currentNavState) { + const auto g4_lvol = g4_pvol->GetLogicalVolume(); + const auto vg_lvol = vg_pvol->GetLogicalVolume(); + + // Push this placed volume into the running navigation state. + currentNavState.Push(vg_pvol); + currentNavState.SetBoundaryState(false); + + // Fill the material-cuts-couple index in the auxiliary array. + int g4mcindex = g4_lvol->GetMaterialCutsCouple()->GetIndex(); + int hepemmcindex = g4tohepmcindex[g4mcindex]; + volAuxData[vg_lvol->id()].fMCIndex = hepemmcindex; + + const int regionId = g4_lvol->GetRegion()->GetInstanceID(); + // Check if the region is a Woodcock tracking region in G4HepEm. + if (hepEmTM->IsWDTRegion(regionId)) { + // Check if this logical volume is one of the declared WDT root LVs for this region. + const int rootIMC = hepEmTM->GetWDTCoupleHepEmIndex(regionId, g4_lvol->GetInstanceID()); + if (rootIMC >= 0) { + // This placed volume belongs to a WDT root LV, so record a WDTRoot. + int idx = (int)wdtRaw.roots.size(); + wdtRaw.roots.push_back(adeptint::WDTRoot{currentNavState, rootIMC}); + wdtRaw.regionToRootIndices[regionId].push_back(idx); + } + } + + // Check if the volume belongs to a GPU region. + if (!trackInAllRegions) { + for (G4Region *gpuRegion : gpuRegions) { + if (g4_lvol->GetRegion() == gpuRegion) { + volAuxData[vg_lvol->id()].fGPUregionId = g4_lvol->GetRegion()->GetInstanceID(); + } + } + } else { + volAuxData[vg_lvol->id()].fGPUregionId = g4_lvol->GetRegion()->GetInstanceID(); + } + + if (g4_lvol->GetSensitiveDetector() != nullptr) { + if (volAuxData[vg_lvol->id()].fSensIndex < 0) { + G4cout << "VecGeom: Making " << vg_lvol->GetName() << " sensitive" << G4endl; + } + volAuxData[vg_lvol->id()].fSensIndex = 1; + } + + // Now do the daughters. + for (size_t id = 0; id < g4_lvol->GetNoDaughters(); ++id) { + visitGeometry(g4_lvol->GetDaughter(id), vg_lvol->GetDaughters()[id], currentNavState); + } + + // Pop the navigation state before returning. + currentNavState.Pop(); + }; + + // Initialize the root navigation state. + vecgeom::NavigationState rootNavState; + visitGeometry(g4world, vecgeomWorld, rootNavState); + + auto findRegionName = [](int rid) -> std::string { + for (auto *r : *G4RegionStore::GetInstance()) { + if (r && r->GetInstanceID() == rid) return r->GetName(); + } + return std::string(""); + }; + + std::cout << "\n=== Woodcock tracking summary (host) ===\n"; + std::cout << "KineticEnergyLimit = " << wdtRaw.ekinMin << " [G4 units]\n"; + std::cout << "Total WDT roots found: " << wdtRaw.roots.size() << std::endl; + std::cout << "Regions with WDT: " << wdtRaw.regionToRootIndices.size() << std::endl; + + if (wdtRaw.regionToRootIndices.empty()) { + std::cout << " (none)\n"; + } else { + for (const auto &[regionId, indices] : wdtRaw.regionToRootIndices) { + std::cout << "\nRegionID " << regionId << " (" << findRegionName(regionId) << "): " << indices.size() + << " root placed-volume(s)\n"; + + for (std::size_t i = 0; i < indices.size(); ++i) { + const auto &root = wdtRaw.roots[indices[i]]; + std::cout << " [" << i << "] hepemIMC=" << root.hepemIMC << "\n"; + std::cout << " NavState (level=" << root.root.GetLevel() << "):\n"; + // vecgeom::NavigationState::Print() prints the full stack. + root.root.Print(); + } + } + } + std::cout << "=== End Woodcock tracking summary ===\n\n"; +} + +/// @brief Resolve the Geant4 placed volume associated with a VecGeom placed volume. +G4VPhysicalVolume const *AdePTGeometryBridge::GetG4PhysicalVolume(vecgeom::VPlacedVolume const *placedVolume) +{ + if (placedVolume == nullptr) { + throw std::runtime_error("AdePTGeometryBridge::GetG4PhysicalVolume: Input VecGeom placed volume is nullptr"); + } + if (placedVolume->id() >= fGlobalVecGeomPvToG4Map.size()) { + throw std::runtime_error( + "AdePTGeometryBridge::GetG4PhysicalVolume: VecGeom placed volume id is outside the lookup table"); + } + + auto *g4Volume = fGlobalVecGeomPvToG4Map[placedVolume->id()]; + if (g4Volume == nullptr) { + throw std::runtime_error( + "AdePTGeometryBridge::GetG4PhysicalVolume: VecGeom placed volume not found in Geant4 mapping"); + } + return g4Volume; +} diff --git a/src/AdePTTrackingManager.cc b/src/AdePTTrackingManager.cc index 888d2226..4ee07016 100644 --- a/src/AdePTTrackingManager.cc +++ b/src/AdePTTrackingManager.cc @@ -2,6 +2,7 @@ // SPDX-License-Identifier: Apache-2.0 #include +#include #include "G4Threading.hh" #include "G4Track.hh" @@ -103,14 +104,14 @@ void AdePTTrackingManager::InitializeSharedAdePTTransport() auto adeptG4HepEmState = std::make_unique(fHepEmTrackingManager->GetConfig()); // Check VecGeom geometry matches Geant4 before deriving any geometry metadata for transport. - fGeant4Integration.CheckGeometry(adeptG4HepEmState->GetData()); + AdePTGeometryBridge::CheckGeometry(adeptG4HepEmState->GetData()); // Initialize auxiliary per-LV data and collect the raw WDT metadata on the Geant4 side. auto *auxData = new adeptint::VolAuxData[vecgeom::GeoManager::Instance().GetRegisteredVolumesCount()]; adeptint::WDTHostRaw wdtRaw; - fGeant4Integration.InitVolAuxData(auxData, adeptG4HepEmState->GetData(), fHepEmTrackingManager.get(), - fAdePTConfiguration->GetTrackInAllRegions(), - fAdePTConfiguration->GetGPURegionNames(), wdtRaw); + AdePTGeometryBridge::InitVolAuxData(auxData, adeptG4HepEmState->GetData(), fHepEmTrackingManager.get(), + fAdePTConfiguration->GetTrackInAllRegions(), + fAdePTConfiguration->GetGPURegionNames(), wdtRaw); adeptint::WDTHostPacked wdtPacked = adeptint::PackWDT(wdtRaw); // Move the fully prepared host-side package into the shared transport. The @@ -165,11 +166,11 @@ void AdePTTrackingManager::InitializeAdePT() auto *tman = G4TransportationManager::GetTransportationManager(); auto *world = tman->GetNavigatorForTracking()->GetWorldVolume(); std::cout << "Loading geometry via G4VG\n"; - AdePTGeant4Integration::CreateVecGeomWorld(world); + AdePTGeometryBridge::CreateVecGeomWorld(world); } else { #ifdef VECGEOM_GDML_SUPPORT std::cout << "Loading geometry via VGDML\n"; - AdePTGeant4Integration::CreateVecGeomWorld(fAdePTConfiguration->GetVecGeomGDML()); + AdePTGeometryBridge::CreateVecGeomWorld(fAdePTConfiguration->GetVecGeomGDML()); #endif }