geo_vectors_utils.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_GEO_VECTORS_UTILS_H
#define LARCOREALG_GEOMETRY_GEO_VECTORS_UTILS_H

// LArSoft libraries
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"

// ROOT libraries
#include "Math/GenVector/DisplacementVector2D.h"
#include "Math/GenVector/DisplacementVector3D.h"
#include "Math/GenVector/LorentzVector.h"
#include "Math/GenVector/PositionVector2D.h"
#include "Math/GenVector/PositionVector3D.h"

// C/C++ standard library
#include <array>
#include <cassert>
#include <functional>  // std::mem_fn()
#include <iterator>    // std::back_inserter()
#include <type_traits> // std::declval(), std::is_same<>, ...
#include <vector>

namespace lar::util::simple_geo {
  template <typename T>
  struct Point2D;
  template <typename T>
  struct Point3D;
}

namespace geo {

  namespace vect {

    //--------------------------------------------------------------------------
    namespace details {
      //------------------------------------------------------------------------
      template <typename Op, typename... T>
      struct AccumulateImpl;

      template <typename Op, typename First, typename Second, typename... Others>
      struct AccumulateImpl<Op, First, Second, Others...> {
        static auto compute(Op op, First&& a, Second&& b, Others&&... others) -> decltype(auto)
        {
          return op(a,
                    AccumulateImpl<Op, Second, Others...>::compute(
                      op, std::forward<Second>(b), std::forward<Others>(others)...));
        }
      }; // AccumulateImpl<>

      template <typename Op, typename T>
      struct AccumulateImpl<Op, T> {
        static auto compute(Op, T&& v) -> decltype(auto) { return std::forward<T>(v); }
      };

      template <typename Op, typename... T>
      auto extended_accumulate(Op op, T&&... args)
      {
        return AccumulateImpl<Op, T...>::compute(op, std::forward<T>(args)...);
      }

      template <typename... T>
      auto extended_and(T... args) -> decltype(auto)
      {
        auto and_op = [](auto&& a, auto&& b) { return a && b; };
        return extended_accumulate(and_op, std::forward<T>(args)...);
      }

      //------------------------------------------------------------------------
      //
      // These two pages of metaprogramming madness are aimed to have objects
      // that can provide a uniform coordinate get/set for any type of vector,
      // with the minimal run-time overhead.
      // Complications arise from the fact that TVector3 setter returns void,
      // while GenVector vectors return a reference to the vector itself
      // (hence the need of the additional SetterResult template parameter,
      // and of quite some code to autodetect its value).
      // The makeXxxx() functions are aimed to enclose the additional
      // autodetection overlay, which should become unnecessary with C++17.
      //

      template <typename Vector>
      struct VectorScalar {
        using type = typename Vector::Scalar;
      };

      template <typename Vector>
      using VectorScalar_t = typename VectorScalar<Vector>::type;

      //------------------------------------------------------------------------
      template <typename Vector>
      struct HasGetter {
      private:
        template <typename Class>
        static constexpr bool TestX(decltype(std::declval<Class>().X())*)
        {
          return true;
        }
        template <typename Class>
        static constexpr bool TestX(...)
        {
          return false;
        }
        template <typename Class>
        static constexpr bool TestY(decltype(std::declval<Class>().Y())*)
        {
          return true;
        }
        template <typename Class>
        static constexpr bool TestY(...)
        {
          return false;
        }
        template <typename Class>
        static constexpr bool TestZ(decltype(std::declval<Class>().Z())*)
        {
          return true;
        }
        template <typename Class>
        static constexpr bool TestZ(...)
        {
          return false;
        }
        template <typename Class>
        static constexpr bool TestT(decltype(std::declval<Class>().T())*)
        {
          return true;
        }
        template <typename Class>
        static constexpr bool TestT(...)
        {
          return false;
        }

      public:
        static constexpr bool X = TestX<Vector>(nullptr);
        static constexpr bool Y = TestY<Vector>(nullptr);
        static constexpr bool Z = TestZ<Vector>(nullptr);
        static constexpr bool T = TestT<Vector>(nullptr);
      }; // struct HasGetter<>

      template <typename Vector>
      constexpr bool HasX()
      {
        return HasGetter<Vector>::X;
      }
      template <typename Vector>
      constexpr bool HasY()
      {
        return HasGetter<Vector>::Y;
      }
      template <typename Vector>
      constexpr bool HasZ()
      {
        return HasGetter<Vector>::Z;
      }
      template <typename Vector>
      constexpr bool HasT()
      {
        return HasGetter<Vector>::T;
      }

      template <typename Vector, typename = void>
      struct DimensionImpl;

      template <typename Vector>
      constexpr unsigned int dimension()
      {
        return DimensionImpl<Vector>::value;
      }

      template <typename Vector>
      using CoordinateArray_t = std::array<VectorScalar_t<Vector>, dimension<Vector>()>;

      template <typename T>
      struct MemberFuncReturnType {
        // with C++17 `result_type` is deprecated,
        // and a different way will be needed
        using type = typename decltype(std::mem_fn(std::declval<T>()))::result_type;
      }; // MemberFuncReturnType

      template <typename T>
      using MemberFuncReturn_t = typename MemberFuncReturnType<T>::type;

      template <typename T>
      struct MemberFuncClassType;

      template <typename Class, typename Func>
      struct MemberFuncClassType<Func Class::*> {
        using type = Class;
      };

      template <typename T>
      using MemberFuncClass_t = typename MemberFuncClassType<T>::type;

      template <typename Vector, typename SetterType = void>
      struct BaseCoordTypes;

      template <typename Vector>
      struct BaseCoordTypes<Vector, void> {
        using Vector_t = std::decay_t<Vector>;
        using Scalar_t = VectorScalar_t<Vector>;
        using Getter_t = Scalar_t (Vector_t::*)() const;
      }; // struct BaseCoordTypes<void>

      template <typename Vector, typename SetterType>
      struct BaseCoordTypes {
      private:
        using BaseTypes_t = BaseCoordTypes<Vector, void>;

      public:
        using Scalar_t = typename BaseTypes_t::Scalar_t;
        using Vector_t = typename BaseTypes_t::Vector_t;
        using Getter_t = typename BaseTypes_t::Getter_t;
        using Setter_t = SetterType;
        using SetterReturn_t = MemberFuncReturn_t<Setter_t>;
        static_assert(std::is_same<Setter_t, SetterReturn_t (Vector_t::*)(Scalar_t)>(),
                      "Invalid setter type");
      }; // struct BaseCoordTypes<>

      template <typename Vector>
      struct CoordGetterTraits {
      private:
        using BaseTypes_t = BaseCoordTypes<Vector>;

      public:
        using Vector_t = typename BaseTypes_t::Vector_t;
        using Scalar_t = typename BaseTypes_t::Scalar_t;
        using Getter_t = typename BaseTypes_t::Getter_t;
      }; // struct CoordGetterTraits

      template <typename Vector>
      class CoordGetter {
        using Traits_t = CoordGetterTraits<Vector>;

      public:
        using Vector_t = typename Traits_t::Vector_t;
        using Scalar_t = typename Traits_t::Scalar_t;
        using Getter_t = typename Traits_t::Getter_t;

        constexpr CoordGetter(Getter_t getter) : fGetter(getter) {}

        Scalar_t operator()(Vector_t const& v) const { return get(v); }

        Scalar_t get(Vector_t const& v) const { return (v.*fGetter)(); }

      private:
        Getter_t fGetter; 

      }; // class CoordGetter<>

      template <typename Getter>
      constexpr auto makeCoordReader(Getter getter)
      {
        using Vector_t = std::remove_reference_t<MemberFuncClass_t<Getter>>;
        return CoordGetter<Vector_t>{getter};
      }

      template <typename Vector, typename SetterType>
      struct CoordManagerTraits {
      private:
        using BaseTypes_t = BaseCoordTypes<Vector, SetterType>;

      public:
        using Vector_t = typename BaseTypes_t::Vector_t;
        using Scalar_t = typename BaseTypes_t::Scalar_t;
        using Getter_t = typename BaseTypes_t::Getter_t;
        using Setter_t = typename BaseTypes_t::Setter_t;
      }; // struct VectorCoordManagerTraits<>

      template <typename Vector, typename SetterType>
      class CoordManager : public CoordGetter<Vector> {
        using Base_t = CoordGetter<Vector>;
        using Traits_t = CoordManagerTraits<Vector, SetterType>;

      public:
        using Vector_t = typename Traits_t::Vector_t; // this is not constant
        using Scalar_t = typename Traits_t::Scalar_t;
        using Getter_t = typename Traits_t::Getter_t;
        using Setter_t = typename Traits_t::Setter_t;

        constexpr CoordManager(Getter_t getter, Setter_t setter) : Base_t(getter), fSetter(setter)
        {}

        void operator()(Vector_t& v, Scalar_t c) const { set(v, c); }

        void set(Vector_t& v, Scalar_t c) const { (v.*fSetter)(c); }

        void incr(Vector_t& v, Scalar_t c) const { set(v, Base_t::get(v) + c); }

        void decr(Vector_t& v, Scalar_t c) const { set(v, Base_t::get(v) - c); }

        void mult(Vector_t& v, Scalar_t f) const { set(v, Base_t::get(v) * f); }

        void div(Vector_t& v, Scalar_t f) const { set(v, Base_t::get(v) / f); }

      private:
        Setter_t fSetter; 
      };                  // class CoordManager<>

      template <typename Getter, typename Setter>
      constexpr auto makeCoordManager(Getter getter, Setter setter)
      {
        using Vector_t = std::remove_reference_t<MemberFuncClass_t<Getter>>;
        return CoordManager<Vector_t, Setter>{getter, setter};
      }

      template <typename CoordHelper, typename StoredVector>
      class BoundCoordGetter {

      public:
        using Stored_t = StoredVector;

        using CoordHelper_t = CoordHelper;
        using Vector_t = typename CoordHelper_t::Vector_t;
        using Scalar_t = typename CoordHelper_t::Scalar_t;
        using Getter_t = typename CoordHelper_t::Getter_t;

        BoundCoordGetter(Stored_t& v, CoordHelper_t coordManager) : fCoord(coordManager), fVector(v)
        {}

        BoundCoordGetter(Stored_t& v, Getter_t getter) : fCoord(getter), fVector(v) {}

        Scalar_t get() const { return manager().get(vector()); }

        Scalar_t operator()() const { return get(); }

        operator Scalar_t() const { return manager().get(vector()); }

      protected:
        CoordHelper_t const& manager() const { return fCoord; }
        Stored_t& vector() const { return fVector; }

      private:
        CoordHelper_t fCoord; 
        Stored_t& fVector;    
      };                      // class VectorCoordGetter<>

      template <typename CoordHelper, typename StoredVector>
      class BoundCoordManager : public BoundCoordGetter<CoordHelper, StoredVector> {
        using Base_t = BoundCoordGetter<CoordHelper, StoredVector>;

      public:
        using typename Base_t::Stored_t;

        using CoordHelper_t = CoordHelper;
        using Vector_t = typename CoordHelper_t::Vector_t;
        using Scalar_t = typename CoordHelper_t::Scalar_t;
        using Getter_t = typename CoordHelper_t::Getter_t;
        using Setter_t = typename CoordHelper_t::Setter_t;

        BoundCoordManager(Stored_t& v, CoordHelper_t coordManager) : Base_t(v, coordManager) {}

        BoundCoordManager(Stored_t& v, Getter_t getter, Setter_t setter)
          : Base_t(v, CoordHelper_t(getter, setter))
        {}

        BoundCoordManager& operator=(Scalar_t c)
        {
          Base_t::manager().set(Base_t::vector(), c);
          return *this;
        }

        BoundCoordManager& operator+=(Scalar_t c)
        {
          Base_t::manager().incr(Base_t::vector(), c);
          return *this;
        }

        BoundCoordManager& operator-=(Scalar_t c)
        {
          Base_t::manager().decr(Base_t::vector(), c);
          return *this;
        }

        BoundCoordManager& operator*=(Scalar_t f)
        {
          Base_t::manager().mult(Base_t::vector(), f);
          return *this;
        }

        BoundCoordManager& operator/=(Scalar_t f)
        {
          Base_t::manager().div(Base_t::vector(), f);
          return *this;
        }

      }; // class BoundCoordManager

      //------------------------------------------------------------------------

    } // namespace details

    // BEGIN Geometry group ------------------------------------------------------

    namespace extra {

      template <typename T>
      constexpr T roundValue0(T value, T tol)
      {
        return (std::abs(value) < tol) ? T{} : value;
      } // roundValue0()

      template <typename T>
      constexpr T roundValue01(T value, T tol)
      {
        if (std::abs(value) < tol) return 0.;
        if (std::abs(std::abs(value) - 1.) < tol) return (value > 0.) ? 1. : -1.;
        return value;
      } // roundValue01()

    } // namespace extra

    // --- BEGIN Vector coordinate access abstraction --------------------------

    template <typename Vector>
    constexpr unsigned int dimension()
    {
      return details::dimension<Vector>();
    }
    template <typename Vector>
    constexpr unsigned int dimension(Vector&&)
    {
      return dimension<Vector>();
    }

    template <typename Vector>
    constexpr std::array<std::size_t, vect::dimension<Vector>()> indices();
    template <typename Vector>
    constexpr auto indices(Vector const&) -> decltype(indices<Vector>());

    template <typename Vector>
    using coordinate_t = details::VectorScalar_t<Vector>;

    template <typename Vector, typename Coords>
    constexpr Vector makeFromCoords(Coords&& coords);

    template <typename Vector, typename Coords>
    unsigned int fillCoords(Coords& dest, Vector const& src);

    template <typename Vector>
    using CoordReader_t = decltype(details::makeCoordReader(&Vector::X));

    template <typename Vector>
    using CoordManager_t = decltype(details::makeCoordManager(&Vector::X, &Vector::SetX));

    template <typename Vector>
    static constexpr auto XcoordManager = details::makeCoordManager(&Vector::X, &Vector::SetX);

    template <typename Vector>
    static constexpr auto XcoordManager<Vector const> = details::makeCoordReader(&Vector::X);

    template <typename Vector>
    static constexpr auto const YcoordManager =
      details::makeCoordManager(&Vector::Y, &Vector::SetY);

    template <typename Vector>
    static constexpr auto YcoordManager<Vector const> = details::makeCoordReader(&Vector::Y);

    template <typename Vector>
    static constexpr auto ZcoordManager = details::makeCoordManager(&Vector::Z, &Vector::SetZ);

    template <typename Vector>
    static constexpr auto ZcoordManager<Vector const> = details::makeCoordReader(&Vector::Z);

    template <typename Vector>
    static constexpr auto TcoordManager = details::makeCoordManager(&Vector::T, &Vector::SetT);

    template <typename Vector>
    static constexpr auto TcoordManager<Vector const> = details::makeCoordReader(&Vector::T);

    template <typename Vector>
    constexpr auto coordManager(unsigned int n);

    template <typename Vector>
    constexpr auto coordManager(unsigned int n, Vector& v);

    template <typename Vector>
    constexpr auto coordManagers();
    template <typename Vector>
    constexpr auto coordManagers(Vector&&);

    template <typename Vector>
    constexpr auto coordReaders();
    template <typename Vector>
    constexpr auto coordReaders(Vector&&);

    template <typename Vector>
    constexpr auto bindCoord(Vector const& v, CoordReader_t<Vector> helper)
    {
      return details::BoundCoordGetter<CoordReader_t<Vector>, Vector const>(v, helper);
    }

    template <typename Vector>
    auto bindCoord(Vector& v, CoordManager_t<Vector> helper)
      -> details::BoundCoordManager<CoordManager_t<Vector>, Vector>
    {
      return {v, helper};
    }

    template <typename Vector>
    auto Xcoord(Vector& v)
    {
      return bindCoord(v, XcoordManager<Vector>);
    }

    template <typename Vector>
    auto Ycoord(Vector& v)
    {
      return bindCoord<Vector>(v, YcoordManager<Vector>);
    }

    template <typename Vector>
    auto Zcoord(Vector& v)
    {
      return bindCoord<Vector>(v, ZcoordManager<Vector>);
    }

    template <typename Vector>
    auto Tcoord(Vector& v)
    {
      return bindCoord<Vector>(v, TcoordManager<Vector>);
    }

    template <typename Vector>
    auto coord(Vector& v, unsigned int n) noexcept;

    template <typename Vector>
    auto coord(Vector& v, Coordinate n) noexcept;

    template <typename Vector>
    constexpr auto bindCoordManagers(Vector& v);

    template <typename Vector>
    constexpr auto bindCoordReaders(Vector const& v);

    template <typename Vector>
    class CoordConstIterator {

    public:
      using iterator_t = CoordConstIterator<Vector>; 

    private:
      using Vector_t = Vector; 
      using Coord_t = vect::coordinate_t<Vector_t>;

      Vector* v = nullptr;
      unsigned int index = std::numeric_limits<unsigned int>::max();

      decltype(auto) access(unsigned int c) const { return vect::bindCoord(v, index); }

      iterator_t copy() const { return *this; }

    public:
      using value_type = std::remove_cv_t<Coord_t>;
      using difference_type = std::ptrdiff_t;
      using reference = Coord_t;
      using pointer = Coord_t const*;
      using iterator_category = std::random_access_iterator_tag;

      CoordConstIterator() = default;

      CoordConstIterator(Vector& v, unsigned int index = 0) : v(&v), index(index) {}

      // --- BEGIN Access ----------------------------------------------------------

      reference operator*() const { return vect::coord(*v, index); }

      reference operator[](difference_type n) const { return vect::coord(*v, index + n); }

      // --- END Access ------------------------------------------------------------

      // --- BEGIN Moving ----------------------------------------------------------

      iterator_t& operator++()
      {
        ++index;
        return *this;
      }

      iterator_t operator++(int)
      {
        iterator_t it(*this);
        ++index;
        return it;
      }

      iterator_t& operator--()
      {
        --index;
        return *this;
      }

      iterator_t operator--(int)
      {
        iterator_t it(*this);
        --index;
        return it;
      }

      iterator_t& operator+=(difference_type d)
      {
        index += d;
        return *this;
      }

      iterator_t operator+(difference_type d) const { return (copy() += d); }

      iterator_t& operator-=(difference_type d)
      {
        index -= d;
        return *this;
      }

      iterator_t operator-(difference_type d) const { return (copy() -= d); }

      iterator_t operator-(iterator_t const& other) const { return index - other.index; }

      // --- END Moving ------------------------------------------------------------

      // --- BEGIN Comparison operators --------------------------------------------

      bool operator==(iterator_t const& other) const
      {
        return (v == other.v) && (index == other.index);
      }

      bool operator!=(iterator_t const& other) const
      {
        return (v != other.v) || (index != other.index);
      }

      bool operator<=(iterator_t const& other) const
      {
        return (v == other.v) && (index <= other.index);
      }

      bool operator>=(iterator_t const& other) const
      {
        return (v == other.v) && (index >= other.index);
      }

      bool operator<(iterator_t const& other) const
      {
        return (v == other.v) && (index < other.index);
      }

      bool operator>(iterator_t const& other) const
      {
        return (v == other.v) && (index > other.index);
      }

      // --- END Comparison operators ----------------------------------------------

    }; // class CoordConstIterator

    template <typename Vector>
    CoordConstIterator<Vector> operator+(typename CoordConstIterator<Vector>::difference_type n,
                                         CoordConstIterator<Vector> const& v)
    {
      return v + n;
    }

    template <typename Vector>
    auto vector_cbegin(Vector const& v)
    {
      return CoordConstIterator(v, 0);
    }

    template <typename Vector>
    auto vector_cend(Vector const& v)
    {
      return CoordConstIterator(v, vect::dimension(v));
    }

    template <typename Dest, typename Source>
    Dest convertTo(Source const& v);

    template <typename Dest, typename Source>
    std::vector<Dest> convertCollTo(std::vector<Source> const& coll);

    template <typename Vector, typename Pred>
    Vector transformCoords(Vector const& v, Pred&& pred);

    // --- END Vector coordinate access abstraction ----------------------------

    // --- BEGIN Functions for common vector operations ------------------------

    template <typename Vector, typename Scalar>
    Vector rounded0(Vector const& v, Scalar tol)
    {
      return transformCoords(v, [tol](auto c) { return extra::roundValue0(c, tol); });
    }

    template <typename Vector, typename Scalar>
    void round0(Vector& v, Scalar tol)
    {
      v = rounded0(v, tol);
    }

    template <typename Vector, typename Scalar>
    Vector rounded01(Vector const& v, Scalar tol)
    {
      return transformCoords(v, [tol](auto c) { return extra::roundValue01(c, tol); });
    }

    template <typename Vector, typename Scalar>
    void round01(Vector& v, Scalar tol)
    {
      v = rounded01(v, tol);
    }

    template <typename Vector>
    bool isfinite(Vector const& v);

    template <typename Vector>
    Vector normalize(Vector const& v)
    {
      return v.Unit();
    }

    template <typename Vector>
    Vector cross(Vector const& a, Vector const& b)
    {
      return a.Cross(b);
    }

    template <typename Vector, typename OtherVector>
    constexpr auto dot(Vector const& a, OtherVector const& b)
    {
      return a.Dot(b);
    }

    template <typename Vector>
    auto mag2(Vector const& v)
    {
      return v.Mag2();
    }

    template <typename Vector>
    auto norm(Vector const& v)
    {
      return v.Mag();
    }

    template <typename Vector>
    auto mixedProduct(Vector const& a, Vector const& b, Vector const& c)
    {
      return dot(cross(a, b), c);
    }

    // --- END Functions for common vector operations --------------------------

    template <unsigned int N = 3U>
    class MiddlePointAccumulatorDim {
      static constexpr unsigned int Dim = N; 
      std::array<Length_t, Dim> fSums;       
      double fW = 0.0;                       

    public:
      MiddlePointAccumulatorDim() { fSums.fill(0.); }

      template <typename BeginIter, typename EndIter>
      MiddlePointAccumulatorDim(BeginIter begin, EndIter end) : MiddlePointAccumulatorDim()
      {
        add(begin, end);
      }

      // --- BEGIN Result query ------------------------------------------------

      bool empty() const { return fW == 0.0; }

      double weight() const { return fW; }

      template <typename Point>
      Point middlePointAs() const
      //  { return { fSums[0] / fW, fSums[1] / fW, fSums[2] / fW }; }
      {
        return makeWeightedPoint<Point>(1.0 / fW);
      }

      Point_t middlePoint() const { return middlePointAs<Point_t>(); }

      // --- END Result query --------------------------------------------------

      // --- BEGIN Addition of points ------------------------------------------

      template <typename Point>
      void add(Point const& p)
      {
        std::size_t ic = 0U;
        for (auto c : vect::bindCoordManagers(p))
          fSums[ic++] += c();
        fW += 1.0;
      }

      template <typename Point>
      void add(Point const& p, double weight)
      {
        std::size_t ic = 0U;
        for (auto c : vect::bindCoordManagers(p))
          fSums[ic++] += weight * c();
        fW += weight;
      }

      template <typename BeginIter, typename EndIter>
      void add(BeginIter begin, EndIter end)
      {
        std::for_each(begin, end, [this](auto const& p) { this->add(p); });
      }

      void clear()
      {
        fSums.fill(0.);
        fW = 0.0;
      }

      // --- END Addition of points --------------------------------------------

    private:
      using IndexSequence_t = std::make_index_sequence<Dim>;

      template <typename Point, std::size_t... I>
      Point makePointImpl(std::index_sequence<I...>) const
      {
        return {fSums.operator[](I)...};
      }

      template <typename Point, std::size_t... I>
      Point makeWeightedPointImpl(double w, std::index_sequence<I...>) const
      {
        return {(fSums.operator[](I) * w)...};
      }

      template <typename Point>
      Point makePoint() const
      {
        return vect::makeFromCoords<Point>(fSums);
      }

      template <typename Point>
      Point makeWeightedPoint(double w) const
      {
        return makeWeightedPointImpl<Point>(w, IndexSequence_t{});
      }

    }; // MiddlePointAccumulatorDim()

    using MiddlePointAccumulator = MiddlePointAccumulatorDim<3U>;

    // --- BEGIN Middle point functions ----------------------------------------

    template <typename Point, typename BeginIter, typename EndIter>
    Point middlePointAs(BeginIter begin, EndIter end)
    {
      constexpr auto Dim = vect::dimension<Point>();
      return MiddlePointAccumulatorDim<Dim>(begin, end).template middlePointAs<Point>();
    }

    template <typename BeginIter, typename EndIter>
    Point_t middlePoint(BeginIter begin, EndIter end)
    {
      return middlePointAs<Point_t>(begin, end);
    }

    template <typename Point>
    Point middlePoint(std::initializer_list<Point> points)
    {
      constexpr auto Dim = vect::dimension<Point>();
      return MiddlePointAccumulatorDim<Dim>(points.begin(), points.end())
        .template middlePointAs<Point>();
    }

    // --- END Middle point functions ------------------------------------------

    // --- BEGIN Support for LArSoft geometry vectors --------------------------

    template <typename Point>
    Point_t toPoint(Point const& p)
    {
      return vect::convertTo<Point_t>(p);
    }

    template <typename Vector>
    Vector_t toVector(Vector const& v)
    {
      return vect::convertTo<Vector_t>(v);
    }

    template <typename Point>
    std::vector<Point_t> convertCollToPoint(std::vector<Point> const& coll)
    {
      return convertCollTo<Point_t>(coll);
    }

    template <typename Vector>
    std::vector<Vector_t> convertCollToVector(std::vector<Vector> const& coll)
    {
      return convertCollTo<Vector_t>(coll);
    }

    template <typename Coords>
    GENVECTOR_CONSTEXPR Point_t makePointFromCoords(Coords&& coords)
    {
      return makeFromCoords<Point_t>(std::forward<Coords>(coords));
    }

    template <typename Coords>
    GENVECTOR_CONSTEXPR Vector_t makeVectorFromCoords(Coords&& coords)
    {
      return makeFromCoords<Vector_t>(std::forward<Coords>(coords));
    }

    // --- END Support for LArSoft geometry vectors ----------------------------

    // END Geometry group ------------------------------------------------------

  } // namespace vect

} // namespace geo

//------------------------------------------------------------------------------
//---  STL specialization for ROOT GenVector vectors
//------------------------------------------------------------------------------
namespace ROOT::Math {


  // --- BEGIN 2D vectors ------------------------------------------------------
  template <class CoordSystem, class Tag>
  decltype(auto) begin(ROOT::Math::PositionVector2D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) cbegin(ROOT::Math::PositionVector2D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) end(ROOT::Math::PositionVector2D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) cend(ROOT::Math::PositionVector2D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) begin(ROOT::Math::DisplacementVector2D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) cbegin(ROOT::Math::DisplacementVector2D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) end(ROOT::Math::DisplacementVector2D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) cend(ROOT::Math::DisplacementVector2D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  // --- END 2D vectors --------------------------------------------------------

  // --- BEGIN 3D vectors ------------------------------------------------------
  template <class CoordSystem, class Tag>
  decltype(auto) begin(ROOT::Math::PositionVector3D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) cbegin(ROOT::Math::PositionVector3D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) end(ROOT::Math::PositionVector3D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) cend(ROOT::Math::PositionVector3D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) begin(ROOT::Math::DisplacementVector3D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) cbegin(ROOT::Math::DisplacementVector3D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) end(ROOT::Math::DisplacementVector3D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  template <class CoordSystem, class Tag>
  decltype(auto) cend(ROOT::Math::DisplacementVector3D<CoordSystem, Tag> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  // --- END 3D vectors --------------------------------------------------------

  // --- BEGIN 4D vectors ------------------------------------------------------
  template <class CoordSystem>
  decltype(auto) begin(ROOT::Math::LorentzVector<CoordSystem> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem>
  decltype(auto) cbegin(ROOT::Math::LorentzVector<CoordSystem> const& v)
  {
    return geo::vect::vector_cbegin(v);
  }

  template <class CoordSystem>
  decltype(auto) end(ROOT::Math::LorentzVector<CoordSystem> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  template <class CoordSystem>
  decltype(auto) cend(ROOT::Math::LorentzVector<CoordSystem> const& v)
  {
    return geo::vect::vector_cend(v);
  }

  // --- END 4D vectors --------------------------------------------------------


} // namespace ROOT::Math

//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
//---  template specializations for standard geometry vectors
//---
namespace geo::vect {

  //----------------------------------------------------------------------------
  template <>
  inline auto norm(geo::Vector_t const& v)
  {
    return v.R();
  }

  //----------------------------------------------------------------------------

} // namespace geo::vect

//------------------------------------------------------------------------------
//---  template implementation

namespace geo::vect::details {

  //----------------------------------------------------------------------------
  template <typename>
  struct AlwaysFalse : std::false_type {};

  // constexpr variant of forward (from StackExchange)
  template <typename T>
  constexpr T&& constexpr_forward(std::remove_reference_t<T>& t)
  {
    return static_cast<T&&>(t);
  }

  template <typename T>
  constexpr T&& constexpr_forward(std::remove_reference_t<T>&& t)
  {
    static_assert(!std::is_lvalue_reference<T>(),
                  "template argument substituting T is an lvalue reference type");
    return static_cast<T&&>(t);
  }

  //----------------------------------------------------------------------------
  template <typename Vector, typename /* = void */>
  struct DimensionImpl : public std::integral_constant<unsigned int,
                                                       HasT<Vector>() ? 4U :
                                                       HasZ<Vector>() ? 3U :
                                                       HasY<Vector>() ? 2U :
                                                       HasX<Vector>() ? 1U :
                                                                        0U> {};

  template <typename Array>
  struct DimensionImpl<Array, std::enable_if_t<(std::extent_v < Array >> 0)>>
    : public std::integral_constant<unsigned int, std::extent_v<Array>> {};

  template <typename T, std::size_t Dim>
  struct DimensionImpl<std::array<T, Dim>, void>
    : public std::integral_constant<unsigned int, Dim> {};

  template <typename T>
  struct DimensionImpl<lar::util::simple_geo::Point2D<T>, void>
    : public std::integral_constant<unsigned int, 2> {};

  template <typename T>
  struct DimensionImpl<lar::util::simple_geo::Point3D<T>, void>
    : public std::integral_constant<unsigned int, 3> {};

  //----------------------------------------------------------------------------
  template <typename Vector>
  using VectorIndices_t = std::make_index_sequence<dimension<Vector>()>;

  template <typename Vector>
  constexpr auto makeVectorIndices()
  {
    return VectorIndices_t<Vector>{};
  }

  template <typename Vector>
  constexpr auto makeVectorIndices(Vector&&)
  {
    return makeVectorIndices<Vector>();
  }

  template <typename T, T... Indices>
  constexpr auto makeIndexSeqImpl(std::integer_sequence<T, Indices...>)
  // BUG the double brace syntax is required to work around clang bug 21629
  // (https://bugs.llvm.org/show_bug.cgi?id=21629)
  {
    return std::array<T, sizeof...(Indices)>{{Indices...}};
  }

  // fill a sequence object with the first `N` values of type `T`.
  template <typename T, T N>
  constexpr auto makeIndexSeq()
  {
    return makeIndexSeqImpl<T>(std::make_integer_sequence<T, N>{});
  }

  //----------------------------------------------------------------------------
  template <std::size_t I, typename Data>
  constexpr auto accessElement(Data&& data)
  {
    return data[I];
  }

  template <typename Vector, typename Coords, std::size_t... Indices>
  constexpr Vector makeFromCoordsImpl(Coords&& coords, std::index_sequence<Indices...>)
  {
    return {accessElement<Indices>(constexpr_forward<Coords>(coords))...};
  }

  //----------------------------------------------------------------------------
  template <typename Vector>
  constexpr CoordManager_t<Vector> NoCoordManager{nullptr, nullptr};

  template <typename Vector, unsigned int Dim = dimension<Vector>()>
  struct CoordManagerImpl;

  template <typename Vector>
  struct CoordManagerImpl<Vector, 1U> {
    static auto get(unsigned int n) noexcept
    {
      return (n == 0) ? XcoordManager<Vector> : NoCoordManager<Vector>;
    }
  };

  template <typename Vector>
  struct CoordManagerImpl<Vector, 2U> {
    static auto get(unsigned int n) noexcept
    {
      return (n == 1) ? YcoordManager<Vector> : CoordManagerImpl<Vector, 1U>::get(n);
    }
  };

  template <typename Vector>
  struct CoordManagerImpl<Vector, 3U> {
    static auto get(unsigned int n) noexcept
    {
      return (n == 2) ? ZcoordManager<Vector> : CoordManagerImpl<Vector, 2U>::get(n);
    }
  };

  template <typename Vector>
  struct CoordManagerImpl<Vector, 4U> {
    static auto get(unsigned int n) noexcept
    {
      return (n == 3) ? TcoordManager<Vector> : CoordManagerImpl<Vector, 3U>::get(n);
    }
  };

  //----------------------------------------------------------------------------
  template <typename Vector, unsigned int N>
  struct CoordManagersImplBase {
    static constexpr unsigned int Dim = N;

    using Manager_t = decltype(XcoordManager<Vector>);
    using Return_t = std::array<Manager_t, Dim>;

    static_assert(dimension<Vector>() == Dim, "Inconsistent vector size.");
  };

  template <typename Vector, unsigned int N>
  struct CoordManagersImpl;

  template <typename Vector>
  struct CoordManagersImpl<Vector, 2U> : private CoordManagersImplBase<Vector, 2U> {
    using Base_t = CoordManagersImplBase<Vector, 2U>;
    using typename Base_t::Return_t;
    static constexpr Return_t get()
    {
      // BUG the double brace syntax is required to work around clang bug 21629
      // (https://bugs.llvm.org/show_bug.cgi?id=21629)
      return {{XcoordManager<Vector>, YcoordManager<Vector>}};
    }
  };

  template <typename Vector>
  struct CoordManagersImpl<Vector, 3U> : private CoordManagersImplBase<Vector, 3U> {
    using Base_t = CoordManagersImplBase<Vector, 3U>;
    using typename Base_t::Return_t;
    static constexpr Return_t get()
    {
      // BUG the double brace syntax is required to work around clang bug 21629
      // (https://bugs.llvm.org/show_bug.cgi?id=21629)
      return {{XcoordManager<Vector>, YcoordManager<Vector>, ZcoordManager<Vector>}};
    }
  };

  template <typename Vector>
  struct CoordManagersImpl<Vector, 4U> : private CoordManagersImplBase<Vector, 4U> {
    using Base_t = CoordManagersImplBase<Vector, 4U>;
    using typename Base_t::Return_t;
    static constexpr Return_t get()
    {
      // BUG the double brace syntax is required to work around clang bug 21629
      // (https://bugs.llvm.org/show_bug.cgi?id=21629)
      return {{XcoordManager<Vector>,
               YcoordManager<Vector>,
               ZcoordManager<Vector>,
               TcoordManager<Vector>}};
    }
  };

  //----------------------------------------------------------------------------
  template <typename Vector, unsigned int N>
  struct BindCoordManagersImplBase {
    static constexpr unsigned int Dim = N;

    using Manager_t = decltype(Xcoord(std::declval<Vector>()));
    using Return_t = std::array<Manager_t, Dim>;

    static_assert(dimension<Vector>() == Dim, "Inconsistent vector size.");
  };

  template <typename Vector, unsigned int N>
  struct BindCoordManagersImpl;

  template <typename Vector>
  struct BindCoordManagersImpl<Vector, 2U> : private BindCoordManagersImplBase<Vector, 2U> {
    using Base_t = CoordManagersImplBase<Vector, 2U>;
    using typename Base_t::Return_t;
    static Return_t bind(Vector& v)
    {
      // BUG the double brace syntax is required to work around clang bug 21629
      // (https://bugs.llvm.org/show_bug.cgi?id=21629)
      return {{Xcoord(v), Ycoord(v)}};
    }
  };

  template <typename Vector>
  struct BindCoordManagersImpl<Vector, 3U> : private BindCoordManagersImplBase<Vector, 3U> {
    using Base_t = BindCoordManagersImplBase<Vector, 3U>;
    using typename Base_t::Return_t;
    static Return_t bind(Vector& v)
    {
      // BUG the double brace syntax is required to work around clang bug 21629
      // (https://bugs.llvm.org/show_bug.cgi?id=21629)
      return {{Xcoord(v), Ycoord(v), Zcoord(v)}};
    }
  };

  template <typename Vector>
  struct BindCoordManagersImpl<Vector, 4U> : private BindCoordManagersImplBase<Vector, 4U> {
    using Base_t = BindCoordManagersImplBase<Vector, 4U>;
    using typename Base_t::Return_t;
    static Return_t bind(Vector& v)
    {
      // BUG the double brace syntax is required to work around clang bug 21629
      // (https://bugs.llvm.org/show_bug.cgi?id=21629)
      return {{Xcoord(v), Ycoord(v), Zcoord(v), Tcoord(v)}};
    }
  };

  //----------------------------------------------------------------------------
  template <typename Dest, typename Source>
  struct ConvertToImplBase {
    static_assert(dimension<Source>() == dimension<Dest>(),
                  "Source and destination vectors must have the same dimension.");
  };

  template <typename Dest, typename Source, unsigned int Dim>
  struct ConvertArrayTo : private ConvertToImplBase<Dest, Source> {
    static_assert(std::is_arithmetic_v<std::decay_t<decltype(std::declval<Source>()[0])>>);
    static Dest convert(Source const& v) { return geo::vect::makeFromCoords<Dest>(v); }
  };

  // will handle cases with specific dimensionality
  template <typename Dest, typename Source, unsigned int Dim>
  struct ConvertToImplDim;

  template <typename Dest, typename Source>
  struct ConvertToImpl : public ConvertToImplDim<Dest, Source, dimension<Source>()> {};

  template <typename Dest, typename T, std::size_t Dim>
  struct ConvertToImpl<Dest, std::array<T, Dim>>
    : public ConvertArrayTo<Dest, std::array<T, Dim>, Dim> {};

  template <typename Dest, typename T, std::size_t Dim>
  struct ConvertToImpl<Dest, T[Dim]> : public ConvertArrayTo<Dest, T[Dim], Dim> {};

  // special special implementation for T*, since we can't guess its dimension
  template <typename Dest, typename T>
  struct ConvertToImpl<Dest, T*> {
    static_assert(std::is_arithmetic_v<T>);
    static Dest convert(T* v) { return geo::vect::makeFromCoords<Dest>(v); }
  }; // struct ConvertToImpl<T*>

  // handled cases with specific dimensionality
  template <typename Dest, typename Source, unsigned int Dim>
  struct ConvertToImplDim {
    // trivial to do: open a feature request!
    static_assert(AlwaysFalse<Dest>(), "This vector dimensionality is not implemented yet.");
  };

  template <typename Dest, typename Source>
  struct ConvertToImplDim<Dest, Source, 2U> : private ConvertToImplBase<Dest, Source> {
    static Dest convert(Source const& v) { return {Xcoord(v)(), Ycoord(v)()}; }
  };

  template <typename Dest, typename Source>
  struct ConvertToImplDim<Dest, Source, 3U> : private ConvertToImplBase<Dest, Source> {
    static Dest convert(Source const& v) { return {Xcoord(v)(), Ycoord(v)(), Zcoord(v)()}; }
  };

  template <typename Dest, typename Source>
  struct ConvertToImplDim<Dest, Source, 4U> : private ConvertToImplBase<Dest, Source> {
    static Dest convert(Source const& v)
    {
      return {Xcoord(v)(), Ycoord(v)(), Zcoord(v)(), Tcoord(v)()};
    }
  };

  template <typename Dest, typename Source>
  struct ConvertToDispatcher : public ConvertToImpl<Dest, Source> {};

  // special pass-through case
  template <typename Vector>
  struct ConvertToDispatcher<Vector, Vector> {
    static constexpr decltype(auto) convert(Vector const& v) { return v; }
  };

  //----------------------------------------------------------------------------
  template <typename Point, std::size_t... I>
  bool isfiniteImpl(Point const& point, std::index_sequence<I...>)
  {
    return extended_and(std::isfinite(coord(point, I).get())...);
  }

  //----------------------------------------------------------------------------

} // namespace geo::vect::details

//------------------------------------------------------------------------------
template <typename Vector>
constexpr std::array<std::size_t, geo::vect::dimension<Vector>()> geo::vect::indices()
{
  return details::makeIndexSeq<std::size_t, dimension<Vector>()>();
}

template <typename Vector>
constexpr auto geo::vect::indices(Vector const&) -> decltype(indices<Vector>())
{
  return indices<Vector>();
}

//------------------------------------------------------------------------
template <typename Vector, typename Coords>
constexpr Vector geo::vect::makeFromCoords(Coords&& coords)
{
  using namespace geo::vect::details;
  return makeFromCoordsImpl<Vector>(constexpr_forward<Coords>(coords), makeVectorIndices<Vector>());
}

//------------------------------------------------------------------------------
template <typename Vector>
constexpr auto geo::vect::coordManager(unsigned int n)
{
  return details::CoordManagerImpl<Vector>::get(n);
}

//------------------------------------------------------------------------------
template <typename Vector>
auto geo::vect::coord(Vector& v, unsigned int n) noexcept
{
  return vect::bindCoord<Vector>(v, coordManager<Vector>(n));
}

//------------------------------------------------------------------------------
template <typename Vector>
auto geo::vect::coord(Vector& v, Coordinate coord) noexcept
{
  return coord(v, to_int(coord));
}

//------------------------------------------------------------------------
template <typename Vector, typename Coords>
unsigned int geo::vect::fillCoords(Coords& dest, Vector const& src)
{
  // this is simpler than makeFromCoords() because doesn't attempt constexpr
  for (std::size_t i = 0; i < vect::dimension(src); ++i)
    dest[i] = vect::coord(src, i);
  return vect::dimension(src);
}

//------------------------------------------------------------------------------
template <typename Vector>
constexpr auto geo::vect::coordManagers()
{
  using PlainVector = std::remove_reference_t<Vector>;
  return details::CoordManagersImpl<PlainVector, dimension<PlainVector>()>::get();
}

template <typename Vector>
constexpr auto geo::vect::coordManagers(Vector&&)
{
  return coordManagers<Vector>();
}

template <typename Vector>
constexpr auto geo::vect::coordReaders()
{
  using ConstVector = std::add_const_t<std::remove_reference_t<Vector>>;
  return details::CoordManagersImpl<ConstVector, dimension<ConstVector>()>::get();
}

template <typename Vector>
constexpr auto geo::vect::coordReaders(Vector&&)
{
  return coordReaders<Vector>();
}

//------------------------------------------------------------------------------
template <typename Vector>
constexpr auto geo::vect::bindCoordManagers(Vector& v)
{
  return details::BindCoordManagersImpl<Vector, dimension<Vector>()>::bind(v);
}

template <typename Vector>
constexpr auto geo::vect::bindCoordReaders(Vector const& v)
{
  using ConstVector = std::add_const_t<std::remove_reference_t<Vector>>;
  return details::BindCoordManagersImpl<ConstVector, dimension<ConstVector>()>::bind(v);
}

//------------------------------------------------------------------------------
template <typename Dest, typename Source>
Dest geo::vect::convertTo(Source const& v)
{
  return details::ConvertToDispatcher<Dest, Source>::convert(v);
}

//----------------------------------------------------------------------------
template <typename Dest, typename Source>
std::vector<Dest> geo::vect::convertCollTo(std::vector<Source> const& coll)
{
  std::vector<Dest> dest;
  dest.reserve(coll.size());
  std::transform(
    coll.begin(), coll.end(), std::back_inserter(dest), geo::vect::convertTo<Dest, Source>);
  return dest;
}

//------------------------------------------------------------------------
template <typename Vector, typename Pred>
Vector geo::vect::transformCoords(Vector const& v, Pred&& pred)
{
  details::CoordinateArray_t<Vector> values;
  unsigned int i = 0;
  for (auto c : bindCoordReaders(v))
    values[i++] = pred(c());
  return makeFromCoords<Vector>(values);
}

//------------------------------------------------------------------------------
template <typename Vector>
bool geo::vect::isfinite(Vector const& v)
{
  return details::isfiniteImpl(v, details::makeVectorIndices<Vector>());
}

//------------------------------------------------------------------------------

#endif // LARCOREALG_GEOMETRY_GEO_VECTORS_UTILS_H