Belle II Software  release-05-01-25
Materials.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 <geometry/Materials.h>
12 #include <framework/logging/Logger.h>
13 #include <framework/gearbox/GearDir.h>
14 #include <framework/gearbox/Unit.h>
15 #include <geometry/dbobjects/GeoMaterial.h>
16 #include <geometry/dbobjects/GeoOpticalSurface.h>
17 
18 #include <memory>
19 #include <limits>
20 #include <boost/algorithm/string.hpp>
21 
22 #include <G4Material.hh>
23 #include <G4Element.hh>
24 #include <G4NistManager.hh>
25 #include <G4OpticalSurface.hh>
26 #include <G4MaterialPropertiesTable.hh>
27 
28 #include "CLHEP/Units/PhysicalConstants.h"
29 
30 using namespace std;
31 
32 namespace Belle2 {
37  namespace geometry {
38  namespace {
43  template<class T> void addProperties(T& object, const gearbox::Interface& parameters)
44  {
45  // check if we have properties defined for this material
46  if (parameters.getNumberNodes("Property") > 0) {
47  //Apparantly we have properties, so lets add them to the Material
48  for (const GearDir& property : parameters.getNodes("Property")) {
49  string name;
50  try {
51  name = property.getString("@name");
52  } catch (gearbox::PathEmptyError&) {
53  B2ERROR("Property at " << property.getPath() << " does not have a name attribute, ignoring");
54  continue;
55  }
56  const double conversionFactor = Unit::convertValue(1, property.getString("@unit", "GeV"));
57  unsigned int nValues = property.getNumberNodes("value");
58  std::vector<double> energies(nValues, 0);
59  std::vector<double> values(nValues, 0);
60  for (unsigned int i = 0; i < nValues; ++i) {
61  GearDir value(property, "value", i + 1);
62  energies[i] = value.getDouble("@energy") * conversionFactor;
63  values[i] = value.getDouble();
64  }
65  object.addProperty({name, energies, values});
66  }
67  }
68  }
69  }
70 
71  G4MaterialPropertiesTable* Materials::createProperties(const std::vector<GeoMaterialProperty>& props)
72  {
73  if (!props.size()) return nullptr;
74  auto* table = new G4MaterialPropertiesTable();
75  m_PropTables.push_back(table);
76  for (const GeoMaterialProperty& prop : props) {
77  std::vector<double> energies = prop.getEnergies();
78  std::vector<double> values = prop.getValues();
79  const size_t n = values.size();
80  // convert energies to Geant4 Units
81  for (double& energy : energies) energy *= CLHEP::GeV / Unit::GeV;
82  table->AddProperty(prop.getName().c_str(), energies.data(), values.data(), n);
83  }
84  return table;
85  }
86 
87  Materials& Materials::getInstance()
88  {
89  static unique_ptr<Materials> instance(new Materials());
90  return *instance;
91  }
92 
93  Materials::~Materials()
94  {
95  clear();
96  }
97 
98  void Materials::initBuilders()
99  {
100  if (!m_nistMaterialBuilder) {
101  B2DEBUG(50, "Creating new Nist Builder instances");
102  m_nistElementBuilder = new G4NistElementBuilder(0);
103  m_nistMaterialBuilder = new G4NistMaterialBuilder(m_nistElementBuilder, 0);
104  }
105  }
106 
107  G4Material* Materials::findMaterial(const string& name)
108  {
109  G4Material* mat{nullptr};
110  if (m_materialCache.retrieve(name, mat)) {
111  return mat;
112  }
113  initBuilders();
114  mat = m_nistMaterialBuilder->FindOrBuildMaterial(name);
115 
116  //Try different combinations of the Material name to fallback to predefined G4 Elements
117  if (!mat && name.substr(0, 3) != "G4_") {
118  //Mainly for materials from single elements, e.g. G4_Al, G4_Si
119  mat = m_nistMaterialBuilder->FindOrBuildMaterial("G4_" + name);
120  //For predefined materials, e.g. G4_AIR, G4_TEFLON
121  if (!mat) mat = m_nistMaterialBuilder->FindOrBuildMaterial("G4_" + boost::to_upper_copy(name));
122  }
123 
124  //Insert into cache
125  m_materialCache.insert(name, mat);
126  return mat;
127  }
128 
129  G4Material* Materials::getMaterial(const string& name, bool showErrors)
130  {
131  G4Material* mat = findMaterial(name);
132  if (!mat && showErrors) B2ERROR("Material '" << name << "' could not be found");
133  // There is a chance that we want to scale densities of all materials for studies.
134  // Now we have the original so we can look to the scaled one or create it if it's missing
135  if (m_densityScale) {
136  if (m_ignoreScaling.count(name) > 0) return mat;
137  if (name.substr(0, 7) == "scaled:")
138  B2WARNING("Requested an already scaled material ... scaling once more" << LogVar("material", name));
139  std::string scaled = "scaled:" + std::to_string(*m_densityScale) + ":" + name;
140  G4Material* scaledMat = findMaterial(scaled);
141  if (!scaledMat) {
142  scaledMat = new G4Material(scaled, mat->GetDensity() * *m_densityScale, mat,
143  mat->GetState(), mat->GetTemperature(), mat->GetPressure());
144  m_materialCache.insert(scaled, scaledMat);
145  }
146  return scaledMat;
147  }
148  return mat;
149  }
150 
151  G4Element* Materials::getElement(const string& name)
152  {
153  initBuilders();
154  G4Element* elm = m_nistElementBuilder->FindOrBuildElement(name);
155  if (!elm) B2ERROR("Element '" << name << "' could not be found");
156  return elm;
157  }
158 
159  G4Material* Materials::createMaterial(const gearbox::Interface& parameters)
160  {
161  GeoMaterial material = createMaterialConfig(parameters);
162  return createMaterial(material);
163  }
164 
165  G4Material* Materials::createMaterial(const GeoMaterial& parameters)
166  {
167  G4Material* oldmat = getMaterial(parameters.getName(), false);
168  if (oldmat) {
169  //We cannot modify or delete materials, so we just use existing ones
170  B2ERROR("Material with name " << parameters.getName() << " already existing");
171  return oldmat;
172  }
173  B2DEBUG(10, "Creating Material " << parameters.getName());
174  //If density is negative or smaller than epsilon we should calculate the
175  //density from the used materials
176  double density = parameters.getDensity() * CLHEP::g / CLHEP::cm3;
177  if (density < 1e-25) {
178  density = 0;
179  for (const GeoMaterialComponent& component : parameters.getComponents()) {
180  //Check if we can calculate density, only works if we just add materials,
181  //elements do not have density
182  if (component.getIselement()) {
183  B2ERROR("createMaterial " << parameters.getName()
184  << ": Cannot calculate density when adding elements, please provde a density");
185  return nullptr;
186  }
187  G4Material* mat = getMaterial(component.getName());
188  if (!mat) {
189  B2ERROR("createMaterial " << parameters.getName() << ": Material '" << component.getName() << "' not found");
190  return nullptr;
191  }
192  density += mat->GetDensity() * component.getFraction();
193  }
194  }
195 
196  //Finally, create Material and add all components
197  G4Material* mat = new G4Material(parameters.getName(), density, parameters.getComponents().size(),
198  (G4State)parameters.getState(), parameters.getTemperature(), parameters.getPressure() * CLHEP::pascal);
199 
200  for (const GeoMaterialComponent& component : parameters.getComponents()) {
201  if (component.getIselement()) {
202  G4Element* cmp = getElement(component.getName());
203  if (!cmp) {
204  B2ERROR("Cannot create material " << parameters.getName() << ": element " << component.getName() << " not found");
205  return nullptr;
206  }
207  mat->AddElement(cmp, component.getFraction());
208  } else {
209  G4Material* cmp = getMaterial(component.getName());
210  if (!cmp) {
211  B2ERROR("Cannot create material " << parameters.getName() << ": material " << component.getName() << " not found");
212  return nullptr;
213  }
214  mat->AddMaterial(cmp, component.getFraction());
215  }
216  }
217 
218  //add properties
219  mat->SetMaterialPropertiesTable(createProperties(parameters.getProperties()));
220 
221  //Insert into cache and return
222  m_materialCache.insert(parameters.getName(), mat);
223  return mat;
224  }
225 
226  GeoMaterial Materials::createMaterialConfig(const gearbox::Interface& parameters)
227  {
228  //Get straightforward parameters from gearbox
229  string name = parameters.getString("@name");
230  B2DEBUG(10, "Creating Material Config " << name);
231  string stateStr = parameters.getString("state", "undefined");
232  double density = parameters.getDensity("density", 0); // * CLHEP::g / CLHEP::cm3;
233  double temperature = parameters.getDouble("temperature", CLHEP::STP_Temperature);
234  double pressure = parameters.getDouble("pressure", CLHEP::STP_Pressure / CLHEP::pascal);
235 
236  // set the simple ones like they are
237  GeoMaterial material;
238  material.setName(name);
239  material.setPressure(pressure);
240  material.setTemperature(temperature);
241  material.setDensity(density);
242 
243  //Sum up fractions to see if they are normalized
244  double sumFractions(0);
245 
246  // Add all component materials
247  for (const GearDir& component : parameters.getNodes("Components/Material")) {
248  const string componentName = component.getString();
249  double fraction = component.getDouble("@fraction", 1.0);
250  material.addComponent({componentName, false, fraction});
251  sumFractions += fraction;
252  }
253  // Add all component elements
254  for (const GearDir& element : parameters.getNodes("Components/Element")) {
255  const std::string elementName = element.getString();
256  double fraction = element.getDouble("@fraction", 1.0);
257  material.addComponent({elementName, true, fraction});
258  sumFractions += fraction;
259  }
260 
261  //Warn if the fractions are not normalized
262  if (abs(sumFractions - 1) > numeric_limits<double>::epsilon()) {
263  B2WARNING("createMaterial " << name << ": Fractions not normalized, scaling by 1/" << sumFractions);
264  for (GeoMaterialComponent& cmp : material.getComponents()) {
265  cmp.setFraction(cmp.getFraction() / sumFractions);
266  }
267  }
268 
269  //Determine material state
270  boost::to_lower(stateStr);
271  G4State state = kStateUndefined;
272  if (stateStr == "solid") {
273  state = kStateSolid;
274  } else if (stateStr == "liquid") {
275  state = kStateLiquid;
276  } else if (stateStr == "gas") {
277  state = kStateGas;
278  } else if (stateStr != "undefined") {
279  B2WARNING("createMaterial " << name << ": Unknown state '" << stateStr << "', using undefined");
280  }
281  material.setState(state);
282 
283  addProperties(material, parameters);
284  // all done, return material description
285  return material;
286  }
287 
288  G4OpticalSurface* Materials::createOpticalSurface(const gearbox::Interface& parameters)
289  {
290  GeoOpticalSurface surface = createOpticalSurfaceConfig(parameters);
291  return createOpticalSurface(surface);
292  }
293 
294  G4OpticalSurface* Materials::createOpticalSurface(const GeoOpticalSurface& surface)
295  {
296  G4OpticalSurface* optSurf = new G4OpticalSurface(surface.getName(),
297  (G4OpticalSurfaceModel) surface.getModel(),
298  (G4OpticalSurfaceFinish) surface.getFinish(),
299  (G4SurfaceType) surface.getType(),
300  surface.getValue());
301  optSurf->SetMaterialPropertiesTable(createProperties(surface.getProperties()));
302  //m_OptSurfaces.push_back(optSurf);
303  return optSurf;
304  }
305 
306  GeoOpticalSurface Materials::createOpticalSurfaceConfig(const gearbox::Interface& parameters)
307  {
308  string name = parameters.getString("@name", "OpticalSurface");
309  string modelString = parameters.getString("Model", "glisur");
310  string finishString = parameters.getString("Finish", "polished");
311  string typeString = parameters.getString("Type", "dielectric_dielectric");
312  double value = parameters.getDouble("Value", 1.0);
313 
314 #define CHECK_ENUM_VALUE(name,value) if (name##String == #value) { name = value; }
315  G4OpticalSurfaceModel model;
316  boost::to_lower(modelString);
317  CHECK_ENUM_VALUE(model, glisur)
318  else CHECK_ENUM_VALUE(model, unified)
319  else CHECK_ENUM_VALUE(model, LUT)
320  else {
321  B2FATAL("Unknown Optical Surface Model: " << modelString);
322  }
323 
324  G4OpticalSurfaceFinish finish;
325  boost::to_lower(finishString);
326  CHECK_ENUM_VALUE(finish, polished)
327  else CHECK_ENUM_VALUE(finish, polishedfrontpainted)
328  else CHECK_ENUM_VALUE(finish, polishedbackpainted)
329  else CHECK_ENUM_VALUE(finish, ground)
330  else CHECK_ENUM_VALUE(finish, groundfrontpainted)
331  else CHECK_ENUM_VALUE(finish, groundbackpainted)
332  else CHECK_ENUM_VALUE(finish, polishedlumirrorair)
333  else CHECK_ENUM_VALUE(finish, polishedlumirrorglue)
334  else CHECK_ENUM_VALUE(finish, polishedair)
335  else CHECK_ENUM_VALUE(finish, polishedteflonair)
336  else CHECK_ENUM_VALUE(finish, polishedtioair)
337  else CHECK_ENUM_VALUE(finish, polishedtyvekair)
338  else CHECK_ENUM_VALUE(finish, polishedvm2000air)
339  else CHECK_ENUM_VALUE(finish, polishedvm2000glue)
340  else CHECK_ENUM_VALUE(finish, etchedlumirrorair)
341  else CHECK_ENUM_VALUE(finish, etchedlumirrorglue)
342  else CHECK_ENUM_VALUE(finish, etchedair)
343  else CHECK_ENUM_VALUE(finish, etchedteflonair)
344  else CHECK_ENUM_VALUE(finish, etchedtioair)
345  else CHECK_ENUM_VALUE(finish, etchedtyvekair)
346  else CHECK_ENUM_VALUE(finish, etchedvm2000air)
347  else CHECK_ENUM_VALUE(finish, etchedvm2000glue)
348  else CHECK_ENUM_VALUE(finish, groundlumirrorair)
349  else CHECK_ENUM_VALUE(finish, groundlumirrorglue)
350  else CHECK_ENUM_VALUE(finish, groundair)
351  else CHECK_ENUM_VALUE(finish, groundteflonair)
352  else CHECK_ENUM_VALUE(finish, groundtioair)
353  else CHECK_ENUM_VALUE(finish, groundtyvekair)
354  else CHECK_ENUM_VALUE(finish, groundvm2000air)
355  else CHECK_ENUM_VALUE(finish, groundvm2000glue)
356  else {
357  B2FATAL("Unknown Optical Surface Finish: " << finishString);
358  }
359 
360  G4SurfaceType type;
361  boost::to_lower(typeString);
362  CHECK_ENUM_VALUE(type, dielectric_metal)
363  else CHECK_ENUM_VALUE(type, dielectric_dielectric)
364  else CHECK_ENUM_VALUE(type, dielectric_LUT)
365  else CHECK_ENUM_VALUE(type, firsov)
366  else CHECK_ENUM_VALUE(type, x_ray)
367  else {
368  B2FATAL("Unknown Optical Surface Type: " << typeString);
369  }
370 #undef CHECK_ENUM_VALUE
371 
372  GeoOpticalSurface surface(name, model, finish, type, value);
373  addProperties(surface, GearDir(parameters, "Properties"));
374  return surface;
375  }
376 
377  void Materials::clear()
378  {
379  // delete all property tables we have
380  B2DEBUG(50, "Cleaning G4MaterialPropertiesTable");
381  for (G4MaterialPropertiesTable* prop : m_PropTables) delete prop;
382  m_PropTables.clear();
383  // and last but not least: get rid of all materials
384  B2DEBUG(50, "Cleaning G4Materials");
385  G4MaterialTable& materials = *G4Material::GetMaterialTable();
386  for (G4Material* mat : materials) delete mat;
387  materials.clear();
388  B2DEBUG(50, "Cleaning G4Elements");
389  G4ElementTable& elements = *G4Element::GetElementTable();
390  for (G4Element* elm : elements) delete elm;
391  elements.clear();
392  B2DEBUG(50, "Cleaning G4Isotopes");
393  auto& isotopes = const_cast<G4IsotopeTable&>(*G4Isotope::GetIsotopeTable());
394  for (G4Isotope* iso : isotopes) delete iso;
395  isotopes.clear();
396  // delete material and element builder as they keep indices to materials they created :/
397  delete m_nistMaterialBuilder;
398  delete m_nistElementBuilder;
399  m_nistMaterialBuilder = nullptr;
400  m_nistElementBuilder = nullptr;
401  // finally, get rid of the cache, it's invalid now anyway
402  B2DEBUG(50, "Clean up material cache");
403  m_materialCache.clear();
404  }
405 
406  } //geometry namespace
408 } //Belle2 namespace
Belle2::GeoMaterial::addComponent
void addComponent(const GeoMaterialComponent &component)
add a component to the material.
Definition: GeoMaterial.h:47
Belle2::GeoOpticalSurface::getProperties
const std::vector< GeoMaterialProperty > & getProperties() const
get all properties
Definition: GeoOpticalSurface.h:61
Belle2::GeoMaterial::setName
void setName(const std::string &name)
set the name of the material
Definition: GeoMaterial.h:37
Belle2::GeoMaterial
Class to represent a material informaion in the Database.
Definition: GeoMaterial.h:32
Belle2::GeoMaterial::setTemperature
void setTemperature(double temperature)
set the temperature of the material (in default framework units
Definition: GeoMaterial.h:43
Belle2::GeoOpticalSurface::getName
const std::string & getName() const
get name of the optical surface
Definition: GeoOpticalSurface.h:51
Belle2::GeoMaterialComponent
Component of a material.
Definition: GeoMaterialComponent.h:29
Belle2::GeoMaterialProperty
Property of a material.
Definition: GeoMaterialProperty.h:29
Belle2::GeoOpticalSurface
Represent an optical finish of a surface.
Definition: GeoOpticalSurface.h:31
Belle2::GeoMaterial::setState
void setState(int state)
set the state of the material
Definition: GeoMaterial.h:39
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::GeoMaterial::setPressure
void setPressure(double pressure)
set the pressure of the material (in default framework units
Definition: GeoMaterial.h:45
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
LogVar
Class to store variables with their name which were sent to the logging service.
Definition: LogVariableStream.h:24
Belle2::geometry::Materials
Thin wrapper around the Geant4 Material system.
Definition: Materials.h:50
Belle2::GeoOpticalSurface::getValue
double getValue() const
get value for the surface condition
Definition: GeoOpticalSurface.h:59
Belle2::gearbox::Interface
Exception to be thrown in case of an empty result.
Definition: Interface.h:37
Belle2::GeoMaterial::getComponents
std::vector< GeoMaterialComponent > & getComponents()
get all components
Definition: GeoMaterial.h:61
Belle2::GeoOpticalSurface::getModel
int getModel() const
get model for the surface
Definition: GeoOpticalSurface.h:53
Belle2::GeoOpticalSurface::getFinish
int getFinish() const
get finish of the surface
Definition: GeoOpticalSurface.h:55
Belle2::GeoOpticalSurface::getType
int getType() const
get type of the surface
Definition: GeoOpticalSurface.h:57
Belle2::GeoMaterial::setDensity
void setDensity(double density)
set the density of the material (in default framework units
Definition: GeoMaterial.h:41