Belle II Software  release-05-02-19
GeometryManager.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Ritter *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <framework/logging/Logger.h>
12 #include <framework/gearbox/GearDir.h>
13 #include <framework/database/IntervalOfValidity.h>
14 #include <geometry/GeometryManager.h>
15 #include <geometry/Materials.h>
16 #include <geometry/CreatorManager.h>
17 #include <geometry/CreatorBase.h>
18 #include <geometry/utilities.h>
19 #include <geometry/dbobjects/GeoConfiguration.h>
20 #include <geometry/bfieldmap/BFieldMap.h>
21 #include <geometry/bfieldmap/BFieldFrameworkInterface.h>
22 #include <framework/dbobjects/MagneticField.h>
23 #include <framework/database/DBStore.h>
24 
25 #include "G4Box.hh"
26 #include "G4ThreeVector.hh"
27 #include "G4LogicalVolume.hh"
28 #include "G4PVPlacement.hh"
29 
30 #include "G4RunManager.hh"
31 #include "G4GeometryManager.hh"
32 #include "G4PhysicalVolumeStore.hh"
33 #include "G4LogicalVolumeStore.hh"
34 #include "G4SolidStore.hh"
35 #include "G4RegionStore.hh"
36 #include "G4SurfaceProperty.hh"
37 #include "G4LogicalBorderSurface.hh"
38 #include "G4LogicalSkinSurface.hh"
39 #include "G4VisAttributes.hh"
40 #include "G4VoxelLimits.hh"
41 
42 //VGM stuff
43 #include "Geant4GM/volumes/Factory.h"
44 #include "RootGM/volumes/Factory.h"
45 #include "VGM/volumes/IPlacement.h"
46 #include "TGeoManager.h"
47 
48 #include <memory>
49 
50 using namespace std;
51 
52 namespace Belle2 {
57  using namespace gearbox;
58 
59  namespace {
69  double getTopMinSize(const std::string& name, EAxis axis, G4LogicalVolume* volume, double current)
70  {
71  G4VoxelLimits dummyLimits;
72  double extent{0};
73  for (int i = 0; i < volume->GetNoDaughters(); ++i) {
74  G4VPhysicalVolume* daughter = volume->GetDaughter(i);
75  G4AffineTransform trans(daughter->GetRotation(), daughter->GetTranslation());
76  double vmin{0}, vmax{0};
77  daughter->GetLogicalVolume()->GetSolid()->CalculateExtent((EAxis) axis, dummyLimits, trans, vmin, vmax);
78  extent = std::max({extent, std::fabs(vmin), std::fabs(vmax)});
79  }
80  B2DEBUG(100, "Current global volume size in " << name << ": " << std::fixed
81  << std::setprecision(std::numeric_limits<double>::max_digits10)
82  << current << ", needed: " << extent);
83  if (current < extent && (extent - current) > Unit::um) {
84  if (current > 0) {
85  B2WARNING("Global volume not large enough in " << name << " direction, enlarging from "
86  << current << " mm to " << extent << " mm");
87  } else {
88  B2DEBUG(10, "Setting global volume size to " << name << "= +-" << extent << " mm");
89  }
90  return extent;
91  }
92  return current;
93  }
94  }
95 
96  namespace geometry {
97 
98  GeometryManager& GeometryManager::getInstance()
99  {
100  static unique_ptr<GeometryManager> instance(new GeometryManager());
101  return *instance;
102  }
103 
104  void GeometryManager::clear()
105  {
106  B2DEBUG(50, "Cleaning up Geometry");
107  for (G4VisAttributes* visAttr : m_VisAttributes) delete visAttr;
108  m_VisAttributes.clear();
109  for (CreatorBase* creator : m_creators) delete creator;
110  m_creators.clear();
111  m_topVolume = nullptr;
112  //Clean up existing Geometry
113  G4GeometryManager::GetInstance()->OpenGeometry();
114  G4PhysicalVolumeStore::Clean();
115  G4LogicalVolumeStore::Clean();
116  G4SolidStore::Clean();
117  G4LogicalBorderSurface::CleanSurfaceTable();
118  G4LogicalSkinSurface::CleanSurfaceTable();
119  G4SurfaceProperty::CleanSurfacePropertyTable();
120  //And finally clean up materials
121  Materials::getInstance().clear();
122  }
123 
124  GeoConfiguration GeometryManager::createGeometryConfig(const GearDir& detectorDir, const IntervalOfValidity& iov)
125  {
126  std::string detectorName;
127  try {
128  detectorName = detectorDir.getString("Name");
129  B2DEBUG(50, "Creating geometry for detector: " << detectorName);
130  } catch (gearbox::PathEmptyError& e) {
131  B2FATAL("Could not read detector name, make sure gearbox is connected and "
132  << detectorDir.getPath() << " points to the geometry description");
133  }
134 
135  const double globalWidth = detectorDir.getLength("Global/width", 0);
136  const double globalHeight = detectorDir.getLength("Global/height", 0);
137  const double globalLength = detectorDir.getLength("Global/length", 0);
138  const std::string globalMaterial = detectorDir.getString("Global/material", "Air");
139 
140  GeoConfiguration config(detectorName, globalWidth, globalHeight, globalLength, globalMaterial);
141 
142  // Add materials
143  Materials& materials = Materials::getInstance();
144  for (const GearDir& matlist : detectorDir.getNodes("Materials")) {
145  for (const GearDir& mat : matlist.getNodes("Material")) {
146  GeoMaterial material = materials.createMaterialConfig(mat);
147  config.addMaterial(material);
148  }
149  }
150 
151  std::set<std::string> componentNames = m_components;
152  std::set<std::string> excludedNames = m_excluded;
153  std::set<std::string> additionalNames = m_additional;
154  if (!m_components.empty() && !m_additional.empty()) {
155  B2WARNING("Additional components are ignored when a list of components is provided.");
156  }
157 
158  //Now create all subcomponents
159  for (const GearDir& component : detectorDir.getNodes("DetectorComponent")) {
160  string name;
161  string creatorName;
162  try {
163  name = component.getString("@name");
164  creatorName = component.getString("Creator");
165  } catch (gearbox::PathEmptyError& e) {
166  B2ERROR("Could not find required element Name or Creator for " << component.getPath());
167  continue;
168  }
169  const bool isDefault = component.getBool("@isDefault", true);
170  //Remove this name from the list of selected, excluded or additional components. At
171  //the end there should be nothing left in these lists
172  componentNames.erase(name);
173  excludedNames.erase(name);
174  additionalNames.erase(name);
175  if (!m_components.empty() && m_components.count(name) == 0) {
176  B2DEBUG(50, "DetectorComponent " << name << " not in list of components, skipping");
177  continue;
178  }
179  if (m_components.empty() && !isDefault) {
180  if (m_additional.count(name) == 0) {
181  B2DEBUG(50, "DectorComponent " << name << " is not enabled by default, skipping");
182  continue;
183  } else {
184  B2DEBUG(50, "DectorComponent " << name << " is enabled in addition to the default components");
185  }
186  }
187  if (!m_excluded.empty() && m_excluded.count(name) > 0) {
188  B2DEBUG(10, "DetectorComponent " << name << " in list of excluded components, skipping");
189  continue;
190  }
191 
192  string libraryName = component.getString("Creator/@library", "");
193  if (!iov.empty()) {
194  CreatorBase* creator = CreatorManager::getCreator(creatorName, libraryName);
195  if (creator) {
196  creator->createPayloads(GearDir(component, "Content"), iov);
197  } else {
198  B2ERROR("Could not load creator " << creatorName << " from " << libraryName);
199  }
200  }
201  config.addComponent({name, creatorName, libraryName});
202  }
203 
204  //If there are still names left in the componentNames, excludedNames or
205  //additionalNames there is probably an typo in the respective component
206  //list. Throw an error for each name left using a small lambda function
207  auto checkRemaining = [](const std::string & type, const std::set<std::string>& remainingNames) {
208  for (const std::string& name : remainingNames) {
209  B2ERROR("Geometry '" << name << "' is specified in list of "
210  << type << " but could not be found");
211  }
212  };
213 
214  checkRemaining("components", componentNames);
215  checkRemaining("excluded components", excludedNames);
216  checkRemaining("additional components", additionalNames);
217 
218  return config;
219  }
220 
221  void GeometryManager::createGeometry(const GearDir& detectorDir, GeometryTypes type)
222  {
223  GeoConfiguration config = createGeometryConfig(detectorDir, IntervalOfValidity());
224  createGeometry(config, type, false);
225  }
226 
227  void GeometryManager::createGeometry(const GeoConfiguration& config, GeometryTypes type, bool useDB)
228  {
229  // remove the old geometry
230  clear();
231 
232  // if we don't use the DB make sure the magnetic field is properly set
233  // by adding it as a fake database payload
234  BFieldMap::Instance().clear();
235  if (!useDB) {
236  auto* fieldmap = new MagneticField();
237  fieldmap->addComponent(new BFieldFrameworkInterface());
238  DBStore::Instance().addConstantOverride("MagneticField", fieldmap, false);
239  }
240 
241  //Let Geant4 know that we "modified" the geometry
242  G4RunManager* runManager = G4RunManager::GetRunManager();
243  if (runManager) runManager->ReinitializeGeometry(true, true);
244 
245  // create new geometry
246  B2DEBUG(10, "Creating geometry for detector: " << config.getName());
247 
248  // create Materials
249  Materials& materials = Materials::getInstance();
250  for (const GeoMaterial& mat : config.getMaterials()) {
251  materials.createMaterial(mat);
252  }
253 
254  //Now set Top volume. Be aware that Geant4 uses "half size" so the size
255  //will be in each direction even though the member name suggests a total size
256  G4Material* top_mat = Materials::get(config.getGlobalMaterial());
257  double xHalfLength = config.getGlobalWidth() / Unit::cm * CLHEP::cm;
258  double yHalfLength = config.getGlobalHeight() / Unit::cm * CLHEP::cm;
259  double zHalfLength = config.getGlobalLength() / Unit::cm * CLHEP::cm;
260  //Geant4 doesn't like negative dimensions so lets make sure it gets at
261  //least a positive value. We'll rezize the box anyway if the length is
262  //zero so this is not important now
263  G4Box* top_box = new G4Box("Top", xHalfLength > 0 ? xHalfLength : 1,
264  yHalfLength > 0 ? yHalfLength : 1,
265  zHalfLength > 0 ? zHalfLength : 1);
266  G4LogicalVolume* top_log = new G4LogicalVolume(top_box, top_mat, "Top", nullptr, nullptr, nullptr);
267  setVisibility(*top_log, false);
268  m_topVolume = new G4PVPlacement(nullptr, G4ThreeVector(), top_log, "Top", nullptr, false, 0);
269  B2DEBUG(10, "Created top volume with x= +-" << config.getGlobalWidth() << " cm, y= +-"
270  << config.getGlobalHeight() << " cm, z= +-" << config.getGlobalLength() << " cm");
271 
272 
273  auto getDensityScale = [this](const std::string & name) {
274  std::optional<double> scale;
275  // cppcheck-suppress stlIfFind ; cppcheck doesn't like if with initializer ...
276  if (auto it = m_densityScaling.find(name); it != m_densityScaling.end()) {
277  scale = it->second;
278  }
279  return scale;
280  };
281  auto globalScale = getDensityScale("*");
282 
283  for (const GeoComponent& component : config.getComponents()) {
284  CreatorBase* creator = CreatorManager::getCreator(component.getCreator(), component.getLibrary());
285  if (creator) {
286  // Do we want to scale the density for this or all components?
287  auto componentScale = getDensityScale(component.getName());
288  if (componentScale or globalScale) {
289  double scale = globalScale ? *globalScale : 1.0;
290  if (componentScale) scale *= *componentScale;
291  Materials::getInstance().setDensityScale(scale);
292  }
293  // remember how many we had to print the difference
294  int oldSolids = G4SolidStore::GetInstance()->size();
295  int oldLogical = G4LogicalVolumeStore::GetInstance()->size();
296  int oldPhysical = G4PhysicalVolumeStore::GetInstance()->size();
297  try {
298  if (!useDB) throw CreatorBase::DBNotImplemented();
299  creator->createFromDB(component.getName(), *top_log, type);
300  B2DEBUG(50, "called creator " << component.getCreator() << " to create component " << component.getName() << " from DB");
301 
302  } catch (CreatorBase::DBNotImplemented& e) {
303  GearDir parameters = Gearbox::getInstance().getDetectorComponent(component.getName());
304  creator->create(parameters, *top_log, type);
305  B2DEBUG(50, "called creator " << component.getCreator() << " to create component " << component.getName() << " from Gearbox");
306  }
307  int newSolids = G4SolidStore::GetInstance()->size() - oldSolids;
308  int newLogical = G4LogicalVolumeStore::GetInstance()->size() - oldLogical;
309  int newPhysical = G4PhysicalVolumeStore::GetInstance()->size() - oldPhysical;
310  B2DEBUG(50, "DetectorComponent " << component.getName() << " created " << newSolids
311  << " solids, " << newLogical << " logical volumes and "
312  << newPhysical << " physical volumes");
313  m_creators.push_back(creator);
314  //Done creating things, please reset density scaling
315  Materials::getInstance().resetDensityScale();
316  if (m_assignRegions) {
317  //Automatically assign a region with the creator name to all volumes
318  G4Region* region {nullptr};
319  //We loop over all children of the top volume and check whether they
320  //have the correct region assigned
321  for (int i = 0; i < top_log->GetNoDaughters(); ++i) {
322  G4LogicalVolume* vol = top_log->GetDaughter(i)->GetLogicalVolume();
323  //We only assign a region if there is not already one
324  if (!vol->GetRegion()) {
325  //Ok, new region, create or get one if not already done and assign it
326  if (!region) region = G4RegionStore::GetInstance()->FindOrCreateRegion(component.getName());
327  vol->SetRegion(region);
328  //And propagate the region to all child volumes
329  region->AddRootLogicalVolume(vol);
330  }
331  }
332  }
333  }
334  }
335 
336  int newSolids = G4SolidStore::GetInstance()->size();
337  int newLogical = G4LogicalVolumeStore::GetInstance()->size();
338  int newPhysical = G4PhysicalVolumeStore::GetInstance()->size();
339  B2DEBUG(50, "Created a total of " << newSolids << " solids, " << newLogical
340  << " logical volumes and " << newPhysical << " physical volumes");
341 
342  //Enlarge top volume to always encompass all volumes
343  top_box->SetXHalfLength(getTopMinSize("x", kXAxis, top_log, xHalfLength));
344  top_box->SetYHalfLength(getTopMinSize("y", kYAxis, top_log, yHalfLength));
345  top_box->SetZHalfLength(getTopMinSize("z", kZAxis, top_log, zHalfLength));
346 
347  B2DEBUG(50, "Optimizing geometry and creating lookup tables ...");
348  G4GeometryManager::GetInstance()->CloseGeometry(true, LogSystem::Instance().isLevelEnabled(LogConfig::c_Debug, 200, PACKAGENAME()));
349  }
350 
351  void GeometryManager::createTGeoRepresentation()
352  {
353  if (!m_topVolume) {
354  B2ERROR("No Geometry found, please create a geometry before converting it to ROOT::TGeo");
355  return;
356  }
357 
358  Geant4GM::Factory g4Factory;
359  if (LogSystem::Instance().isLevelEnabled(LogConfig::c_Debug, 200, PACKAGENAME())) {
360  g4Factory.SetDebug(1);
361  }
362  g4Factory.Import(m_topVolume);
363  RootGM::Factory rtFactory;
364  rtFactory.SetIgnore(1);
365  if (LogSystem::Instance().isLevelEnabled(LogConfig::c_Debug, 200, PACKAGENAME())) {
366  rtFactory.SetDebug(1);
367  }
368  g4Factory.Export(&rtFactory);
369  gGeoManager->CloseGeometry();
370  delete g4Factory.Top();
371  delete rtFactory.Top();
372  }
373 
374  G4VisAttributes* GeometryManager::newVisAttributes()
375  {
376  m_VisAttributes.push_back(new G4VisAttributes());
377  return m_VisAttributes.back();
378  }
379  }
381 } //Belle2 namespace
Belle2::IntervalOfValidity
A class that describes the interval of experiments/runs for which an object in the database is valid.
Definition: IntervalOfValidity.h:35
Belle2::GeoMaterial
Class to represent a material informaion in the Database.
Definition: GeoMaterial.h:32
Belle2::GeoConfiguration
configuration of the geometry
Definition: GeoConfiguration.h:33
Belle2::GeoComponent
Describe one component of the Geometry.
Definition: GeoComponent.h:29
Belle2::geometry::CreatorBase::create
virtual void create(const GearDir &content, G4LogicalVolume &topVolume, GeometryTypes type)=0
Function to actually create the geometry, has to be overridden by derived classes.
Belle2::geometry::Materials::createMaterial
G4Material * createMaterial(const gearbox::Interface &parameters)
Create a material from the parameters specified by parameters.
Definition: Materials.cc:159
Belle2::gearbox::Interface::getPath
std::string getPath() const
Return path of the current interface.
Definition: Interface.h:72
Belle2::geometry::Materials::createMaterialConfig
GeoMaterial createMaterialConfig(const gearbox::Interface &parameters)
Create Material from XML description.
Definition: Materials.cc:226
Belle2::geometry::GeometryManager
Class to manage the creation and conversion of the geometry.
Definition: GeometryManager.h:50
Belle2::geometry::CreatorBase::createFromDB
virtual void createFromDB(const std::string &name, G4LogicalVolume &topVolume, GeometryTypes type)
Function to create the geometry from the Database.
Definition: CreatorBase.cc:27
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::geometry::setVisibility
void setVisibility(G4LogicalVolume &volume, bool visible)
Helper function to quickly set the visibility of a given volume.
Definition: utilities.cc:115
Belle2::BFieldFrameworkInterface
Simple BFieldComponent to just wrap the existing BFieldMap with the new BFieldManager.
Definition: BFieldFrameworkInterface.h:25
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
Belle2::geometry::CreatorBase::createPayloads
virtual void createPayloads(const GearDir &content, const IntervalOfValidity &iov)
Function to create the geometry database.
Definition: CreatorBase.cc:34
Belle2::geometry::Materials
Thin wrapper around the Geant4 Material system.
Definition: Materials.h:50
Belle2::MagneticField
Magnetic field map.
Definition: MagneticField.h:43
Belle2::gearbox::Interface::getLength
double getLength(const std::string &path="") const noexcept(false)
Get the parameter path as a double converted to the standard length unit.
Definition: Interface.h:261
Belle2::GearDir::getString
virtual std::string getString(const std::string &path="") const noexcept(false) override
Get the parameter path as a string.
Definition: GearDir.h:79
Belle2::geometry::CreatorBase
Pure virtual base class for all geometry creators.
Definition: CreatorBase.h:31
Belle2::gearbox::Interface::getNodes
std::vector< GearDir > getNodes(const std::string &path="") const
Get vector of GearDirs which point to all the nodes the given path evaluates to.
Definition: Interface.cc:31
Belle2::geometry::GeometryTypes
GeometryTypes
Flag indiciating the type of geometry to be used.
Definition: GeometryManager.h:39