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