GeometryHierarchyMap.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GeometryHierarchyMap.hpp

// This file is part of the Acts project.
//
// Copyright (C) 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 "Acts/Geometry/GeometryIdentifier.hpp"

#include <algorithm>
#include <cassert>
#include <initializer_list>
#include <iterator>
#include <stdexcept>
#include <utility>
#include <vector>

namespace Acts {

template <typename value_t>
class GeometryHierarchyMap {
 public:
  using InputElement = typename std::pair<GeometryIdentifier, value_t>;
  using Iterator = typename std::vector<value_t>::const_iterator;
  using Size = typename std::vector<value_t>::size_type;
  using Value = value_t;

  GeometryHierarchyMap(std::vector<InputElement> elements);
  GeometryHierarchyMap(std::initializer_list<InputElement> elements);

  // defaulted constructors and assignment operators
  GeometryHierarchyMap() = default;
  GeometryHierarchyMap(const GeometryHierarchyMap&) = default;
  GeometryHierarchyMap(GeometryHierarchyMap&&) = default;
  ~GeometryHierarchyMap() = default;
  GeometryHierarchyMap& operator=(const GeometryHierarchyMap&) = default;
  GeometryHierarchyMap& operator=(GeometryHierarchyMap&&) = default;

  Iterator begin() const { return m_values.begin(); }
  Iterator end() const { return m_values.end(); }
  bool empty() const { return m_values.empty(); }
  Size size() const { return m_values.size(); }

  GeometryIdentifier idAt(Size index) const { return m_ids.at(index); }
  const Value& valueAt(Size index) const { return m_values.at(index); }

  Iterator find(GeometryIdentifier id) const;

 private:
  // NOTE this class assumes that it knows the ordering of the levels within
  //      the geometry id. if the geometry id changes, this code has to be
  //      adapted too. the asserts ensure that such a change is caught.
  static_assert(GeometryIdentifier().setVolume(1).value() <
                    GeometryIdentifier().setVolume(1).setBoundary(1).value(),
                "Incompatible GeometryIdentifier hierarchy");
  static_assert(GeometryIdentifier().setBoundary(1).value() <
                    GeometryIdentifier().setBoundary(1).setLayer(1).value(),
                "Incompatible GeometryIdentifier hierarchy");
  static_assert(GeometryIdentifier().setLayer(1).value() <
                    GeometryIdentifier().setLayer(1).setApproach(1).value(),
                "Incompatible GeometryIdentifier hierarchy");
  static_assert(GeometryIdentifier().setApproach(1).value() <
                    GeometryIdentifier().setApproach(1).setSensitive(1).value(),
                "Incompatible GeometryIdentifier hierarchy");

  using Identifier = GeometryIdentifier::Value;

  // encoded ids for all elements for faster lookup.
  std::vector<Identifier> m_ids;
  // validity bit masks for the ids: which parts to use for comparison
  std::vector<Identifier> m_masks;
  std::vector<Value> m_values;

  static constexpr Identifier makeLeadingLevelsMask(GeometryIdentifier id) {
    // construct id from encoded value with all bits set
    auto allSet = GeometryIdentifier(~GeometryIdentifier::Value(0u));
    // manually iterate over identifier levels starting from the lowest
    if (id.sensitive() != 0u) {
      // all levels are valid; keep all bits set.
      return allSet.value();
    }
    if (id.approach() != 0u) {
      return allSet.setSensitive(0u).value();
    }
    if (id.layer() != 0u) {
      return allSet.setSensitive(0u).setApproach(0u).value();
    }
    if (id.boundary() != 0u) {
      return allSet.setSensitive(0u).setApproach(0u).setLayer(0u).value();
    }
    if (id.volume() != 0u) {
      return allSet.setSensitive(0u)
          .setApproach(0u)
          .setLayer(0u)
          .setBoundary(0u)
          .value();
    }
    // no valid levels; all bits are zero.
    return Identifier(0u);
  }
  static constexpr Identifier makeHighestLevelMask() {
    return makeLeadingLevelsMask(GeometryIdentifier(0u).setVolume(1u));
  }
  static constexpr bool equalWithinMask(Identifier lhs, Identifier rhs,
                                        Identifier mask) {
    return (lhs & mask) == (rhs & mask);
  }
  template <typename iterator_t>
  static void sortAndCheckDuplicates(iterator_t beg, iterator_t end);

  template <typename iterator_t>
  void fill(iterator_t beg, iterator_t end);
};

// implementations

template <typename value_t>
inline GeometryHierarchyMap<value_t>::GeometryHierarchyMap(
    std::vector<InputElement> elements) {
  sortAndCheckDuplicates(elements.begin(), elements.end());
  fill(elements.begin(), elements.end());
}

template <typename value_t>
inline GeometryHierarchyMap<value_t>::GeometryHierarchyMap(
    std::initializer_list<InputElement> elements)
    : GeometryHierarchyMap(
          std::vector<InputElement>(elements.begin(), elements.end())) {}

template <typename value_t>
template <typename iterator_t>
inline void GeometryHierarchyMap<value_t>::sortAndCheckDuplicates(
    iterator_t beg, iterator_t end) {
  // ensure elements are sorted by identifier
  std::sort(beg, end, [=](const auto& lhs, const auto& rhs) {
    return lhs.first < rhs.first;
  });
  // check that all elements have unique identifier
  auto dup = std::adjacent_find(beg, end, [](const auto& lhs, const auto& rhs) {
    return lhs.first == rhs.first;
  });
  if (dup != end) {
    throw std::invalid_argument("Input elements contain duplicates");
  }
}

template <typename value_t>
template <typename iterator_t>
inline void GeometryHierarchyMap<value_t>::fill(iterator_t beg,
                                                iterator_t end) {
  const auto n = std::distance(beg, end);
  m_ids.clear();
  m_ids.reserve(n);
  m_masks.clear();
  m_masks.reserve(n);
  m_values.clear();
  m_values.reserve(n);
  for (; beg != end; ++beg) {
    m_ids.push_back(beg->first.value());
    m_masks.push_back(makeLeadingLevelsMask(beg->first.value()));
    m_values.push_back(std::move(beg->second));
  }
}

template <typename value_t>
inline auto GeometryHierarchyMap<value_t>::find(GeometryIdentifier id) const
    -> Iterator {
  assert((m_ids.size() == m_values.size()) and
         "Inconsistent container state: #ids != # values");
  assert((m_masks.size() == m_values.size()) and
         "Inconsistent container state: #masks != #values");

  // we can not search for the element directly since the relevant one
  // might be stored at a higher level. ids for higher levels would always
  // be sorted before the requested id. searching for the first element
  // after the requested ensures that we include the full hierarchy.
  const auto it = std::upper_bound(m_ids.begin(), m_ids.end(), id.value());
  auto i = std::distance(m_ids.begin(), it);

  // now go up the hierarchy to find the first matching element.
  // example: the container stores four identifiers
  //
  //     2|x|x (volume-only)
  //     2|2|1 (volume, layer, and sensitive)
  //     2|3|x (volume and layer)
  //     2|3|4 (volume, layer, and sensitive)
  //
  // where | marks level boundaries. searching for either 2|3|4, 2|3|7, or
  // 2|4|x would first point to 2|3|4 and thus needs to go up the hierarchy.
  while (0 < i) {
    // index always starts after item of interest due to upper bound search
    --i;

    // if the input id does not even match at the highest hierarchy level
    // with the current comparison id, then have reached the end of this
    // hierarchy. having a special check for the highest level avoids an
    // unbounded search window all the way to the beginning of the container for
    // the global default entry.
    if (not equalWithinMask(id.value(), m_ids[i], makeHighestLevelMask())) {
      // check if a global default entry exists
      if (m_ids.front() == Identifier(0u)) {
        return begin();
      } else {
        return end();
      }
    }

    // since the search is going backwards in the sorted container, it
    // progresses from more specific to less specific elements. the first
    // match is automatically the appropriate one.
    if (equalWithinMask(id.value(), m_ids[i], m_masks[i])) {
      return std::next(begin(), i);
    }
  }

  // all options are exhausted and no matching element was found.
  return end();
}

}  // namespace Acts