Skip to content

Commit

Permalink
Using g4vg-style placed volume map when loading geometry from GDML (#343
Browse files Browse the repository at this point in the history 
)

For consistency, we now use a vector for mapping vecgeom PlacedVolume
IDs to Geant4 Physical volumes, as provided by g4vg.

For the time being we also generate this mapping on AdePT side rather
than using the G4VG-provided map as there are a small number of volumes
with different mappings when using either method.
  • Loading branch information
@JuanGonzalezCaminero
JuanGonzalezCaminero authored Feb 6, 2025
1 parent b1bda50 commit a0ac9bb
Show file tree
Hide file tree
Showing 3 changed files with 50 additions and 41 deletions.
2 changes: 0 additions & 2 deletions include/AdePT/core/AdePTTransport.icc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,8 +186,6 @@ void AdePTTransport<IntegrationLayer>::Initialize(bool common_data)
return;
}

fIntegrationLayer.InitScoringData();

std::cout << "=== AdePTTransport: initializing transport engine for thread: " << fIntegrationLayer.GetThreadID()
<< std::endl;

Expand Down
7 changes: 3 additions & 4 deletions include/AdePT/integration/AdePTGeant4Integration.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,8 +48,8 @@ public:
static void InitVolAuxData(adeptint::VolAuxData *volAuxData, G4HepEmState *hepEmState, bool trackInAllRegions,
std::vector<std::string> const *gpuRegionNames);

/// @brief Initializes the mapping of VecGeom to G4 volumes for sensitive volumes and their parents
void InitScoringData();
/// @brief Returns a mapping of VecGeom placed volume IDs to Geant4 physical volumes
static void GetPhysicalVolumeMap(std::vector<G4VPhysicalVolume const *> &vecgeomToG4Map);

/// @brief Reconstructs GPU hits on host and calls the user-defined sensitive detector code
void ProcessGPUHit(GPUHit const &hit);
Expand Down Expand Up @@ -83,8 +83,7 @@ private:

void ReturnTrack(adeptint::TrackData const &track, unsigned int trackIndex, int debugLevel) const;

std::unordered_map<size_t,
const G4VPhysicalVolume *> fglobal_vecgeom_to_g4_map; ///< Maps Vecgeom PV IDs to G4 PV IDs
static std::vector<G4VPhysicalVolume const *> fglobal_vecgeom_to_g4_map;
std::unique_ptr<AdePTGeant4Integration_detail::ScoringObjects, AdePTGeant4Integration_detail::Deleter>
fScoringObjects{nullptr};
};
Expand Down
82 changes: 47 additions & 35 deletions src/AdePTGeant4Integration.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,8 +84,42 @@ void Deleter::operator()(ScoringObjects *ptr)

} // namespace AdePTGeant4Integration_detail

std::vector<G4VPhysicalVolume const *> AdePTGeant4Integration::fglobal_vecgeom_to_g4_map;

AdePTGeant4Integration::~AdePTGeant4Integration() {}

void AdePTGeant4Integration::GetPhysicalVolumeMap(std::vector<G4VPhysicalVolume const *> &vecgeomToG4Map)
{
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<void(G4VPhysicalVolume const *, vecgeom::VPlacedVolume const *)> 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 IDs
vecgeomToG4Map.resize(std::max<std::size_t>(vecgeomToG4Map.size(), vg_pvol->id() + 1), nullptr);
vecgeomToG4Map[vg_pvol->id()] = g4_pvol;

// 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];

// VecGeom does not strip pointers from logical volume names
if (std::string(pvol_d->GetLogicalVolume()->GetName()).rfind(g4pvol_d->GetLogicalVolume()->GetName(), 0) != 0)
throw std::runtime_error("Fatal: CheckGeometry: Volume names " +
std::string(pvol_d->GetLogicalVolume()->GetName()) + " and " +
std::string(g4pvol_d->GetLogicalVolume()->GetName()) + " mismatch");
visitGeometry(g4pvol_d, pvol_d);
}
};
visitGeometry(g4world, vecgeomWorld);
}

void AdePTGeant4Integration::CreateVecGeomWorld(std::string filename)
{
// Import the gdml file into VecGeom
Expand All @@ -97,6 +131,9 @@ void AdePTGeant4Integration::CreateVecGeomWorld(std::string filename)
return;
}

// Generate the mapping of VecGeom volume IDs to Geant4 physical volumes
GetPhysicalVolumeMap(fglobal_vecgeom_to_g4_map);

const vecgeom::VPlacedVolume *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld();
if (vecgeomWorld == nullptr) {
std::cerr << "GeoManager vecgeomWorld volume is nullptr" << G4endl;
Expand All @@ -115,6 +152,14 @@ void AdePTGeant4Integration::CreateVecGeomWorld(G4VPhysicalVolume const *physvol
auto conversion = g4vg::convert(physvol);
vecgeom::GeoManager::Instance().SetWorldAndClose(conversion.world);

// Get the mapping of VecGeom volume IDs to Geant4 physical volumes from g4vg
// fglobal_vecgeom_to_g4_map = conversion.physical_volumes;

// TODO: We generate the volume map using a visitor since the one provided by the converter
// appears to have a small number of wrongly mapped IDs
// Generate the mapping of VecGeom volume IDs to Geant4 physical volumes
GetPhysicalVolumeMap(fglobal_vecgeom_to_g4_map);

// EXPECT: we finish with a non-null VecGeom host geometry
vecgeom::VPlacedVolume const *vecgeomWorld = vecgeom::GeoManager::Instance().GetWorld();
if (vecgeomWorld == nullptr) {
Expand Down Expand Up @@ -267,38 +312,6 @@ void AdePTGeant4Integration::InitVolAuxData(adeptint::VolAuxData *volAuxData, G4
visitGeometry(g4world, vecgeomWorld);
}

void AdePTGeant4Integration::InitScoringData()
{
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<void(G4VPhysicalVolume const *, vecgeom::VPlacedVolume const *)> 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 IDs
// Though we only record and reconstruct hits for sensitive volumes, this map needs to store every
// volume in the geometry, as a step may begin in a sensitive volume and end in a non-sensitive one
fglobal_vecgeom_to_g4_map.insert(std::pair<int, const G4VPhysicalVolume *>(vg_pvol->id(), g4_pvol));
// 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];

// VecGeom does not strip pointers from logical volume names
if (std::string(pvol_d->GetLogicalVolume()->GetName()).rfind(g4pvol_d->GetLogicalVolume()->GetName(), 0) != 0)
throw std::runtime_error("Fatal: CheckGeometry: Volume names " +
std::string(pvol_d->GetLogicalVolume()->GetName()) + " and " +
std::string(g4pvol_d->GetLogicalVolume()->GetName()) + " mismatch");
visitGeometry(g4pvol_d, pvol_d);
}
};
visitGeometry(g4world, vecgeomWorld);
}

void AdePTGeant4Integration::ProcessGPUHit(GPUHit const &hit)
{
if (!fScoringObjects) {
Expand Down Expand Up @@ -363,9 +376,8 @@ void AdePTGeant4Integration::FillG4NavigationHistory(vecgeom::NavigationState aN
for (aLevel = 0; aLevel <= aVecGeomLevel; aLevel++) {
// While we are in levels shallower than the history depth, it may be that we already
// have the correct volume in the history
const auto item = fglobal_vecgeom_to_g4_map.find(aNavState.At(aLevel)->id());
pnewvol = item == fglobal_vecgeom_to_g4_map.end() ? nullptr : const_cast<G4VPhysicalVolume *>(item->second);

pnewvol = const_cast<G4VPhysicalVolume *>(fglobal_vecgeom_to_g4_map[aNavState.At(aLevel)->id()]);
assert(pnewvol != nullptr);
if (!pnewvol) throw std::runtime_error("VecGeom volume not found in G4 mapping!");

if (aG4HistoryDepth && (aLevel <= aG4HistoryDepth)) {
Expand Down

0 comments on commit a0ac9bb

Please sign in to comment.