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