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>
20 #include <boost/algorithm/string.hpp>
22 #include <G4Material.hh>
23 #include <G4Element.hh>
24 #include <G4NistManager.hh>
25 #include <G4OpticalSurface.hh>
26 #include <G4MaterialPropertiesTable.hh>
28 #include "CLHEP/Units/PhysicalConstants.h"
43 template<
class T>
void addProperties(T&
object,
const gearbox::Interface& parameters)
46 if (parameters.getNumberNodes(
"Property") > 0) {
48 for (
const GearDir& property : parameters.getNodes(
"Property")) {
51 name =
property.getString(
"@name");
52 }
catch (gearbox::PathEmptyError&) {
53 B2ERROR(
"Property at " << property.getPath() <<
" does not have a name attribute, ignoring");
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();
65 object.addProperty({name, energies, values});
71 G4MaterialPropertiesTable* Materials::createProperties(
const std::vector<GeoMaterialProperty>& props)
73 if (!props.size())
return nullptr;
74 auto* table =
new G4MaterialPropertiesTable();
75 m_PropTables.push_back(table);
77 std::vector<double> energies = prop.getEnergies();
78 std::vector<double> values = prop.getValues();
79 const size_t n = values.size();
81 for (
double& energy : energies) energy *= CLHEP::GeV / Unit::GeV;
82 table->AddProperty(prop.getName().c_str(), energies.data(), values.data(), n);
89 static unique_ptr<Materials> instance(
new Materials());
93 Materials::~Materials()
98 void Materials::initBuilders()
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);
107 G4Material* Materials::findMaterial(
const string& name)
109 G4Material* mat{
nullptr};
110 if (m_materialCache.retrieve(name, mat)) {
114 mat = m_nistMaterialBuilder->FindOrBuildMaterial(name);
117 if (!mat && name.substr(0, 3) !=
"G4_") {
119 mat = m_nistMaterialBuilder->FindOrBuildMaterial(
"G4_" + name);
121 if (!mat) mat = m_nistMaterialBuilder->FindOrBuildMaterial(
"G4_" + boost::to_upper_copy(name));
125 m_materialCache.insert(name, mat);
129 G4Material* Materials::getMaterial(
const string& name,
bool showErrors)
131 G4Material* mat = findMaterial(name);
132 if (!mat && showErrors) B2ERROR(
"Material '" << name <<
"' could not be found");
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);
142 scaledMat =
new G4Material(scaled, mat->GetDensity() * *m_densityScale, mat,
143 mat->GetState(), mat->GetTemperature(), mat->GetPressure());
144 m_materialCache.insert(scaled, scaledMat);
151 G4Element* Materials::getElement(
const string& name)
154 G4Element* elm = m_nistElementBuilder->FindOrBuildElement(name);
155 if (!elm) B2ERROR(
"Element '" << name <<
"' could not be found");
161 GeoMaterial material = createMaterialConfig(parameters);
162 return createMaterial(material);
165 G4Material* Materials::createMaterial(
const GeoMaterial& parameters)
167 G4Material* oldmat = getMaterial(parameters.getName(),
false);
170 B2ERROR(
"Material with name " << parameters.getName() <<
" already existing");
173 B2DEBUG(10,
"Creating Material " << parameters.getName());
176 double density = parameters.getDensity() * CLHEP::g / CLHEP::cm3;
177 if (density < 1e-25) {
182 if (component.getIselement()) {
183 B2ERROR(
"createMaterial " << parameters.getName()
184 <<
": Cannot calculate density when adding elements, please provde a density");
187 G4Material* mat = getMaterial(component.getName());
189 B2ERROR(
"createMaterial " << parameters.getName() <<
": Material '" << component.getName() <<
"' not found");
192 density += mat->GetDensity() * component.getFraction();
197 G4Material* mat =
new G4Material(parameters.getName(), density, parameters.getComponents().size(),
198 (G4State)parameters.getState(), parameters.getTemperature(), parameters.getPressure() * CLHEP::pascal);
201 if (component.getIselement()) {
202 G4Element* cmp = getElement(component.getName());
204 B2ERROR(
"Cannot create material " << parameters.getName() <<
": element " << component.getName() <<
" not found");
207 mat->AddElement(cmp, component.getFraction());
209 G4Material* cmp = getMaterial(component.getName());
211 B2ERROR(
"Cannot create material " << parameters.getName() <<
": material " << component.getName() <<
" not found");
214 mat->AddMaterial(cmp, component.getFraction());
219 mat->SetMaterialPropertiesTable(createProperties(parameters.getProperties()));
222 m_materialCache.insert(parameters.getName(), mat);
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);
233 double temperature = parameters.getDouble(
"temperature", CLHEP::STP_Temperature);
234 double pressure = parameters.getDouble(
"pressure", CLHEP::STP_Pressure / CLHEP::pascal);
244 double sumFractions(0);
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;
254 for (
const GearDir& element : parameters.getNodes(
"Components/Element")) {
255 const std::string elementName = element.getString();
256 double fraction = element.getDouble(
"@fraction", 1.0);
258 sumFractions += fraction;
262 if (abs(sumFractions - 1) > numeric_limits<double>::epsilon()) {
263 B2WARNING(
"createMaterial " << name <<
": Fractions not normalized, scaling by 1/" << sumFractions);
265 cmp.setFraction(cmp.getFraction() / sumFractions);
270 boost::to_lower(stateStr);
271 G4State state = kStateUndefined;
272 if (stateStr ==
"solid") {
274 }
else if (stateStr ==
"liquid") {
275 state = kStateLiquid;
276 }
else if (stateStr ==
"gas") {
278 }
else if (stateStr !=
"undefined") {
279 B2WARNING(
"createMaterial " << name <<
": Unknown state '" << stateStr <<
"', using undefined");
283 addProperties(material, parameters);
291 return createOpticalSurface(surface);
296 G4OpticalSurface* optSurf =
new G4OpticalSurface(surface.
getName(),
297 (G4OpticalSurfaceModel) surface.
getModel(),
298 (G4OpticalSurfaceFinish) surface.
getFinish(),
299 (G4SurfaceType) surface.
getType(),
301 optSurf->SetMaterialPropertiesTable(createProperties(surface.
getProperties()));
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);
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)
321 B2FATAL(
"Unknown Optical Surface Model: " << modelString);
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)
357 B2FATAL(
"Unknown Optical Surface Finish: " << finishString);
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)
368 B2FATAL(
"Unknown Optical Surface Type: " << typeString);
370 #undef CHECK_ENUM_VALUE
373 addProperties(surface,
GearDir(parameters,
"Properties"));
377 void Materials::clear()
380 B2DEBUG(50,
"Cleaning G4MaterialPropertiesTable");
381 for (G4MaterialPropertiesTable* prop : m_PropTables)
delete prop;
382 m_PropTables.clear();
384 B2DEBUG(50,
"Cleaning G4Materials");
385 G4MaterialTable& materials = *G4Material::GetMaterialTable();
386 for (G4Material* mat : materials)
delete mat;
388 B2DEBUG(50,
"Cleaning G4Elements");
389 G4ElementTable& elements = *G4Element::GetElementTable();
390 for (G4Element* elm : elements)
delete elm;
392 B2DEBUG(50,
"Cleaning G4Isotopes");
393 auto& isotopes =
const_cast<G4IsotopeTable&
>(*G4Isotope::GetIsotopeTable());
394 for (G4Isotope* iso : isotopes)
delete iso;
397 delete m_nistMaterialBuilder;
398 delete m_nistElementBuilder;
399 m_nistMaterialBuilder =
nullptr;
400 m_nistElementBuilder =
nullptr;
402 B2DEBUG(50,
"Clean up material cache");
403 m_materialCache.clear();