8 #include "GeoModelKernel/GeoIntrusivePtr.h"
9 #include "GeoModelKernel/GeoMaterial.h"
10 #include "GeoModelKernel/GeoElement.h"
11 #include "GeoModelKernel/Units.h"
12 #include "GaudiKernel/SystemOfUnits.h"
27 m_managerName(managerName),
28 m_extraFunctionality(false),
29 m_athenaComps(nullptr) {
37 const std::string& space,
38 bool extraFunctionality)
40 m_managerName(managerName),
41 m_extraFunctionality(extraFunctionality),
42 m_athenaComps(nullptr) {
57 const std::string& space)
59 m_managerName(managerName),
60 m_extraFunctionality(true),
61 m_athenaComps(nullptr) {
71 m_managerName(managerName),
72 m_extraFunctionality(true),
73 m_athenaComps(athenaComps) {
88 std::string errorMessage(
"Null pointer to Stored Material Manager!");
90 throw std::runtime_error(errorMessage);
113 std::string errorMessage(
"Null pointer to Stored Material Manager!");
115 throw std::runtime_error(errorMessage);
125 MaterialStore::const_iterator iter;
135 const double volumeTot,
136 const double volume1,
const std::string& matName1,
137 const double volume2,
const std::string& matName2
139 std::vector<std::string> baseMaterials;
140 std::vector<double> fracWeight;
141 baseMaterials.reserve(2);
142 fracWeight.reserve(2);
146 ATH_MSG_DEBUG(
"Composite material : " << matName1 <<
" " << matName2);
148 double density1, density2;
150 MaterialWeightMap::const_iterator iter;
153 density1 = mat1->getDensity();
157 density1 = mat1->getDensity();
163 density2 = mat2->getDensity();
167 density2 = mat2->getDensity();
171 double weight1 = density1 * volume1;
172 double weight2 = density2 * volume2;
173 double invWeightTot = 1.0 / (weight1 + weight2);
175 double density = (weight1 + weight2) / volumeTot;
177 double frac1 = weight1 / (weight1 + weight2);
178 double frac2 = weight2 / (weight1 + weight2);
179 double density_2 = 1.0 / (frac1 / density1 + frac2 / density2);
180 double density_3 = (weight1 + weight2) / (volume1 + volume2);
186 baseMaterials.push_back(matName1);
187 baseMaterials.push_back(matName2);
188 fracWeight.push_back(weight1 * invWeightTot);
189 fracWeight.push_back(weight2 * invWeightTot);
191 return getMaterial(newMatName, baseMaterials, fracWeight, density);
217 std::string newName2 =
newName;
218 bool newNameProvided = !newName2.empty();
219 if (!newNameProvided) {
220 newName2 = origMaterialName +
"Modified";
223 const GeoMaterial* newMaterial =
nullptr;
233 newMaterial = material;
236 newMaterial = origMaterial;
240 if (newNameProvided || !
compareDensity(origMaterial->getDensity(), density)) {
242 GeoMaterial* newMaterialTmp =
new GeoMaterial(newName2, density);
243 newMaterialTmp->add(
const_cast<GeoMaterial*
>(origMaterial), 1.);
245 newMaterial = newMaterialTmp;
265 if (scaleFactor > 1000 || scaleFactor < 0.001) {
266 ATH_MSG_ERROR(
"Scale factor must be between 0.001 and 1000.");
274 if (
newName.empty() && scaleFactor == 1.)
return origMaterial;
276 const GeoMaterial* newMaterial =
nullptr;
279 double density = origMaterial->getDensity() * scaleFactor;
280 std::string newName2 =
newName;
281 if (newName2.empty()) {
283 int scaleInt =
static_cast<int>(scaleFactor * 10000);
284 int scale1 = scaleInt / 10000;
285 int scale2 = scaleInt % 10000;
287 std::ostringstream
os;
288 os << origMaterialName <<
scale1 <<
"_" << std::setw(4) << std::setfill(
'0') <<
scale2;
300 GeoIntrusivePtr<const GeoMaterial> matPtr{material};
301 std::string
name(material->getName());
310 ATH_MSG_DEBUG(
"Created new material: " <<
name <<
", " << material->getDensity() /
317 return(std::abs(
d1 /
d2 - 1.) < 1
e-5);
322 ATH_MSG_DEBUG(
"Reading in weight table: " << weightTable->nodeName());
325 ATH_MSG_DEBUG(
"GeometryDBSvc not available. Using old version.");
330 std::string materialName =
db()->
getString(weightTable,
"MATERIAL",
i);
331 if (!space.empty()) {
332 materialName = space +
"::" + materialName;
334 std::string materialBase;
335 if (
db()->testField(weightTable,
"BASEMATERIAL",
i)) {
336 materialBase =
db()->
getString(weightTable,
"BASEMATERIAL",
i);
341 bool linearWeightFlag =
false;
343 linearWeightFlag =
db()->
getInt(weightTable,
"LINWEIGHTFLAG",
i);
347 ATH_MSG_WARNING(
"Material: " << materialName <<
" already exists in weight table");
351 <<
" linearWeightFlag " << linearWeightFlag
352 <<
" raw weight " <<
db()->getDouble(weightTable,
"WEIGHT",
i)
354 <<
" to weight table");
362 int linearWeightFlag) {
367 ATH_MSG_WARNING(
"Material: " << materialName <<
" already exists in weight table");
371 <<
" linearWeightFlag " << linearWeightFlag
372 <<
" to weight table");
379 for (
unsigned int i = 0;
i < weightTable->size();
i++) {
381 std::string materialName = record->
getString(
"MATERIAL");
382 if (!space.empty()) {
383 materialName = space +
"::" + materialName;
385 std::string materialBase;
387 materialBase = record->
getString(
"BASEMATERIAL");
392 bool linearWeightFlag =
false;
394 linearWeightFlag = record->
getInt(
"LINWEIGHTFLAG");
398 ATH_MSG_WARNING(
"Material: " << materialName <<
" already exists in weight table");
407 ATH_MSG_DEBUG(
"Reading in composition table: " << compositionTable->nodeName());
410 ATH_MSG_ERROR(
"GeometryDBSvc not available. Unable to read in composition table.");
413 std::string materialName =
db()->
getString(compositionTable,
"MATERIAL",
i);
414 if (!space.empty()) {
415 materialName = space +
"::" + materialName;
418 std::string componentName =
db()->
getString(compositionTable,
"COMPONENT",
i);
420 double factor =
db()->
getDouble(compositionTable,
"FACTOR",
i);
421 double actualLength =
db()->
getDouble(compositionTable,
"ACTUALLENGTH",
i);
432 if (!scalingTable)
return;
434 if (
db()->getTableSize(scalingTable) == 0)
return;
436 ATH_MSG_DEBUG(
"Reading in extra material scaling table: " << scalingTable->nodeName());
438 ATH_MSG_ERROR(
"GeometryDBSvc not available. Unable to read in scaling table.");
441 std::string materialName =
db()->
getString(scalingTable,
"MATERIAL",
i);
442 double scalingFactor =
db()->
getDouble(scalingTable,
"FACTOR",
i);
445 if (scalingFactor >= 0 || scalingFactor == 1) {
446 msg(
MSG::DEBUG) <<
"Material " << materialName <<
" will be scaled by: " << scalingFactor <<
endmsg;
453 ATH_MSG_WARNING(
"Overriding material: " << materialName <<
" which already exists in scaling table");
477 MaterialWeightMap::const_iterator iter;
479 const std::string& materialBase = iter->second.name;
480 double weight = iter->second.weight;
481 double density =
weight / volume;
482 if (iter->second.linearWeightFlag) {
483 ATH_MSG_ERROR(
"Material defined by linear weight cannot be created with getMaterialForVolume method: " <<
487 ATH_MSG_VERBOSE(
"Found material in weight table - name, base, weight(g), volume(cm3), density(g/cm3): "
488 << materialName <<
", "
489 << materialBase <<
", "
494 if (materialBase.empty()) {
498 return getMaterial(materialBase, density, materialName);
506 ATH_MSG_VERBOSE(
"Material not in weight table, using standard material: "
532 if (volume <= 0 ||
length <= 0) {
538 std::pair<MaterialCompositionMap::const_iterator, MaterialCompositionMap::const_iterator> iterRange;
541 ATH_MSG_VERBOSE(
"Found material in material composition table:" << materialName);
543 std::vector<double> factors;
544 std::vector<std::string> components;
545 for (MaterialCompositionMap::const_iterator iter = iterRange.first; iter != iterRange.second; ++iter) {
546 double factorTmp = iter->second.factor;
547 if (iter->second.actualLength > 0) factorTmp *= iter->second.actualLength /
length;
548 factors.push_back(factorTmp);
549 components.push_back(iter->second.name);
555 MaterialWeightMap::const_iterator iter;
557 const std::string& materialBase = iter->second.name;
558 double weight = iter->second.weight;
559 double density =
weight / volume;
562 if (materialBase.empty()) {
570 ATH_MSG_VERBOSE(
"Material not in weight table, using standard material: "
579 const std::string& materialComponent,
583 std::vector<std::string> tmpMaterialComponents(1, materialComponent);
584 std::vector<double> tmpFactors(1, factor);
590 const std::vector<std::string>& materialComponents,
591 const std::vector<double>& factors,
595 if (volume <= 0 ||
length <= 0) {
600 if (!factors.empty() && factors.size() < materialComponents.size()) {
601 ATH_MSG_WARNING(
"getMaterialForVolumeLength: factor vector size too small. Setting remaining factors to 1.");
604 std::vector<std::string> baseMaterials;
605 std::vector<double> fracWeight;
606 baseMaterials.reserve(materialComponents.size());
607 fracWeight.reserve(materialComponents.size());
609 double totWeight = 0;
610 for (
unsigned int iComp = 0; iComp < materialComponents.size(); ++iComp) {
611 const std::string& materialName = materialComponents[iComp];
614 MaterialWeightMap::const_iterator iter;
616 const std::string& materialBase = iter->second.name;
617 double weight = iter->second.weight;
619 if (iComp < factors.size()) {
624 <<
" factor " << factors[iComp]
626 <<
" basMat " << materialBase
627 <<
" linear? " << iter->second.linearWeightFlag);
630 if (materialBase.empty()) {
632 baseMaterials.push_back(materialName);
634 baseMaterials.push_back(materialBase);
636 fracWeight.push_back(
weight);
655 baseMaterials.push_back(materialName);
656 fracWeight.push_back(
weight);
663 for (
unsigned int i = 0;
i < materialComponents.size(); ++
i) {
669 for (
unsigned int i = 0;
i < fracWeight.size(); ++
i) {
670 fracWeight[
i] /= totWeight;
672 double density = totWeight / volume;
719 const std::vector<std::string>& materialComponents,
720 const std::vector<double>& fracWeight,
727 const std::vector<std::string>& materialComponents,
728 const std::vector<double>& fracWeight,
730 const GeoMaterial* newMaterial =
nullptr;
739 newMaterial = material;
741 GeoMaterial* newMaterialTmp =
new GeoMaterial(
name, density);
742 for (
unsigned int i = 0;
i < materialComponents.size();
i++) {
745 newMaterialTmp->add(
const_cast<GeoMaterial*
>(origMaterial), fracWeight[
i]);
747 ATH_MSG_ERROR(
"Material component missing " << materialComponents[
i]);
751 newMaterial = newMaterialTmp;
765 const std::string materialTable =
"ExtraMaterials";
766 const std::string componentsTable =
"ExtraMatComponents";
771 if (!
db() || !
db()->testField(
"",
"TableSize:" + materialTable) || !
db()->getTableSize(materialTable)
772 || !
db()->testField(
"",
"TableSize:" + componentsTable) || !
db()->getTableSize(componentsTable))
return;
775 ATH_MSG_INFO(
"Extra materials being read in from text file.");
777 using MatMap = std::map<std::string, MaterialDef>;
781 for (
unsigned int iMat = 0; iMat <
db()->
getTableSize(materialTable); iMat++) {
782 std::string materialName =
db()->
getString(materialTable,
"NAME", iMat);
784 materials[materialName] =
MaterialDef(materialName, density);
788 for (
unsigned int iComp = 0; iComp <
db()->
getTableSize(componentsTable); iComp++) {
789 std::string materialName =
db()->
getString(componentsTable,
"NAME", iComp);
791 double fracWeight =
db()->
getDouble(componentsTable,
"FRACTION", iComp);
793 if (iter != materials.end()) {
794 iter->second.addComponent(
compName, fracWeight);
803 int matCountLast = -1;
804 bool someUndefined =
true;
808 while (someUndefined && matCount != matCountLast) {
809 matCountLast = matCount;
810 someUndefined =
false;
811 for (
MatMap::iterator iter = materials.begin(); iter != materials.end(); ++iter) {
816 bool compsDefined =
true;
817 for (
unsigned int iComp = 0; iComp < tmpMat.
numComponents(); ++iComp) {
820 if (iter2 != materials.end()) {
821 if (!iter2->second.isCreated()) {
822 compsDefined =
false;
832 someUndefined =
true;
840 ATH_MSG_ERROR(
"Not all materials could be defined due to cyclic definitions");
853 bool byAtomicRatio =
false;
854 if (totWeight > 1.1) {
855 byAtomicRatio =
true;
857 if (material.
compName(
i).find(
"::") != std::string::npos) {
860 ", is assumed to be defined by atomic ratio (due to total fraction > 1)"
861 " but component is not an element: " << material.
compName(
i));
866 ATH_MSG_ERROR(
"Error making material " << material.
name() <<
". Element not found: " <<
870 totWeight += material.
fraction(
i) * element->getA();
874 if (std::abs(totWeight - 1) > 0.01) {
875 ATH_MSG_WARNING(
"Total fractional weight does not sum to 1. Will renormalize. Total = " << totWeight);
879 GeoIntrusivePtr<GeoMaterial> newMaterial{
new GeoMaterial(material.
name(), material.
density())};
883 double fracWeight = material.
fraction(
i) / totWeight;
884 if (material.
compName(
i).find(
"::") == std::string::npos) {
892 fracWeight = material.
fraction(
i) * element->getA() / totWeight;
894 newMaterial->add(
const_cast<GeoElement*
>(element), fracWeight);
907 newMaterial->add(
const_cast<GeoMaterial*
>(materialTmp.get()), fracWeight);
929 m_fractions.push_back(fraction);
936 for (
unsigned int i = 0;
i < m_fractions.size();
i++) {
937 sum += m_fractions[
i];
952 const GeoMaterial* origMaterial) {
962 if (!origMaterial)
throw std::runtime_error(std::string(
"Invalid material: ") + materialName);
967 if (scaleFactor < 0 || scaleFactor == 1 || materialName.find(
"Ether") != std::string::npos)
return origMaterial;
971 std::string
newName = materialName +
"_ExtraScaling";
977 if (newMaterial)
return newMaterial;
980 double density = origMaterial->getDensity() * scaleFactor;
983 GeoMaterial* newMaterialTmp =
new GeoMaterial(
newName, density);
984 newMaterialTmp->add(
const_cast<GeoMaterial*
>(origMaterial), 1.);
986 newMaterial = newMaterialTmp;
999 ExtraScaleFactorMap::const_iterator iter =
m_scalingMap.find(materialName);
1001 return iter->second;
1006 if (iter !=
m_scalingMap.end() && materialName !=
"std::Air" && materialName !=
"std::Vacuum") {
1007 return iter->second;