ActsExamples/EventData/GeometryContainers.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 "Acts/Geometry/GeometryIdentifier.hpp"
0012 #include "Acts/Utilities/PointerTraits.hpp"
0013 #include "ActsExamples/Utilities/GroupBy.hpp"
0014 #include "ActsExamples/Utilities/Range.hpp"
0015
0016 #include <algorithm>
0017 #include <cassert>
0018 #include <utility>
0019
0020 #include <boost/bimap.hpp>
0021 #include <boost/container/flat_map.hpp>
0022 #include <boost/container/flat_set.hpp>
0023
0024 namespace ActsExamples {
0025 namespace detail {
0026 /// @brief Concept to define objects that have a geometryId getter method
0027 template <typename ObjType>
0028 concept GeometryIdObj = requires(const ObjType& obj) {
0029   { obj.geometryId() } -> std::same_as<Acts::GeometryIdentifier>;
0030 };
0031 /// @brief Concept to define objects to pointer that have a geometryId getter method
0032 template <typename ObjType>
0033 concept GeometryIdPtrObj =
0034     Acts::PointerConcept<ObjType> &&
0035     GeometryIdObj<typename Acts::RemovePointer<ObjType>::type>;
0036
0037 /// @brief Concept to define an object which carries time
0038 template <typename ObjType>
0039 concept TimedObj = requires(const ObjType& obj) {
0040   { obj.time() };
0041 };
0042 /// @brief Concept to define a pointer to an object which carries time
0043 template <typename ObjType>
0044 concept TimedObjPtr = Acts::PointerConcept<ObjType> &&
0045                       TimedObj<typename Acts::RemovePointer<ObjType>::type>;
0046
0047 // extract the geometry identifier from a variety of types
0048 struct GeometryIdGetter {
0049   // explicit geometry identifier are just forwarded
0050   constexpr Acts::GeometryIdentifier operator()(
0051       Acts::GeometryIdentifier geometryId) const {
0052     return geometryId;
0053   }
0054   // encoded geometry ids are converted back to geometry identifiers.
0055   constexpr Acts::GeometryIdentifier operator()(
0056       Acts::GeometryIdentifier::Value encoded) const {
0057     return Acts::GeometryIdentifier(encoded);
0058   }
0059   // support elements in map-like structures.
0060   template <typename T>
0061   constexpr Acts::GeometryIdentifier operator()(
0062       const std::pair<Acts::GeometryIdentifier, T>& mapItem) const {
0063     return mapItem.first;
0064   }
0065   // Support pointer to object that implement `.geometryId()`.
0066   template <GeometryIdPtrObj T>
0067   constexpr Acts::GeometryIdentifier operator()(const T& thing) const {
0068     return thing->geometryId();
0069   }
0070   // support elements that implement `.geometryId()`.
0071   template <GeometryIdObj T>
0072   inline auto operator()(const T& thing) const {
0073     return thing.geometryId();
0074   }
0075   // support reference_wrappers around such types as well
0076   template <GeometryIdObj T>
0077   inline auto operator()(std::reference_wrapper<T> thing) const
0078       -> decltype(thing.get().geometryId(), Acts::GeometryIdentifier()) {
0079     return thing.get().geometryId();
0080   }
0081 };
0082
0083 struct CompareGeometryId {
0084   // indicate that comparisons between keys and full objects are allowed.
0085   using is_transparent = void;
0086
0087   /// @brief Comparator that sorts objects with the same geometry id by time
0088   template <TimedObj Left, TimedObj Right>
0089   constexpr bool operator()(const Left& lhs, const Right& rhs) const {
0090     const auto lId = GeometryIdGetter{}(lhs);
0091     const auto rId = GeometryIdGetter{}(rhs);
0092     if (lId == rId) {
0093       return lhs.time() < rhs.time();
0094     }
0095     return lId < rId;
0096   }
0097   template <TimedObjPtr Left, TimedObjPtr Right>
0098   constexpr bool operator()(const Left& lhs, const Right& rhs) const {
0099     return (*this)(*lhs, *rhs);
0100   }
0101   // compare two elements using the automatic key extraction.
0102   template <typename Left, typename Right>
0103   constexpr bool operator()(Left&& lhs, Right&& rhs) const {
0104     return GeometryIdGetter()(lhs) < GeometryIdGetter()(rhs);
0105   }
0106 };
0107
0108 }  // namespace detail
0109
0110 /// Store elements that know their detector geometry id, e.g. simulation hits.
0111 ///
0112 /// @tparam T type to be stored, must be compatible with `CompareGeometryId`
0113 ///
0114 /// The container stores an arbitrary number of elements for any geometry
0115 /// id. Elements can be retrieved via the geometry id; elements can be selected
0116 /// for a specific geometry id or for a larger range, e.g. a volume or a layer
0117 /// within the geometry hierarchy using the helper functions below. Elements can
0118 /// also be accessed by index that uniquely identifies each element regardless
0119 /// of geometry id.
0120 template <typename T>
0121 using GeometryIdMultiset =
0122     boost::container::flat_multiset<T, detail::CompareGeometryId>;
0123 /// Store elements indexed by an geometry id.
0124 ///
0125 /// @tparam T type to be stored
0126 ///
0127 /// The behaviour is the same as for the `GeometryIdMultiset` except that the
0128 /// stored elements do not know their geometry id themself. When iterating
0129 /// the iterator elements behave as for the `std::map`, i.e.
0130 ///
0131 ///     for (const auto& entry: elements) {
0132 ///         auto id = entry.first; // geometry id
0133 ///         const auto& el = entry.second; // stored element
0134 ///     }
0135 ///
0136 template <typename T>
0137 using GeometryIdMultimap =
0138     GeometryIdMultiset<std::pair<Acts::GeometryIdentifier, T>>;
0139
0140 /// Select all elements within the given volume.
0141 template <typename T>
0142 inline Range<typename GeometryIdMultiset<T>::const_iterator> selectVolume(
0143     const GeometryIdMultiset<T>& container,
0144     Acts::GeometryIdentifier::Value volume) {
0145   auto cmp = Acts::GeometryIdentifier().withVolume(volume);
0146   auto beg = std::lower_bound(container.begin(), container.end(), cmp,
0147                               detail::CompareGeometryId{});
0148   // WARNING overflows to volume==0 if the input volume is the last one
0149   cmp = Acts::GeometryIdentifier().withVolume(volume + 1u);
0150   // optimize search by using the lower bound as start point. also handles
0151   // volume overflows since the geo id would be located before the start of
0152   // the upper edge search window.
0153   auto end =
0154       std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
0155   return makeRange(beg, end);
0156 }
0157
0158 /// Select all elements within the given volume.
0159 template <typename T>
0160 inline auto selectVolume(const GeometryIdMultiset<T>& container,
0161                          Acts::GeometryIdentifier id) {
0162   return selectVolume(container, id.volume());
0163 }
0164
0165 /// Select all elements within the given layer.
0166 template <typename T>
0167 inline Range<typename GeometryIdMultiset<T>::const_iterator> selectLayer(
0168     const GeometryIdMultiset<T>& container,
0169     Acts::GeometryIdentifier::Value volume,
0170     Acts::GeometryIdentifier::Value layer) {
0171   auto cmp = Acts::GeometryIdentifier().withVolume(volume).withLayer(layer);
0172   auto beg = std::lower_bound(container.begin(), container.end(), cmp,
0173                               detail::CompareGeometryId{});
0174   // WARNING resets to layer==0 if the input layer is the last one
0175   cmp = Acts::GeometryIdentifier().withVolume(volume).withLayer(layer + 1u);
0176   // optimize search by using the lower bound as start point. also handles
0177   // volume overflows since the geo id would be located before the start of
0178   // the upper edge search window.
0179   auto end =
0180       std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
0181   return makeRange(beg, end);
0182 }
0183
0184 // Select all elements within the given layer.
0185 template <typename T>
0186 inline auto selectLayer(const GeometryIdMultiset<T>& container,
0187                         Acts::GeometryIdentifier id) {
0188   return selectLayer(container, id.volume(), id.layer());
0189 }
0190
0191 /// Select all elements for the given module / sensitive surface.
0192 template <typename T>
0193 inline Range<typename GeometryIdMultiset<T>::const_iterator> selectModule(
0194     const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier geoId) {
0195   // module is the lowest level and defines a single geometry id value
0196   return makeRange(container.equal_range(geoId));
0197 }
0198
0199 /// Select all elements for the given module / sensitive surface.
0200 template <typename T>
0201 inline auto selectModule(const GeometryIdMultiset<T>& container,
0202                          Acts::GeometryIdentifier::Value volume,
0203                          Acts::GeometryIdentifier::Value layer,
0204                          Acts::GeometryIdentifier::Value sensitive) {
0205   return selectModule(container, Acts::GeometryIdentifier()
0206                                      .withVolume(volume)
0207                                      .withLayer(layer)
0208                                      .withSensitive(sensitive));
0209 }
0210
0211 /// Select all elements for the lowest non-zero identifier component.
0212 ///
0213 /// Zero values of lower components are interpreted as wildcard search patterns
0214 /// that select all element at the given geometry hierarchy and below. This only
0215 /// applies to the lower components and not to intermediate zeros.
0216 ///
0217 /// Examples:
0218 /// - volume=2,layer=0,sensitive=3 -> select all elements in the sensitive
0219 /// - volume=1,layer=2,sensitive=0 -> select all elements in the layer
0220 /// - volume=3,layer=0,sensitive=0 -> select all elements in the volume
0221 ///
0222 /// @note An identifier with all components set to zero selects the whole input
0223 ///   container.
0224 /// @note Boundary and approach surfaces do not really fit into the geometry
0225 ///   hierarchy and must be set to zero for the selection. If they are set on an
0226 ///   input identifier, the behaviour of this search method is undefined.
0227 template <typename T>
0228 inline Range<typename GeometryIdMultiset<T>::const_iterator>
0229 selectLowestNonZeroGeometryObject(const GeometryIdMultiset<T>& container,
0230                                   Acts::GeometryIdentifier geoId) {
0231   assert((geoId.boundary() == 0u) && "Boundary component must be zero");
0232   assert((geoId.approach() == 0u) && "Approach component must be zero");
0233
0234   if (geoId.sensitive() != 0u) {
0235     return selectModule(container, geoId);
0236   } else if (geoId.layer() != 0u) {
0237     return selectLayer(container, geoId);
0238   } else if (geoId.volume() != 0u) {
0239     return selectVolume(container, geoId);
0240   } else {
0241     return makeRange(container.begin(), container.end());
0242   }
0243 }
0244
0245 /// Iterate over groups of elements belonging to each module/ sensitive surface.
0246 template <typename T>
0247 inline GroupBy<typename GeometryIdMultiset<T>::const_iterator,
0248                detail::GeometryIdGetter>
0249 groupByModule(const GeometryIdMultiset<T>& container) {
0250   return makeGroupBy(container, detail::GeometryIdGetter());
0251 }
0252
0253 /// The accessor for the GeometryIdMultiset container
0254 ///
0255 /// It wraps up a few lookup methods to be used in the Combinatorial Kalman
0256 /// Filter
0257 template <typename T>
0258 struct GeometryIdMultisetAccessor {
0259   using Container = GeometryIdMultiset<T>;
0260   using Key = Acts::GeometryIdentifier;
0261   using Value = typename GeometryIdMultiset<T>::value_type;
0262   using Iterator = typename GeometryIdMultiset<T>::const_iterator;
0263
0264   // pointer to the container
0265   const Container* container = nullptr;
0266 };
0267
0268 /// A map that allows mapping back and forth between ACTS and Athena Geometry
0269 /// Ids
0270 using GeometryIdMapActsAthena =
0271     boost::bimap<std::uint64_t, Acts::GeometryIdentifier>;
0272
0273 }  // namespace ActsExamples