JugReco/JugReco/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 "JugReco/Utilities/GroupBy.hpp"
0012 #include "JugReco/Utilities/Range.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014
0015 #include <algorithm>
0016 #include <cstddef>
0017 #include <utility>
0018
0019 #include <boost/container/flat_map.hpp>
0020 #include <boost/container/flat_set.hpp>
0021
0022 namespace Jug {
0023 namespace detail {
0024 // extract the geometry identifier from a variety of types
0025 struct GeometryIdGetter {
0026   // explicit geometry identifier are just forwarded
0027   constexpr Acts::GeometryIdentifier operator()(Acts::GeometryIdentifier geometryId) const {
0028     return geometryId;
0029   }
0030   // encoded geometry ids are converted back to geometry identifiers.
0031   constexpr Acts::GeometryIdentifier operator()(Acts::GeometryIdentifier::Value encoded) const {
0032     return Acts::GeometryIdentifier(encoded);
0033   }
0034   // support elements in map-like structures.
0035   template <typename T>
0036   constexpr Acts::GeometryIdentifier operator()(
0037       const std::pair<Acts::GeometryIdentifier, T>& mapItem) const {
0038     return mapItem.first;
0039   }
0040   // support elements that implement `.geometryId()`.
0041   template <typename T>
0042   inline auto operator()(const T& thing) const
0043       -> decltype(thing.geometryId(), Acts::GeometryIdentifier()) {
0044     return thing.geometryId();
0045   }
0046 };
0047
0048 struct CompareGeometryId {
0049   // indicate that comparisons between keys and full objects are allowed.
0050   using is_transparent = void;
0051   // compare two elements using the automatic key extraction.
0052   template <typename Left, typename Right>
0053   constexpr bool operator()(Left&& lhs, Right&& rhs) const {
0054     return GeometryIdGetter()(lhs) < GeometryIdGetter()(rhs);
0055   }
0056 };
0057 }  // namespace detail
0058
0059 /// Store elements that know their detector geometry id, e.g. simulation hits.
0060 ///
0061 /// @tparam T type to be stored, must be compatible with `CompareGeometryId`
0062 ///
0063 /// The container stores an arbitrary number of elements for any geometry
0064 /// id. Elements can be retrieved via the geometry id; elements can be selected
0065 /// for a specific geometry id or for a larger range, e.g. a volume or a layer
0066 /// within the geometry hierachy using the helper functions below. Elements can
0067 /// also be accessed by index that uniquely identifies each element regardless
0068 /// of geometry id.
0069 template <typename T>
0070 using GeometryIdMultiset =
0071     boost::container::flat_multiset<T, detail::CompareGeometryId>;
0072
0073 /// Store elements indexed by an geometry id.
0074 ///
0075 /// @tparam T type to be stored
0076 ///
0077 /// The behaviour is the same as for the `GeometryIdMultiset` except that the
0078 /// stored elements do not know their geometry id themself. When iterating
0079 /// the iterator elements behave as for the `std::map`, i.e.
0080 ///
0081 ///     for (const auto& entry: elements) {
0082 ///         auto id = entry.first; // geometry id
0083 ///         const auto& el = entry.second; // stored element
0084 ///     }
0085 ///
0086 template <typename T>
0087 using GeometryIdMultimap = GeometryIdMultiset<std::pair<Acts::GeometryIdentifier, T>>;
0088
0089 /// Select all elements within the given volume.
0090 template <typename T>
0091 inline Range<typename GeometryIdMultiset<T>::const_iterator> selectVolume(
0092     const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier::Value volume) {
0093   auto cmp = Acts::GeometryIdentifier().setVolume(volume);
0094   auto beg = std::lower_bound(container.begin(), container.end(), cmp,
0095                               detail::CompareGeometryId{});
0096   // WARNING overflows to volume==0 if the input volume is the last one
0097   cmp = Acts::GeometryIdentifier().setVolume(volume + 1u);
0098   // optimize search by using the lower bound as start point. also handles
0099   // volume overflows since the geo id would be located before the start of
0100   // the upper edge search window.
0101   auto end =
0102       std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
0103   return makeRange(beg, end);
0104 }
0105 template <typename T>
0106 inline auto selectVolume(const GeometryIdMultiset<T>& container,
0107                          Acts::GeometryIdentifier id) {
0108   return selectVolume(container, id.volume());
0109 }
0110
0111 /// Select all elements within the given layer.
0112 template <typename T>
0113 inline Range<typename GeometryIdMultiset<T>::const_iterator> selectLayer(
0114     const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier::Value volume,
0115     Acts::GeometryIdentifier::Value layer) {
0116   auto cmp = Acts::GeometryIdentifier().setVolume(volume).setLayer(layer);
0117   auto beg = std::lower_bound(container.begin(), container.end(), cmp,
0118                               detail::CompareGeometryId{});
0119   // WARNING resets to layer==0 if the input layer is the last one
0120   cmp = Acts::GeometryIdentifier().setVolume(volume).setLayer(layer + 1u);
0121   // optimize search by using the lower bound as start point. also handles
0122   // volume overflows since the geo id would be located before the start of
0123   // the upper edge search window.
0124   auto end =
0125       std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
0126   return makeRange(beg, end);
0127 }
0128 template <typename T>
0129 inline auto selectLayer(const GeometryIdMultiset<T>& container,
0130                         Acts::GeometryIdentifier id) {
0131   return selectLayer(container, id.volume(), id.layer());
0132 }
0133
0134 /// Select all elements for the given module / sensitive surface.
0135 template <typename T>
0136 inline Range<typename GeometryIdMultiset<T>::const_iterator> selectModule(
0137     const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier geoId) {
0138   // module is the lowest level and defines a single geometry id value
0139   return makeRange(container.equal_range(geoId));
0140 }
0141 template <typename T>
0142 inline auto selectModule(const GeometryIdMultiset<T>& container,
0143                          Acts::GeometryIdentifier::Value volume,
0144                          Acts::GeometryIdentifier::Value layer,
0145                          Acts::GeometryIdentifier::Value module) {
0146   return selectModule(
0147       container,
0148       Acts::GeometryIdentifier().setVolume(volume).setLayer(layer).setSensitive(
0149           module));
0150 }
0151
0152 /// Iterate over groups of elements belonging to each module/ sensitive surface.
0153 template <typename T>
0154 inline GroupBy<typename GeometryIdMultiset<T>::const_iterator,
0155                detail::GeometryIdGetter>
0156 groupByModule(const GeometryIdMultiset<T>& container) {
0157   return makeGroupBy(container, detail::GeometryIdGetter());
0158 }
0159
0160 }  // namespace FW