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