Acts/Material/MaterialComposition.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018-2020 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 http://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<uint8_t>(e)),
0042         m_fraction(static_cast<uint8_t>(f * UINT8_MAX)) {
0043     assert((0u < e) && ("The atomic number must be positive"));
0044     assert((0.0f <= f) && (f <= 1.0f) && "Relative fraction must be in [0,1]");
0045   }
0046   /// Construct from atomic number and integer weight.
0047   ///
0048   /// @param e is the atomic number of the element
0049   /// @param w is the integer weight and must be a value in [0,256)
0050   constexpr explicit ElementFraction(unsigned int e, unsigned int w)
0051       : m_element(static_cast<uint8_t>(e)),
0052         m_fraction(static_cast<uint8_t>(w)) {
0053     assert((0u < e) && ("The atomic number must be positive"));
0054     assert((w < 256u) && "Integer weight must be in [0,256)");
0055   }
0056
0057   /// Must always be created with valid data.
0058   ElementFraction() = delete;
0059   ElementFraction(ElementFraction&&) = default;
0060   ElementFraction(const ElementFraction&) = default;
0061   ~ElementFraction() = default;
0062   ElementFraction& operator=(ElementFraction&&) = default;
0063   ElementFraction& operator=(const ElementFraction&) = default;
0064
0065   /// The element atomic number.
0066   constexpr uint8_t element() const { return m_element; }
0067   /// The relative fraction of this element.
0068   constexpr float fraction() const {
0069     return static_cast<float>(m_fraction) / UINT8_MAX;
0070   }
0071
0072  private:
0073   // element atomic number
0074   uint8_t m_element;
0075   // element fraction in the compound scaled to the [0,256) range.
0076   uint8_t m_fraction;
0077
0078   friend constexpr bool operator==(ElementFraction lhs, ElementFraction rhs) {
0079     return (lhs.m_fraction == rhs.m_fraction) &&
0080            (lhs.m_element == rhs.m_element);
0081   }
0082   /// Sort by fraction for fastest access to the most probable element.
0083   friend constexpr bool operator<(ElementFraction lhs, ElementFraction rhs) {
0084     return lhs.m_fraction < rhs.m_fraction;
0085   }
0086   friend class MaterialComposition;
0087 };
0088
0089 /// Material composed from multiple elements with varying factions.
0090 ///
0091 /// @see ElementFraction for details.
0092 class MaterialComposition {
0093  public:
0094   /// Construct an empty composition corresponding to vacuum.
0095   MaterialComposition() = default;
0096   /// Constructor from element fractions.
0097   ///
0098   /// Rescales the fractions so they all add up to unity within the accuracy.
0099   MaterialComposition(std::vector<ElementFraction> elements)
0100       : m_elements(std::move(elements)) {
0101     std::sort(m_elements.begin(), m_elements.end());
0102     // compute the total weight first
0103     unsigned total = 0u;
0104     for (auto element : m_elements) {
0105       total += element.m_fraction;
0106     }
0107     // compute scale factor into the [0, 256) range
0108     float scale = float{std::numeric_limits<std::uint8_t>::max()} / total;
0109     for (auto& element : m_elements) {
0110       element.m_fraction =
0111           static_cast<std::uint8_t>(element.m_fraction * scale);
0112     }
0113   }
0114
0115   MaterialComposition(MaterialComposition&&) = default;
0116   MaterialComposition(const MaterialComposition&) = default;
0117   ~MaterialComposition() = default;
0118   MaterialComposition& operator=(MaterialComposition&&) = default;
0119   MaterialComposition& operator=(const MaterialComposition&) = default;
0120
0121   // Support range-based iteration over contained elements.
0122   auto begin() const { return m_elements.begin(); }
0123   auto end() const { return m_elements.end(); }
0124
0125   /// Check if the composed material is valid, i.e. it is not vacuum.
0126   operator bool() const { return !m_elements.empty(); }
0127   /// Return the number of elements.
0128   std::size_t size() const { return m_elements.size(); }
0129
0130  private:
0131   std::vector<ElementFraction> m_elements;
0132
0133   friend inline bool operator==(const MaterialComposition& lhs,
0134                                 const MaterialComposition& rhs) {
0135     return (lhs.m_elements == rhs.m_elements);
0136   }
0137 };
0138
0139 }  // namespace Acts