File indexing completed on 2026-07-16 08:07:47
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::copy(gapVolumes.begin(), gapVolumes.end(),
0208 std::back_inserter(volumeTuples));
0209
0210 ACTS_VERBOSE("*** Volume configuration after "
0211 << axisDirectionName(m_direction) << " attachment:");
0212 printVolumeSequence(volumeTuples, logger, Acts::Logging::VERBOSE);
0213
0214 ACTS_VERBOSE("Synchronizing bounds in " << axisDirectionName(m_dirOrth1)
0215 << "/"
0216 << axisDirectionName(m_dirOrth2));
0217 const auto [hl1, hl2] = synchronizeBounds(volumeTuples, logger);
0218
0219 for (auto& vt : volumeTuples) {
0220 ACTS_VERBOSE("Updated bounds for volume at "
0221 << axisDirectionName(m_direction) << ": "
0222 << vt.localTransform.translation()[dirIdx]);
0223 ACTS_VERBOSE(*vt.updatedBounds);
0224
0225 vt.commit(gctx, logger);
0226 }
0227
0228 ACTS_VERBOSE("*** Volume configuration after "
0229 << axisDirectionName(m_dirOrth1) << "/"
0230 << axisDirectionName(m_dirOrth2) << " synchronization:");
0231 printVolumeSequence(volumeTuples, logger, Acts::Logging::VERBOSE);
0232
0233 std::ranges::sort(volumeTuples, {},
0234 [*this](const auto& v) { return v.mid(m_direction); });
0235
0236 m_volumes.clear();
0237 for (const auto& vt : volumeTuples) {
0238 m_volumes.push_back(vt.volume);
0239 }
0240
0241 ACTS_DEBUG("*** Volume configuration after final "
0242 << axisDirectionName(m_direction) << " sorting:");
0243 printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG);
0244
0245 double min = volumeTuples.front().min(m_direction);
0246 double max = volumeTuples.back().max(m_direction);
0247
0248 double mid = std::midpoint(min, max);
0249 double hl = std::midpoint(max, -min);
0250
0251 Translation3 translation(Vector3::Unit(dirIdx) * mid);
0252 auto bounds = std::make_shared<CuboidVolumeBounds>(
0253 std::initializer_list<std::pair<CuboidVolumeBounds::BoundValues, double>>{
0254 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction), hl},
0255 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1), hl1},
0256 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2), hl2}});
0257 Volume::update(gctx, bounds, m_groupTransform * translation, logger);
0258 ACTS_DEBUG("Outer bounds are:\n" << volumeBounds());
0259 ACTS_DEBUG("Outer transform / new group transform is:\n"
0260 << toString(localToGlobalTransform(gctx)));
0261
0262
0263
0264 m_groupTransform = localToGlobalTransform(gctx);
0265 }
0266
0267 void CuboidVolumeStack::overlapPrint(const CuboidVolumeStack::VolumeTuple& a,
0268 const CuboidVolumeStack::VolumeTuple& b,
0269 const Logger& logger) {
0270 if (logger().doPrint(Acts::Logging::DEBUG)) {
0271 std::stringstream ss;
0272 ss << std::fixed;
0273 ss << std::setprecision(3);
0274 ss << std::setfill(' ');
0275
0276 int w = 9;
0277
0278 ACTS_VERBOSE("Checking overlap between");
0279 ss << " - " << " " << axisDirectionName(m_direction) << ": [ "
0280 << std::setw(w) << a.min(m_direction) << " <- " << std::setw(w)
0281 << a.mid(m_direction) << " -> " << std::setw(w) << a.max(m_direction)
0282 << " ]";
0283 ACTS_VERBOSE(ss.str());
0284
0285 ss.str("");
0286 ss << " - " << " " << axisDirectionName(m_direction) << ": [ "
0287 << std::setw(w) << b.min(m_direction) << " <- " << std::setw(w)
0288 << b.mid(m_direction) << " -> " << std::setw(w) << b.max(m_direction)
0289 << " ]";
0290 ACTS_VERBOSE(ss.str());
0291 }
0292 }
0293
0294 std::vector<CuboidVolumeStack::VolumeTuple>
0295 CuboidVolumeStack::checkOverlapAndAttach(const GeometryContext& gctx,
0296 std::vector<VolumeTuple>& volumes,
0297 VolumeAttachmentStrategy strategy,
0298 const Logger& logger) {
0299
0300 auto dirIdx = axisToIndex(m_direction);
0301 auto dirBoundIdx = CuboidVolumeBounds::boundsFromAxisDirection(m_direction);
0302
0303 std::vector<VolumeTuple> gapVolumes;
0304 for (std::size_t i = 0; i < volumes.size() - 1; i++) {
0305 std::size_t j = i + 1;
0306 auto& a = volumes.at(i);
0307 auto& b = volumes.at(j);
0308
0309 overlapPrint(a, b, logger);
0310
0311
0312 constexpr auto tolerance = s_onSurfaceTolerance;
0313 if (a.max(m_direction) - tolerance > b.min(m_direction)) {
0314 ACTS_ERROR(" -> Overlap in " << axisDirectionName(m_direction));
0315 throw std::invalid_argument("Volumes overlap in " +
0316 axisDirectionName(m_direction));
0317 } else {
0318 ACTS_VERBOSE(" -> No overlap");
0319 }
0320
0321 if (std::abs(a.max(m_direction) - b.min(m_direction)) < tolerance) {
0322 ACTS_VERBOSE("No gap between volumes, no attachment needed");
0323 } else {
0324 double gapWidth = b.min(m_direction) - a.max(m_direction);
0325 ACTS_VERBOSE("Gap width: " << gapWidth);
0326
0327 ACTS_VERBOSE("Synchronizing bounds in "
0328 << axisDirectionName(m_direction)
0329 << " with strategy: " << strategy);
0330 switch (strategy) {
0331 case VolumeAttachmentStrategy::Midpoint: {
0332 ACTS_VERBOSE(" -> Strategy: Expand both volumes to midpoint");
0333
0334 double aMidNew =
0335 (a.min(m_direction) + a.max(m_direction)) / 2.0 + gapWidth / 4.0;
0336 double aHlNew = a.halfLength(m_direction) + gapWidth / 4.0;
0337 ACTS_VERBOSE(" - New halflength for first volume: " << aHlNew);
0338 ACTS_VERBOSE(" - New bounds for first volume: ["
0339 << (aMidNew - aHlNew) << " <- " << aMidNew << " -> "
0340 << (aMidNew + aHlNew) << "]");
0341
0342 assert(std::abs(a.min(m_direction) - (aMidNew - aHlNew)) < 1e-9 &&
0343 "Volume shrunk");
0344 assert(aHlNew >= a.halfLength(m_direction) && "Volume shrunk");
0345
0346 double bMidNew =
0347 (b.min(m_direction) + b.max(m_direction)) / 2.0 - gapWidth / 4.0;
0348 double bHlNew = b.halfLength(m_direction) + gapWidth / 4.0;
0349 ACTS_VERBOSE(" - New halflength for second volume: " << bHlNew);
0350 ACTS_VERBOSE(" - New bounds for second volume: ["
0351 << (bMidNew - bHlNew) << " <- " << bMidNew << " -> "
0352 << (bMidNew + bHlNew) << "]");
0353
0354 assert(bHlNew >= b.halfLength(m_direction) && "Volume shrunk");
0355 assert(std::abs(b.max(m_direction) - (bMidNew + bHlNew)) < 1e-9 &&
0356 "Volume shrunk");
0357
0358 Translation3 translationA(Vector3::Unit(dirIdx) * aMidNew);
0359 a.setLocalTransform(Transform3{translationA}, m_groupTransform);
0360 a.updatedBounds->set(dirBoundIdx, aHlNew);
0361
0362 Translation3 translationB(Vector3::Unit(dirIdx) * bMidNew);
0363 b.setLocalTransform(Transform3{translationB}, m_groupTransform);
0364 b.updatedBounds->set(dirBoundIdx, bHlNew);
0365
0366 break;
0367 }
0368 case VolumeAttachmentStrategy::First: {
0369 ACTS_VERBOSE(" -> Strategy: Expand first volume");
0370 double aMidNew = (a.min(m_direction) + b.min(m_direction)) / 2.0;
0371 double aHlNew = (b.min(m_direction) - a.min(m_direction)) / 2.0;
0372 ACTS_VERBOSE(" - Gap width: " << gapWidth);
0373 ACTS_VERBOSE(" - New bounds for first volume: ["
0374 << (aMidNew - aHlNew) << " <- " << aMidNew << " -> "
0375 << (aMidNew + aHlNew) << "]");
0376
0377 assert(std::abs(a.min(m_direction) - (aMidNew - aHlNew)) < 1e-9 &&
0378 "Volume shrunk");
0379 assert(aHlNew >= a.halfLength(m_direction) && "Volume shrunk");
0380
0381 Translation3 translationA(Vector3::Unit(dirIdx) * aMidNew);
0382 a.setLocalTransform(Transform3{translationA}, m_groupTransform);
0383 a.updatedBounds->set(dirBoundIdx, aHlNew);
0384
0385 break;
0386 }
0387 case VolumeAttachmentStrategy::Second: {
0388 ACTS_VERBOSE(" -> Strategy: Expand second volume");
0389 double bMidNew = (a.max(m_direction) + b.max(m_direction)) / 2.0;
0390 double bHlNew = (b.max(m_direction) - a.max(m_direction)) / 2.0;
0391 ACTS_VERBOSE(" - New halflength for second volume: " << bHlNew);
0392 ACTS_VERBOSE(" - New bounds for second volume: ["
0393 << (bMidNew - bHlNew) << " <- " << bMidNew << " -> "
0394 << (bMidNew + bHlNew) << "]");
0395
0396 assert(bHlNew >= b.halfLength(m_direction) && "Volume shrunk");
0397 assert(std::abs(b.max(m_direction) - (bMidNew + bHlNew)) < 1e-9 &&
0398 "Volume shrunk");
0399
0400 Translation3 translationB(Vector3::Unit(dirIdx) * bMidNew);
0401 b.setLocalTransform(Transform3{translationB}, m_groupTransform);
0402 b.updatedBounds->set(dirBoundIdx, bHlNew);
0403 break;
0404 }
0405 case VolumeAttachmentStrategy::Gap: {
0406 ACTS_VERBOSE(" -> Strategy: Create a gap volume");
0407 double gapHl = (b.min(m_direction) - a.max(m_direction)) / 2.0;
0408 double gapMid = (b.min(m_direction) + a.max(m_direction)) / 2.0;
0409
0410 ACTS_VERBOSE(" - Gap half length: " << gapHl << " at "
0411 << axisDirectionName(m_direction)
0412 << ": " << gapMid);
0413
0414 Translation3 gapTranslation(Vector3::Unit(dirIdx) * gapMid);
0415
0416 double min1 = std::min(a.min(m_dirOrth1), b.min(m_dirOrth1));
0417 double max1 = std::max(a.max(m_dirOrth1), b.max(m_dirOrth1));
0418
0419 double min2 = std::min(a.min(m_dirOrth2), b.min(m_dirOrth2));
0420 double max2 = std::max(a.max(m_dirOrth2), b.max(m_dirOrth2));
0421
0422 Transform3 gapLocalTransform{gapTranslation};
0423 Transform3 gapGlobalTransform = m_groupTransform * gapLocalTransform;
0424
0425 auto gapBounds = std::make_shared<CuboidVolumeBounds>(
0426 std::initializer_list<
0427 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0428 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0429 gapHl},
0430 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0431 (max1 - min1) / 2},
0432 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0433 (max2 - min2) / 2}});
0434 auto gap = addGapVolume(gapGlobalTransform, gapBounds);
0435 gapVolumes.emplace_back(gctx, *gap, m_groupTransform);
0436
0437 break;
0438 }
0439 default:
0440 ACTS_ERROR("Attachment strategy " << strategy << " not implemented");
0441 std::stringstream ss;
0442 ss << strategy;
0443 throw std::invalid_argument("Attachment strategy " + ss.str() +
0444 " not implemented");
0445 }
0446 }
0447 }
0448
0449 return gapVolumes;
0450 }
0451
0452 void CuboidVolumeStack::printVolumeSequence(
0453 const std::vector<VolumeTuple>& volumes, const Logger& logger,
0454 Acts::Logging::Level lvl) {
0455 if (!logger().doPrint(lvl)) {
0456 return;
0457 }
0458 for (const auto& vt : volumes) {
0459 std::stringstream ss;
0460 ss << std::fixed;
0461 ss << std::setprecision(3);
0462 ss << std::setfill(' ');
0463
0464 int w = 9;
0465
0466 for (const auto& axis :
0467 {AxisDirection::AxisX, AxisDirection::AxisY, AxisDirection::AxisZ}) {
0468 ss << axisDirectionName(axis) << ": [ " << std::setw(w) << vt.min(axis)
0469 << " <- " << std::setw(w) << vt.mid(axis) << " -> " << std::setw(w)
0470 << vt.max(axis) << " ]\n";
0471 }
0472 logger().log(lvl, ss.str());
0473 }
0474 }
0475
0476 void CuboidVolumeStack::checkVolumeAlignment(
0477 const std::vector<VolumeTuple>& volumes, const Logger& logger) const {
0478 std::size_t n = 0;
0479 auto dirIdx = axisToIndex(m_direction);
0480 auto dirOrth1Idx = axisToIndex(m_dirOrth1);
0481 auto dirOrth2Idx = axisToIndex(m_dirOrth2);
0482
0483 for (auto& vt : volumes) {
0484 ACTS_VERBOSE("Checking volume #"
0485 << n << " at " << axisDirectionName(m_direction) << ": "
0486 << vt.localTransform.translation()[dirIdx]);
0487 ACTS_VERBOSE("- Local transform is:\n" << vt.localTransform.matrix());
0488
0489
0490 constexpr auto tolerance = s_onSurfaceTolerance;
0491
0492
0493
0494
0495 if ((vt.localTransform.rotation().matrix() - RotationMatrix3::Identity())
0496 .norm() > tolerance) {
0497 ACTS_ERROR("Volumes are not aligned: rotation is different");
0498 throw std::invalid_argument(
0499 "Volumes are not aligned: rotation is different");
0500 }
0501
0502 ACTS_VERBOSE(" -> Rotation is ok!");
0503
0504
0505 if (std::abs(vt.localTransform.translation()[dirOrth1Idx]) > tolerance ||
0506 std::abs(vt.localTransform.translation()[dirOrth2Idx]) > tolerance) {
0507 ACTS_ERROR("Volumes are not aligned: translation in "
0508 << axisDirectionName(m_dirOrth1) << " or "
0509 << axisDirectionName(m_dirOrth2));
0510 throw std::invalid_argument("Volumes are not aligned: translation in " +
0511 axisDirectionName(m_dirOrth1) + " or " +
0512 axisDirectionName(m_dirOrth2));
0513 }
0514 ACTS_VERBOSE(" -> Translation in " << axisDirectionName(m_dirOrth1) << "/"
0515 << axisDirectionName(m_dirOrth2)
0516 << " is ok!");
0517
0518 n++;
0519 }
0520 }
0521
0522 std::pair<double, double> CuboidVolumeStack::synchronizeBounds(
0523 std::vector<VolumeTuple>& volumes, const Logger& logger) {
0524 auto boundDirOrth1 = CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1);
0525 auto boundDirOrth2 = CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2);
0526
0527 const double maxHl1 =
0528 std::max_element(volumes.begin(), volumes.end(),
0529 [boundDirOrth1](const auto& a, const auto& b) {
0530 return a.bounds->get(boundDirOrth1) <
0531 b.bounds->get(boundDirOrth1);
0532 })
0533 ->bounds->get(boundDirOrth1);
0534 const double maxHl2 =
0535 std::max_element(volumes.begin(), volumes.end(),
0536 [boundDirOrth2](const auto& a, const auto& b) {
0537 return a.bounds->get(boundDirOrth2) <
0538 b.bounds->get(boundDirOrth2);
0539 })
0540 ->bounds->get(boundDirOrth2);
0541 ACTS_VERBOSE("Found: half length " << axisDirectionName(m_dirOrth1) << ":"
0542 << maxHl1 << ", half length "
0543 << axisDirectionName(m_dirOrth2) << ":"
0544 << maxHl2);
0545
0546 for (auto& vt : volumes) {
0547 vt.set({
0548 {boundDirOrth1, maxHl1},
0549 {boundDirOrth2, maxHl2},
0550 });
0551 }
0552
0553 return {maxHl1, maxHl2};
0554 }
0555
0556 void CuboidVolumeStack::update(const GeometryContext& gctx,
0557 std::shared_ptr<VolumeBounds> volbounds,
0558 std::optional<Transform3> transform,
0559 const Logger& logger) {
0560 ACTS_DEBUG(
0561 "Resizing CuboidVolumeStack with strategy: " << m_resizeStrategies.first);
0562 ACTS_DEBUG("Currently have " << m_volumes.size() << " children");
0563 ACTS_DEBUG(m_gaps.size() << " gaps");
0564 for (const auto& v : m_volumes) {
0565 ACTS_DEBUG(" - volume bounds: \n" << v->volumeBounds());
0566 ACTS_DEBUG(" transform: \n"
0567 << v->localToGlobalTransform(gctx).matrix());
0568 }
0569
0570 ACTS_DEBUG("New bounds are: \n" << *volbounds);
0571
0572 auto bounds = std::dynamic_pointer_cast<CuboidVolumeBounds>(volbounds);
0573 if (bounds == nullptr) {
0574 throw std::invalid_argument(
0575 "CuboidVolumeStack requires CuboidVolumeBounds");
0576 }
0577
0578 if (*bounds == volumeBounds()) {
0579 ACTS_VERBOSE("Bounds are the same, no resize needed");
0580 return;
0581 }
0582
0583 ACTS_VERBOSE("Group transform is:\n" << toString(m_groupTransform));
0584 ACTS_VERBOSE("Current transform is:\n"
0585 << toString(localToGlobalTransform(gctx)));
0586 if (transform.has_value()) {
0587 ACTS_VERBOSE("Input transform:\n" << toString(transform.value()));
0588 }
0589
0590 VolumeTuple oldVolume{gctx, *this, localToGlobalTransform(gctx)};
0591 VolumeTuple newVolume{gctx, *this, localToGlobalTransform(gctx)};
0592 newVolume.updatedBounds = std::make_shared<CuboidVolumeBounds>(*bounds);
0593 newVolume.globalTransform = transform.value_or(localToGlobalTransform(gctx));
0594 newVolume.localTransform =
0595 globalToLocalTransform(gctx) * newVolume.globalTransform;
0596
0597 if (!transform.has_value()) {
0598 ACTS_VERBOSE("Local transform does not change");
0599 } else {
0600 ACTS_VERBOSE("Local transform changes from\n"
0601 << m_groupTransform.matrix() << "\nto\n"
0602 << newVolume.localTransform.matrix());
0603 ACTS_VERBOSE("Checking transform consistency");
0604
0605 std::vector<VolumeTuple> volTemp{newVolume};
0606 checkVolumeAlignment(volTemp, logger);
0607 }
0608
0609 constexpr auto tolerance = s_onSurfaceTolerance;
0610 auto same = [](double a, double b) { return std::abs(a - b) < tolerance; };
0611
0612 for (const auto& dir : {m_direction, m_dirOrth1, m_dirOrth2}) {
0613 const double newMin = newVolume.min(dir);
0614 const double newMax = newVolume.max(dir);
0615 const double newMid = newVolume.mid(dir);
0616 const double newHl = newVolume.halfLength(dir);
0617
0618 const double oldMin = oldVolume.min(dir);
0619 const double oldMax = oldVolume.max(dir);
0620 const double oldMid = oldVolume.mid(dir);
0621 const double oldHl = oldVolume.halfLength(dir);
0622
0623 ACTS_VERBOSE("Previous bounds are: " << axisDirectionName(dir) << ": [ "
0624 << oldMin << " <- " << oldMid << " -> "
0625 << oldMax << " ] (" << oldHl << ")\n");
0626 ACTS_VERBOSE("New bounds are: " << axisDirectionName(dir) << ": [ "
0627 << newMin << " <- " << newMid << " -> "
0628 << newMax << " ] (" << newHl << ")\n");
0629
0630 if (!same(newMin, oldMin) && newMin > oldMin) {
0631 ACTS_ERROR("Shrinking the stack size in "
0632 << axisDirectionName(dir) << " is not supported: " << newMin
0633 << " -> " << oldMin);
0634 throw std::invalid_argument("Shrinking the stack in " +
0635 axisDirectionName(dir) + " is not supported");
0636 }
0637
0638 if (!same(newMax, oldMax) && newMax < oldMax) {
0639 ACTS_ERROR("Shrinking the stack size in "
0640 << axisDirectionName(dir) << " is not supported: " << newMax
0641 << " -> " << oldMax);
0642 throw std::invalid_argument("Shrinking the stack in " +
0643 axisDirectionName(dir) + " is not supported");
0644 }
0645 }
0646 auto isGap = [this](const Volume* vol) {
0647 return std::ranges::any_of(
0648 m_gaps, [&](const auto& gap) { return vol == gap.get(); });
0649 };
0650 ACTS_VERBOSE("Stack direction is " << axisDirectionName(m_direction));
0651
0652 std::vector<VolumeTuple> volumeTuples;
0653 volumeTuples.reserve(m_volumes.size());
0654 std::ranges::transform(m_volumes, std::back_inserter(volumeTuples),
0655 [this, &gctx](const auto& volume) {
0656 return VolumeTuple{gctx, *volume, m_groupTransform};
0657 });
0658
0659 ACTS_VERBOSE("*** Initial volume configuration:");
0660 printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG);
0661 for (const auto& dir : {m_dirOrth1, m_dirOrth2}) {
0662 if (!same(newVolume.min(dir), oldVolume.min(dir)) ||
0663 !same(newVolume.max(dir), oldVolume.max(dir))) {
0664 ACTS_VERBOSE("Resize all volumes to new " << axisDirectionName(dir)
0665 << " bounds");
0666 for (auto& volume : volumeTuples) {
0667 volume.set({{CuboidVolumeBounds::boundsFromAxisDirection(dir),
0668 newVolume.halfLength(dir)}});
0669 }
0670 ACTS_VERBOSE("*** Volume configuration after " << axisDirectionName(dir)
0671 << " resizing:");
0672 printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG);
0673 } else {
0674 ACTS_VERBOSE(axisDirectionName(dir)
0675 << " bounds are the same, no " << axisDirectionName(dir)
0676 << " resize needed");
0677 }
0678 }
0679
0680 if (same(newVolume.halfLength(m_direction),
0681 oldVolume.halfLength(m_direction))) {
0682 ACTS_VERBOSE("Halflength "
0683 << axisDirectionName(m_direction) << "is the same, no "
0684 << axisDirectionName(m_direction) << "resize needed");
0685 } else {
0686 auto dirIdx = axisToIndex(m_direction);
0687 auto boundDirIdx = CuboidVolumeBounds::boundsFromAxisDirection(m_direction);
0688 auto [firstStrategy, secondStrategy] = m_resizeStrategies;
0689 if (firstStrategy == VolumeResizeStrategy::Expand) {
0690 if (newVolume.min(m_direction) < oldVolume.min(m_direction)) {
0691 ACTS_VERBOSE("Expanding first volume to new "
0692 << axisDirectionName(m_direction) << "bounds");
0693
0694 auto& first = volumeTuples.front();
0695 double newMinFirst = newVolume.min(m_direction);
0696 double newMidFirst = (newMinFirst + first.max(m_direction)) / 2.0;
0697 double newHlFirst = (first.max(m_direction) - newMinFirst) / 2.0;
0698
0699 ACTS_VERBOSE(" -> first " << axisDirectionName(m_direction) << ": [ "
0700 << newMinFirst << " <- " << newMidFirst
0701 << " -> " << first.max(m_direction)
0702 << " ] (hl: " << newHlFirst << ")");
0703
0704 Translation3 translation(Vector3::Unit(dirIdx) * newMidFirst);
0705 first.set({{boundDirIdx, newHlFirst}});
0706 first.setLocalTransform(Transform3{translation}, m_groupTransform);
0707 }
0708
0709 if (newVolume.max(m_direction) > oldVolume.max(m_direction)) {
0710 ACTS_VERBOSE("Expanding last volume to new "
0711 << axisDirectionName(m_direction) << " bounds");
0712
0713 auto& last = volumeTuples.back();
0714 double newMaxLast = newVolume.max(m_direction);
0715 double newMidLast = (last.min(m_direction) + newMaxLast) / 2.0;
0716 double newHlLast = (newMaxLast - last.min(m_direction)) / 2.0;
0717
0718 ACTS_VERBOSE(" -> last " << axisDirectionName(m_direction) << ": [ "
0719 << last.min(m_direction) << " <- "
0720 << newMidLast << " -> " << newMaxLast
0721 << " ] (hl: " << newHlLast << ")");
0722
0723 Translation3 translation(Vector3::Unit(dirIdx) * newMidLast);
0724 last.set({{boundDirIdx, newHlLast}});
0725 last.setLocalTransform(Transform3{translation}, m_groupTransform);
0726 }
0727 } else if (firstStrategy == VolumeResizeStrategy::Gap) {
0728 ACTS_VERBOSE("Creating gap volumes to fill the new "
0729 << axisDirectionName(m_direction) << " bounds");
0730
0731 auto printGapDimensions = [&](const VolumeTuple& gap,
0732 const std::string& prefix = "") {
0733 for (const auto& dir : {m_direction, m_dirOrth1, m_dirOrth2}) {
0734 ACTS_VERBOSE(" -> gap" << prefix << ": " << axisDirectionName(dir)
0735 << ": [ " << gap.min(m_direction) << " <- "
0736 << gap.mid(dir) << " -> " << gap.max(dir)
0737 << " ]");
0738 }
0739 };
0740
0741 if (!same(newVolume.min(m_direction), oldVolume.min(m_direction)) &&
0742 newVolume.min(m_direction) < oldVolume.min(m_direction)) {
0743 double gap1Min = newVolume.min(m_direction);
0744 double gap1Max = oldVolume.min(m_direction);
0745 double gap1Hl = (gap1Max - gap1Min) / 2.0;
0746 double gap1P = (gap1Max + gap1Min) / 2.0;
0747
0748
0749 auto& candidate = volumeTuples.front();
0750 if (isGap(candidate.volume)) {
0751 ACTS_VERBOSE("~> Reusing existing gap volume at negative "
0752 << axisDirectionName(m_direction));
0753
0754 gap1Hl =
0755 candidate.bounds->get(
0756 CuboidVolumeBounds::boundsFromAxisDirection(m_direction)) +
0757 gap1Hl;
0758 gap1Max = gap1Min + gap1Hl * 2;
0759 gap1P = (gap1Max + gap1Min) / 2.0;
0760
0761 printGapDimensions(candidate, " before");
0762
0763 auto gap1Bounds = std::make_shared<CuboidVolumeBounds>(
0764 std::initializer_list<
0765 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0766 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0767 gap1Hl},
0768 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0769 newVolume.halfLength(m_dirOrth1)},
0770 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0771 newVolume.halfLength(m_dirOrth2)}});
0772 Translation3 gap1Translation(Vector3::Unit(dirIdx) * gap1P);
0773 Transform3 gap1Transform = m_groupTransform * gap1Translation;
0774 candidate.volume->update(gctx, std::move(gap1Bounds), gap1Transform);
0775 candidate = VolumeTuple{gctx, *candidate.volume, m_groupTransform};
0776 ACTS_VERBOSE("After:");
0777 printGapDimensions(candidate, " after ");
0778
0779 } else {
0780 ACTS_VERBOSE("~> Creating new gap volume at negative ");
0781 auto gap1Bounds = std::make_shared<CuboidVolumeBounds>(
0782 std::initializer_list<
0783 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0784 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0785 gap1Hl},
0786 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0787 newVolume.halfLength(m_dirOrth1)},
0788 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0789 newVolume.halfLength(m_dirOrth2)}});
0790 Translation3 gap1Translation(Vector3::Unit(dirIdx) * gap1P);
0791 Transform3 gap1Transform = m_groupTransform * gap1Translation;
0792 auto gap1 = addGapVolume(gap1Transform, std::move(gap1Bounds));
0793 volumeTuples.insert(volumeTuples.begin(),
0794 VolumeTuple{gctx, *gap1, m_groupTransform});
0795 printGapDimensions(volumeTuples.front());
0796 }
0797 }
0798
0799 if (!same(newVolume.max(m_direction), oldVolume.max(m_direction)) &&
0800 newVolume.max(m_direction) > oldVolume.max(m_direction)) {
0801 double gap2Min = oldVolume.max(m_direction);
0802 double gap2Max = newVolume.max(m_direction);
0803 double gap2Hl = (gap2Max - gap2Min) / 2.0;
0804 double gap2P = (gap2Max + gap2Min) / 2.0;
0805
0806
0807 auto& candidate = volumeTuples.back();
0808 if (isGap(candidate.volume)) {
0809 ACTS_VERBOSE("~> Reusing existing gap volume at positive ");
0810
0811 gap2Hl =
0812 candidate.bounds->get(
0813 CuboidVolumeBounds::boundsFromAxisDirection(m_direction)) +
0814 gap2Hl;
0815 gap2Min = newVolume.max(m_direction) - gap2Hl * 2;
0816 gap2P = (gap2Max + gap2Min) / 2.0;
0817
0818 auto gap2Bounds = std::make_shared<CuboidVolumeBounds>(
0819 std::initializer_list<
0820 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0821 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0822 gap2Hl},
0823 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0824 newVolume.halfLength(m_dirOrth1)},
0825 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0826 newVolume.halfLength(m_dirOrth2)}});
0827 Translation3 gap2Translation(Vector3::Unit(dirIdx) * gap2P);
0828 Transform3 gap2Transform = m_groupTransform * gap2Translation;
0829
0830 candidate.volume->update(gctx, std::move(gap2Bounds), gap2Transform);
0831 candidate = VolumeTuple{gctx, *candidate.volume, m_groupTransform};
0832 printGapDimensions(candidate, " after ");
0833 } else {
0834 ACTS_VERBOSE("~> Creating new gap volume at positive ");
0835 auto gap2Bounds = std::make_shared<CuboidVolumeBounds>(
0836 std::initializer_list<
0837 std::pair<CuboidVolumeBounds::BoundValues, double>>{
0838 {CuboidVolumeBounds::boundsFromAxisDirection(m_direction),
0839 gap2Hl},
0840 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth1),
0841 newVolume.halfLength(m_dirOrth1)},
0842 {CuboidVolumeBounds::boundsFromAxisDirection(m_dirOrth2),
0843 newVolume.halfLength(m_dirOrth2)}});
0844 Translation3 gap2Translation(Vector3::Unit(dirIdx) * gap2P);
0845 Transform3 gap2Transform = m_groupTransform * gap2Translation;
0846 auto gap2 = addGapVolume(gap2Transform, std::move(gap2Bounds));
0847 volumeTuples.emplace_back(gctx, *gap2, m_groupTransform);
0848 printGapDimensions(volumeTuples.back());
0849 }
0850 }
0851 }
0852
0853 ACTS_VERBOSE("*** Volume configuration after "
0854 << axisDirectionName(m_direction) << " resizing:");
0855 printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG);
0856 }
0857
0858 ACTS_VERBOSE("Commit and update outer vector of volumes");
0859 m_volumes.clear();
0860 for (auto& vt : volumeTuples) {
0861 vt.commit(gctx, logger);
0862 m_volumes.push_back(vt.volume);
0863 }
0864
0865 Volume::update(gctx, std::move(bounds), newVolume.globalTransform, logger);
0866
0867 m_groupTransform = localToGlobalTransform(gctx);
0868 }
0869
0870 }
