ActsExamples/EventData/GeometryContainers.hpp

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