Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:27:38

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 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 "Acts/Geometry/GeometryIdentifier.hpp"
0012 
0013 #include <algorithm>
0014 #include <cassert>
0015 #include <initializer_list>
0016 #include <iterator>
0017 #include <stdexcept>
0018 #include <utility>
0019 #include <vector>
0020 
0021 namespace Acts {
0022 
0023 /// Store values mapped into the geometry hierarchy.
0024 ///
0025 /// @tparam value_t stored value type
0026 ///
0027 /// The core functionality is to find an equivalent element, i.e. an
0028 /// identifier-value pair, for a given geometry identifier via
0029 ///
0030 ///     auto it = container.find(GeometryIdentifier(...));
0031 ///     if (it != container.end()) {
0032 ///         ...
0033 ///     }
0034 ///
0035 /// Trailing zero levels of stored geometry identifiers are used as
0036 /// broadcast values to refer to higher-level objects within the geometry, e.g.
0037 /// a geometry identifier with vanishing approach and sensitive index identifies
0038 /// a layer. An entry will all geometry identifier levels set to zero acts as
0039 /// the global default value.
0040 ///
0041 /// The container also supports range-based iteration over all stored elements
0042 ///
0043 ///     for (const auto& element : container) {
0044 ///         ...
0045 ///     }
0046 ///
0047 /// and index-based access to stored elements and associated geometry
0048 /// identifiers
0049 ///
0050 ///     GeometryIdentifier id3 = container.idAt(3);
0051 ///     const auto& element4 = container.valueAt(4);
0052 ///
0053 /// @note No guarantees are given for the element order when using range-based
0054 ///   or index-based access. Any apparent ordering must be considered an
0055 ///   implementation detail and might change.
0056 ///
0057 /// Adding elements is potentially expensive as the internal lookup structure
0058 /// must be updated. In addition, modifying an element in-place could change its
0059 /// identifier which would also break the lookup. Thus, the container can not be
0060 /// modified after construction to prevent misuse.
0061 template <typename value_t>
0062 class GeometryHierarchyMap {
0063  public:
0064   /// Combined geometry identifier and value element. Only used for input.
0065   using InputElement = typename std::pair<GeometryIdentifier, value_t>;
0066   using Iterator = typename std::vector<value_t>::const_iterator;
0067   using Size = typename std::vector<value_t>::size_type;
0068   using Value = value_t;
0069 
0070   /// Construct the container from the given elements.
0071   ///
0072   /// @param elements input elements (must be unique with respect to identifier)
0073   GeometryHierarchyMap(std::vector<InputElement> elements);
0074   /// Construct the container from an initializer list.
0075   ///
0076   /// @param elements input initializer list
0077   GeometryHierarchyMap(std::initializer_list<InputElement> elements);
0078 
0079   // defaulted constructors and assignment operators
0080   GeometryHierarchyMap() = default;
0081   GeometryHierarchyMap(const GeometryHierarchyMap&) = default;
0082   GeometryHierarchyMap(GeometryHierarchyMap&&) = default;
0083   ~GeometryHierarchyMap() = default;
0084   GeometryHierarchyMap& operator=(const GeometryHierarchyMap&) = default;
0085   GeometryHierarchyMap& operator=(GeometryHierarchyMap&&) = default;
0086 
0087   /// Return an iterator pointing to the beginning of the stored values.
0088   Iterator begin() const { return m_values.begin(); }
0089   /// Return an iterator pointing to the end of the stored values.
0090   Iterator end() const { return m_values.end(); }
0091   /// Check if any elements are stored.
0092   bool empty() const { return m_values.empty(); }
0093   /// Return the number of stored elements.
0094   Size size() const { return m_values.size(); }
0095 
0096   /// Access the geometry identifier for the i-th element with bounds check.
0097   ///
0098   /// @throws std::out_of_range for invalid indices
0099   GeometryIdentifier idAt(Size index) const { return m_ids.at(index); }
0100   /// Access the value of the i-th element in the container with bounds check.
0101   ///
0102   /// @throws std::out_of_range for invalid indices
0103   const Value& valueAt(Size index) const { return m_values.at(index); }
0104 
0105   /// Find the most specific value for a given geometry identifier.
0106   ///
0107   /// This can be either from the element matching exactly to the given geometry
0108   /// id, if it exists, or from the element for the next available higher level
0109   /// within the geometry hierarchy.
0110   ///
0111   /// @param id geometry identifier for which information is requested
0112   /// @retval iterator to an existing value
0113   /// @retval `.end()` iterator if no matching element exists
0114   Iterator find(GeometryIdentifier id) const;
0115 
0116  private:
0117   // NOTE this class assumes that it knows the ordering of the levels within
0118   //      the geometry id. if the geometry id changes, this code has to be
0119   //      adapted too. the asserts ensure that such a change is caught.
0120   static_assert(GeometryIdentifier().setVolume(1).value() <
0121                     GeometryIdentifier().setVolume(1).setBoundary(1).value(),
0122                 "Incompatible GeometryIdentifier hierarchy");
0123   static_assert(GeometryIdentifier().setBoundary(1).value() <
0124                     GeometryIdentifier().setBoundary(1).setLayer(1).value(),
0125                 "Incompatible GeometryIdentifier hierarchy");
0126   static_assert(GeometryIdentifier().setLayer(1).value() <
0127                     GeometryIdentifier().setLayer(1).setApproach(1).value(),
0128                 "Incompatible GeometryIdentifier hierarchy");
0129   static_assert(GeometryIdentifier().setApproach(1).value() <
0130                     GeometryIdentifier().setApproach(1).setSensitive(1).value(),
0131                 "Incompatible GeometryIdentifier hierarchy");
0132 
0133   using Identifier = GeometryIdentifier::Value;
0134 
0135   // encoded ids for all elements for faster lookup.
0136   std::vector<Identifier> m_ids;
0137   // validity bit masks for the ids: which parts to use for comparison
0138   std::vector<Identifier> m_masks;
0139   std::vector<Value> m_values;
0140 
0141   /// Construct a mask where all leading non-zero levels are set.
0142   static constexpr Identifier makeLeadingLevelsMask(GeometryIdentifier id) {
0143     // construct id from encoded value with all bits set
0144     auto allSet = GeometryIdentifier(~GeometryIdentifier::Value{0u});
0145     // manually iterate over identifier levels starting from the lowest
0146     if (id.sensitive() != 0u) {
0147       // all levels are valid; keep all bits set.
0148       return allSet.setExtra(0u).value();
0149     }
0150     if (id.approach() != 0u) {
0151       return allSet.setExtra(0u).setSensitive(0u).value();
0152     }
0153     if (id.layer() != 0u) {
0154       return allSet.setExtra(0u).setSensitive(0u).setApproach(0u).value();
0155     }
0156     if (id.boundary() != 0u) {
0157       return allSet.setExtra(0u)
0158           .setSensitive(0u)
0159           .setApproach(0u)
0160           .setLayer(0u)
0161           .value();
0162     }
0163     if (id.volume() != 0u) {
0164       return allSet.setExtra(0u)
0165           .setSensitive(0u)
0166           .setApproach(0u)
0167           .setLayer(0u)
0168           .setBoundary(0u)
0169           .value();
0170     }
0171     // no valid levels; all bits are zero.
0172     return Identifier{0u};
0173   }
0174   /// Construct a mask where only the highest level is set.
0175   static constexpr Identifier makeHighestLevelMask() {
0176     return makeLeadingLevelsMask(GeometryIdentifier(0u).setVolume(1u));
0177   }
0178   /// Compare the two identifiers only within the masked bits.
0179   static constexpr bool equalWithinMask(Identifier lhs, Identifier rhs,
0180                                         Identifier mask) {
0181     return (lhs & mask) == (rhs & mask);
0182   }
0183   /// Ensure identifier ordering and uniqueness.
0184   template <typename iterator_t>
0185   static void sortAndCheckDuplicates(iterator_t beg, iterator_t end);
0186 
0187   /// Fill the container from the input elements.
0188   ///
0189   /// This assumes that the elements are ordered and unique with respect to
0190   /// their identifiers.
0191   template <typename iterator_t>
0192   void fill(iterator_t beg, iterator_t end);
0193 };
0194 
0195 // implementations
0196 
0197 template <typename value_t>
0198 inline GeometryHierarchyMap<value_t>::GeometryHierarchyMap(
0199     std::vector<InputElement> elements) {
0200   sortAndCheckDuplicates(elements.begin(), elements.end());
0201   fill(elements.begin(), elements.end());
0202 }
0203 
0204 template <typename value_t>
0205 inline GeometryHierarchyMap<value_t>::GeometryHierarchyMap(
0206     std::initializer_list<InputElement> elements)
0207     : GeometryHierarchyMap(
0208           std::vector<InputElement>(elements.begin(), elements.end())) {}
0209 
0210 template <typename value_t>
0211 template <typename iterator_t>
0212 inline void GeometryHierarchyMap<value_t>::sortAndCheckDuplicates(
0213     iterator_t beg, iterator_t end) {
0214   // ensure elements are sorted by identifier
0215   std::sort(beg, end, [=](const auto& lhs, const auto& rhs) {
0216     return lhs.first < rhs.first;
0217   });
0218   // check that all elements have unique identifier
0219   auto dup = std::adjacent_find(beg, end, [](const auto& lhs, const auto& rhs) {
0220     return lhs.first == rhs.first;
0221   });
0222   if (dup != end) {
0223     throw std::invalid_argument("Input elements contain duplicates");
0224   }
0225 }
0226 
0227 template <typename value_t>
0228 template <typename iterator_t>
0229 inline void GeometryHierarchyMap<value_t>::fill(iterator_t beg,
0230                                                 iterator_t end) {
0231   const auto n = std::distance(beg, end);
0232   m_ids.clear();
0233   m_ids.reserve(n);
0234   m_masks.clear();
0235   m_masks.reserve(n);
0236   m_values.clear();
0237   m_values.reserve(n);
0238   for (; beg != end; ++beg) {
0239     m_ids.push_back(beg->first.value());
0240     m_masks.push_back(makeLeadingLevelsMask(beg->first.value()));
0241     m_values.push_back(std::move(beg->second));
0242   }
0243 }
0244 
0245 template <typename value_t>
0246 inline auto GeometryHierarchyMap<value_t>::find(GeometryIdentifier id) const
0247     -> Iterator {
0248   assert((m_ids.size() == m_values.size()) &&
0249          "Inconsistent container state: #ids != # values");
0250   assert((m_masks.size() == m_values.size()) &&
0251          "Inconsistent container state: #masks != #values");
0252 
0253   // we can not search for the element directly since the relevant one
0254   // might be stored at a higher level. ids for higher levels would always
0255   // be sorted before the requested id. searching for the first element
0256   // after the requested ensures that we include the full hierarchy.
0257   const auto it = std::upper_bound(m_ids.begin(), m_ids.end(), id.value());
0258   auto i = std::distance(m_ids.begin(), it);
0259 
0260   // now go up the hierarchy to find the first matching element.
0261   // example: the container stores four identifiers
0262   //
0263   //     2|x|x (volume-only)
0264   //     2|2|1 (volume, layer, and sensitive)
0265   //     2|3|x (volume and layer)
0266   //     2|3|4 (volume, layer, and sensitive)
0267   //
0268   // where | marks level boundaries. searching for either 2|3|4, 2|3|7, or
0269   // 2|4|x would first point to 2|3|4 and thus needs to go up the hierarchy.
0270   while (0 < i) {
0271     // index always starts after item of interest due to upper bound search
0272     --i;
0273 
0274     // if the input id does not even match at the highest hierarchy level
0275     // with the current comparison id, then have reached the end of this
0276     // hierarchy. having a special check for the highest level avoids an
0277     // unbounded search window all the way to the beginning of the container for
0278     // the global default entry.
0279     if (!equalWithinMask(id.value(), m_ids[i], makeHighestLevelMask())) {
0280       // check if a global default entry exists
0281       if (m_ids.front() == Identifier{0u}) {
0282         return begin();
0283       } else {
0284         return end();
0285       }
0286     }
0287 
0288     // since the search is going backwards in the sorted container, it
0289     // progresses from more specific to less specific elements. the first
0290     // match is automatically the appropriate one.
0291     if (equalWithinMask(id.value(), m_ids[i], m_masks[i])) {
0292       return std::next(begin(), i);
0293     }
0294   }
0295 
0296   // all options are exhausted and no matching element was found.
0297   return end();
0298 }
0299 
0300 }  // namespace Acts