src/Geometry/Portal.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/Geometry/Portal.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Tolerance.hpp"
0013 #include "Acts/Geometry/CompositePortalLink.hpp"
0014 #include "Acts/Geometry/GeometryContext.hpp"
0015 #include "Acts/Geometry/GridPortalLink.hpp"
0016 #include "Acts/Geometry/PortalLinkBase.hpp"
0017 #include "Acts/Geometry/TrivialPortalLink.hpp"
0018 #include "Acts/Surfaces/RegularSurface.hpp"
0019 #include "Acts/Utilities/BinningType.hpp"
0020 #include "Acts/Utilities/Zip.hpp"
0021
0022 #include <cstdlib>
0023 #include <memory>
0024 #include <stdexcept>
0025
0026 namespace Acts {
0027
0028 const char* PortalMergingException::what() const noexcept {
0029   return "Failure to merge portals";
0030 }
0031
0032 const char* PortalFusingException::what() const noexcept {
0033   return "Failure to fuse portals";
0034 }
0035
0036 Portal::Portal(Direction direction, std::unique_ptr<PortalLinkBase> link) {
0037   if (link == nullptr) {
0038     throw std::invalid_argument("Link must not be null");
0039   }
0040
0041   m_surface = link->surfacePtr();
0042
0043   if (direction == Direction::AlongNormal()) {
0044     m_alongNormal = std::move(link);
0045   } else {
0046     m_oppositeNormal = std::move(link);
0047   }
0048 }
0049
0050 Portal::Portal(Direction direction, std::shared_ptr<RegularSurface> surface,
0051                TrackingVolume& volume)
0052     : Portal(direction,
0053              std::make_unique<TrivialPortalLink>(std::move(surface), volume)) {}
0054
0055 Portal::Portal(const GeometryContext& gctx,
0056                std::unique_ptr<PortalLinkBase> alongNormal,
0057                std::unique_ptr<PortalLinkBase> oppositeNormal) {
0058   if (alongNormal == nullptr && oppositeNormal == nullptr) {
0059     throw std::invalid_argument("At least one link must be provided");
0060   }
0061
0062   if (alongNormal != nullptr) {
0063     setLink(gctx, Direction::AlongNormal(), std::move(alongNormal));
0064   }
0065   if (oppositeNormal != nullptr) {
0066     setLink(gctx, Direction::OppositeNormal(), std::move(oppositeNormal));
0067   }
0068 }
0069
0070 Portal::Portal(const GeometryContext& gctx, Arguments&& args) {
0071   if (!args.alongNormal.surface && !args.oppositeNormal.surface) {
0072     throw std::invalid_argument("At least one link must be provided");
0073   }
0074
0075   if (args.alongNormal.surface) {
0076     setLink(gctx, Direction::AlongNormal(),
0077             std::make_unique<TrivialPortalLink>(
0078                 std::move(args.alongNormal.surface), *args.alongNormal.volume));
0079   }
0080   if (args.oppositeNormal.surface) {
0081     setLink(gctx, Direction::OppositeNormal(),
0082             std::make_unique<TrivialPortalLink>(
0083                 std::move(args.oppositeNormal.surface),
0084                 *args.oppositeNormal.volume));
0085   }
0086 }
0087
0088 void Portal::setLink(const GeometryContext& gctx, Direction direction,
0089                      std::unique_ptr<PortalLinkBase> link) {
0090   if (link == nullptr) {
0091     throw std::invalid_argument("Link must not be null");
0092   }
0093
0094   auto& target =
0095       direction == Direction::AlongNormal() ? m_alongNormal : m_oppositeNormal;
0096   const auto& other =
0097       direction == Direction::AlongNormal() ? m_oppositeNormal : m_alongNormal;
0098
0099   // check if surfaces are identical
0100   if (m_surface != nullptr &&
0101       !isSameSurface(gctx, link->surface(), *m_surface)) {
0102     throw PortalFusingException();
0103   }
0104
0105   // check if they both have material but are not the same surface
0106   if (m_surface != nullptr && (m_surface.get() != &link->surface()) &&
0107       link->surface().surfaceMaterial() != nullptr &&
0108       m_surface->surfaceMaterial() != nullptr) {
0109     throw PortalFusingException();
0110   }
0111
0112   target = std::move(link);
0113
0114   if (other == nullptr) {
0115     // We don't have an existing surface, take the one we just got
0116     m_surface = target->surfacePtr();
0117     return;
0118   }
0119
0120   if (target->surface().surfaceMaterial() != nullptr) {
0121     // new link has material: assign that to existing link
0122     m_surface = target->surfacePtr();
0123     other->setSurface(m_surface);
0124   } else {
0125     // none have material, or the existing surface had material: assign the
0126     // existing surface by convention
0127     target->setSurface(m_surface);
0128   }
0129 }
0130
0131 void Portal::setLink(const GeometryContext& gctx, Direction direction,
0132                      std::shared_ptr<RegularSurface> surface,
0133                      TrackingVolume& volume) {
0134   setLink(gctx, direction,
0135           std::make_unique<TrivialPortalLink>(std::move(surface), volume));
0136 }
0137
0138 const PortalLinkBase* Portal::getLink(Direction direction) const {
0139   if (direction == Direction::AlongNormal()) {
0140     return m_alongNormal.get();
0141   } else {
0142     return m_oppositeNormal.get();
0143   }
0144 }
0145
0146 Result<const TrackingVolume*> Portal::resolveVolume(
0147     const GeometryContext& gctx, const Vector3& position,
0148     const Vector3& direction) const {
0149   assert(m_surface != nullptr);
0150   const Vector3 normal = m_surface->normal(gctx, position);
0151   Direction side = Direction::fromScalarZeroAsPositive(normal.dot(direction));
0152
0153   const PortalLinkBase* link = side == Direction::AlongNormal()
0154                                    ? m_alongNormal.get()
0155                                    : m_oppositeNormal.get();
0156
0157   if (link == nullptr) {
0158     // no link is attached in this direction => this is the end of the world as
0159     // we know it. (i feel fine)
0160     return nullptr;
0161   } else {
0162     auto res = link->resolveVolume(gctx, position);
0163     if (!res.ok()) {
0164       return res.error();
0165     }
0166     return *res;
0167   }
0168 }
0169
0170 bool Portal::isValid() const {
0171   return m_alongNormal != nullptr || m_oppositeNormal != nullptr;
0172 }
0173
0174 const RegularSurface& Portal::surface() const {
0175   assert(m_surface != nullptr);
0176   return *m_surface;
0177 }
0178
0179 RegularSurface& Portal::surface() {
0180   assert(m_surface != nullptr);
0181   return *m_surface;
0182 }
0183
0184 Portal Portal::merge(const GeometryContext& gctx, Portal& aPortal,
0185                      Portal& bPortal, AxisDirection direction,
0186                      const Logger& logger) {
0187   ACTS_VERBOSE("Merging two portals along " << direction);
0188
0189   if (&aPortal == &bPortal) {
0190     ACTS_ERROR("Cannot merge a portal with itself");
0191     throw PortalMergingException{};
0192   }
0193
0194   if (aPortal.m_surface->surfaceMaterial() != nullptr ||
0195       bPortal.m_surface->surfaceMaterial() != nullptr) {
0196     ACTS_ERROR("Cannot merge portals with material");
0197     throw PortalMergingException{};
0198   }
0199
0200   std::unique_ptr<PortalLinkBase> mergedAlongNormal = nullptr;
0201   std::unique_ptr<PortalLinkBase> mergedOppositeNormal = nullptr;
0202
0203   bool aHasAlongNormal = aPortal.m_alongNormal != nullptr;
0204   bool aHasOppositeNormal = aPortal.m_oppositeNormal != nullptr;
0205   bool bHasAlongNormal = bPortal.m_alongNormal != nullptr;
0206   bool bHasOppositeNormal = bPortal.m_oppositeNormal != nullptr;
0207
0208   if (aHasAlongNormal != bHasAlongNormal ||
0209       aHasOppositeNormal != bHasOppositeNormal) {
0210     ACTS_ERROR("Portals do not have the same links attached");
0211     throw PortalMergingException();
0212   }
0213
0214   if (aPortal.m_alongNormal != nullptr) {
0215     if (bPortal.m_alongNormal == nullptr) {
0216       ACTS_ERROR(
0217           "Portal A has link along normal, while b does not. This is not "
0218           "supported");
0219       throw PortalMergingException();
0220     }
0221
0222     ACTS_VERBOSE("Portals have links along normal, merging");
0223     mergedAlongNormal = PortalLinkBase::merge(std::move(aPortal.m_alongNormal),
0224                                               std::move(bPortal.m_alongNormal),
0225                                               direction, logger);
0226   }
0227
0228   if (aPortal.m_oppositeNormal != nullptr) {
0229     if (bPortal.m_oppositeNormal == nullptr) {
0230       ACTS_ERROR(
0231           "Portal A has link opposite normal, while b does not. This is not "
0232           "supported");
0233       throw PortalMergingException();
0234     }
0235
0236     ACTS_VERBOSE("Portals have links opposite normal, merging");
0237     mergedOppositeNormal = PortalLinkBase::merge(
0238         std::move(aPortal.m_oppositeNormal),
0239         std::move(bPortal.m_oppositeNormal), direction, logger);
0240   }
0241
0242   aPortal.m_surface.reset();
0243   bPortal.m_surface.reset();
0244   return Portal{gctx, std::move(mergedAlongNormal),
0245                 std::move(mergedOppositeNormal)};
0246 }
0247
0248 Portal Portal::fuse(const GeometryContext& gctx, Portal& aPortal,
0249                     Portal& bPortal, const Logger& logger) {
0250   ACTS_VERBOSE("Fusing two portals");
0251   if (&aPortal == &bPortal) {
0252     ACTS_ERROR("Cannot fuse a portal with itself");
0253     throw PortalMergingException{};
0254   }
0255
0256   bool aHasAlongNormal = aPortal.m_alongNormal != nullptr;
0257   bool aHasOppositeNormal = aPortal.m_oppositeNormal != nullptr;
0258   bool bHasAlongNormal = bPortal.m_alongNormal != nullptr;
0259   bool bHasOppositeNormal = bPortal.m_oppositeNormal != nullptr;
0260
0261   if (aPortal.m_surface == nullptr || bPortal.m_surface == nullptr) {
0262     ACTS_ERROR("Portals have no surface");
0263     throw PortalFusingException();
0264   }
0265
0266   if (aPortal.m_surface->associatedDetectorElement() != nullptr ||
0267       bPortal.m_surface->associatedDetectorElement() != nullptr) {
0268     ACTS_ERROR("Cannot fuse portals with detector elements");
0269     throw PortalFusingException();
0270   }
0271
0272   if (!isSameSurface(gctx, *aPortal.m_surface, *bPortal.m_surface)) {
0273     ACTS_ERROR("Portals have different surfaces");
0274     ACTS_ERROR("A: " << aPortal.m_surface->bounds());
0275     ACTS_ERROR("\n" << aPortal.m_surface->transform(gctx).matrix());
0276     ACTS_ERROR("B: " << bPortal.m_surface->bounds());
0277     ACTS_ERROR("\n" << bPortal.m_surface->transform(gctx).matrix());
0278     throw PortalFusingException();
0279   }
0280
0281   if (aPortal.m_surface->surfaceMaterial() != nullptr &&
0282       bPortal.m_surface->surfaceMaterial() != nullptr) {
0283     ACTS_ERROR("Cannot fuse portals if both have material");
0284     throw PortalFusingException();
0285   }
0286
0287   if (aHasAlongNormal == bHasAlongNormal ||
0288       aHasOppositeNormal == bHasOppositeNormal) {
0289     ACTS_ERROR("Portals have the same links attached");
0290     throw PortalFusingException();
0291   }
0292
0293   auto maybeConvertToGrid = [&](std::unique_ptr<PortalLinkBase> link)
0294       -> std::unique_ptr<PortalLinkBase> {
0295     auto* composite = dynamic_cast<CompositePortalLink*>(link.get());
0296     if (composite == nullptr) {
0297       return link;
0298     }
0299
0300     ACTS_VERBOSE("Converting composite to grid during portal fusing");
0301     return composite->makeGrid(gctx, logger);
0302   };
0303
0304   aPortal.m_surface.reset();
0305   bPortal.m_surface.reset();
0306   if (aHasAlongNormal) {
0307     ACTS_VERBOSE("Taking along normal from lhs, opposite normal from rhs");
0308     return Portal{gctx, maybeConvertToGrid(std::move(aPortal.m_alongNormal)),
0309                   maybeConvertToGrid(std::move(bPortal.m_oppositeNormal))};
0310   } else {
0311     ACTS_VERBOSE("Taking along normal from rhs, opposite normal from lhs");
0312     return Portal{gctx, maybeConvertToGrid(std::move(bPortal.m_alongNormal)),
0313                   maybeConvertToGrid(std::move(aPortal.m_oppositeNormal))};
0314   }
0315 }
0316
0317 bool Portal::isSameSurface(const GeometryContext& gctx, const Surface& a,
0318                            const Surface& b) {
0319   if (&a == &b) {
0320     return true;
0321   }
0322
0323   if (a.type() != b.type()) {
0324     return false;
0325   }
0326
0327   std::vector<double> aValues = a.bounds().values();
0328   std::vector<double> bValues = b.bounds().values();
0329   bool different = false;
0330   for (auto [aVal, bVal] : zip(aValues, bValues)) {
0331     if (std::abs(aVal - bVal) > s_onSurfaceTolerance) {
0332       different = true;
0333       break;
0334     }
0335   }
0336
0337   if (a.bounds().type() != b.bounds().type() || different) {
0338     return false;
0339   }
0340
0341   if (!a.transform(gctx).linear().isApprox(b.transform(gctx).linear(),
0342                                            s_transformEquivalentTolerance)) {
0343     return false;
0344   }
0345
0346   Vector3 delta =
0347       (a.transform(gctx).translation() - b.transform(gctx).translation())
0348           .cwiseAbs();
0349
0350   if (delta.maxCoeff() > s_onSurfaceTolerance) {
0351     return false;
0352   }
0353
0354   return true;
0355 };
0356
0357 void Portal::fill(TrackingVolume& volume) {
0358   if (m_alongNormal != nullptr && m_oppositeNormal != nullptr) {
0359     throw std::logic_error{"Portal is already filled"};
0360   }
0361
0362   if (m_surface == nullptr) {
0363     throw std::logic_error{"Portal has no existing link set, can't fill"};
0364   }
0365
0366   if (m_alongNormal == nullptr) {
0367     m_alongNormal = std::make_unique<TrivialPortalLink>(m_surface, volume);
0368   } else {
0369     assert(m_oppositeNormal == nullptr);
0370     m_oppositeNormal = std::make_unique<TrivialPortalLink>(m_surface, volume);
0371   }
0372 }
0373
0374 }  // namespace Acts