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