diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoMaterial.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoMaterial.h
index 1a328fd896f261de29ff60b4c914e0f7abe3935b..f7ee8632ba8280b2b65c1d38d3e10b14b7b436cb 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoMaterial.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoMaterial.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
*/
#ifndef GEOMODELKERNEL_GEOMATERIAL_H
@@ -28,12 +28,13 @@
#include "GeoModelKernel/RCBase.h"
#include "GeoModelKernel/GeoElement.h"
+#include "GeoModelKernel/GeoIntrusivePtr.h"
#include <vector>
class GeoMaterial : public RCBase
{
public:
- GeoMaterial (const std::string &Name, double Density);
+ GeoMaterial (const std::string& name, double density);
// Add an element to the material.
void add (const GeoElement* element, double fraction = 1.0);
@@ -63,7 +64,8 @@ class GeoMaterial : public RCBase
// Return the nuclear interaction length, computed from the
// density, the list of constituents, and their properties.
double getIntLength () const;
-
+
+ // Return the number of elements
unsigned int getNumElements () const;
// Gets the ith element.
@@ -83,45 +85,44 @@ class GeoMaterial : public RCBase
const unsigned int& getID () const;
protected:
- virtual ~GeoMaterial();
+ virtual ~GeoMaterial() = default;
private:
- GeoMaterial(const GeoMaterial &right);
- GeoMaterial & operator=(const GeoMaterial &right);
+ GeoMaterial() = delete;
+ GeoMaterial(const GeoMaterial &right) = delete;
+ GeoMaterial & operator=(const GeoMaterial &right) = delete;
- std::string m_name;
- double m_density;
- unsigned int m_iD;
+ std::string m_name{};
+ double m_density{0.};
+ unsigned int m_iD{0};
// A list of the fractional composition of each material.
// Fraction is by mass.
- std::vector<double> m_fraction;
+ std::vector<double> m_fractions;
// The radiation length of the material.
- double m_radLength;
+ double m_radLength{0.};
// The interaction length of the material.
- double m_intLength;
-
- // A static used to assign unique identifiers to new
- // materials.
- static unsigned int s_lastID;
+ double m_intLength{0.};
// The constant term in the formula governing dE/dx.
- double m_dedDxConst;
+ double m_dedDxConst{0.};
// The ionization potential in the formula governing dE/dx.
- double m_deDxI0;
+ double m_deDxI0{0.};
// A flag used to lock the material from further addition
// of elements or other constituents.
- bool m_locked;
-
- static const double s_ionizationPotential[93];
+ bool m_locked{false};
- private:
// The list of GeoElements composing a GeoMaterial.
- std::vector<const GeoElement *> m_element;
+ std::vector<GeoIntrusivePtr<const GeoElement>> m_elements;
+
+ // A static used to assign unique identifiers to new materials.
+ static unsigned int s_lastID;
+
+ static const double s_ionizationPotential[93];
};
inline const std::string& GeoMaterial::getName () const
diff --git a/GeoModelCore/GeoModelKernel/src/GeoMaterial.cxx b/GeoModelCore/GeoModelKernel/src/GeoMaterial.cxx
index 26f3e09e7bd68d418f91269db58660da4ec3269d..1b693c83509ae75b1d4b5386323765fb3ee8fc12 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoMaterial.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoMaterial.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
*/
#include "GeoModelKernel/GeoMaterial.h"
@@ -213,61 +213,37 @@ const double
unsigned int GeoMaterial::s_lastID = 0;
-GeoMaterial::GeoMaterial (const std::string &Name, double Density)
- : m_name(Name)
- , m_density(Density)
+GeoMaterial::GeoMaterial (const std::string& name, double density)
+ : m_name(name)
+ , m_density(density)
, m_iD(s_lastID++)
- , m_radLength(0)
- , m_intLength(0)
- , m_dedDxConst(0)
- , m_deDxI0(0)
- , m_locked(false)
{
}
-GeoMaterial::~GeoMaterial()
-{
- for (size_t i = 0; i < m_element.size (); i++)
- {
- m_element[i]->unref ();
- }
-}
-
void GeoMaterial::add (const GeoElement* element, double fraction)
{
// You can only add materials until you call "lock"...
- if (!m_locked)
- {
- std::vector <const GeoElement *>::iterator e = std::find(m_element.begin(),m_element.end(),element);
- if (e==m_element.end()) {
- m_element.push_back (element);
- m_fraction.push_back (fraction);
- element->ref ();
- }
- else {
- int n = e-m_element.begin();
- m_fraction[n]+=fraction;
- }
- }
- else
- {
- throw std::out_of_range ("Element added after material locked");
- }
+ if(m_locked) throw std::out_of_range ("Element added after material locked");
+
+ auto e = std::find(m_elements.begin(),m_elements.end(),element);
+ if (e==m_elements.end()) {
+ m_elements.push_back (element);
+ m_fractions.push_back (fraction);
+ element->ref ();
+ }
+ else {
+ int n = e-m_elements.begin();
+ m_fractions[n]+=fraction;
+ }
}
void GeoMaterial::add (const GeoMaterial* material, double fraction)
{
- if (!m_locked)
- {
- for (size_t e = 0; e < material->getNumElements (); e++)
- {
- add(material->m_element[e],fraction * material->m_fraction[e]);
- }
- }
- else
- {
- throw std::out_of_range ("Material added after material locked");
- }
+ if(m_locked) throw std::out_of_range ("Material added after material locked");
+
+ for(size_t e = 0; e < material->getNumElements (); ++e) {
+ add(material->m_elements[e],fraction * material->m_fractions[e]);
+ }
}
void GeoMaterial::lock ()
@@ -292,12 +268,7 @@ void GeoMaterial::lock ()
// //
//--------------------------------------------//
- if (getNumElements () == 0)
- {
- throw std::out_of_range ("Attempt to lock a material with no elements "+getName());
- return;
- }
-
+ if(getNumElements() == 0) throw std::out_of_range ("Attempt to lock a material with no elements "+getName());
//--------------------------------------------//
// //
@@ -305,42 +276,44 @@ void GeoMaterial::lock ()
double dEDxConstant = 0, dEDxI0 = 0, NILinv = 0.0, radInv = 0.0;
- { // ===============Renormalization================================
- double wSum=std::accumulate(m_fraction.begin(),m_fraction.end(),0.0);
- if (fabs(wSum-1.0)>FLT_EPSILON) // 'FLT_EPSILON' defined in <cfloat>
- {
- if(std::getenv("GEOMODELKERNEL_VERBOSE")) {
- std::cerr << "Warning in material "
- << m_name
- << ". Mass fractions sum to "
- << wSum << "; renormalizing to 1.0" << std::endl;
- }
+ // ===============Renormalization================================
+ {
+ double wSum=std::accumulate(m_fractions.begin(),m_fractions.end(),0.0);
+ if(fabs(wSum-1.0)>FLT_EPSILON) { // 'FLT_EPSILON' defined in <cfloat>
+ if(std::getenv("GEOMODELKERNEL_VERBOSE")) {
+ std::cerr << "Warning in material "
+ << m_name
+ << ". Mass fractions sum to "
+ << wSum << "; renormalizing to 1.0" << std::endl;
+ }
}
double inv_wSum = 1. / wSum;
- for (size_t e=0;e<getNumElements();e++) {m_fraction[e]*=inv_wSum;}
- } // ==============================================================
+ for(size_t e=0;e<getNumElements();++e) {
+ m_fractions[e]*=inv_wSum;
+ }
+ }
+ // ==============================================================
const double inv_lambda0 = 1. / lambda0;
- for (size_t e = 0; e < getNumElements (); e++)
- {
- double w = getFraction (e); // Weight fraction.
- double Z = m_element[e]->getZ (); // Atomic number
- double A = m_element[e]->getA (); // Atomic mass.
- double N = m_element[e]->getN (); // Number of nucleons.
- double dovera = m_density ? m_density / A : 0; // don't crash if both are 0
- double n = m_fraction[e] * GeoModelKernelUnits::Avogadro * dovera; // Number density.
- int iZ = (int) (m_element[e]->getZ () + 0.5) - 1; // Atomic number index
-
- dEDxConstant += w * C0 * dovera * Z;
- // Make sure we don't overflow the table:
- // the `Ether' special `element' has Z set to 500.
- if (iZ >= 0 && iZ < (std::end(s_ionizationPotential)-std::begin(s_ionizationPotential)))
- dEDxI0 += w * s_ionizationPotential[iZ];
- NILinv += n * pow (N, 2.0 / 3.0) * GeoModelKernelUnits::amu * inv_lambda0;
-
- double nAtomsPerVolume = A ? GeoModelKernelUnits::Avogadro*m_density*m_fraction[e]/A : 0.;
- radInv += (nAtomsPerVolume*m_element[e]->getRadTsai());
- }
+ for(size_t e = 0; e < getNumElements (); ++e) {
+ double w = getFraction (e); // Weight fraction.
+ double Z = m_elements[e]->getZ (); // Atomic number
+ double A = m_elements[e]->getA (); // Atomic mass.
+ double N = m_elements[e]->getN (); // Number of nucleons.
+ double dovera = m_density ? m_density / A : 0; // don't crash if both are 0
+ double n = m_fractions[e] * GeoModelKernelUnits::Avogadro * dovera; // Number density.
+ int iZ = (int) (m_elements[e]->getZ () + 0.5) - 1; // Atomic number index
+
+ dEDxConstant += w * C0 * dovera * Z;
+ // Make sure we don't overflow the table:
+ // the `Ether' special `element' has Z set to 500.
+ if (iZ >= 0 && iZ < (std::end(s_ionizationPotential)-std::begin(s_ionizationPotential)))
+ dEDxI0 += w * s_ionizationPotential[iZ];
+ NILinv += n * pow (N, 2.0 / 3.0) * GeoModelKernelUnits::amu * inv_lambda0;
+
+ double nAtomsPerVolume = A ? GeoModelKernelUnits::Avogadro*m_density*m_fractions[e]/A : 0.;
+ radInv += (nAtomsPerVolume*m_elements[e]->getRadTsai());
+ }
m_dedDxConst = dEDxConstant;
m_deDxI0 = dEDxI0 ;
m_intLength = NILinv ? 1.0 / NILinv : 0;
@@ -382,7 +355,6 @@ double GeoMaterial::getDeDxMin () const
// -----------------------------------------------------------//
return m_dedDxConst * ConstToMin;
-
}
double GeoMaterial::getRadLength () const
@@ -403,19 +375,19 @@ unsigned int GeoMaterial::getNumElements () const
{
if (!m_locked)
throw std::out_of_range ("Material accessed before lock");
- return m_element.size ();
+ return m_elements.size();
}
const GeoElement* GeoMaterial::getElement (unsigned int i) const
{
if (!m_locked)
throw std::out_of_range ("Material accessed before lock");
- return m_element[i];
+ return m_elements[i];
}
double GeoMaterial::getFraction (int i) const
{
if (!m_locked)
throw std::out_of_range ("Material accessed before lock");
- return m_fraction[i];
+ return m_fractions[i];
}