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/Material/MergedMaterialMarker.hpp"
0019 #include "Acts/Surfaces/RegularSurface.hpp"
0020 #include "Acts/Utilities/Zip.hpp"
0021
0022 #include <algorithm>
0023 #include <cstdlib>
0024 #include <memory>
0025 #include <sstream>
0026 #include <stdexcept>
0027
0028 namespace Acts {
0029
0030 PortalMergingException::PortalMergingException(std::string message)
0031     : m_message{std::move(message)} {}
0032
0033 const char* PortalMergingException::what() const noexcept {
0034   return m_message.c_str();
0035 }
0036
0037 const char* PortalFusingException::what() const noexcept {
0038   return "Failure to fuse portals";
0039 }
0040
0041 Portal::Portal(Direction direction, std::unique_ptr<PortalLinkBase> link) {
0042   if (link == nullptr) {
0043     throw std::invalid_argument("Link must not be null");
0044   }
0045
0046   m_surface = link->surfacePtr();
0047
0048   if (direction == Direction::AlongNormal()) {
0049     m_alongNormal = std::move(link);
0050   } else {
0051     m_oppositeNormal = std::move(link);
0052   }
0053 }
0054
0055 Portal::Portal(Direction direction, std::shared_ptr<RegularSurface> surface,
0056                TrackingVolume& volume)
0057     : Portal(direction,
0058              std::make_unique<TrivialPortalLink>(std::move(surface), volume)) {}
0059
0060 Portal::Portal(const GeometryContext& gctx,
0061                std::unique_ptr<PortalLinkBase> alongNormal,
0062                std::unique_ptr<PortalLinkBase> oppositeNormal) {
0063   if (alongNormal == nullptr && oppositeNormal == nullptr) {
0064     throw std::invalid_argument("At least one link must be provided");
0065   }
0066
0067   if (alongNormal != nullptr) {
0068     setLink(gctx, Direction::AlongNormal(), std::move(alongNormal));
0069   }
0070   if (oppositeNormal != nullptr) {
0071     setLink(gctx, Direction::OppositeNormal(), std::move(oppositeNormal));
0072   }
0073 }
0074
0075 Portal::Portal(const GeometryContext& gctx, Arguments&& args) {
0076   if (!args.alongNormal.surface && !args.oppositeNormal.surface) {
0077     throw std::invalid_argument("At least one link must be provided");
0078   }
0079
0080   if (args.alongNormal.surface) {
0081     setLink(gctx, Direction::AlongNormal(),
0082             std::make_unique<TrivialPortalLink>(
0083                 std::move(args.alongNormal.surface), *args.alongNormal.volume));
0084   }
0085   if (args.oppositeNormal.surface) {
0086     setLink(gctx, Direction::OppositeNormal(),
0087             std::make_unique<TrivialPortalLink>(
0088                 std::move(args.oppositeNormal.surface),
0089                 *args.oppositeNormal.volume));
0090   }
0091 }
0092
0093 void Portal::setLink(const GeometryContext& gctx, Direction direction,
0094                      std::unique_ptr<PortalLinkBase> link) {
0095   if (link == nullptr) {
0096     throw std::invalid_argument("Link must not be null");
0097   }
0098
0099   auto& target =
0100       direction == Direction::AlongNormal() ? m_alongNormal : m_oppositeNormal;
0101   const auto& other =
0102       direction == Direction::AlongNormal() ? m_oppositeNormal : m_alongNormal;
0103
0104   // check if surfaces are identical
0105   if (m_surface != nullptr &&
0106       !isSameSurface(gctx, link->surface(), *m_surface)) {
0107     throw PortalFusingException();
0108   }
0109
0110   // check if they both have material but are not the same surface
0111   if (m_surface != nullptr && (m_surface.get() != &link->surface()) &&
0112       link->surface().hasMaterial() && m_surface->hasMaterial()) {
0113     throw PortalFusingException();
0114   }
0115
0116   target = std::move(link);
0117
0118   if (other == nullptr) {
0119     // We don't have an existing surface, take the one we just got
0120     m_surface = target->surfacePtr();
0121     return;
0122   }
0123
0124   if (target->surface().hasMaterial()) {
0125     // new link has material: assign that to existing link
0126     m_surface = target->surfacePtr();
0127     other->setSurface(m_surface);
0128   } else {
0129     // none have material, or the existing surface had material: assign the
0130     // existing surface by convention
0131     target->setSurface(m_surface);
0132   }
0133 }
0134
0135 void Portal::setLink(const GeometryContext& gctx, Direction direction,
0136                      std::shared_ptr<RegularSurface> surface,
0137                      TrackingVolume& volume) {
0138   setLink(gctx, direction,
0139           std::make_unique<TrivialPortalLink>(std::move(surface), volume));
0140 }
0141
0142 const PortalLinkBase* Portal::getLink(Direction direction) const {
0143   if (direction == Direction::AlongNormal()) {
0144     return m_alongNormal.get();
0145   } else {
0146     return m_oppositeNormal.get();
0147   }
0148 }
0149
0150 Result<const TrackingVolume*> Portal::resolveVolume(
0151     const GeometryContext& gctx, const Vector3& position,
0152     const Vector3& direction) const {
0153   assert(m_surface != nullptr);
0154   const Vector3 normal = m_surface->normal(gctx, position);
0155   Direction side = Direction::fromScalarZeroAsPositive(normal.dot(direction));
0156
0157   const PortalLinkBase* link = side == Direction::AlongNormal()
0158                                    ? m_alongNormal.get()
0159                                    : m_oppositeNormal.get();
0160
0161   if (link == nullptr) {
0162     // no link is attached in this direction => this is the end of the world as
0163     // we know it. (i feel fine)
0164     return nullptr;
0165   } else {
0166     auto res = link->resolveVolume(gctx, position);
0167     if (!res.ok()) {
0168       return res.error();
0169     }
0170     return *res;
0171   }
0172 }
0173
0174 bool Portal::isValid() const {
0175   return m_alongNormal != nullptr || m_oppositeNormal != nullptr;
0176 }
0177
0178 const RegularSurface& Portal::surface() const {
0179   assert(m_surface != nullptr);
0180   return *m_surface;
0181 }
0182
0183 RegularSurface& Portal::surface() {
0184   assert(m_surface != nullptr);
0185   return *m_surface;
0186 }
0187
0188 void Portal::addTag(std::string tag) {
0189   if (std::ranges::find(m_tags, tag) != m_tags.end()) {
0190     throw std::invalid_argument("Portal already has tag: " + tag);
0191   }
0192   m_tags.push_back(std::move(tag));
0193 }
0194
0195 std::span<const std::string> Portal::tags() const {
0196   return m_tags;
0197 }
0198
0199 // Note: tags are intentionally *not* propagated through merge/fuse. These
0200 // operations build new portals, and portal tagging is a post-stacking operation
0201 // (applied in the finalize phase of the blueprint construction, after all
0202 // merging and fusing has happened), so there is nothing to carry over here.
0203 Portal Portal::merge(const GeometryContext& gctx, Portal& aPortal,
0204                      Portal& bPortal, AxisDirection direction,
0205                      const Logger& logger,
0206                      PortalMaterialMergePolicy materialPolicy) {
0207   ACTS_VERBOSE("Merging two portals along " << direction);
0208
0209   if (&aPortal == &bPortal) {
0210     ACTS_ERROR("Cannot merge a portal with itself");
0211     throw PortalMergingException{};
0212   }
0213
0214   const bool aHasMaterial = aPortal.m_surface->hasMaterial();
0215   const bool bHasMaterial = bPortal.m_surface->hasMaterial();
0216   const bool hadMaterial = aHasMaterial || bHasMaterial;
0217
0218   if (hadMaterial) {
0219     std::stringstream ss;
0220     ss << "portals with material along " << direction << ": ";
0221     if (aHasMaterial) {
0222       ss << "portal A surface (bounds=" << aPortal.m_surface->bounds()
0223          << ", center=" << aPortal.m_surface->center(gctx).transpose()
0224          << ") carries material";
0225     }
0226     if (bHasMaterial) {
0227       ss << (aHasMaterial ? " and " : "")
0228          << "portal B surface (bounds=" << bPortal.m_surface->bounds()
0229          << ", center=" << bPortal.m_surface->center(gctx).transpose()
0230          << ") carries material";
0231     }
0232
0233     if (materialPolicy == PortalMaterialMergePolicy::eThrow) {
0234       std::stringstream es;
0235       es << "Cannot merge " << ss.str()
0236          << ". This typically means material was designated on a portal face "
0237             "that is subsequently merged during container stacking. Move the "
0238             "material designation to a face that is not merged (e.g. the face "
0239             "of the enclosing container).";
0240       ACTS_ERROR(es.str());
0241       throw PortalMergingException{es.str()};
0242     }
0243
0244     ACTS_WARNING("Merging "
0245                  << ss.str()
0246                  << ". The input material is discarded and the merged surface "
0247                     "is tagged with a MergedMaterialMarker. This is lossy: the "
0248                     "material designation should be moved to a face that is "
0249                     "not merged (e.g. the face of the enclosing container).");
0250   }
0251
0252   std::unique_ptr<PortalLinkBase> mergedAlongNormal = nullptr;
0253   std::unique_ptr<PortalLinkBase> mergedOppositeNormal = nullptr;
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 (aHasAlongNormal != bHasAlongNormal ||
0261       aHasOppositeNormal != bHasOppositeNormal) {
0262     ACTS_ERROR("Portals do not have the same links attached");
0263     throw PortalMergingException();
0264   }
0265
0266   if (aPortal.m_alongNormal != nullptr) {
0267     if (bPortal.m_alongNormal == nullptr) {
0268       ACTS_ERROR(
0269           "Portal A has link along normal, while b does not. This is not "
0270           "supported");
0271       throw PortalMergingException();
0272     }
0273
0274     ACTS_VERBOSE("Portals have links along normal, merging");
0275     mergedAlongNormal = PortalLinkBase::merge(std::move(aPortal.m_alongNormal),
0276                                               std::move(bPortal.m_alongNormal),
0277                                               direction, logger);
0278   }
0279
0280   if (aPortal.m_oppositeNormal != nullptr) {
0281     if (bPortal.m_oppositeNormal == nullptr) {
0282       ACTS_ERROR(
0283           "Portal A has link opposite normal, while b does not. This is not "
0284           "supported");
0285       throw PortalMergingException();
0286     }
0287
0288     ACTS_VERBOSE("Portals have links opposite normal, merging");
0289     mergedOppositeNormal = PortalLinkBase::merge(
0290         std::move(aPortal.m_oppositeNormal),
0291         std::move(bPortal.m_oppositeNormal), direction, logger);
0292   }
0293
0294   aPortal.m_surface.reset();
0295   bPortal.m_surface.reset();
0296   Portal merged{gctx, std::move(mergedAlongNormal),
0297                 std::move(mergedOppositeNormal)};
0298
0299   if (hadMaterial &&
0300       materialPolicy == PortalMaterialMergePolicy::eDiscardAndMark) {
0301     // Tag the merged surface so the lossy merge remains discoverable
0302     // downstream.
0303     merged.m_surface->assignSurfaceMaterial(
0304         std::make_shared<MergedMaterialMarker>());
0305   }
0306
0307   return merged;
0308 }
0309
0310 Portal Portal::fuse(const GeometryContext& gctx, Portal& aPortal,
0311                     Portal& bPortal, const Logger& logger) {
0312   ACTS_VERBOSE("Fusing two portals");
0313   if (&aPortal == &bPortal) {
0314     ACTS_ERROR("Cannot fuse a portal with itself");
0315     throw PortalMergingException{};
0316   }
0317
0318   bool aHasAlongNormal = aPortal.m_alongNormal != nullptr;
0319   bool aHasOppositeNormal = aPortal.m_oppositeNormal != nullptr;
0320   bool bHasAlongNormal = bPortal.m_alongNormal != nullptr;
0321   bool bHasOppositeNormal = bPortal.m_oppositeNormal != nullptr;
0322
0323   if (aPortal.m_surface == nullptr || bPortal.m_surface == nullptr) {
0324     ACTS_ERROR("Portals have no surface");
0325     throw PortalFusingException();
0326   }
0327
0328   if (aPortal.m_surface->isSensitive() || bPortal.m_surface->isSensitive()) {
0329     ACTS_ERROR("Cannot fuse portals with detector elements");
0330     throw PortalFusingException();
0331   }
0332
0333   if (!isSameSurface(gctx, *aPortal.m_surface, *bPortal.m_surface)) {
0334     ACTS_ERROR("Portals have different surfaces");
0335     ACTS_ERROR("A: " << aPortal.m_surface->bounds());
0336     ACTS_ERROR("\n"
0337                << aPortal.m_surface->localToGlobalTransform(gctx).matrix());
0338     ACTS_ERROR("B: " << bPortal.m_surface->bounds());
0339     ACTS_ERROR("\n"
0340                << bPortal.m_surface->localToGlobalTransform(gctx).matrix());
0341     throw PortalFusingException();
0342   }
0343
0344   if (aPortal.m_surface->hasMaterial() && bPortal.m_surface->hasMaterial()) {
0345     ACTS_ERROR("Cannot fuse portals if both have material");
0346     throw PortalFusingException();
0347   }
0348
0349   if (aHasAlongNormal == bHasAlongNormal ||
0350       aHasOppositeNormal == bHasOppositeNormal) {
0351     ACTS_ERROR("Portals have the same links attached");
0352     throw PortalFusingException();
0353   }
0354
0355   auto maybeConvertToGrid = [&](std::unique_ptr<PortalLinkBase> link)
0356       -> std::unique_ptr<PortalLinkBase> {
0357     auto* composite = dynamic_cast<CompositePortalLink*>(link.get());
0358     if (composite == nullptr) {
0359       return link;
0360     }
0361
0362     ACTS_VERBOSE("Converting composite to grid during portal fusing");
0363     return composite->makeGrid(gctx, logger);
0364   };
0365
0366   aPortal.m_surface.reset();
0367   bPortal.m_surface.reset();
0368   if (aHasAlongNormal) {
0369     ACTS_VERBOSE("Taking along normal from lhs, opposite normal from rhs");
0370     return Portal{gctx, maybeConvertToGrid(std::move(aPortal.m_alongNormal)),
0371                   maybeConvertToGrid(std::move(bPortal.m_oppositeNormal))};
0372   } else {
0373     ACTS_VERBOSE("Taking along normal from rhs, opposite normal from lhs");
0374     return Portal{gctx, maybeConvertToGrid(std::move(bPortal.m_alongNormal)),
0375                   maybeConvertToGrid(std::move(aPortal.m_oppositeNormal))};
0376   }
0377 }
0378
0379 bool Portal::isSameSurface(const GeometryContext& gctx, const Surface& a,
0380                            const Surface& b) {
0381   if (&a == &b) {
0382     return true;
0383   }
0384
0385   if (a.type() != b.type()) {
0386     return false;
0387   }
0388
0389   std::vector<double> aValues = a.bounds().values();
0390   std::vector<double> bValues = b.bounds().values();
0391   bool different = false;
0392   for (auto [aVal, bVal] : zip(aValues, bValues)) {
0393     if (std::abs(aVal - bVal) > s_onSurfaceTolerance) {
0394       different = true;
0395       break;
0396     }
0397   }
0398
0399   if (a.bounds().type() != b.bounds().type() || different) {
0400     return false;
0401   }
0402
0403   if (!a.localToGlobalTransform(gctx).linear().isApprox(
0404           b.localToGlobalTransform(gctx).linear(),
0405           s_transformEquivalentTolerance)) {
0406     return false;
0407   }
0408
0409   Vector3 delta = (a.localToGlobalTransform(gctx).translation() -
0410                    b.localToGlobalTransform(gctx).translation())
0411                       .cwiseAbs();
0412
0413   if (delta.maxCoeff() > s_onSurfaceTolerance) {
0414     return false;
0415   }
0416
0417   return true;
0418 };
0419
0420 void Portal::fill(TrackingVolume& volume) {
0421   if (m_alongNormal != nullptr && m_oppositeNormal != nullptr) {
0422     throw std::logic_error{"Portal is already filled"};
0423   }
0424
0425   if (m_surface == nullptr) {
0426     throw std::logic_error{"Portal has no existing link set, can't fill"};
0427   }
0428
0429   if (m_alongNormal == nullptr) {
0430     m_alongNormal = std::make_unique<TrivialPortalLink>(m_surface, volume);
0431   } else {
0432     assert(m_oppositeNormal == nullptr);
0433     m_oppositeNormal = std::make_unique<TrivialPortalLink>(m_surface, volume);
0434   }
0435 }
0436
0437 }  // namespace Acts