Acts/Material/MaterialComposition.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <algorithm>
0012 #include <cassert>
0013 #include <cstdint>
0014 #include <limits>
0015 #include <ostream>
0016 #include <vector>
0017
0018 namespace Acts {
0019
0020 /// Memory-efficient storage of the relative fraction of an element.
0021 ///
0022 /// @ingroup material
0023 ///
0024 /// This can be used to define materials that are compounds of multiple elements
0025 /// with varying fractions. The element is identified by its atomic number
0026 /// stored as a single byte (allows up to 256 elements; more than we need).
0027 /// Its fraction is also stored as a single byte with values between 0 and
0028 /// 255. This gives an accuracy of 1/256 ~ 0.5 %.
0029 ///
0030 /// The element fraction allows you to store element composition in merged
0031 /// materials with a large number of bins. Depending on the
0032 /// detector and the description granularity this can be a lot of information
0033 /// and thus requires the reduced memory footprint. This is really only needed
0034 /// for nuclear interaction in the fast simulation where the reduced fractional
0035 /// accuracy is not a problem. The fractional accuracy should be much better
0036 /// than the parametrization uncertainty for hadronic interactions.
0037 class ElementFraction {
0038  public:
0039   /// Construct from atomic number and relative fraction.
0040   ///
0041   /// @param e is the atomic number of the element
0042   /// @param f is the relative fraction and must be a value in [0,1]
0043   constexpr ElementFraction(unsigned int e, float f)
0044       : m_element(static_cast<std::uint8_t>(e)),
0045         m_fraction(static_cast<std::uint8_t>(
0046             f * std::numeric_limits<std::uint8_t>::max())) {
0047     assert((0u < e) && ("The atomic number must be positive"));
0048     assert((0.0f <= f) && (f <= 1.0f) && "Relative fraction must be in [0,1]");
0049   }
0050   /// Construct from atomic number and integer weight.
0051   ///
0052   /// @param e is the atomic number of the element
0053   /// @param w is the integer weight and must be a value in [0,256)
0054   constexpr explicit ElementFraction(unsigned int e, unsigned int w)
0055       : m_element(static_cast<std::uint8_t>(e)),
0056         m_fraction(static_cast<std::uint8_t>(w)) {
0057     assert((0u < e) && ("The atomic number must be positive"));
0058     assert((w < 256u) && "Integer weight must be in [0,256)");
0059   }
0060
0061   /// Must always be created with valid data.
0062   ElementFraction() = delete;
0063   /// Move constructor
0064   ElementFraction(ElementFraction&&) = default;
0065   /// Copy constructor
0066   ElementFraction(const ElementFraction&) = default;
0067   ~ElementFraction() = default;
0068   /// Move assignment operator
0069   /// @return Reference to this element fraction after move assignment
0070   ElementFraction& operator=(ElementFraction&&) = default;
0071   /// Copy assignment operator
0072   /// @return Reference to this element fraction after copy assignment
0073   ElementFraction& operator=(const ElementFraction&) = default;
0074
0075   /// The element atomic number.
0076   /// @return The atomic number of the element
0077   constexpr std::uint8_t element() const { return m_element; }
0078   /// The relative fraction of this element.
0079   /// @return The relative fraction as a float in [0,1]
0080   constexpr float fraction() const {
0081     return static_cast<float>(m_fraction) /
0082            std::numeric_limits<std::uint8_t>::max();
0083   }
0084
0085  private:
0086   // element atomic number
0087   std::uint8_t m_element;
0088   // element fraction in the compound scaled to the [0,256) range.
0089   std::uint8_t m_fraction;
0090
0091   friend constexpr bool operator==(ElementFraction lhs, ElementFraction rhs) {
0092     return (lhs.m_fraction == rhs.m_fraction) &&
0093            (lhs.m_element == rhs.m_element);
0094   }
0095   /// Sort by fraction for fastest access to the most probable element.
0096   friend constexpr bool operator<(ElementFraction lhs, ElementFraction rhs) {
0097     return lhs.m_fraction < rhs.m_fraction;
0098   }
0099   friend class MaterialComposition;
0100
0101   /// Stream operator for ElementFraction
0102   friend std::ostream& operator<<(std::ostream& os, const ElementFraction& ef) {
0103     os << "ElementFraction(Z=" << static_cast<unsigned int>(ef.m_element)
0104        << ", f=" << ef.fraction() << ")";
0105     return os;
0106   }
0107 };
0108
0109 /// Material composed from multiple elements with varying factions.
0110 ///
0111 /// @see ElementFraction for details.
0112 class MaterialComposition {
0113  public:
0114   /// Construct an empty composition corresponding to vacuum.
0115   MaterialComposition() = default;
0116   /// Constructor from element fractions.
0117   ///
0118   /// Rescales the fractions so they all add up to unity within the accuracy.
0119   /// @param elements Vector of element fractions that define the composition
0120   explicit MaterialComposition(std::vector<ElementFraction> elements)
0121       : m_elements(std::move(elements)) {
0122     std::ranges::sort(m_elements, std::less<ElementFraction>{});
0123     // compute the total weight first
0124     unsigned total = 0u;
0125     for (const auto& element : m_elements) {
0126       total += element.m_fraction;
0127     }
0128     // compute scale factor into the [0, 256) range
0129     float scale = float{std::numeric_limits<std::uint8_t>::max()} / total;
0130     for (auto& element : m_elements) {
0131       element.m_fraction =
0132           static_cast<std::uint8_t>(element.m_fraction * scale);
0133     }
0134   }
0135
0136   /// Move constructor
0137   MaterialComposition(MaterialComposition&&) = default;
0138   /// Copy constructor
0139   MaterialComposition(const MaterialComposition&) = default;
0140   ~MaterialComposition() = default;
0141   /// Move assignment operator
0142   /// @return Reference to this material composition after move assignment
0143   MaterialComposition& operator=(MaterialComposition&&) = default;
0144   /// Copy assignment operator
0145   /// @return Reference to this material composition after copy assignment
0146   MaterialComposition& operator=(const MaterialComposition&) = default;
0147
0148   /// Support range-based iteration over contained elements.
0149   /// @return Iterator to the first element
0150   auto begin() const { return m_elements.begin(); }
0151   /// Get iterator to end of elements
0152   /// @return Iterator past the last element
0153   auto end() const { return m_elements.end(); }
0154
0155   /// Check if the composed material is valid, i.e. it is not vacuum.
0156   explicit operator bool() const { return !m_elements.empty(); }
0157   /// Return the number of elements.
0158   /// @return The number of elements in the composition
0159   std::size_t size() const { return m_elements.size(); }
0160
0161  private:
0162   std::vector<ElementFraction> m_elements;
0163
0164   friend inline bool operator==(const MaterialComposition& lhs,
0165                                 const MaterialComposition& rhs) {
0166     return lhs.m_elements == rhs.m_elements;
0167   }
0168
0169   /// Stream operator for MaterialComposition
0170   friend std::ostream& operator<<(std::ostream& os,
0171                                   const MaterialComposition& mc) {
0172     os << "MaterialComposition(elements=[";
0173     for (std::size_t i = 0; i < mc.m_elements.size(); ++i) {
0174       if (i > 0) {
0175         os << ", ";
0176       }
0177       os << mc.m_elements[i];
0178     }
0179     os << "])";
0180     return os;
0181   }
0182 };
0183
0184 }  // namespace Acts