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