src/Utilities/MultiAxisSpec.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Utilities/MultiAxisSpec.hpp"
0010
0011 #include "Acts/Surfaces/CylinderBounds.hpp"
0012 #include "Acts/Surfaces/RadialBounds.hpp"
0013 #include "Acts/Surfaces/RectangleBounds.hpp"
0014 #include "Acts/Surfaces/Surface.hpp"
0015 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0016
0017 #include <algorithm>
0018 #include <memory>
0019 #include <ostream>
0020 #include <sstream>
0021 #include <stdexcept>
0022
0023 namespace Acts {
0024
0025 namespace {
0026
0027 std::string unsupportedSurfaceMessage(const Surface& surface) {
0028   std::stringstream ss;
0029   ss << "Axis resolution is not supported for this surface:\n"
0030      << surface.toStream(GeometryContext::dangerouslyDefaultConstruct());
0031   return ss.str();
0032 }
0033
0034 AxisSpec::Options resolveCylinder(const CylinderBounds& cBounds,
0035                                   AxisDirection aDir) {
0036   using enum AxisDirection;
0037
0038   double r = cBounds.get(CylinderBounds::eR);
0039   double hz = cBounds.get(CylinderBounds::eHalfLengthZ);
0040   double avgPhi = cBounds.get(CylinderBounds::eAveragePhi);
0041   double halfPhi = cBounds.get(CylinderBounds::eHalfPhiSector);
0042   AxisBoundaryType phiBoundaryType = cBounds.coversFullAzimuth()
0043                                          ? AxisBoundaryType::Closed
0044                                          : AxisBoundaryType::Bound;
0045   switch (aDir) {
0046     case AxisRPhi:
0047       return {.min = r * (avgPhi - halfPhi),
0048               .max = r * (avgPhi + halfPhi),
0049               .boundaryType = phiBoundaryType};
0050     case AxisPhi:
0051       return {.min = avgPhi - halfPhi,
0052               .max = avgPhi + halfPhi,
0053               .boundaryType = phiBoundaryType};
0054     case AxisZ:
0055       return {.min = -hz, .max = hz, .boundaryType = AxisBoundaryType::Bound};
0056     default:
0057       throw std::invalid_argument(
0058           "Cylinder axis resolution must be along rphi, phi or z");
0059   }
0060 }
0061
0062 AxisSpec::Options resolveDisc(const RadialBounds& rBounds, AxisDirection aDir) {
0063   using enum AxisBoundaryType;
0064
0065   double avgPhi = rBounds.get(RadialBounds::eAveragePhi);
0066   double halfPhi = rBounds.get(RadialBounds::eHalfPhiSector);
0067   switch (aDir) {
0068     case AxisDirection::AxisR:
0069       return {.min = rBounds.get(RadialBounds::eMinR),
0070               .max = rBounds.get(RadialBounds::eMaxR),
0071               .boundaryType = Bound};
0072     case AxisDirection::AxisPhi:
0073       return {.min = avgPhi - halfPhi,
0074               .max = avgPhi + halfPhi,
0075               .boundaryType = rBounds.coversFullAzimuth() ? Closed : Bound};
0076     default:
0077       throw std::invalid_argument(
0078           "Disc axis resolution must be along r or phi");
0079   }
0080 }
0081
0082 AxisSpec::Options resolveRectangle(const RectangleBounds& pBounds,
0083                                    AxisDirection aDir) {
0084   using enum AxisBoundaryType;
0085
0086   switch (aDir) {
0087     case AxisDirection::AxisX:
0088       return {.min = pBounds.get(RectangleBounds::eMinX),
0089               .max = pBounds.get(RectangleBounds::eMaxX),
0090               .boundaryType = Bound};
0091     case AxisDirection::AxisY:
0092       return {.min = pBounds.get(RectangleBounds::eMinY),
0093               .max = pBounds.get(RectangleBounds::eMaxY),
0094               .boundaryType = Bound};
0095     default:
0096       throw std::invalid_argument(
0097           "Rectangle axis resolution must be along x or y");
0098   }
0099 }
0100
0101 AxisSpec::Options resolveTrapezoid(const TrapezoidBounds& pBounds,
0102                                    AxisDirection aDir) {
0103   using enum AxisBoundaryType;
0104
0105   switch (aDir) {
0106     case AxisDirection::AxisX: {
0107       double halfX = std::max(pBounds.get(TrapezoidBounds::eHalfLengthXnegY),
0108                               pBounds.get(TrapezoidBounds::eHalfLengthXposY));
0109       return {.min = -halfX, .max = halfX, .boundaryType = Bound};
0110     }
0111     case AxisDirection::AxisY: {
0112       double halfY = pBounds.get(TrapezoidBounds::eHalfLengthY);
0113       return {.min = -halfY, .max = halfY, .boundaryType = Bound};
0114     }
0115     default:
0116       throw std::invalid_argument(
0117           "Trapezoid axis resolution must be along x or y");
0118   }
0119 }
0120
0121 AxisSpec::Options resolveAxisOptions(const Surface& surface,
0122                                      AxisDirection aDir) {
0123   const SurfaceBounds& bounds = surface.bounds();
0124   AxisSpec::Options options = [&]() -> AxisSpec::Options {
0125     switch (bounds.type()) {
0126       case SurfaceBounds::eCylinder:
0127         return resolveCylinder(static_cast<const CylinderBounds&>(bounds),
0128                                aDir);
0129       case SurfaceBounds::eDisc:
0130         return resolveDisc(static_cast<const RadialBounds&>(bounds), aDir);
0131       case SurfaceBounds::eRectangle:
0132         return resolveRectangle(static_cast<const RectangleBounds&>(bounds),
0133                                 aDir);
0134       case SurfaceBounds::eTrapezoid:
0135         return resolveTrapezoid(static_cast<const TrapezoidBounds&>(bounds),
0136                                 aDir);
0137       default:
0138         throw std::invalid_argument(unsupportedSurfaceMessage(surface));
0139     }
0140   }();
0141   options.direction = aDir;
0142   return options;
0143 }
0144
0145 }  // namespace
0146
0147 MultiAxisSpec::MultiAxisSpec(std::vector<AxisSpec> axisSpecs)
0148     : m_axisSpecs(std::move(axisSpecs)) {
0149   if (m_axisSpecs.empty()) {
0150     throw std::invalid_argument(
0151         "MultiAxisSpec: at least one axis spec is required");
0152   }
0153 }
0154
0155 std::size_t MultiAxisSpec::size() const {
0156   return m_axisSpecs.size();
0157 }
0158
0159 const AxisSpec& MultiAxisSpec::axisSpec(std::size_t i) const {
0160   return m_axisSpecs.at(i);
0161 }
0162
0163 std::span<const AxisSpec> MultiAxisSpec::axisSpecs() const {
0164   return m_axisSpecs;
0165 }
0166
0167 bool MultiAxisSpec::isDeferred() const {
0168   return std::ranges::any_of(
0169       m_axisSpecs, [](const AxisSpec& af) { return af.isDeferred(); });
0170 }
0171
0172 std::unique_ptr<IMultiAxis> MultiAxisSpec::buildMultiAxis(
0173     const Options& options) const {
0174   return buildMultiAxisImpl(options);
0175 }
0176
0177 std::unique_ptr<IMultiAxis> MultiAxisSpec::buildMultiAxisImpl(
0178     std::span<const AxisSpec::Options> options) const {
0179   if (!options.empty() && options.size() != size()) {
0180     throw std::invalid_argument(
0181         "MultiAxisSpec: either one option set per axis or none at all");
0182   }
0183
0184   std::vector<std::unique_ptr<IAxis>> axes;
0185   axes.reserve(size());
0186   for (std::size_t i = 0; i < size(); ++i) {
0187     axes.push_back(m_axisSpecs[i].buildAxis(
0188         options.empty() ? AxisSpec::Options{} : options[i]));
0189   }
0190
0191   switch (axes.size()) {
0192     case 1:
0193       return IMultiAxis::create(*axes[0]);
0194     case 2:
0195       return IMultiAxis::create(*axes[0], *axes[1]);
0196     case 3:
0197       return IMultiAxis::create(*axes[0], *axes[1], *axes[2]);
0198     default:
0199       throw std::domain_error(
0200           "MultiAxisSpec: multi-axes support at most 3 dimensions");
0201   }
0202 }
0203
0204 std::string MultiAxisSpec::toString() const {
0205   std::stringstream ss;
0206   ss << "MultiAxisSpec: " << size() << " axes [";
0207   for (std::size_t i = 0; i < size(); ++i) {
0208     ss << (i > 0 ? "; " : "") << axisSpec(i);
0209   }
0210   ss << "]";
0211   return ss.str();
0212 }
0213
0214 std::unique_ptr<IMultiAxis2D> resolveMultiAxis(
0215     const MultiAxisSpec2D& multiAxisSpec, const Surface& surface) {
0216   std::span<const AxisSpec> axisSpecs = multiAxisSpec.axisSpecs();
0217   std::array<AxisDirection, 2> canonical = surface.localAxes();
0218
0219   const std::size_t nDirected = std::ranges::count_if(
0220       axisSpecs, [](const AxisSpec& af) { return af.direction().has_value(); });
0221   if (nDirected == 1) {
0222     throw std::invalid_argument(
0223         "resolveMultiAxis: either both or neither of the axes must carry a "
0224         "direction");
0225   }
0226
0227   // Bind a canonical direction to its spec, positionally without stored
0228   // directions and by direction otherwise
0229   auto slot = [&](std::size_t i) -> const AxisSpec& {
0230     if (nDirected == 0) {
0231       return axisSpecs[i];
0232     }
0233     AxisDirection aDir = canonical[i];
0234     auto it = std::ranges::find_if(axisSpecs, [aDir](const AxisSpec& af) {
0235       return af.direction() == aDir;
0236     });
0237     if (it == axisSpecs.end()) {
0238       std::stringstream ss;
0239       ss << "resolveMultiAxis: binning directions do not match the surface "
0240             "axes ("
0241          << axisDirectionName(canonical[0]) << ", "
0242          << axisDirectionName(canonical[1]) << ")";
0243       throw std::invalid_argument(ss.str());
0244     }
0245     return *it;
0246   };
0247
0248   return IMultiAxis::create(
0249       *slot(0).buildAxis(resolveAxisOptions(surface, canonical[0])),
0250       *slot(1).buildAxis(resolveAxisOptions(surface, canonical[1])));
0251 }
0252
0253 }  // namespace Acts