File indexing completed on 2026-09-20 08:20:14
0001
0002
0003
0004
0005
0006
0007
0008
0009 #include "Acts/Geometry/CuboidVolumeStack.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Tolerance.hpp"
0013 #include "Acts/Geometry/CuboidVolumeBounds.hpp"
0014 #include "Acts/Utilities/AxisDefinitions.hpp"
0015 #include "Acts/Utilities/Logger.hpp"
0016 #include "Acts/Utilities/StringHelpers.hpp"
0017
0018 #include <algorithm>
0019 #include <cstddef>
0020 #include <initializer_list>
0021 #include <iomanip>
0022 #include <memory>
0023 #include <numeric>
0024 #include <sstream>
0025 #include <stdexcept>
0026 #include <utility>
0027
0028 namespace Acts {
0029
0030 struct CuboidVolumeStack::VolumeTuple {
0031 Volume* volume{};
0032 const CuboidVolumeBounds* bounds{};
0033 std::shared_ptr<CuboidVolumeBounds> updatedBounds{};
0034 Transform3 localTransform = Transform3::Identity();
0035 Transform3 globalTransform = Transform3::Identity();
0036
0037 bool transformDirty = false;
0038
0039 explicit VolumeTuple(const GeometryContext& gctx, Volume& volume_,
0040 const Transform3& groupTransform)
0041 : volume{&volume_},
0042 localTransform{groupTransform.inverse() *
0043 volume_.localToGlobalTransform(gctx)},
0044 globalTransform{volume_.localToGlobalTransform(gctx)} {
0045 bounds = dynamic_cast<const CuboidVolumeBounds*>(&volume_.volumeBounds());
0046 assert(bounds != nullptr);
0047 updatedBounds = std::make_shared<CuboidVolumeBounds>(*bounds);
0048 }
0049
0050 double mid(AxisDirection direction) const {
0051 return localTransform.translation()[axisToIndex(direction)];
0052 }
0053 double halfLength(AxisDirection direction) const {
0054 return updatedBounds->get(
0055 CuboidVolumeBounds::boundsFromAxisDirection(direction));
0056 }
0057 double min(AxisDirection direction) const {
0058 return mid(direction) - halfLength(direction);
0059 }
0060 double max(AxisDirection direction) const {
0061 return mid(direction) + halfLength(direction);
0062 }
0063
0064 void set(
0065 std::initializer_list<std::pair<CuboidVolumeBounds::BoundValues, double>>
0066 keyValues) {
0067 updatedBounds->set(keyValues);
0068 }
0069
0070 void setLocalTransform(const Transform3& transform,
0071 const Transform3& groupTransform) {
0072 localTransform = transform;
0073 globalTransform = groupTransform * localTransform;
0074 transformDirty = true;
0075 }
0076
0077 void commit(const GeometryContext& gctx, const Logger& logger) {
0078
0079 auto copy = std::make_shared<CuboidVolumeBounds>(*updatedBounds);
0080
0081 std::optional<Transform3> transform = std::nullopt;
0082 if (transformDirty) {
0083 transform = globalTransform;
0084 }
0085
0086 volume->update(gctx, std::move(updatedBounds), transform, logger);
0087 bounds = copy.get();
0088 updatedBounds = std::move(copy);
0089 transformDirty = false;
0090 }
0091 };
0092
0093 std::size_t CuboidVolumeStack::axisToIndex(AxisDirection direction) {
0094 switch (direction) {
0095 case AxisDirection::AxisX:
0096 return 0;
0097 break;
0098 case AxisDirection::AxisY:
0099 return 1;
0100 break;
0101 case AxisDirection::AxisZ:
0102 return 2;
0103 break;
0104 default:
0105 throw std::invalid_argument("Invalid axis direction");
0106 }
0107 }
0108
0109 std::pair<AxisDirection, AxisDirection> CuboidVolumeStack::getOrthogonalAxes(
0110 AxisDirection direction) {
0111 switch (direction) {
0112 case AxisDirection::AxisX:
0113 return {AxisDirection::AxisY, AxisDirection::AxisZ};
0114 break;
0115 case AxisDirection::AxisY:
0116 return {AxisDirection::AxisZ, AxisDirection::AxisX};
0117 break;
0118 case AxisDirection::AxisZ:
0119 return {AxisDirection::AxisX, AxisDirection::AxisY};
0120 break;
0121 default:
0122 throw std::invalid_argument("Invalid axis direction");
0123 }
0124 }
0125
0126 CuboidVolumeStack::CuboidVolumeStack(const GeometryContext& gctx,
0127 std::vector<Volume*>& volumes,
0128 AxisDirection direction,
0129 VolumeAttachmentStrategy strategy,
0130 VolumeResizeStrategy resizeStrategy,
0131 const Logger& logger)
0132 : VolumeStack(volumes, direction, {resizeStrategy, resizeStrategy}) {
0133 std::tie(m_dirOrth1, m_dirOrth2) = getOrthogonalAxes(m_direction);
0134
0135 initializeOuterVolume(gctx, strategy, logger);
0136 }
0137
0138 void CuboidVolumeStack::initializeOuterVolume(const GeometryContext& gctx,
0139 VolumeAttachmentStrategy strategy,
0140 const Logger& logger) {
0141 ACTS_DEBUG("Creating CuboidVolumeStack from "
0142 << m_volumes.size() << " volumes in direction "
0143 << axisDirectionName(m_direction));
0144 if (m_volumes.empty()) {
0145 throw std::invalid_argument(
0146 "CuboidVolumeStack requires at least one volume");
0147 }
0148
0149 if (m_direction != Acts::AxisDirection::AxisX &&
0150 m_direction != Acts::AxisDirection::AxisY &&
0151 m_direction != Acts::AxisDirection::AxisZ) {
0152 throw std::invalid_argument(axisDirectionName(m_direction) +
0153 " is not supported ");
0154 }
0155
0156
0157 m_groupTransform = m_volumes.front()->localToGlobalTransform(gctx);
0158 ACTS_VERBOSE("Initial group transform is:\n" << m_groupTransform.matrix());
0159
0160 std::vector<VolumeTuple> volumeTuples;
0161 volumeTuples.reserve(m_volumes.size());
0162
0163 for (const auto& volume : m_volumes) {
0164 const auto* cuboidBounds =
0165 dynamic_cast<const CuboidVolumeBounds*>(&volume->volumeBounds());
0166 if (cuboidBounds == nullptr) {
0167 throw std::invalid_argument{
0168 "CuboidVolumeStack requires all volumes to "
0169 "have CuboidVolumeBounds"};
0170 }
0171
0172 volumeTuples.emplace_back(gctx, *volume, m_groupTransform);
0173 }
0174
0175 ACTS_DEBUG("*** Initial volume configuration:");
0176 printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG);
0177
0178 if (m_volumes.size() == 1) {
0179 ACTS_VERBOSE("Only one volume, returning");
0180 setTransform(m_volumes.front()->localToGlobalTransform(gctx));
0181 const auto* bounds = dynamic_cast<const CuboidVolumeBounds*>(
0182 &m_volumes.front()->volumeBounds());
0183 assert(bounds != nullptr && "Volume bounds are not cuboid bounds");
0184 Volume::update(gctx, std::make_shared<CuboidVolumeBounds>(*bounds),
0185 std::nullopt, logger);
0186 ACTS_VERBOSE(
0187 "Transform is now: " << toString(localToGlobalTransform(gctx)));
0188 return;
0189 }
0190
0191 ACTS_VERBOSE("Checking volume alignment");
0192 checkVolumeAlignment(volumeTuples, logger);
0193
0194 auto dirIdx = axisToIndex(m_direction);
0195 ACTS_VERBOSE("Sorting by volume " << axisDirectionName(m_direction)
0196 << " position");
0197 std::ranges::sort(volumeTuples, {}, [dirIdx](const auto& v) {
0198 return v.localTransform.translation()[dirIdx];
0199 });
0200 ACTS_VERBOSE("Checking for overlaps and attaching volumes in "
0201 << axisDirectionName(m_direction));
0202 std::vector<VolumeTuple> gapVolumes =
0203 checkOverlapAndAttach(gctx, volumeTuples, strategy, logger);
0204
0205 ACTS_VERBOSE("Appending " << gapVolumes.size()
0206 << " gap volumes to the end of the volume vector");
0207 std::ranges::copy(gapVolumes, std::back_inserter(volumeTuples));
0208
0209 ACTS_VERBOSE("*** Volume configuration after "
0210 << axisDirectionName(m_direction) << " attachment:");
0211 printVolumeSequence(volumeTuples, logger, Acts::Logging::VERBOSE);
0212
0213 ACTS_VERBOSE("Synchronizing bounds in " << axisDirectionName(m_dirOrth1)
0214 << "/"
0215 << axisDirectionName(m_dirOrth2));
0216 const auto [hl1, hl2] = synchronizeBounds(volumeTuples, logger);
0217
0218 for (auto& vt : volumeTuples) {
0219 ACTS_VERBOSE("Updated bounds for volume at "
0220 << axisDirectionName(m_direction) << ": "
0221 << vt.localTransform.translation()[dirIdx]);
0222 ACTS_VERBOSE(*vt.updatedBounds);
0223
0224 vt.commit(gctx, logger);
0225 }
0226
0227 ACTS_VERBOSE("*** Volume configuration after "
0228 << axisDirectionName(m_dirOrth1) << "/"
0229 << axisDirectionName(m_dirOrth2) << " synchronization:");
0230 printVolumeSequence(volumeTuples, logger, Acts::Logging::VERBOSE);
0231
0232 std::ranges::sort(volumeTuples, {},
0233 [*this](const auto& v) { return v.mid(m_direction); });
0234
0235 m_volumes.clear();
0236 for (const auto& vt : volumeTuples) {
0237 m_volumes.push_back(vt.volume);
0238 }
0239
0240 ACTS_DEBUG("*** Volume configuration after final "
0241 << axisDirectionName(m_direction) << " sorting:");
0242 printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG);
0243
0244 double min = volumeTuples.front().min(m_direction);
0245 double max = volumeTuples.back().max(m_direction);
0246
0247 double mid = std::midpoint(min, max);
0248 double hl = std::midpoint(max, -min);
0249
0250 Translation3 translation(Vector3::Unit(dirIdx) * mid);
0251 auto bounds = std::make_shared<CuboidVolumeBounds>(
0252 std::initializer_list<std::pair<CuboidVolumeBounds::BoundValues, double>>{
0253 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction), hl},
0254 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1), hl1},
0255 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2), hl2}});
0256 Volume::update(gctx, bounds, m_groupTransform * translation, logger);
0257 ACTS_DEBUG("Outer bounds are:\n" << volumeBounds());
0258 ACTS_DEBUG("Outer transform / new group transform is:\n"
0259 << toString(localToGlobalTransform(gctx)));
0260
0261
0262
0263 m_groupTransform = localToGlobalTransform(gctx);
0264 }
0265
0266 void CuboidVolumeStack::overlapPrint(const CuboidVolumeStack::VolumeTuple& a,
0267 const CuboidVolumeStack::VolumeTuple& b,
0268 const Logger& logger) {
0269 if (logger().doPrint(Acts::Logging::DEBUG)) {
0270 std::stringstream ss;
0271 ss << std::fixed;
0272 ss << std::setprecision(3);
0273 ss << std::setfill(' ');
0274
0275 int w = 9;
0276
0277 ACTS_VERBOSE("Checking overlap between");
0278 ss << " - " << " " << axisDirectionName(m_direction) << ": [ "
0279 << std::setw(w) << a.min(m_direction) << " <- " << std::setw(w)
0280 << a.mid(m_direction) << " -> " << std::setw(w) << a.max(m_direction)
0281 << " ]";
0282 ACTS_VERBOSE(ss.str());
0283
0284 ss.str("");
0285 ss << " - " << " " << axisDirectionName(m_direction) << ": [ "
0286 << std::setw(w) << b.min(m_direction) << " <- " << std::setw(w)
0287 << b.mid(m_direction) << " -> " << std::setw(w) << b.max(m_direction)
0288 << " ]";
0289 ACTS_VERBOSE(ss.str());
0290 }
0291 }
0292
0293 std::vector<CuboidVolumeStack::VolumeTuple>
0294 CuboidVolumeStack::checkOverlapAndAttach(const GeometryContext& gctx,
0295 std::vector<VolumeTuple>& volumes,
0296 VolumeAttachmentStrategy strategy,
0297 const Logger& logger) {
0298
0299 auto dirIdx = axisToIndex(m_direction);
0300 auto dirBoundIdx = CuboidVolumeBounds::boundsFromAxisDirection(m_direction);
0301
0302 std::vector<VolumeTuple> gapVolumes;
0303 for (std::size_t i = 0; i < volumes.size() - 1; i++) {
0304 std::size_t j = i + 1;
0305 auto& a = volumes.at(i);
0306 auto& b = volumes.at(j);
0307
0308 overlapPrint(a, b, logger);
0309
0310
0311 constexpr auto tolerance = s_onSurfaceTolerance;
0312 if (a.max(m_direction) - tolerance > b.min(m_direction)) {
0313 ACTS_ERROR(" -> Overlap in " << axisDirectionName(m_direction));
0314 throw std::invalid_argument("Volumes overlap in " +
0315 axisDirectionName(m_direction));
0316 } else {
0317 ACTS_VERBOSE(" -> No overlap");
0318 }
0319
0320 if (std::abs(a.max(m_direction) - b.min(m_direction)) < tolerance) {
0321 ACTS_VERBOSE("No gap between volumes, no attachment needed");
0322 } else {
0323 double gapWidth = b.min(m_direction) - a.max(m_direction);
0324 ACTS_VERBOSE("Gap width: " << gapWidth);
0325
0326 ACTS_VERBOSE("Synchronizing bounds in "
0327 << axisDirectionName(m_direction)
0328 << " with strategy: " << strategy);
0329 switch (strategy) {
0330 case VolumeAttachmentStrategy::Midpoint: {
0331 ACTS_VERBOSE(" -> Strategy: Expand both volumes to midpoint");
0332
0333 double aMidNew =
0334 (a.min(m_direction) + a.max(m_direction)) / 2.0 + gapWidth / 4.0;
0335 double aHlNew = a.halfLength(m_direction) + gapWidth / 4.0;
0336 ACTS_VERBOSE(" - New halflength for first volume: " << aHlNew);
0337 ACTS_VERBOSE(" - New bounds for first volume: ["
0338 << (aMidNew - aHlNew) << " <- " << aMidNew << " -> "
0339 << (aMidNew + aHlNew) << "]");
0340
0341 assert(std::abs(a.min(m_direction) - (aMidNew - aHlNew)) < 1e-9 &&
0342 "Volume shrunk");
0343 assert(aHlNew >= a.halfLength(m_direction) && "Volume shrunk");
0344
0345 double bMidNew =
0346 (b.min(m_direction) + b.max(m_direction)) / 2.0 - gapWidth / 4.0;
0347 double bHlNew = b.halfLength(m_direction) + gapWidth / 4.0;
0348 ACTS_VERBOSE(" - New halflength for second volume: " << bHlNew);
0349 ACTS_VERBOSE(" - New bounds for second volume: ["
0350 << (bMidNew - bHlNew) << " <- " << bMidNew << " -> "
0351 << (bMidNew + bHlNew) << "]");
0352
0353 assert(bHlNew >= b.halfLength(m_direction) && "Volume shrunk");
0354 assert(std::abs(b.max(m_direction) - (bMidNew + bHlNew)) < 1e-9 &&
0355 "Volume shrunk");
0356
0357 Translation3 translationA(Vector3::Unit(dirIdx) * aMidNew);
0358 a.setLocalTransform(Transform3{translationA}, m_groupTransform);
0359 a.updatedBounds->set(dirBoundIdx, aHlNew);
0360
0361 Translation3 translationB(Vector3::Unit(dirIdx) * bMidNew);
0362 b.setLocalTransform(Transform3{translationB}, m_groupTransform);
0363 b.updatedBounds->set(dirBoundIdx, bHlNew);
0364
0365 break;
0366 }
0367 case VolumeAttachmentStrategy::First: {
0368 ACTS_VERBOSE(" -> Strategy: Expand first volume");
0369 double aMidNew = (a.min(m_direction) + b.min(m_direction)) / 2.0;
0370 double aHlNew = (b.min(m_direction) - a.min(m_direction)) / 2.0;
0371 ACTS_VERBOSE(" - Gap width: " << gapWidth);
0372 ACTS_VERBOSE(" - New bounds for first volume: ["
0373 << (aMidNew - aHlNew) << " <- " << aMidNew << " -> "
0374 << (aMidNew + aHlNew) << "]");
0375
0376 assert(std::abs(a.min(m_direction) - (aMidNew - aHlNew)) < 1e-9 &&
0377 "Volume shrunk");
0378 assert(aHlNew >= a.halfLength(m_direction) && "Volume shrunk");
0379
0380 Translation3 translationA(Vector3::Unit(dirIdx) * aMidNew);
0381 a.setLocalTransform(Transform3{translationA}, m_groupTransform);
0382 a.updatedBounds->set(dirBoundIdx, aHlNew);
0383
0384 break;
0385 }
0386 case VolumeAttachmentStrategy::Second: {
0387 ACTS_VERBOSE(" -> Strategy: Expand second volume");
0388 double bMidNew = (a.max(m_direction) + b.max(m_direction)) / 2.0;
0389 double bHlNew = (b.max(m_direction) - a.max(m_direction)) / 2.0;
0390 ACTS_VERBOSE(" - New halflength for second volume: " << bHlNew);
0391 ACTS_VERBOSE(" - New bounds for second volume: ["
0392 << (bMidNew - bHlNew) << " <- " << bMidNew << " -> "
0393 << (bMidNew + bHlNew) << "]");
0394
0395 assert(bHlNew >= b.halfLength(m_direction) && "Volume shrunk");
0396 assert(std::abs(b.max(m_direction) - (bMidNew + bHlNew)) < 1e-9 &&
0397 "Volume shrunk");
0398
0399 Translation3 translationB(Vector3::Unit(dirIdx) * bMidNew);
0400 b.setLocalTransform(Transform3{translationB}, m_groupTransform);
0401 b.updatedBounds->set(dirBoundIdx, bHlNew);
0402 break;
0403 }
0404 case VolumeAttachmentStrategy::Gap: {
0405 ACTS_VERBOSE(" -> Strategy: Create a gap volume");
0406 double gapHl = (b.min(m_direction) - a.max(m_direction)) / 2.0;
0407 double gapMid = (b.min(m_direction) + a.max(m_direction)) / 2.0;
0408
0409 ACTS_VERBOSE(" - Gap half length: " << gapHl << " at "
0410 << axisDirectionName(m_direction)
0411 << ": " << gapMid);
0412
0413 Translation3 gapTranslation(Vector3::Unit(dirIdx) * gapMid);
0414
0415 double min1 = std::min(a.min(m_dirOrth1), b.min(m_dirOrth1));
0416 double max1 = std::max(a.max(m_dirOrth1), b.max(m_dirOrth1));
0417
0418 double min2 = std::min(a.min(m_dirOrth2), b.min(m_dirOrth2));
0419 double max2 = std::max(a.max(m_dirOrth2), b.max(m_dirOrth2));
0420
0421 Transform3 gapLocalTransform{gapTranslation};
0422 Transform3 gapGlobalTransform = m_groupTransform * gapLocalTransform;
0423
0424 auto gapBounds = std::make_shared<CuboidVolumeBounds>(
0425 std::initializer_list<
0426 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0427 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0428 gapHl},
0429 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0430 (max1 - min1) / 2},
0431 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0432 (max2 - min2) / 2}});
0433 auto gap = addGapVolume(gapGlobalTransform, gapBounds);
0434 gapVolumes.emplace_back(gctx, *gap, m_groupTransform);
0435
0436 break;
0437 }
0438 default:
0439 ACTS_ERROR("Attachment strategy " << strategy << " not implemented");
0440 std::stringstream ss;
0441 ss << strategy;
0442 throw std::invalid_argument("Attachment strategy " + ss.str() +
0443 " not implemented");
0444 }
0445 }
0446 }
0447
0448 return gapVolumes;
0449 }
0450
0451 void CuboidVolumeStack::printVolumeSequence(
0452 const std::vector<VolumeTuple>& volumes, const Logger& logger,
0453 Acts::Logging::Level lvl) {
0454 if (!logger().doPrint(lvl)) {
0455 return;
0456 }
0457 for (const auto& vt : volumes) {
0458 std::stringstream ss;
0459 ss << std::fixed;
0460 ss << std::setprecision(3);
0461 ss << std::setfill(' ');
0462
0463 int w = 9;
0464
0465 for (const auto& axis :
0466 {AxisDirection::AxisX, AxisDirection::AxisY, AxisDirection::AxisZ}) {
0467 ss << axisDirectionName(axis) << ": [ " << std::setw(w) << vt.min(axis)
0468 << " <- " << std::setw(w) << vt.mid(axis) << " -> " << std::setw(w)
0469 << vt.max(axis) << " ]\n";
0470 }
0471 logger().log(lvl, ss.str());
0472 }
0473 }
0474
0475 void CuboidVolumeStack::checkVolumeAlignment(
0476 const std::vector<VolumeTuple>& volumes, const Logger& logger) const {
0477 std::size_t n = 0;
0478 auto dirIdx = axisToIndex(m_direction);
0479 auto dirOrth1Idx = axisToIndex(m_dirOrth1);
0480 auto dirOrth2Idx = axisToIndex(m_dirOrth2);
0481
0482 for (auto& vt : volumes) {
0483 ACTS_VERBOSE("Checking volume #"
0484 << n << " at " << axisDirectionName(m_direction) << ": "
0485 << vt.localTransform.translation()[dirIdx]);
0486 ACTS_VERBOSE("- Local transform is:\n" << vt.localTransform.matrix());
0487
0488
0489 constexpr auto tolerance = s_onSurfaceTolerance;
0490
0491
0492
0493
0494 if ((vt.localTransform.rotation().matrix() - RotationMatrix3::Identity())
0495 .norm() > tolerance) {
0496 ACTS_ERROR("Volumes are not aligned: rotation is different");
0497 throw std::invalid_argument(
0498 "Volumes are not aligned: rotation is different");
0499 }
0500
0501 ACTS_VERBOSE(" -> Rotation is ok!");
0502
0503
0504 if (std::abs(vt.localTransform.translation()[dirOrth1Idx]) > tolerance ||
0505 std::abs(vt.localTransform.translation()[dirOrth2Idx]) > tolerance) {
0506 ACTS_ERROR("Volumes are not aligned: translation in "
0507 << axisDirectionName(m_dirOrth1) << " or "
0508 << axisDirectionName(m_dirOrth2));
0509 throw std::invalid_argument("Volumes are not aligned: translation in " +
0510 axisDirectionName(m_dirOrth1) + " or " +
0511 axisDirectionName(m_dirOrth2));
0512 }
0513 ACTS_VERBOSE(" -> Translation in " << axisDirectionName(m_dirOrth1) << "/"
0514 << axisDirectionName(m_dirOrth2)
0515 << " is ok!");
0516
0517 n++;
0518 }
0519 }
0520
0521 std::pair<double, double> CuboidVolumeStack::synchronizeBounds(
0522 std::vector<VolumeTuple>& volumes, const Logger& logger) {
0523 auto boundDirOrth1 = CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1);
0524 auto boundDirOrth2 = CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2);
0525
0526 const double maxHl1 =
0527 std::ranges::max_element(volumes, [boundDirOrth1](const auto& a,
0528 const auto& b) {
0529 return a.bounds->get(boundDirOrth1) < b.bounds->get(boundDirOrth1);
0530 })->bounds->get(boundDirOrth1);
0531 const double maxHl2 =
0532 std::ranges::max_element(volumes, [boundDirOrth2](const auto& a,
0533 const auto& b) {
0534 return a.bounds->get(boundDirOrth2) < b.bounds->get(boundDirOrth2);
0535 })->bounds->get(boundDirOrth2);
0536 ACTS_VERBOSE("Found: half length " << axisDirectionName(m_dirOrth1) << ":"
0537 << maxHl1 << ", half length "
0538 << axisDirectionName(m_dirOrth2) << ":"
0539 << maxHl2);
0540
0541 for (auto& vt : volumes) {
0542 vt.set({
0543 {boundDirOrth1, maxHl1},
0544 {boundDirOrth2, maxHl2},
0545 });
0546 }
0547
0548 return {maxHl1, maxHl2};
0549 }
0550
0551 void CuboidVolumeStack::update(const GeometryContext& gctx,
0552 std::shared_ptr<VolumeBounds> volbounds,
0553 std::optional<Transform3> transform,
0554 const Logger& logger) {
0555 ACTS_DEBUG(
0556 "Resizing CuboidVolumeStack with strategy: " << m_resizeStrategies.first);
0557 ACTS_DEBUG("Currently have " << m_volumes.size() << " children");
0558 ACTS_DEBUG(m_gaps.size() << " gaps");
0559 for (const auto& v : m_volumes) {
0560 ACTS_DEBUG(" - volume bounds: \n" << v->volumeBounds());
0561 ACTS_DEBUG(" transform: \n"
0562 << v->localToGlobalTransform(gctx).matrix());
0563 }
0564
0565 ACTS_DEBUG("New bounds are: \n" << *volbounds);
0566
0567 auto bounds = std::dynamic_pointer_cast<CuboidVolumeBounds>(volbounds);
0568 if (bounds == nullptr) {
0569 throw std::invalid_argument(
0570 "CuboidVolumeStack requires CuboidVolumeBounds");
0571 }
0572
0573 if (*bounds == volumeBounds()) {
0574 ACTS_VERBOSE("Bounds are the same, no resize needed");
0575 return;
0576 }
0577
0578 ACTS_VERBOSE("Group transform is:\n" << toString(m_groupTransform));
0579 ACTS_VERBOSE("Current transform is:\n"
0580 << toString(localToGlobalTransform(gctx)));
0581 if (transform.has_value()) {
0582 ACTS_VERBOSE("Input transform:\n" << toString(transform.value()));
0583 }
0584
0585 VolumeTuple oldVolume{gctx, *this, localToGlobalTransform(gctx)};
0586 VolumeTuple newVolume{gctx, *this, localToGlobalTransform(gctx)};
0587 newVolume.updatedBounds = std::make_shared<CuboidVolumeBounds>(*bounds);
0588 newVolume.globalTransform = transform.value_or(localToGlobalTransform(gctx));
0589 newVolume.localTransform =
0590 globalToLocalTransform(gctx) * newVolume.globalTransform;
0591
0592 if (!transform.has_value()) {
0593 ACTS_VERBOSE("Local transform does not change");
0594 } else {
0595 ACTS_VERBOSE("Local transform changes from\n"
0596 << m_groupTransform.matrix() << "\nto\n"
0597 << newVolume.localTransform.matrix());
0598 ACTS_VERBOSE("Checking transform consistency");
0599
0600 std::vector<VolumeTuple> volTemp{newVolume};
0601 checkVolumeAlignment(volTemp, logger);
0602 }
0603
0604 constexpr auto tolerance = s_onSurfaceTolerance;
0605 auto same = [](double a, double b) { return std::abs(a - b) < tolerance; };
0606
0607 for (const auto& dir : {m_direction, m_dirOrth1, m_dirOrth2}) {
0608 const double newMin = newVolume.min(dir);
0609 const double newMax = newVolume.max(dir);
0610 const double newMid = newVolume.mid(dir);
0611 const double newHl = newVolume.halfLength(dir);
0612
0613 const double oldMin = oldVolume.min(dir);
0614 const double oldMax = oldVolume.max(dir);
0615 const double oldMid = oldVolume.mid(dir);
0616 const double oldHl = oldVolume.halfLength(dir);
0617
0618 ACTS_VERBOSE("Previous bounds are: " << axisDirectionName(dir) << ": [ "
0619 << oldMin << " <- " << oldMid << " -> "
0620 << oldMax << " ] (" << oldHl << ")\n");
0621 ACTS_VERBOSE("New bounds are: " << axisDirectionName(dir) << ": [ "
0622 << newMin << " <- " << newMid << " -> "
0623 << newMax << " ] (" << newHl << ")\n");
0624
0625 if (!same(newMin, oldMin) && newMin > oldMin) {
0626 ACTS_ERROR("Shrinking the stack size in "
0627 << axisDirectionName(dir) << " is not supported: " << newMin
0628 << " -> " << oldMin);
0629 throw std::invalid_argument("Shrinking the stack in " +
0630 axisDirectionName(dir) + " is not supported");
0631 }
0632
0633 if (!same(newMax, oldMax) && newMax < oldMax) {
0634 ACTS_ERROR("Shrinking the stack size in "
0635 << axisDirectionName(dir) << " is not supported: " << newMax
0636 << " -> " << oldMax);
0637 throw std::invalid_argument("Shrinking the stack in " +
0638 axisDirectionName(dir) + " is not supported");
0639 }
0640 }
0641 auto isGap = [this](const Volume* vol) {
0642 return std::ranges::any_of(
0643 m_gaps, [&](const auto& gap) { return vol == gap.get(); });
0644 };
0645 ACTS_VERBOSE("Stack direction is " << axisDirectionName(m_direction));
0646
0647 std::vector<VolumeTuple> volumeTuples;
0648 volumeTuples.reserve(m_volumes.size());
0649 std::ranges::transform(m_volumes, std::back_inserter(volumeTuples),
0650 [this, &gctx](const auto& volume) {
0651 return VolumeTuple{gctx, *volume, m_groupTransform};
0652 });
0653
0654 ACTS_VERBOSE("*** Initial volume configuration:");
0655 printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG);
0656 for (const auto& dir : {m_dirOrth1, m_dirOrth2}) {
0657 if (!same(newVolume.min(dir), oldVolume.min(dir)) ||
0658 !same(newVolume.max(dir), oldVolume.max(dir))) {
0659 ACTS_VERBOSE("Resize all volumes to new " << axisDirectionName(dir)
0660 << " bounds");
0661 for (auto& volume : volumeTuples) {
0662 volume.set({{CuboidVolumeBounds::boundsFromAxisDirection(dir),
0663 newVolume.halfLength(dir)}});
0664 }
0665 ACTS_VERBOSE("*** Volume configuration after " << axisDirectionName(dir)
0666 << " resizing:");
0667 printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG);
0668 } else {
0669 ACTS_VERBOSE(axisDirectionName(dir)
0670 << " bounds are the same, no " << axisDirectionName(dir)
0671 << " resize needed");
0672 }
0673 }
0674
0675 if (same(newVolume.halfLength(m_direction),
0676 oldVolume.halfLength(m_direction))) {
0677 ACTS_VERBOSE("Halflength "
0678 << axisDirectionName(m_direction) << "is the same, no "
0679 << axisDirectionName(m_direction) << "resize needed");
0680 } else {
0681 auto dirIdx = axisToIndex(m_direction);
0682 auto boundDirIdx = CuboidVolumeBounds::boundsFromAxisDirection(m_direction);
0683 auto [firstStrategy, secondStrategy] = m_resizeStrategies;
0684 if (firstStrategy == VolumeResizeStrategy::Expand) {
0685 if (newVolume.min(m_direction) < oldVolume.min(m_direction)) {
0686 ACTS_VERBOSE("Expanding first volume to new "
0687 << axisDirectionName(m_direction) << "bounds");
0688
0689 auto& first = volumeTuples.front();
0690 double newMinFirst = newVolume.min(m_direction);
0691 double newMidFirst = (newMinFirst + first.max(m_direction)) / 2.0;
0692 double newHlFirst = (first.max(m_direction) - newMinFirst) / 2.0;
0693
0694 ACTS_VERBOSE(" -> first " << axisDirectionName(m_direction) << ": [ "
0695 << newMinFirst << " <- " << newMidFirst
0696 << " -> " << first.max(m_direction)
0697 << " ] (hl: " << newHlFirst << ")");
0698
0699 Translation3 translation(Vector3::Unit(dirIdx) * newMidFirst);
0700 first.set({{boundDirIdx, newHlFirst}});
0701 first.setLocalTransform(Transform3{translation}, m_groupTransform);
0702 }
0703
0704 if (newVolume.max(m_direction) > oldVolume.max(m_direction)) {
0705 ACTS_VERBOSE("Expanding last volume to new "
0706 << axisDirectionName(m_direction) << " bounds");
0707
0708 auto& last = volumeTuples.back();
0709 double newMaxLast = newVolume.max(m_direction);
0710 double newMidLast = (last.min(m_direction) + newMaxLast) / 2.0;
0711 double newHlLast = (newMaxLast - last.min(m_direction)) / 2.0;
0712
0713 ACTS_VERBOSE(" -> last " << axisDirectionName(m_direction) << ": [ "
0714 << last.min(m_direction) << " <- "
0715 << newMidLast << " -> " << newMaxLast
0716 << " ] (hl: " << newHlLast << ")");
0717
0718 Translation3 translation(Vector3::Unit(dirIdx) * newMidLast);
0719 last.set({{boundDirIdx, newHlLast}});
0720 last.setLocalTransform(Transform3{translation}, m_groupTransform);
0721 }
0722 } else if (firstStrategy == VolumeResizeStrategy::Gap) {
0723 ACTS_VERBOSE("Creating gap volumes to fill the new "
0724 << axisDirectionName(m_direction) << " bounds");
0725
0726 auto printGapDimensions = [&](const VolumeTuple& gap,
0727 const std::string& prefix = "") {
0728 for (const auto& dir : {m_direction, m_dirOrth1, m_dirOrth2}) {
0729 ACTS_VERBOSE(" -> gap" << prefix << ": " << axisDirectionName(dir)
0730 << ": [ " << gap.min(m_direction) << " <- "
0731 << gap.mid(dir) << " -> " << gap.max(dir)
0732 << " ]");
0733 }
0734 };
0735
0736 if (!same(newVolume.min(m_direction), oldVolume.min(m_direction)) &&
0737 newVolume.min(m_direction) < oldVolume.min(m_direction)) {
0738 double gap1Min = newVolume.min(m_direction);
0739 double gap1Max = oldVolume.min(m_direction);
0740 double gap1Hl = (gap1Max - gap1Min) / 2.0;
0741 double gap1P = (gap1Max + gap1Min) / 2.0;
0742
0743
0744 auto& candidate = volumeTuples.front();
0745 if (isGap(candidate.volume)) {
0746 ACTS_VERBOSE("~> Reusing existing gap volume at negative "
0747 << axisDirectionName(m_direction));
0748
0749 gap1Hl =
0750 candidate.bounds->get(
0751 CuboidVolumeBounds::boundsFromAxisDirection(m_direction)) +
0752 gap1Hl;
0753 gap1Max = gap1Min + gap1Hl * 2;
0754 gap1P = (gap1Max + gap1Min) / 2.0;
0755
0756 printGapDimensions(candidate, " before");
0757
0758 auto gap1Bounds = std::make_shared<CuboidVolumeBounds>(
0759 std::initializer_list<
0760 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0761 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0762 gap1Hl},
0763 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0764 newVolume.halfLength(m_dirOrth1)},
0765 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0766 newVolume.halfLength(m_dirOrth2)}});
0767 Translation3 gap1Translation(Vector3::Unit(dirIdx) * gap1P);
0768 Transform3 gap1Transform = m_groupTransform * gap1Translation;
0769 candidate.volume->update(gctx, std::move(gap1Bounds), gap1Transform);
0770 candidate = VolumeTuple{gctx, *candidate.volume, m_groupTransform};
0771 ACTS_VERBOSE("After:");
0772 printGapDimensions(candidate, " after ");
0773
0774 } else {
0775 ACTS_VERBOSE("~> Creating new gap volume at negative ");
0776 auto gap1Bounds = std::make_shared<CuboidVolumeBounds>(
0777 std::initializer_list<
0778 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0779 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0780 gap1Hl},
0781 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0782 newVolume.halfLength(m_dirOrth1)},
0783 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0784 newVolume.halfLength(m_dirOrth2)}});
0785 Translation3 gap1Translation(Vector3::Unit(dirIdx) * gap1P);
0786 Transform3 gap1Transform = m_groupTransform * gap1Translation;
0787 auto gap1 = addGapVolume(gap1Transform, std::move(gap1Bounds));
0788 volumeTuples.insert(volumeTuples.begin(),
0789 VolumeTuple{gctx, *gap1, m_groupTransform});
0790 printGapDimensions(volumeTuples.front());
0791 }
0792 }
0793
0794 if (!same(newVolume.max(m_direction), oldVolume.max(m_direction)) &&
0795 newVolume.max(m_direction) > oldVolume.max(m_direction)) {
0796 double gap2Min = oldVolume.max(m_direction);
0797 double gap2Max = newVolume.max(m_direction);
0798 double gap2Hl = (gap2Max - gap2Min) / 2.0;
0799 double gap2P = (gap2Max + gap2Min) / 2.0;
0800
0801
0802 auto& candidate = volumeTuples.back();
0803 if (isGap(candidate.volume)) {
0804 ACTS_VERBOSE("~> Reusing existing gap volume at positive ");
0805
0806 gap2Hl =
0807 candidate.bounds->get(
0808 CuboidVolumeBounds::boundsFromAxisDirection(m_direction)) +
0809 gap2Hl;
0810 gap2Min = newVolume.max(m_direction) - gap2Hl * 2;
0811 gap2P = (gap2Max + gap2Min) / 2.0;
0812
0813 auto gap2Bounds = std::make_shared<CuboidVolumeBounds>(
0814 std::initializer_list<
0815 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0816 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0817 gap2Hl},
0818 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0819 newVolume.halfLength(m_dirOrth1)},
0820 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0821 newVolume.halfLength(m_dirOrth2)}});
0822 Translation3 gap2Translation(Vector3::Unit(dirIdx) * gap2P);
0823 Transform3 gap2Transform = m_groupTransform * gap2Translation;
0824
0825 candidate.volume->update(gctx, std::move(gap2Bounds), gap2Transform);
0826 candidate = VolumeTuple{gctx, *candidate.volume, m_groupTransform};
0827 printGapDimensions(candidate, " after ");
0828 } else {
0829 ACTS_VERBOSE("~> Creating new gap volume at positive ");
0830 auto gap2Bounds = std::make_shared<CuboidVolumeBounds>(
0831 std::initializer_list<
0832 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0833 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0834 gap2Hl},
0835 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0836 newVolume.halfLength(m_dirOrth1)},
0837 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0838 newVolume.halfLength(m_dirOrth2)}});
0839 Translation3 gap2Translation(Vector3::Unit(dirIdx) * gap2P);
0840 Transform3 gap2Transform = m_groupTransform * gap2Translation;
0841 auto gap2 = addGapVolume(gap2Transform, std::move(gap2Bounds));
0842 volumeTuples.emplace_back(gctx, *gap2, m_groupTransform);
0843 printGapDimensions(volumeTuples.back());
0844 }
0845 }
0846 }
0847
0848 ACTS_VERBOSE("*** Volume configuration after "
0849 << axisDirectionName(m_direction) << " resizing:");
0850 printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG);
0851 }
0852
0853 ACTS_VERBOSE("Commit and update outer vector of volumes");
0854 m_volumes.clear();
0855 for (auto& vt : volumeTuples) {
0856 vt.commit(gctx, logger);
0857 m_volumes.push_back(vt.volume);
0858 }
0859
0860 Volume::update(gctx, std::move(bounds), newVolume.globalTransform, logger);
0861
0862 m_groupTransform = localToGlobalTransform(gctx);
0863 }
0864
0865 }
