GeometryContainers.h

// This file is part of the Acts project.
//
// Copyright (C) 2017-2020 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
#pragma once
 
#include <algorithm>
#include <boost/container/flat_set.hpp>
#include <cassert>
#include <utility>
 
#include "Acts/Geometry/GeometryIdentifier.hpp"
#include "Tracking/Sim/GroupBy.h"
#include "Tracking/Sim/Range.h"
 
namespace acts_examples {
namespace detail {
 
// extract the geometry identifier from a variety of types
struct GeometryIdGetter {
  // explicit geometry identifier are just forwarded
  constexpr Acts::GeometryIdentifier operator()(
      Acts::GeometryIdentifier geometryId) const {
    return geometryId;
  }
  // encoded geometry ids are converted back to geometry identifiers.
  constexpr Acts::GeometryIdentifier operator()(
      Acts::GeometryIdentifier::Value encoded) const {
    return Acts::GeometryIdentifier(encoded);
  }
  // support elements in map-like structures.
  template <typename T>
  constexpr Acts::GeometryIdentifier operator()(
      const std::pair<Acts::GeometryIdentifier, T>& mapItem) const {
    return mapItem.first;
  }
  // support elements that implement `.geometryId()`.
  template <typename T>
  inline auto operator()(const T& thing) const
      -> decltype(thing.geometryId(), Acts::GeometryIdentifier()) {
    return thing.geometryId();
  }
  // support reference_wrappers around such types as well
  template <typename T>
  inline auto operator()(std::reference_wrapper<T> thing) const
      -> decltype(thing.get().geometryId(), Acts::GeometryIdentifier()) {
    return thing.get().geometryId();
  }
};
 
struct CompareGeometryId {
  // indicate that comparisons between keys and full objects are allowed.
  using is_transparent = void;
  // compare two elements using the automatic key extraction.
  template <typename Left, typename Right>
  constexpr bool operator()(Left&& lhs, Right&& rhs) const {
    return GeometryIdGetter()(lhs) < GeometryIdGetter()(rhs);
  }
};
 
}  // namespace detail
 
template <typename T>
using GeometryIdMultiset =
    boost::container::flat_multiset<T, detail::CompareGeometryId>;
 
template <typename T>
using GeometryIdMultimap =
    GeometryIdMultiset<std::pair<Acts::GeometryIdentifier, T>>;
 
template <typename T>
inline acts_examples::Range<typename GeometryIdMultiset<T>::const_iterator>
selectVolume(const GeometryIdMultiset<T>& container,
             Acts::GeometryIdentifier::Value volume) {
  auto cmp = Acts::GeometryIdentifier().withVolume(volume);
  auto beg = std::lower_bound(container.begin(), container.end(), cmp,
                              detail::CompareGeometryId{});
  // WARNING overflows to volume==0 if the input volume is the last one
  cmp = Acts::GeometryIdentifier().withVolume(volume + 1u);
  // optimize search by using the lower bound as start point. also handles
  // volume overflows since the geo id would be located before the start of
  // the upper edge search window.
  auto end =
      std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
  return makeRange(beg, end);
}
 
template <typename T>
inline auto selectVolume(const GeometryIdMultiset<T>& container,
                         Acts::GeometryIdentifier id) {
  return selectVolume(container, id.volume());
}
 
template <typename T>
inline acts_examples::Range<typename GeometryIdMultiset<T>::const_iterator>
selectLayer(const GeometryIdMultiset<T>& container,
            Acts::GeometryIdentifier::Value volume,
            Acts::GeometryIdentifier::Value layer) {
  auto cmp = Acts::GeometryIdentifier().withVolume(volume).withLayer(layer);
  auto beg = std::lower_bound(container.begin(), container.end(), cmp,
                              detail::CompareGeometryId{});
  // WARNING resets to layer==0 if the input layer is the last one
  cmp = Acts::GeometryIdentifier().withVolume(volume).withLayer(layer + 1u);
  // optimize search by using the lower bound as start point. also handles
  // volume overflows since the geo id would be located before the start of
  // the upper edge search window.
  auto end =
      std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
  return makeRange(beg, end);
}
 
// Select all elements within the given layer.
template <typename T>
inline auto selectLayer(const GeometryIdMultiset<T>& container,
                        Acts::GeometryIdentifier id) {
  return selectLayer(container, id.volume(), id.layer());
}
 
template <typename T>
inline acts_examples::Range<typename GeometryIdMultiset<T>::const_iterator>
selectModule(const GeometryIdMultiset<T>& container,
             Acts::GeometryIdentifier geoId) {
  // module is the lowest level and defines a single geometry id value
  return makeRange(container.equal_range(geoId));
}
 
template <typename T>
inline auto selectModule(const GeometryIdMultiset<T>& container,
                         Acts::GeometryIdentifier::Value volume,
                         Acts::GeometryIdentifier::Value layer,
                         Acts::GeometryIdentifier::Value module_) {
  return selectModule(container, Acts::GeometryIdentifier()
                                     .withVolume(volume)
                                     .withLayer(layer)
                                     .withSensitive(module_));
}
 
template <typename T>
inline acts_examples::Range<typename GeometryIdMultiset<T>::const_iterator>
selectLowestNonZeroGeometryObject(const GeometryIdMultiset<T>& container,
                                  Acts::GeometryIdentifier geoId) {
  assert((geoId.boundary() == 0u) and "Boundary component must be zero");
  assert((geoId.approach() == 0u) and "Approach component must be zero");
 
  if (geoId.sensitive() != 0u) {
    return selectModule(container, geoId);
  } else if (geoId.layer() != 0u) {
    return selectLayer(container, geoId);
  } else if (geoId.volume() != 0u) {
    return selectVolume(container, geoId);
  } else {
    return makeRange(container.begin(), container.end());
  }
}
 
template <typename T>
inline acts_examples::GroupBy<typename GeometryIdMultiset<T>::const_iterator,
                              detail::GeometryIdGetter>
groupByModule(const GeometryIdMultiset<T>& container) {
  return makeGroupBy(container, detail::GeometryIdGetter());
}
 
template <typename T>
struct GeometryIdMultisetAccessor {
  using Container = GeometryIdMultiset<T>;
  using Key = Acts::GeometryIdentifier;
  using Value = typename GeometryIdMultiset<T>::value_type;
  using Iterator = typename GeometryIdMultiset<T>::const_iterator;
 
  // pointer to the container
  const Container* container_ = nullptr;
};
 
}  // namespace acts_examples