Add support for multiple mother volumes to vertex confinement

include/RMGNavigationTools.hh

@@ -16,6 +16,7 @@
 #ifndef _RMG_NAVIGATION_TOOLS_HH_
 #define _RMG_NAVIGATION_TOOLS_HH_
 
+#include <set>
 #include <string>
 
 #include "G4LogicalVolume.hh"
@@ -30,6 +31,7 @@ namespace RMGNavigationTools {
   G4LogicalVolume* FindLogicalVolume(std::string name);
   G4VPhysicalVolume* FindPhysicalVolume(std::string name, int copy_nr = 0);
   G4VPhysicalVolume* FindDirectMother(G4VPhysicalVolume* volume);
+  std::set<G4VPhysicalVolume*> FindDirectMothers(G4VPhysicalVolume* volume);
 
   void PrintListOfLogicalVolumes();
   void PrintListOfPhysicalVolumes();

include/RMGVertexConfinement.hh

@@ -18,9 +18,11 @@
 
 #include <chrono>
 #include <optional>
+#include <queue>
 #include <regex>
 #include <vector>
 
+#include "G4GenericMessenger.hh"
 #include "G4RotationMatrix.hh"
 #include "G4ThreeVector.hh"
 #include "G4Transform3D.hh"
@@ -30,7 +32,6 @@
 
 class G4VPhysicalVolume;
 class G4VSolid;
-class G4GenericMessenger;
 class RMGVertexConfinement : public RMGVVertexGenerator {
 
   public:
@@ -85,7 +86,8 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
     struct SampleableObject {
 
         SampleableObject() = default;
-        SampleableObject(G4VPhysicalVolume* v, G4RotationMatrix r, G4ThreeVector t, G4VSolid* s);
+        SampleableObject(G4VPhysicalVolume* v, G4RotationMatrix r, G4ThreeVector t, G4VSolid* s,
+            bool cc = true);
         // NOTE: G4 volume/solid pointers should be fully owned by G4, avoid trying to delete them
         ~SampleableObject() = default;
 
@@ -111,9 +113,10 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
         [[nodiscard]] size_t size() const { return data.size(); }
         SampleableObject& at(size_t i) { return data.at(i); }
         void emplace_back(G4VPhysicalVolume* v, const G4RotationMatrix& r, const G4ThreeVector& t,
-            G4VSolid* s);
+            G4VSolid* s, bool cc = true);
         inline void push_back(const SampleableObject& obj) {
-          this->emplace_back(obj.physical_volume, obj.rotation, obj.translation, obj.sampling_solid);
+          this->emplace_back(obj.physical_volume, obj.rotation, obj.translation, obj.sampling_solid,
+              obj.containment_check);
         }
         [[nodiscard]] inline bool empty() const { return data.empty(); }
         inline SampleableObject& back() { return data.back(); }
@@ -129,6 +132,15 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
 
   private:
 
+    struct VolumeTreeEntry {
+        G4VPhysicalVolume* physvol;
+
+        G4ThreeVector vol_global_translation; // origin
+        G4RotationMatrix vol_global_rotation; // identity
+        std::vector<G4RotationMatrix> partial_rotations;
+        std::vector<G4ThreeVector> partial_translations;
+    };
+
     void InitializePhysicalVolumes();
     void InitializeGeometricalVolumes();
     bool ActualGenerateVertex(G4ThreeVector& v);
@@ -142,6 +154,7 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
 
     SamplingMode fSamplingMode = kUnionAll;
     bool fOnSurface = false;
+    bool fForceContainmentCheck = false;
 
     // counters used for the current run.
     long fTrials = 0;

src/RMGNavigationTools.cc

@@ -53,6 +53,18 @@ G4LogicalVolume* RMGNavigationTools::FindLogicalVolume(std::string name) {
 
 G4VPhysicalVolume* RMGNavigationTools::FindDirectMother(G4VPhysicalVolume* volume) {
 
+  auto ancestors = RMGNavigationTools::FindDirectMothers(volume);
+
+  if (ancestors.size() != 1) {
+    RMGLog::Out(RMGLog::fatal, "Could not find a unique direct mother volume, ",
+        "this cannot be! Returning first ancestor in the list");
+  }
+
+  return *ancestors.begin();
+}
+
+std::set<G4VPhysicalVolume*> RMGNavigationTools::FindDirectMothers(G4VPhysicalVolume* volume) {
+
   std::set<G4VPhysicalVolume*> ancestors;
   for (const auto& v : *G4PhysicalVolumeStore::GetInstance()) {
     if (v->GetLogicalVolume()->IsAncestor(volume)) ancestors.insert(v);
@@ -78,12 +90,7 @@ G4VPhysicalVolume* RMGNavigationTools::FindDirectMother(G4VPhysicalVolume* volum
     else it++;
   }
 
-  if (ancestors.size() != 1) {
-    RMGLog::Out(RMGLog::error, "Could not find a unique direct mother volume, ",
-        "this cannot be! Returning first ancestor in the list");
-  }
-
-  return *ancestors.begin();
+  return ancestors;
 }
 
 void RMGNavigationTools::PrintListOfLogicalVolumes() {

src/RMGVertexConfinement.cc

@@ -44,10 +44,8 @@
 //  - physical volumes get always a bounding box assigned, but at later time
 //  - purely geometrical volumes only have the sampling_solid member defined
 RMGVertexConfinement::SampleableObject::SampleableObject(G4VPhysicalVolume* v, G4RotationMatrix r,
-    G4ThreeVector t, G4VSolid* s)
-    :
-
-      rotation(r), translation(t) {
+    G4ThreeVector t, G4VSolid* s, bool cc)
+    : rotation(r), translation(t), containment_check(cc) {
 
   // at least one volume must be specified
   if (!v and !s) RMGLog::Out(RMGLog::error, "Invalid pointers given to constructor");
@@ -115,9 +113,9 @@ bool RMGVertexConfinement::SampleableObjectCollection::IsInside(const G4ThreeVec
 }
 
 void RMGVertexConfinement::SampleableObjectCollection::emplace_back(G4VPhysicalVolume* v,
-    const G4RotationMatrix& r, const G4ThreeVector& t, G4VSolid* s) {
+    const G4RotationMatrix& r, const G4ThreeVector& t, G4VSolid* s, bool cc) {
 
-  this->data.emplace_back(v, r, t, s);
+  this->data.emplace_back(v, r, t, s, cc);
 
   const auto& _v = this->data.back().volume;
   const auto& _s = this->data.back().surface;
@@ -164,8 +162,9 @@ void RMGVertexConfinement::InitializePhysicalVolumes() {
         // do not specify a bounding solid at this stage
         fPhysicalVolumes.emplace_back(*it, G4RotationMatrix(), G4ThreeVector(), nullptr);
 
-        RMGLog::OutFormat(RMGLog::detail, " · '{}[{}]', mass = {}", (*it)->GetName().c_str(),
-            (*it)->GetCopyNo(), std::string(G4BestUnit(fPhysicalVolumes.data.back().volume, "Mass")));
+        RMGLog::OutFormat(RMGLog::detail, " · '{}[{}]', volume = {}", (*it)->GetName().c_str(),
+            (*it)->GetCopyNo(),
+            std::string(G4BestUnit(fPhysicalVolumes.data.back().volume, "Volume")));
 
         found = true;
       }
@@ -193,6 +192,7 @@ void RMGVertexConfinement::InitializePhysicalVolumes() {
   // now inspect the solids the physical volumes refer to and configure the
   // appropriate sampling strategy, i.e. setting the sampling solid and
   // containment check flag
+  std::vector<SampleableObject> new_obj_from_inspection;
   for (auto&& el : fPhysicalVolumes.data) {
 
     RMGLog::OutFormatDev(RMGLog::debug, "Inspecting volume '{}'", el.physical_volume->GetName());
@@ -260,44 +260,83 @@ void RMGVertexConfinement::InitializePhysicalVolumes() {
               (lim_max.getZ() - lim_min.getZ()) / 2);
     } // sampling_solid and containment_check must hold a valid value at this point
 
+
     // determine solid transformation w.r.t. world volume reference
 
-    G4ThreeVector vol_global_translation; // origin
-    G4RotationMatrix vol_global_rotation; // identity
-    std::vector<G4RotationMatrix> partial_rotations;
-    std::vector<G4ThreeVector> partial_translations;
+    // found paths to the mother volume.
+    std::vector<VolumeTreeEntry> trees;
+
+    // queue for paths to the mother volume that still have to be searched.
+    std::queue<VolumeTreeEntry> q;
+    q.push({el.physical_volume});
 
-    for (auto v = el.physical_volume; v != world_volume; v = RMGNavigationTools::FindDirectMother(v)) {
+    for (; !q.empty(); q.pop()) {
+      auto v = q.front();
 
-      if (!v)
-        RMGLog::Out(RMGLog::fatal, "nullptr detected in loop condition, this is unexpected. ",
-            "Blame RMGPhysVolNavigator::FindDirectMother?");
+      if (!v.physvol)
+        RMGLog::OutDev(RMGLog::fatal, "nullptr detected in loop condition, this is unexpected. ",
+            "Blame RMGNavigationTools::FindDirectMother?");
 
-      partial_rotations.push_back(v->GetObjectRotationValue());
-      partial_translations.push_back(v->GetObjectTranslation());
+      v.partial_rotations.push_back(v.physvol->GetObjectRotationValue());
+      v.partial_translations.push_back(v.physvol->GetObjectTranslation());
 
-      vol_global_rotation = partial_rotations.back() * vol_global_rotation;
+      v.vol_global_rotation = v.partial_rotations.back() * v.vol_global_rotation;
+
+      for (auto m : RMGNavigationTools::FindDirectMothers(v.physvol)) {
+        if (m != world_volume) {
+          auto v_m = VolumeTreeEntry(v); // create a copy of the current helper object.
+          v_m.physvol = m;
+          q.push(v_m);
+        } else { // we finished that branch!
+          trees.push_back(v);
+        }
+      }
     }
-    // world volume not included in loop
-    partial_translations.emplace_back(); // origin
-    partial_rotations.emplace_back();    // identity
-
-    // partial_rotations[0] and partial_translations[0] refer to the target
-    // volume partial_rotations[1] and partial_translations[1], to the direct
-    // mother, etc. It is necessary to rotate with respect to the frame of the
-    // mother. If there are no rotations (or only the target volume is
-    // rotated): rotations are identity matrices and vol_global_translation =
-    // sum(partial_translations)
-    for (size_t i = 0; i < partial_translations.size() - 1; i++) {
-      G4ThreeVector tmp = partial_translations[i];
-      for (size_t j = i + 1; j < partial_rotations.size() - 1; j++) { tmp *= partial_rotations[j]; }
-      vol_global_translation += tmp;
+
+    RMGLog::OutFormatDev(RMGLog::debug, "Found {} ways to reach world volume from {}", trees.size(),
+        el.physical_volume->GetName());
+
+    // finalize all found paths to the mother volume.
+    for (auto&& v : trees) {
+      // world volume not included in loop
+      v.partial_translations.emplace_back(); // origin
+      v.partial_rotations.emplace_back();    // identity
+
+      // partial_rotations[0] and partial_translations[0] refer to the target
+      // volume partial_rotations[1] and partial_translations[1], to the direct
+      // mother, etc. It is necessary to rotate with respect to the frame of the
+      // mother. If there are no rotations (or only the target volume is
+      // rotated): rotations are identity matrices and vol_global_translation =
+      // sum(partial_translations)
+      for (size_t i = 0; i < v.partial_translations.size() - 1; i++) {
+        G4ThreeVector tmp = v.partial_translations[i];
+        for (size_t j = i + 1; j < v.partial_rotations.size() - 1; j++) {
+          tmp *= v.partial_rotations[j];
+        }
+        v.vol_global_translation += tmp;
+      }
     }
 
-    // assign transformation to sampling solid
-    el.rotation = vol_global_rotation;
-    el.translation = vol_global_translation;
+    if (trees.empty())
+      RMGLog::OutDev(RMGLog::fatal, "No path to world volume found, that should not be!");
+    // assign first found transformation to current sampling solid
+    el.rotation = trees[0].vol_global_rotation;
+    el.translation = trees[0].vol_global_translation;
+
+    // we found more than one path to the mother volume. Append them to the list of sampling
+    // solids separately.
+    for (size_t i = 1; i < trees.size(); ++i) {
+      const auto& v = trees.at(i);
+      new_obj_from_inspection.emplace_back(el.physical_volume, v.vol_global_rotation,
+          v.vol_global_translation, el.sampling_solid, el.containment_check);
+    }
   }
+
+  // finally add all newly found sampling solids (this os not done directly above as this
+  // invalidates old iterators).
+  for (const auto& s : new_obj_from_inspection) { fPhysicalVolumes.push_back(s); }
+
+  RMGLog::Out(RMGLog::detail, "Sampling from ", fPhysicalVolumes.size(), " physical volumes");
 }
 
 void RMGVertexConfinement::InitializeGeometricalVolumes() {
@@ -350,6 +389,7 @@ void RMGVertexConfinement::Reset() {
   fGeomVolumeSolids.clear();
   fSamplingMode = RMGVertexConfinement::kUnionAll;
   fOnSurface = false;
+  fForceContainmentCheck = false;
 }
 
 bool RMGVertexConfinement::GenerateVertex(G4ThreeVector& vertex) {
@@ -421,6 +461,14 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) {
         } else {
           vertex = choice.translation +
                    choice.rotation * RMGGeneratorUtil::rand(choice.sampling_solid, fOnSurface);
+          if (fForceContainmentCheck) {
+            auto is_inside = physical_first ? fPhysicalVolumes.IsInside(vertex)
+                                            : fGeomVolumeSolids.IsInside(vertex);
+            if (!is_inside)
+              RMGLog::OutDev(RMGLog::error,
+                  "Generated vertex not inside sampling volumes (forced containment check): ",
+                  vertex / CLHEP::cm, " cm");
+          }
         }
 
         // is it also in the other volume class (geometrical/physical)?
@@ -493,6 +541,11 @@ bool RMGVertexConfinement::ActualGenerateVertex(G4ThreeVector& vertex) {
           vertex = choice.translation +
                    choice.rotation * RMGGeneratorUtil::rand(choice.sampling_solid, fOnSurface);
           RMGLog::OutDev(RMGLog::debug, "Generated vertex: ", vertex / CLHEP::cm, " cm");
+          if (fForceContainmentCheck && !all_volumes.IsInside(vertex)) {
+            RMGLog::OutDev(RMGLog::error,
+                "Generated vertex not inside sampling volumes (forced containment check): ",
+                vertex / CLHEP::cm, " cm");
+          }
         }
 
         RMGLog::OutDev(RMGLog::debug, "Found good vertex ", vertex / CLHEP::cm, " cm", " after ",
@@ -609,6 +662,15 @@ void RMGVertexConfinement::DefineCommands() {
       .SetStates(G4State_PreInit, G4State_Idle)
       .SetToBeBroadcasted(true);
 
+  fMessengers.back()
+      ->DeclareProperty("ForceContainmentCheck", fForceContainmentCheck)
+      .SetGuidance("If true (or omitted argument), perform a containment check even after sampling "
+                   "from a natively sampleable object. This is only an extra sanity check that does"
+                   "not alter the behaviour.")
+      .SetParameterName("flag", true)
+      .SetStates(G4State_PreInit, G4State_Idle)
+      .SetToBeBroadcasted(false);
+
   fMessengers.push_back(
       std::make_unique<G4GenericMessenger>(this, "/RMG/Generator/Confinement/Physical/",
           "Commands for setting physical volumes up for primary confinement"));

tests/confinement/CMakeLists.txt

@@ -20,7 +20,8 @@ foreach(_mac ${_macros})
                                                test-out.root -- macros/${_mac})
   set_tests_properties(confinement-mt/${_mac} PROPERTIES LABELS mt)
 
-  add_test(NAME confinement-vis/${_mac} COMMAND remage-cli -g gdml/geometry.gdml -o test-out.root
-                                                -- macros/vis.mac macros/${_mac})
+  add_test(NAME confinement-vis/${_mac}
+           COMMAND remage-cli -g gdml/geometry.gdml -o test-out.root -- macros/_vis.mac
+                   macros/${_mac} macros/_vis-export.mac)
   set_tests_properties(confinement-vis/${_mac} PROPERTIES LABELS vis)
 endforeach()

tests/confinement/macros/_vis-export.mac

@@ -0,0 +1 @@
+/vis/ogl/export {export-fn}

tests/confinement/macros/complex-volume.mac

@@ -1,6 +1,7 @@
-/control/execute macros/init.mac
+/control/execute macros/_init.mac
 
 /RMG/Generator/Confine Volume
+/RMG/Generator/Confinement/ForceContainmentCheck true
 
 /RMG/Generator/Confinement/Physical/AddVolume BoxWithHole
 /RMG/Generator/Confinement/Physical/AddVolume BoxAndOrb
@@ -14,4 +15,4 @@
 
 /run/beamOn 5000
 
-/vis/ogl/export complex-volume.pdf
+/control/alias export-fn complex-volume.pdf

tests/confinement/macros/geometrical-and-physical.mac

@@ -1,6 +1,7 @@
-/control/execute macros/init.mac
+/control/execute macros/_init.mac
 
 /RMG/Generator/Confine Volume
+/RMG/Generator/Confinement/ForceContainmentCheck true
 
 /RMG/Generator/Confinement/SamplingMode IntersectPhysicalWithGeometrical
 
@@ -19,4 +20,4 @@
 
 /run/beamOn 2000
 
-/vis/ogl/export geometrical-and-physical.pdf
+/control/alias export-fn eometrical-and-physical.pdf

tests/confinement/macros/geometrical-or-physical.mac

@@ -1,6 +1,7 @@
-/control/execute macros/init.mac
+/control/execute macros/_init.mac
 
 /RMG/Generator/Confine Volume
+/RMG/Generator/Confinement/ForceContainmentCheck true
 
 # /RMG/Generator/Confinement/SampleOnSurface
 
@@ -20,4 +21,4 @@
 
 /run/beamOn 2000
 
-/vis/ogl/export geometrical-or-physical.pdf
+/control/alias export-fn geometrical-or-physical.pdf

tests/confinement/macros/geometrical.mac

@@ -1,6 +1,7 @@
-/control/execute macros/init.mac
+/control/execute macros/_init.mac
 
 /RMG/Generator/Confine Volume
+/RMG/Generator/Confinement/ForceContainmentCheck true
 
 # /RMG/Generator/Confinement/SampleOnSurface
 
@@ -35,4 +36,4 @@
 
 /run/beamOn 5000
 
-/vis/ogl/export geometrical.pdf
+/control/alias export-fn geometrical.pdf