diff --git a/Core/include/Acts/Geometry/CompositePortalLink.hpp b/Core/include/Acts/Geometry/CompositePortalLink.hpp new file mode 100644 index 00000000000..ccf7e8bcf17 --- /dev/null +++ b/Core/include/Acts/Geometry/CompositePortalLink.hpp @@ -0,0 +1,99 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/Tolerance.hpp" +#include "Acts/Geometry/PortalLinkBase.hpp" + +#include + +#include + +namespace Acts { + +/// Composite portal links can graft together other portal link instances, for +/// example grids that could not be merged due to invalid binnings. +/// +/// +-------+ +-------+ +/// | | | | +/// | | | | +/// | | | | +/// +-------+ | | +/// | | | | +/// | | + +-------+ +/// | | | | +/// +-------+ | | +/// | | | | +/// | | +-------+ +/// | | | | +/// +-------+ +-------+ +/// +/// During resolution, it will consult each of it's children and return +/// the result on the first surface where the lookup position is within +/// bounds. +class CompositePortalLink final : public PortalLinkBase { + public: + /// Construct a composite portal from two arbitrary other portal links. The + /// only requirement is that the portal link surfaces are mergeable. + /// @param a The first portal link + /// @param b The second portal link + /// @param direction The binning direction + /// @param flatten If true, the composite will flatten any nested composite + CompositePortalLink(std::unique_ptr a, + std::unique_ptr b, BinningValue direction, + bool flatten = true); + + /// Print the composite portal link + /// @param os The output stream + void toStream(std::ostream& os) const override; + + /// Resolve the volume for a 2D position + /// @note This will transform the position to global coordinates before + /// consulting its children. + /// @note @p position is assumed to be on surface + /// @param gctx The geometry context + /// @param position The 2D position + /// @param tolerance The on-surface tolerance + Result resolveVolume( + const GeometryContext& gctx, const Vector2& position, + double tolerance = s_onSurfaceTolerance) const override; + + /// Resolve the volume for a 3D position + /// @note @p position is assumed to be on surface + /// @param gctx The geometry context + /// @param position The 3D position + /// @param tolerance The tolerance + Result resolveVolume( + const GeometryContext& gctx, const Vector3& position, + double tolerance = s_onSurfaceTolerance) const override; + + /// Get the depth of the composite tree + /// @return The depth + std::size_t depth() const; + + /// Get the number of children + /// @return The number of children + std::size_t size() const; + + private: + /// Helper function to construct a merged surface from two portal links along + /// a given direction + /// @param a The first portal link + /// @param b The second portal link + /// @param direction The merging direction + /// @return The merged surface + static std::shared_ptr mergedSurface(const PortalLinkBase* a, + const PortalLinkBase* b, + BinningValue direction); + + boost::container::small_vector, 4> + m_children{}; +}; + +} // namespace Acts diff --git a/Core/include/Acts/Geometry/GridPortalLink.hpp b/Core/include/Acts/Geometry/GridPortalLink.hpp new file mode 100644 index 00000000000..400b86e6b82 --- /dev/null +++ b/Core/include/Acts/Geometry/GridPortalLink.hpp @@ -0,0 +1,615 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/Tolerance.hpp" +#include "Acts/Geometry/PortalLinkBase.hpp" +#include "Acts/Geometry/TrackingVolume.hpp" +#include "Acts/Surfaces/BoundaryTolerance.hpp" +#include "Acts/Surfaces/CylinderSurface.hpp" +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Surfaces/PlaneSurface.hpp" +#include "Acts/Utilities/Grid.hpp" +#include "Acts/Utilities/Logger.hpp" + +#include + +namespace Acts { + +class IGrid; + +template + requires(sizeof...(Axes) <= 2) +class GridPortalLinkT; + +/// GridPortalLink implements a subdivided surface where the target volume +/// depends on the position on the surface. The link can be in two states: +/// 1. One-dimensional binning with an associated *direction* to determine which +/// local coordinate to use for the lookup +/// 2. Two-dimensional binning +/// @note The grid dimensions and boundaries are **required** (and checked) to be +/// consistent with the surface bounds. +class GridPortalLink : public PortalLinkBase { + protected: + /// Constructor from a surface and a direction, for initialization by derived + /// class + /// @param surface The surface + /// @param direction The binning direction + GridPortalLink(std::shared_ptr surface, + BinningValue direction) + : PortalLinkBase(std::move(surface)), m_direction(direction) {} + + public: + /// Factory function for a one-dimensional grid portal link, which allows + /// using template deduction to figure out the right type + /// @tparam axis_t The axis type + /// @param surface The surface + /// @param direction The binning direction + /// @param axis The axis to use for the binning + /// @note The axis boundaries are checked against the bounds of @p surface. + /// @return A unique pointer to the grid portal link + template + static std::unique_ptr> make( + std::shared_ptr surface, BinningValue direction, + axis_t&& axis) { + using enum BinningValue; + if (dynamic_cast(surface.get()) != nullptr) { + if (direction != binZ && direction != binRPhi) { + throw std::invalid_argument{"Invalid binning direction"}; + } + } else if (dynamic_cast(surface.get()) != nullptr && + direction != binR && direction != binPhi) { + throw std::invalid_argument{"Invalid binning direction"}; + } + + return std::make_unique>( + surface, direction, std::forward(axis)); + } + + /// Factory function for a two-dimensional grid portal link, which allows + /// using template deduction to figure out the right type. + /// @tparam axis_1_t The first axis type + /// @tparam axis_2_t The second axis type + /// @param surface The surface + /// @param axis1 The first axis to use for the binning + /// @param axis2 The second axis to use for the binning + /// @note The axis boundaries are checked against the bounds of @p surface. + /// @return A unique pointer to the grid portal link + template + static std::unique_ptr> make( + std::shared_ptr surface, axis_1_t axis1, axis_2_t axis2) { + std::optional direction; + if (dynamic_cast(surface.get()) != nullptr) { + direction = BinningValue::binRPhi; + } else if (dynamic_cast(surface.get()) != nullptr) { + direction = BinningValue::binR; + } + + return std::make_unique>( + surface, direction.value(), std::move(axis1), std::move(axis2)); + } + + /// Factory function for an automatically sized one-dimensional grid. This produces a single-bin grid with boundaries taken from the bounds of @p surface along @p direction. + /// @param surface The surface + /// @param volume The tracking volume + /// @param direction The binning direction + /// @return A unique pointer to the grid portal link + static std::unique_ptr make( + const std::shared_ptr& surface, TrackingVolume& volume, + BinningValue direction); + + /// Merge two grid portal links into a single one. The routine can merge + /// one-dimenaional, tow-dimensional and mixed links. The merge will try to + /// preserve equidistant binning in case bin widths match. Otherwise, the + /// merge falls back to variable binning. + /// + /// 1D merge scenarios: + /// + /// +----------------------------------------------------------+ + /// |Colinear | + /// | | + /// | +-------+-------+-------+ +-------+-------+-------+ | + /// | | | | | | | | | | + /// | | | | | + | | | | | + /// | | | | | | | | | | + /// | +-------+-------+-------+ +-------+-------+-------+ | + /// | +-------+-------+-------+-------+-------+-------+ | + /// | | | | | | | | | + /// | = | | | | | | | | + /// | | | | | | | | | + /// | +-------+-------+-------+-------+-------+-------+ | + /// | | + /// +----------------------------------------------------------+ + /// + /// Two grid along a shared direction are merged along their shared direction + /// + /// +-------------------------------------------------+ + /// |Parallel | + /// | | + /// | +-------+ +-------+ +-------+-------+ | + /// | | | | | | | | | + /// | | | | | | | | | + /// | | | | | | | | | + /// | +-------+ +-------+ +-------+-------+ | + /// | | | | | | | | | + /// | | | + | | = | | | | + /// | | | | | | | | | + /// | +-------+ +-------+ +-------+-------+ | + /// | | | | | | | | | + /// | | | | | | | | | + /// | | | | | | | | | + /// | +-------+ +-------+ +-------+-------+ | + /// | | + /// +-------------------------------------------------+ + /// + /// Two grids along a shared direction a merged in the direction that is + /// orthogonal to their shared direction. + /// + /// +-------------------------------------------+ + /// |Perpendicular | + /// | | + /// | +-------+ | + /// | | | | + /// | | | | + /// | | | | + /// | +-------+ +-------+-------+-------+ | + /// | | | | | | | | + /// | | | + | | | | | + /// | | | | | | | | + /// | +-------+ +-------+-------+-------+ | + /// | | | | + /// | | | | + /// | | | +-------+-------+-------+ | + /// | +-------+ | | | | | + /// | | | | | | + /// | | | | | | + /// | +-------+-------+-------+ | + /// | | | | | | + /// | = | | | | | + /// | | | | | | + /// | +-------+-------+-------+ | + /// | | | | | | + /// | | | | | | + /// | | | | | | + /// | +-------+-------+-------+ | + /// | | + /// +-------------------------------------------+ + /// + /// Two grids whose directions are not shared are merged (ordering does not + /// matter here). The routine will expand one of the grids to match the + /// other's binning, by subdividing the grid in the as-of-yet unbinned + /// direction, while filling all bins with the original bin contents. + /// Afterwards, a conventional mixed-dimension merge is performed. + /// + /// Mixed merge scenarios + /// The order is normalized by always taking the 2D grid as the left hand + /// side. The 1D grid is expanded to match the binning in the as-of-yet + /// unbinned direction with the binning taken from the 2D grid. + /// + /// +-----------------------------------------+ + /// |2D + 1D | + /// | | + /// | +-------+-------+ +-------+-------+ | + /// | | | | | | | | + /// | | | | | | | | + /// | | | | | | | | + /// | +-------+-------+ +-------+-------+ | + /// | | | | | | | | + /// | | | | | | | | + /// | | | | | | | | + /// | +-------+-------+ = +-------+-------+ | + /// | | | | | | | | + /// | | | | | | | | + /// | | | | | | | | + /// | +-------+-------+ +-------+-------+ | + /// | + | | | | + /// | +-------+-------+ | | | | + /// | | | | | | | | + /// | | | | +-------+-------+ | + /// | | | | | + /// | +-------+-------+ | + /// +-----------------------------------------+ + /// + /// +--------------------------------------------------------------+ + /// |2D + 1D | + /// | | + /// | +-------+-------+ +-------+ +-------+-------+-------+ | + /// | | | | | | | | | | | + /// | | | | | | | | | | | + /// | | | | | | | | | | | + /// | +-------+-------+ +-------+ +-------+-------+-------+ | + /// | | | | | | | | | | | + /// | | | | + | | = | | | | | + /// | | | | | | | | | | | + /// | +-------+-------+ +-------+ +-------+-------+-------+ | + /// | | | | | | | | | | | + /// | | | | | | | | | | | + /// | | | | | | | | | | | + /// | +-------+-------+ +-------+ +-------+-------+-------+ | + /// | | + /// +--------------------------------------------------------------+ + /// + /// 2D merges + /// The grids need to already share a common axis. If that is not the case, + /// the merge routine returns a `nullptr`. This can be handled by composite + /// merging as a fallback if needed. + /// Ordering and direction does not matter here. + /// + /// +-----------------------------------------+ + /// |2D + 2D | + /// | | + /// | +-------+-------+ +-------+-------+ | + /// | | | | | | | | + /// | | | | | | | | + /// | | | | | | | | + /// | +-------+-------+ +-------+-------+ | + /// | | | | | | | | + /// | | | | + | | | | + /// | | | | | | | | + /// | +-------+-------+ +-------+-------+ | + /// | | | | | | | | + /// | | | | | | | | + /// | | | | | | | | + /// | +-------+-------+ +-------+-------+ | + /// | | + /// | +-------+-------+-------+-------+ | + /// | | | | | | | + /// | | | | | | | + /// | | | | | | | + /// | +-------+-------+-------+-------+ | + /// | | | | | | | + /// | = | | | | | | + /// | | | | | | | + /// | +-------+-------+-------+-------+ | + /// | | | | | | | + /// | | | | | | | + /// | | | | | | | + /// | +-------+-------+-------+-------+ | + /// | | + /// +-----------------------------------------+ + /// + /// @param a The first grid portal link + /// @param b The second grid portal link + /// @param direction The merging direction + /// @param logger The logger to use for messages + /// @return The merged grid portal link or nullptr if grid merging did not succeed + /// @note The returned pointer can be nullptr, if the grids + /// given are not mergeable, because their binnings are + /// not compatible. This is not an error case, and needs + /// to be handled by th caller! Invalid input is handled + /// via exceptions. + static std::unique_ptr merge( + const GridPortalLink& a, const GridPortalLink& b, BinningValue direction, + const Logger& logger = getDummyLogger()); + + /// Return the associated grid in a type-erased form + /// @return The grid + virtual const IGrid& grid() const = 0; + + /// Set the volume on all grid bins + /// @param volume The volume to set + virtual void setVolume(TrackingVolume* volume) = 0; + + /// Get the number of dimensions of the grid + /// @return The number of dimensions + virtual unsigned int dim() const = 0; + + /// Expand a 1D grid to a 2D one, by using the provided axis along the + /// *missing* direction. + /// @param other The axis to use for the missing direction, + /// can be null for auto determination + /// @return A unique pointer to the 2D grid portal link + virtual std::unique_ptr extendTo2d( + const IAxis* other) const = 0; + + /// The binning direction of the grid + /// @note For 2D grids, this will always be the loc0 + /// direction, depending on the surface type. + /// @return The binning direction + BinningValue direction() const { return m_direction; } + + /// Helper function to fill the bin contents after merging. + /// This called by the merging routine, and requires access to the internal + /// grid state. + /// @param a The first grid portal link + /// @param b The second grid portal link + /// @param merged The merged grid portal link + /// @param direction The merging direction + /// @param logger The logger to use for messages + static void fillMergedGrid(const GridPortalLink& a, const GridPortalLink& b, + GridPortalLink& merged, BinningValue direction, + const Logger& logger); + + /// Helper function that prints a textual representation of the grid with the + /// volume names. + /// @param os The output stream + void printContents(std::ostream& os) const; + + protected: + /// Helper function to check consistency for grid on a cylinder surface + /// @param cyl The cylinder surface + void checkConsistency(const CylinderSurface& cyl) const; + + /// Helper function to check consistency for grid on a disc surface + /// @param disc The disc surface + void checkConsistency(const DiscSurface& disc) const; + + /// Expand a 1D grid to a 2D one for a cylinder surface + /// @param surface The cylinder surface + /// @param other The axis to use for the missing direction, + /// can be null for auto determination + /// @return A unique pointer to the 2D grid portal link + std::unique_ptr extendTo2dImpl( + const std::shared_ptr& surface, + const IAxis* other) const; + + /// Expand a 1D grid to a 2D one for a disc surface + /// @param surface The disc surface + /// @param other The axis to use for the missing direction, + /// can be null for auto determination + /// @return A unique pointer to the 2D grid portal link + std::unique_ptr extendTo2dImpl( + const std::shared_ptr& surface, const IAxis* other) const; + + /// Helper enum to declare which local direction to fill + enum class FillDirection { + loc0, + loc1, + }; + + /// Helper function to fill a 2D grid from a 1D grid, by extending all + /// bins along the different direction. + /// @param dir The direction to fill + /// @param grid1d The 1D grid + /// @param grid2d The 2D grid + static void fillGrid1dTo2d(FillDirection dir, const GridPortalLink& grid1d, + GridPortalLink& grid2d); + + /// Index type for type earsed (**slow**) bin access + using IndexType = boost::container::static_vector; + + /// Helper function to get grid bin count in type-eraased way. + /// @return The number of bins in each direction + virtual IndexType numLocalBins() const = 0; + + /// Helper function to get grid bin content in type-eraased way. + /// @param indices The bin indices + /// @return The tracking volume at the bin + virtual TrackingVolume*& atLocalBins(IndexType indices) = 0; + + /// Helper function to get grid bin content in type-eraased way. + /// @param indices The bin indices + /// @return The tracking volume at the bin + virtual TrackingVolume* atLocalBins(IndexType indices) const = 0; + + private: + BinningValue m_direction; +}; + +/// Concrete class deriving from @c GridPortalLink that boxes a concrete grid for lookup. +/// @tparam Axes The axis types of the grid +template + requires(sizeof...(Axes) <= 2) +class GridPortalLinkT final : public GridPortalLink { + public: + /// The internal grid type + using GridType = Grid; + + /// The dimension of the grid + static constexpr std::size_t DIM = sizeof...(Axes); + + /// Constructor from a surface, axes and direction + /// @param surface The surface + /// @param direction The binning direction + /// @param axes The axes for the grid + /// @note The axes are checked for consistency with the bounds of @p surface. + GridPortalLinkT(std::shared_ptr surface, + BinningValue direction, Axes&&... axes) + : GridPortalLink(std::move(surface), direction), + m_grid(std::tuple{std::move(axes)...}) { + using enum BinningValue; + + if (const auto* cylinder = + dynamic_cast(m_surface.get())) { + checkConsistency(*cylinder); + + if (direction == binRPhi) { + m_projection = &projection; + } else if (direction == binZ) { + m_projection = &projection; + } else { + throw std::invalid_argument{"Invalid binning direction"}; + } + + } else if (const auto* disc = + dynamic_cast(m_surface.get())) { + checkConsistency(*disc); + + if (direction == binR) { + m_projection = &projection; + } else if (direction == BinningValue::binPhi) { + m_projection = &projection; + } else { + throw std::invalid_argument{"Invalid binning direction"}; + } + + } else { + throw std::logic_error{"Surface type is not supported"}; + } + } + + /// Get the grid + /// @return The grid + const GridType& grid() const override { return m_grid; } + + /// Get the grid + /// @return The grid + GridType& grid() { return m_grid; } + + /// Get the number of dimensions of the grid + /// @return The number of dimensions + unsigned int dim() const override { return DIM; } + + /// Prints an identification to the output stream + /// @param os The output stream + void toStream(std::ostream& os) const override { + os << "GridPortalLink"; + } + + /// Makes a 2D grid from a 1D grid by extending it using an optional axis + /// @param other The axis to use for the missing direction, + /// can be null for auto determination + /// @return A unique pointer to the 2D grid portal link + std::unique_ptr extendTo2d( + const IAxis* other) const override { + if constexpr (DIM == 2) { + return std::make_unique>(*this); + } else { + if (auto cylinder = + std::dynamic_pointer_cast(m_surface)) { + return extendTo2dImpl(cylinder, other); + } else if (auto disc = + std::dynamic_pointer_cast(m_surface)) { + return extendTo2dImpl(disc, other); + } else { + throw std::logic_error{ + "Surface type is not supported (this should not happen)"}; + } + } + } + + /// Set the volume on all grid bins + /// @param volume The volume to set + void setVolume(TrackingVolume* volume) override { + auto loc = m_grid.numLocalBins(); + if constexpr (GridType::DIM == 1) { + for (std::size_t i = 1; i <= loc[0]; i++) { + m_grid.atLocalBins({i}) = volume; + } + } else { + for (std::size_t i = 1; i <= loc[0]; i++) { + for (std::size_t j = 1; j <= loc[1]; j++) { + m_grid.atLocalBins({i, j}) = volume; + } + } + } + } + + /// Lookup the tracking volume at a position on the surface + /// @param gctx The geometry context + /// @param position The local position + /// @param tolerance The tolerance for the lookup + /// @note The position is required to be on the associated surface + /// @return The tracking volume (can be null) + Result resolveVolume( + const GeometryContext& gctx, const Vector3& position, + double tolerance = s_onSurfaceTolerance) const override { + auto res = m_surface->globalToLocal(gctx, position, tolerance); + if (!res.ok()) { + return res.error(); + } + + const Vector2& local = *res; + return resolveVolume(gctx, local, tolerance); + } + + /// Lookup the tracking volume at a position on the surface + /// @param position The local position + /// @note The position is required to be on the associated surface + /// @return The tracking volume (can be null) + Result resolveVolume( + const GeometryContext& /*gctx*/, const Vector2& position, + double /*tolerance*/ = s_onSurfaceTolerance) const override { + assert(surface().insideBounds(position, BoundaryTolerance::None())); + return m_grid.atPosition(m_projection(position)); + } + + protected: + /// Type erased access to the number of bins + /// @return The number of bins in each direction + IndexType numLocalBins() const override { + typename GridType::index_t idx = m_grid.numLocalBins(); + IndexType result; + for (std::size_t i = 0; i < DIM; i++) { + result.push_back(idx[i]); + } + return result; + } + + /// Type erased local bin access + /// @param indices The bin indices + /// @return The tracking volume at the bin + TrackingVolume*& atLocalBins(IndexType indices) override { + throw_assert(indices.size() == DIM, "Invalid number of indices"); + typename GridType::index_t idx; + for (std::size_t i = 0; i < DIM; i++) { + idx[i] = indices[i]; + } + return m_grid.atLocalBins(idx); + } + + /// Type erased local bin access + /// @param indices The bin indices + /// @return The tracking volume at the bin + TrackingVolume* atLocalBins(IndexType indices) const override { + throw_assert(indices.size() == DIM, "Invalid number of indices"); + typename GridType::index_t idx; + for (std::size_t i = 0; i < DIM; i++) { + idx[i] = indices[i]; + } + return m_grid.atLocalBins(idx); + } + + private: + /// Helper function that's assigned to project from the 2D local position to a + /// possible 1D grid. + template + static ActsVector projection(const Vector2& position) { + using enum BinningValue; + if constexpr (DIM == 2) { + return position; + } else { + if constexpr (std::is_same_v) { + static_assert(direction == binRPhi || direction == binZ, + "Invalid binning direction"); + + if constexpr (direction == binRPhi) { + return ActsVector<1>{position[0]}; + } else if constexpr (direction == binZ) { + return ActsVector<1>{position[1]}; + } + } else if constexpr (std::is_same_v) { + static_assert(direction == binR || direction == binPhi, + "Invalid binning direction"); + + if constexpr (direction == binR) { + return ActsVector<1>{position[0]}; + } else if constexpr (direction == binPhi) { + return ActsVector<1>{position[1]}; + } + } else if constexpr (std::is_same_v) { + static_assert(direction == binX || direction == binY, + "Invalid binning direction"); + + if constexpr (direction == binX) { + return ActsVector<1>{position[0]}; + } else if constexpr (direction == binY) { + return ActsVector<1>{position[1]}; + } + } + } + } + + GridType m_grid; + + /// Stores a function pointer that can project from 2D to 1D if needed + ActsVector (*m_projection)(const Vector2& position); +}; + +} // namespace Acts diff --git a/Core/include/Acts/Geometry/PortalError.hpp b/Core/include/Acts/Geometry/PortalError.hpp new file mode 100644 index 00000000000..0c704b1e9e5 --- /dev/null +++ b/Core/include/Acts/Geometry/PortalError.hpp @@ -0,0 +1,29 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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 +#include + +namespace Acts { + +enum class PortalError { + // ensure all values are non-zero + PositionNotOnAnyChildPortalLink = 1, +}; + +std::error_code make_error_code(Acts::PortalError e); + +} // namespace Acts + +namespace std { +// register with STL +template <> +struct is_error_code_enum : std::true_type {}; +} // namespace std diff --git a/Core/include/Acts/Geometry/PortalLinkBase.hpp b/Core/include/Acts/Geometry/PortalLinkBase.hpp new file mode 100644 index 00000000000..a445190af05 --- /dev/null +++ b/Core/include/Acts/Geometry/PortalLinkBase.hpp @@ -0,0 +1,119 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/Definitions/Tolerance.hpp" +#include "Acts/Utilities/BinningType.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "Acts/Utilities/Result.hpp" + +#include + +namespace Acts { + +class RegularSurface; +class TrackingVolume; +class GeometryContext; + +/// PortalLinkBase is the abstract base class for all portal links. +/// A portal link is a mapping between a surface and a point on the surface and +/// a destination tracking volum. +/// The derived classes implement different ways to resolve a volume +class PortalLinkBase { + protected: + /// Constructor from a surface. This constructor is only + /// called from derived classes + /// @param surface The surface + explicit PortalLinkBase(std::shared_ptr surface) + : m_surface(std::move(surface)) { + if (!m_surface) { + throw std::invalid_argument("Surface pointer must not be null"); + } + } + + public: + /// Virtual destructor in case the object is held as a derived + virtual ~PortalLinkBase() = default; + + /// Resolve a volume given a global position. Depending on the derived class, + /// the global position might be converted to a local position before lookup. + /// @param gctx The geometry context + /// @param position The global position + /// @param tolerance The tolerance for the lookup + /// + /// @return The tracking volume or null if no connection was found + virtual Result resolveVolume( + const GeometryContext& gctx, const Vector3& position, + double tolerance = s_onSurfaceTolerance) const = 0; + + /// Resolve a volume given a local position. The local position is assumed to + /// be on surface. + /// @param gctx The geometry context + /// @param position The local position + /// @param tolerance The tolerance for the lookup + /// + /// @return The tracking volume or null if no connection was found + virtual Result resolveVolume( + const GeometryContext& gctx, const Vector2& position, + double tolerance = s_onSurfaceTolerance) const = 0; + + //// Merge two portal link into a single one. The merge can resolve + /// combinations of difference derived classes, and will try to flatten and + /// deep merge given links if possible. + /// @param a The first portal link + /// @param b The second portal link + /// @param direction The binning direction in which to merge. Valid values are + /// depend on the surface types associated with the links. + /// @param logger The logger to use for messages + /// @return The merged portal link + static std::unique_ptr merge( + std::unique_ptr a, std::unique_ptr b, + BinningValue direction, const Logger& logger = getDummyLogger()); + + /// Stream output function + /// @param os The output stream + virtual void toStream(std::ostream& os) const = 0; + + /// Stream output operator + /// @param os The output stream + /// @param link The portal link + friend std::ostream& operator<<(std::ostream& os, + const PortalLinkBase& link) { + link.toStream(os); + return os; + } + + /// Getter for the associated surface + /// @return The surface + const RegularSurface& surface() const { return *m_surface; } + + /// Setter for the surface + /// @param surface The surface + void setSurface(std::shared_ptr surface) { + m_surface = std::move(surface); + } + + /// Getter for the underlying shared pointer + /// @return The shared pointer to the surface + const std::shared_ptr& surfacePtr() const { + return m_surface; + } + + protected: + /// Helper function to check a number of preconditions before merging is + /// executed. + static void checkMergePreconditions(const PortalLinkBase& a, + const PortalLinkBase& b, + BinningValue direction); + + std::shared_ptr m_surface; +}; + +} // namespace Acts diff --git a/Core/include/Acts/Geometry/TrivialPortalLink.hpp b/Core/include/Acts/Geometry/TrivialPortalLink.hpp new file mode 100644 index 00000000000..73199028feb --- /dev/null +++ b/Core/include/Acts/Geometry/TrivialPortalLink.hpp @@ -0,0 +1,65 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/Tolerance.hpp" +#include "Acts/Geometry/PortalLinkBase.hpp" + +namespace Acts { + +class TrackingVolume; +class GridPortalLink; + +/// Trivial portal link links to a single target volume on every point on a +/// surface +class TrivialPortalLink final : public PortalLinkBase { + public: + /// Construct a trivial portal link from a surface and a volume + /// @param surface is the surface + /// @param volume is the target + TrivialPortalLink(std::shared_ptr surface, + TrackingVolume& volume) + : PortalLinkBase(std::move(surface)), m_volume{&volume} {} + + /// Make a 1D grid portal link from this trivial portal link + /// The grid size is automatically determined from the surface bounds. + /// @param direction The binning direction + /// @return A grid + std::unique_ptr makeGrid(BinningValue direction) const; + + /// Print the portal link to a stream + /// @param os output stream + void toStream(std::ostream& os) const override; + + /// Resolve the volume for a 2D position + /// @note Always returns the single target volume + /// @param gctx is the geometry context + /// @param position is the 2D position + /// @param tolerance is the tolerance + /// @return The target volume (can be null) + Result resolveVolume( + const GeometryContext& gctx, const Vector2& position, + double tolerance = s_onSurfaceTolerance) const override; + + /// Resolve the volume for a 3D position + /// @note Always returns the single target volume + /// @param gctx is the geometry context + /// @param position is the 2D position + /// @param tolerance is the tolerance + /// @return The target volume (can be null) + /// @note The position is assumed to be on the associated surface. + Result resolveVolume( + const GeometryContext& gctx, const Vector3& position, + double tolerance = s_onSurfaceTolerance) const override; + + private: + TrackingVolume* m_volume; +}; + +} // namespace Acts diff --git a/Core/include/Acts/Utilities/IAxis.hpp b/Core/include/Acts/Utilities/IAxis.hpp index f61e3310d3b..62b241a9b19 100644 --- a/Core/include/Acts/Utilities/IAxis.hpp +++ b/Core/include/Acts/Utilities/IAxis.hpp @@ -75,6 +75,8 @@ class IAxis { case Variable: return callable( dynamic_cast&>(*this)); + default: + throw std::logic_error("Unknown axis type"); } }; @@ -86,6 +88,8 @@ class IAxis { return switchOnType(AxisBound); case Closed: return switchOnType(AxisClosed); + default: + throw std::logic_error("Unknown axis type"); } } diff --git a/Core/src/Geometry/CMakeLists.txt b/Core/src/Geometry/CMakeLists.txt index 0391a99b2fb..00d9618f062 100644 --- a/Core/src/Geometry/CMakeLists.txt +++ b/Core/src/Geometry/CMakeLists.txt @@ -35,4 +35,10 @@ target_sources( Volume.cpp VolumeBounds.cpp CylinderVolumeStack.cpp + GridPortalLink.cpp + GridPortalLinkMerging.cpp + TrivialPortalLink.cpp + CompositePortalLink.cpp + PortalLinkBase.cpp + PortalError.cpp ) diff --git a/Core/src/Geometry/CompositePortalLink.cpp b/Core/src/Geometry/CompositePortalLink.cpp new file mode 100644 index 00000000000..fa21653146f --- /dev/null +++ b/Core/src/Geometry/CompositePortalLink.cpp @@ -0,0 +1,128 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/. + +#include "Acts/Geometry/CompositePortalLink.hpp" + +#include "Acts/Geometry/PortalError.hpp" +#include "Acts/Surfaces/CylinderSurface.hpp" +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Surfaces/PlaneSurface.hpp" +#include "Acts/Surfaces/RegularSurface.hpp" + +#include +#include + +namespace Acts { + +CompositePortalLink::CompositePortalLink(std::unique_ptr a, + std::unique_ptr b, + BinningValue direction, bool flatten) + : PortalLinkBase(mergedSurface(a.get(), b.get(), direction)) { + if (!flatten) { + m_children.push_back(std::move(a)); + m_children.push_back(std::move(b)); + + } else { + auto handle = [&](std::unique_ptr link) { + if (auto* composite = dynamic_cast(link.get()); + composite != nullptr) { + m_children.insert( + m_children.end(), + std::make_move_iterator(composite->m_children.begin()), + std::make_move_iterator(composite->m_children.end())); + } else { + m_children.push_back(std::move(link)); + } + }; + + handle(std::move(a)); + handle(std::move(b)); + } +} + +Result CompositePortalLink::resolveVolume( + const GeometryContext& gctx, const Vector2& position, + double tolerance) const { + // In this overload, we have to go back to global, because the children have + // their own local coordinate systems + Vector3 global = m_surface->localToGlobal(gctx, position); + auto res = resolveVolume(gctx, global, tolerance); + if (!res.ok()) { + return res.error(); + } + return *res; +} + +Result CompositePortalLink::resolveVolume( + const GeometryContext& gctx, const Vector3& position, + double tolerance) const { + assert(m_surface->isOnSurface(gctx, position, BoundaryTolerance::None(), + tolerance)); + + for (const auto& child : m_children) { + if (child->surface().isOnSurface(gctx, position, BoundaryTolerance::None(), + tolerance)) { + return child->resolveVolume(gctx, position); + } + } + + return PortalError::PositionNotOnAnyChildPortalLink; +} + +std::shared_ptr CompositePortalLink::mergedSurface( + const PortalLinkBase* a, const PortalLinkBase* b, BinningValue direction) { + assert(a != nullptr); + assert(b != nullptr); + if (a->surface().type() != b->surface().type()) { + throw std::invalid_argument{"Cannot merge surfaces of different types"}; + } + + if (const auto* cylinderA = + dynamic_cast(&a->surface()); + cylinderA != nullptr) { + const auto& cylinderB = dynamic_cast(b->surface()); + + auto [merged, reversed] = cylinderA->mergedWith(cylinderB, direction, true); + return merged; + } else if (const auto* discA = + dynamic_cast(&a->surface()); + discA != nullptr) { + const auto& discB = dynamic_cast(b->surface()); + auto [merged, reversed] = discA->mergedWith(discB, direction, true); + return merged; + } else if (const auto* planeA = + dynamic_cast(&a->surface()); + planeA != nullptr) { + throw std::logic_error{"Plane surfaces not implemented yet"}; + } else { + throw std::invalid_argument{"Unsupported surface type"}; + } +} + +std::size_t CompositePortalLink::depth() const { + std::size_t result = 1; + + for (const auto& child : m_children) { + if (const auto* composite = + dynamic_cast(child.get())) { + result = std::max(result, composite->depth() + 1); + } + } + + return result; +} + +std::size_t CompositePortalLink::size() const { + return m_children.size(); +} + +void CompositePortalLink::toStream(std::ostream& os) const { + os << "CompositePortalLink"; +} + +} // namespace Acts diff --git a/Core/src/Geometry/GridPortalLink.cpp b/Core/src/Geometry/GridPortalLink.cpp new file mode 100644 index 00000000000..d7c11941374 --- /dev/null +++ b/Core/src/Geometry/GridPortalLink.cpp @@ -0,0 +1,416 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/. + +#include "Acts/Geometry/GridPortalLink.hpp" + +#include "Acts/Surfaces/PlaneSurface.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" + +namespace Acts { + +std::unique_ptr GridPortalLink::make( + const std::shared_ptr& surface, TrackingVolume& volume, + BinningValue direction) { + std::unique_ptr grid; + + if (const auto* cylinder = + dynamic_cast(surface.get()); + cylinder != nullptr) { + if (direction == BinningValue::binRPhi) { + ActsScalar r = cylinder->bounds().get(CylinderBounds::eR); + if (cylinder->bounds().coversFullAzimuth()) { + grid = GridPortalLink::make(surface, direction, + Axis{AxisClosed, -M_PI * r, M_PI * r, 1}); + } else { + ActsScalar hlPhi = + cylinder->bounds().get(CylinderBounds::eHalfPhiSector); + + grid = GridPortalLink::make(surface, direction, + Axis{AxisBound, -hlPhi * r, hlPhi * r, 1}); + } + } else if (direction == BinningValue::binZ) { + ActsScalar hlZ = cylinder->bounds().get(CylinderBounds::eHalfLengthZ); + grid = GridPortalLink::make(surface, direction, + Axis{AxisBound, -hlZ, hlZ, 1}); + } else { + throw std::invalid_argument{"Invalid binning direction"}; + } + } else if (const auto* disc = dynamic_cast(surface.get()); + disc != nullptr) { + const auto& bounds = dynamic_cast(disc->bounds()); + if (direction == BinningValue::binR) { + ActsScalar minR = bounds.get(RadialBounds::eMinR); + ActsScalar maxR = bounds.get(RadialBounds::eMaxR); + grid = GridPortalLink::make(surface, direction, + Axis{AxisBound, minR, maxR, 1}); + } else if (direction == BinningValue::binPhi) { + if (bounds.coversFullAzimuth()) { + grid = GridPortalLink::make(surface, direction, + Axis{AxisClosed, -M_PI, M_PI, 1}); + } else { + ActsScalar hlPhi = bounds.get(RadialBounds::eHalfPhiSector); + grid = GridPortalLink::make(surface, direction, + Axis{AxisBound, -hlPhi, hlPhi, 1}); + } + } else { + throw std::invalid_argument{"Invalid binning direction"}; + } + } else if (const auto* plane = + dynamic_cast(surface.get()); + plane != nullptr) { + throw std::invalid_argument{"Plane surface is not implemented yet"}; + + } else { + throw std::invalid_argument{"Surface type is not supported"}; + } + + assert(grid != nullptr); + grid->setVolume(&volume); + + return grid; +} + +void GridPortalLink::checkConsistency(const CylinderSurface& cyl) const { + if (cyl.bounds().get(CylinderBounds::eAveragePhi) != 0) { + throw std::invalid_argument( + "GridPortalLink: CylinderBounds: only average phi == 0 is " + "supported. Rotate the cylinder surface."); + } + + constexpr auto tolerance = s_onSurfaceTolerance; + auto same = [](auto a, auto b) { return std::abs(a - b) < tolerance; }; + + auto checkZ = [&cyl, same](const IAxis& axis) { + ActsScalar hlZ = cyl.bounds().get(CylinderBounds::eHalfLengthZ); + if (!same(axis.getMin(), -hlZ) || !same(axis.getMax(), hlZ)) { + throw std::invalid_argument( + "GridPortalLink: CylinderBounds: invalid length setup."); + } + }; + auto checkRPhi = [&cyl, same](const IAxis& axis) { + ActsScalar hlPhi = cyl.bounds().get(CylinderBounds::eHalfPhiSector); + ActsScalar r = cyl.bounds().get(CylinderBounds::eR); + if (ActsScalar hlRPhi = r * hlPhi; + !same(axis.getMin(), -hlRPhi) || !same(axis.getMax(), hlRPhi)) { + throw std::invalid_argument( + "GridPortalLink: CylinderBounds: invalid phi sector setup: axes " + "don't match bounds"); + } + + // If full cylinder, make sure axis wraps around + if (same(hlPhi, M_PI)) { + if (axis.getBoundaryType() != AxisBoundaryType::Closed) { + throw std::invalid_argument( + "GridPortalLink: CylinderBounds: invalid phi sector setup: " + "axis is not closed."); + } + } else { + if (axis.getBoundaryType() == AxisBoundaryType::Closed) { + throw std::invalid_argument( + "GridPortalLink: CylinderBounds: invalid phi sector setup: " + "axis is closed."); + } + } + }; + + if (dim() == 1) { + const IAxis& axisLoc0 = *grid().axes().front(); + if (direction() == BinningValue::binRPhi) { + checkRPhi(axisLoc0); + } else { + checkZ(axisLoc0); + } + } else { // DIM == 2 + const auto& axisLoc0 = *grid().axes().front(); + const auto& axisLoc1 = *grid().axes().back(); + checkRPhi(axisLoc0); + checkZ(axisLoc1); + } +} + +void GridPortalLink::checkConsistency(const DiscSurface& disc) const { + constexpr auto tolerance = s_onSurfaceTolerance; + auto same = [](auto a, auto b) { return std::abs(a - b) < tolerance; }; + + const auto* bounds = dynamic_cast(&disc.bounds()); + if (bounds == nullptr) { + throw std::invalid_argument( + "GridPortalLink: DiscBounds: invalid bounds type."); + } + + if (bounds->get(RadialBounds::eAveragePhi) != 0) { + throw std::invalid_argument( + "GridPortalLink: DiscBounds: only average phi == 0 is supported. " + "Rotate the disc surface."); + } + + auto checkR = [&bounds, same](const IAxis& axis) { + ActsScalar minR = bounds->get(RadialBounds::eMinR); + ActsScalar maxR = bounds->get(RadialBounds::eMaxR); + if (!same(axis.getMin(), minR) || !same(axis.getMax(), maxR)) { + throw std::invalid_argument( + "GridPortalLink: DiscBounds: invalid radius setup."); + } + }; + + auto checkPhi = [&bounds, same](const IAxis& axis) { + ActsScalar hlPhi = bounds->get(RadialBounds::eHalfPhiSector); + if (!same(axis.getMin(), -hlPhi) || !same(axis.getMax(), hlPhi)) { + throw std::invalid_argument( + "GridPortalLink: DiscBounds: invalid phi sector setup."); + } + // If full disc, make sure axis wraps around + if (same(hlPhi, M_PI)) { + if (axis.getBoundaryType() != Acts::AxisBoundaryType::Closed) { + throw std::invalid_argument( + "GridPortalLink: DiscBounds: invalid phi sector setup: axis is " + "not closed."); + } + } else { + if (axis.getBoundaryType() == Acts::AxisBoundaryType::Closed) { + throw std::invalid_argument( + "GridPortalLink: DiscBounds: invalid phi sector setup: axis " + "is closed."); + } + } + }; + + if (dim() == 1) { + const IAxis& axisLoc0 = *grid().axes().front(); + if (direction() == BinningValue::binR) { + checkR(axisLoc0); + } else { + checkPhi(axisLoc0); + } + } else { // DIM == 2 + const auto& axisLoc0 = *grid().axes().front(); + const auto& axisLoc1 = *grid().axes().back(); + checkR(axisLoc0); + checkPhi(axisLoc1); + } +} + +void GridPortalLink::printContents(std::ostream& os) const { + std::size_t dim = grid().axes().size(); + os << "----- GRID " << dim << "d -----" << std::endl; + os << grid() << " along " << direction() << std::endl; + + auto lpad = [](const std::string& s, std::size_t n) { + return s.size() < n ? std::string(n - s.size(), ' ') + s : s; + }; + + auto rpad = [](const std::string& s, std::size_t n) { + return s.size() < n ? s + std::string(n - s.size(), ' ') : s; + }; + + std::string loc0; + std::string loc1; + + bool flipped = false; + + if (surface().type() == Surface::Cylinder) { + loc0 = "rPhi"; + loc1 = "z"; + flipped = direction() != BinningValue::binRPhi; + } else if (surface().type() == Surface::Disc) { + loc0 = "r"; + loc1 = "phi"; + flipped = direction() != BinningValue::binR; + } else if (surface().type() == Surface::Plane) { + loc0 = "x"; + loc1 = "y"; + flipped = direction() != BinningValue::binX; + } else { + throw std::invalid_argument{"Unsupported surface type"}; + } + + if (dim == 1) { + auto loc = numLocalBins(); + + if (flipped) { + os << lpad(loc1, 4) << " > " << lpad("i=0", 10) << " "; + for (std::size_t i = 1; i <= loc.at(0) + 1; i++) { + os << lpad("i=" + std::to_string(i), 13) + " "; + } + os << std::endl; + + os << std::string(4, ' '); + for (std::size_t i = 0; i <= loc.at(0) + 1; i++) { + std::string name = "0x0"; + if (const auto* v = atLocalBins({i}); v != nullptr) { + name = v->volumeName(); + } + name = name.substr(0, std::min(name.size(), std::size_t{13})); + os << lpad(name, 13) << " "; + } + os << std::endl; + + } else { + os << "v " << loc0 << std::endl; + for (std::size_t i = 0; i <= loc.at(0) + 1; i++) { + os << "i=" << i << " "; + std::string name = "0x0"; + if (const auto* v = atLocalBins({i}); v != nullptr) { + name = v->volumeName(); + } + name = name.substr(0, std::min(name.size(), std::size_t{13})); + os << lpad(name, 13) << " "; + os << std::endl; + } + } + + } else { + auto loc = numLocalBins(); + os << rpad("v " + loc0 + "|" + loc1 + " >", 14) + "j=0 "; + for (std::size_t j = 1; j <= loc.at(1) + 1; j++) { + os << lpad("j=" + std::to_string(j), 13) << " "; + } + os << std::endl; + for (std::size_t i = 0; i <= loc.at(0) + 1; i++) { + os << "i=" << i << " "; + for (std::size_t j = 0; j <= loc.at(1) + 1; j++) { + std::string name = "0x0"; + if (const auto* v = atLocalBins({i, j}); v != nullptr) { + name = v->volumeName(); + } + name = name.substr(0, std::min(name.size(), std::size_t{13})); + os << lpad(name, 13) << " "; + } + os << std::endl; + } + } +} + +void GridPortalLink::fillGrid1dTo2d(FillDirection dir, + const GridPortalLink& grid1d, + GridPortalLink& grid2d) { + const auto locSource = grid1d.numLocalBins(); + const auto locDest = grid2d.numLocalBins(); + assert(locSource.size() == 1); + assert(locDest.size() == 2); + + for (std::size_t i = 0; i <= locSource[0] + 1; ++i) { + TrackingVolume* source = grid1d.atLocalBins({i}); + + if (dir == FillDirection::loc1) { + for (std::size_t j = 0; j <= locDest[1] + 1; ++j) { + grid2d.atLocalBins({i, j}) = source; + } + } else if (dir == FillDirection::loc0) { + for (std::size_t j = 0; j <= locDest[0] + 1; ++j) { + grid2d.atLocalBins({j, i}) = source; + } + } + } +} + +std::unique_ptr GridPortalLink::extendTo2dImpl( + const std::shared_ptr& surface, const IAxis* other) const { + assert(dim() == 1); + if (direction() == BinningValue::binRPhi) { + const auto& axisRPhi = *grid().axes().front(); + // 1D direction is binRPhi, so add a Z axis + ActsScalar hlZ = surface->bounds().get(CylinderBounds::eHalfLengthZ); + + auto grid = axisRPhi.visit([&](const auto& axis0) { + Axis axisZ{AxisBound, -hlZ, hlZ, 1}; + const IAxis* axis = other != nullptr ? other : &axisZ; + return axis->visit( + [&surface, + &axis0](const auto& axis1) -> std::unique_ptr { + return GridPortalLink::make(surface, axis0, axis1); + }); + }); + + fillGrid1dTo2d(FillDirection::loc1, *this, *grid); + return grid; + + } else { + const auto& axisZ = *grid().axes().front(); + // 1D direction is binZ, so add an rPhi axis + ActsScalar r = surface->bounds().get(CylinderBounds::eR); + ActsScalar hlPhi = surface->bounds().get(CylinderBounds::eHalfPhiSector); + ActsScalar hlRPhi = r * hlPhi; + + auto makeGrid = [&](auto bdt) { + auto grid = axisZ.visit([&](const auto& axis1) { + Axis axisRPhi{bdt, -hlRPhi, hlRPhi, 1}; + const IAxis* axis = other != nullptr ? other : &axisRPhi; + return axis->visit( + [&](const auto& axis0) -> std::unique_ptr { + return GridPortalLink::make(surface, axis0, axis1); + }); + }); + + fillGrid1dTo2d(FillDirection::loc0, *this, *grid); + return grid; + }; + + if (surface->bounds().coversFullAzimuth()) { + return makeGrid(AxisClosed); + } else { + return makeGrid(AxisBound); + } + } +} + +std::unique_ptr GridPortalLink::extendTo2dImpl( + const std::shared_ptr& surface, const IAxis* other) const { + assert(dim() == 1); + + const auto* bounds = dynamic_cast(&surface->bounds()); + if (bounds == nullptr) { + throw std::invalid_argument( + "GridPortalLink: DiscBounds: invalid bounds type."); + } + + if (direction() == BinningValue::binR) { + const auto& axisR = *grid().axes().front(); + // 1D direction is binR, so add a phi axis + ActsScalar hlPhi = bounds->get(RadialBounds::eHalfPhiSector); + + auto makeGrid = [&](auto bdt) { + auto grid = axisR.visit([&](const auto& axis0) { + Axis axisPhi{bdt, -hlPhi, hlPhi, 1}; + const IAxis* axis = other != nullptr ? other : &axisPhi; + return axis->visit( + [&](const auto& axis1) -> std::unique_ptr { + return GridPortalLink::make(surface, axis0, axis1); + }); + }); + + fillGrid1dTo2d(FillDirection::loc1, *this, *grid); + return grid; + }; + + if (bounds->coversFullAzimuth()) { + return makeGrid(AxisClosed); + } else { + return makeGrid(AxisBound); + } + } else { + const auto& axisPhi = *grid().axes().front(); + // 1D direction is binPhi, so add an R axis + ActsScalar rMin = bounds->get(RadialBounds::eMinR); + ActsScalar rMax = bounds->get(RadialBounds::eMaxR); + + auto grid = axisPhi.visit([&](const auto& axis1) { + Axis axisR{AxisBound, rMin, rMax, 1}; + const IAxis* axis = other != nullptr ? other : &axisR; + return axis->visit( + [&surface, + &axis1](const auto& axis0) -> std::unique_ptr { + return GridPortalLink::make(surface, axis0, axis1); + }); + }); + fillGrid1dTo2d(FillDirection::loc0, *this, *grid); + return grid; + } +} + +} // namespace Acts diff --git a/Core/src/Geometry/GridPortalLinkMerging.cpp b/Core/src/Geometry/GridPortalLinkMerging.cpp new file mode 100644 index 00000000000..51dfd89548b --- /dev/null +++ b/Core/src/Geometry/GridPortalLinkMerging.cpp @@ -0,0 +1,572 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/. + +#include "Acts/Definitions/Tolerance.hpp" +#include "Acts/Geometry/CompositePortalLink.hpp" +#include "Acts/Geometry/GridPortalLink.hpp" +#include "Acts/Surfaces/CylinderSurface.hpp" +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Surfaces/PlaneSurface.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" + +#include +#include +#include +#include +#include + +namespace Acts { + +namespace { + +template +std::unique_ptr makeGrid( + const std::shared_ptr& surface, BinningValue direction, + const Logger& logger, std::tuple args, const IAxis* otherAxis, + bool prepend) { + // @TODO: PlaneSurface support + + ACTS_VERBOSE("Make resulting merged grid"); + + // This is to make it possible to construct Axis with the tuple arguments + auto axisFactory = [](Ts&&... axisArgs) { + return Axis{std::forward(axisArgs)...}; + }; + + // Avoid copy-pasting identical code twice below + auto makeGrid = [&](auto boundaryType) -> std::unique_ptr { + auto axisArgs = std::tuple_cat(std::tuple{boundaryType}, std::move(args)); + auto merged = std::apply(axisFactory, std::move(axisArgs)); + + if (otherAxis == nullptr) { + // No other axis + ACTS_VERBOSE(" ~> merged axis: " << merged); + return GridPortalLink::make(surface, direction, std::move(merged)); + } else { + return otherAxis->visit( + [&](const auto& axis) -> std::unique_ptr { + if (prepend) { + ACTS_VERBOSE(" ~> other axis (prepend): " << axis); + ACTS_VERBOSE(" ~> merged axis: " << merged); + return GridPortalLink::make(surface, axis, std::move(merged)); + } else { + ACTS_VERBOSE(" ~> merged axis: " << merged); + ACTS_VERBOSE(" ~> other axis (append): " << axis); + return GridPortalLink::make(surface, std::move(merged), axis); + } + }); + } + }; + + // Check if we're in the cylinder or disc case, and the resulting bounds wrap + // around and should have closed binning + if (direction == BinningValue::binPhi || direction == BinningValue::binRPhi) { + if (const auto* cylinder = + dynamic_cast(surface.get()); + cylinder != nullptr) { + if (cylinder->bounds().coversFullAzimuth()) { + return makeGrid(AxisClosed); + } + } else if (const auto* disc = + dynamic_cast(surface.get()); + disc != nullptr) { + const auto& radialBounds = + dynamic_cast(disc->bounds()); + if (radialBounds.coversFullAzimuth()) { + return makeGrid(AxisClosed); + } + } + } + + return makeGrid(AxisBound); +} + +std::unique_ptr mergeVariable( + const std::shared_ptr& mergedSurface, const IAxis& axisA, + const IAxis& axisB, ActsScalar /*tolerance*/, BinningValue direction, + const Logger& logger, const IAxis* otherAxis, bool prepend) { + ACTS_VERBOSE("Variable merge: direction is " << direction); + + ACTS_VERBOSE("~> axis a: " << axisA); + ACTS_VERBOSE("~> axis b: " << axisB); + + std::vector binEdges; + + binEdges.reserve(axisA.getNBins() + axisB.getNBins() + 1); + + auto edgesA = axisA.getBinEdges(); + + if (direction == BinningValue::binR) { + ACTS_VERBOSE("Performing asymmetric merge"); + std::ranges::copy(edgesA, std::back_inserter(binEdges)); + + } else { + ACTS_VERBOSE("Performing symmetrized merge"); + ActsScalar halfWidth = + (axisA.getMax() - axisA.getMin() + axisB.getMax() - axisB.getMin()) / + 2.0; + ACTS_VERBOSE(" ~> half width: " << halfWidth); + + ActsScalar shift = axisA.getMax() - halfWidth; + ACTS_VERBOSE(" ~> shift: " << shift); + + std::ranges::transform(edgesA, std::back_inserter(binEdges), + [&](ActsScalar edge) { return edge + shift; }); + } + + ActsScalar stitchPoint = binEdges.back(); + auto edgesB = axisB.getBinEdges(); + std::transform( + std::next(edgesB.begin()), edgesB.end(), std::back_inserter(binEdges), + [&](ActsScalar edge) { return edge - axisB.getMin() + stitchPoint; }); + + return makeGrid(mergedSurface, direction, logger, + std::tuple{std::move(binEdges)}, otherAxis, prepend); +} + +std::unique_ptr mergeEquidistant( + const std::shared_ptr& mergedSurface, const IAxis& axisA, + const IAxis& axisB, ActsScalar tolerance, BinningValue direction, + const Logger& logger, const IAxis* otherAxis, bool prepend) { + ACTS_VERBOSE("===> potentially equidistant merge: checking bin widths"); + + ACTS_VERBOSE("~> axis a: " << axisA); + ACTS_VERBOSE("~> axis b: " << axisB); + + ActsScalar binsWidthA = (axisA.getMax() - axisA.getMin()) / + static_cast(axisA.getNBins()); + ActsScalar binsWidthB = (axisB.getMax() - axisB.getMin()) / + static_cast(axisB.getNBins()); + + ACTS_VERBOSE(" ~> binWidths: " << binsWidthA << " vs " << binsWidthB); + + if (std::abs(binsWidthA - binsWidthB) < tolerance) { + ACTS_VERBOSE("==> binWidths same: " << binsWidthA); + + ActsScalar min = std::numeric_limits::signaling_NaN(); + ActsScalar max = std::numeric_limits::signaling_NaN(); + + if (direction == BinningValue::binR) { + ACTS_VERBOSE("Performing asymmetric merge"); + min = axisA.getMin(); + max = axisB.getMax(); + } else { + ACTS_VERBOSE("Performing symmetrized merge"); + + ActsScalar halfWidth = + (axisA.getMax() - axisA.getMin() + axisB.getMax() - axisB.getMin()) / + 2.0; + + min = -halfWidth; + max = halfWidth; + } + + return makeGrid(mergedSurface, direction, logger, + std::tuple{min, max, axisA.getNBins() + axisB.getNBins()}, + otherAxis, prepend); + + } else { + ACTS_VERBOSE("==> binWidths differ: " << binsWidthA << " vs " << binsWidthB + << " ~> variable merge"); + + std::unique_ptr mergedPortalLink = + mergeVariable(mergedSurface, axisA, axisB, tolerance, direction, logger, + otherAxis, prepend); + return mergedPortalLink; + } +} + +std::unique_ptr colinearMerge( + const std::shared_ptr& mergedSurface, const IAxis& axisA, + const IAxis& axisB, ActsScalar tolerance, BinningValue direction, + const Logger& logger, const IAxis* otherAxis, bool prepend) { + AxisType aType = axisA.getType(); + AxisType bType = axisB.getType(); + if (axisA.getBoundaryType() != axisB.getBoundaryType()) { + ACTS_WARNING("AxisBoundaryTypes are different"); + return nullptr; + } + + if (axisA.getBoundaryType() != AxisBoundaryType::Bound) { + ACTS_WARNING("AxisBoundaryType is not Bound, cannot do colinear merge"); + return nullptr; + } + + // convenience only + auto mergeVariableLocal = [&] { + return mergeVariable(mergedSurface, axisA, axisB, tolerance, direction, + logger, otherAxis, prepend); + }; + + if (aType == AxisType::Equidistant && bType == AxisType::Equidistant) { + auto mergedPortalLink = + mergeEquidistant(mergedSurface, axisA, axisB, tolerance, direction, + logger, otherAxis, prepend); + return mergedPortalLink; + } else if (aType == AxisType::Variable && bType == AxisType::Variable) { + ACTS_VERBOSE("===> variable merge"); + auto mergedPortalLink = mergeVariableLocal(); + return mergedPortalLink; + } else if (aType == AxisType::Equidistant && bType == AxisType::Variable) { + ACTS_WARNING("===> mixed merged"); + auto mergedPortalLink = mergeVariableLocal(); + return mergedPortalLink; + } else { + ACTS_WARNING("===> mixed merged"); + auto mergedPortalLink = mergeVariableLocal(); + return mergedPortalLink; + } +} + +std::unique_ptr mergeGridPortals( + const GridPortalLink* a, const GridPortalLink* b, + const RegularSurface* surfaceA, const RegularSurface* surfaceB, + BinningValue loc0, BinningValue loc1, BinningValue direction, + const Logger& logger) { + assert(surfaceA != nullptr); + assert(surfaceB != nullptr); + assert(surfaceA->type() == surfaceB->type()); + assert(a->dim() == 2 || a->dim() == 1); + assert(a->dim() == b->dim()); + + ACTS_VERBOSE(" - a: " << a->surface().bounds()); + ACTS_VERBOSE(" - b: " << b->surface().bounds()); + + // This tolerance is used checking bin equivalence. It's not intended to be + // user configurable, as we don't foresee cases where the tolerance would have + // to be adjusteed to the input geometry. + constexpr auto tolerance = s_onSurfaceTolerance; + + std::shared_ptr mergedSurface = nullptr; + bool reversed = false; + + if (const auto* cylinderA = dynamic_cast(surfaceA); + cylinderA != nullptr) { + std::tie(mergedSurface, reversed) = + cylinderA->mergedWith(dynamic_cast(*surfaceB), + direction, true, logger); + } else if (const auto* discA = dynamic_cast(surfaceA); + discA != nullptr) { + std::tie(mergedSurface, reversed) = discA->mergedWith( + dynamic_cast(*surfaceB), direction, true, logger); + } else if (const auto* planeA = dynamic_cast(surfaceA); + planeA != nullptr) { + throw std::logic_error{"Plane surfaces not implemented yet"}; + } else { + throw std::invalid_argument{"Unsupported surface type"}; + } + + ACTS_VERBOSE("Merged surface: " << mergedSurface->bounds()); + ACTS_VERBOSE("~> reversed? " << std::boolalpha << reversed); + + // Normalize ordering of grid portals and surfaces: a is always at lower + // range than b + if (reversed) { + ACTS_VERBOSE("Swapping grid portals and surfaces after merging"); + std::swap(surfaceA, surfaceB); + std::swap(a, b); + } + + // We do this here, because this is the only place where we've normalized + // the ordering of grids just above: a is always "before" b, so we can + // iterate over the bins in a unified way to fill the merged grid. + auto fillGrid = [&](std::unique_ptr merged) { + if (merged != nullptr) { + ACTS_VERBOSE("Post processing merged grid: " << merged->grid()); + GridPortalLink::fillMergedGrid(*a, *b, *merged, direction, logger); + } + return merged; + }; + + if (a->dim() == 1) { + ACTS_VERBOSE("Merge two 1D GridPortalLinks on " << surfaceA->name() + << "s in " << direction); + if (a->direction() != b->direction()) { + ACTS_VERBOSE("GridPortalLinks have different directions"); + ACTS_VERBOSE("=> cross merge"); + + // Find the one which is already binned in the merging direction + const auto* aligned = a->direction() == direction ? a : b; + const auto* alignedSurface = + a->direction() == direction ? surfaceA : surfaceB; + const auto* other = a->direction() == direction ? b : a; + const auto* otherSurface = + a->direction() == direction ? surfaceB : surfaceA; + + // Extend the aligned one by the other one's axis + auto aligned2D = aligned->extendTo2d(other->grid().axes().front()); + // Edtend the other one with a single bin + auto other2D = other->extendTo2d(nullptr); + + assert(aligned2D != nullptr); + assert(other2D != nullptr); + + ACTS_VERBOSE("Expanded grids:"); + ACTS_VERBOSE(" - aligned: " << aligned2D->grid()); + ACTS_VERBOSE(" - other: " << other2D->grid()); + + // @Fix ordering (a,b) should already be in good order, to save us one additional roundtrip + if (a->direction() == direction) { + return mergeGridPortals(aligned2D.get(), other2D.get(), alignedSurface, + otherSurface, loc0, loc1, direction, logger); + } else { + return mergeGridPortals(other2D.get(), aligned2D.get(), otherSurface, + alignedSurface, loc0, loc1, direction, logger); + } + } + + const auto& axisA = *a->grid().axes().front(); + const auto& axisB = *b->grid().axes().front(); + + if (direction == loc0) { + ACTS_VERBOSE("Grids are binned along " << a->direction()); + if (a->direction() == loc0) { + ACTS_VERBOSE("=> colinear merge"); + + return fillGrid(colinearMerge(mergedSurface, axisA, axisB, tolerance, + loc0, logger, nullptr, false)); + + } else { + ACTS_VERBOSE("=> parallel merge"); + + auto a2D = a->extendTo2d(nullptr); + auto b2D = b->extendTo2d(nullptr); + assert(a2D != nullptr); + assert(b2D != nullptr); + + ACTS_VERBOSE("Expanded grids:"); + ACTS_VERBOSE(" - a: " << a2D->grid()); + ACTS_VERBOSE(" - b: " << b2D->grid()); + + return mergeGridPortals(a2D.get(), b2D.get(), surfaceA, surfaceB, loc0, + loc1, direction, logger); + } + } else if (direction == loc1) { + ACTS_VERBOSE("Grids are binned along " << a->direction()); + if (a->direction() == loc1) { + ACTS_VERBOSE("=> colinear merge"); + + return fillGrid(colinearMerge(mergedSurface, axisA, axisB, tolerance, + loc1, logger, nullptr, false)); + + } else { + ACTS_VERBOSE("=> parallel merge"); + auto a2D = a->extendTo2d(nullptr); + auto b2D = b->extendTo2d(nullptr); + assert(a2D != nullptr); + assert(b2D != nullptr); + + ACTS_VERBOSE("Expanded grids:"); + ACTS_VERBOSE(" - a: " << a2D->grid()); + ACTS_VERBOSE(" - b: " << b2D->grid()); + return mergeGridPortals(a2D.get(), b2D.get(), surfaceA, surfaceB, loc0, + loc1, direction, logger); + } + + } else { + ACTS_ERROR("Invalid binning direction: " << direction); + throw std::invalid_argument{"Invalid binning direction"}; + } + } else { + ACTS_VERBOSE("Merging two 2D GridPortalLinks on " + << surfaceA->name() << "s in " << direction << " direction"); + + const auto& loc0AxisA = *a->grid().axes().front(); + const auto& loc1AxisA = *a->grid().axes().back(); + const auto& loc0AxisB = *b->grid().axes().front(); + const auto& loc1AxisB = *b->grid().axes().back(); + + if (direction == loc0) { + ACTS_VERBOSE("=> colinear merge along " << loc0); + ACTS_VERBOSE("--> Checking if " << loc1 << " axes are identical"); + + if (loc1AxisA != loc1AxisB) { + ACTS_WARNING(" ~> " + << loc1 + << " axes are not identical, falling back to composite " + "merging"); + return nullptr; + } + ACTS_VERBOSE(" ~> they are!"); + + ACTS_VERBOSE(" ~> " << loc1 << " axis: " << loc1AxisA); + return fillGrid(colinearMerge(mergedSurface, loc0AxisA, loc0AxisB, + tolerance, loc0, logger, &loc1AxisA, + false)); + + } else if (direction == loc1) { + ACTS_VERBOSE("=> colinear merge along " << loc1); + ACTS_VERBOSE("--> Checking if " << loc0 << " axes are identical"); + + if (loc0AxisA != loc0AxisB) { + ACTS_WARNING(" ~> " + << loc0 + << " axes are not identical, falling back to composite " + "merging"); + return nullptr; + } + ACTS_VERBOSE(" ~> they are!"); + + ACTS_VERBOSE(" ~> " << loc0 << " axis: " << loc0AxisA); + return fillGrid(colinearMerge(mergedSurface, loc1AxisA, loc1AxisB, + tolerance, loc1, logger, &loc0AxisA, true)); + + } else { + ACTS_ERROR("Invalid binning direction: " << a->direction()); + throw std::invalid_argument{"Invalid binning direction"}; + } + } +} + +std::unique_ptr mergeGridPortals(const GridPortalLink* a, + const GridPortalLink* b, + BinningValue direction, + const Logger& logger) { + using enum BinningValue; + assert(a->dim() == 2 || a->dim() == 1); + assert(b->dim() == 2 || b->dim() == 1); + + if (a->dim() < b->dim()) { + return mergeGridPortals(b, a, direction, logger); + } + + ACTS_VERBOSE("Merging GridPortalLinks along " << direction << ":"); + ACTS_VERBOSE(" - a: " << a->grid() << " along: " << a->direction()); + ACTS_VERBOSE(" - b: " << b->grid() << " along: " << b->direction()); + + const auto* cylinder = dynamic_cast(&a->surface()); + const auto* disc = dynamic_cast(&a->surface()); + + if (a->dim() == b->dim()) { + ACTS_VERBOSE("Grid both have same dimension: " << a->dim()); + + if (cylinder != nullptr) { + return mergeGridPortals( + a, b, cylinder, &dynamic_cast(b->surface()), + binRPhi, binZ, direction, logger); + } else if (disc != nullptr) { + return mergeGridPortals(a, b, disc, + &dynamic_cast(b->surface()), + binR, binPhi, direction, logger); + } else { + // @TODO: Support PlaneSurface + ACTS_VERBOSE("Surface type is not supported here, falling back"); + return nullptr; + } + } else { + ACTS_VERBOSE("Grids have different dimension, extending rhs to 2D"); + const IAxis* otherAxis = nullptr; + + if (b->direction() == direction) { + ACTS_VERBOSE("1D grid is binned in merging direction " << direction); + ACTS_VERBOSE("~> Adding complementary axis"); + + if (cylinder != nullptr) { + otherAxis = direction == binRPhi ? a->grid().axes().back() + : a->grid().axes().front(); + } else if (disc != nullptr) { + otherAxis = direction == binR ? a->grid().axes().back() + : a->grid().axes().front(); + } else { + ACTS_VERBOSE("Surface type is not supported here, falling back"); + return nullptr; + } + } else { + ACTS_VERBOSE("1D grid is binned in complementary direction"); + } + + auto b2D = b->extendTo2d(otherAxis); + ACTS_VERBOSE("-> new grid: " << b2D->grid()); + return mergeGridPortals(a, b2D.get(), direction, logger); + } +} + +} // namespace + +void GridPortalLink::fillMergedGrid(const GridPortalLink& a, + const GridPortalLink& b, + GridPortalLink& merged, + BinningValue direction, + const Logger& logger) { + ACTS_VERBOSE("Filling merged grid"); + assert(a.dim() == b.dim()); + assert(a.direction() == b.direction()); + const auto locBinsA = a.numLocalBins(); + const auto locBinsB = b.numLocalBins(); + + ACTS_VERBOSE("a: " << a.grid()); + ACTS_VERBOSE("b: " << b.grid()); + ACTS_VERBOSE("merged: " << merged.grid()); + + if (a.dim() == 1) { + std::size_t nBinsA = locBinsA.at(0); + std::size_t nBinsB = locBinsB.at(0); + + for (std::size_t i = 1; i <= nBinsA; ++i) { + merged.atLocalBins({i}) = a.atLocalBins({i}); + } + for (std::size_t i = 1; i <= nBinsB; ++i) { + merged.atLocalBins({nBinsA + i}) = b.atLocalBins({i}); + } + } else { + if (a.direction() == direction) { + std::size_t nBinsA = locBinsA.at(0); + std::size_t nBinsB = locBinsB.at(0); + std::size_t nBinsCommon = locBinsB.at(1); + + for (std::size_t i = 1; i <= nBinsA; ++i) { + for (std::size_t j = 1; j <= nBinsCommon; ++j) { + merged.atLocalBins({i, j}) = a.atLocalBins({i, j}); + } + } + + for (std::size_t i = 1; i <= nBinsB; ++i) { + for (std::size_t j = 1; j <= nBinsCommon; ++j) { + std::size_t ti = i + nBinsA; + merged.atLocalBins({ti, j}) = b.atLocalBins({i, j}); + } + } + } else { + std::size_t nBinsA = locBinsA.at(1); + std::size_t nBinsB = locBinsB.at(1); + std::size_t nBinsCommon = locBinsB.at(0); + + for (std::size_t i = 1; i <= nBinsCommon; ++i) { + for (std::size_t j = 1; j <= nBinsA; ++j) { + merged.atLocalBins({i, j}) = a.atLocalBins({i, j}); + } + } + + for (std::size_t i = 1; i <= nBinsCommon; ++i) { + for (std::size_t j = 1; j <= nBinsB; ++j) { + std::size_t tj = j + nBinsA; + merged.atLocalBins({i, tj}) = b.atLocalBins({i, j}); + } + } + } + } +} + +std::unique_ptr GridPortalLink::merge(const GridPortalLink& a, + const GridPortalLink& b, + BinningValue direction, + const Logger& logger) { + ACTS_VERBOSE("Merging two GridPortalLinks"); + + checkMergePreconditions(a, b, direction); + + auto merged = mergeGridPortals(&a, &b, direction, logger); + if (merged == nullptr) { + ACTS_WARNING("Grid merging failed returning nullptr"); + return nullptr; + } + + return merged; +} + +} // namespace Acts diff --git a/Core/src/Geometry/PortalError.cpp b/Core/src/Geometry/PortalError.cpp new file mode 100644 index 00000000000..bbedd3010e5 --- /dev/null +++ b/Core/src/Geometry/PortalError.cpp @@ -0,0 +1,39 @@ +// 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/. + +#include "Acts/Surfaces/SurfaceError.hpp" + +#include + +namespace { + +class SurfaceErrorCategory : public std::error_category { + public: + // Return a short descriptive name for the category. + const char* name() const noexcept final { return "SurfaceError"; } + + // Return what each enum means in text. + std::string message(int c) const final { + using Acts::SurfaceError; + + switch (static_cast(c)) { + case SurfaceError::GlobalPositionNotOnSurface: + return "Global to local transformation failed: position not on " + "surface."; + default: + return "unknown"; + } + } +}; + +} // namespace + +std::error_code Acts::make_error_code(Acts::SurfaceError e) { + static SurfaceErrorCategory c; + return {static_cast(e), c}; +} diff --git a/Core/src/Geometry/PortalLinkBase.cpp b/Core/src/Geometry/PortalLinkBase.cpp new file mode 100644 index 00000000000..9a6be75794c --- /dev/null +++ b/Core/src/Geometry/PortalLinkBase.cpp @@ -0,0 +1,173 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/. + +#include "Acts/Geometry/PortalLinkBase.hpp" + +#include "Acts/Geometry/CompositePortalLink.hpp" +#include "Acts/Geometry/GridPortalLink.hpp" +#include "Acts/Geometry/TrivialPortalLink.hpp" +#include "Acts/Surfaces/CylinderSurface.hpp" +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" +#include "Acts/Surfaces/RegularSurface.hpp" +#include "Acts/Utilities/ThrowAssert.hpp" + +namespace Acts { + +void PortalLinkBase::checkMergePreconditions(const PortalLinkBase& a, + const PortalLinkBase& b, + BinningValue direction) { + const auto& surfaceA = a.surface(); + const auto& surfaceB = b.surface(); + + throw_assert(&surfaceA != &surfaceB, + "Cannot merge portals to the same surface"); + + throw_assert(surfaceA.type() == surfaceB.type(), + "Cannot merge portals of different surface types"); + + throw_assert(surfaceA.bounds().type() == surfaceB.bounds().type(), + "Cannot merge portals of different surface bounds"); + + if (const auto* cylA = dynamic_cast(&surfaceA); + cylA != nullptr) { + const auto* cylB = dynamic_cast(&surfaceB); + throw_assert(cylB != nullptr, + "Cannot merge CylinderSurface with " + "non-CylinderSurface"); + throw_assert( + direction == BinningValue::binZ || direction == BinningValue::binRPhi, + "Invalid binning direction: " + binningValueName(direction)); + } else if (const auto* discA = dynamic_cast(&surfaceA); + discA != nullptr) { + const auto* discB = dynamic_cast(&surfaceB); + throw_assert(discB != nullptr, + "Cannot merge DiscSurface with non-DiscSurface"); + throw_assert( + direction == BinningValue::binR || direction == BinningValue::binPhi, + "Invalid binning direction: " + binningValueName(direction)); + + throw_assert(dynamic_cast(&discA->bounds()) && + dynamic_cast(&discB->bounds()), + "DiscSurface bounds must be RadialBounds"); + + } else { + throw std::logic_error{"Surface type is not supported"}; + } +} + +std::unique_ptr PortalLinkBase::merge( + std::unique_ptr a, std::unique_ptr b, + BinningValue direction, const Logger& logger) { + ACTS_DEBUG("Merging two arbitrary portals"); + + ACTS_VERBOSE(" - a: " << *a); + ACTS_VERBOSE(" - b: " << *b); + + checkMergePreconditions(*a, *b, direction); + + // Three options: + // 1. Grid + // 2. Trivial + // 3. Composite + + // Grid Grid + // Grid Trivial + // Grid Composite + // Trivial Grid + // Trivial Trivial + // Trivial Composite + // Composite Grid + // Composite Trivial + // Composite Composite + + auto gridMerge = + [&](const GridPortalLink& aGrid, + const GridPortalLink& bGrid) -> std::unique_ptr { + assert(a != nullptr); + assert(b != nullptr); + auto merged = GridPortalLink::merge(aGrid, bGrid, direction, logger); + if (merged == nullptr) { + return std::make_unique(std::move(a), std::move(b), + direction); + } + return merged; + }; + + if (const auto* aGrid = dynamic_cast(a.get()); + aGrid != nullptr) { + if (const auto* bGrid = dynamic_cast(b.get()); + bGrid != nullptr) { + ACTS_VERBOSE("Merging two grid portals"); + return gridMerge(*aGrid, *bGrid); + + } else if (const auto* bTrivial = + dynamic_cast(b.get()); + bTrivial != nullptr) { + ACTS_VERBOSE("Merging a grid portal with a trivial portal"); + return gridMerge(*aGrid, *bTrivial->makeGrid(direction)); + + } else if (dynamic_cast(b.get()) != nullptr) { + ACTS_WARNING("Merging a grid portal with a composite portal"); + return std::make_unique(std::move(a), std::move(b), + direction); + + } else { + throw std::logic_error{"Portal link type is not supported"}; + } + + } else if (const auto* aTrivial = + dynamic_cast(a.get()); + aTrivial != nullptr) { + if (const auto* bGrid = dynamic_cast(b.get()); + bGrid) { + ACTS_VERBOSE("Merging a trivial portal with a grid portal"); + return gridMerge(*aTrivial->makeGrid(direction), *bGrid); + + } else if (const auto* bTrivial = + dynamic_cast(b.get()); + bTrivial != nullptr) { + ACTS_VERBOSE("Merging two trivial portals"); + return gridMerge(*aTrivial->makeGrid(direction), + *bTrivial->makeGrid(direction)); + + } else if (dynamic_cast(b.get()) != nullptr) { + ACTS_WARNING("Merging a trivial portal with a composite portal"); + return std::make_unique(std::move(a), std::move(b), + direction); + + } else { + throw std::logic_error{"Portal link type is not supported"}; + } + + } else if (dynamic_cast(a.get()) != nullptr) { + if (dynamic_cast(b.get()) != nullptr) { + ACTS_WARNING("Merging a composite portal with a grid portal"); + return std::make_unique(std::move(a), std::move(b), + direction); + + } else if (dynamic_cast(b.get()) != nullptr) { + ACTS_WARNING("Merging a composite portal with a trivial portal"); + return std::make_unique(std::move(a), std::move(b), + direction); + + } else if (dynamic_cast(b.get()) != nullptr) { + ACTS_WARNING("Merging two composite portals"); + return std::make_unique(std::move(a), std::move(b), + direction); + + } else { + throw std::logic_error{"Portal link type is not supported"}; + } + + } else { + throw std::logic_error{"Portal link type is not supported"}; + } +} + +} // namespace Acts diff --git a/Core/src/Geometry/TrivialPortalLink.cpp b/Core/src/Geometry/TrivialPortalLink.cpp new file mode 100644 index 00000000000..6c699b9a320 --- /dev/null +++ b/Core/src/Geometry/TrivialPortalLink.cpp @@ -0,0 +1,45 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/. + +#include "Acts/Geometry/TrivialPortalLink.hpp" + +#include "Acts/Geometry/GridPortalLink.hpp" +#include "Acts/Geometry/TrackingVolume.hpp" + +#include + +namespace Acts { + +std::unique_ptr TrivialPortalLink::makeGrid( + BinningValue direction) const { + return GridPortalLink::make(m_surface, *m_volume, direction); +} + +Result TrivialPortalLink::resolveVolume( + const GeometryContext& /*gctx*/, const Vector2& /*position*/, + double /*tolerance*/) const { + return m_volume; +} + +Result TrivialPortalLink::resolveVolume( + const GeometryContext& gctx, const Vector3& position, + double tolerance) const { + static_cast(gctx); + static_cast(position); + static_cast(tolerance); + assert(m_surface->isOnSurface(gctx, position, BoundaryTolerance::None(), + tolerance) && + "Trivial portal lookup point should be on surface"); + return m_volume; +} + +void TrivialPortalLink::toStream(std::ostream& os) const { + os << "TrivialPortalLink"; +} + +} // namespace Acts diff --git a/Core/src/Surfaces/SurfaceError.cpp b/Core/src/Surfaces/SurfaceError.cpp index bbedd3010e5..a9b7ad7e4ba 100644 --- a/Core/src/Surfaces/SurfaceError.cpp +++ b/Core/src/Surfaces/SurfaceError.cpp @@ -1,30 +1,29 @@ // This file is part of the Acts project. // -// Copyright (C) 2020 CERN for the benefit of the Acts project +// Copyright (C) 2024 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/. -#include "Acts/Surfaces/SurfaceError.hpp" +#include "Acts/Geometry/PortalError.hpp" #include namespace { -class SurfaceErrorCategory : public std::error_category { +class PortalErrorCategory : public std::error_category { public: // Return a short descriptive name for the category. - const char* name() const noexcept final { return "SurfaceError"; } + const char* name() const noexcept final { return "PortalError"; } // Return what each enum means in text. std::string message(int c) const final { - using Acts::SurfaceError; + using Acts::PortalError; - switch (static_cast(c)) { - case SurfaceError::GlobalPositionNotOnSurface: - return "Global to local transformation failed: position not on " - "surface."; + switch (static_cast(c)) { + case PortalError::PositionNotOnAnyChildPortalLink: + return "Position not on any of the composite child portal links"; default: return "unknown"; } @@ -33,7 +32,7 @@ class SurfaceErrorCategory : public std::error_category { } // namespace -std::error_code Acts::make_error_code(Acts::SurfaceError e) { - static SurfaceErrorCategory c; +std::error_code Acts::make_error_code(Acts::PortalError e) { + static PortalErrorCategory c; return {static_cast(e), c}; } diff --git a/Tests/UnitTests/Core/Geometry/CMakeLists.txt b/Tests/UnitTests/Core/Geometry/CMakeLists.txt index 00ee68e3c80..3b058a882b4 100644 --- a/Tests/UnitTests/Core/Geometry/CMakeLists.txt +++ b/Tests/UnitTests/Core/Geometry/CMakeLists.txt @@ -32,3 +32,4 @@ add_unittest(TrapezoidVolumeBounds TrapezoidVolumeBoundsTests.cpp) add_unittest(VolumeBounds VolumeBoundsTests.cpp) add_unittest(Volume VolumeTests.cpp) add_unittest(CylinderVolumeStack CylinderVolumeStackTests.cpp) +add_unittest(PortalLink PortalLinkTests.cpp) diff --git a/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp b/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp new file mode 100644 index 00000000000..727e01c9f69 --- /dev/null +++ b/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp @@ -0,0 +1,2322 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 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/. + +#include +#include +#include +#include + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Units.hpp" +#include "Acts/Geometry/CompositePortalLink.hpp" +#include "Acts/Geometry/CylinderVolumeBounds.hpp" +#include "Acts/Geometry/GridPortalLink.hpp" +#include "Acts/Geometry/TrackingVolume.hpp" +#include "Acts/Geometry/TrivialPortalLink.hpp" +#include "Acts/Surfaces/CylinderSurface.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" +#include "Acts/Surfaces/SurfaceMergingException.hpp" +#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp" +#include "Acts/Utilities/AxisFwd.hpp" +#include "Acts/Utilities/ThrowAssert.hpp" + +#include +#include +#include +#include + +using namespace Acts::UnitLiterals; + +namespace Acts::Test { + +auto logger = Acts::getDefaultLogger("UnitTests", Acts::Logging::VERBOSE); + +struct Fixture { + Logging::Level m_level; + Fixture() { + m_level = Acts::Logging::getFailureThreshold(); + Acts::Logging::setFailureThreshold(Acts::Logging::FATAL); + } + + ~Fixture() { Acts::Logging::setFailureThreshold(m_level); } +}; + +std::shared_ptr makeDummyVolume() { + return std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); +} + +GeometryContext gctx; + +template +std::unique_ptr copy(const std::unique_ptr& p) { + return std::make_unique(*p); +} + +template +void visitBins(const link_t& link, + const std::function& func) { + auto& grid = link.grid(); + auto loc = grid.numLocalBins(); + if constexpr (std::decay_t::DIM == 1) { + for (std::size_t i = 1; i <= loc[0]; i++) { + func(grid.atLocalBins({i})); + } + } else { + for (std::size_t i = 1; i <= loc[0]; i++) { + for (std::size_t j = 1; j <= loc[1]; j++) { + func(grid.atLocalBins({i, j})); + } + } + } +} + +BOOST_FIXTURE_TEST_SUITE(Geometry, Fixture) + +BOOST_AUTO_TEST_SUITE(GridConstruction) + +BOOST_AUTO_TEST_CASE(Cylinder) { + BOOST_TEST_CONTEXT("1D") { + auto cyl = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm); + + // Volume for bin testing + auto vol = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + // Incompatible binning + BOOST_CHECK_THROW( + GridPortalLink::make(cyl, BinningValue::binZ, Axis{AxisBound, 0, 5, 5}), + std::invalid_argument); + + auto grid1dCyl = GridPortalLink::make(cyl, BinningValue::binZ, + Axis{AxisBound, -100_mm, 100_mm, 10}); + BOOST_REQUIRE(grid1dCyl); + grid1dCyl->setVolume(vol.get()); + + // Throws because non-closed axis + BOOST_CHECK_THROW(GridPortalLink::make(cyl, BinningValue::binRPhi, + Axis{AxisBound, -180_degree * 30_mm, + 180_degree * 30_mm, 10}), + std::invalid_argument); + + auto grid1dCylRPhi = GridPortalLink::make( + cyl, BinningValue::binRPhi, + Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 10}); + BOOST_REQUIRE_NE(grid1dCylRPhi, nullptr); + grid1dCylRPhi->setVolume(vol.get()); + + Axis axisExpected{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 10}; + BOOST_CHECK_EQUAL(grid1dCylRPhi->grid().axes().size(), 1); + const auto& axis = *grid1dCylRPhi->grid().axes().front(); + BOOST_CHECK_EQUAL(axis, axisExpected); + + // Another cylinder, shifted in z + auto cyl2 = Surface::makeShared( + Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm); + + auto grid1dCyl2 = GridPortalLink::make(cyl2, BinningValue::binZ, + Axis{AxisBound, -50_mm, 50_mm, 5}); + + // Test exception on cylinder with non-zero average phi + auto cylNonZeroAverage = Surface::makeShared( + Transform3::Identity(), 30_mm, 100_mm, 20_degree, 45_degree); + BOOST_CHECK_THROW( + GridPortalLink::make(cylNonZeroAverage, BinningValue::binZ, + Axis{AxisBound, -100_mm, 100_mm, 10}), + std::invalid_argument); + + auto checkAllBins = [&](const auto& link) { + visitBins(link, [&](const TrackingVolume* content) { + BOOST_CHECK_EQUAL(content, vol.get()); + }); + }; + + checkAllBins(*grid1dCyl); + checkAllBins(*grid1dCylRPhi); + + // Extend to a 2D grid with auto phi binning + + auto grid2dCyl1 = grid1dCyl->extendTo2d(nullptr); + BOOST_REQUIRE(grid2dCyl1); + BOOST_CHECK_EQUAL(grid2dCyl1->grid().axes().size(), 2); + BOOST_CHECK_EQUAL(grid2dCyl1->surface().bounds(), cyl->bounds()); + const auto* axis1 = grid2dCyl1->grid().axes().front(); + const auto* axis2 = grid2dCyl1->grid().axes().back(); + + Axis axis1Expected{AxisClosed, -M_PI * 30_mm, M_PI * 30_mm, 1}; + BOOST_CHECK_EQUAL(*axis1, axis1Expected); + Axis axis2Expected{AxisBound, -100_mm, 100_mm, 10}; + BOOST_CHECK_EQUAL(*axis2, axis2Expected); + + auto& concrete = dynamic_cast< + GridPortalLinkT&>( + *grid2dCyl1); + + checkAllBins(concrete); + + Axis axis1Explicit{AxisClosed, -M_PI * 30_mm, M_PI * 30_mm, 13}; + auto grid2dCyl1Explicit = grid1dCyl->extendTo2d(&axis1Explicit); + BOOST_REQUIRE(grid2dCyl1Explicit); + BOOST_CHECK_EQUAL(grid2dCyl1Explicit->grid().axes().size(), 2); + axis1 = grid2dCyl1Explicit->grid().axes().front(); + axis2 = grid2dCyl1Explicit->grid().axes().back(); + + BOOST_CHECK_EQUAL(*axis1, axis1Explicit); + BOOST_CHECK_EQUAL(*axis2, axis2Expected); + + auto& concrete2 = dynamic_cast< + GridPortalLinkT&>( + *grid2dCyl1Explicit); + + checkAllBins(concrete2); + + auto cylPhi = Surface::makeShared( + Transform3::Identity(), 30_mm, 100_mm, 45_degree); + std::unique_ptr grid1dCylPhi = GridPortalLink::make( + cylPhi, BinningValue::binZ, Axis{AxisBound, -100_mm, 100_mm, 10}); + + grid1dCylPhi->setVolume(vol.get()); + + // Check that phi sector portal does not accept closed axis + BOOST_CHECK_THROW(GridPortalLink::make(cylPhi, BinningValue::binRPhi, + Axis{AxisClosed, -45_degree * 30_mm, + 45_degree * 30_mm, 10}), + std::invalid_argument); + + auto grid2dCylPhi = grid1dCylPhi->extendTo2d(nullptr); + BOOST_CHECK_EQUAL(grid2dCylPhi->grid().axes().size(), 2); + BOOST_CHECK_EQUAL(grid2dCylPhi->surface().bounds(), cylPhi->bounds()); + const auto* axis1Phi = grid2dCylPhi->grid().axes().front(); + const auto* axis2Phi = grid2dCylPhi->grid().axes().back(); + + Axis axis1PhiExpected{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 1}; + BOOST_CHECK_EQUAL(*axis1Phi, axis1PhiExpected); + Axis axis2PhiExpected{AxisBound, -100_mm, 100_mm, 10}; + BOOST_CHECK_EQUAL(*axis2Phi, axis2PhiExpected); + + auto& concrete3 = + dynamic_cast&>( + *grid2dCylPhi); + + checkAllBins(concrete3); + + Axis axis1PhiExplicit{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 13}; + auto grid2dCylPhiExplicit = grid1dCylPhi->extendTo2d(&axis1PhiExplicit); + BOOST_REQUIRE(grid2dCylPhiExplicit); + BOOST_CHECK_EQUAL(grid2dCylPhiExplicit->grid().axes().size(), 2); + axis1Phi = grid2dCylPhiExplicit->grid().axes().front(); + axis2Phi = grid2dCylPhiExplicit->grid().axes().back(); + BOOST_CHECK_EQUAL(*axis1Phi, axis1PhiExplicit); + BOOST_CHECK_EQUAL(*axis2Phi, axis2PhiExpected); + + auto& concrete4 = + dynamic_cast&>( + *grid2dCylPhiExplicit); + + checkAllBins(concrete4); + } + + BOOST_TEST_CONTEXT("2D") { + auto cyl = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 45_degree); + + // z bad, rphi bad + BOOST_CHECK_THROW(GridPortalLink::make(cyl, Axis{AxisBound, 1, 2, 5}, + Axis{AxisBound, 3_mm, 4_mm, 5}), + std::invalid_argument); + + // z good, rphi bad + BOOST_CHECK_THROW(GridPortalLink::make(cyl, Axis{AxisBound, 3_mm, 4_mm, 5}, + Axis{AxisBound, -100_mm, 100_m, 5}), + std::invalid_argument); + + // z bad, rphi good + BOOST_CHECK_THROW( + GridPortalLink::make( + cyl, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5}, + Axis{AxisBound, -80_mm, 100_mm, 5}), + std::invalid_argument); + + auto grid1 = GridPortalLink::make( + cyl, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5}, + Axis{AxisBound, -100_mm, 100_mm, 5}); + + auto cylFull = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm); + + // Throws because non-closed axis + BOOST_CHECK_THROW(GridPortalLink::make(cylFull, + Axis{AxisBound, -180_degree * 30_mm, + 180_degree * 30_mm, 5}, + Axis{AxisBound, -100_mm, 100_mm, 5}), + std::invalid_argument); + + auto gridFull = GridPortalLink::make( + cylFull, Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5}, + Axis{AxisBound, -100_mm, 100_mm, 5}); + + BOOST_CHECK_EQUAL(gridFull->grid().axes().size(), 2); + BOOST_CHECK_EQUAL(gridFull->grid().axes().size(), 2); + Axis axisFullExpected{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, + 5}; + const auto& axisFull = *gridFull->grid().axes().front(); + BOOST_CHECK_EQUAL(axisFull, axisFullExpected); + } +} + +BOOST_AUTO_TEST_CASE(Disc) { + using enum AxisType; + using enum AxisBoundaryType; + BOOST_TEST_CONTEXT("1D") { + auto disc1 = + Surface::makeShared(Transform3::Identity(), 30_mm, 100_mm); + + // Volume for bin testing + auto vol = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + BOOST_CHECK_THROW(GridPortalLink::make(disc1, BinningValue::binZ, + Axis{AxisBound, 30_mm, 100_mm, 3}), + std::invalid_argument); + + // Check exception for full disc and non-closed phi axis + BOOST_CHECK_THROW( + GridPortalLink::make(disc1, BinningValue::binPhi, + Axis{AxisBound, -180_degree, 180_degree, 3}), + std::invalid_argument); + + auto grid1 = GridPortalLink::make(disc1, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 3}); + BOOST_REQUIRE_NE(grid1, nullptr); + BOOST_CHECK_EQUAL(grid1->grid().axes().size(), 1); + const auto& axis = *grid1->grid().axes().front(); + Axis axis1Expected{AxisBound, 30_mm, 100_mm, 3}; + BOOST_CHECK_EQUAL(axis, axis1Expected); + + grid1->setVolume(vol.get()); + + auto discPhi = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 45_degree); + + // Check thet disc with phi sector does not accept closed axis + BOOST_CHECK_THROW( + GridPortalLink::make(discPhi, BinningValue::binPhi, + Axis{AxisClosed, -45_degree, 45_degree, 3}), + std::invalid_argument); + + auto gridPhi = + GridPortalLink::make(discPhi, BinningValue::binPhi, + Axis{AxisBound, -45_degree, 45_degree, 3}); + BOOST_REQUIRE_NE(gridPhi, nullptr); + gridPhi->setVolume(vol.get()); + + // Test exception on disc with non-zero average phi + auto discNonZeroAverage = Surface::makeShared( + Transform3::Identity(), + std::make_shared(30_mm, 100_mm, 45_degree, 75_degree)); + BOOST_CHECK_THROW( + GridPortalLink::make(discNonZeroAverage, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 3}), + std::invalid_argument); + + BOOST_CHECK_EQUAL(gridPhi->grid().axes().size(), 1); + const auto& axisPhi = *gridPhi->grid().axes().front(); + Axis axisPhi1Expected{AxisBound, -45_degree, 45_degree, 3}; + BOOST_CHECK_EQUAL(axisPhi, axisPhi1Expected); + + auto checkAllBins = [&](const auto& grid) { + visitBins(grid, [&](const TrackingVolume* content) { + BOOST_CHECK_EQUAL(content, vol.get()); + }); + }; + + checkAllBins(*grid1); + checkAllBins(*gridPhi); + + // Test making 2D grids from the 1D ones + auto grid2d = grid1->extendTo2d(nullptr); + BOOST_REQUIRE(grid2d); + BOOST_CHECK_EQUAL(grid2d->grid().axes().size(), 2); + const auto* axis1 = grid2d->grid().axes().front(); + BOOST_CHECK_EQUAL(*axis1, axis1Expected); + const auto* axis2 = grid2d->grid().axes().back(); + BOOST_CHECK_CLOSE(axis2->getMin(), -180_degree, 1e-6); + BOOST_CHECK_CLOSE(axis2->getMax(), 180_degree, 1e-6); + BOOST_CHECK_EQUAL(axis2->getNBins(), 1); + BOOST_CHECK_EQUAL(axis2->getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis2->getBoundaryType(), AxisBoundaryType::Closed); + + checkAllBins( + dynamic_cast>&>(*grid2d)); + + Axis axis2Explicit{AxisClosed, -180_degree, 180_degree, 3}; + auto grid2dExplicit = grid1->extendTo2d(&axis2Explicit); + BOOST_REQUIRE(grid2dExplicit); + BOOST_CHECK_EQUAL(grid2dExplicit->grid().axes().size(), 2); + axis1 = grid2dExplicit->grid().axes().front(); + axis2 = grid2dExplicit->grid().axes().back(); + BOOST_CHECK_EQUAL(*axis1, axis1Expected); + BOOST_CHECK_EQUAL(*axis2, axis2Explicit); + + checkAllBins(dynamic_cast&>( + *grid2dExplicit)); + + auto gridPhiBinnedInR = GridPortalLink::make( + discPhi, BinningValue::binR, Axis{AxisBound, 30_mm, 100_mm, 3}); + gridPhiBinnedInR->setVolume(vol.get()); + auto grid2dPhiNonClosed = gridPhiBinnedInR->extendTo2d(nullptr); + BOOST_REQUIRE(grid2dPhiNonClosed); + BOOST_CHECK_EQUAL(grid2dPhiNonClosed->grid().axes().size(), 2); + Axis gridPhiBinnedInRExpected{AxisBound, 30_mm, 100_mm, 3}; + BOOST_CHECK_EQUAL(*grid2dPhiNonClosed->grid().axes().front(), + gridPhiBinnedInRExpected); + const auto* axisPhiNonClosed = grid2dPhiNonClosed->grid().axes().back(); + BOOST_CHECK_CLOSE(axisPhiNonClosed->getMin(), -45_degree, 1e-6); + BOOST_CHECK_CLOSE(axisPhiNonClosed->getMax(), 45_degree, 1e-6); + BOOST_CHECK_EQUAL(axisPhiNonClosed->getNBins(), 1); + BOOST_CHECK_EQUAL(axisPhiNonClosed->getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axisPhiNonClosed->getBoundaryType(), + AxisBoundaryType::Bound); + + checkAllBins( + dynamic_cast>&>( + *grid2dPhiNonClosed)); + + Axis axisPhiNonClosedExplicit{AxisBound, -45_degree, 45_degree, 3}; + auto grid2dPhiNonClosedExplicit = + gridPhiBinnedInR->extendTo2d(&axisPhiNonClosedExplicit); + BOOST_REQUIRE(grid2dPhiNonClosedExplicit); + BOOST_CHECK_EQUAL(grid2dPhiNonClosedExplicit->grid().axes().size(), 2); + axisPhiNonClosed = grid2dPhiNonClosedExplicit->grid().axes().back(); + BOOST_CHECK_EQUAL(*axisPhiNonClosed, axisPhiNonClosedExplicit); + BOOST_CHECK_EQUAL(*grid2dPhiNonClosedExplicit->grid().axes().front(), + gridPhiBinnedInRExpected); + + checkAllBins( + dynamic_cast&>( + *grid2dPhiNonClosedExplicit)); + + auto grid2dPhi = gridPhi->extendTo2d(nullptr); + BOOST_REQUIRE(grid2dPhi); + BOOST_CHECK_EQUAL(grid2dPhi->grid().axes().size(), 2); + Axis axis2dPhiExpected{AxisBound, 30_mm, 100_mm, 1}; + BOOST_CHECK_EQUAL(*grid2dPhi->grid().axes().front(), axis2dPhiExpected); + BOOST_CHECK_EQUAL(*grid2dPhi->grid().axes().back(), axisPhi1Expected); + + checkAllBins( + dynamic_cast&>(*grid2dPhi)); + + Axis axis2dPhiExplicit{AxisBound, 30_mm, 100_mm, 3}; + auto grid2dPhiExplicit = gridPhi->extendTo2d(&axis2dPhiExplicit); + BOOST_REQUIRE(grid2dPhiExplicit); + BOOST_CHECK_EQUAL(grid2dPhiExplicit->grid().axes().size(), 2); + BOOST_CHECK_EQUAL(*grid2dPhiExplicit->grid().axes().front(), + axis2dPhiExplicit); + BOOST_CHECK_EQUAL(*grid2dPhiExplicit->grid().axes().back(), + axisPhi1Expected); + + checkAllBins(dynamic_cast&>( + *grid2dPhiExplicit)); + } + + BOOST_TEST_CONTEXT("2D") { + auto discPhi = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 45_degree); + + Axis rBad{AxisBound, 1, 2, 5}; + Axis rGood{AxisBound, 30_mm, 100_mm, 5}; + Axis phiBad{AxisBound, 1, 2, 5}; + Axis phiGood{AxisBound, -45_degree, 45_degree, 5}; + + // r bad, phi bad + BOOST_CHECK_THROW(GridPortalLink::make(discPhi, rBad, phiBad), + std::invalid_argument); + // r bad, phi good + BOOST_CHECK_THROW(GridPortalLink::make(discPhi, rBad, phiGood), + std::invalid_argument); + // r good, phi bad + BOOST_CHECK_THROW(GridPortalLink::make(discPhi, rGood, phiBad), + std::invalid_argument); + // r good phi good + auto grid = GridPortalLink::make(discPhi, rGood, phiGood); + BOOST_REQUIRE_NE(grid, nullptr); + } +} + +BOOST_AUTO_TEST_CASE(Plane) { + // @TODO: Add plane tests +} + +BOOST_AUTO_TEST_CASE(FromTrivial) { + BOOST_TEST_CONTEXT("Cylinder") { + auto cyl = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm); + + auto vol = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto trivial = std::make_unique(cyl, *vol); + BOOST_REQUIRE(trivial); + + BOOST_CHECK_EQUAL(trivial->resolveVolume(gctx, Vector2{1, 2}).value(), + vol.get()); + + auto gridZ = trivial->makeGrid(BinningValue::binZ); + BOOST_REQUIRE(gridZ); + + BOOST_CHECK_EQUAL(gridZ->grid().axes().size(), 1); + BOOST_CHECK_EQUAL(gridZ->surface().bounds(), cyl->bounds()); + Axis axisZExpected{AxisBound, -100_mm, 100_mm, 1}; + BOOST_CHECK_EQUAL(*gridZ->grid().axes().front(), axisZExpected); + + BOOST_CHECK_EQUAL( + gridZ->resolveVolume(gctx, Vector2{20_degree * 30_mm, 90_mm}).value(), + vol.get()); + + auto gridRPhi = trivial->makeGrid(BinningValue::binRPhi); + BOOST_REQUIRE(gridRPhi); + + BOOST_CHECK_EQUAL(gridRPhi->grid().axes().size(), 1); + BOOST_CHECK_EQUAL(gridRPhi->surface().bounds(), cyl->bounds()); + Axis axisRPhiExpected{AxisClosed, -M_PI * 30_mm, M_PI * 30_mm, 1}; + BOOST_CHECK_EQUAL(*gridRPhi->grid().axes().front(), axisRPhiExpected); + + auto cylPhi = Surface::makeShared( + Transform3::Identity(), 30_mm, 100_mm, 30_degree); + + auto trivialPhi = std::make_unique(cylPhi, *vol); + BOOST_REQUIRE(trivialPhi); + + BOOST_CHECK_EQUAL(trivialPhi->resolveVolume(gctx, Vector2{1, 2}).value(), + vol.get()); + + auto gridRPhiSector = trivialPhi->makeGrid(BinningValue::binRPhi); + BOOST_REQUIRE(gridRPhiSector); + + BOOST_CHECK_EQUAL( + gridRPhiSector->resolveVolume(gctx, Vector2{20_degree * 30_mm, 90_mm}) + .value(), + vol.get()); + + BOOST_CHECK_EQUAL(gridRPhiSector->grid().axes().size(), 1); + Axis axisRPhiSectorExpected{AxisBound, -30_degree * 30_mm, + 30_degree * 30_mm, 1}; + BOOST_CHECK_EQUAL(*gridRPhiSector->grid().axes().front(), + axisRPhiSectorExpected); + } + + BOOST_TEST_CONTEXT("Disc") { + auto disc = + Surface::makeShared(Transform3::Identity(), 30_mm, 100_mm); + + auto vol = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto trivial = std::make_unique(disc, *vol); + BOOST_REQUIRE(trivial); + + // Doesn't matter which position + BOOST_CHECK_EQUAL(trivial->resolveVolume(gctx, Vector2{1, 2}).value(), + vol.get()); + + auto gridR = trivial->makeGrid(BinningValue::binR); + BOOST_REQUIRE(gridR); + + BOOST_CHECK_EQUAL(gridR->grid().axes().size(), 1); + BOOST_CHECK_EQUAL(gridR->surface().bounds(), disc->bounds()); + Axis axisRExpected{AxisBound, 30_mm, 100_mm, 1}; + BOOST_CHECK_EQUAL(*gridR->grid().axes().front(), axisRExpected); + + BOOST_CHECK_EQUAL( + gridR->resolveVolume(gctx, Vector2{90_mm, 10_degree}).value(), + vol.get()); + + auto gridPhi = trivial->makeGrid(BinningValue::binPhi); + BOOST_REQUIRE(gridPhi); + + BOOST_CHECK_EQUAL(gridPhi->grid().axes().size(), 1); + BOOST_CHECK_EQUAL(gridPhi->surface().bounds(), disc->bounds()); + Axis axisPhiExpected{AxisClosed, -M_PI, M_PI, 1}; + BOOST_CHECK_EQUAL(*gridPhi->grid().axes().front(), axisPhiExpected); + + BOOST_CHECK_EQUAL( + gridPhi->resolveVolume(gctx, Vector2{90_mm, 10_degree}).value(), + vol.get()); + } +} + +BOOST_AUTO_TEST_SUITE_END() // GridConstruction + +BOOST_AUTO_TEST_SUITE(GridMerging) + +BOOST_AUTO_TEST_SUITE(Merging1dCylinder) + +BOOST_AUTO_TEST_SUITE(ZDirection) + +BOOST_AUTO_TEST_CASE(ColinearMerge) { + auto vol1 = makeDummyVolume(); + auto vol2 = makeDummyVolume(); + + auto cyl = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm); + + auto grid1dCyl = GridPortalLink::make(cyl, BinningValue::binZ, + Axis{AxisBound, -100_mm, 100_mm, 10}); + grid1dCyl->setVolume(vol1.get()); + + // Another cylinder, shifted in z + auto cyl2 = Surface::makeShared( + Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm); + + auto grid1dCyl2 = GridPortalLink::make(cyl2, BinningValue::binZ, + Axis{AxisBound, -50_mm, 50_mm, 5}); + + grid1dCyl2->setVolume(vol2.get()); + + // Completely invalid + BOOST_CHECK_THROW(GridPortalLink::merge(*grid1dCyl, *grid1dCyl2, + BinningValue::binPhi, *logger), + AssertionFailureException); + // Invalid direction, as the cylinders are shifted in z, and can't be merged + // in r x phi + BOOST_CHECK_THROW(GridPortalLink::merge(*grid1dCyl, *grid1dCyl2, + BinningValue::binRPhi, *logger), + SurfaceMergingException); + + BOOST_TEST_CONTEXT("Consistent equidistant") { + auto mergedPtr = GridPortalLink::merge(*grid1dCyl, *grid1dCyl2, + BinningValue::binZ, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_EQUAL(axis.getMin(), -150_mm); + BOOST_CHECK_EQUAL(axis.getMax(), 150_mm); + BOOST_CHECK_EQUAL(axis.getNBins(), 15); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + } + + BOOST_TEST_CONTEXT("Inconsistent equidistant") { + auto grid1dCyl2BinWidthChanged = GridPortalLink::make( + cyl2, BinningValue::binZ, Axis{AxisBound, -50_mm, 50_mm, 6}); + + auto mergedPtr = GridPortalLink::merge( + *grid1dCyl, *grid1dCyl2BinWidthChanged, BinningValue::binZ, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_EQUAL(axis.getMin(), -150_mm); + BOOST_CHECK_EQUAL(axis.getMax(), 150_mm); + BOOST_CHECK_EQUAL(axis.getNBins(), 16); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + } + + BOOST_TEST_CONTEXT("Right Variable") { + auto gridLeft = GridPortalLink::make(cyl, BinningValue::binZ, + Axis{AxisBound, -100_mm, 100_mm, 10}); + + auto gridRight = + GridPortalLink::make(cyl2, BinningValue::binZ, + Axis{AxisBound, {-50_mm, -10_mm, 10_mm, 50_mm}}); + + auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight, + BinningValue::binZ, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_EQUAL(axis.getMin(), -150_mm); + BOOST_CHECK_EQUAL(axis.getMax(), 150_mm); + BOOST_CHECK_EQUAL(axis.getNBins(), 13); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + } + + BOOST_TEST_CONTEXT("Left Variable") { + auto gridLeft = + GridPortalLink::make(cyl, BinningValue::binZ, + Axis{AxisBound, {-100_mm, -80_mm, 10_mm, 100_mm}}); + + auto gridRight = GridPortalLink::make(cyl2, BinningValue::binZ, + Axis{AxisBound, -50_mm, 50_mm, 8}); + + auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight, + BinningValue::binZ, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_EQUAL(axis.getMin(), -150_mm); + BOOST_CHECK_EQUAL(axis.getMax(), 150_mm); + BOOST_CHECK_EQUAL(axis.getNBins(), 11); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + } + + BOOST_TEST_CONTEXT("Both Variable") { + auto gridLeft = + GridPortalLink::make(cyl, BinningValue::binZ, + Axis{AxisBound, {-100_mm, -80_mm, 10_mm, 100_mm}}); + + auto gridRight = + GridPortalLink::make(cyl2, BinningValue::binZ, + Axis{AxisBound, {-50_mm, -10_mm, 10_mm, 50_mm}}); + + auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight, + BinningValue::binZ, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_EQUAL(axis.getMin(), -150_mm); + BOOST_CHECK_EQUAL(axis.getMax(), 150_mm); + BOOST_CHECK_EQUAL(axis.getNBins(), 6); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + } + + BOOST_TEST_CONTEXT("Non bound axis") { + std::unique_ptr gridLeft = + GridPortalLink::make(cyl, BinningValue::binZ, + Axis{AxisBound, {-100_mm, -80_mm, 10_mm, 100_mm}}); + std::unique_ptr gridRightClosed = + GridPortalLink::make(cyl2, BinningValue::binZ, + Axis{AxisClosed, {-50_mm, -10_mm, 10_mm, 50_mm}}); + std::unique_ptr gridRightOpen = + GridPortalLink::make(cyl2, BinningValue::binZ, + Axis{AxisOpen, {-50_mm, -10_mm, 10_mm, 50_mm}}); + + auto compositeLR = PortalLinkBase::merge( + copy(gridLeft), copy(gridRightClosed), BinningValue::binZ, *logger); + BOOST_CHECK_NE(dynamic_cast(compositeLR.get()), + nullptr); + auto compositeRL = PortalLinkBase::merge( + copy(gridLeft), copy(gridRightOpen), BinningValue::binZ, *logger); + BOOST_CHECK_NE(dynamic_cast(compositeRL.get()), + nullptr); + } +} + +BOOST_AUTO_TEST_CASE(ParallelMerge) { + // Merge in z direction with phi sectors + auto cyl1 = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 35_degree); + + auto grid1 = GridPortalLink::make( + cyl1, BinningValue::binRPhi, + + Axis{AxisBound, -35_degree * 30_mm, 35_degree * 30_mm, 3}); + auto cyl2 = Surface::makeShared( + Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm, + 35_degree); + + auto grid2 = GridPortalLink::make( + cyl2, BinningValue::binRPhi, + Axis{AxisBound, -35_degree * 30_mm, 35_degree * 30_mm, 3}); + + auto merged12Ptr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binZ, *logger); + BOOST_REQUIRE_NE(merged12Ptr, nullptr); + auto merged12 = dynamic_cast(merged12Ptr.get()); + BOOST_REQUIRE_NE(merged12, nullptr); + + BOOST_REQUIRE_EQUAL(merged12->grid().axes().size(), 2); + + const auto& axis1 = *merged12->grid().axes().front(); + const auto& axis2 = *merged12->grid().axes().back(); + Axis axis1Expected{AxisBound, -35_degree * 30_mm, 35_degree * 30_mm, 3}; + BOOST_CHECK_EQUAL(axis1, axis1Expected); + Axis axis2Expected{AxisBound, {-150_mm, 50_mm, 150_mm}}; + BOOST_CHECK_EQUAL(axis2, axis2Expected); +} + +BOOST_AUTO_TEST_SUITE_END() // ZDirection + +BOOST_AUTO_TEST_SUITE(RPhiDirection) + +BOOST_AUTO_TEST_CASE(ColinearMerge) { + auto cyl = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm); + BOOST_CHECK_THROW(GridPortalLink::make(cyl, BinningValue::binRPhi, + Axis{AxisBound, 0, 5, 5}), + std::invalid_argument); + + auto cylNonZeroAverage = Surface::makeShared( + Transform3::Identity(), 30_mm, 100_mm, 20_degree, 45_degree); + + BOOST_CHECK_THROW( + GridPortalLink::make( + cylNonZeroAverage, BinningValue::binRPhi, + Axis{AxisBound, -20_degree * 30_mm, 20_degree * 30_mm, 5}), + std::invalid_argument); + + BOOST_TEST_CONTEXT("Colinear merge in rPhi") { + auto cylPhi1 = Surface::makeShared( + Transform3::Identity() * AngleAxis3(45_degree, Vector3::UnitZ()), 30_mm, + 100_mm, 20_degree, 0_degree); + + auto cylPhi2 = Surface::makeShared( + Transform3::Identity() * AngleAxis3(105_degree, Vector3::UnitZ()), + 30_mm, 100_mm, 40_degree, 0_degree); + + auto portalPhi1 = GridPortalLink::make( + cylPhi1, BinningValue::binRPhi, + Axis{AxisBound, -20_degree * 30_mm, 20_degree * 30_mm, 5}); + + auto portalPhi2 = GridPortalLink::make( + cylPhi2, BinningValue::binRPhi, + Axis{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 10}); + + auto cylPhi3 = Surface::makeShared( + Transform3::Identity() * AngleAxis3(45_degree, Vector3::UnitZ()), 30_mm, + 100_mm, 90_degree, 0_degree); + + auto cylPhi4 = Surface::makeShared( + Transform3::Identity() * AngleAxis3(-135_degree, Vector3::UnitZ()), + 30_mm, 100_mm, 90_degree, 0_degree); + + auto portalPhi3 = GridPortalLink::make( + cylPhi3, BinningValue::binRPhi, + Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 2}); + + auto portalPhi4 = GridPortalLink::make( + cylPhi4, BinningValue::binRPhi, + Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 2}); + + BOOST_TEST_CONTEXT("Consistent equidistant") { + auto portalMerged = GridPortalLink::merge(*portalPhi1, *portalPhi2, + BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(portalMerged, nullptr); + + const auto* merged = + dynamic_cast(portalMerged.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis.getNBins(), 15); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + + // Test that if you merge half-circles, we get a closed axis + auto portalMerged34 = GridPortalLink::merge( + *portalPhi3, *portalPhi4, BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(portalMerged34, nullptr); + + const auto* merged34 = + dynamic_cast(portalMerged34.get()); + BOOST_REQUIRE_NE(merged34, nullptr); + BOOST_CHECK_EQUAL(merged34->grid().axes().size(), 1); + const auto& axis34 = *merged34->grid().axes().front(); + BOOST_CHECK_CLOSE(axis34.getMin(), -180_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis34.getMax(), 180_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis34.getNBins(), 4); + BOOST_CHECK_EQUAL(axis34.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis34.getBoundaryType(), AxisBoundaryType::Closed); + } + + BOOST_TEST_CONTEXT("Inconsistent equidistant") { + auto portalPhi2Mod = GridPortalLink::make( + cylPhi2, BinningValue::binRPhi, + Axis{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 3}); + + auto portalMergedMod = GridPortalLink::merge( + *portalPhi1, *portalPhi2Mod, BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(portalMergedMod, nullptr); + + const auto* merged12 = + dynamic_cast(portalMergedMod.get()); + BOOST_REQUIRE_NE(merged12, nullptr); + BOOST_CHECK_EQUAL(merged12->grid().axes().size(), 1); + const auto& axis12 = *merged12->grid().axes().front(); + BOOST_CHECK_CLOSE(axis12.getMin(), -60_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis12.getMax(), 60_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis12.getNBins(), 8); + BOOST_CHECK_EQUAL(axis12.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis12.getBoundaryType(), AxisBoundaryType::Bound); + + std::vector expected12 = {-31.4159, -17.4533, -3.49066, + 10.472, 14.6608, 18.8496, + 23.0383, 27.2271, 31.4159}; + CHECK_CLOSE_OR_SMALL(axis12.getBinEdges(), expected12, 1e-4, 10e-10); + + auto portalPhi4Mod = GridPortalLink::make( + cylPhi4, BinningValue::binRPhi, + Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 1}); + + auto portalMerged34 = GridPortalLink::merge( + *portalPhi3, *portalPhi4Mod, BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(portalMerged34, nullptr); + + const auto* merged34 = + dynamic_cast(portalMerged34.get()); + BOOST_REQUIRE_NE(merged34, nullptr); + BOOST_CHECK_EQUAL(merged34->grid().axes().size(), 1); + const auto& axis34 = *merged34->grid().axes().front(); + BOOST_CHECK_CLOSE(axis34.getMin(), -180_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis34.getMax(), 180_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis34.getNBins(), 3); + BOOST_CHECK_EQUAL(axis34.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis34.getBoundaryType(), AxisBoundaryType::Closed); + + // Caution: for full-azimuth cases, the ordering is preserved, you get + // in what you get out. -> this can flip + std::vector expected34 = {-94.2478, -47.1239, 0, 94.2478}; + CHECK_CLOSE_OR_SMALL(axis34.getBinEdges(), expected34, 1e-4, 10e-10); + } + + BOOST_TEST_CONTEXT("Left variable") { + BOOST_TEST_CONTEXT("Non-closed") { + auto gridLeft = + GridPortalLink::make(cylPhi1, BinningValue::binRPhi, + Axis{AxisBound, + {-20_degree * 30_mm, -10_degree * 30_mm, + 10_degree * 30_mm, 20_degree * 30_mm}}); + auto gridRight = GridPortalLink::make( + cylPhi2, BinningValue::binRPhi, + Axis{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 3}); + + auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight, + BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + + const auto* merged = + dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis.getNBins(), 6); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + + std::vector expected = { + -31.4159, -17.4533, -3.49066, 10.472, 15.708, 26.1799, 31.4159}; + CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10); + } + + BOOST_TEST_CONTEXT("Closed") { + auto gridLeft = GridPortalLink::make( + cylPhi4, BinningValue::binRPhi, + Axis{AxisBound, + {-90_degree * 30_mm, 25_degree * 30_mm, 90_degree * 30_mm}}); + + auto gridRight = GridPortalLink::make( + cylPhi3, BinningValue::binRPhi, + Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 3}); + + auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight, + BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + + const auto* merged = + dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_CLOSE(axis.getMin(), -180_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis.getMax(), 180_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis.getNBins(), 5); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Closed); + + // Caution: for full-azimuth cases, the ordering is preserved, you get + // in what you get out. -> this can flip + std::vector expected = {-94.2478, -34.0339, 0, + 31.4159, 62.8319, 94.2478}; + CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10); + } + } + + BOOST_TEST_CONTEXT("Right variable") { + BOOST_TEST_CONTEXT("Non-closed") { + auto gridLeft = GridPortalLink::make( + cylPhi1, BinningValue::binRPhi, + Axis{AxisBound, -20_degree * 30_mm, 20_degree * 30_mm, 3}); + auto gridRight = + GridPortalLink::make(cylPhi2, BinningValue::binRPhi, + Axis{AxisBound, + {-40_degree * 30_mm, -10_degree * 30_mm, + 10_degree * 30_mm, 40_degree * 30_mm}}); + + auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight, + BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + + const auto* merged = + dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis.getNBins(), 6); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + + std::vector expected = {-31.4159, -15.708, -5.23599, 10.472, + 17.4533, 24.4346, 31.4159}; + CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10); + } + + BOOST_TEST_CONTEXT("Closed") { + auto gridLeft = GridPortalLink::make( + cylPhi4, BinningValue::binRPhi, + Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 3}); + + auto gridRight = GridPortalLink::make( + cylPhi3, BinningValue::binRPhi, + Axis{AxisBound, + {-90_degree * 30_mm, 25_degree * 30_mm, 90_degree * 30_mm}}); + + auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight, + BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + + const auto* merged = + dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_CLOSE(axis.getMin(), -180_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis.getMax(), 180_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis.getNBins(), 5); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Closed); + + // Caution: for full-azimuth cases, the ordering is preserved, you get + // in what you get out. -> this can flip + std::vector expected = {-94.2478, -62.8319, -31.4159, + 0, 60.2139, 94.2478}; + CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10); + } + } + + BOOST_TEST_CONTEXT("Both variable") { + BOOST_TEST_CONTEXT("Non-closed") { + auto gridLeft = + GridPortalLink::make(cylPhi1, BinningValue::binRPhi, + Axis{AxisBound, + {-20_degree * 30_mm, -10_degree * 30_mm, + 10_degree * 30_mm, 20_degree * 30_mm}}); + auto gridRight = GridPortalLink::make( + cylPhi2, BinningValue::binRPhi, + Axis{AxisBound, + {-40_degree * 30_mm, -5_degree * 30_mm, 40_degree * 30_mm}}); + + auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight, + BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + + const auto* merged = + dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis.getNBins(), 5); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + + std::vector expected = {-31.4159, -13.09, 10.472, + 15.708, 26.1799, 31.4159}; + CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10); + } + + BOOST_TEST_CONTEXT("Closed") { + auto gridLeft = GridPortalLink::make( + cylPhi4, BinningValue::binRPhi, + Axis{AxisBound, + {-90_degree * 30_mm, 25_degree * 30_mm, 90_degree * 30_mm}}); + + auto gridRight = + GridPortalLink::make(cylPhi3, BinningValue::binRPhi, + Axis{AxisBound, + {-90_degree * 30_mm, -10_degree * 30_mm, + 10_degree * 30_mm, 90_degree * 30_mm}}); + + auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight, + BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + + const auto* merged = + dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_CLOSE(axis.getMin(), -180_degree * 30_mm, 1e-9); + BOOST_CHECK_CLOSE(axis.getMax(), 180_degree * 30_mm, 1e-9); + BOOST_CHECK_EQUAL(axis.getNBins(), 5); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Closed); + + // Caution: for full-azimuth cases, the ordering is preserved, you get + // in what you get out. -> this can flip + std::vector expected = {-94.2478, -34.0339, 0, + 41.8879, 52.3599, 94.2478}; + CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10); + } + } + } +} + +BOOST_AUTO_TEST_CASE(ParallelMerge) { + // Merge in phi direction with z binning + auto cylPhi1 = Surface::makeShared( + Transform3::Identity() * AngleAxis3(45_degree, Vector3::UnitZ()), 30_mm, + 100_mm, 20_degree, 0_degree); + + auto cylPhi2 = Surface::makeShared( + Transform3::Identity() * AngleAxis3(85_degree, Vector3::UnitZ()), 30_mm, + 100_mm, 20_degree, 0_degree); + + auto portalPhi1 = GridPortalLink::make(cylPhi1, BinningValue::binZ, + Axis{AxisBound, -100_mm, 100_mm, 5}); + + auto portalPhi2 = GridPortalLink::make(cylPhi2, BinningValue::binZ, + Axis{AxisBound, -100_mm, 100_mm, 5}); + + auto merged12Ptr = GridPortalLink::merge(*portalPhi1, *portalPhi2, + BinningValue::binRPhi, *logger); + BOOST_REQUIRE_NE(merged12Ptr, nullptr); + auto merged12 = dynamic_cast(merged12Ptr.get()); + BOOST_REQUIRE_NE(merged12, nullptr); + + const auto& axis1 = *merged12->grid().axes().front(); + const auto& axis2 = *merged12->grid().axes().back(); + // Phi sectors were same size, should give equidistant binning + Axis axis1Expected{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 2}; + BOOST_CHECK_EQUAL(axis1, axis1Expected); + Axis axis2Expected{AxisBound, -100_mm, 100_mm, 5}; + BOOST_CHECK_EQUAL(axis2, axis2Expected); +} + +BOOST_AUTO_TEST_SUITE_END() // RPhiDirection + +BOOST_AUTO_TEST_SUITE_END() // Merging1dCylinder + +BOOST_AUTO_TEST_SUITE(Merging2dCylinder) + +BOOST_AUTO_TEST_CASE(ZDirection) { + BOOST_TEST_CONTEXT("Phi sector") { + auto cyl1 = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 45_degree); + + // z good, rphi good + auto grid1 = GridPortalLink::make( + cyl1, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5}, + Axis{AxisBound, -100_mm, 100_mm, 5}); + + auto trf2 = Transform3{Translation3{Vector3::UnitZ() * 150_mm}}; + auto cyl2 = + Surface::makeShared(trf2, 30_mm, 50_mm, 45_degree); + + // Second grid portal with compatible phi binning + auto grid2 = GridPortalLink::make( + cyl2, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5}, + Axis{AxisBound, -50_mm, 50_mm, 5}); + + // We're merging in z direction, so the phi binnings need to be the same + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binZ, *logger); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + + const auto& axis1 = *merged->grid().axes().front(); + const auto& axis2 = *merged->grid().axes().back(); + + BOOST_CHECK_EQUAL(axis1.getMin(), -45_degree * 30_mm); + BOOST_CHECK_EQUAL(axis1.getMax(), 45_degree * 30_mm); + BOOST_CHECK_EQUAL(axis1.getNBins(), 5); + BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound); + + BOOST_CHECK_EQUAL(axis2.getMin(), -150_mm); + BOOST_CHECK_EQUAL(axis2.getMax(), 150_mm); + BOOST_CHECK_EQUAL(axis2.getNBins(), 10); + BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound); + + auto grid3 = GridPortalLink::make( + cyl2, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 3}, + Axis{AxisBound, -50_mm, 50_mm, 5}); + + auto composite = PortalLinkBase::merge(copy(grid1), copy(grid3), + BinningValue::binZ, *logger); + BOOST_CHECK_NE(dynamic_cast(composite.get()), + nullptr); + } + + BOOST_TEST_CONTEXT("Check wraparound for full circle in phi") { + auto cyl1 = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 180_degree); + + // z good, rphi good + auto grid1 = GridPortalLink::make( + cyl1, Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5}, + Axis{AxisBound, -100_mm, 100_mm, 5}); + + auto trf2 = Transform3{Translation3{Vector3::UnitZ() * 150_mm}}; + auto cyl2 = + Surface::makeShared(trf2, 30_mm, 50_mm, 180_degree); + + // Second grid portal with compatible phi binning + auto grid2 = GridPortalLink::make( + cyl2, Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5}, + Axis{AxisBound, -50_mm, 50_mm, 5}); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binZ, *logger); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + + const auto& axis1 = *merged->grid().axes().front(); + const auto& axis2 = *merged->grid().axes().back(); + + Axis axis1Expected{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5}; + BOOST_CHECK_EQUAL(axis1, axis1Expected); + Axis axis2Expected{AxisBound, + {-150, -110, -70, -30, 10, 50, 70, 90, 110, 130, 150}}; + BOOST_CHECK_EQUAL(axis2, axis2Expected); + } +} + +BOOST_AUTO_TEST_CASE(RPhiDirection) { + auto cyl1 = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 45_degree); + + // z good, rphi good + auto grid1 = GridPortalLink::make( + cyl1, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5}, + Axis{AxisBound, -100_mm, 100_mm, 5}); + BOOST_REQUIRE_NE(grid1, nullptr); + + auto trf2 = Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}; + auto cyl2 = + Surface::makeShared(trf2, 30_mm, 100_mm, 45_degree); + + // Second grid portal with compatible phi binning + auto grid2 = GridPortalLink::make( + cyl2, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5}, + Axis{AxisBound, -100_mm, 100_mm, 5}); + BOOST_REQUIRE_NE(grid2, nullptr); + + // We're merging in z direction, so the phi binnings need to be the same + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binRPhi, *logger); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(mergedPtr, nullptr); + + const auto& axis1 = *merged->grid().axes().front(); + const auto& axis2 = *merged->grid().axes().back(); + BOOST_CHECK_CLOSE(axis1.getMin(), -90_degree * 30_mm, 1e-8); + BOOST_CHECK_CLOSE(axis1.getMax(), 90_degree * 30_mm, 1e-8); + BOOST_CHECK_EQUAL(axis1.getNBins(), 10); + BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound); + Axis axis2Expected{AxisBound, -100_mm, 100_mm, 5}; + BOOST_CHECK_EQUAL(axis2, axis2Expected); +} + +BOOST_AUTO_TEST_SUITE_END() // Merging2dCylinder + +BOOST_AUTO_TEST_SUITE(Merging1dDisc) + +BOOST_AUTO_TEST_SUITE(RDirection) + +BOOST_AUTO_TEST_CASE(ColinearMerge) { + // Without phi sector + auto disc1 = + Surface::makeShared(Transform3::Identity(), 30_mm, 100_mm); + + auto grid1 = GridPortalLink::make(disc1, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 7}); + + auto disc2 = + Surface::makeShared(Transform3::Identity(), 100_mm, 150_mm); + + auto grid2 = GridPortalLink::make(disc2, BinningValue::binR, + Axis{AxisBound, 100_mm, 150_mm, 5}); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binR, *logger); + BOOST_REQUIRE(mergedPtr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + Axis axisExpected{AxisBound, 30_mm, 150_mm, 12}; + BOOST_CHECK_EQUAL(*merged->grid().axes().front(), axisExpected); + + // With phi sector + auto discPhi1 = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 30_degree); + + auto discPhiGrid1 = GridPortalLink::make(discPhi1, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 7}); + + auto discPhi2 = Surface::makeShared(Transform3::Identity(), + 100_mm, 150_mm, 30_degree); + + auto discPhiGrid2 = GridPortalLink::make(discPhi2, BinningValue::binR, + Axis{AxisBound, 100_mm, 150_mm, 5}); + + auto mergedPhiPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2, + BinningValue::binR, *logger); + BOOST_REQUIRE(mergedPhiPtr); + const auto* mergedPhi = + dynamic_cast(mergedPhiPtr.get()); + BOOST_REQUIRE_NE(mergedPhi, nullptr); + + BOOST_CHECK_EQUAL(mergedPhi->grid().axes().size(), 1); + BOOST_CHECK_EQUAL(*mergedPhi->grid().axes().front(), axisExpected); +} + +BOOST_AUTO_TEST_CASE(ParallelMerge) { + auto disc1 = + Surface::makeShared(Transform3::Identity(), 30_mm, 100_mm); + + auto grid1 = + GridPortalLink::make(disc1, BinningValue::binPhi, + Axis{AxisClosed, -180_degree, 180_degree, 5}); + + auto disc2 = + Surface::makeShared(Transform3::Identity(), 100_mm, 150_mm); + + auto grid2 = + GridPortalLink::make(disc2, BinningValue::binPhi, + Axis{AxisClosed, -180_degree, 180_degree, 5}); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binR, *logger); + BOOST_REQUIRE(mergedPtr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2); + const auto& axis1 = *merged->grid().axes().front(); + const auto& axis2 = *merged->grid().axes().back(); + Axis axis1Expected{AxisBound, {30_mm, 100_mm, 150_mm}}; + BOOST_CHECK_EQUAL(axis1, axis1Expected); + Axis axis2Expected{AxisClosed, -180_degree, 180_degree, 5}; + BOOST_CHECK_EQUAL(axis2, axis2Expected); +} + +BOOST_AUTO_TEST_SUITE_END() // RDirection + +BOOST_AUTO_TEST_SUITE(PhiDirection) + +BOOST_AUTO_TEST_CASE(ColinearMerge) { + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm, 30_degree); + + auto grid1 = GridPortalLink::make(disc1, BinningValue::binPhi, + Axis{AxisBound, -30_degree, 30_degree, 3}); + + auto disc2 = Surface::makeShared( + Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 60_degree); + + auto grid2 = GridPortalLink::make(disc2, BinningValue::binPhi, + Axis{AxisBound, -60_degree, 60_degree, 6}); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binPhi, *logger); + BOOST_REQUIRE(mergedPtr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1); + const auto& axis = *merged->grid().axes().front(); + BOOST_CHECK_CLOSE(axis.getMin(), -90_degree, 1e-6); + BOOST_CHECK_CLOSE(axis.getMax(), 90_degree, 1e-6); + BOOST_CHECK_EQUAL(axis.getNBins(), 9); + + // Check wrapping + + auto disc1Half = Surface::makeShared( + Transform3{AngleAxis3{15_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 90_degree); + + auto grid1Half = + GridPortalLink::make(disc1Half, BinningValue::binPhi, + Axis{AxisBound, -90_degree, 90_degree, 3}); + + auto disc2Half = Surface::makeShared( + Transform3{AngleAxis3{-165_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 90_degree); + + auto grid2Half = + GridPortalLink::make(disc2Half, BinningValue::binPhi, + Axis{AxisBound, -90_degree, 90_degree, 3}); + + auto mergedHalfPtr = GridPortalLink::merge(*grid1Half, *grid2Half, + BinningValue::binPhi, *logger); + BOOST_REQUIRE(mergedHalfPtr); + const auto* mergedHalf = + dynamic_cast(mergedHalfPtr.get()); + BOOST_REQUIRE_NE(mergedHalf, nullptr); + + BOOST_CHECK_EQUAL(mergedHalf->grid().axes().size(), 1); + Axis axisHalfExpected{AxisClosed, -180_degree, 180_degree, 6}; + BOOST_CHECK_EQUAL(axisHalfExpected, *mergedHalf->grid().axes().front()); +} + +BOOST_AUTO_TEST_CASE(ParallelMerge) { + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm, 30_degree); + + auto grid1 = GridPortalLink::make(disc1, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 3}); + + auto disc2 = Surface::makeShared( + Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 60_degree); + + auto grid2 = GridPortalLink::make(disc2, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 3}); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binPhi, *logger); + BOOST_REQUIRE(mergedPtr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2); + const auto& axis1 = *merged->grid().axes().front(); + const auto& axis2 = *merged->grid().axes().back(); + + Axis axis1Expected{AxisBound, 30_mm, 100_mm, 3}; + BOOST_CHECK_EQUAL(axis1, axis1Expected); + + BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6); + BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6); + BOOST_CHECK_EQUAL(axis2.getNBins(), 2); + BOOST_CHECK_CLOSE(axis2.getBinEdges().at(1), 30_degree, 1e-6); + BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound); +} + +BOOST_AUTO_TEST_SUITE_END() // PhiDirection + +BOOST_AUTO_TEST_CASE(BinFilling) { + // Volumes for bin content checking + + auto vol1 = makeDummyVolume(); + auto vol2 = makeDummyVolume(); + + BOOST_TEST_CONTEXT("RDirection") { + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 60_mm, 30_degree); + + auto grid1 = GridPortalLink::make(disc1, BinningValue::binR, + Axis{AxisBound, 30_mm, 60_mm, 2}); + + grid1->setVolume(vol1.get()); + + auto disc2 = Surface::makeShared(Transform3::Identity(), 60_mm, + 90_mm, 30_degree); + + auto grid2 = GridPortalLink::make(disc2, BinningValue::binR, + Axis{AxisBound, 60_mm, 90_mm, 2}); + + grid2->setVolume(vol2.get()); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binR, *logger); + + using merged_type = + GridPortalLinkT>; + + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE(merged); + + grid1->printContents(std::cout); + grid2->printContents(std::cout); + merged->printContents(std::cout); + + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({1}), vol1.get()); + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({2}), vol1.get()); + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({3}), vol2.get()); + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({4}), vol2.get()); + } + + BOOST_TEST_CONTEXT("PhiDirection") { + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm, 30_degree); + + auto grid1 = GridPortalLink::make( + disc1, BinningValue::binPhi, Axis{AxisBound, -30_degree, 30_degree, 2}); + + grid1->setVolume(vol1.get()); + + auto disc2 = Surface::makeShared( + Transform3{AngleAxis3{60_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 30_degree); + + auto grid2 = GridPortalLink::make( + disc2, BinningValue::binPhi, Axis{AxisBound, -30_degree, 30_degree, 2}); + + grid2->setVolume(vol2.get()); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binPhi, *logger); + BOOST_REQUIRE(mergedPtr); + + using merged_type = + GridPortalLinkT>; + + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE(merged); + + grid1->printContents(std::cout); + grid2->printContents(std::cout); + merged->printContents(std::cout); + + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({1}), vol2.get()); + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({2}), vol2.get()); + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({3}), vol1.get()); + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({4}), vol1.get()); + } +} + +BOOST_AUTO_TEST_SUITE_END() // Merging1dDisc + +BOOST_AUTO_TEST_SUITE(Merging2dDisc) + +BOOST_AUTO_TEST_CASE(RDirection) { + // Basic, because the parallel 1D case already tests this to some degree + auto discPhi1 = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 30_degree); + + auto discPhiGrid1 = + GridPortalLink::make(discPhi1, Axis{AxisBound, 30_mm, 100_mm, 7}, + Axis{AxisBound, -30_degree, 30_degree, 3}); + + auto discPhi2 = Surface::makeShared(Transform3::Identity(), + 100_mm, 150_mm, 30_degree); + + auto discPhiGrid2 = + GridPortalLink::make(discPhi2, Axis{AxisBound, 100_mm, 150_mm, 5}, + Axis{AxisBound, -30_degree, 30_degree, 3}); + + auto mergedPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2, + BinningValue::binR, *logger); + BOOST_REQUIRE(mergedPtr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2); + const auto& axis1 = *merged->grid().axes().front(); + const auto& axis2 = *merged->grid().axes().back(); + + BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm); + BOOST_CHECK_EQUAL(axis1.getMax(), 150_mm); + BOOST_CHECK_EQUAL(axis1.getNBins(), 12); + BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound); + BOOST_CHECK_EQUAL(axis2.getMin(), -30_degree); + BOOST_CHECK_EQUAL(axis2.getMax(), 30_degree); + BOOST_CHECK_EQUAL(axis2.getNBins(), 3); + BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant); +} + +BOOST_AUTO_TEST_CASE(PhiDirection) { + // Basic, because the parallel 1D case already tests this to some degree + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm, 30_degree); + + auto grid1 = GridPortalLink::make(disc1, Axis{AxisBound, 30_mm, 100_mm, 3}, + Axis{AxisBound, -30_degree, 30_degree, 3}); + + auto disc2 = Surface::makeShared( + Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 60_degree); + + auto grid2 = GridPortalLink::make(disc2, Axis{AxisBound, 30_mm, 100_mm, 3}, + Axis{AxisBound, -60_degree, 60_degree, 6}); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binPhi, *logger); + BOOST_REQUIRE(mergedPtr); + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2); + const auto& axis1 = *merged->grid().axes().front(); + const auto& axis2 = *merged->grid().axes().back(); + + BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm); + BOOST_CHECK_EQUAL(axis1.getMax(), 100_mm); + BOOST_CHECK_EQUAL(axis1.getNBins(), 3); + BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound); + BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6); + BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6); + BOOST_CHECK_EQUAL(axis2.getNBins(), 9); + BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound); +} + +BOOST_AUTO_TEST_CASE(BinFilling) { + // Volumes for bin content checking + // Volume shape/transform is irrelevant, only used for pointer identity + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto fillCheckerBoard = [&](auto& grid) { + auto loc = grid.numLocalBins(); + for (std::size_t i = 1; i <= loc[0]; ++i) { + for (std::size_t j = 1; j <= loc[1]; ++j) { + grid.atLocalBins({i, j}) = (i + j) % 2 == 0 ? vol1.get() : vol2.get(); + } + } + }; + + auto checkCheckerBoard = [&](const auto& grid) { + auto loc = grid.numLocalBins(); + for (std::size_t i = 1; i <= loc[0]; ++i) { + for (std::size_t j = 1; j <= loc[1]; ++j) { + const auto* vol = grid.atLocalBins({i, j}); + if (vol != ((i + j) % 2 == 0 ? vol1.get() : vol2.get())) { + BOOST_ERROR("Is not a checkerboard pattern"); + return; + } + } + } + }; + + BOOST_TEST_CONTEXT("RDirection") { + auto discPhi1 = Surface::makeShared(Transform3::Identity(), + 30_mm, 60_mm, 30_degree); + + auto discPhiGrid1 = + GridPortalLink::make(discPhi1, Axis{AxisBound, 30_mm, 60_mm, 2}, + Axis{AxisBound, -30_degree, 30_degree, 2}); + + fillCheckerBoard(discPhiGrid1->grid()); + checkCheckerBoard(discPhiGrid1->grid()); + + auto discPhi2 = Surface::makeShared(Transform3::Identity(), + 60_mm, 90_mm, 30_degree); + + auto discPhiGrid2 = + GridPortalLink::make(discPhi2, Axis{AxisBound, 60_mm, 90_mm, 2}, + Axis{AxisBound, -30_degree, 30_degree, 2}); + + fillCheckerBoard(discPhiGrid2->grid()); + checkCheckerBoard(discPhiGrid2->grid()); + + auto mergedPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2, + BinningValue::binR, *logger); + + using merged_type = + GridPortalLinkT, + Axis>; + + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE(merged); + checkCheckerBoard(merged->grid()); + + // Fill a / b + discPhiGrid1->setVolume(vol1.get()); + discPhiGrid2->setVolume(vol2.get()); + + mergedPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2, + BinningValue::binR, *logger); + + merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE(merged); + + const auto* v1 = vol1.get(); + const auto* v2 = vol2.get(); + + std::vector> locations = { + {{40_mm, -20_degree}, v1}, {{40_mm, 20_degree}, v1}, + {{50_mm, -20_degree}, v1}, {{50_mm, 20_degree}, v1}, + {{70_mm, -20_degree}, v2}, {{70_mm, 20_degree}, v2}, + {{80_mm, -20_degree}, v2}, {{80_mm, 20_degree}, v2}, + }; + + for (const auto& [loc, vol] : locations) { + BOOST_TEST_CONTEXT(loc.transpose()) + BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol); + } + + std::vector> contents = { + {v1, v1}, + {v1, v1}, + {v2, v2}, + {v2, v2}, + }; + + for (std::size_t i = 0; i < 4; ++i) { + for (std::size_t j = 0; j < 2; ++j) { + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({i + 1, j + 1}), + contents.at(i).at(j)); + } + } + } + + BOOST_TEST_CONTEXT("PhiDirection") { + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm, 30_degree); + + auto grid1 = + GridPortalLink::make(disc1, Axis{AxisBound, 30_mm, 100_mm, 2}, + Axis{AxisBound, -30_degree, 30_degree, 2}); + fillCheckerBoard(grid1->grid()); + checkCheckerBoard(grid1->grid()); + + auto disc2 = Surface::makeShared( + Transform3{AngleAxis3{60_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 30_degree); + + auto grid2 = + GridPortalLink::make(disc2, Axis{AxisBound, 30_mm, 100_mm, 2}, + Axis{AxisBound, -30_degree, 30_degree, 2}); + fillCheckerBoard(grid2->grid()); + checkCheckerBoard(grid2->grid()); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binPhi, *logger); + + BOOST_REQUIRE(mergedPtr); + + using merged_type = + GridPortalLinkT, + Axis>; + + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE(merged); + + checkCheckerBoard(merged->grid()); + + // Fill a / b + grid1->setVolume(vol1.get()); + grid2->setVolume(vol2.get()); + + mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binPhi, *logger); + merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE(merged); + + const auto* v1 = vol1.get(); + const auto* v2 = vol2.get(); + + std::vector> locations = { + {{40_mm, -50_degree}, v2}, {{40_mm, -10_degree}, v2}, + {{50_mm, -50_degree}, v2}, {{50_mm, -10_degree}, v2}, + {{40_mm, 10_degree}, v1}, {{50_mm, 50_degree}, v1}, + {{50_mm, 10_degree}, v1}, {{50_mm, 50_degree}, v1}, + }; + + for (const auto& [loc, vol] : locations) { + BOOST_TEST_CONTEXT(loc.transpose()) + BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol); + } + + std::vector> contents = { + {v2, v2, v1, v1}, + {v2, v2, v1, v1}, + }; + + for (std::size_t i = 0; i < 2; ++i) { + for (std::size_t j = 0; j < 4; ++j) { + BOOST_CHECK_EQUAL(merged->grid().atLocalBins({i + 1, j + 1}), + contents.at(i).at(j)); + } + } + } +} + +BOOST_AUTO_TEST_SUITE_END() // Merging2dDisc + +BOOST_AUTO_TEST_SUITE(MergingMixedDisc) + +BOOST_AUTO_TEST_CASE(RDirection) { + auto discPhi1 = Surface::makeShared(Transform3::Identity(), + 30_mm, 100_mm, 30_degree); + + auto discPhiGrid1 = + GridPortalLink::make(discPhi1, Axis{AxisBound, 30_mm, 100_mm, 7}, + Axis{AxisBound, -30_degree, 30_degree, 3}); + + auto discPhi2 = Surface::makeShared(Transform3::Identity(), + 100_mm, 150_mm, 30_degree); + + auto discPhiGrid21dPhi = + GridPortalLink::make(discPhi2, BinningValue::binPhi, + Axis{AxisBound, -30_degree, 30_degree, 3}); + + auto merged12PhiPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid21dPhi, + BinningValue::binR, *logger); + BOOST_REQUIRE(merged12PhiPtr); + const auto* merged12Phi = + dynamic_cast(merged12PhiPtr.get()); + BOOST_REQUIRE_NE(merged12Phi, nullptr); + + auto merged21PhiPtr = GridPortalLink::merge(*discPhiGrid21dPhi, *discPhiGrid1, + BinningValue::binR, *logger); + BOOST_REQUIRE(merged21PhiPtr); + const auto* merged21Phi = + dynamic_cast(merged21PhiPtr.get()); + BOOST_REQUIRE_NE(merged21Phi, nullptr); + + BOOST_CHECK_EQUAL(merged12Phi->grid(), merged21Phi->grid()); + + BOOST_CHECK_EQUAL(merged12Phi->grid().axes().size(), 2); + const auto& axis1 = *merged12Phi->grid().axes().front(); + const auto& axis2 = *merged12Phi->grid().axes().back(); + + BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm); + BOOST_CHECK_EQUAL(axis1.getMax(), 150_mm); + BOOST_CHECK_EQUAL(axis1.getNBins(), 8); + BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound); + BOOST_CHECK_EQUAL(axis2.getMin(), -30_degree); + BOOST_CHECK_EQUAL(axis2.getMax(), 30_degree); + BOOST_CHECK_EQUAL(axis2.getNBins(), 3); + BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound); + + auto discPhiGrid21dR = GridPortalLink::make( + discPhi2, BinningValue::binR, Axis{AxisBound, 100_mm, 150_mm, 5}); + + auto merged12RPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid21dR, + BinningValue::binR, *logger); + BOOST_REQUIRE(merged12RPtr); + const auto* merged12R = + dynamic_cast(merged12RPtr.get()); + BOOST_REQUIRE_NE(merged12R, nullptr); + + auto merged21RPtr = GridPortalLink::merge(*discPhiGrid21dR, *discPhiGrid1, + BinningValue::binR, *logger); + BOOST_REQUIRE(merged21RPtr); + const auto* merged21R = + dynamic_cast(merged21RPtr.get()); + BOOST_REQUIRE_NE(merged21R, nullptr); + BOOST_CHECK_EQUAL(merged12R->grid(), merged21R->grid()); + + BOOST_CHECK_EQUAL(merged12R->grid().axes().size(), 2); + const auto& axis1R = *merged12R->grid().axes().front(); + const auto& axis2R = *merged12R->grid().axes().back(); + + BOOST_CHECK_EQUAL(axis1R.getMin(), 30_mm); + BOOST_CHECK_EQUAL(axis1R.getMax(), 150_mm); + BOOST_CHECK_EQUAL(axis1R.getNBins(), 12); + BOOST_CHECK_EQUAL(axis1R.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis1R.getBoundaryType(), AxisBoundaryType::Bound); + BOOST_CHECK_EQUAL(axis2R.getMin(), -30_degree); + BOOST_CHECK_EQUAL(axis2R.getMax(), 30_degree); + BOOST_CHECK_EQUAL(axis2R.getNBins(), 3); + BOOST_CHECK_EQUAL(axis2R.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis2R.getBoundaryType(), AxisBoundaryType::Bound); +} + +BOOST_AUTO_TEST_CASE(PhiDirection) { + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm, 30_degree); + + auto grid1 = GridPortalLink::make(disc1, Axis{AxisBound, 30_mm, 100_mm, 3}, + Axis{AxisBound, -30_degree, 30_degree, 3}); + + auto disc2 = Surface::makeShared( + Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 60_degree); + + auto grid21dPhi = + GridPortalLink::make(disc2, BinningValue::binPhi, + + Axis{AxisBound, -60_degree, 60_degree, 6}); + + auto merged12PhiPtr = + GridPortalLink::merge(*grid1, *grid21dPhi, BinningValue::binPhi, *logger); + BOOST_REQUIRE(merged12PhiPtr); + const auto* merged12Phi = + dynamic_cast(merged12PhiPtr.get()); + BOOST_REQUIRE_NE(merged12Phi, nullptr); + + BOOST_CHECK_EQUAL(merged12Phi->grid().axes().size(), 2); + const auto& axis1 = *merged12Phi->grid().axes().front(); + const auto& axis2 = *merged12Phi->grid().axes().back(); + + BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm); + BOOST_CHECK_EQUAL(axis1.getMax(), 100_mm); + BOOST_CHECK_EQUAL(axis1.getNBins(), 3); + BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound); + BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6); + BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6); + BOOST_CHECK_EQUAL(axis2.getNBins(), 9); + BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound); + + auto grid21dR = GridPortalLink::make(disc2, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 3}); + + auto merged12RPtr = + GridPortalLink::merge(*grid1, *grid21dR, BinningValue::binPhi, *logger); + BOOST_REQUIRE(merged12RPtr); + const auto* merged12R = + dynamic_cast(merged12RPtr.get()); + BOOST_REQUIRE_NE(merged12R, nullptr); + + BOOST_CHECK_EQUAL(merged12R->grid().axes().size(), 2); + const auto& axis1R = *merged12R->grid().axes().front(); + const auto& axis2R = *merged12R->grid().axes().back(); + + BOOST_CHECK_EQUAL(axis1R.getMin(), 30_mm); + BOOST_CHECK_EQUAL(axis1R.getMax(), 100_mm); + BOOST_CHECK_EQUAL(axis1R.getNBins(), 3); + BOOST_CHECK_EQUAL(axis1R.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis1R.getBoundaryType(), AxisBoundaryType::Bound); + BOOST_CHECK_CLOSE(axis2R.getMin(), -90_degree, 1e-6); + BOOST_CHECK_CLOSE(axis2R.getMax(), 90_degree, 1e-6); + BOOST_CHECK_EQUAL(axis2R.getNBins(), 4); + BOOST_CHECK_EQUAL(axis2R.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis2R.getBoundaryType(), AxisBoundaryType::Bound); +} + +BOOST_AUTO_TEST_SUITE_END() // MergingMixedDisc + +BOOST_AUTO_TEST_SUITE(MergingCrossDisc) + +BOOST_AUTO_TEST_CASE(RDirection) { + // Volumes for bin content checking + auto vol1 = makeDummyVolume(); + auto vol2 = makeDummyVolume(); + auto vol3 = makeDummyVolume(); + auto vol4 = makeDummyVolume(); + + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm, 30_degree); + + auto grid1 = GridPortalLink::make(disc1, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 2}); + grid1->grid().atLocalBins({1}) = vol1.get(); + grid1->grid().atLocalBins({2}) = vol2.get(); + + auto disc2 = Surface::makeShared(Transform3::Identity(), 100_mm, + 150_mm, 30_degree); + + auto grid2 = GridPortalLink::make(disc2, BinningValue::binPhi, + Axis{AxisBound, -30_degree, 30_degree, 2}); + + grid2->grid().atLocalBins({1}) = vol3.get(); + grid2->grid().atLocalBins({2}) = vol4.get(); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binR, *logger); + + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2); + const auto& axis1 = *merged->grid().axes().front(); + const auto& axis2 = *merged->grid().axes().back(); + Axis axisExpectedR{AxisBound, {30_mm, 65_mm, 100_mm, 150_mm}}; + Axis axisExpectedPhi{AxisBound, -30_degree, 30_degree, 2}; + BOOST_CHECK_EQUAL(axis1, axisExpectedR); + BOOST_CHECK_EQUAL(axis2, axisExpectedPhi); + + std::vector> locations = { + {{40_mm, -15_degree}, vol1.get()}, {{40_mm, 15_degree}, vol1.get()}, + {{90_mm, -15_degree}, vol2.get()}, {{90_mm, 15_degree}, vol2.get()}, + + {{110_mm, -15_degree}, vol3.get()}, {{110_mm, 15_degree}, vol4.get()}, + {{140_mm, -15_degree}, vol3.get()}, {{140_mm, 15_degree}, vol4.get()}, + }; + + for (const auto& [loc, vol] : locations) { + BOOST_TEST_CONTEXT(loc.transpose()) + BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol); + } + + grid1->printContents(std::cout); + grid2->printContents(std::cout); + merged->printContents(std::cout); +} + +BOOST_AUTO_TEST_CASE(PhiDirection) { + // Volumes for bin content checking + + auto vol1 = makeDummyVolume(); + auto vol2 = makeDummyVolume(); + auto vol3 = makeDummyVolume(); + auto vol4 = makeDummyVolume(); + + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm, 30_degree); + + auto grid1 = GridPortalLink::make(disc1, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 2}); + + grid1->grid().atLocalBins({1}) = vol1.get(); + grid1->grid().atLocalBins({2}) = vol2.get(); + + auto disc2 = Surface::makeShared( + Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 60_degree); + + auto grid2 = GridPortalLink::make(disc2, BinningValue::binPhi, + Axis{AxisBound, -60_degree, 60_degree, 2}); + + grid2->grid().atLocalBins({1}) = vol3.get(); + grid2->grid().atLocalBins({2}) = vol4.get(); + + auto mergedPtr = + GridPortalLink::merge(*grid1, *grid2, BinningValue::binPhi, *logger); + + using merged_type = + GridPortalLinkT, + Axis>; + + const auto* merged = dynamic_cast(mergedPtr.get()); + BOOST_REQUIRE_NE(merged, nullptr); + + BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2); + const auto& axis1 = *merged->grid().axes().front(); + const auto& axis2 = *merged->grid().axes().back(); + Axis axisExpectedR{AxisBound, 30_mm, 100_mm, 2}; + BOOST_CHECK_EQUAL(axis1, axisExpectedR); + + BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6); + BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6); + BOOST_CHECK_EQUAL(axis2.getNBins(), 3); + BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant); + BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound); + + std::vector> locations = { + {{40_mm, 45_degree}, vol1.get()}, {{40_mm, 0_degree}, vol4.get()}, + {{40_mm, -80_degree}, vol3.get()}, {{90_mm, 45_degree}, vol2.get()}, + {{90_mm, 0_degree}, vol4.get()}, {{90_mm, -80_degree}, vol3.get()}, + }; + + grid1->printContents(std::cout); + grid2->printContents(std::cout); + merged->printContents(std::cout); + + for (const auto& [loc, vol] : locations) { + BOOST_TEST_CONTEXT((Vector2{loc[0], loc[1] / 1_degree}.transpose())) + BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol); + } +} + +BOOST_AUTO_TEST_SUITE_END() // MergeCrossDisc + +BOOST_AUTO_TEST_SUITE_END() // GridMerging + +BOOST_AUTO_TEST_CASE(CompositeConstruction) { + // Arbitrary volumes for testing only + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol3 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto disc1 = + Surface::makeShared(Transform3::Identity(), 30_mm, 60_mm); + + auto trivial1 = std::make_unique(disc1, *vol1); + + auto disc2 = + Surface::makeShared(Transform3::Identity(), 60_mm, 90_mm); + auto trivial2 = std::make_unique(disc2, *vol2); + + auto composite = std::make_unique( + copy(trivial1), copy(trivial2), BinningValue::binR); + + auto compositeCopy = std::make_unique( + copy(trivial1), copy(trivial2), BinningValue::binR); + + BOOST_CHECK_EQUAL( + composite->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(), + vol1.get()); + BOOST_CHECK_EQUAL( + composite->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(), + vol2.get()); + + BOOST_CHECK_EQUAL(composite->depth(), 1); + + // Test exception on different surface types + auto cyl = Surface::makeShared(Transform3::Identity(), 30_mm, + 40_mm); + auto trivialCyl = std::make_unique(cyl, *vol3); + BOOST_CHECK_THROW( + CompositePortalLink(copy(trivial1), copy(trivialCyl), BinningValue::binR), + std::invalid_argument); + + auto disc3 = + Surface::makeShared(Transform3::Identity(), 90_mm, 120_mm); + auto trivial3 = std::make_unique(disc3, *vol3); + + // Test exception on un-mergable surfaces + BOOST_CHECK_THROW( + CompositePortalLink(copy(trivial1), copy(trivial3), BinningValue::binR), + SurfaceMergingException); + + // Composite with a composite (this should work regardless of flattening) + + CompositePortalLink composite2(std::move(composite), copy(trivial3), + BinningValue::binR, false); + + BOOST_CHECK_EQUAL( + composite2.resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(), + vol1.get()); + BOOST_CHECK_EQUAL( + composite2.resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(), + vol2.get()); + BOOST_CHECK_EQUAL( + composite2.resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(), + vol3.get()); + + // Two without flattening + BOOST_CHECK_EQUAL(composite2.depth(), 2); + + CompositePortalLink composite2Flat(std::move(compositeCopy), copy(trivial3), + BinningValue::binR, true); + + // One because of flattening + BOOST_CHECK_EQUAL(composite2Flat.depth(), 1); + + BOOST_CHECK_EQUAL( + composite2Flat.resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(), + vol1.get()); + BOOST_CHECK_EQUAL( + composite2Flat.resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(), + vol2.get()); + BOOST_CHECK_EQUAL( + composite2Flat.resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(), + vol3.get()); +} + +BOOST_AUTO_TEST_SUITE(PortalMerging) + +BOOST_AUTO_TEST_CASE(TrivialTrivial) {} + +BOOST_AUTO_TEST_CASE(TrivialGridR) { + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto disc1 = + Surface::makeShared(Transform3::Identity(), 30_mm, 60_mm); + + auto disc2 = + Surface::makeShared(Transform3::Identity(), 60_mm, 90_mm); + + auto trivial = std::make_unique(disc2, *vol2); + BOOST_REQUIRE(trivial); + + auto gridPhi = GridPortalLink::make(disc1, BinningValue::binPhi, + Axis{AxisClosed, -M_PI, M_PI, 2}); + gridPhi->setVolume(vol1.get()); + + auto gridR = GridPortalLink::make(disc1, BinningValue::binR, + Axis{AxisBound, 30_mm, 60_mm, 2}); + gridR->setVolume(vol1.get()); + + BOOST_TEST_CONTEXT("Colinear") { + auto merged = PortalLinkBase::merge(copy(trivial), copy(gridR), + BinningValue::binR, *logger); + BOOST_REQUIRE(merged); + + auto* mergedGrid = dynamic_cast(merged.get()); + BOOST_REQUIRE(mergedGrid); + + mergedGrid->printContents(std::cout); + + BOOST_CHECK_EQUAL(mergedGrid->grid().axes().size(), 1); + Axis axisExpected{AxisBound, {30_mm, 45_mm, 60_mm, 90_mm}}; + BOOST_CHECK_EQUAL(*mergedGrid->grid().axes().front(), axisExpected); + } + + BOOST_TEST_CONTEXT("Orthogonal") { + auto merged = PortalLinkBase::merge(copy(gridPhi), copy(trivial), + BinningValue::binR, *logger); + BOOST_REQUIRE(merged); + + auto* mergedGrid = dynamic_cast(merged.get()); + BOOST_REQUIRE(mergedGrid); + + mergedGrid->printContents(std::cout); + + BOOST_CHECK_EQUAL(mergedGrid->grid().axes().size(), 2); + Axis axis1Expected{AxisBound, 30_mm, 90_mm, 2}; + Axis axis2Expected{AxisClosed, -M_PI, M_PI, 2}; + BOOST_CHECK_EQUAL(*mergedGrid->grid().axes().front(), axis1Expected); + BOOST_CHECK_EQUAL(*mergedGrid->grid().axes().back(), axis2Expected); + } +} + +BOOST_AUTO_TEST_CASE(TrivialGridPhi) { + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto disc1 = Surface::makeShared(Transform3::Identity(), 30_mm, + 100_mm, 30_degree); + + auto disc2 = Surface::makeShared( + Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm, + 60_degree); + + auto trivial = std::make_unique(disc2, *vol2); + BOOST_REQUIRE(trivial); + + auto gridPhi = GridPortalLink::make( + disc1, BinningValue::binPhi, Axis{AxisBound, -30_degree, 30_degree, 2}); + gridPhi->setVolume(vol1.get()); + + auto gridR = GridPortalLink::make(disc1, BinningValue::binR, + Axis{AxisBound, 30_mm, 100_mm, 2}); + gridR->setVolume(vol1.get()); + + BOOST_TEST_CONTEXT("Colinear") { + auto merged = PortalLinkBase::merge(copy(trivial), copy(gridPhi), + BinningValue::binPhi, *logger); + BOOST_REQUIRE(merged); + + auto* mergedGrid = dynamic_cast(merged.get()); + BOOST_REQUIRE(mergedGrid); + + mergedGrid->printContents(std::cout); + + BOOST_CHECK_EQUAL(mergedGrid->grid().axes().size(), 1); + const auto& axis = *mergedGrid->grid().axes().front(); + BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); + BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); + std::vector expectedBins{-90_degree, 30_degree, 60_degree, + 90_degree}; + CHECK_CLOSE_REL(axis.getBinEdges(), expectedBins, 1e-6); + } + + BOOST_TEST_CONTEXT("Orthogonal") { + auto merged = PortalLinkBase::merge(copy(gridR), copy(trivial), + BinningValue::binPhi, *logger); + BOOST_REQUIRE(merged); + + auto* mergedGrid = dynamic_cast(merged.get()); + BOOST_REQUIRE(mergedGrid); + + mergedGrid->printContents(std::cout); + + BOOST_CHECK_EQUAL(mergedGrid->grid().axes().size(), 2); + const auto& axis1 = *mergedGrid->grid().axes().front(); + const auto& axis2 = *mergedGrid->grid().axes().back(); + Axis axis1Expected{AxisBound, 30_mm, 100_mm, 2}; + BOOST_CHECK_EQUAL(axis1, axis1Expected); + std::vector expectedBins{-90_degree, 30_degree, 90_degree}; + CHECK_CLOSE_REL(axis2.getBinEdges(), expectedBins, 1e-6); + } +} + +BOOST_AUTO_TEST_CASE(CompositeOther) { + auto vol1 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol3 = std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto disc1 = + Surface::makeShared(Transform3::Identity(), 30_mm, 60_mm); + auto disc2 = + Surface::makeShared(Transform3::Identity(), 60_mm, 90_mm); + auto disc3 = + Surface::makeShared(Transform3::Identity(), 90_mm, 120_mm); + + auto grid1 = GridPortalLink::make(disc1, *vol1, BinningValue::binR); + auto trivial2 = std::make_unique(disc2, *vol2); + + auto composite12 = std::make_unique( + std::move(grid1), std::move(trivial2), BinningValue::binR); + + BOOST_CHECK_EQUAL(composite12->depth(), 1); + BOOST_CHECK_EQUAL(composite12->size(), 2); + + auto trivial3 = std::make_unique(disc3, *vol3); + + auto composite123Ptr = PortalLinkBase::merge( + std::move(composite12), std::move(trivial3), BinningValue::binR, *logger); + + const auto* composite123 = + dynamic_cast(composite123Ptr.get()); + BOOST_REQUIRE(composite123); + + BOOST_CHECK_EQUAL(composite123->depth(), 1); + + BOOST_CHECK_EQUAL( + composite123->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(), + vol1.get()); + BOOST_CHECK_EQUAL( + composite123->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(), + vol2.get()); + BOOST_CHECK_EQUAL( + composite123->resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(), + vol3.get()); + + BOOST_CHECK_EQUAL(composite123->size(), 3); +} + +BOOST_AUTO_TEST_SUITE_END() // PortalMerging + +BOOST_AUTO_TEST_SUITE_END() // Geometry + +} // namespace Acts::Test diff --git a/Tests/UnitTests/Core/Utilities/AxesTests.cpp b/Tests/UnitTests/Core/Utilities/AxesTests.cpp index b155bec7a43..6ec592746d4 100644 --- a/Tests/UnitTests/Core/Utilities/AxesTests.cpp +++ b/Tests/UnitTests/Core/Utilities/AxesTests.cpp @@ -438,6 +438,10 @@ BOOST_AUTO_TEST_CASE(AxisVisit) { BOOST_CHECK( (std::is_same_v, Axis>)); }); + + std::vector edges = + varClosed.visit([](const auto& axis) { return axis.getBinEdges(); }); + BOOST_CHECK_EQUAL(edges.size(), varClosed.getBinEdges().size()); } BOOST_AUTO_TEST_CASE(Output) {