File indexing completed on 2025-01-18 09:11:17
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "Acts/Detector/detail/CylindricalDetectorHelper.hpp"
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/Definitions/Tolerance.hpp"
0013 #include "Acts/Detector/DetectorVolume.hpp"
0014 #include "Acts/Detector/Portal.hpp"
0015 #include "Acts/Detector/detail/DetectorVolumeConsistency.hpp"
0016 #include "Acts/Detector/detail/PortalHelper.hpp"
0017 #include "Acts/Geometry/CutoutCylinderVolumeBounds.hpp"
0018 #include "Acts/Surfaces/CylinderBounds.hpp"
0019 #include "Acts/Surfaces/CylinderSurface.hpp"
0020 #include "Acts/Surfaces/DiscSurface.hpp"
0021 #include "Acts/Surfaces/PlanarBounds.hpp"
0022 #include "Acts/Surfaces/PlaneSurface.hpp"
0023 #include "Acts/Surfaces/RadialBounds.hpp"
0024 #include "Acts/Surfaces/RectangleBounds.hpp"
0025 #include "Acts/Surfaces/Surface.hpp"
0026 #include "Acts/Surfaces/SurfaceBounds.hpp"
0027 #include "Acts/Utilities/Axis.hpp"
0028 #include "Acts/Utilities/BinningType.hpp"
0029 #include "Acts/Utilities/Enumerate.hpp"
0030 #include "Acts/Utilities/Grid.hpp"
0031 #include "Acts/Utilities/Helpers.hpp"
0032 #include "Acts/Utilities/StringHelpers.hpp"
0033
0034 #include <algorithm>
0035 #include <cmath>
0036 #include <cstddef>
0037 #include <iterator>
0038 #include <map>
0039 #include <numbers>
0040 #include <ostream>
0041 #include <stdexcept>
0042 #include <string>
0043 #include <tuple>
0044 #include <utility>
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072 namespace {
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082 Acts::Experimental::PortalReplacement createDiscReplacement(
0083 const Acts::Transform3& transform, const std::vector<double>& rBoundaries,
0084 const std::vector<double>& phiBoundaries, unsigned int index,
0085 Acts::Direction dir) {
0086
0087 Acts::AxisDirection stitchValue = phiBoundaries.size() == 2u
0088 ? Acts::AxisDirection::AxisR
0089 : Acts::AxisDirection::AxisPhi;
0090
0091 auto [minRit, maxRit] = std::ranges::minmax_element(rBoundaries);
0092 auto [sectorPhi, avgPhi] = Acts::range_medium(phiBoundaries);
0093
0094
0095 auto bounds = std::make_unique<Acts::RadialBounds>(*minRit, *maxRit,
0096 0.5 * sectorPhi, avgPhi);
0097
0098 auto surface = Acts::Surface::makeShared<Acts::DiscSurface>(
0099 transform, std::move(bounds));
0100
0101 const auto& stitchBoundaries =
0102 (stitchValue == Acts::AxisDirection::AxisR) ? rBoundaries : phiBoundaries;
0103 return Acts::Experimental::PortalReplacement(
0104 std::make_shared<Acts::Experimental::Portal>(surface), index, dir,
0105 stitchBoundaries, stitchValue);
0106 }
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118 Acts::Experimental::PortalReplacement createCylinderReplacement(
0119 const Acts::Transform3& transform, double r,
0120 const std::vector<double>& zBoundaries,
0121 const std::vector<double>& phiBoundaries, unsigned int index,
0122 Acts::Direction dir) {
0123
0124 Acts::AxisDirection stitchValue = phiBoundaries.size() == 2u
0125 ? Acts::AxisDirection::AxisZ
0126 : Acts::AxisDirection::AxisPhi;
0127 auto [lengthZ, medZ] = Acts::range_medium(zBoundaries);
0128 auto [sectorPhi, avgPhi] = Acts::range_medium(phiBoundaries);
0129
0130
0131 auto bounds = std::make_unique<Acts::CylinderBounds>(r, 0.5 * lengthZ,
0132 0.5 * sectorPhi, avgPhi);
0133
0134 auto surface = Acts::Surface::makeShared<Acts::CylinderSurface>(
0135 transform, std::move(bounds));
0136
0137
0138 const auto& stitchBoundaries =
0139 (stitchValue == Acts::AxisDirection::AxisZ) ? zBoundaries : phiBoundaries;
0140 return Acts::Experimental::PortalReplacement(
0141 std::make_shared<Acts::Experimental::Portal>(surface), index, dir,
0142 stitchBoundaries, stitchValue);
0143 }
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155 Acts::Experimental::PortalReplacement createSectorReplacement(
0156 const Acts::GeometryContext& gctx, const Acts::Vector3& volumeCenter,
0157 const Acts::Surface& refSurface, const std::vector<double>& boundaries,
0158 Acts::AxisDirection binning, unsigned int index, Acts::Direction dir) {
0159
0160 const auto& refTransform = refSurface.transform(gctx);
0161 auto refRotation = refTransform.rotation();
0162
0163 const auto& boundValues = refSurface.bounds().values();
0164 std::unique_ptr<Acts::PlanarBounds> bounds = nullptr;
0165
0166
0167 Acts::Transform3 transform = Acts::Transform3::Identity();
0168 if (binning == Acts::AxisDirection::AxisR) {
0169
0170 auto [range, medium] = Acts::range_medium(boundaries);
0171
0172
0173
0174 Acts::Vector3 pCenter = volumeCenter + medium * refRotation.col(1u);
0175 transform.pretranslate(pCenter);
0176
0177 double halfX =
0178 0.5 * (boundValues[Acts::RectangleBounds::BoundValues::eMaxX] -
0179 boundValues[Acts::RectangleBounds::BoundValues::eMinX]);
0180
0181 bounds = std::make_unique<Acts::RectangleBounds>(halfX, 0.5 * range);
0182 } else if (binning == Acts::AxisDirection::AxisZ) {
0183
0184 auto [range, medium] = Acts::range_medium(boundaries);
0185
0186 const auto& surfaceCenter = refSurface.center(gctx);
0187 Acts::Vector3 centerDiffs = (surfaceCenter - volumeCenter);
0188 double centerR = centerDiffs.dot(refRotation.col(2));
0189
0190 Acts::Vector3 pCenter = volumeCenter + centerR * refRotation.col(2);
0191 transform.pretranslate(pCenter);
0192
0193 double halfY =
0194 0.5 * (boundValues[Acts::RectangleBounds::BoundValues::eMaxY] -
0195 boundValues[Acts::RectangleBounds::BoundValues::eMinY]);
0196 bounds = std::make_unique<Acts::RectangleBounds>(0.5 * range, halfY);
0197 }
0198
0199 transform.prerotate(refRotation);
0200
0201 auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0202 transform, std::move(bounds));
0203
0204 Acts::Experimental::PortalReplacement pRep = {
0205 std::make_shared<Acts::Experimental::Portal>(surface), index, dir,
0206 boundaries, binning};
0207 return pRep;
0208 }
0209
0210
0211
0212
0213
0214
0215
0216 void checkVolumes(
0217 const Acts::GeometryContext& gctx,
0218 const std::vector<std::shared_ptr<Acts::Experimental::DetectorVolume>>&
0219 volumes) {
0220
0221
0222 std::string message = "CylindricalDetectorHelper: ";
0223 if (volumes.size() < 2u) {
0224 message += std::string("not enough volume given (") +
0225 std::to_string(volumes.size());
0226 message += std::string(" ), when required >=2.");
0227 throw std::invalid_argument(message.c_str());
0228 }
0229
0230 for (auto [iv, v] : Acts::enumerate(volumes)) {
0231
0232 if (v == nullptr) {
0233 message += "nullptr detector instead of volume " + std::to_string(iv);
0234 throw std::invalid_argument(message.c_str());
0235 }
0236
0237 if (v->volumeBounds().type() != Acts::VolumeBounds::BoundsType::eCylinder) {
0238 message += "non-cylindrical volume bounds detected for volume " +
0239 std::to_string(iv);
0240 throw std::invalid_argument(message.c_str());
0241 }
0242 }
0243
0244 Acts::Experimental::detail::DetectorVolumeConsistency::checkRotationAlignment(
0245 gctx, volumes);
0246 }
0247
0248
0249
0250
0251
0252
0253
0254
0255 void checkBounds(
0256 [[maybe_unused]] const Acts::GeometryContext& gctx,
0257 const std::vector<std::shared_ptr<Acts::Experimental::DetectorVolume>>&
0258 volumes,
0259 const std::vector<std::array<unsigned int, 2u>>& refCur) {
0260
0261 auto refValues = volumes[0u]->volumeBounds().values();
0262 for (auto [iv, v] : Acts::enumerate(volumes)) {
0263 if (iv > 0u) {
0264 auto curValues = v->volumeBounds().values();
0265 for (auto [r, c] : refCur) {
0266 if (std::abs(refValues[r] - curValues[c]) >
0267 Acts::s_onSurfaceTolerance) {
0268 std::string message = "CylindricalDetectorHelper: '";
0269 message += volumes[iv - 1]->name();
0270 if (r != c) {
0271 message += "' does not attach to '";
0272 } else {
0273 message += "' does not match with '";
0274 }
0275 message += volumes[iv]->name();
0276 message += "'\n";
0277 message += " - at bound values ";
0278 message += std::to_string(refValues[r]);
0279 message += " / " + std::to_string(curValues[c]);
0280 throw std::runtime_error(message.c_str());
0281 }
0282 }
0283 refValues = curValues;
0284 }
0285 }
0286 }
0287
0288 }
0289
0290 Acts::Experimental::DetectorComponent::PortalContainer
0291 Acts::Experimental::detail::CylindricalDetectorHelper::connectInR(
0292 const GeometryContext& gctx,
0293 std::vector<std::shared_ptr<Experimental::DetectorVolume>>& volumes,
0294 const std::vector<unsigned int>& selectedOnly,
0295 Acts::Logging::Level logLevel) {
0296
0297 checkVolumes(gctx, volumes);
0298
0299
0300
0301
0302 std::vector<std::array<unsigned int, 2u>> checks = {
0303 {1u, 0u}, {3u, 3u}, {4u, 4u}};
0304
0305 if (selectedOnly.empty()) {
0306 checks.push_back({2u, 2u});
0307 }
0308 checkBounds(gctx, volumes, checks);
0309
0310
0311 ACTS_LOCAL_LOGGER(getDefaultLogger("CylindricalDetectorHelper", logLevel));
0312
0313 ACTS_DEBUG("Connect " << volumes.size() << " detector volumes in R.");
0314
0315
0316 DetectorComponent::PortalContainer dShell;
0317
0318
0319 std::vector<double> rBoundaries = {};
0320 auto refValues = volumes[0u]->volumeBounds().values();
0321
0322
0323 rBoundaries.push_back(refValues[CylinderVolumeBounds::BoundValues::eMinR]);
0324 rBoundaries.push_back(refValues[CylinderVolumeBounds::BoundValues::eMaxR]);
0325
0326
0327 bool connectR = selectedOnly.empty() || rangeContainsValue(selectedOnly, 2u);
0328
0329
0330 double phiSector =
0331 refValues[CylinderVolumeBounds::BoundValues::eHalfPhiSector];
0332 double avgPhi = refValues[CylinderVolumeBounds::BoundValues::eAveragePhi];
0333
0334
0335 for (unsigned int iv = 1; iv < volumes.size(); ++iv) {
0336 refValues = volumes[iv]->volumeBounds().values();
0337
0338 rBoundaries.push_back(refValues[CylinderVolumeBounds::BoundValues::eMaxR]);
0339
0340 if (connectR) {
0341 ACTS_VERBOSE("Connect volume '" << volumes[iv - 1]->name() << "' to "
0342 << volumes[iv]->name() << "'.");
0343
0344
0345 auto innerCylinder = volumes[iv - 1]->portalPtrs()[2u];
0346 auto outerCylinder = volumes[iv]->portalPtrs()[3u];
0347 auto fusedCylinder = Portal::fuse(innerCylinder, outerCylinder);
0348 volumes[iv - 1]->updatePortal(fusedCylinder, 2u);
0349 volumes[iv]->updatePortal(fusedCylinder, 3u);
0350 }
0351 }
0352
0353
0354 if (connectR) {
0355
0356 if (volumes[0u]->portals().size() == 4u ||
0357 volumes[0u]->portals().size() == 6u) {
0358 dShell[3u] = volumes[0u]->portalPtrs()[3u];
0359 }
0360 dShell[2u] = volumes[volumes.size() - 1u]->portalPtrs()[2u];
0361 }
0362
0363
0364
0365
0366 bool sectorsPresent = volumes[volumes.size() - 1u]->portals().size() > 4u;
0367
0368
0369
0370 std::vector<PortalReplacement> pReplacements = {};
0371
0372
0373 std::vector<Acts::Direction> discDirs = {Acts::Direction::Forward(),
0374 Acts::Direction::Backward()};
0375 for (const auto [iu, idir] : enumerate(discDirs)) {
0376 if (selectedOnly.empty() || rangeContainsValue(selectedOnly, iu)) {
0377 const Surface& refSurface = volumes[0u]->portals()[iu]->surface();
0378 const Transform3& refTransform = refSurface.transform(gctx);
0379 pReplacements.push_back(createDiscReplacement(
0380 refTransform, rBoundaries, {avgPhi - phiSector, avgPhi + phiSector},
0381 iu, idir));
0382 }
0383 }
0384
0385
0386 if (sectorsPresent) {
0387 ACTS_VERBOSE("Sector planes are present, they need replacement.");
0388
0389 std::vector<Acts::Direction> sectorDirs = {Acts::Direction::Forward(),
0390 Acts::Direction::Backward()};
0391 Acts::Vector3 vCenter = volumes[0u]->transform(gctx).translation();
0392 for (const auto [iu, idir] : enumerate(sectorDirs)) {
0393
0394
0395 if (selectedOnly.empty() || rangeContainsValue(selectedOnly, iu + 4u)) {
0396
0397 const Surface& refSurface =
0398 volumes[volumes.size() - 1u]->portals()[iu + 4u]->surface();
0399 pReplacements.push_back(createSectorReplacement(
0400 gctx, vCenter, refSurface, rBoundaries, Acts::AxisDirection::AxisR,
0401 iu + 4ul, idir));
0402 }
0403 }
0404 } else {
0405 ACTS_VERBOSE(
0406 "No sector planes present, full 2 * std::numbers::pi cylindrical "
0407 "geometry.");
0408 }
0409
0410
0411 PortalHelper::attachExternalNavigationDelegates(gctx, volumes, pReplacements);
0412
0413
0414 ACTS_VERBOSE("Portals of " << volumes.size() << " volumes need updating.");
0415
0416 for (auto& iv : volumes) {
0417 ACTS_VERBOSE("- update portals of volume '" << iv->name() << "'.");
0418 for (auto& [p, i, dir, boundaries, binning] : pReplacements) {
0419
0420 dShell[i] = p;
0421
0422
0423
0424
0425
0426 std::size_t nPortals = iv->portals().size();
0427 bool innerPresent = (nPortals == 3u || nPortals == 5u);
0428 int iOffset = (innerPresent && i > 2u) ? -1 : 0;
0429 ACTS_VERBOSE("-- update portal with index "
0430 << i + iOffset << " (including offset " << iOffset << ")");
0431 iv->updatePortal(p, static_cast<unsigned int>(i + iOffset));
0432 }
0433 }
0434
0435 return dShell;
0436 }
0437
0438 Acts::Experimental::DetectorComponent::PortalContainer
0439 Acts::Experimental::detail::CylindricalDetectorHelper::connectInZ(
0440 const Acts::GeometryContext& gctx,
0441 std::vector<std::shared_ptr<Acts::Experimental::DetectorVolume>>& volumes,
0442 const std::vector<unsigned int>& selectedOnly,
0443 Acts::Logging::Level logLevel) {
0444
0445 checkVolumes(gctx, volumes);
0446
0447
0448 std::vector<std::array<unsigned int, 2u>> checks = {{3u, 3u}, {4u, 4u}};
0449
0450
0451 if (selectedOnly.empty()) {
0452 checks.push_back({0u, 0u});
0453 checks.push_back({1u, 1u});
0454 }
0455 checkBounds(gctx, volumes, checks);
0456
0457
0458 ACTS_LOCAL_LOGGER(getDefaultLogger("CylindricalDetectorHelper", logLevel));
0459
0460 ACTS_DEBUG("Connect " << volumes.size() << " detector volumes in Z.");
0461
0462
0463 DetectorComponent::PortalContainer dShell;
0464
0465
0466
0467
0468 const bool connectZ =
0469 selectedOnly.empty() || rangeContainsValue(selectedOnly, 1u);
0470
0471 const auto rotation = volumes[0u]->transform(gctx).rotation();
0472
0473 std::vector<Vector3> zBoundaries3D = {};
0474
0475
0476
0477
0478
0479 auto addZboundary3D = [&](const Experimental::DetectorVolume& volume,
0480 int side) -> void {
0481 const auto boundValues = volume.volumeBounds().values();
0482 double halflengthZ =
0483 boundValues[CylinderVolumeBounds::BoundValues::eHalfLengthZ];
0484 zBoundaries3D.push_back(volume.transform(gctx).translation() +
0485 side * halflengthZ * rotation.col(2));
0486 };
0487
0488
0489 addZboundary3D(*volumes[0u].get(), -1);
0490 addZboundary3D(*volumes[0u].get(), 1);
0491 for (unsigned int iv = 1; iv < volumes.size(); ++iv) {
0492
0493 addZboundary3D(*volumes[iv].get(), 1u);
0494
0495 if (connectZ) {
0496 ACTS_VERBOSE("Connect volume '" << volumes[iv - 1]->name() << "' to "
0497 << volumes[iv]->name() << "'.");
0498
0499 auto& pDisc = volumes[iv - 1]->portalPtrs()[1u];
0500 auto& nDisc = volumes[iv]->portalPtrs()[0u];
0501
0502 Vector3 pPosition = pDisc->surface().center(gctx);
0503 Vector3 nPosition = nDisc->surface().center(gctx);
0504 if (!pPosition.isApprox(nPosition)) {
0505 std::string message = "CylindricalDetectorHelper: '";
0506 message += volumes[iv - 1]->name();
0507 message += "' does not attach to '";
0508 message += volumes[iv]->name();
0509 message += "'\n";
0510 message += " - along z with values ";
0511 message += Acts::toString(pPosition);
0512 message += " / " + Acts::toString(nPosition);
0513 throw std::runtime_error(message.c_str());
0514 }
0515 auto fusedDisc = Portal::fuse(pDisc, nDisc);
0516 volumes[iv - 1]->updatePortal(fusedDisc, 1u);
0517 volumes[iv]->updatePortal(fusedDisc, 0u);
0518 }
0519 }
0520
0521
0522 if (connectZ) {
0523 dShell[0u] = volumes[0u]->portalPtrs()[0u];
0524 dShell[1u] = volumes[volumes.size() - 1u]->portalPtrs()[1u];
0525 }
0526
0527
0528 Vector3 combinedCenter =
0529 0.5 * (zBoundaries3D[zBoundaries3D.size() - 1u] + zBoundaries3D[0u]);
0530
0531 ACTS_VERBOSE("New combined center calculated at "
0532 << toString(combinedCenter));
0533
0534
0535 std::vector<double> zBoundaries = {};
0536 for (const auto& zb3D : zBoundaries3D) {
0537 auto proj3D = (zb3D - combinedCenter).dot(rotation.col(2));
0538 double zBoundary = std::copysign((zb3D - combinedCenter).norm(), proj3D);
0539 zBoundaries.push_back(zBoundary);
0540 }
0541
0542 Transform3 combinedTransform = Transform3::Identity();
0543 combinedTransform.pretranslate(combinedCenter);
0544 combinedTransform.rotate(rotation);
0545
0546
0547 const auto& refVolume = volumes[0u];
0548 const auto refValues = refVolume->volumeBounds().values();
0549
0550
0551 double minR = refValues[CylinderVolumeBounds::BoundValues::eMinR];
0552 double maxR = refValues[CylinderVolumeBounds::BoundValues::eMaxR];
0553 double phiSector =
0554 refValues[CylinderVolumeBounds::BoundValues::eHalfPhiSector];
0555 double avgPhi = refValues[CylinderVolumeBounds::BoundValues::eAveragePhi];
0556
0557
0558 std::size_t nPortals = volumes[volumes.size() - 1u]->portals().size();
0559 bool innerPresent = (nPortals != 3u && nPortals != 5u);
0560 bool sectorsPresent = nPortals > 4u;
0561
0562
0563
0564 std::vector<PortalReplacement> pReplacements = {};
0565
0566
0567 std::vector<Acts::Direction> cylinderDirs = {Acts::Direction::Backward()};
0568
0569 std::vector<double> cylinderR = {maxR};
0570 if (innerPresent) {
0571 ACTS_VERBOSE("Inner surface present, tube geometry detected.");
0572 cylinderDirs.push_back(Direction::Forward());
0573 cylinderR.push_back(minR);
0574 } else {
0575 ACTS_VERBOSE("No inner surface present, solid cylinder geometry detected.");
0576 }
0577
0578 unsigned int iSecOffset = innerPresent ? 4u : 3u;
0579
0580 for (const auto [iu, idir] : enumerate(cylinderDirs)) {
0581 if (selectedOnly.empty() || rangeContainsValue(selectedOnly, iu + 2u)) {
0582 pReplacements.push_back(createCylinderReplacement(
0583 combinedTransform, cylinderR[iu], zBoundaries,
0584 {avgPhi - phiSector, avgPhi + phiSector}, iu + 2u, idir));
0585 }
0586 }
0587
0588
0589 if (sectorsPresent) {
0590 ACTS_VERBOSE("Sector planes are present, they need replacement.");
0591
0592 std::vector<Acts::Direction> sectorDirs = {Acts::Direction::Forward(),
0593 Acts::Direction::Backward()};
0594 for (const auto [iu, idir] : enumerate(sectorDirs)) {
0595
0596 if (selectedOnly.empty() || rangeContainsValue(selectedOnly, iu + 4u)) {
0597 const Surface& refSurface =
0598 volumes[0u]->portals()[iu + iSecOffset]->surface();
0599 pReplacements.push_back(createSectorReplacement(
0600 gctx, combinedCenter, refSurface, zBoundaries,
0601 Acts::AxisDirection::AxisZ, iu + 4ul, idir));
0602 }
0603 }
0604 } else {
0605 ACTS_VERBOSE(
0606 "No sector planes present, full 2 * std::numbers::pi cylindrical "
0607 "geometry.");
0608 }
0609
0610
0611 PortalHelper::attachExternalNavigationDelegates(gctx, volumes, pReplacements);
0612
0613
0614 ACTS_VERBOSE("Portals of " << volumes.size() << " volumes need updating.");
0615 for (auto& iv : volumes) {
0616 ACTS_VERBOSE("- update portals of volume '" << iv->name() << "'.");
0617 for (auto& [p, i, dir, boundaries, binning] : pReplacements) {
0618
0619
0620
0621 int iOffset = (i > 2u && !innerPresent) ? -1 : 0;
0622 ACTS_VERBOSE("-- update portal with index " << i);
0623 iv->updatePortal(p, static_cast<unsigned int>(i + iOffset));
0624
0625 dShell[i] = p;
0626 }
0627 }
0628
0629 return dShell;
0630 }
0631
0632 Acts::Experimental::DetectorComponent::PortalContainer
0633 Acts::Experimental::detail::CylindricalDetectorHelper::connectInPhi(
0634 const Acts::GeometryContext& gctx,
0635 std::vector<std::shared_ptr<Acts::Experimental::DetectorVolume>>& volumes,
0636 const std::vector<unsigned int>& ,
0637 Acts::Logging::Level logLevel) {
0638
0639 checkVolumes(gctx, volumes);
0640
0641
0642 ACTS_LOCAL_LOGGER(getDefaultLogger("CylindricalDetectorHelper", logLevel));
0643
0644 ACTS_DEBUG("Connect " << volumes.size() << " detector volumes in phi.");
0645
0646
0647 DetectorComponent::PortalContainer dShell;
0648
0649
0650 std::size_t nPortals = volumes[volumes.size() - 1u]->portals().size();
0651 bool innerPresent = (nPortals != 3u && nPortals != 5u);
0652
0653 Transform3 refTransform = volumes[0u]->transform(gctx);
0654
0655
0656 unsigned int iSecOffset = innerPresent ? 4u : 3u;
0657 std::vector<double> phiBoundaries = {};
0658 auto refValues = volumes[0u]->volumeBounds().values();
0659 phiBoundaries.push_back(
0660 refValues[CylinderVolumeBounds::BoundValues::eAveragePhi] -
0661 refValues[CylinderVolumeBounds::BoundValues::eHalfPhiSector]);
0662 phiBoundaries.push_back(
0663 refValues[CylinderVolumeBounds::BoundValues::eAveragePhi] +
0664 refValues[CylinderVolumeBounds::BoundValues::eHalfPhiSector]);
0665
0666 for (unsigned int iv = 1; iv < volumes.size(); ++iv) {
0667 ACTS_VERBOSE("Connect volume '" << volumes[iv - 1]->name() << "' to "
0668 << volumes[iv]->name() << "'.");
0669
0670
0671 auto& rSector = volumes[iv - 1]->portalPtrs()[iSecOffset + 1u];
0672 auto& lSector = volumes[iv]->portalPtrs()[iSecOffset];
0673 auto fusedSector = Portal::fuse(rSector, lSector);
0674 volumes[iv - 1]->updatePortal(fusedSector, iSecOffset + 1u);
0675 volumes[iv]->updatePortal(fusedSector, iSecOffset);
0676
0677 auto curValues = volumes[iv]->volumeBounds().values();
0678
0679 double lowPhi =
0680 curValues[CylinderVolumeBounds::BoundValues::eAveragePhi] -
0681 curValues[CylinderVolumeBounds::BoundValues::eHalfPhiSector];
0682 double highPhi =
0683 curValues[CylinderVolumeBounds::BoundValues::eAveragePhi] +
0684 curValues[CylinderVolumeBounds::BoundValues::eHalfPhiSector];
0685
0686 if (std::abs(phiBoundaries[phiBoundaries.size() - 1u] - lowPhi) >
0687 Acts::s_onSurfaceTolerance) {
0688 std::string message = "CylindricalDetectorHelper: '";
0689 message += volumes[iv - 1]->name();
0690 message += "' does not attach to '";
0691 message += volumes[iv]->name();
0692 message += "'\n";
0693 message += " - within phi sectors ";
0694 message += std::to_string(lowPhi);
0695 message +=
0696 " / " + std::to_string(phiBoundaries[phiBoundaries.size() - 1u]);
0697 throw std::runtime_error(message.c_str());
0698 }
0699
0700 phiBoundaries.push_back(highPhi);
0701
0702 refValues = curValues;
0703 }
0704
0705
0706
0707 std::vector<PortalReplacement> pReplacements = {};
0708
0709 pReplacements.push_back(createDiscReplacement(
0710 refTransform,
0711 {refValues[CylinderVolumeBounds::BoundValues::eMinR],
0712 refValues[CylinderVolumeBounds::BoundValues::eMaxR]},
0713 phiBoundaries, 0u, Acts::Direction::Forward()));
0714
0715
0716 pReplacements.push_back(createDiscReplacement(
0717 refTransform,
0718 {refValues[CylinderVolumeBounds::BoundValues::eMinR],
0719 refValues[CylinderVolumeBounds::BoundValues::eMaxR]},
0720 phiBoundaries, 1u, Acts::Direction::Backward()));
0721
0722
0723 pReplacements.push_back(createCylinderReplacement(
0724 refTransform, refValues[CylinderVolumeBounds::BoundValues::eMaxR],
0725 {-refValues[CylinderVolumeBounds::BoundValues::eHalfLengthZ],
0726 refValues[CylinderVolumeBounds::BoundValues::eHalfLengthZ]},
0727 phiBoundaries, 2u, Acts::Direction::Backward()));
0728
0729
0730
0731 if (refValues[CylinderVolumeBounds::BoundValues::eMinR] > 0.) {
0732
0733 pReplacements.push_back(createCylinderReplacement(
0734 refTransform, refValues[CylinderVolumeBounds::BoundValues::eMinR],
0735 {-refValues[CylinderVolumeBounds::BoundValues::eHalfLengthZ],
0736 refValues[CylinderVolumeBounds::BoundValues::eHalfLengthZ]},
0737 phiBoundaries, 3u, Acts::Direction::Forward()));
0738 }
0739
0740
0741 PortalHelper::attachExternalNavigationDelegates(gctx, volumes, pReplacements);
0742
0743 ACTS_VERBOSE("Portals of " << volumes.size() << " volumes need updating.");
0744 for (auto& iv : volumes) {
0745 ACTS_VERBOSE("- update portals of volume '" << iv->name() << "'.");
0746 for (auto& [p, i, dir, boundaries, binning] : pReplacements) {
0747 ACTS_VERBOSE("-- update portal with index " << i);
0748 iv->updatePortal(p, static_cast<unsigned int>(i));
0749
0750 dShell[i] = p;
0751 }
0752 }
0753
0754
0755 return dShell;
0756 }
0757
0758 Acts::Experimental::DetectorComponent::PortalContainer
0759 Acts::Experimental::detail::CylindricalDetectorHelper::wrapInZR(
0760 const Acts::GeometryContext& gctx,
0761 std::vector<std::shared_ptr<Acts::Experimental::DetectorVolume>>& volumes,
0762 Acts::Logging::Level logLevel) {
0763
0764 ACTS_LOCAL_LOGGER(getDefaultLogger("CylindricalDetectorHelper", logLevel));
0765
0766 ACTS_DEBUG("Wrapping volumes in Z-R.");
0767
0768
0769 if (volumes.size() != 2u) {
0770 throw std::invalid_argument(
0771 "Wrapping the detector volume requires exactly 2 volumes.");
0772 }
0773
0774
0775 DetectorComponent::PortalContainer dShell;
0776
0777
0778 dShell[0u] = volumes[1u]->portalPtrs()[0u];
0779 dShell[1u] = volumes[1u]->portalPtrs()[1u];
0780 dShell[2u] = volumes[1u]->portalPtrs()[2u];
0781
0782
0783 auto& outerCover = volumes[0u]->portalPtrs()[2u];
0784 auto& innerCover = volumes[1u]->portalPtrs()[3u];
0785 auto fusedCover = Portal::fuse(outerCover, innerCover);
0786 volumes[0u]->updatePortal(fusedCover, 2u);
0787 volumes[1u]->updatePortal(fusedCover, 3u);
0788
0789
0790 auto& firstDiscN = volumes[1u]->portalPtrs()[4u];
0791 auto& secondDiscN = volumes[0u]->portalPtrs()[0u];
0792 auto fusedDiscN = Portal::fuse(firstDiscN, secondDiscN);
0793 volumes[1u]->updatePortal(fusedDiscN, 4u);
0794 volumes[0u]->updatePortal(fusedDiscN, 0u);
0795
0796
0797 auto& firstDiscP = volumes[0u]->portalPtrs()[1u];
0798 auto& secondDiscP = volumes[1u]->portalPtrs()[5u];
0799 auto fusedDiscP = Portal::fuse(firstDiscP, secondDiscP);
0800 volumes[0u]->updatePortal(fusedDiscP, 1u);
0801 volumes[1u]->updatePortal(fusedDiscP, 5u);
0802
0803
0804 if (volumes[0u]->portalPtrs().size() == 4u &&
0805 volumes[1u]->portalPtrs().size() == 8u) {
0806 const auto* cylVolBounds =
0807 dynamic_cast<const CylinderVolumeBounds*>(&volumes[0u]->volumeBounds());
0808 const auto* ccylVolBounds = dynamic_cast<const CutoutCylinderVolumeBounds*>(
0809 &volumes[1u]->volumeBounds());
0810 if (cylVolBounds == nullptr || ccylVolBounds == nullptr) {
0811 throw std::invalid_argument(
0812 "Wrapping the detector volume requires a cylinder and a cutout "
0813 "cylinder volume.");
0814 }
0815
0816 double hlZ = cylVolBounds->get(
0817 Acts::CylinderVolumeBounds::BoundValues::eHalfLengthZ);
0818 double HlZ = ccylVolBounds->get(
0819 Acts::CutoutCylinderVolumeBounds::BoundValues::eHalfLengthZ);
0820 double innerR = cylVolBounds->get(CylinderVolumeBounds::BoundValues::eMinR);
0821
0822 std::vector<PortalReplacement> pReplacements;
0823 pReplacements.push_back(createCylinderReplacement(
0824 volumes[0u]->transform(gctx), innerR, {-HlZ, -hlZ, hlZ, HlZ},
0825 {-std::numbers::pi, std::numbers::pi}, 3u, Direction::Forward()));
0826 std::vector<std::shared_ptr<DetectorVolume>> zVolumes = {
0827 volumes[1u], volumes[0u], volumes[1u]};
0828
0829 PortalHelper::attachExternalNavigationDelegates(gctx, zVolumes,
0830 pReplacements);
0831 auto& [p, i, dir, boundaries, binning] = pReplacements[0u];
0832
0833 volumes[1u]->updatePortal(p, 6u);
0834 volumes[0u]->updatePortal(p, 3u);
0835 volumes[1u]->updatePortal(p, 7u);
0836
0837 dShell[3u] = p;
0838 }
0839
0840 return dShell;
0841 }
0842
0843 Acts::Experimental::DetectorComponent::PortalContainer
0844 Acts::Experimental::detail::CylindricalDetectorHelper::connectInR(
0845 const GeometryContext& gctx,
0846 const std::vector<DetectorComponent::PortalContainer>& containers,
0847 const std::vector<unsigned int>& selectedOnly,
0848 Acts::Logging::Level logLevel) noexcept(false) {
0849
0850 ACTS_LOCAL_LOGGER(getDefaultLogger("CylindricalDetectorHelper", logLevel));
0851
0852 ACTS_DEBUG("Connect " << containers.size() << " proto containers in R.");
0853
0854
0855 DetectorComponent::PortalContainer dShell;
0856
0857
0858 for (unsigned int ic = 1; ic < containers.size(); ++ic) {
0859 auto& formerContainer = containers[ic - 1];
0860 auto& currentContainer = containers[ic];
0861
0862 if (!formerContainer.contains(2u)) {
0863 throw std::invalid_argument(
0864 "CylindricalDetectorHelper: proto container has no outer cover, can "
0865 "not be connected in R");
0866 }
0867 if (!currentContainer.contains(3u)) {
0868 throw std::invalid_argument(
0869 "CylindricalDetectorHelper: proto container has no inner cover, can "
0870 "not be connected in R");
0871 }
0872
0873
0874 std::shared_ptr<Portal> innerCylinder = containers[ic - 1].find(2u)->second;
0875
0876 auto innerAttachedVolumes =
0877 innerCylinder->attachedDetectorVolumes()[Direction::Backward().index()];
0878 std::shared_ptr<Portal> outerCylinder = containers[ic].find(3u)->second;
0879 auto outerAttachedVolume =
0880 outerCylinder->attachedDetectorVolumes()[Direction::Forward().index()];
0881 auto fusedCylinder = Portal::fuse(innerCylinder, outerCylinder);
0882
0883
0884 std::for_each(innerAttachedVolumes.begin(), innerAttachedVolumes.end(),
0885 [&](std::shared_ptr<DetectorVolume>& av) {
0886 av->updatePortal(fusedCylinder, 2u);
0887 });
0888 std::for_each(outerAttachedVolume.begin(), outerAttachedVolume.end(),
0889 [&](std::shared_ptr<DetectorVolume>& av) {
0890 av->updatePortal(fusedCylinder, 3u);
0891 });
0892 }
0893
0894
0895 if (containers[0u].contains(3u)) {
0896 dShell[3u] = containers[0u].find(3u)->second;
0897 }
0898 dShell[2u] = containers[containers.size() - 1u].find(2u)->second;
0899
0900 auto sideVolumes = PortalHelper::stripSideVolumes(
0901 containers, {0u, 1u, 4u, 5u}, selectedOnly, logLevel);
0902
0903 for (auto [s, volumes] : sideVolumes) {
0904 auto pR = connectInR(gctx, volumes, {s});
0905 if (pR.contains(s)) {
0906 dShell[s] = pR.find(s)->second;
0907 }
0908 }
0909
0910
0911 return dShell;
0912 }
0913
0914 Acts::Experimental::DetectorComponent::PortalContainer
0915 Acts::Experimental::detail::CylindricalDetectorHelper::connectInZ(
0916 const GeometryContext& gctx,
0917 const std::vector<DetectorComponent::PortalContainer>& containers,
0918 const std::vector<unsigned int>& selectedOnly,
0919 Acts::Logging::Level logLevel) noexcept(false) {
0920
0921 ACTS_LOCAL_LOGGER(getDefaultLogger("CylindricalDetectorHelper", logLevel));
0922
0923 ACTS_DEBUG("Connect " << containers.size() << " proto containers in Z.");
0924
0925
0926 DetectorComponent::PortalContainer dShell;
0927
0928 for (unsigned int ic = 1; ic < containers.size(); ++ic) {
0929 auto& formerContainer = containers[ic - 1];
0930 auto& currentContainer = containers[ic];
0931
0932 if (!formerContainer.contains(1u)) {
0933 throw std::invalid_argument(
0934 "CylindricalDetectorHelper: proto container has no negative disc, "
0935 "can not be connected in Z");
0936 }
0937 if (!currentContainer.contains(0u)) {
0938 throw std::invalid_argument(
0939 "CylindricalDetectorHelper: proto container has no positive disc, "
0940 "can not be connected in Z");
0941 }
0942
0943 std::shared_ptr<Portal> pDisc = formerContainer.find(1u)->second;
0944 auto pAttachedVolumes =
0945 pDisc->attachedDetectorVolumes()[Direction::Backward().index()];
0946
0947 std::shared_ptr<Portal> nDisc = currentContainer.find(0u)->second;
0948 auto nAttachedVolumes =
0949 nDisc->attachedDetectorVolumes()[Direction::Forward().index()];
0950
0951 auto fusedDisc = Portal::fuse(pDisc, nDisc);
0952
0953 std::for_each(pAttachedVolumes.begin(), pAttachedVolumes.end(),
0954 [&](std::shared_ptr<DetectorVolume>& av) {
0955 av->updatePortal(fusedDisc, 1u);
0956 });
0957 std::for_each(nAttachedVolumes.begin(), nAttachedVolumes.end(),
0958 [&](std::shared_ptr<DetectorVolume>& av) {
0959 av->updatePortal(fusedDisc, 0u);
0960 });
0961 }
0962
0963
0964 dShell[0u] = containers[0u].find(0u)->second;
0965 dShell[1u] = containers[containers.size() - 1u].find(1u)->second;
0966
0967
0968 std::vector<unsigned int> nominalSides = {2u, 4u, 5u};
0969 if (containers[0u].contains(3u)) {
0970 nominalSides.push_back(3u);
0971 }
0972
0973
0974 auto sideVolumes = PortalHelper::stripSideVolumes(containers, nominalSides,
0975 selectedOnly, logLevel);
0976
0977 ACTS_VERBOSE("There remain " << sideVolumes.size()
0978 << " side volume packs to be connected");
0979 for (auto [s, volumes] : sideVolumes) {
0980 ACTS_VERBOSE(" - connect " << volumes.size() << " at selected side " << s);
0981 auto pR = connectInZ(gctx, volumes, {s}, logLevel);
0982 if (pR.contains(s)) {
0983 dShell[s] = pR.find(s)->second;
0984 }
0985 }
0986
0987
0988 return dShell;
0989 }
0990
0991 Acts::Experimental::DetectorComponent::PortalContainer
0992 Acts::Experimental::detail::CylindricalDetectorHelper::connectInPhi(
0993 [[maybe_unused]] const GeometryContext& gctx,
0994 [[maybe_unused]] const std::vector<DetectorComponent::PortalContainer>&
0995 containers,
0996 [[maybe_unused]] const std::vector<unsigned int>& selectedOnly,
0997 [[maybe_unused]] Acts::Logging::Level logLevel) noexcept(false) {
0998 throw std::invalid_argument(
0999 "CylindricalDetectorHelper: container connection in phi not implemented "
1000 "yet.");
1001 DetectorComponent::PortalContainer dShell;
1002
1003 return dShell;
1004 }
1005
1006 Acts::Experimental::DetectorComponent::PortalContainer
1007 Acts::Experimental::detail::CylindricalDetectorHelper::wrapInZR(
1008 [[maybe_unused]] const GeometryContext& gctx,
1009 const std::vector<DetectorComponent::PortalContainer>& containers,
1010 Acts::Logging::Level logLevel) {
1011 if (containers.size() != 2u) {
1012 throw std::invalid_argument(
1013 "CylindricalDetectorHelper: wrapping must take exactly two "
1014 "containers.");
1015 }
1016
1017
1018 auto innerContainer = containers.front();
1019
1020 auto outerContainer = containers.back();
1021 std::shared_ptr<DetectorVolume> wrappingVolume = nullptr;
1022 for (const auto& [key, value] : outerContainer) {
1023 auto attachedVolumes = value->attachedDetectorVolumes();
1024 for (const auto& ava : attachedVolumes) {
1025 for (const auto& av : ava) {
1026 if (wrappingVolume == nullptr && av != nullptr) {
1027 wrappingVolume = av;
1028 } else if (wrappingVolume != nullptr && av != wrappingVolume) {
1029 throw std::invalid_argument(
1030 "CylindricalDetectorHelper: wrapping container must represent a "
1031 "single volume.");
1032 }
1033 }
1034 }
1035 }
1036 if (wrappingVolume == nullptr) {
1037 throw std::invalid_argument(
1038 "CylindricalDetectorHelper: wrapping volume could not be "
1039 "determined.");
1040 }
1041
1042
1043 ACTS_LOCAL_LOGGER(getDefaultLogger("CylindricalDetectorHelper", logLevel));
1044
1045 ACTS_DEBUG("Wrapping a container with volume `" << wrappingVolume->name()
1046 << "'.");
1047
1048 DetectorComponent::PortalContainer dShell;
1049
1050
1051 dShell[0u] = wrappingVolume->portalPtrs()[0u];
1052 dShell[1u] = wrappingVolume->portalPtrs()[1u];
1053 dShell[2u] = wrappingVolume->portalPtrs()[2u];
1054
1055
1056 auto& innerCover = innerContainer[2u];
1057 auto innerAttachedVolumes =
1058 innerCover->attachedDetectorVolumes()[Direction::Backward().index()];
1059 auto& innerTube = wrappingVolume->portalPtrs()[3u];
1060 auto fusedCover = Portal::fuse(innerCover, innerTube);
1061
1062 std::for_each(innerAttachedVolumes.begin(), innerAttachedVolumes.end(),
1063 [&](std::shared_ptr<DetectorVolume>& av) {
1064 av->updatePortal(fusedCover, 2u);
1065 });
1066 wrappingVolume->updatePortal(fusedCover, 3u);
1067
1068
1069
1070 auto& firstDiscN = innerContainer[0u];
1071
1072 auto firstNAttachedVolumes =
1073 firstDiscN->attachedDetectorVolumes()[Direction::Forward().index()];
1074
1075 auto& secondDiscN = wrappingVolume->portalPtrs()[4u];
1076 auto fusedDiscN = Portal::fuse(firstDiscN, secondDiscN);
1077
1078 std::for_each(firstNAttachedVolumes.begin(), firstNAttachedVolumes.end(),
1079 [&](std::shared_ptr<DetectorVolume>& av) {
1080 av->updatePortal(fusedDiscN, 0u);
1081 });
1082 wrappingVolume->updatePortal(fusedDiscN, 4u);
1083
1084
1085 auto& firstDiscP = innerContainer[1u];
1086 auto firstPAttachedVolumes =
1087 firstDiscP->attachedDetectorVolumes()[Direction::Backward().index()];
1088
1089 auto& secondDiscP = wrappingVolume->portalPtrs()[5u];
1090 auto fusedDiscP = Portal::fuse(firstDiscP, secondDiscP);
1091
1092 std::for_each(firstPAttachedVolumes.begin(), firstPAttachedVolumes.end(),
1093 [&](std::shared_ptr<DetectorVolume>& av) {
1094 av->updatePortal(fusedDiscP, 1u);
1095 });
1096
1097 wrappingVolume->updatePortal(fusedDiscP, 5u);
1098
1099
1100 if (innerContainer.size() == 4u &&
1101 wrappingVolume->portalPtrs().size() == 8u) {
1102
1103 auto& centralSegment = innerContainer[3u];
1104 auto centralValues = centralSegment->surface().bounds().values();
1105 double centralHalfLengthZ =
1106 centralValues[CylinderBounds::BoundValues::eHalfLengthZ];
1107
1108 auto& nSegment = wrappingVolume->portalPtrs()[6u];
1109 auto nValues = nSegment->surface().bounds().values();
1110 double nHalfLengthZ = nValues[CylinderBounds::BoundValues::eHalfLengthZ];
1111 auto& pSegment = wrappingVolume->portalPtrs()[7u];
1112 auto pValues = pSegment->surface().bounds().values();
1113 double pHalfLengthZ = pValues[CylinderBounds::BoundValues::eHalfLengthZ];
1114
1115 auto sideVolumes =
1116 PortalHelper::stripSideVolumes({innerContainer}, {3u}, {3u}, logLevel);
1117
1118
1119 std::vector<std::shared_ptr<DetectorVolume>> innerVolumes = {
1120 wrappingVolume->getSharedPtr()};
1121
1122 std::vector<double> zBoundaries = {-centralHalfLengthZ - 2 * nHalfLengthZ,
1123 centralHalfLengthZ};
1124
1125 for (auto& svs : sideVolumes) {
1126 for (auto& v : svs.second) {
1127 const auto* cylVolBounds =
1128 dynamic_cast<const CylinderVolumeBounds*>(&v->volumeBounds());
1129 if (cylVolBounds == nullptr) {
1130 throw std::invalid_argument(
1131 "CylindricalDetectorHelper: side volume must be a cylinder.");
1132 }
1133 double hlZ =
1134 cylVolBounds->get(CylinderVolumeBounds::BoundValues::eHalfLengthZ);
1135 zBoundaries.push_back(zBoundaries.back() + 2 * hlZ);
1136 innerVolumes.push_back(v);
1137 }
1138 }
1139
1140 zBoundaries.push_back(zBoundaries.back() + 2 * pHalfLengthZ);
1141 innerVolumes.push_back(wrappingVolume);
1142 }
1143
1144
1145 return dShell;
1146 }