From 6327d3f0e0fb9230180f5d800406ed383a422aa6 Mon Sep 17 00:00:00 2001
From: Johannes Junggeburth <johannes.josef.junggeburth@cern.ch>
Date: Fri, 24 Jan 2025 15:19:15 +0100
Subject: [PATCH] GeoMaterial - Minor clean-up
---
.../GeoModelKernel/GeoElement.h | 67 +++------
.../GeoModelKernel/GeoMaterial.h | 71 +++++----
.../GeoModelKernel/src/GeoElement.cxx | 27 ++--
.../GeoModelKernel/src/GeoMaterial.cxx | 136 +++++++++---------
4 files changed, 131 insertions(+), 170 deletions(-)
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoElement.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoElement.h
index f8cfeee60..b236915ca 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoElement.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoElement.h
@@ -26,42 +26,45 @@ class GeoElement : public RCBase
int operator!=(const GeoElement &right) const;
// Returns the effective number of nucleons in the element.
- double getN () const;
+ double getN () const{
+ return m_a * (GeoModelKernelUnits::mole / GeoModelKernelUnits::gram);
+ }
// The name of the element, e.g. "Carbon".
- const std::string& getName () const;
+ const std::string& getName () const {
+ return m_name;
+ }
// The chemical symbol for the element, e.g. C, O, S, Na....
- const std::string& getSymbol () const;
+ const std::string& getSymbol () const {
+ return m_symbol;
+ }
// The atomic number Z for the material.
- const double& getZ () const;
+ double getZ() const{
+ return m_z;
+ }
// The average atomic mass for the element.
- const double& getA () const;
+ double getA() const {
+ return m_a;
+ }
// Tsai formula for the radiation length
- double getRadTsai () const;
+ double getRadTsai() const;
protected:
- virtual ~GeoElement();
+ virtual ~GeoElement() = default;
private:
- GeoElement(const GeoElement &right);
- GeoElement & operator=(const GeoElement &right);
-
- std::string m_name;
- std::string m_symbol;
-
- double m_z;
- double m_a;
+ std::string m_name{};
+ std::string m_symbol{};
+ double m_z{0.};
+ double m_a{0.};
};
-inline int GeoElement::operator==(const GeoElement &right) const
-{
- return
- m_name ==
- right.m_name && m_symbol == right.m_symbol && m_z == right.m_z && m_a == right.m_a;
+inline int GeoElement::operator==(const GeoElement &right) const{
+ return m_name == right.m_name && m_symbol == right.m_symbol && m_z == right.m_z && m_a == right.m_a;
}
inline int GeoElement::operator!=(const GeoElement &right) const
@@ -70,29 +73,5 @@ inline int GeoElement::operator!=(const GeoElement &right) const
m_symbol != right.m_symbol || m_z != right.m_z || m_a != right.m_a;
}
-inline double GeoElement::getN () const
-{
- return m_a * (GeoModelKernelUnits::mole / GeoModelKernelUnits::gram);
-}
-
-inline const std::string& GeoElement::getName () const
-{
- return m_name;
-}
-
-inline const std::string& GeoElement::getSymbol () const
-{
- return m_symbol;
-}
-
-inline const double& GeoElement::getZ () const
-{
- return m_z;
-}
-
-inline const double& GeoElement::getA () const
-{
- return m_a;
-}
#endif
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoMaterial.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoMaterial.h
index 78a1e6e38..8050cbf9f 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoMaterial.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoMaterial.h
@@ -31,95 +31,92 @@
#include "GeoModelKernel/GeoIntrusivePtr.h"
#include <vector>
-class GeoMaterial : public RCBase
-{
+class GeoMaterial : public RCBase {
public:
GeoMaterial (const std::string& name, double density);
- // Add an element to the material.
+ // Add an element to the material.
void add (const GeoElement* element, double fraction = 1.0);
- // Add another material to the material (this copies all of
- // the element from the other material into this one).
- // Fraction is by mass.
+ // Add another material to the material (this copies all of
+ // the element from the other material into this one).
+ // Fraction is by mass.
void add (const GeoMaterial* material, double fraction);
- // Lock the material against the addition of other
- // materials or elements.
+ // Lock the material against the addition of other
+ // materials or elements.
void lock ();
- // Constant dEdx term.
+ // Constant dEdx term.
double getDeDxConstant () const;
- // I0 term (ionization potential). From tables. Units: eV.
+ // I0 term (ionization potential). From tables. Units: eV.
double getDeDxI0 () const;
- // Get dEdx_min. dEdxConstant*11.528 Paul Avery, CBX-92-39.
+ // Get dEdx_min. dEdxConstant*11.528 Paul Avery, CBX-92-39.
double getDeDxMin () const;
- // Returns the radiation length, computed from the density
- // and the list of constituents, and their properties.
+ // Returns the radiation length, computed from the density
+ // and the list of constituents, and their properties.
double getRadLength () const;
- // Return the nuclear interaction length, computed from the
- // density, the list of constituents, and their properties.
+ // Return the nuclear interaction length, computed from the
+ // density, the list of constituents, and their properties.
double getIntLength () const;
- // Return the number of elements
+ // Return the number of elements
unsigned int getNumElements () const;
- // Gets the ith element.
+ // Gets the ith element.
const GeoElement* getElement (unsigned int i) const;
- // Gets the fraction by weight of the ith element
+ // Gets the fraction by weight of the ith element
double getFraction (int i) const;
- // The name of the material.
+ // The name of the material.
const std::string& getName () const;
- // The density of the material.
+ // The density of the material.
const double& getDensity () const;
- // Gives an integral identifier for the material. For
- // convenience.
+ // Gives an integral identifier for the material. For
+ // convenience.
const unsigned int& getID () const;
protected:
virtual ~GeoMaterial() = default;
private:
- GeoMaterial() = delete;
- GeoMaterial(const GeoMaterial &right) = delete;
- GeoMaterial & operator=(const GeoMaterial &right) = delete;
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_fractions;
- // The radiation length of the material.
+ // The radiation length of the material.
double m_radLength{0.};
- // The interaction length of the material.
+ // The interaction length of the material.
double m_intLength{0.};
- // The constant term in the formula governing dE/dx.
+ // The constant term in the formula governing dE/dx.
double m_dedDxConst{0.};
- // The ionization potential in the formula governing dE/dx.
+ // The ionization potential in the formula governing dE/dx.
double m_deDxI0{0.};
- // A flag used to lock the material from further addition
- // of elements or other constituents.
+ // A flag used to lock the material from further addition
+ // of elements or other constituents.
bool m_locked{false};
- // The list of GeoElements composing a GeoMaterial.
- std::vector<GeoIntrusivePtr<const GeoElement>> m_elements;
+ // A list of the fractional composition of each material.
+ // Fraction is by mass.
+
+ /// The list of GeoElements composing a GeoMaterial with the corresponding fraction.
+ using ElementWithFrac = std::pair<GeoIntrusivePtr<const GeoElement>, double>;
+ std::vector<ElementWithFrac> m_elements{};
- // A static used to assign unique identifiers to new materials.
+ // A static used to assign unique identifiers to new materials.
static std::atomic<unsigned int> s_lastID;
};
diff --git a/GeoModelCore/GeoModelKernel/src/GeoElement.cxx b/GeoModelCore/GeoModelKernel/src/GeoElement.cxx
index d2225aeb4..c73fe32e6 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoElement.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoElement.cxx
@@ -4,32 +4,26 @@
#include "GeoModelKernel/GeoElement.h"
#include <cmath>
+#include <array>
-GeoElement::GeoElement (const std::string &Name, const std::string &Symbol, double Z, double A)
- : m_name (Name)
- , m_symbol (Symbol)
- , m_z (Z)
- , m_a (A)
-{
-}
-
-GeoElement::~GeoElement()
-{
-}
-
+GeoElement::GeoElement (const std::string &Name, const std::string &Symbol, double Z, double A):
+ m_name{Name},
+ m_symbol{Symbol},
+ m_z {Z},
+ m_a {A} {}
double GeoElement::getRadTsai() const
{
// Compute Coulomb correction factor (Phys Rev. D50 3-1 (1994) page 1254)
- const double k1 = 0.0083, k2 = 0.20206, k3 = 0.0020, k4 = 0.0369;
+ constexpr double k1 = 0.0083, k2 = 0.20206, k3 = 0.0020, k4 = 0.0369;
double az2 = (GeoModelKernelUnits::fine_structure_const*m_z)*(GeoModelKernelUnits::fine_structure_const*m_z);
double az4 = az2 * az2;
double coulomb = (k1*az4 + k2 + 1./(1.+az2))*az2 - (k3*az4 + k4)*az4;
// Compute Tsai's Expression for the Radiation Length (Phys Rev. D50 3-1 (1994) page 1254)
- const double Lrad_light[] = {5.31, 4.79, 4.74, 4.71};
- const double Lprad_light[] = {6.144, 5.621, 5.805, 5.924};
+ constexpr std::array<double,4> Lrad_light {5.31, 4.79, 4.74, 4.71};
+ constexpr std::array<double,4> Lprad_light {6.144, 5.621, 5.805, 5.924};
const double logZ3 = std::log(m_z) * (1./3);
@@ -43,6 +37,5 @@ double GeoElement::getRadTsai() const
Lprad = std::log(1194.) - 2*logZ3;
}
- double radTsai = 4*GeoModelKernelUnits::alpha_rcl2*m_z*(m_z*(Lrad-coulomb) + Lprad);
- return radTsai;
+ return 4.*GeoModelKernelUnits::alpha_rcl2*m_z*(m_z*(Lrad-coulomb) + Lprad);
}
diff --git a/GeoModelCore/GeoModelKernel/src/GeoMaterial.cxx b/GeoModelCore/GeoModelKernel/src/GeoMaterial.cxx
index af6582782..e790a1c61 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoMaterial.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoMaterial.cxx
@@ -3,6 +3,7 @@
*/
#include "GeoModelKernel/GeoMaterial.h"
+#include "GeoModelKernel/throwExcept.h"
#include <stdexcept>
#include <cmath>
#include <numeric>
@@ -213,41 +214,36 @@ constexpr std::array<double, 93> s_ionizationPotential{
std::atomic<unsigned int> GeoMaterial::s_lastID = 0;
-GeoMaterial::GeoMaterial (const std::string& name, double density)
- : m_name(name)
- , m_density(density)
- , m_iD(s_lastID++)
-{
-}
+GeoMaterial::GeoMaterial (const std::string& name, double density):
+ m_name{name},
+ m_density{density},
+ m_iD{s_lastID++}{}
-void GeoMaterial::add (const GeoElement* element, double fraction)
-{
+void GeoMaterial::add (const GeoElement* element, double fraction) {
// You can only add materials until you call "lock"...
- if(m_locked) throw std::out_of_range ("Element added after material locked");
+ if(m_locked) THROW_EXCEPTION("Element added after material locked");
GeoIntrusivePtr<const GeoElement> elementPtr{element};
- auto e = std::find(m_elements.begin(),m_elements.end(),element);
+ std::vector<ElementWithFrac>::iterator e = std::find_if(m_elements.begin(),m_elements.end(),
+ [elementPtr](ElementWithFrac& known){
+ return known.first == elementPtr;
+ });
if (e==m_elements.end()) {
- m_elements.push_back (elementPtr);
- m_fractions.push_back (fraction);
- }
- else {
- int n = e-m_elements.begin();
- m_fractions[n]+=fraction;
+ m_elements.emplace_back(std::make_pair(elementPtr, fraction));
+ } else {
+ e->second+=fraction;
}
}
-void GeoMaterial::add (const GeoMaterial* material, double fraction)
-{
- if(m_locked) throw std::out_of_range ("Material added after material locked");
+void GeoMaterial::add (const GeoMaterial* material, double fraction) {
+ if(m_locked) THROW_EXCEPTION ("Material added after material locked");
- for(size_t e = 0; e < material->getNumElements (); ++e) {
- add(material->m_elements[e],fraction * material->m_fractions[e]);
+ for(const auto& [element, weight] : material->m_elements) {
+ add(element, fraction * weight);
}
}
-void GeoMaterial::lock ()
-{
+void GeoMaterial::lock () {
if(m_locked) return;
m_locked = true;
@@ -261,48 +257,52 @@ void GeoMaterial::lock ()
// //
// For energy loss //
// //
- const double C0 = 0.00307 * GeoModelKernelUnits::cm3 / GeoModelKernelUnits::gram; //
+ constexpr double C0 = 0.00307 * GeoModelKernelUnits::cm3 / GeoModelKernelUnits::gram; //
// //
// For nuclear absorption length. //
- const double lambda0 = 35 * GeoModelKernelUnits::gram / GeoModelKernelUnits::cm2; //
+ constexpr double lambda0 = 35 * GeoModelKernelUnits::gram / GeoModelKernelUnits::cm2; //
// //
//--------------------------------------------//
- if(getNumElements() == 0) throw std::out_of_range ("Attempt to lock a material with no elements "+getName());
+ if(getNumElements() == 0) THROW_EXCEPTION("Attempt to lock a material with no elements "<<getName());
//--------------------------------------------//
// //
// -------------------------------------------//
- double dEDxConstant = 0, dEDxI0 = 0, NILinv = 0.0, radInv = 0.0;
+ double dEDxConstant{0.}, dEDxI0{0.}, NILinv{0.}, radInv{0.};
+// std::sort(m_elements.begin(), m_elements.end(),[](const ElementWithFrac& a, const ElementWithFrac& b){
+// return a.second > b.second;
+// });
// ===============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>
+
+ const double wSum=std::accumulate(m_elements.begin(),m_elements.end(),0.0,
+ [](const double w, const ElementWithFrac& a) {
+ return a.second +w;
+ });
+ if(std::abs(wSum-1.0)>std::numeric_limits<float>::epsilon()) {
if(std::getenv("GEOMODELKERNEL_VERBOSE")) {
- std::cerr << "Warning in material "
- << m_name
- << ". Mass fractions sum to "
- << wSum << "; renormalizing to 1.0" << std::endl;
+ 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_fractions[e]*=inv_wSum;
+ const double inv_wSum = 1. / wSum;
+ for(auto& [element, weight] : m_elements) {
+ weight*=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_elements[e]->getZ (); // Atomic number
- double A = m_elements[e]->getA (); // Atomic mass.
- double N = m_elements[e]->getN (); // Number of nucleons.
+ for(const auto& [element, w] : m_elements) {
+ double Z = element->getZ(); // Atomic number
+ double A = element->getA(); // Atomic mass.
+ double N = element->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
+ double n = element * GeoModelKernelUnits::Avogadro * dovera; // Number density.
+ int iZ = (int) (element->getZ () + 0.5) - 1; // Atomic number index
dEDxConstant += w * C0 * dovera * Z;
// Make sure we don't overflow the table:
@@ -311,8 +311,8 @@ void GeoMaterial::lock ()
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());
+ double nAtomsPerVolume = A ? GeoModelKernelUnits::Avogadro*m_density*w/A : 0.;
+ radInv += (nAtomsPerVolume*element->getRadTsai());
}
m_dedDxConst = dEDxConstant;
m_deDxI0 = dEDxI0 ;
@@ -326,28 +326,25 @@ void GeoMaterial::lock ()
m_radLength = radInv ? 1.0 / radInv : 0;
}
-double GeoMaterial::getDeDxConstant () const
-{
+double GeoMaterial::getDeDxConstant () const {
if (!m_locked)
- throw std::out_of_range ("Material accessed before lock");
+ THROW_EXCEPTION ("Material accessed before lock");
return m_dedDxConst;
}
-double GeoMaterial::getDeDxI0 () const
-{
+double GeoMaterial::getDeDxI0 () const {
if (!m_locked)
- throw std::out_of_range ("Material accessed before lock");
+ THROW_EXCEPTION ("Material accessed before lock");
return m_deDxI0;
}
-double GeoMaterial::getDeDxMin () const
-{
+double GeoMaterial::getDeDxMin () const {
//------------------------------------------------------------//
- static const double ConstToMin = 11.528; //
+ constexpr double ConstToMin = 11.528; //
//------------------------------------------------------------//
if (!m_locked)
- throw std::out_of_range ("Material accessed before lock");
+ THROW_EXCEPTION ("Material accessed before lock");
// -----------------------------------------------------------//
// See: Paul Avery's notes on track fitting, CBX 92-39. //
@@ -357,37 +354,32 @@ double GeoMaterial::getDeDxMin () const
return m_dedDxConst * ConstToMin;
}
-double GeoMaterial::getRadLength () const
-{
+double GeoMaterial::getRadLength () const {
if (!m_locked)
- throw std::out_of_range ("Material accessed before lock");
+ THROW_EXCEPTION ("Material accessed before lock");
return m_radLength;
}
-double GeoMaterial::getIntLength () const
-{
+double GeoMaterial::getIntLength () const {
if (!m_locked)
- throw std::out_of_range ("Material accessed before lock");
+ THROW_EXCEPTION ("Material accessed before lock");
return m_intLength;
}
-unsigned int GeoMaterial::getNumElements () const
-{
+unsigned int GeoMaterial::getNumElements() const {
if (!m_locked)
- throw std::out_of_range ("Material accessed before lock");
+ THROW_EXCEPTION ("Material accessed before lock");
return m_elements.size();
}
-const GeoElement* GeoMaterial::getElement (unsigned int i) const
-{
+const GeoElement* GeoMaterial::getElement(unsigned int i) const {
if (!m_locked)
- throw std::out_of_range ("Material accessed before lock");
- return m_elements[i];
+ THROW_EXCEPTION ("Material accessed before lock");
+ return m_elements[i].first;
}
-double GeoMaterial::getFraction (int i) const
-{
+double GeoMaterial::getFraction (int i) const {
if (!m_locked)
- throw std::out_of_range ("Material accessed before lock");
- return m_fractions[i];
+ THROW_EXCEPTION ("Material accessed before lock");
+ return m_elements[i].second;
}
--
GitLab