Belle II Software light-2406-ragdoll
OverlapCheckerModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <geometry/modules/overlapchecker/OverlapCheckerModule.h>
10#include <geometry/GeometryManager.h>
11#include <framework/gearbox/Unit.h>
12
13#include "G4VPhysicalVolume.hh"
14#include "G4LogicalVolume.hh"
15#include "G4UnitsTable.hh"
16#include "G4VExceptionHandler.hh"
17#include "G4StateManager.hh"
18#include "G4AffineTransform.hh"
19#include "G4VSolid.hh"
20
21#include <regex>
22#include <utility>
23
24using namespace Belle2;
25
26REG_MODULE(OverlapChecker);
27
28namespace {
30 class OverlapHandler: public G4VExceptionHandler {
31 public:
33 explicit OverlapHandler(std::function<void(const std::string&)> callback): m_callback(std::move(callback)) {}
35 bool Notify(const char*, const char*, G4ExceptionSeverity, const char* description) override
36 {
37 m_callback(description);
38 return false;
39 }
40 private:
42 std::function<void(const std::string&)> m_callback;
43 };
44}
45
47{
48 //Set module properties and the description
49 setDescription("Checks the geometry for overlaps.");
50
51 //Parameter definition
52 addParam("points", m_points, R"DOC(
53Number of test points we will generate on randomly on the surface of each geometry volume and then check for all of them that
54
551. The points are inside the mother volume
562. The points are not inside any neighbor volume
57
58The higher the number the more precise the check for overlaps becomes, and the slower it gets.
59See also https://questions.belle2.org/question/7264/ )DOC",
60 m_points);
61 addParam("tolerance", m_tolerance, "Tolerance of overlap check.", m_tolerance);
62 addParam("maxErrors", m_maxErrors, "Number of overlap errors per volume before continuing with next volume", m_maxErrors);
63 addParam("maxDepth", m_maxDepth, "Maximum depth to go into the geometry tree, 0 means no maximum", m_maxDepth);
64 addParam("prefix", m_prefix, "Prefix path, only volumes starting with the given path are checked", m_prefix);
65}
66
68{
69 m_displayData.registerInDataStore();
70 if (m_maxErrors <= 0) m_maxErrors = m_points + 1;
71}
72
74{
75 if (!m_displayData.isValid()) m_displayData.create();
76 //Check the geometry tree for overlaps
77 G4UnitDefinition::BuildUnitsTable();
78 G4VPhysicalVolume* volume = geometry::GeometryManager::getInstance().getTopVolume();
79 if (!volume) {
80 B2ERROR("No geometry found. => Add the Geometry module to the path before the OverlapChecker module.");
81 return;
82 }
83 m_seen.clear();
84 // reset the navigation history to be in the top level volume
85 m_nav.Reset();
86 // remember the existing G4Exception handler and set our own one to find the position of the intersection
87 G4VExceptionHandler* old = G4StateManager::GetStateManager()->GetExceptionHandler();
88 // set our own exception handler
89 OverlapHandler handler([&](const std::string & message) { handleOverlap(message); });
90 G4StateManager::GetStateManager()->SetExceptionHandler(&handler);
91 // now check for overlaps
92 checkVolume(volume, "");
93 // and reset the exception handler
94 G4StateManager::GetStateManager()->SetExceptionHandler(old);
95
96 //Print the list of found overlaps
97 for (const auto& m_overlap : m_overlaps) {
98 B2ERROR("Overlaps detected for " << m_overlap);
99 }
100}
101
102void OverlapCheckerModule::handleOverlap(const std::string& geant4Message)
103{
104 // ok, let's handle the Overlap message
105 G4VPhysicalVolume* volume = m_nav.GetTopVolume();
106 m_nav.BackLevel();
107 B2ERROR(geant4Message);
108 std::regex r(R"((mother)?\s*local point \‍(([-+0-9eE.]+),([-+0-9eE.]+),([-+0-9eE.]+)\))");
109 std::smatch m;
110 if (std::regex_search(geant4Message, m, r)) {
111 G4ThreeVector posLocal(std::atof(m[2].str().c_str()), std::atof(m[3].str().c_str()), std::atof(m[4].str().c_str()));
112 if (m[1].length() == 0) {
113 // Damn you Geant4, giving me the coordinates in the coordinate system of
114 // the sister volume which I don't know exactly is almost useless. Aaah
115 // darn it, let's check all sisters, see if they have the correct name
116 // and if so check if the point transformed to the mother system is
117 // inside the one were checking. So let's remember the transformation
118 // relative to the top volume of our current volume.
119 G4AffineTransform trans_volume(volume->GetRotation(), volume->GetTranslation());
120 trans_volume.Invert();
121 // And the solid
122 G4VSolid* solid = volume->GetLogicalVolume()->GetSolid();
123 // Now look for the name of the intersecting volume in the exception message (jaaaay)
124 std::regex nameRegex("with (.*) volume's");
125 std::smatch nameMatch;
126 // And if we can find the name we can draw the intersection
127 if (std::regex_search(geant4Message, nameMatch, nameRegex)) {
128 const std::string& name = nameMatch[1].str();
129 // By looping over all sisters
130 for (size_t i = 0; i < volume->GetMotherLogical()->GetNoDaughters(); ++i) {
131 G4VPhysicalVolume* sister = volume->GetMotherLogical()->GetDaughter(i);
132 // ignoring the ones which don't match the name
133 std::string sisterName{(sister->GetName()).data()};
134 if (name != sisterName) continue;
135 // now transform the point into the coordinate system of the volume we actually look at
136 G4AffineTransform trans_sister(sister->GetRotation(), sister->GetTranslation());
137 G4ThreeVector posMother = trans_sister.TransformPoint(posLocal);
138 G4ThreeVector posSister = trans_volume.TransformPoint(posMother);
139 // and check if the point is inside our volume
140 if (solid->Inside(posSister) != kOutside) {
141 // if so this is the right volume which is intersecting.
142 B2INFO("Found intersecting volume " << sister->GetName() << "." << sister->GetCopyNo());
143 // so add the point to the display data
144 G4AffineTransform t = m_nav.GetTopTransform().Inverse();
145 G4ThreeVector global = t.TransformPoint(posMother);
146 m_displayData->addPoint(geant4Message, ROOT::Math::XYZVector(global[0], global[1], global[2]) * Unit::mm);
147 }
148 }
149 } else {
150 B2ERROR("Could not find name of intersecting volume");
151 }
152 } else {
153 // in case of overlap with mother the life is so much simpler: just convert the point to global and save it.
154 G4AffineTransform t = m_nav.GetTopTransform().Inverse();
155 G4ThreeVector global = t.TransformPoint(posLocal);
156 m_displayData->addPoint(geant4Message, ROOT::Math::XYZVector(global[0], global[1], global[2]) * Unit::mm);
157 }
158 }
159 m_nav.NewLevel(volume);
160}
161
162bool OverlapCheckerModule::checkVolume(G4VPhysicalVolume* volume, const std::string& path, int depth)
163{
164 // check if we exceeded maximum recursion depth
165 if (m_maxDepth > 0 && depth >= m_maxDepth) return false;
166 // add the volume to the navigation history
167 m_nav.NewLevel(volume);
168 // remember the path to the volume
169 std::string volumePath = path + "/" + volume->GetName();
170 bool result{false};
171 // check if we have a prefix we have to match. if not prefix or prefix matches, check for overlaps
172 if (m_prefix.empty() || volumePath.substr(0, m_prefix.size()) == m_prefix) {
173 result = volume->CheckOverlaps(m_points, m_tolerance, true, m_maxErrors);
174 if (result) {
175 m_overlaps.push_back(volumePath);
176 }
177 }
178
179 //Check the daughter volumes for overlaps
180 G4LogicalVolume* logicalVolume = volume->GetLogicalVolume();
181 for (size_t iDaughter = 0; iDaughter < logicalVolume->GetNoDaughters(); iDaughter++) {
182 G4VPhysicalVolume* daughter = logicalVolume->GetDaughter(iDaughter);
183 // check if we already checked this particular volume, if so skip it
184 auto it = m_seen.insert(daughter);
185 if (!it.second) continue;
186 result |= checkVolume(daughter, volumePath, depth + 1);
187 }
188 m_nav.BackLevel();
189 return result;
190}
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
OverlapCheckerModule()
Constructor of the module.
double m_tolerance
tolerance of overlap check
void initialize() override
Initialize the module.
std::string m_prefix
check only volumes beginning with prefix
void event() override
event function: this runs the overlap checker for each event
bool checkVolume(G4VPhysicalVolume *volume, const std::string &path, int depth=0)
Check a volume for overlaps.
int m_points
number of test points
int m_maxErrors
maximum number of errors before skipping current volume
void handleOverlap(const std::string &geant4Message)
Handle a G4Exception with the overlap message issued by Geant4.
std::vector< std::string > m_overlaps
list of overlapping volumes
int m_maxDepth
maximum depth to check
std::set< G4VPhysicalVolume * > m_seen
set of logical volumes we already checked
G4NavigationHistory m_nav
navigation history to remember coordinate transformations
StoreObjPtr< DisplayData > m_displayData
Pointer to the DisplayData where we add the overlap points for rendering.
static const double mm
[millimeters]
Definition: Unit.h:70
G4VPhysicalVolume * getTopVolume()
Return a pointer to the top volume.
static GeometryManager & getInstance()
Return a reference to the instance.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.