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