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