GeoMaterial - Minor clean-up

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

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;
   
 };

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);
 }

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;
 }