Support sensitive detectors in replica/parameterized volumes (#1649)

* Tweak function signature * Define a unique instance class to enable replicas * Add is_replica test * Update example calo * Add replica test for geant4 * Add tracking tolerance and bumps, and avoid coincident ray * Add reconstruction of replica in test * Support reconstructing in replicas * Add test for level touchable updater * Revert unnecessary test changes and fix units * Fix typo
celeritas-project · Mar 4, 2025 · b02445f · b02445f
1 parent 2348fa6
commit b02445f
Show file tree

Hide file tree

Showing 30 changed files with 2,425 additions and 150 deletions.
diff --git a/src/celeritas/ext/detail/LevelTouchableUpdater.cc b/src/celeritas/ext/detail/LevelTouchableUpdater.cc
@@ -11,6 +11,7 @@
 #include <G4VTouchable.hh>
 
 #include "geocel/GeantGeoUtils.hh"
+#include "geocel/GeoParamsInterface.hh"
 #include "geocel/GeoTraits.hh"
 #include "celeritas/geo/GeoParams.hh"  // IWYU pragma: keep
 #include "celeritas/user/DetectorSteps.hh"
@@ -60,23 +61,23 @@ bool LevelTouchableUpdater::operator()(SpanVolInst ids,
     CELER_EXPECT(touchable);
     CELER_EXPECT(!ids.empty() && ids.front());
 
-    // Update phys_vol_ from geometry and volume instance id
-    phys_vol_.clear();
+    // Update phys_inst_ from geometry and volume instance id
+    phys_inst_.clear();
     for (auto vi_id : ids)
     {
         if (!vi_id)
             break;
-        auto pv = geo_->id_to_pv(vi_id);
+        auto phys_inst = geo_->id_to_geant(vi_id);
         CELER_VALIDATE(
-            pv,
+            phys_inst,
             << "no Geant4 physical volume is attached to volume instance "
             << vi_id.get() << "='" << geo_->volume_instances().at(vi_id)
             << "' (geometry type: " << GeoTraits<GeoParams>::name << ')');
-        phys_vol_.push_back(pv);
+        phys_inst_.push_back(phys_inst);
     }
-    CELER_ASSERT(!phys_vol_.empty());
+    CELER_ASSERT(!phys_inst_.empty());
 
-    set_history(make_span(phys_vol_), nav_hist_.get());
+    set_history(make_span(phys_inst_), nav_hist_.get());
     touchable->UpdateYourself(nav_hist_->GetTopVolume(), nav_hist_.get());
     return true;
 }

diff --git a/src/celeritas/ext/detail/LevelTouchableUpdater.hh b/src/celeritas/ext/detail/LevelTouchableUpdater.hh
@@ -60,7 +60,7 @@ class LevelTouchableUpdater final : public TouchableUpdaterInterface
     // Geometry for doing G4PV translation
     SPConstGeo geo_;
     // Temporary storage for physical volumes
-    std::vector<G4VPhysicalVolume const*> phys_vol_;
+    std::vector<GeantPhysicalInstance> phys_inst_;
     // Temporary history
     std::unique_ptr<G4NavigationHistory> nav_hist_;
 };

diff --git a/src/geocel/GeantGeoUtils.cc b/src/geocel/GeantGeoUtils.cc
@@ -22,10 +22,12 @@
 #include <G4PhysicalVolumeStore.hh>
 #include <G4ReflectionFactory.hh>
 #include <G4RegionStore.hh>
+#include <G4ReplicaNavigation.hh>
 #include <G4SolidStore.hh>
 #include <G4Threading.hh>
 #include <G4TouchableHistory.hh>
 #include <G4TransportationManager.hh>
+#include <G4VPVParameterisation.hh>
 #include <G4VPhysicalVolume.hh>
 #include <G4Version.hh>
 #include <G4ios.hh>
@@ -38,6 +40,7 @@
 #include "corecel/math/Algorithms.hh"
 #include "orange/g4org/Converter.hh"
 
+#include "GeoParamsInterface.hh"
 #include "g4/VisitVolumes.hh"
 
 #include "detail/MakeLabelVector.hh"
@@ -51,6 +54,8 @@ static_assert(G4VERSION_NUMBER
               "CMake-reported Geant4 version does not match installed "
               "<G4Version.hh>: compare to 'corecel/Config.hh'");
 
+using ReplicaId = celeritas::GeantPhysicalInstance::ReplicaId;
+
 namespace celeritas
 {
 namespace
@@ -74,6 +79,25 @@ void free_and_clear(std::vector<T*>* table)
     table->clear();
 }
 
+void update_replica(G4VPhysicalVolume* pv, ReplicaId replica)
+{
+    CELER_EXPECT(pv);
+    CELER_EXPECT(replica);
+    static G4ThreadLocal G4ReplicaNavigation nav_;
+    nav_.ComputeTransformation(replica.get(), pv);
+    pv->SetCopyNo(replica.get());
+}
+
+void update_parameterised(G4VPhysicalVolume* pv, ReplicaId replica)
+{
+    CELER_EXPECT(pv);
+    CELER_EXPECT(replica);
+    auto* param = pv->GetParameterisation();
+    CELER_ASSERT(param);
+    param->ComputeTransformation(replica.get(), pv);
+    pv->SetCopyNo(replica.get());
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace
 
@@ -194,6 +218,16 @@ G4VPhysicalVolume const* geant_world_volume()
     return nav->GetWorldVolume();
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Whether a physical volume is parameterized or replicated.
+ */
+bool is_replica(G4VPhysicalVolume const& pv)
+{
+    auto vt = pv.VolumeType();
+    return (vt == EVolume::kReplica || vt == EVolume::kParameterised);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Find Geant4 logical volumes corresponding to a list of names.
@@ -261,9 +295,9 @@ std::vector<Label> make_logical_vol_labels(G4VPhysicalVolume const& world)
         },
         world);
 
-    return detail::make_label_vector<G4LogicalVolume>(
+    return detail::make_label_vector<G4LogicalVolume const*>(
         std::move(names),
-        [](G4LogicalVolume const& lv) { return lv.GetInstanceID(); });
+        [](G4LogicalVolume const* lv) { return lv->GetInstanceID(); });
 }
 
 //---------------------------------------------------------------------------//
@@ -298,21 +332,20 @@ std::vector<Label> make_physical_vol_labels(G4VPhysicalVolume const& world)
         },
         world);
 
-    return detail::make_label_vector<G4VPhysicalVolume>(
+    return detail::make_label_vector<G4VPhysicalVolume const*>(
         std::move(names),
-        [](G4VPhysicalVolume const& lv) { return lv.GetInstanceID(); });
+        [](G4VPhysicalVolume const* pv) { return pv->GetInstanceID(); });
 }
 
 //---------------------------------------------------------------------------//
 /*!
  * Update a nav history to match the given pv stack.
  *
- * \warning The stack cannot have a parameterized/replicated volume.
- *
  * \note The stack should have the same semantics as \c LevelId, i.e. the
  * initial entry is the "most global" level.
  */
-void set_history(Span<G4VPhysicalVolume const*> stack, G4NavigationHistory* nav)
+void set_history(Span<GeantPhysicalInstance const> stack,
+                 G4NavigationHistory* nav)
 {
     CELER_EXPECT(!stack.empty());
     CELER_EXPECT(std::all_of(stack.begin(), stack.end(), LogicalTrue{}));
@@ -327,7 +360,8 @@ void set_history(Span<G4VPhysicalVolume const*> stack, G4NavigationHistory* nav)
          level != end_level;
          ++level)
     {
-        if (nav->GetVolume(level) != stack[level])
+        auto* nav_pv = nav->GetVolume(level);
+        if (nav_pv != stack[level].pv || stack[level].replica)
         {
             break;
         }
@@ -338,7 +372,7 @@ void set_history(Span<G4VPhysicalVolume const*> stack, G4NavigationHistory* nav)
         // Top level disagrees (rare? should always be world):
         // reset to top level
         nav->Reset();
-        nav->SetFirstEntry(const_cast<G4VPhysicalVolume*>(stack[0]));
+        nav->SetFirstEntry(const_cast<G4VPhysicalVolume*>(stack[0].pv));
         ++level;
     }
     else if (level < nav_stack_size())
@@ -348,18 +382,32 @@ void set_history(Span<G4VPhysicalVolume const*> stack, G4NavigationHistory* nav)
         CELER_ASSERT(nav_stack_size() == level);
     }
 
-    // Add all remaining levels
+    // Add all remaining levels: see G4Navigator::LocateGLobalPoint
     for (auto end_level = stack.size(); level != end_level; ++level)
     {
-        G4VPhysicalVolume const* pv = stack[level];
-        constexpr auto volume_type = EVolume::kNormal;
-        CELER_VALIDATE(pv->VolumeType() == volume_type,
-                       << "sensitive detectors inside of "
-                          "replica/parameterized volumes are not supported: '"
-                       << pv->GetName() << "' inside "
-                       << PrintableNavHistory{nav});
-        nav->NewLevel(
-            const_cast<G4VPhysicalVolume*>(pv), volume_type, pv->GetCopyNo());
+        auto* pv = const_cast<G4VPhysicalVolume*>(stack[level].pv);
+        auto vol_type = pv->VolumeType();
+        auto replica = stack[level].replica;
+        switch (vol_type)
+        {
+            case EVolume::kNormal:
+                CELER_ASSERT(!replica);
+                replica = id_cast<ReplicaId>(pv->GetCopyNo());
+                break;
+            case EVolume::kReplica:
+                update_replica(pv, replica);
+                break;
+            case EVolume::kParameterised:
+                update_parameterised(pv, replica);
+                break;
+            default:
+                CELER_LOG_LOCAL(error)
+                    << R"(Encountered abnormal Geant4 volume inside navigation history: )"
+                    << pv->GetName() << "' inside "
+                    << PrintableNavHistory{nav};
+                break;
+        }
+        nav->NewLevel(pv, vol_type, replica.get());
     }
 
     CELER_ENSURE(nav_stack_size() == stack.size());

diff --git a/src/geocel/GeantGeoUtils.hh b/src/geocel/GeantGeoUtils.hh
@@ -32,6 +32,7 @@ class G4VTouchable;
 
 namespace celeritas
 {
+struct GeantPhysicalInstance;
 //---------------------------------------------------------------------------//
 #if CELERITAS_GEANT4_VERSION >= 0x0b0200
 //! Version-independent typedef to Geant4 touchable history
@@ -78,6 +79,10 @@ Span<G4LogicalVolume*> geant_logical_volumes();
 // Get the world volume if the geometry has been set up
 G4VPhysicalVolume const* geant_world_volume();
 
+//---------------------------------------------------------------------------//
+// Whether the volume is a replica/parameterization
+bool is_replica(G4VPhysicalVolume const&);
+
 //---------------------------------------------------------------------------//
 // Find Geant4 logical volumes corresponding to a list of names
 std::unordered_set<G4LogicalVolume const*>
@@ -90,8 +95,8 @@ std::vector<Label> make_logical_vol_labels(G4VPhysicalVolume const& world);
 std::vector<Label> make_physical_vol_labels(G4VPhysicalVolume const& world);
 
 //---------------------------------------------------------------------------//
-// Update a nav history to match the given pv stack
-void set_history(Span<G4VPhysicalVolume const*> stack,
+// Update a nav history to match the given volume instance stack
+void set_history(Span<GeantPhysicalInstance const> stack,
                  G4NavigationHistory* nav);
 
 //---------------------------------------------------------------------------//

diff --git a/src/geocel/GeoParamsInterface.hh b/src/geocel/GeoParamsInterface.hh
@@ -6,6 +6,8 @@
 //---------------------------------------------------------------------------//
 #pragma once
 
+#include <utility>
+
 #include "corecel/Macros.hh"
 #include "corecel/cont/LabelIdMultiMap.hh"  // IWYU pragma: export
 #include "corecel/cont/Span.hh"  // IWYU pragma: export
@@ -19,6 +21,31 @@ class G4VPhysicalVolume;
 
 namespace celeritas
 {
+//---------------------------------------------------------------------------//
+/*!
+ * Unique placement/replica of a Geant4 physical volume.
+ *
+ * This should correspond to a \c VolumeInstanceId and be a unique
+ * instantiation of a Geant4 physical volume (PV). Some Geant4 PVs are
+ * "parameterised" or "replica" types, which allows a single instance to be
+ * mutated at runtime to form a sort of array.
+ * If the pointed-to physical volume is *not* a replica/parameterised volume,
+ * \c replica is false. Otherwise, it corresponds to the PV's copy number,
+ * which can be used to reconstruct the placed volume instance.
+ */
+struct GeantPhysicalInstance
+{
+    using ReplicaId = OpaqueId<struct Replica_>;
+
+    //! Geant4 physical volume
+    G4VPhysicalVolume const* pv{nullptr};
+    //! Replica/parameterisation instance
+    ReplicaId replica;
+
+    //! False if no PV is associated
+    constexpr explicit operator bool() const { return static_cast<bool>(pv); }
+};
+
 //---------------------------------------------------------------------------//
 /*!
  * Interface class for accessing host geometry metadata.
@@ -47,7 +74,7 @@ class GeoParamsInterface
     //! Outer bounding box of geometry
     virtual BBox const& bbox() const = 0;
 
-    //! Maximum nested scene/volume depth
+    //! Maximum nested volume instance depth
     virtual LevelId::size_type max_depth() const = 0;
 
     //// VOLUMES ////
@@ -62,7 +89,7 @@ class GeoParamsInterface
     virtual VolumeId find_volume(G4LogicalVolume const* volume) const = 0;
 
     //! Get the Geant4 PV corresponding to a volume instance
-    virtual G4VPhysicalVolume const* id_to_pv(VolumeInstanceId id) const = 0;
+    virtual GeantPhysicalInstance id_to_geant(VolumeInstanceId id) const = 0;
 
     //// DEPRECATED: remove in v0.6 ////
 

diff --git a/src/geocel/detail/MakeLabelVector.hh b/src/geocel/detail/MakeLabelVector.hh
@@ -25,9 +25,9 @@ namespace detail
  * ID.
  */
 template<class T, class GetId>
-std::vector<Label> make_label_vector(
-    std::unordered_map<std::string, std::vector<T const*>> const& names,
-    GetId&& get_id)
+std::vector<Label>
+make_label_vector(std::unordered_map<std::string, std::vector<T>> const& names,
+                  GetId&& get_id)
 {
     std::vector<Label> result;
     auto add_label = [&result](std::size_t id, Label&& label) {
@@ -45,15 +45,15 @@ std::vector<Label> make_label_vector(
         {
             // Label is just the name since this item is unique
             CELER_ASSERT(items.front());
-            add_label(get_id(*items.front()), Label{name});
+            add_label(get_id(items.front()), Label{name});
             continue;
         }
 
         // Set labels in increasing order
         for (auto idx : range(items.size()))
         {
             CELER_ASSERT(items[idx]);
-            add_label(get_id(*items[idx]), Label{name, std::to_string(idx)});
+            add_label(get_id(items[idx]), Label{name, std::to_string(idx)});
         }
     }
     return result;

diff --git a/src/geocel/g4/GeantGeoParams.cc b/src/geocel/g4/GeantGeoParams.cc
@@ -165,20 +165,39 @@ VolumeId GeantGeoParams::find_volume(G4LogicalVolume const* volume) const
 /*!
  * Get the Geant4 physical volume corresponding to a volume instance ID.
  *
- * If the input ID is false, a null pointer will be returned.
+ * \warning For Geant4 parameterised/replicated volumes, external state (e.g.
+ * the local navigation) \em must be used in concert with this method: i.e.,
+ * navigation on the current thread needs to have just "visited" the instance.
+ *
+ * \todo Create our own volume instance mapping for Geant4, where
+ * VolumeInstanceId corresponds to a (G4PV*, int) pair rather than just G4PV*
+ * (which in Geant4 is not sufficiently unique).
  */
-G4VPhysicalVolume const* GeantGeoParams::id_to_pv(VolumeInstanceId id) const
+GeantPhysicalInstance GeantGeoParams::id_to_geant(VolumeInstanceId id) const
 {
     CELER_EXPECT(!id || id < vol_instances_.size());
     if (!id)
     {
-        return nullptr;
+        return {};
     }
 
     G4PhysicalVolumeStore* pv_store = G4PhysicalVolumeStore::GetInstance();
     auto index = id.unchecked_get() - pv_offset_;
     CELER_ASSERT(index < pv_store->size());
-    return (*pv_store)[index];
+
+    GeantPhysicalInstance result;
+    result.pv = (*pv_store)[index];
+    CELER_ASSERT(result.pv);
+    if (result.pv->VolumeType() != EVolume::kNormal)
+    {
+        auto copy_no = result.pv->GetCopyNo();
+        // NOTE: if this line fails, Geant4 may be returning uninitialized
+        // memory on the local thread.
+        CELER_ASSERT(copy_no >= 0 && copy_no < result.pv->GetMultiplicity());
+        result.replica = id_cast<GeantPhysicalInstance::ReplicaId>(copy_no);
+    }
+
+    return result;
 }
 
 //---------------------------------------------------------------------------//
@@ -187,7 +206,7 @@ G4VPhysicalVolume const* GeantGeoParams::id_to_pv(VolumeInstanceId id) const
  *
  * If the input volume ID is unassigned, a null pointer will be returned.
  */
-G4LogicalVolume const* GeantGeoParams::id_to_lv(VolumeId id) const
+G4LogicalVolume const* GeantGeoParams::id_to_geant(VolumeId id) const
 {
     CELER_EXPECT(!id || id < volumes_.size());
     if (!id)