From 05c2168ba470269b2cd01735b3f78b8b9351ca1a Mon Sep 17 00:00:00 2001 From: Matthew Leary Harris <44590310+Matthewharri@users.noreply.github.com> Date: Fri, 11 Oct 2024 06:26:30 -0400 Subject: [PATCH 1/5] feat: Add GeoSimplePolygonBrep to GeoModelToDetectorVolume (#3713) This PR adds back in the conversion method for `GeoSimplePolygonBrep` to `DetectorVolume`. This somehow got removed during a refactoring. Also included is a small bug fix that didn't propagate the volume transformation correctly for `GeoShapeShift` objects. Tagging @junggjo9 @paulgessinger @noemina --- Plugins/GeoModel/src/GeoModelToDetectorVolume.cpp | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/Plugins/GeoModel/src/GeoModelToDetectorVolume.cpp b/Plugins/GeoModel/src/GeoModelToDetectorVolume.cpp index 6303dec2c26..f784a12a990 100644 --- a/Plugins/GeoModel/src/GeoModelToDetectorVolume.cpp +++ b/Plugins/GeoModel/src/GeoModelToDetectorVolume.cpp @@ -46,6 +46,13 @@ Volume convertVolume(const Transform3& trf, const GeoShape& shape) { const GeoBox* box = dynamic_cast(&shape); bounds = std::make_shared( box->getXHalfLength(), box->getYHalfLength(), box->getZHalfLength()); + } else if (shape.typeID() == GeoSimplePolygonBrep::getClassTypeID()) { + const GeoSimplePolygonBrep* brep = + dynamic_cast(&shape); + double xmin{0}, xmax{0}, ymin{0}, ymax{0}, zmin{0}, zmax{0}; + brep->extent(xmin, ymin, zmin, xmax, ymax, zmax); + bounds = std::make_shared( + (xmax - xmin) / 2, (ymax - ymin) / 2, (zmax - zmin) / 2); } else if (shape.typeID() == GeoTrd::getClassTypeID()) { const GeoTrd* trd = dynamic_cast(&shape); double x1 = trd->getXHalfLength1(); @@ -117,9 +124,10 @@ Volume convertVolume(const Transform3& trf, const GeoShape& shape) { dynamic_cast(&shape); const GeoShape* shapeOp = shiftShape->getOp(); newTrf = trf * shiftShape->getX(); - return convertVolume(trf, *shapeOp); + return convertVolume(newTrf, *shapeOp); } else { - throw std::runtime_error("FATAL: Unsupported GeoModel shape"); + throw std::runtime_error("FATAL: Unsupported GeoModel shape: " + + shape.type()); } return Volume(newTrf, bounds); } From 16f6341afd79d63c9250bae48a86f0eaf6fde3fc Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Fri, 11 Oct 2024 15:21:25 +0200 Subject: [PATCH 2/5] fix: CylVolStack resizing issue (#3715) This fixes a bug where the local transform was not correctly synced after resizing. Also improves a number of assertions to not rely on floating point identity anymore Related to #3502 --- This pull request includes several changes to the `CylinderVolumeStack` and `Volume` classes to improve logging, add new parameters, and enhance the overlap checking and volume updating mechanisms. The most important changes include adding a logger parameter to several methods, updating method signatures, and improving overlap checking logic. ### Enhancements to logging and method signatures: * [`Core/include/Acts/Geometry/CylinderVolumeStack.hpp`](diffhunk://#diff-86daf525fefbaa344566ea4fc6493ca85afd3627922f04c3f16758fca8717a6eL88-R93): Added a `logger` parameter to the `update` method and modified its signature to include this parameter. [[1]](diffhunk://#diff-86daf525fefbaa344566ea4fc6493ca85afd3627922f04c3f16758fca8717a6eL88-R93) [[2]](diffhunk://#diff-86daf525fefbaa344566ea4fc6493ca85afd3627922f04c3f16758fca8717a6eR125-R130) * [`Core/include/Acts/Geometry/Volume.hpp`](diffhunk://#diff-beb0194dfdd77fcfe13c0b73c9e7790cbe7bae313f1723f7979a4e9beb5d16d5R16): Added a `logger` parameter to the `update` method and included the necessary import for the `Logger` class. [[1]](diffhunk://#diff-beb0194dfdd77fcfe13c0b73c9e7790cbe7bae313f1723f7979a4e9beb5d16d5R16) [[2]](diffhunk://#diff-beb0194dfdd77fcfe13c0b73c9e7790cbe7bae313f1723f7979a4e9beb5d16d5R87-R90) * [`Core/src/Geometry/CylinderVolumeStack.cpp`](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL70-R70): Updated the `commit` method to include a `logger` parameter and modified the `initializeOuterVolume` method to pass the `logger` parameter to the `update` method. [[1]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL70-R70) [[2]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL154-R156) [[3]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL188-R190) [[4]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL212-R215) [[5]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL244-R247) [[6]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL268-R272) ### Improvements to overlap checking: * [`Core/src/Geometry/CylinderVolumeStack.cpp`](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL286-R301): Enhanced the `overlapPrint` method to include a `direction` parameter and updated the overlap checking logic to handle different directions more robustly. [[1]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL286-R301) [[2]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cR312-R323) [[3]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL319-R338) [[4]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL347-R368) [[5]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL357-R386) [[6]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL379-R405) [[7]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL397-R424) [[8]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cR436-R442) [[9]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL446-R472) ### Volume update improvements: * [`Core/src/Geometry/CylinderVolumeStack.cpp`](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL625-R689): Updated the `update` method to provide detailed logging of the current and new bounds, and added checks to prevent shrinking the stack size. [[1]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL625-R689) [[2]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL676-R703) [[3]](diffhunk://#diff-39babeb41776345f10ffd90e329b67faf7a8a4f47f0abe7533f665442d52a44cL693-R754) --- .../Acts/Geometry/CylinderVolumeStack.hpp | 22 +- Core/include/Acts/Geometry/Volume.hpp | 5 +- Core/src/Geometry/CylinderVolumeStack.cpp | 202 +++++++++++------- Core/src/Geometry/GridPortalLink.cpp | 5 +- Core/src/Geometry/TrivialPortalLink.cpp | 2 +- Core/src/Geometry/Volume.cpp | 3 +- Core/src/Surfaces/CylinderBounds.cpp | 6 +- .../Geometry/CylinderVolumeStackTests.cpp | 51 +++++ 8 files changed, 193 insertions(+), 103 deletions(-) diff --git a/Core/include/Acts/Geometry/CylinderVolumeStack.hpp b/Core/include/Acts/Geometry/CylinderVolumeStack.hpp index 5e708a5c3c9..3fcb8f1039b 100644 --- a/Core/include/Acts/Geometry/CylinderVolumeStack.hpp +++ b/Core/include/Acts/Geometry/CylinderVolumeStack.hpp @@ -85,23 +85,12 @@ class CylinderVolumeStack : public Volume { /// construction. /// @param volbounds is the new bounds /// @param transform is the new transform - /// @pre The volume bounds need to be of type - /// @c CylinderVolumeBounds. - void update(std::shared_ptr volbounds, - std::optional transform = std::nullopt) override; - - /// Update the volume bounds and transform. This - /// will update the bounds of all volumes in the stack - /// to accommodate the new bounds and optionally create - /// gap volumes according to the resize strategy set during - /// construction. - /// @param newBounds is the new bounds - /// @param transform is the new transform /// @param logger is the logger /// @pre The volume bounds need to be of type /// @c CylinderVolumeBounds. - void update(std::shared_ptr newBounds, - std::optional transform, const Logger& logger); + void update(std::shared_ptr volbounds, + std::optional transform = std::nullopt, + const Logger& logger = getDummyLogger()) override; /// Access the gap volume that were created during attachment or resizing. /// @return the vector of gap volumes @@ -133,11 +122,12 @@ class CylinderVolumeStack : public Volume { Acts::Logging::Level lvl); /// Helper function that prints output helping in debugging overlaps + /// @param direction is the overlap check direction /// @param a is the first volume /// @param b is the second volume /// @param logger is the logger - static void overlapPrint(const VolumeTuple& a, const VolumeTuple& b, - const Logger& logger); + static void overlapPrint(BinningValue direction, const VolumeTuple& a, + const VolumeTuple& b, const Logger& logger); /// Helper function that checks if volumes are properly aligned /// for attachment. diff --git a/Core/include/Acts/Geometry/Volume.hpp b/Core/include/Acts/Geometry/Volume.hpp index 9217e758444..da8932b1eab 100644 --- a/Core/include/Acts/Geometry/Volume.hpp +++ b/Core/include/Acts/Geometry/Volume.hpp @@ -13,6 +13,7 @@ #include "Acts/Geometry/GeometryObject.hpp" #include "Acts/Utilities/BinningType.hpp" #include "Acts/Utilities/BoundingBox.hpp" +#include "Acts/Utilities/Logger.hpp" #include #include @@ -83,8 +84,10 @@ class Volume : public GeometryObject { /// Set the volume bounds and optionally also update the volume transform /// @param volbounds The volume bounds to be assigned /// @param transform The transform to be assigned, can be optional + /// @param logger A logger object to log messages virtual void update(std::shared_ptr volbounds, - std::optional transform = std::nullopt); + std::optional transform = std::nullopt, + const Logger& logger = Acts::getDummyLogger()); /// Construct bounding box for this shape /// @param envelope Optional envelope to add / subtract from min/max diff --git a/Core/src/Geometry/CylinderVolumeStack.cpp b/Core/src/Geometry/CylinderVolumeStack.cpp index cb1cfdeb822..c0346e5fe2f 100644 --- a/Core/src/Geometry/CylinderVolumeStack.cpp +++ b/Core/src/Geometry/CylinderVolumeStack.cpp @@ -67,7 +67,7 @@ struct CylinderVolumeStack::VolumeTuple { transformDirty = true; } - void commit() { + void commit(const Logger& logger) { // make a copy so we can't accidentally modify in-place auto copy = std::make_shared(*updatedBounds); @@ -76,7 +76,7 @@ struct CylinderVolumeStack::VolumeTuple { transform = globalTransform; } - volume->update(std::move(updatedBounds), transform); + volume->update(std::move(updatedBounds), transform, logger); bounds = copy.get(); updatedBounds = std::move(copy); transformDirty = false; @@ -151,7 +151,9 @@ void CylinderVolumeStack::initializeOuterVolume(BinningValue direction, const auto* cylBounds = dynamic_cast( &m_volumes.front()->volumeBounds()); assert(cylBounds != nullptr && "Volume bounds are not cylinder bounds"); - Volume::update(std::make_shared(*cylBounds)); + Volume::update(std::make_shared(*cylBounds), + std::nullopt, logger); + ACTS_VERBOSE("Transform is now: " << m_transform.matrix()); return; } @@ -185,7 +187,7 @@ void CylinderVolumeStack::initializeOuterVolume(BinningValue direction, << vt.localTransform.translation()[eZ]); ACTS_VERBOSE(*vt.updatedBounds); - vt.commit(); + vt.commit(logger); } ACTS_VERBOSE("*** Volume configuration after r synchronization:"); @@ -209,7 +211,8 @@ void CylinderVolumeStack::initializeOuterVolume(BinningValue direction, m_transform = m_groupTransform * Translation3{0, 0, midZ}; - Volume::update(std::make_shared(minR, maxR, hlZ)); + Volume::update(std::make_shared(minR, maxR, hlZ), + std::nullopt, logger); ACTS_DEBUG("Outer bounds are:\n" << volumeBounds()); ACTS_DEBUG("Outer transform / new group transform is:\n" << m_transform.matrix()); @@ -241,7 +244,7 @@ void CylinderVolumeStack::initializeOuterVolume(BinningValue direction, for (auto& vt : volumeTuples) { ACTS_VERBOSE("Updated bounds for volume at r: " << vt.midR()); ACTS_VERBOSE(*vt.updatedBounds); - vt.commit(); + vt.commit(logger); } ACTS_VERBOSE("*** Volume configuration after z synchronization:"); @@ -265,7 +268,8 @@ void CylinderVolumeStack::initializeOuterVolume(BinningValue direction, m_transform = m_groupTransform * Translation3{0, 0, midZ}; - Volume::update(std::make_shared(minR, maxR, hlZ)); + Volume::update(std::make_shared(minR, maxR, hlZ), + std::nullopt, logger); ACTS_DEBUG("Outer bounds are:\n" << volumeBounds()); ACTS_DEBUG("Outer transform is:\n" << m_transform.matrix()); @@ -283,25 +287,40 @@ void CylinderVolumeStack::initializeOuterVolume(BinningValue direction, } void CylinderVolumeStack::overlapPrint( - const CylinderVolumeStack::VolumeTuple& a, + BinningValue direction, const CylinderVolumeStack::VolumeTuple& a, const CylinderVolumeStack::VolumeTuple& b, const Logger& logger) { if (logger().doPrint(Acts::Logging::DEBUG)) { std::stringstream ss; ss << std::fixed; ss << std::setprecision(3); ss << std::setfill(' '); - ACTS_VERBOSE("Checking overlap between"); + int w = 9; - ss << " - " - << " z: [ " << std::setw(w) << a.minZ() << " <- " << std::setw(w) - << a.midZ() << " -> " << std::setw(w) << a.maxZ() << " ]"; - ACTS_VERBOSE(ss.str()); - ss.str(""); - ss << " - " - << " z: [ " << std::setw(w) << b.minZ() << " <- " << std::setw(w) - << b.midZ() << " -> " << std::setw(w) << b.maxZ() << " ]"; - ACTS_VERBOSE(ss.str()); + ACTS_VERBOSE("Checking overlap between"); + if (direction == BinningValue::binZ) { + ss << " - " + << " z: [ " << std::setw(w) << a.minZ() << " <- " << std::setw(w) + << a.midZ() << " -> " << std::setw(w) << a.maxZ() << " ]"; + ACTS_VERBOSE(ss.str()); + + ss.str(""); + ss << " - " + << " z: [ " << std::setw(w) << b.minZ() << " <- " << std::setw(w) + << b.midZ() << " -> " << std::setw(w) << b.maxZ() << " ]"; + ACTS_VERBOSE(ss.str()); + } else { + ss << " - " + << " r: [ " << std::setw(w) << a.minR() << " <-> " << std::setw(w) + << a.maxR() << " ]"; + ACTS_VERBOSE(ss.str()); + + ss.str(""); + ss << " - " + << " r: [ " << std::setw(w) << b.minR() << " <-> " << std::setw(w) + << b.maxR() << " ]"; + ACTS_VERBOSE(ss.str()); + } } } @@ -316,7 +335,7 @@ CylinderVolumeStack::checkOverlapAndAttachInZ( auto& a = volumes.at(i); auto& b = volumes.at(j); - overlapPrint(a, b, logger); + overlapPrint(BinningValue::binZ, a, b, logger); if (a.maxZ() > b.minZ()) { ACTS_ERROR(" -> Overlap in z"); @@ -344,8 +363,9 @@ CylinderVolumeStack::checkOverlapAndAttachInZ( << (aZMidNew - aHlZNew) << " <- " << aZMidNew << " -> " << (aZMidNew + aHlZNew) << "]"); - assert(a.minZ() == aZMidNew - aHlZNew && "Volume shrunk"); - assert(a.maxZ() <= aZMidNew + aHlZNew && "Volume shrunk"); + assert(std::abs(a.minZ() - (aZMidNew - aHlZNew)) < 1e-9 && + "Volume shrunk"); + assert(aHlZNew >= a.halfLengthZ() && "Volume shrunk"); ActsScalar bZMidNew = (b.minZ() + b.maxZ()) / 2.0 - gapWidth / 4.0; ActsScalar bHlZNew = b.halfLengthZ() + gapWidth / 4.0; @@ -354,15 +374,16 @@ CylinderVolumeStack::checkOverlapAndAttachInZ( << (bZMidNew - bHlZNew) << " <- " << bZMidNew << " -> " << (bZMidNew + bHlZNew) << "]"); - assert(b.minZ() >= bZMidNew - bHlZNew && "Volume shrunk"); - assert(b.maxZ() == bZMidNew + bHlZNew && "Volume shrunk"); + assert(bHlZNew >= b.halfLengthZ() && "Volume shrunk"); + assert(std::abs(b.maxZ() - (bZMidNew + bHlZNew)) < 1e-9 && + "Volume shrunk"); - a.localTransform = Translation3{0, 0, aZMidNew}; - a.volume->setTransform(m_groupTransform * a.localTransform); + a.setLocalTransform(Transform3{Translation3{0, 0, aZMidNew}}, + m_groupTransform); a.updatedBounds->set(CylinderVolumeBounds::eHalfLengthZ, aHlZNew); - b.localTransform = Translation3{0, 0, bZMidNew}; - b.volume->setTransform(m_groupTransform * b.localTransform); + b.setLocalTransform(Transform3{Translation3{0, 0, bZMidNew}}, + m_groupTransform); b.updatedBounds->set(CylinderVolumeBounds::eHalfLengthZ, bHlZNew); break; @@ -376,11 +397,12 @@ CylinderVolumeStack::checkOverlapAndAttachInZ( << (aZMidNew - aHlZNew) << " <- " << aZMidNew << " -> " << (aZMidNew + aHlZNew) << "]"); - assert(a.minZ() == aZMidNew - aHlZNew && "Volume shrunk"); - assert(a.maxZ() <= aZMidNew + aHlZNew && "Volume shrunk"); + assert(std::abs(a.minZ() - (aZMidNew - aHlZNew)) < 1e-9 && + "Volume shrunk"); + assert(aHlZNew >= a.halfLengthZ() && "Volume shrunk"); - a.localTransform = Translation3{0, 0, aZMidNew}; - a.volume->setTransform(m_groupTransform * a.localTransform); + a.setLocalTransform(Transform3{Translation3{0, 0, aZMidNew}}, + m_groupTransform); a.updatedBounds->set(CylinderVolumeBounds::eHalfLengthZ, aHlZNew); break; @@ -394,11 +416,12 @@ CylinderVolumeStack::checkOverlapAndAttachInZ( << (bZMidNew - bHlZNew) << " <- " << bZMidNew << " -> " << (bZMidNew + bHlZNew) << "]"); - assert(b.minZ() >= bZMidNew - bHlZNew && "Volume shrunk"); - assert(b.maxZ() == bZMidNew + bHlZNew && "Volume shrunk"); + assert(bHlZNew >= b.halfLengthZ() && "Volume shrunk"); + assert(std::abs(b.maxZ() - (bZMidNew + bHlZNew)) < 1e-9 && + "Volume shrunk"); - b.localTransform = Translation3{0, 0, bZMidNew}; - b.volume->setTransform(m_groupTransform * b.localTransform); + b.setLocalTransform(Transform3{Translation3{0, 0, bZMidNew}}, + m_groupTransform); b.updatedBounds->set(CylinderVolumeBounds::eHalfLengthZ, bHlZNew); break; } @@ -410,10 +433,13 @@ CylinderVolumeStack::checkOverlapAndAttachInZ( ACTS_VERBOSE(" - Gap half length: " << gapHlZ << " at z: " << gapMidZ); + ActsScalar minR = std::min(a.minR(), b.minR()); + ActsScalar maxR = std::max(a.maxR(), b.maxR()); + Transform3 gapLocalTransform{Translation3{0, 0, gapMidZ}}; Transform3 gapGlobalTransform = m_groupTransform * gapLocalTransform; - auto gapBounds = std::make_shared( - a.minR(), b.maxR(), gapHlZ); + auto gapBounds = + std::make_shared(minR, maxR, gapHlZ); auto gap = addGapVolume(gapGlobalTransform, gapBounds); gapVolumes.emplace_back(*gap, m_groupTransform); @@ -443,7 +469,7 @@ CylinderVolumeStack::checkOverlapAndAttachInR( auto& a = volumes.at(i); auto& b = volumes.at(j); - overlapPrint(a, b, logger); + overlapPrint(BinningValue::binR, a, b, logger); if (a.maxR() > b.minR()) { ACTS_ERROR(" -> Overlap in r"); @@ -622,44 +648,45 @@ std::pair CylinderVolumeStack::synchronizeZBounds( } void CylinderVolumeStack::update(std::shared_ptr volbounds, - std::optional transform) { - if (volbounds == nullptr) { - throw std::invalid_argument("New bounds are nullptr"); + std::optional transform, + const Logger& logger) { + ACTS_DEBUG( + "Resizing CylinderVolumeStack with strategy: " << m_resizeStrategy); + ACTS_DEBUG("Currently have " << m_volumes.size() << " children"); + ACTS_DEBUG(m_gaps.size() << " gaps"); + for (const auto& v : m_volumes) { + ACTS_DEBUG(" - volume bounds: \n" << v->volumeBounds()); + ACTS_DEBUG(" transform: \n" << v->transform().matrix()); } + + ACTS_DEBUG("New bounds are: \n" << *volbounds); + auto cylBounds = std::dynamic_pointer_cast(volbounds); if (cylBounds == nullptr) { throw std::invalid_argument( "CylinderVolumeStack requires CylinderVolumeBounds"); } - update(std::move(cylBounds), transform, - *Acts::getDefaultLogger("CYLSTACK", Logging::VERBOSE)); -} - -void CylinderVolumeStack::update( - std::shared_ptr newBounds, - std::optional transform, const Logger& logger) { - ACTS_DEBUG( - "Resizing CylinderVolumeStack with strategy: " << m_resizeStrategy); - if (newBounds == nullptr) { + if (cylBounds == nullptr) { throw std::invalid_argument("New bounds are nullptr"); } - if (*newBounds == volumeBounds()) { + if (*cylBounds == volumeBounds()) { ACTS_VERBOSE("Bounds are the same, no resize needed"); return; } + ACTS_VERBOSE("Group transform is:\n" << m_groupTransform.matrix()); + ACTS_VERBOSE("Current transform is:\n" << m_transform.matrix()); if (transform.has_value()) { - ACTS_VERBOSE(transform.value().matrix()); + ACTS_VERBOSE("Input transform:\n" << transform.value().matrix()); } VolumeTuple oldVolume{*this, m_transform}; VolumeTuple newVolume{*this, m_transform}; - newVolume.updatedBounds = std::make_shared(*newBounds); + newVolume.updatedBounds = std::make_shared(*cylBounds); newVolume.globalTransform = transform.value_or(m_transform); - newVolume.localTransform = - m_groupTransform.inverse() * newVolume.globalTransform; + newVolume.localTransform = m_transform.inverse() * newVolume.globalTransform; if (!transform.has_value()) { ACTS_VERBOSE("Local transform does not change"); @@ -673,7 +700,7 @@ void CylinderVolumeStack::update( checkVolumeAlignment(volTemp, logger); } - checkNoPhiOrBevel(*newBounds, logger); + checkNoPhiOrBevel(*cylBounds, logger); const ActsScalar newMinR = newVolume.minR(); const ActsScalar newMaxR = newVolume.maxR(); @@ -690,30 +717,41 @@ void CylinderVolumeStack::update( const ActsScalar oldHlZ = oldVolume.halfLengthZ(); ACTS_VERBOSE("Previous bounds are: z: [ " - << oldMinZ << " <- " << oldMidZ << " -> " << oldMaxZ - << " ], r: [ " << oldMinR << " <-> " << oldMaxR << " ]"); + << oldMinZ << " <- " << oldMidZ << " -> " << oldMaxZ << " ] (" + << oldHlZ << "), r: [ " << oldMinR << " <-> " << oldMaxR + << " ]"); ACTS_VERBOSE("New bounds are: z: [ " - << newMinZ << " <- " << newMidZ << " -> " << newMaxZ - << " ], r: [ " << newMinR << " <-> " << newMaxR << " ]"); - - if (newMinZ > oldMinZ) { - ACTS_ERROR("Shrinking the stack size in z is not supported"); - throw std::invalid_argument("Shrinking the stack is not supported"); + << newMinZ << " <- " << newMidZ << " -> " << newMaxZ << " ] (" + << newHlZ << "), r: [ " << newMinR << " <-> " << newMaxR + << " ]"); + + constexpr auto tolerance = s_onSurfaceTolerance; + auto same = [](ActsScalar a, ActsScalar b) { + return std::abs(a - b) < tolerance; + }; + + if (!same(newMinZ, oldMinZ) && newMinZ > oldMinZ) { + ACTS_ERROR("Shrinking the stack size in z is not supported: " + << newMinZ << " -> " << oldMinZ); + throw std::invalid_argument("Shrinking the stack in z is not supported"); } - if (newMaxZ < oldMaxZ) { - ACTS_ERROR("Shrinking the stack size in z is not supported"); - throw std::invalid_argument("Shrinking the stack is not supported"); + if (!same(newMaxZ, oldMaxZ) && newMaxZ < oldMaxZ) { + ACTS_ERROR("Shrinking the stack size in z is not supported: " + << newMaxZ << " -> " << oldMaxZ); + throw std::invalid_argument("Shrinking the stack in z is not supported"); } - if (newMinR > oldMinR) { - ACTS_ERROR("Shrinking the stack size in r is not supported"); - throw std::invalid_argument("Shrinking the stack is not supported"); + if (!same(newMinR, oldMinR) && newMinR > oldMinR) { + ACTS_ERROR("Shrinking the stack size in r is not supported: " + << newMinR << " -> " << oldMinR); + throw std::invalid_argument("Shrinking the stack in r is not supported"); } - if (newMaxR < oldMaxR) { - ACTS_ERROR("Shrinking the stack size in r is not supported"); - throw std::invalid_argument("Shrinking the stack is not supported"); + if (!same(newMaxR, oldMaxR) && newMaxR < oldMaxR) { + ACTS_ERROR("Shrinking the stack size in r is not supported: " + << newMaxR << " -> " << oldMaxR); + throw std::invalid_argument("Shrinking the stack is r in not supported"); } if (m_direction == BinningValue::binZ) { @@ -730,7 +768,7 @@ void CylinderVolumeStack::update( ACTS_VERBOSE("*** Initial volume configuration:"); printVolumeSequence(volumeTuples, logger, Acts::Logging::DEBUG); - if (newMinR != oldMinR || newMaxR != oldMaxR) { + if (!same(newMinR, oldMinR) || !same(newMaxR, oldMaxR)) { ACTS_VERBOSE("Resize all volumes to new r bounds"); for (auto& volume : volumeTuples) { volume.set({ @@ -744,7 +782,7 @@ void CylinderVolumeStack::update( ACTS_VERBOSE("R bounds are the same, no r resize needed"); } - if (newHlZ == oldHlZ) { + if (same(newHlZ, oldHlZ)) { ACTS_VERBOSE("Halflength z is the same, no z resize needed"); } else { if (m_resizeStrategy == ResizeStrategy::Expand) { @@ -784,7 +822,7 @@ void CylinderVolumeStack::update( } else if (m_resizeStrategy == ResizeStrategy::Gap) { ACTS_VERBOSE("Creating gap volumes to fill the new z bounds"); - if (newMinZ < oldMinZ) { + if (!same(newMinZ, oldMinZ) && newMinZ < oldMinZ) { ACTS_VERBOSE("Adding gap volume at negative z"); ActsScalar gap1MinZ = newVolume.minZ(); @@ -804,7 +842,7 @@ void CylinderVolumeStack::update( VolumeTuple{*gap1, m_groupTransform}); } - if (newMaxZ > oldMaxZ) { + if (!same(newMaxZ, oldMaxZ) && newMaxZ > oldMaxZ) { ACTS_VERBOSE("Adding gap volume at positive z"); ActsScalar gap2MinZ = oldVolume.maxZ(); @@ -831,7 +869,7 @@ void CylinderVolumeStack::update( ACTS_VERBOSE("Commit and update outer vector of volumes"); m_volumes.clear(); for (auto& vt : volumeTuples) { - vt.commit(); + vt.commit(logger); m_volumes.push_back(vt.volume); } @@ -923,13 +961,15 @@ void CylinderVolumeStack::update( ACTS_VERBOSE("Commit and update outer vector of volumes"); m_volumes.clear(); for (auto& vt : volumeTuples) { - vt.commit(); + vt.commit(logger); m_volumes.push_back(vt.volume); } } m_transform = newVolume.globalTransform; - Volume::update(std::move(newBounds)); + // @TODO: We probably can reuse m_transform + m_groupTransform = m_transform; + Volume::update(std::move(cylBounds), std::nullopt, logger); } void CylinderVolumeStack::checkNoPhiOrBevel(const CylinderVolumeBounds& bounds, diff --git a/Core/src/Geometry/GridPortalLink.cpp b/Core/src/Geometry/GridPortalLink.cpp index 2ad7ea11ae5..5e724b63120 100644 --- a/Core/src/Geometry/GridPortalLink.cpp +++ b/Core/src/Geometry/GridPortalLink.cpp @@ -89,7 +89,10 @@ void GridPortalLink::checkConsistency(const CylinderSurface& cyl) const { ActsScalar hlZ = cyl.bounds().get(CylinderBounds::eHalfLengthZ); if (!same(axis.getMin(), -hlZ) || !same(axis.getMax(), hlZ)) { throw std::invalid_argument( - "GridPortalLink: CylinderBounds: invalid length setup."); + "GridPortalLink: CylinderBounds: invalid length setup: " + + std::to_string(axis.getMin()) + " != " + std::to_string(-hlZ) + + " or " + std::to_string(axis.getMax()) + + " != " + std::to_string(hlZ)); } }; auto checkRPhi = [&cyl, same](const IAxis& axis) { diff --git a/Core/src/Geometry/TrivialPortalLink.cpp b/Core/src/Geometry/TrivialPortalLink.cpp index abcdb44929c..d0e9fd70b8a 100644 --- a/Core/src/Geometry/TrivialPortalLink.cpp +++ b/Core/src/Geometry/TrivialPortalLink.cpp @@ -39,7 +39,7 @@ Result TrivialPortalLink::resolveVolume( } void TrivialPortalLink::toStream(std::ostream& os) const { - os << "TrivialPortalLink"; + os << "TrivialPortalLinkvolumeName() << ">"; } } // namespace Acts diff --git a/Core/src/Geometry/Volume.cpp b/Core/src/Geometry/Volume.cpp index c1aefc96186..97a026b039f 100644 --- a/Core/src/Geometry/Volume.cpp +++ b/Core/src/Geometry/Volume.cpp @@ -79,7 +79,8 @@ void Volume::assignVolumeBounds(std::shared_ptr volbounds) { } void Volume::update(std::shared_ptr volbounds, - std::optional transform) { + std::optional transform, + const Logger& /*logger*/) { if (volbounds) { m_volumeBounds = std::move(volbounds); } diff --git a/Core/src/Surfaces/CylinderBounds.cpp b/Core/src/Surfaces/CylinderBounds.cpp index e7d5032c993..27123d29ec0 100644 --- a/Core/src/Surfaces/CylinderBounds.cpp +++ b/Core/src/Surfaces/CylinderBounds.cpp @@ -153,10 +153,12 @@ std::vector Acts::CylinderBounds::circleVertices( void Acts::CylinderBounds::checkConsistency() noexcept(false) { if (get(eR) <= 0.) { - throw std::invalid_argument("CylinderBounds: invalid radial setup."); + throw std::invalid_argument( + "CylinderBounds: invalid radial setup: radius is negative"); } if (get(eHalfLengthZ) <= 0.) { - throw std::invalid_argument("CylinderBounds: invalid length setup."); + throw std::invalid_argument( + "CylinderBounds: invalid length setup: half length is negative"); } if (get(eHalfPhiSector) <= 0. || get(eHalfPhiSector) > M_PI) { throw std::invalid_argument("CylinderBounds: invalid phi sector setup."); diff --git a/Tests/UnitTests/Core/Geometry/CylinderVolumeStackTests.cpp b/Tests/UnitTests/Core/Geometry/CylinderVolumeStackTests.cpp index ab3c23a3dfc..cdd56a22899 100644 --- a/Tests/UnitTests/Core/Geometry/CylinderVolumeStackTests.cpp +++ b/Tests/UnitTests/Core/Geometry/CylinderVolumeStackTests.cpp @@ -813,6 +813,57 @@ BOOST_DATA_TEST_CASE( } } +BOOST_AUTO_TEST_CASE(ResizeReproduction1) { + Transform3 trf1 = + Transform3::Identity() * Translation3{Vector3::UnitZ() * -2000}; + auto bounds1 = std::make_shared(70, 100, 100.0); + Volume vol1{trf1, bounds1}; + + std::vector volumes = {&vol1}; + CylinderVolumeStack stack(volumes, BinningValue::binZ, + CylinderVolumeStack::AttachmentStrategy::Gap, + CylinderVolumeStack::ResizeStrategy::Gap, *logger); + + Transform3 trf2 = + Transform3::Identity() * Translation3{Vector3::UnitZ() * -1500}; + stack.update(std::make_shared(30.0, 100, 600), trf2, + *logger); + + std::cout << stack.volumeBounds() << std::endl; + std::cout << stack.transform().matrix() << std::endl; + + Transform3 trf3 = + Transform3::Identity() * Translation3{Vector3::UnitZ() * -1600}; + stack.update(std::make_shared(30.0, 100, 700), trf3, + *logger); +} + +BOOST_AUTO_TEST_CASE(ResizeReproduction2) { + // The numbers are tuned a bit to reproduce the faulty behavior + Transform3 trf1 = + Transform3::Identity() * Translation3{Vector3::UnitZ() * 263}; + auto bounds1 = std::make_shared(30, 100, 4.075); + Volume vol1{trf1, bounds1}; + + std::vector volumes = {&vol1}; + CylinderVolumeStack stack(volumes, BinningValue::binZ, + CylinderVolumeStack::AttachmentStrategy::Gap, + CylinderVolumeStack::ResizeStrategy::Gap, *logger); + + Transform3 trf2 = + Transform3::Identity() * Translation3{Vector3::UnitZ() * 260.843}; + stack.update(std::make_shared(30.0, 100, 6.232), trf2, + *logger); + + std::cout << stack.volumeBounds() << std::endl; + std::cout << stack.transform().matrix() << std::endl; + + Transform3 trf3 = + Transform3::Identity() * Translation3{Vector3::UnitZ() * 1627.31}; + stack.update(std::make_shared(30.0, 100, 1372.699), + trf3, *logger); +} + BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE(RDirection) From 2b9d1ef8f0c66f4d3edc8d8ef7bce8511b77dc2b Mon Sep 17 00:00:00 2001 From: Carlo Varni <75478407+CarloVarni@users.noreply.github.com> Date: Fri, 11 Oct 2024 17:02:23 +0200 Subject: [PATCH 3/5] refactor: z and r axis in grid for seeding are Open instead of Bound (#3712) Not clear why these were `Bound` instead of `Open`. According to documentation (i.e. https://github.com/acts-project/acts/blob/main/Core/include/Acts/Utilities/AxisFwd.hpp) ```cpp enum class AxisBoundaryType { /// Default behaviour: out of bounds /// positions are filled into the over or underflow bins Open, /// Out-of-bounds positions resolve to first/last bin /// respectively Bound, /// Out-of-bounds positions resolve to the outermost /// bin on the opposite side Closed, }; ``` For the grid case, `Open` seems the proper option --- .../include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp | 4 ++-- .../include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp index 5a72449373a..7b020010771 100644 --- a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp +++ b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp @@ -23,8 +23,8 @@ template using CylindricalSpacePointGrid = Acts::Grid< std::vector, Acts::Axis, - Acts::Axis, - Acts::Axis>; + Acts::Axis, + Acts::Axis>; /// Cylindrical Binned Group template diff --git a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp index a21c08abd63..062aa7595ef 100644 --- a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp +++ b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp @@ -139,8 +139,8 @@ Acts::CylindricalSpacePointGridCreator::createGrid( config.rBinEdges.end()); } - Axis zAxis(std::move(zValues)); - Axis rAxis(std::move(rValues)); + Axis zAxis(std::move(zValues)); + Axis rAxis(std::move(rValues)); return Acts::CylindricalSpacePointGrid( std::make_tuple(std::move(phiAxis), std::move(zAxis), std::move(rAxis))); } From 5dd00138a848cc14e5ac917f5931eafebe19ae85 Mon Sep 17 00:00:00 2001 From: Carlo Varni <75478407+CarloVarni@users.noreply.github.com> Date: Fri, 11 Oct 2024 18:39:59 +0200 Subject: [PATCH 4/5] refactor: Only compute middle range once per bin (#3714) Seeding is done one bin at a time. Where the bin corresponds to the area where the middle space point candidate may be. That means that all middle space points processed by each call of `createSeedsForGroup` belong to the same bin and thus have the same z-bin. In case a `rRangeMiddleSP` vector is provided by the user (that is the case where the radius validity range changes according to the z-bin), we have to retrieve from the middle space point the z-bin and then the right validity range. In the current code this is done for every middle space point in the bin! We can do it only once. --- Core/include/Acts/Seeding/SeedFinder.hpp | 10 ++++ Core/include/Acts/Seeding/SeedFinder.ipp | 61 ++++++++++++------------ 2 files changed, 41 insertions(+), 30 deletions(-) diff --git a/Core/include/Acts/Seeding/SeedFinder.hpp b/Core/include/Acts/Seeding/SeedFinder.hpp index ecad26c540f..120610b36a1 100644 --- a/Core/include/Acts/Seeding/SeedFinder.hpp +++ b/Core/include/Acts/Seeding/SeedFinder.hpp @@ -107,6 +107,16 @@ class SeedFinder { const Acts::Range1D& rMiddleSPRange) const; private: + /// Given a middle space point candidate, get the proper radius validity range + /// In case the radius range changes according to the z-bin we need to + /// retrieve the proper range. We can do this computation only once, since + /// all the middle space point candidates belong to the same z-bin + /// @param spM space point candidate to be used as middle SP in a seed + /// @param rMiddleSPRange range object containing the minimum and maximum r for middle SP for a certain z bin. + std::pair retrieveRadiusRangeForMiddle( + const external_spacepoint_t& spM, + const Acts::Range1D& rMiddleSPRange) const; + /// Iterates over dublets and tests the compatibility between them by applying /// a series of cuts that can be tested with only two SPs /// @param options frequently changing configuration (like beam position) diff --git a/Core/include/Acts/Seeding/SeedFinder.ipp b/Core/include/Acts/Seeding/SeedFinder.ipp index 9f10c356259..f05e14fa0fb 100644 --- a/Core/include/Acts/Seeding/SeedFinder.ipp +++ b/Core/include/Acts/Seeding/SeedFinder.ipp @@ -104,40 +104,20 @@ void SeedFinder::createSeedsForGroup( return; } + // we compute this here since all middle space point candidates belong to the + // same z-bin + auto [minRadiusRangeForMiddle, maxRadiusRangeForMiddle] = + retrieveRadiusRangeForMiddle(*middleSPs.front(), rMiddleSPRange); for (const external_spacepoint_t* spM : middleSPs) { const float rM = spM->radius(); // check if spM is outside our radial region of interest - if (m_config.useVariableMiddleSPRange) { - if (rM < rMiddleSPRange.min()) { - continue; - } - if (rM > rMiddleSPRange.max()) { - // break because SPs are sorted in r - break; - } - } else if (!m_config.rRangeMiddleSP.empty()) { - /// get zBin position of the middle SP - auto pVal = std::lower_bound(m_config.zBinEdges.begin(), - m_config.zBinEdges.end(), spM->z()); - int zBin = std::distance(m_config.zBinEdges.begin(), pVal); - /// protects against zM at the limit of zBinEdges - zBin == 0 ? zBin : --zBin; - if (rM < m_config.rRangeMiddleSP[zBin][0]) { - continue; - } - if (rM > m_config.rRangeMiddleSP[zBin][1]) { - // break because SPs are sorted in r - break; - } - } else { - if (rM < m_config.rMinMiddle) { - continue; - } - if (rM > m_config.rMaxMiddle) { - // break because SPs are sorted in r - break; - } + if (rM < minRadiusRangeForMiddle) { + continue; + } + if (rM > maxRadiusRangeForMiddle) { + // break because SPs are sorted in r + break; } const float zM = spM->z(); @@ -827,4 +807,25 @@ SeedFinder::filterCandidates( } // loop on bottoms } +template +std::pair SeedFinder:: + retrieveRadiusRangeForMiddle( + const external_spacepoint_t& spM, + const Acts::Range1D& rMiddleSPRange) const { + if (m_config.useVariableMiddleSPRange) { + return std::make_pair(rMiddleSPRange.min(), rMiddleSPRange.max()); + } + if (!m_config.rRangeMiddleSP.empty()) { + /// get zBin position of the middle SP + auto pVal = std::lower_bound(m_config.zBinEdges.begin(), + m_config.zBinEdges.end(), spM.z()); + int zBin = std::distance(m_config.zBinEdges.begin(), pVal); + /// protects against zM at the limit of zBinEdges + zBin == 0 ? zBin : --zBin; + return std::make_pair(m_config.rRangeMiddleSP[zBin][0], + m_config.rRangeMiddleSP[zBin][1]); + } + return std::make_pair(m_config.rMinMiddle, m_config.rMaxMiddle); +} + } // namespace Acts From cf50c3cc931b747351c2a0e28fda36b6615ce1a3 Mon Sep 17 00:00:00 2001 From: Stephen Nicholas Swatman Date: Fri, 11 Oct 2024 19:53:13 +0200 Subject: [PATCH 5/5] fix: Make charge smearing optional in digi config (#3710) Currently, the JSON reading for digitization configurations unconditionally reads the charge smearing data, but this need not necessarily be present; the geometric config for the ODD, for example, does not have this data. As a result, parser was failing to parse the geometric ODD digitization config. This commit makes the reading conditional and adds tests to ensure the ODD config reading remains functional. --- .../Io/Json/src/JsonDigitizationConfig.cpp | 8 +++-- Examples/Python/tests/test_examples.py | 34 +++++++++++++++++++ 2 files changed, 40 insertions(+), 2 deletions(-) diff --git a/Examples/Io/Json/src/JsonDigitizationConfig.cpp b/Examples/Io/Json/src/JsonDigitizationConfig.cpp index 524fd4d32ef..489c86087bf 100644 --- a/Examples/Io/Json/src/JsonDigitizationConfig.cpp +++ b/Examples/Io/Json/src/JsonDigitizationConfig.cpp @@ -139,7 +139,9 @@ void ActsExamples::to_json(nlohmann::json& j, j["thickness"] = gdc.thickness; j["threshold"] = gdc.threshold; j["digital"] = gdc.digital; - to_json(j["charge-smearing"], gdc.chargeSmearer); + if (j.find("charge-smearing") != j.end()) { + to_json(j["charge-smearing"], gdc.chargeSmearer); + } } void ActsExamples::from_json(const nlohmann::json& j, @@ -161,7 +163,9 @@ void ActsExamples::from_json(const nlohmann::json& j, gdc.varianceMap[idx] = vars; } } - from_json(j["charge-smearing"], gdc.chargeSmearer); + if (j.find("charge-smearing") != j.end()) { + from_json(j["charge-smearing"], gdc.chargeSmearer); + } } void ActsExamples::to_json(nlohmann::json& j, diff --git a/Examples/Python/tests/test_examples.py b/Examples/Python/tests/test_examples.py index e396e17d0f7..46044a9513c 100644 --- a/Examples/Python/tests/test_examples.py +++ b/Examples/Python/tests/test_examples.py @@ -1000,6 +1000,40 @@ def test_digitization_example(trk_geo, tmp_path, assert_root_hash, digi_config_f assert_root_hash(root_file.name, root_file) +@pytest.mark.parametrize( + "digi_config_file", + [ + DIGI_SHARE_DIR / "default-smearing-config-generic.json", + DIGI_SHARE_DIR / "default-geometric-config-generic.json", + pytest.param( + ( + getOpenDataDetectorDirectory() + / "config" + / "odd-digi-smearing-config.json" + ), + marks=[ + pytest.mark.odd, + ], + ), + pytest.param( + ( + getOpenDataDetectorDirectory() + / "config" + / "odd-digi-geometric-config.json" + ), + marks=[ + pytest.mark.odd, + ], + ), + ], + ids=["smeared", "geometric", "odd-smeared", "odd-geometric"], +) +def test_digitization_example_input_parsing(digi_config_file): + from acts.examples import readDigiConfigFromJson + + acts.examples.readDigiConfigFromJson(str(digi_config_file)) + + @pytest.mark.parametrize( "digi_config_file", [