SurfaceArray.hpp

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

// 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 "Acts/Definitions/Algebra.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Utilities/BinningType.hpp"
#include "Acts/Utilities/IAxis.hpp"
#include "Acts/Utilities/detail/Axis.hpp"
#include "Acts/Utilities/detail/Grid.hpp"

#include <iostream>
#include <type_traits>
#include <vector>

namespace Acts {

using SurfaceVector = std::vector<const Surface*>;

class SurfaceArray {
 public:
  struct ISurfaceGridLookup {

    virtual void fill(const GeometryContext& gctx,
                      const SurfaceVector& surfaces) = 0;


    virtual size_t completeBinning(const GeometryContext& gctx,
                                   const SurfaceVector& surfaces) = 0;

    virtual SurfaceVector& lookup(const Vector3& position) = 0;

    virtual const SurfaceVector& lookup(const Vector3& position) const = 0;

    virtual SurfaceVector& lookup(size_t bin) = 0;

    virtual const SurfaceVector& lookup(size_t bin) const = 0;

    virtual const SurfaceVector& neighbors(const Vector3& position) const = 0;

    virtual size_t size() const = 0;

    virtual Vector3 getBinCenter(size_t bin) const = 0;

    virtual std::vector<const IAxis*> getAxes() const = 0;

    virtual size_t dimensions() const = 0;

    virtual bool isValidBin(size_t bin) const = 0;

    virtual std::vector<BinningValue> binningValues() const { return {}; };

    virtual ~ISurfaceGridLookup() = 0;
  };

  template <class... Axes>
  struct SurfaceGridLookup : ISurfaceGridLookup {
    static constexpr size_t DIM = sizeof...(Axes);

   public:
    using point_t =
        std::conditional_t<DIM == 1, std::array<double, 1>, ActsVector<DIM>>;
    using Grid_t = detail::Grid<SurfaceVector, Axes...>;

    SurfaceGridLookup(std::function<point_t(const Vector3&)> globalToLocal,
                      std::function<Vector3(const point_t&)> localToGlobal,
                      std::tuple<Axes...> axes,
                      std::vector<BinningValue> bValues = {})
        : m_globalToLocal(std::move(globalToLocal)),
          m_localToGlobal(std::move(localToGlobal)),
          m_grid(std::move(axes)),
          m_binValues(std::move(bValues)) {
      m_neighborMap.resize(m_grid.size());
    }

    void fill(const GeometryContext& gctx,
              const SurfaceVector& surfaces) override {
      for (const auto& srf : surfaces) {
        Vector3 pos = srf->binningPosition(gctx, binR);
        lookup(pos).push_back(srf);
      }

      populateNeighborCache();
    }

    size_t completeBinning(const GeometryContext& gctx,
                           const SurfaceVector& surfaces) override {
      size_t binCompleted = 0;
      size_t nBins = size();
      double minPath = 0;
      double curPath = 0;
      const Surface* minSrf = nullptr;

      for (size_t b = 0; b < nBins; ++b) {
        if (!isValidBin(b)) {
          continue;
        }
        std::vector<const Surface*>& binContent = lookup(b);
        // only complete if we have an empty bin
        if (!binContent.empty()) {
          continue;
        }

        Vector3 binCtr = getBinCenter(b);
        minPath = std::numeric_limits<double>::max();
        for (const auto& srf : surfaces) {
          curPath = (binCtr - srf->binningPosition(gctx, binR)).norm();

          if (curPath < minPath) {
            minPath = curPath;
            minSrf = srf;
          }
        }

        binContent.push_back(minSrf);
        ++binCompleted;
      }

      // recreate neighborcache
      populateNeighborCache();
      return binCompleted;
    }

    SurfaceVector& lookup(const Vector3& position) override {
      return m_grid.atPosition(m_globalToLocal(position));
    }

    const SurfaceVector& lookup(const Vector3& position) const override {
      return m_grid.atPosition(m_globalToLocal(position));
    }

    SurfaceVector& lookup(size_t bin) override { return m_grid.at(bin); }

    const SurfaceVector& lookup(size_t bin) const override {
      return m_grid.at(bin);
    }

    const SurfaceVector& neighbors(const Vector3& position) const override {
      auto lposition = m_globalToLocal(position);
      return m_neighborMap.at(m_grid.globalBinFromPosition(lposition));
    }

    size_t size() const override { return m_grid.size(); }

    std::vector<BinningValue> binningValues() const override {
      return m_binValues;
    }

    Vector3 getBinCenter(size_t bin) const override {
      return getBinCenterImpl(bin);
    }

    std::vector<const IAxis*> getAxes() const override {
      auto arr = m_grid.axes();
      return std::vector<const IAxis*>(arr.begin(), arr.end());
    }

    size_t dimensions() const override { return DIM; }

    bool isValidBin(size_t bin) const override {
      std::array<size_t, DIM> indices = m_grid.localBinsFromGlobalBin(bin);
      std::array<size_t, DIM> nBins = m_grid.numLocalBins();
      for (size_t i = 0; i < indices.size(); ++i) {
        size_t idx = indices.at(i);
        if (idx <= 0 || idx >= nBins.at(i) + 1) {
          return false;
        }
      }

      return true;
    }

   private:
    void populateNeighborCache() {
      // calculate neighbors for every bin and store in map
      for (size_t i = 0; i < m_grid.size(); i++) {
        if (!isValidBin(i)) {
          continue;
        }
        typename Grid_t::index_t loc = m_grid.localBinsFromGlobalBin(i);
        auto neighborIdxs = m_grid.neighborHoodIndices(loc, 1u);
        std::vector<const Surface*>& neighbors = m_neighborMap.at(i);
        neighbors.clear();

        for (const auto idx : neighborIdxs) {
          const std::vector<const Surface*>& binContent = m_grid.at(idx);
          std::copy(binContent.begin(), binContent.end(),
                    std::back_inserter(neighbors));
        }
      }
    }

    template <size_t D = DIM, std::enable_if_t<D != 1, int> = 0>
    Vector3 getBinCenterImpl(size_t bin) const {
      return m_localToGlobal(ActsVector<DIM>(
          m_grid.binCenter(m_grid.localBinsFromGlobalBin(bin)).data()));
    }

    template <size_t D = DIM, std::enable_if_t<D == 1, int> = 0>
    Vector3 getBinCenterImpl(size_t bin) const {
      point_t pos = m_grid.binCenter(m_grid.localBinsFromGlobalBin(bin));
      return m_localToGlobal(pos);
    }

    std::function<point_t(const Vector3&)> m_globalToLocal;
    std::function<Vector3(const point_t&)> m_localToGlobal;
    Grid_t m_grid;
    std::vector<BinningValue> m_binValues;
    std::vector<SurfaceVector> m_neighborMap;
  };

  struct SingleElementLookup : ISurfaceGridLookup {
    SingleElementLookup(SurfaceVector::value_type element)
        : m_element({element}) {}

    SingleElementLookup(const SurfaceVector& elements) : m_element(elements) {}

    SurfaceVector& lookup(const Vector3& position) override {
      (void)position;
      return m_element;
    }

    const SurfaceVector& lookup(const Vector3& position) const override {
      (void)position;
      return m_element;
    }

    SurfaceVector& lookup(size_t bin) override {
      (void)bin;
      return m_element;
    }

    const SurfaceVector& lookup(size_t bin) const override {
      (void)bin;
      return m_element;
    }

    const SurfaceVector& neighbors(const Vector3& position) const override {
      (void)position;
      return m_element;
    }

    size_t size() const override { return 1; }

    Vector3 getBinCenter(size_t bin) const override {
      (void)bin;
      return Vector3(0, 0, 0);
    }

    std::vector<const IAxis*> getAxes() const override { return {}; }

    size_t dimensions() const override { return 0; }

    void fill(const GeometryContext& /*gctx*/,
              const SurfaceVector& /*surfaces*/) override {}

    size_t completeBinning(const GeometryContext& /*gctx*/,
                           const SurfaceVector& /*surfaces*/) override {
      return 0;
    }

    bool isValidBin(size_t bin) const override {
      (void)bin;
      return true;
    }

   private:
    SurfaceVector m_element;
  };

  SurfaceArray(std::unique_ptr<ISurfaceGridLookup> gridLookup,
               std::vector<std::shared_ptr<const Surface>> surfaces,
               const Transform3& transform = Transform3::Identity());

  SurfaceArray(std::shared_ptr<const Surface> srf);

  SurfaceVector& at(const Vector3& position) {
    return p_gridLookup->lookup(position);
  }

  const SurfaceVector& at(const Vector3& position) const {
    return p_gridLookup->lookup(position);
  }

  SurfaceVector& at(size_t bin) { return p_gridLookup->lookup(bin); }

  const SurfaceVector& at(size_t bin) const {
    return p_gridLookup->lookup(bin);
  }

  const SurfaceVector& neighbors(const Vector3& position) const {
    return p_gridLookup->neighbors(position);
  }

  size_t size() const { return p_gridLookup->size(); }

  Vector3 getBinCenter(size_t bin) { return p_gridLookup->getBinCenter(bin); }

  const SurfaceVector& surfaces() const { return m_surfacesRawPointers; }

  std::vector<const IAxis*> getAxes() const { return p_gridLookup->getAxes(); }

  bool isValidBin(size_t bin) const { return p_gridLookup->isValidBin(bin); }

  const Transform3& transform() const { return m_transform; }

  std::vector<BinningValue> binningValues() const {
    return p_gridLookup->binningValues();
  };

  std::ostream& toStream(const GeometryContext& gctx, std::ostream& sl) const;

 private:
  std::unique_ptr<ISurfaceGridLookup> p_gridLookup;
  // this vector makes sure we have shared ownership over the surfaces
  std::vector<std::shared_ptr<const Surface>> m_surfaces;
  // this vector is returned, so that (expensive) copying of the shared_ptr
  // vector does not happen by default
  SurfaceVector m_surfacesRawPointers;
  // this is only used to keep info on transform applied
  // by l2g and g2l
  Transform3 m_transform;
};

}  // namespace Acts