diff --git a/Core/include/Acts/Geometry/CompositePortalLink.hpp b/Core/include/Acts/Geometry/CompositePortalLink.hpp index 5ae6d76a106..6586d875b79 100644 --- a/Core/include/Acts/Geometry/CompositePortalLink.hpp +++ b/Core/include/Acts/Geometry/CompositePortalLink.hpp @@ -17,6 +17,9 @@ namespace Acts { +class GridPortalLink; +class Surface; + /// Composite portal links can graft together other portal link instances, for /// example grids that could not be merged due to invalid binnings. /// @@ -49,6 +52,15 @@ class CompositePortalLink final : public PortalLinkBase { std::unique_ptr b, BinningValue direction, bool flatten = true); + /// Construct a composite portal from any number of arbitrary other portal + /// links. The only requirement is that the portal link surfaces are + /// mergeable. + /// @param links The portal links + /// @param direction The binning direction + /// @param flatten If true, the composite will flatten any nested composite + CompositePortalLink(std::vector> links, + BinningValue direction, bool flatten = true); + /// Print the composite portal link /// @param os The output stream void toStream(std::ostream& os) const override; @@ -81,19 +93,22 @@ class CompositePortalLink final : public PortalLinkBase { /// @return The number of children std::size_t size() const; - private: - /// Helper function to construct a merged surface from two portal links along - /// a given direction - /// @param a The first portal link - /// @param b The second portal link - /// @param direction The merging direction - /// @return The merged surface - static std::shared_ptr mergedSurface(const PortalLinkBase* a, - const PortalLinkBase* b, - BinningValue direction); + /// (Potentially) create a grid portal link that represents this composite + /// portal link. + /// @note This only works, if the composite is **flat** and only contains + /// **trivial portal links**. If these preconditions are not met, this + /// function returns a nullptr. + /// @param gctx The geometry context + /// @param logger The logger + /// @return The grid portal link + std::unique_ptr makeGrid(const GeometryContext& gctx, + const Logger& logger) const; + private: boost::container::small_vector, 4> m_children{}; + + BinningValue m_direction; }; } // namespace Acts diff --git a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp index c92e001fa9f..1951d96b2da 100644 --- a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp +++ b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Definitions/Algebra.hpp" +#include "Acts/Geometry/BoundarySurfaceFace.hpp" #include "Acts/Geometry/Volume.hpp" #include "Acts/Geometry/VolumeBounds.hpp" #include "Acts/Utilities/BinningType.hpp" @@ -80,6 +81,18 @@ class CylinderVolumeBounds : public VolumeBounds { eSize }; + /// Enum describing the possible faces of a cylinder volume + /// @note These values are synchronized with the BoundarySurfaceFace enum. + /// Once Gen1 is removed, this can be changed. + enum class Face : unsigned int { + PositiveDisc = BoundarySurfaceFace::positiveFaceXY, + NegativeDisc = BoundarySurfaceFace::negativeFaceXY, + OuterCylinder = BoundarySurfaceFace::tubeOuterCover, + InnerCylinder = BoundarySurfaceFace::tubeInnerCover, + NegativePhiPlane = BoundarySurfaceFace::tubeSectorNegativePhi, + PositivePhiPlane = BoundarySurfaceFace::tubeSectorPositivePhi + }; + CylinderVolumeBounds() = delete; /// Constructor 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/GridPortalLink.hpp b/Core/include/Acts/Geometry/GridPortalLink.hpp index 0f84d496920..411e6bb9962 100644 --- a/Core/include/Acts/Geometry/GridPortalLink.hpp +++ b/Core/include/Acts/Geometry/GridPortalLink.hpp @@ -381,12 +381,12 @@ class GridPortalLink : public PortalLinkBase { /// Helper function to get grid bin content in type-eraased way. /// @param indices The bin indices /// @return The tracking volume at the bin - virtual TrackingVolume*& atLocalBins(IndexType indices) = 0; + virtual const TrackingVolume*& atLocalBins(IndexType indices) = 0; /// Helper function to get grid bin content in type-eraased way. /// @param indices The bin indices /// @return The tracking volume at the bin - virtual TrackingVolume* atLocalBins(IndexType indices) const = 0; + virtual const TrackingVolume* atLocalBins(IndexType indices) const = 0; private: BinningValue m_direction; @@ -399,7 +399,7 @@ template class GridPortalLinkT final : public GridPortalLink { public: /// The internal grid type - using GridType = Grid; + using GridType = Grid; /// The dimension of the grid static constexpr std::size_t DIM = sizeof...(Axes); @@ -530,7 +530,6 @@ class GridPortalLinkT final : public GridPortalLink { return m_grid.atPosition(m_projection(position)); } - protected: /// Type erased access to the number of bins /// @return The number of bins in each direction IndexType numLocalBins() const override { @@ -545,7 +544,7 @@ class GridPortalLinkT final : public GridPortalLink { /// Type erased local bin access /// @param indices The bin indices /// @return The tracking volume at the bin - TrackingVolume*& atLocalBins(IndexType indices) override { + const TrackingVolume*& atLocalBins(IndexType indices) override { throw_assert(indices.size() == DIM, "Invalid number of indices"); typename GridType::index_t idx; for (std::size_t i = 0; i < DIM; i++) { @@ -557,7 +556,7 @@ class GridPortalLinkT final : public GridPortalLink { /// Type erased local bin access /// @param indices The bin indices /// @return The tracking volume at the bin - TrackingVolume* atLocalBins(IndexType indices) const override { + const TrackingVolume* atLocalBins(IndexType indices) const override { throw_assert(indices.size() == DIM, "Invalid number of indices"); typename GridType::index_t idx; for (std::size_t i = 0; i < DIM; i++) { diff --git a/Core/include/Acts/Geometry/PortalShell.hpp b/Core/include/Acts/Geometry/PortalShell.hpp index 39e60b8774b..7405a64ad0a 100644 --- a/Core/include/Acts/Geometry/PortalShell.hpp +++ b/Core/include/Acts/Geometry/PortalShell.hpp @@ -8,7 +8,7 @@ #pragma once -#include "Acts/Geometry/BoundarySurfaceFace.hpp" +#include "Acts/Geometry/CylinderVolumeBounds.hpp" #include "Acts/Utilities/BinningType.hpp" #include "Acts/Utilities/Logger.hpp" @@ -68,19 +68,9 @@ class PortalShellBase { /// volumes class CylinderPortalShell : public PortalShellBase { public: - /// Enum describing the possible faces of a cylinder portal shell - /// @note These values are synchronized with the BoundarySurfaceFace enum. - /// Once Gen1 is removed, this can be changed. - enum class Face : unsigned int { - PositiveDisc = BoundarySurfaceFace::positiveFaceXY, - NegativeDisc = BoundarySurfaceFace::negativeFaceXY, - OuterCylinder = BoundarySurfaceFace::tubeOuterCover, - InnerCylinder = BoundarySurfaceFace::tubeInnerCover, - NegativePhiPlane = BoundarySurfaceFace::tubeSectorNegativePhi, - PositivePhiPlane = BoundarySurfaceFace::tubeSectorPositivePhi - }; - - using enum Face; + using Face = CylinderVolumeBounds::Face; + + using enum CylinderVolumeBounds::Face; /// Retrieve the portal associated to the given face. Can be nullptr if unset. /// @param face The face to retrieve the portal for diff --git a/Core/include/Acts/Geometry/TrivialPortalLink.hpp b/Core/include/Acts/Geometry/TrivialPortalLink.hpp index 9f99881fdc7..f68a7e326c6 100644 --- a/Core/include/Acts/Geometry/TrivialPortalLink.hpp +++ b/Core/include/Acts/Geometry/TrivialPortalLink.hpp @@ -58,6 +58,10 @@ class TrivialPortalLink final : public PortalLinkBase { const GeometryContext& gctx, const Vector3& position, double tolerance = s_onSurfaceTolerance) const override; + /// Get the single volume that this trivial portal link is associated with + /// @return The target volume + const TrackingVolume& volume() const; + private: TrackingVolume* m_volume; }; 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/include/Acts/Seeding/SeedFinder.hpp b/Core/include/Acts/Seeding/SeedFinder.hpp index 02047b88740..cfe9908574c 100644 --- a/Core/include/Acts/Seeding/SeedFinder.hpp +++ b/Core/include/Acts/Seeding/SeedFinder.hpp @@ -122,6 +122,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 b27dd710f76..2896fba925f 100644 --- a/Core/include/Acts/Seeding/SeedFinder.ipp +++ b/Core/include/Acts/Seeding/SeedFinder.ipp @@ -106,40 +106,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(); @@ -829,4 +809,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 diff --git a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp index 026a8f330a1..4410f1c125b 100644 --- a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp +++ b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp @@ -32,8 +32,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 4c6ea40541b..523d470eae5 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))); } diff --git a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp index afb74140b8d..05e236e1e29 100644 --- a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp +++ b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp @@ -34,6 +34,7 @@ #include #include #include +#include #include namespace Acts { @@ -579,8 +580,14 @@ class CombinatorialKalmanFilter { ACTS_ERROR("Error in filter: " << res.error()); result.lastError = res.error(); } + + if (result.finished) { + return; + } } + assert(!result.activeBranches.empty() && "No active branches"); + const bool isEndOfWorldReached = endOfWorldReached.checkAbort(state, stepper, navigator, logger()); const bool isVolumeConstraintReached = volumeConstraintAborter.checkAbort( @@ -589,9 +596,8 @@ class CombinatorialKalmanFilter { state, stepper, navigator, logger()); const bool isTargetReached = targetReached.checkAbort(state, stepper, navigator, logger()); - const bool allBranchesStopped = result.activeBranches.empty(); - if (isEndOfWorldReached || isPathLimitReached || isTargetReached || - allBranchesStopped || isVolumeConstraintReached) { + if (isEndOfWorldReached || isVolumeConstraintReached || + isPathLimitReached || isTargetReached) { if (isEndOfWorldReached) { ACTS_VERBOSE("End of world reached"); } else if (isVolumeConstraintReached) { @@ -620,26 +626,12 @@ class CombinatorialKalmanFilter { stepper.releaseStepSize(state.stepping, ConstrainedStep::actor); } - if (!allBranchesStopped) { - // Record the active branch and remove it from the list - storeLastActiveBranch(result); - result.activeBranches.pop_back(); - } else { - // This can happen if we stopped all branches in the filter step - ACTS_VERBOSE("All branches stopped"); - } + // Record the active branch and remove it from the list + storeLastActiveBranch(result); + result.activeBranches.pop_back(); - // If no more active branches, done with filtering; Otherwise, reset - // propagation state to track state at next active branch - if (!result.activeBranches.empty()) { - ACTS_VERBOSE("Propagation jumps to branch with tip = " - << result.activeBranches.back().tipIndex()); - reset(state, stepper, navigator, result); - } else { - ACTS_VERBOSE("Stop Kalman filtering with " - << result.collectedTracks.size() << " found tracks"); - result.finished = true; - } + // Reset propagation state to track state at next active branch + reset(state, stepper, navigator, result); } } @@ -665,9 +657,20 @@ class CombinatorialKalmanFilter { typename navigator_t> void reset(propagator_state_t& state, const stepper_t& stepper, const navigator_t& navigator, result_type& result) const { + if (result.activeBranches.empty()) { + ACTS_VERBOSE("Stop CKF with " << result.collectedTracks.size() + << " found tracks"); + result.finished = true; + + return; + } + auto currentBranch = result.activeBranches.back(); auto currentState = currentBranch.outermostTrackState(); + ACTS_VERBOSE("Propagation jumps to branch with tip = " + << currentBranch.tipIndex()); + // Reset the stepping state stepper.resetState(state.stepping, currentState.filtered(), currentState.filteredCovariance(), @@ -713,40 +716,67 @@ class CombinatorialKalmanFilter { using PM = TrackStatePropMask; bool isSensitive = surface->associatedDetectorElement() != nullptr; - bool isMaterial = surface->surfaceMaterial() != nullptr; - - std::size_t nBranchesOnSurface = 0; + bool hasMaterial = surface->surfaceMaterial() != nullptr; + bool isMaterialOnly = hasMaterial && !isSensitive; + bool expectMeasurements = isSensitive; - if (auto [slBegin, slEnd] = m_sourceLinkAccessor(*surface); - slBegin != slEnd) { - // Screen output message + if (isSensitive) { ACTS_VERBOSE("Measurement surface " << surface->geometryId() << " detected."); + } else if (isMaterialOnly) { + ACTS_VERBOSE("Material surface " << surface->geometryId() + << " detected."); + } else { + ACTS_VERBOSE("Passive surface " << surface->geometryId() + << " detected."); + return Result::success(); + } - // Transport the covariance to the surface + using SourceLinkRange = decltype(m_sourceLinkAccessor(*surface)); + std::optional slRange = std::nullopt; + bool hasMeasurements = false; + if (isSensitive) { + slRange = m_sourceLinkAccessor(*surface); + hasMeasurements = slRange->first != slRange->second; + } + bool isHole = isSensitive && !hasMeasurements; + + if (isHole) { + ACTS_VERBOSE("Detected hole before measurement selection on surface " + << surface->geometryId()); + } + + // Transport the covariance to the surface + if (isHole || isMaterialOnly) { + stepper.transportCovarianceToCurvilinear(state.stepping); + } else { stepper.transportCovarianceToBound(state.stepping, *surface); + } - // Update state and stepper with pre material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::PreUpdate); + // Update state and stepper with pre material effects + materialInteractor(surface, state, stepper, navigator, + MaterialUpdateStage::PreUpdate); - // Bind the transported state to the current surface - auto boundStateRes = - stepper.boundState(state.stepping, *surface, false); - if (!boundStateRes.ok()) { - return boundStateRes.error(); - } - auto& boundState = *boundStateRes; - auto& [boundParams, jacobian, pathLength] = boundState; - boundParams.covariance() = state.stepping.cov; + // Bind the transported state to the current surface + auto boundStateRes = stepper.boundState(state.stepping, *surface, false); + if (!boundStateRes.ok()) { + return boundStateRes.error(); + } + auto& boundState = *boundStateRes; + auto& [boundParams, jacobian, pathLength] = boundState; + boundParams.covariance() = state.stepping.cov; + + auto currentBranch = result.activeBranches.back(); + TrackIndexType prevTip = currentBranch.tipIndex(); - auto currentBranch = result.activeBranches.back(); - TrackIndexType prevTip = currentBranch.tipIndex(); + // Create trackstates for all source links (will be filtered later) + using TrackStatesResult = + Acts::Result>; + TrackStatesResult tsRes = TrackStatesResult::success({}); + if (hasMeasurements) { + auto [slBegin, slEnd] = *slRange; - // Create trackstates for all source links (will be filtered later) - using TrackStatesResult = - Acts::Result>; - TrackStatesResult tsRes = trackStateCandidateCreator( + tsRes = trackStateCandidateCreator( state.geoContext, *calibrationContextPtr, *surface, boundState, slBegin, slEnd, prevTip, *result.trackStates, result.trackStateCandidates, *result.trackStates, logger()); @@ -755,175 +785,98 @@ class CombinatorialKalmanFilter { "Processing of selected track states failed: " << tsRes.error()); return tsRes.error(); } - const CkfTypes::BranchVector& newTrackStateList = - *tsRes; - - if (newTrackStateList.empty()) { - ACTS_VERBOSE("Detected hole after measurement selection on surface " - << surface->geometryId()); - - // Setting the number of branches on the surface to 1 as the hole - // still counts as a branch - nBranchesOnSurface = 1; - - auto stateMask = PM::Predicted | PM::Jacobian; - - // Add a hole track state to the multitrajectory - TrackIndexType currentTip = addNonSourcelinkState( - stateMask, boundState, result, true, prevTip); - auto nonSourcelinkState = - result.trackStates->getTrackState(currentTip); - currentBranch.tipIndex() = currentTip; - currentBranch.nHoles()++; + } + const CkfTypes::BranchVector& newTrackStateList = *tsRes; + + if (!newTrackStateList.empty()) { + Result procRes = + processNewTrackStates(state.geoContext, newTrackStateList, result); + if (!procRes.ok()) { + ACTS_ERROR("Processing of selected track states failed: " + << procRes.error()); + return procRes.error(); + } + unsigned int nBranchesOnSurface = *procRes; - BranchStopperResult branchStopperResult = - m_extensions.branchStopper(currentBranch, nonSourcelinkState); + if (nBranchesOnSurface == 0) { + ACTS_VERBOSE("All branches on surface " << surface->geometryId() + << " have been stopped"); - // Check the branch - if (branchStopperResult == BranchStopperResult::Continue) { - // Remembered the active branch and its state - } else { - // No branch on this surface - nBranchesOnSurface = 0; - if (branchStopperResult == BranchStopperResult::StopAndKeep) { - storeLastActiveBranch(result); - } - // Remove the branch from list - result.activeBranches.pop_back(); - } - } else { - Result procRes = processNewTrackStates( - state.geoContext, newTrackStateList, result); - if (!procRes.ok()) { - ACTS_ERROR("Processing of selected track states failed: " - << procRes.error()); - return procRes.error(); - } - nBranchesOnSurface = *procRes; + reset(state, stepper, navigator, result); - if (nBranchesOnSurface == 0) { - // All branches on the surface have been stopped. Reset happens at - // the end of the function - ACTS_VERBOSE("All branches on surface " << surface->geometryId() - << " have been stopped"); - } else { - // `currentBranch` is invalidated after `processNewTrackStates` - currentBranch = result.activeBranches.back(); - prevTip = currentBranch.tipIndex(); - - auto currentState = currentBranch.outermostTrackState(); - - if (currentState.typeFlags().test(TrackStateFlag::OutlierFlag)) { - // We don't need to update the stepper given an outlier state - ACTS_VERBOSE("Outlier state detected on surface " - << surface->geometryId()); - } else { - // If there are measurement track states on this surface - ACTS_VERBOSE("Filtering step successful with " - << nBranchesOnSurface << " branches"); - // Update stepping state using filtered parameters of last track - // state on this surface - stepper.update(state.stepping, - MultiTrajectoryHelpers::freeFiltered( - state.options.geoContext, currentState), - currentState.filtered(), - currentState.filteredCovariance(), *surface); - ACTS_VERBOSE( - "Stepping state is updated with filtered parameter:"); - ACTS_VERBOSE("-> " << currentState.filtered().transpose() - << " of track state with tip = " - << currentState.index()); - } - } + return Result::success(); } - // Update state and stepper with post material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::PostUpdate); - } else if (isSensitive || isMaterial) { - ACTS_VERBOSE("Handle " << (isSensitive ? "sensitive" : "passive") - << " surface: " << surface->geometryId() - << " without measurements"); - - // No splitting on the surface without source links. Set it to one - // first, but could be changed later - nBranchesOnSurface = 1; - - auto currentBranch = result.activeBranches.back(); - TrackIndexType prevTip = currentBranch.tipIndex(); + // `currentBranch` is invalidated after `processNewTrackStates` + currentBranch = result.activeBranches.back(); + prevTip = currentBranch.tipIndex(); + } else { + if (expectMeasurements) { + ACTS_VERBOSE("Detected hole after measurement selection on surface " + << surface->geometryId()); + } - // No source links on surface, add either hole or passive material - // TrackState. No storage allocation for uncalibrated/calibrated - // measurement and filtered parameter auto stateMask = PM::Predicted | PM::Jacobian; - // Transport the covariance to a curvilinear surface - stepper.transportCovarianceToCurvilinear(state.stepping); - - // Update state and stepper with pre material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::PreUpdate); - - // Transport & bind the state to the current surface - auto boundStateRes = - stepper.boundState(state.stepping, *surface, false); - if (!boundStateRes.ok()) { - return boundStateRes.error(); - } - auto& boundState = *boundStateRes; - auto& [boundParams, jacobian, pathLength] = boundState; - boundParams.covariance() = state.stepping.cov; - // Add a hole or material track state to the multitrajectory TrackIndexType currentTip = addNonSourcelinkState( - stateMask, boundState, result, isSensitive, prevTip); - auto nonSourcelinkState = result.trackStates->getTrackState(currentTip); + stateMask, boundState, result, expectMeasurements, prevTip); currentBranch.tipIndex() = currentTip; - - if (isSensitive) { + auto currentState = currentBranch.outermostTrackState(); + if (expectMeasurements) { currentBranch.nHoles()++; } BranchStopperResult branchStopperResult = - m_extensions.branchStopper(currentBranch, nonSourcelinkState); + m_extensions.branchStopper(currentBranch, currentState); // Check the branch if (branchStopperResult == BranchStopperResult::Continue) { // Remembered the active branch and its state } else { // No branch on this surface - nBranchesOnSurface = 0; if (branchStopperResult == BranchStopperResult::StopAndKeep) { storeLastActiveBranch(result); } // Remove the branch from list result.activeBranches.pop_back(); - } - // Update state and stepper with post material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::PostUpdate); - } else { - // Neither measurement nor material on surface, this branch is still - // valid. Count the branch on current surface - nBranchesOnSurface = 1; - } + // Branch on the surface has been stopped - reset + ACTS_VERBOSE("Branch on surface " << surface->geometryId() + << " has been stopped"); - // Reset current branch if there is no branch on current surface - if (nBranchesOnSurface == 0) { - ACTS_DEBUG("Branch on surface " << surface->geometryId() - << " is stopped"); - if (!result.activeBranches.empty()) { - ACTS_VERBOSE("Propagation jumps to branch with tip = " - << result.activeBranches.back().tipIndex()); reset(state, stepper, navigator, result); - } else { - ACTS_VERBOSE("Stop Kalman filtering with " - << result.collectedTracks.size() << " found tracks"); - result.finished = true; + + return Result::success(); } } + auto currentState = currentBranch.outermostTrackState(); + + if (currentState.typeFlags().test(TrackStateFlag::OutlierFlag)) { + // We don't need to update the stepper given an outlier state + ACTS_VERBOSE("Outlier state detected on surface " + << surface->geometryId()); + } else if (currentState.typeFlags().test( + TrackStateFlag::MeasurementFlag)) { + // If there are measurement track states on this surface + // Update stepping state using filtered parameters of last track + // state on this surface + stepper.update(state.stepping, + MultiTrajectoryHelpers::freeFiltered( + state.options.geoContext, currentState), + currentState.filtered(), + currentState.filteredCovariance(), *surface); + ACTS_VERBOSE("Stepping state is updated with filtered parameter:"); + ACTS_VERBOSE("-> " << currentState.filtered().transpose() + << " of track state with tip = " + << currentState.index()); + } + + // Update state and stepper with post material effects + materialInteractor(surface, state, stepper, navigator, + MaterialUpdateStage::PostUpdate); + return Result::success(); } @@ -947,12 +900,13 @@ class CombinatorialKalmanFilter { auto rootBranch = result.activeBranches.back(); - // Build the new branches by forking the root branch. Reverse the order to - // process the best candidate first + // Build the new branches by forking the root branch. Reverse the order + // to process the best candidate first CkfTypes::BranchVector newBranches; for (auto it = newTrackStateList.rbegin(); it != newTrackStateList.rend(); ++it) { - // Keep the root branch as the first branch, make a copy for the others + // Keep the root branch as the first branch, make a copy for the + // others auto newBranch = (it == newTrackStateList.rbegin()) ? rootBranch : rootBranch.shallowCopy(); diff --git a/Core/include/Acts/TrackFinding/MeasurementSelector.ipp b/Core/include/Acts/TrackFinding/MeasurementSelector.ipp index 3e32c884585..30d09c002be 100644 --- a/Core/include/Acts/TrackFinding/MeasurementSelector.ipp +++ b/Core/include/Acts/TrackFinding/MeasurementSelector.ipp @@ -21,9 +21,9 @@ MeasurementSelector::select( ACTS_VERBOSE("Invoked MeasurementSelector"); - // Return error if no measurement + // Return if no measurement if (candidates.empty()) { - return CombinatorialKalmanFilterError::MeasurementSelectionFailed; + return Result::success(std::pair(candidates.begin(), candidates.end())); } // Get geoID of this surface @@ -92,20 +92,19 @@ MeasurementSelector::select( "No measurement candidate. Return an outlier measurement chi2=" << minChi2); isOutlier = true; - return Result::success(std::make_pair(candidates.begin() + minIndex, - candidates.begin() + minIndex + 1)); + return Result::success(std::pair(candidates.begin() + minIndex, + candidates.begin() + minIndex + 1)); } else { ACTS_VERBOSE("No measurement candidate. Return empty chi2=" << minChi2); - return Result::success( - std::make_pair(candidates.begin(), candidates.begin())); + return Result::success(std::pair(candidates.begin(), candidates.begin())); } } if (passedCandidates <= 1ul) { // return single item range, no sorting necessary ACTS_VERBOSE("Returning only 1 element chi2=" << minChi2); - return Result::success(std::make_pair(candidates.begin() + minIndex, - candidates.begin() + minIndex + 1)); + return Result::success(std::pair(candidates.begin() + minIndex, + candidates.begin() + minIndex + 1)); } std::sort( @@ -115,7 +114,7 @@ MeasurementSelector::select( ACTS_VERBOSE("Number of selected measurements: " << passedCandidates << ", max: " << cuts.numMeasurements); - return Result::success(std::make_pair( + return Result::success(std::pair( candidates.begin(), candidates.begin() + std::min(cuts.numMeasurements, passedCandidates))); } diff --git a/Core/src/Detector/detail/BlueprintDrawer.cpp b/Core/src/Detector/detail/BlueprintDrawer.cpp index bc17c4c5922..b518266dc4d 100644 --- a/Core/src/Detector/detail/BlueprintDrawer.cpp +++ b/Core/src/Detector/detail/BlueprintDrawer.cpp @@ -57,28 +57,35 @@ std::string labelStr( void Acts::Experimental::detail::BlueprintDrawer::dotStream( std::ostream& ss, const Acts::Experimental::Blueprint::Node& node, const Options& options) { + // Replace the "/" in node names + std::string nodeName = node.name; + std::replace(nodeName.begin(), nodeName.end(), '/', '_'); + // Root / leaf or branch if (node.isRoot()) { ss << "digraph " << options.graphName << " {" << '\n'; - ss << node.name << " " << labelStr(options.root, node.name, node.auxiliary) + ss << nodeName << " " << labelStr(options.root, nodeName, node.auxiliary) << '\n'; - ss << node.name << " " << shapeStr(options.root) << '\n'; + ss << nodeName << " " << shapeStr(options.root) << '\n'; } else if (node.isLeaf()) { - ss << node.name << " " << labelStr(options.leaf, node.name, node.auxiliary) + ss << nodeName << " " << labelStr(options.leaf, nodeName, node.auxiliary) << '\n'; - ss << node.name << " " + ss << nodeName << " " << ((node.internalsBuilder != nullptr) ? shapeStr(options.leaf) : shapeStr(options.gap)) << '\n'; } else { - ss << node.name << " " - << labelStr(options.branch, node.name, node.auxiliary) << '\n'; - ss << node.name << " " << shapeStr(options.branch) << '\n'; + ss << nodeName << " " << labelStr(options.branch, nodeName, node.auxiliary) + << '\n'; + ss << nodeName << " " << shapeStr(options.branch) << '\n'; } // Recursive for children for (const auto& c : node.children) { - ss << node.name << " -> " << c->name << ";" << '\n'; + // Replace the "/" in node names + std::string childName = c->name; + std::replace(childName.begin(), childName.end(), '/', '_'); + ss << nodeName << " -> " << childName << ";" << '\n'; dotStream(ss, *c, options); } @@ -86,29 +93,29 @@ void Acts::Experimental::detail::BlueprintDrawer::dotStream( Options::Node shape = node.isLeaf() ? options.shape : options.virtualShape; std::stringstream bts; bts << node.boundsType; - ss << node.name + "_shape " << shapeStr(shape) << '\n'; - ss << node.name + "_shape " + ss << nodeName + "_shape " << shapeStr(shape) << '\n'; + ss << nodeName + "_shape " << labelStr(shape, bts.str(), {"t = " + toString(node.transform.translation(), 1), "b = " + toString(node.boundaryValues, 1)}) << '\n'; - ss << node.name << " -> " << node.name + "_shape [ arrowhead = \"none\" ];" + ss << nodeName << " -> " << nodeName + "_shape [ arrowhead = \"none\" ];" << '\n'; // Sub node detection if (node.internalsBuilder != nullptr) { - ss << node.name + "_int " << shapeStr(options.internals) << '\n'; - ss << node.name << " -> " << node.name + "_int;" << '\n'; + ss << nodeName + "_int " << shapeStr(options.internals) << '\n'; + ss << nodeName << " -> " << nodeName + "_int;" << '\n'; } if (node.geoIdGenerator != nullptr) { - ss << node.name + "_geoID " << shapeStr(options.geoID) << '\n'; - ss << node.name << " -> " << node.name + "_geoID;" << '\n'; + ss << nodeName + "_geoID " << shapeStr(options.geoID) << '\n'; + ss << nodeName << " -> " << nodeName + "_geoID;" << '\n'; } if (node.rootVolumeFinderBuilder != nullptr) { - ss << node.name + "_roots " << shapeStr(options.roots) << '\n'; - ss << node.name << " -> " << node.name + "_roots;" << '\n'; + ss << nodeName + "_roots " << shapeStr(options.roots) << '\n'; + ss << nodeName << " -> " << nodeName + "_roots;" << '\n'; } if (node.isRoot()) { diff --git a/Core/src/Geometry/CompositePortalLink.cpp b/Core/src/Geometry/CompositePortalLink.cpp index 7d80e3b331e..d8010bbbce2 100644 --- a/Core/src/Geometry/CompositePortalLink.cpp +++ b/Core/src/Geometry/CompositePortalLink.cpp @@ -8,21 +8,74 @@ #include "Acts/Geometry/CompositePortalLink.hpp" +#include "Acts/Geometry/GridPortalLink.hpp" #include "Acts/Geometry/PortalError.hpp" +#include "Acts/Geometry/TrivialPortalLink.hpp" #include "Acts/Surfaces/CylinderSurface.hpp" #include "Acts/Surfaces/DiscSurface.hpp" #include "Acts/Surfaces/PlaneSurface.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" #include "Acts/Surfaces/RegularSurface.hpp" +#include "Acts/Utilities/Axis.hpp" +#include #include +#include #include +#include + namespace Acts { +namespace { +std::shared_ptr mergedSurface(const Surface& a, + const Surface& b, + BinningValue direction) { + if (a.type() != b.type()) { + throw std::invalid_argument{"Cannot merge surfaces of different types"}; + } + + if (const auto* cylinderA = dynamic_cast(&a); + cylinderA != nullptr) { + const auto& cylinderB = dynamic_cast(b); + + auto [merged, reversed] = cylinderA->mergedWith(cylinderB, direction, true); + return merged; + } else if (const auto* discA = dynamic_cast(&a); + discA != nullptr) { + const auto& discB = dynamic_cast(b); + auto [merged, reversed] = discA->mergedWith(discB, direction, true); + return merged; + } else if (const auto* planeA = dynamic_cast(&a); + planeA != nullptr) { + throw std::logic_error{"Plane surfaces not implemented yet"}; + } else { + throw std::invalid_argument{"Unsupported surface type"}; + } +} + +std::shared_ptr mergePortalLinks( + const std::vector>& links, + BinningValue direction) { + assert(std::ranges::all_of(links, + [](const auto& link) { return link != nullptr; })); + assert(!links.empty()); + + std::shared_ptr result = links.front()->surfacePtr(); + for (auto it = std::next(links.begin()); it != links.end(); ++it) { + assert(result != nullptr); + result = mergedSurface(*result, it->get()->surface(), direction); + } + + return result; +} +} // namespace + CompositePortalLink::CompositePortalLink(std::unique_ptr a, std::unique_ptr b, BinningValue direction, bool flatten) - : PortalLinkBase(mergedSurface(a.get(), b.get(), direction)) { + : PortalLinkBase(mergedSurface(a->surface(), b->surface(), direction)), + m_direction{direction} { if (!flatten) { m_children.push_back(std::move(a)); m_children.push_back(std::move(b)); @@ -45,6 +98,35 @@ CompositePortalLink::CompositePortalLink(std::unique_ptr a, } } +CompositePortalLink::CompositePortalLink( + std::vector> links, BinningValue direction, + bool flatten) + : PortalLinkBase(mergePortalLinks(links, direction)), + m_direction(direction) { + if (!flatten) { + for (auto& child : links) { + m_children.push_back(std::move(child)); + } + + } else { + auto handle = [&](std::unique_ptr link) { + if (auto* composite = dynamic_cast(link.get()); + composite != nullptr) { + m_children.insert( + m_children.end(), + std::make_move_iterator(composite->m_children.begin()), + std::make_move_iterator(composite->m_children.end())); + } else { + m_children.push_back(std::move(link)); + } + }; + + for (auto& child : links) { + handle(std::move(child)); + } + } +} + Result CompositePortalLink::resolveVolume( const GeometryContext& gctx, const Vector2& position, double tolerance) const { @@ -74,36 +156,6 @@ Result CompositePortalLink::resolveVolume( return PortalError::PositionNotOnAnyChildPortalLink; } -std::shared_ptr CompositePortalLink::mergedSurface( - const PortalLinkBase* a, const PortalLinkBase* b, BinningValue direction) { - assert(a != nullptr); - assert(b != nullptr); - if (a->surface().type() != b->surface().type()) { - throw std::invalid_argument{"Cannot merge surfaces of different types"}; - } - - if (const auto* cylinderA = - dynamic_cast(&a->surface()); - cylinderA != nullptr) { - const auto& cylinderB = dynamic_cast(b->surface()); - - auto [merged, reversed] = cylinderA->mergedWith(cylinderB, direction, true); - return merged; - } else if (const auto* discA = - dynamic_cast(&a->surface()); - discA != nullptr) { - const auto& discB = dynamic_cast(b->surface()); - auto [merged, reversed] = discA->mergedWith(discB, direction, true); - return merged; - } else if (const auto* planeA = - dynamic_cast(&a->surface()); - planeA != nullptr) { - throw std::logic_error{"Plane surfaces not implemented yet"}; - } else { - throw std::invalid_argument{"Unsupported surface type"}; - } -} - std::size_t CompositePortalLink::depth() const { std::size_t result = 1; @@ -125,4 +177,122 @@ void CompositePortalLink::toStream(std::ostream& os) const { os << "CompositePortalLink"; } +std::unique_ptr CompositePortalLink::makeGrid( + const GeometryContext& gctx, const Logger& logger) const { + ACTS_VERBOSE("Attempting to make a grid from a composite portal link"); + std::vector trivialLinks; + std::ranges::transform(m_children, std::back_inserter(trivialLinks), + [](const auto& child) { + return dynamic_cast(child.get()); + }); + + if (std::ranges::any_of(trivialLinks, + [](auto link) { return link == nullptr; })) { + ACTS_VERBOSE( + "Failed to make a grid from a composite portal link -> returning " + "nullptr"); + return nullptr; + } + + // we already know all children have surfaces that are compatible, we produced + // a merged surface for the overall dimensions. + + auto printEdges = [](const auto& edges) -> std::string { + std::vector strings; + std::ranges::transform(edges, std::back_inserter(strings), + [](const auto& v) { return std::to_string(v); }); + return boost::algorithm::join(strings, ", "); + }; + + if (surface().type() == Surface::SurfaceType::Cylinder) { + ACTS_VERBOSE("Combining composite into cylinder grid"); + + if (m_direction != BinningValue::binZ) { + ACTS_ERROR("Cylinder grid only supports binning in Z direction"); + throw std::runtime_error{"Unsupported binning direction"}; + } + + std::vector edges; + edges.reserve(m_children.size() + 1); + + const Transform3& groupTransform = m_surface->transform(gctx); + Transform3 itransform = groupTransform.inverse(); + + std::ranges::sort( + trivialLinks, [&itransform, &gctx](const auto& a, const auto& b) { + return (itransform * a->surface().transform(gctx)).translation()[eZ] < + (itransform * b->surface().transform(gctx)).translation()[eZ]; + }); + + for (const auto& [i, child] : enumerate(trivialLinks)) { + const auto& bounds = + dynamic_cast(child->surface().bounds()); + Transform3 ltransform = itransform * child->surface().transform(gctx); + ActsScalar hlZ = bounds.get(CylinderBounds::eHalfLengthZ); + ActsScalar minZ = ltransform.translation()[eZ] - hlZ; + ActsScalar maxZ = ltransform.translation()[eZ] + hlZ; + if (i == 0) { + edges.push_back(minZ); + } + edges.push_back(maxZ); + } + + ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges)); + + Axis axis{AxisBound, edges}; + + auto grid = GridPortalLink::make(m_surface, m_direction, std::move(axis)); + for (const auto& [i, child] : enumerate(trivialLinks)) { + grid->atLocalBins({i + 1}) = &child->volume(); + } + + return grid; + + } else if (surface().type() == Surface::SurfaceType::Disc) { + ACTS_VERBOSE("Combining composite into disc grid"); + + if (m_direction != BinningValue::binR) { + ACTS_ERROR("Disc grid only supports binning in R direction"); + throw std::runtime_error{"Unsupported binning direction"}; + } + + std::vector edges; + edges.reserve(m_children.size() + 1); + + std::ranges::sort(trivialLinks, [](const auto& a, const auto& b) { + const auto& boundsA = + dynamic_cast(a->surface().bounds()); + const auto& boundsB = + dynamic_cast(b->surface().bounds()); + return boundsA.get(RadialBounds::eMinR) < + boundsB.get(RadialBounds::eMinR); + }); + + for (const auto& [i, child] : enumerate(trivialLinks)) { + const auto& bounds = + dynamic_cast(child->surface().bounds()); + + if (i == 0) { + edges.push_back(bounds.get(RadialBounds::eMinR)); + } + edges.push_back(bounds.get(RadialBounds::eMaxR)); + } + + ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges)); + + Axis axis{AxisBound, edges}; + + auto grid = GridPortalLink::make(m_surface, m_direction, std::move(axis)); + for (const auto& [i, child] : enumerate(trivialLinks)) { + grid->atLocalBins({i + 1}) = &child->volume(); + } + + return grid; + } else if (surface().type() == Surface::SurfaceType::Plane) { + throw std::runtime_error{"Plane surfaces not implemented yet"}; + } else { + throw std::invalid_argument{"Unsupported surface type"}; + } +} + } // namespace Acts diff --git a/Core/src/Geometry/CylinderVolumeStack.cpp b/Core/src/Geometry/CylinderVolumeStack.cpp index cb1cfdeb822..9155087d57e 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,32 +717,48 @@ 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"); } + auto isGap = [this](const Volume* vol) { + return std::ranges::any_of( + m_gaps, [&](const auto& gap) { return vol == gap.get(); }); + }; + if (m_direction == BinningValue::binZ) { ACTS_VERBOSE("Stack direction is z"); @@ -730,7 +773,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 +787,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,43 +827,86 @@ void CylinderVolumeStack::update( } else if (m_resizeStrategy == ResizeStrategy::Gap) { ACTS_VERBOSE("Creating gap volumes to fill the new z bounds"); - if (newMinZ < oldMinZ) { - ACTS_VERBOSE("Adding gap volume at negative z"); + auto printGapDimensions = [&](const VolumeTuple& gap, + const std::string& prefix = "") { + ACTS_VERBOSE(" -> gap" << prefix << ": [ " << gap.minZ() << " <- " + << gap.midZ() << " -> " << gap.maxZ() + << " ], r: [ " << gap.minR() << " <-> " + << gap.maxR() << " ]"); + }; + if (!same(newMinZ, oldMinZ) && newMinZ < oldMinZ) { ActsScalar gap1MinZ = newVolume.minZ(); ActsScalar gap1MaxZ = oldVolume.minZ(); ActsScalar gap1HlZ = (gap1MaxZ - gap1MinZ) / 2.0; ActsScalar gap1PZ = (gap1MaxZ + gap1MinZ) / 2.0; - ACTS_VERBOSE(" -> gap1 z: [ " << gap1MinZ << " <- " << gap1PZ - << " -> " << gap1MaxZ - << " ] (hl: " << gap1HlZ << ")"); - - auto gap1Bounds = - std::make_shared(newMinR, newMaxR, gap1HlZ); - auto gap1Transform = m_groupTransform * Translation3{0, 0, gap1PZ}; - auto gap1 = addGapVolume(gap1Transform, std::move(gap1Bounds)); - volumeTuples.insert(volumeTuples.begin(), - VolumeTuple{*gap1, m_groupTransform}); + // // check if we need a new gap volume or reuse an existing one + auto& candidate = volumeTuples.front(); + if (isGap(candidate.volume)) { + ACTS_VERBOSE("~> Reusing existing gap volume at negative z"); + + gap1HlZ = + candidate.bounds->get(CylinderVolumeBounds::eHalfLengthZ) + + gap1HlZ; + gap1MaxZ = gap1MinZ + gap1HlZ * 2; + gap1PZ = (gap1MaxZ + gap1MinZ) / 2.0; + + printGapDimensions(candidate, " before"); + auto gap1Bounds = std::make_shared( + newMinR, newMaxR, gap1HlZ); + auto gap1Transform = m_groupTransform * Translation3{0, 0, gap1PZ}; + candidate.volume->update(std::move(gap1Bounds), gap1Transform); + candidate = VolumeTuple{*candidate.volume, m_groupTransform}; + ACTS_VERBOSE("After:"); + printGapDimensions(candidate, " after "); + + } else { + ACTS_VERBOSE("~> Creating new gap volume at negative z"); + auto gap1Bounds = std::make_shared( + newMinR, newMaxR, gap1HlZ); + auto gap1Transform = m_groupTransform * Translation3{0, 0, gap1PZ}; + auto gap1 = addGapVolume(gap1Transform, std::move(gap1Bounds)); + volumeTuples.insert(volumeTuples.begin(), + VolumeTuple{*gap1, m_groupTransform}); + printGapDimensions(volumeTuples.front()); + } } - if (newMaxZ > oldMaxZ) { - ACTS_VERBOSE("Adding gap volume at positive z"); - + if (!same(newMaxZ, oldMaxZ) && newMaxZ > oldMaxZ) { ActsScalar gap2MinZ = oldVolume.maxZ(); ActsScalar gap2MaxZ = newVolume.maxZ(); ActsScalar gap2HlZ = (gap2MaxZ - gap2MinZ) / 2.0; ActsScalar gap2PZ = (gap2MaxZ + gap2MinZ) / 2.0; - ACTS_VERBOSE(" -> gap2 z: [ " << gap2MinZ << " <- " << gap2PZ - << " -> " << gap2MaxZ - << " ] (hl: " << gap2HlZ << ")"); - - auto gap2Bounds = - std::make_shared(newMinR, newMaxR, gap2HlZ); - auto gap2Transform = m_groupTransform * Translation3{0, 0, gap2PZ}; - auto gap2 = addGapVolume(gap2Transform, std::move(gap2Bounds)); - volumeTuples.emplace_back(*gap2, m_groupTransform); + // check if we need a new gap volume or reuse an existing one + auto& candidate = volumeTuples.back(); + if (isGap(candidate.volume)) { + ACTS_VERBOSE("~> Reusing existing gap volume at positive z"); + + gap2HlZ = + candidate.bounds->get(CylinderVolumeBounds::eHalfLengthZ) + + gap2HlZ; + gap2MinZ = newVolume.maxZ() - gap2HlZ * 2; + gap2PZ = (gap2MaxZ + gap2MinZ) / 2.0; + + printGapDimensions(candidate, " before"); + auto gap2Bounds = std::make_shared( + newMinR, newMaxR, gap2HlZ); + auto gap2Transform = m_groupTransform * Translation3{0, 0, gap2PZ}; + + candidate.volume->update(std::move(gap2Bounds), gap2Transform); + candidate = VolumeTuple{*candidate.volume, m_groupTransform}; + printGapDimensions(candidate, " after "); + } else { + ACTS_VERBOSE("~> Creating new gap volume at positive z"); + auto gap2Bounds = std::make_shared( + newMinR, newMaxR, gap2HlZ); + auto gap2Transform = m_groupTransform * Translation3{0, 0, gap2PZ}; + auto gap2 = addGapVolume(gap2Transform, std::move(gap2Bounds)); + volumeTuples.emplace_back(*gap2, m_groupTransform); + printGapDimensions(volumeTuples.back()); + } } } @@ -831,7 +917,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); } @@ -887,32 +973,56 @@ void CylinderVolumeStack::update( << " ]"); } } else if (m_resizeStrategy == ResizeStrategy::Gap) { + auto printGapDimensions = [&](const VolumeTuple& gap, + const std::string& prefix = "") { + ACTS_VERBOSE(" -> gap" << prefix << ": [ " << gap.minZ() << " <- " + << gap.midZ() << " -> " << gap.maxZ() + << " ], r: [ " << gap.minR() << " <-> " + << gap.maxR() << " ]"); + }; + if (oldMinR > newMinR) { - ACTS_VERBOSE("Creating gap volume at the innermost end"); - auto gapBounds = - std::make_shared(newMinR, oldMinR, newHlZ); - auto gapTransform = m_groupTransform; - auto gapVolume = addGapVolume(gapTransform, gapBounds); - volumeTuples.insert(volumeTuples.begin(), - VolumeTuple{*gapVolume, m_groupTransform}); - auto gap = volumeTuples.front(); - ACTS_VERBOSE(" -> gap inner: [ " - << gap.minZ() << " <- " << gap.midZ() << " -> " - << gap.maxZ() << " ], r: [ " << gap.minR() << " <-> " - << gap.maxR() << " ]"); + auto& candidate = volumeTuples.front(); + if (isGap(candidate.volume)) { + ACTS_VERBOSE("~> Reusing existing gap volume at inner r"); + auto& candidateCylBounds = dynamic_cast( + candidate.volume->volumeBounds()); + printGapDimensions(candidate, " before"); + candidateCylBounds.set(CylinderVolumeBounds::eMinR, newMinR); + candidate = VolumeTuple{*candidate.volume, m_groupTransform}; + printGapDimensions(candidate, " after "); + } else { + ACTS_VERBOSE("~> Creating new gap volume at inner r"); + auto gapBounds = std::make_shared( + newMinR, oldMinR, newHlZ); + auto gapTransform = m_groupTransform; + auto gapVolume = addGapVolume(gapTransform, gapBounds); + volumeTuples.insert(volumeTuples.begin(), + VolumeTuple{*gapVolume, m_groupTransform}); + auto gap = volumeTuples.front(); + printGapDimensions(gap); + } } if (oldMaxR < newMaxR) { - ACTS_VERBOSE("Creating gap volume at the outermost end"); - auto gapBounds = - std::make_shared(oldMaxR, newMaxR, newHlZ); - auto gapTransform = m_groupTransform; - auto gapVolume = addGapVolume(gapTransform, gapBounds); - volumeTuples.emplace_back(*gapVolume, m_groupTransform); - auto gap = volumeTuples.back(); - ACTS_VERBOSE(" -> gap outer: [ " - << gap.minZ() << " <- " << gap.midZ() << " -> " - << gap.maxZ() << " ], r: [ " << gap.minR() << " <-> " - << gap.maxR() << " ]"); + auto& candidate = volumeTuples.back(); + if (isGap(candidate.volume)) { + ACTS_VERBOSE("~> Reusing existing gap volume at outer r"); + auto& candidateCylBounds = dynamic_cast( + candidate.volume->volumeBounds()); + printGapDimensions(candidate, " before"); + candidateCylBounds.set(CylinderVolumeBounds::eMaxR, newMaxR); + candidate = VolumeTuple{*candidate.volume, m_groupTransform}; + printGapDimensions(candidate, " after "); + } else { + ACTS_VERBOSE("~> Creating new gap volume at outer r"); + auto gapBounds = std::make_shared( + oldMaxR, newMaxR, newHlZ); + auto gapTransform = m_groupTransform; + auto gapVolume = addGapVolume(gapTransform, gapBounds); + volumeTuples.emplace_back(*gapVolume, m_groupTransform); + auto gap = volumeTuples.back(); + printGapDimensions(gap); + } } } @@ -923,13 +1033,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..57a3e1d3141 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) { @@ -295,7 +298,7 @@ void GridPortalLink::fillGrid1dTo2d(FillDirection dir, assert(locDest.size() == 2); for (std::size_t i = 0; i <= locSource[0] + 1; ++i) { - TrackingVolume* source = grid1d.atLocalBins({i}); + const TrackingVolume* source = grid1d.atLocalBins({i}); if (dir == FillDirection::loc1) { for (std::size_t j = 0; j <= locDest[1] + 1; ++j) { diff --git a/Core/src/Geometry/Portal.cpp b/Core/src/Geometry/Portal.cpp index 2a9d09b5f68..69cab13b8e0 100644 --- a/Core/src/Geometry/Portal.cpp +++ b/Core/src/Geometry/Portal.cpp @@ -9,13 +9,17 @@ #include "Acts/Geometry/Portal.hpp" #include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Tolerance.hpp" +#include "Acts/Geometry/CompositePortalLink.hpp" #include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/GridPortalLink.hpp" #include "Acts/Geometry/PortalLinkBase.hpp" #include "Acts/Geometry/TrivialPortalLink.hpp" #include "Acts/Surfaces/RegularSurface.hpp" #include "Acts/Utilities/BinningType.hpp" +#include "Acts/Utilities/Zip.hpp" -#include +#include #include #include @@ -245,7 +249,7 @@ Portal Portal::fuse(const GeometryContext& gctx, Portal& aPortal, Portal& bPortal, const Logger& logger) { ACTS_VERBOSE("Fusing two portals"); if (&aPortal == &bPortal) { - ACTS_ERROR("Cannot merge a portal with itself"); + ACTS_ERROR("Cannot fuse a portal with itself"); throw PortalMergingException{}; } @@ -267,6 +271,10 @@ Portal Portal::fuse(const GeometryContext& gctx, Portal& aPortal, if (!isSameSurface(gctx, *aPortal.m_surface, *bPortal.m_surface)) { ACTS_ERROR("Portals have different surfaces"); + ACTS_ERROR("A: " << aPortal.m_surface->bounds()); + ACTS_ERROR("\n" << aPortal.m_surface->transform(gctx).matrix()); + ACTS_ERROR("B: " << bPortal.m_surface->bounds()); + ACTS_ERROR("\n" << bPortal.m_surface->transform(gctx).matrix()); throw PortalFusingException(); } @@ -282,16 +290,27 @@ Portal Portal::fuse(const GeometryContext& gctx, Portal& aPortal, throw PortalFusingException(); } + auto maybeConvertToGrid = [&](std::unique_ptr link) + -> std::unique_ptr { + auto* composite = dynamic_cast(link.get()); + if (composite == nullptr) { + return link; + } + + ACTS_VERBOSE("Converting composite to grid during portal fusing"); + return composite->makeGrid(gctx, logger); + }; + aPortal.m_surface.reset(); bPortal.m_surface.reset(); if (aHasAlongNormal) { ACTS_VERBOSE("Taking along normal from lhs, opposite normal from rhs"); - return Portal{gctx, std::move(aPortal.m_alongNormal), - std::move(bPortal.m_oppositeNormal)}; + return Portal{gctx, maybeConvertToGrid(std::move(aPortal.m_alongNormal)), + maybeConvertToGrid(std::move(bPortal.m_oppositeNormal))}; } else { ACTS_VERBOSE("Taking along normal from rhs, opposite normal from lhs"); - return Portal{gctx, std::move(bPortal.m_alongNormal), - std::move(aPortal.m_oppositeNormal)}; + return Portal{gctx, maybeConvertToGrid(std::move(bPortal.m_alongNormal)), + maybeConvertToGrid(std::move(aPortal.m_oppositeNormal))}; } } @@ -305,12 +324,30 @@ bool Portal::isSameSurface(const GeometryContext& gctx, const Surface& a, return false; } - if (a.bounds() != b.bounds()) { + std::vector aValues = a.bounds().values(); + std::vector bValues = b.bounds().values(); + bool different = false; + for (auto [aVal, bVal] : zip(aValues, bValues)) { + if (std::abs(aVal - bVal) > s_onSurfaceTolerance) { + different = true; + break; + } + } + + if (a.bounds().type() != b.bounds().type() || different) { return false; } - if (!a.transform(gctx).isApprox(b.transform(gctx), - s_transformEquivalentTolerance)) { + if (!a.transform(gctx).linear().isApprox(b.transform(gctx).linear(), + s_transformEquivalentTolerance)) { + return false; + } + + Vector3 delta = + (a.transform(gctx).translation() - b.transform(gctx).translation()) + .cwiseAbs(); + + if (delta.maxCoeff() > s_onSurfaceTolerance) { return false; } diff --git a/Core/src/Geometry/PortalLinkBase.cpp b/Core/src/Geometry/PortalLinkBase.cpp index 459c144272e..bdc65f8a552 100644 --- a/Core/src/Geometry/PortalLinkBase.cpp +++ b/Core/src/Geometry/PortalLinkBase.cpp @@ -109,8 +109,10 @@ std::unique_ptr PortalLinkBase::merge( } else if (const auto* bTrivial = dynamic_cast(b.get()); bTrivial != nullptr) { - ACTS_VERBOSE("Merging a grid portal with a trivial portal"); - return gridMerge(*aGrid, *bTrivial->makeGrid(direction)); + ACTS_VERBOSE( + "Merging a grid portal with a trivial portal (via composite)"); + return std::make_unique(std::move(a), std::move(b), + direction); } else if (dynamic_cast(b.get()) != nullptr) { ACTS_WARNING("Merging a grid portal with a composite portal"); @@ -126,15 +128,17 @@ std::unique_ptr PortalLinkBase::merge( aTrivial != nullptr) { if (const auto* bGrid = dynamic_cast(b.get()); bGrid) { - ACTS_VERBOSE("Merging a trivial portal with a grid portal"); - return gridMerge(*aTrivial->makeGrid(direction), *bGrid); + ACTS_VERBOSE( + "Merging a trivial portal with a grid portal (via composite)"); + return std::make_unique(std::move(a), std::move(b), + direction); } else if (const auto* bTrivial = dynamic_cast(b.get()); bTrivial != nullptr) { - ACTS_VERBOSE("Merging two trivial portals"); - return gridMerge(*aTrivial->makeGrid(direction), - *bTrivial->makeGrid(direction)); + ACTS_VERBOSE("Merging two trivial portals (via composite"); + return std::make_unique(std::move(a), std::move(b), + direction); } else if (dynamic_cast(b.get()) != nullptr) { ACTS_WARNING("Merging a trivial portal with a composite portal"); diff --git a/Core/src/Geometry/PortalShell.cpp b/Core/src/Geometry/PortalShell.cpp index 59c83c1529f..b9baf403812 100644 --- a/Core/src/Geometry/PortalShell.cpp +++ b/Core/src/Geometry/PortalShell.cpp @@ -367,7 +367,7 @@ std::string CylinderStackPortalShell::label() const { std::ostream& operator<<(std::ostream& os, CylinderPortalShell::Face face) { switch (face) { - using enum CylinderPortalShell::Face; + using enum CylinderVolumeBounds::Face; case PositiveDisc: return os << "PositiveDisc"; case NegativeDisc: diff --git a/Core/src/Geometry/TrivialPortalLink.cpp b/Core/src/Geometry/TrivialPortalLink.cpp index abcdb44929c..6976388dbc0 100644 --- a/Core/src/Geometry/TrivialPortalLink.cpp +++ b/Core/src/Geometry/TrivialPortalLink.cpp @@ -39,7 +39,11 @@ Result TrivialPortalLink::resolveVolume( } void TrivialPortalLink::toStream(std::ostream& os) const { - os << "TrivialPortalLink"; + os << "TrivialPortalLinkvolumeName() << ">"; +} + +const TrackingVolume& TrivialPortalLink::volume() const { + return *m_volume; } } // 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/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/src/Base.cpp b/Examples/Python/src/Base.cpp index 4d5eec557bd..cc8a5f67dc1 100644 --- a/Examples/Python/src/Base.cpp +++ b/Examples/Python/src/Base.cpp @@ -364,12 +364,16 @@ void addBinning(Context& ctx) { .value("binY", Acts::BinningValue::binY) .value("binZ", Acts::BinningValue::binZ) .value("binR", Acts::BinningValue::binR) - .value("binPhi", Acts::BinningValue::binPhi); + .value("binPhi", Acts::BinningValue::binPhi) + .value("binRPhi", Acts::BinningValue::binRPhi) + .value("binH", Acts::BinningValue::binH) + .value("binEta", Acts::BinningValue::binEta) + .value("binMag", Acts::BinningValue::binMag); auto boundaryType = py::enum_(m, "AxisBoundaryType") - .value("bound", Acts::AxisBoundaryType::Bound) - .value("closed", Acts::AxisBoundaryType::Closed) - .value("open", Acts::AxisBoundaryType::Open); + .value("Bound", Acts::AxisBoundaryType::Bound) + .value("Closed", Acts::AxisBoundaryType::Closed) + .value("Open", Acts::AxisBoundaryType::Open); auto axisType = py::enum_(m, "AxisType") .value("equidistant", Acts::AxisType::Equidistant) diff --git a/Examples/Python/src/GeoModel.cpp b/Examples/Python/src/GeoModel.cpp index eb9b4406808..58609a573d5 100644 --- a/Examples/Python/src/GeoModel.cpp +++ b/Examples/Python/src/GeoModel.cpp @@ -192,8 +192,9 @@ void addGeoModel(Context& ctx) { .def_readwrite( "topBoundsOverride", &Acts::GeoModelBlueprintCreater::Options::topBoundsOverride) - .def_readwrite("table", - &Acts::GeoModelBlueprintCreater::Options::table); + .def_readwrite("table", &Acts::GeoModelBlueprintCreater::Options::table) + .def_readwrite("dotGraph", + &Acts::GeoModelBlueprintCreater::Options::dotGraph); } gm.def( diff --git a/Examples/Python/src/Geometry.cpp b/Examples/Python/src/Geometry.cpp index 8ec8ed97ef1..d676dbe6243 100644 --- a/Examples/Python/src/Geometry.cpp +++ b/Examples/Python/src/Geometry.cpp @@ -25,6 +25,7 @@ #include "Acts/Detector/interface/IInternalStructureBuilder.hpp" #include "Acts/Detector/interface/IRootVolumeFinderBuilder.hpp" #include "Acts/Geometry/CylinderVolumeBounds.hpp" +#include "Acts/Geometry/CylinderVolumeStack.hpp" #include "Acts/Geometry/Extent.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Geometry/GeometryHierarchyMap.hpp" @@ -36,6 +37,7 @@ #include "Acts/Plugins/Python/Utilities.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Surfaces/SurfaceArray.hpp" +#include "Acts/Utilities/BinningType.hpp" #include "Acts/Utilities/Helpers.hpp" #include "Acts/Utilities/RangeXD.hpp" #include "Acts/Visualization/ViewConfig.hpp" @@ -46,6 +48,7 @@ #include #include +#include #include #include @@ -106,7 +109,12 @@ void addGeometry(Context& ctx) { .def("approach", &Acts::GeometryIdentifier::approach) .def("sensitive", &Acts::GeometryIdentifier::sensitive) .def("extra", &Acts::GeometryIdentifier::extra) - .def("value", &Acts::GeometryIdentifier::value); + .def("value", &Acts::GeometryIdentifier::value) + .def("__str__", [](const Acts::GeometryIdentifier& self) { + std::stringstream ss; + ss << self; + return ss.str(); + }); } { @@ -178,15 +186,31 @@ void addGeometry(Context& ctx) { { py::class_>( - m, "VolumeBounds"); - - py::class_, Acts::VolumeBounds>( - m, "CylinderVolumeBounds") - .def(py::init(), - "rmin"_a, "rmax"_a, "halfz"_a, "halfphi"_a = M_PI, "avgphi"_a = 0., - "bevelMinZ"_a = 0., "bevelMaxZ"_a = 0.); + m, "VolumeBounds") + .def("type", &Acts::VolumeBounds::type) + .def("__str__", [](const Acts::VolumeBounds& self) { + std::stringstream ss; + ss << self; + return ss.str(); + }); + + auto cvb = + py::class_, + Acts::VolumeBounds>(m, "CylinderVolumeBounds") + .def(py::init(), + "rmin"_a, "rmax"_a, "halfz"_a, "halfphi"_a = M_PI, + "avgphi"_a = 0., "bevelMinZ"_a = 0., "bevelMaxZ"_a = 0.); + + py::enum_(cvb, "Face") + .value("PositiveDisc", CylinderVolumeBounds::Face::PositiveDisc) + .value("NegativeDisc", CylinderVolumeBounds::Face::NegativeDisc) + .value("OuterCylinder", CylinderVolumeBounds::Face::OuterCylinder) + .value("InnerCylinder", CylinderVolumeBounds::Face::InnerCylinder) + .value("NegativePhiPlane", CylinderVolumeBounds::Face::NegativePhiPlane) + .value("PositivePhiPlane", + CylinderVolumeBounds::Face::PositivePhiPlane); } { @@ -209,22 +233,72 @@ void addGeometry(Context& ctx) { })); } + py::class_(m, "ExtentEnvelope") + .def(py::init<>()) + .def(py::init()) + .def(py::init([](Envelope x, Envelope y, Envelope z, Envelope r, + Envelope phi, Envelope rPhi, Envelope h, Envelope eta, + Envelope mag) { + return ExtentEnvelope({.x = x, + .y = y, + .z = z, + .r = r, + .phi = phi, + .rPhi = rPhi, + .h = h, + .eta = eta, + .mag = mag}); + }), + py::arg("x") = zeroEnvelope, py::arg("y") = zeroEnvelope, + py::arg("z") = zeroEnvelope, py::arg("r") = zeroEnvelope, + py::arg("phi") = zeroEnvelope, py::arg("rPhi") = zeroEnvelope, + py::arg("h") = zeroEnvelope, py::arg("eta") = zeroEnvelope, + py::arg("mag") = zeroEnvelope) + .def_static("Zero", &ExtentEnvelope::Zero) + .def("__getitem__", [](ExtentEnvelope& self, + BinningValue bValue) { return self[bValue]; }) + .def("__setitem__", [](ExtentEnvelope& self, BinningValue bValue, + const Envelope& value) { self[bValue] = value; }) + .def("__str__", [](const ExtentEnvelope& self) { + std::array values; + + std::stringstream ss; + for (BinningValue val : allBinningValues()) { + ss << val << "=(" << self[val][0] << ", " << self[val][1] << ")"; + values.at(toUnderlying(val)) = ss.str(); + ss.str(""); + } + + ss.str(""); + ss << "ExtentEnvelope("; + ss << boost::algorithm::join(values, ", "); + ss << ")"; + return ss.str(); + }); + + py::class_(m, "Extent") + .def(py::init(), + py::arg("envelope") = ExtentEnvelope::Zero()) + .def("range", + [](const Acts::Extent& self, + Acts::BinningValue bval) -> std::array { + return {self.min(bval), self.max(bval)}; + }) + .def("__str__", &Extent::toString); + { - py::class_(m, "Extent") - .def(py::init( - [](const std::vector>>& - franges) { - Acts::Extent extent; - for (const auto& [bval, frange] : franges) { - extent.set(bval, frange[0], frange[1]); - } - return extent; - })) - .def("range", [](const Acts::Extent& self, Acts::BinningValue bval) { - return std::array{self.min(bval), - self.max(bval)}; - }); + auto cylStack = py::class_(m, "CylinderVolumeStack"); + + py::enum_(cylStack, + "AttachmentStrategy") + .value("Gap", CylinderVolumeStack::AttachmentStrategy::Gap) + .value("First", CylinderVolumeStack::AttachmentStrategy::First) + .value("Second", CylinderVolumeStack::AttachmentStrategy::Second) + .value("Midpoint", CylinderVolumeStack::AttachmentStrategy::Midpoint); + + py::enum_(cylStack, "ResizeStrategy") + .value("Gap", CylinderVolumeStack::ResizeStrategy::Gap) + .value("Expand", CylinderVolumeStack::ResizeStrategy::Expand); } } @@ -307,10 +381,15 @@ void addExperimentalGeometry(Context& ctx) { // Be able to construct a proto binning py::class_(m, "ProtoBinning") .def(py::init&, std::size_t>()) + const std::vector&, std::size_t>(), + "bValue"_a, "bType"_a, "e"_a, "exp"_a = 0u) .def(py::init()); + std::size_t>(), + "bValue"_a, "bType"_a, "minE"_a, "maxE"_a, "nbins"_a, "exp"_a = 0u) + .def(py::init(), + "bValue"_a, "bType"_a, "nbins"_a, "exp"_a = 0u); } { 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", [ diff --git a/Plugins/GeoModel/include/Acts/Plugins/GeoModel/GeoModelBlueprintCreater.hpp b/Plugins/GeoModel/include/Acts/Plugins/GeoModel/GeoModelBlueprintCreater.hpp index df895eb0741..c6358bdbec2 100644 --- a/Plugins/GeoModel/include/Acts/Plugins/GeoModel/GeoModelBlueprintCreater.hpp +++ b/Plugins/GeoModel/include/Acts/Plugins/GeoModel/GeoModelBlueprintCreater.hpp @@ -58,6 +58,8 @@ class GeoModelBlueprintCreater { std::string topEntry; /// Optionally override the top node bounds std::string topBoundsOverride = ""; + /// Export dot graph + std::string dotGraph = ""; }; /// The Blueprint return object diff --git a/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp b/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp index 93712d5d911..0b3d2f66e42 100644 --- a/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp +++ b/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp @@ -10,6 +10,7 @@ #include "Acts/Detector/GeometryIdGenerator.hpp" #include "Acts/Detector/LayerStructureBuilder.hpp" +#include "Acts/Detector/detail/BlueprintDrawer.hpp" #include "Acts/Detector/detail/BlueprintHelper.hpp" #include "Acts/Detector/interface/IGeometryIdGenerator.hpp" #include "Acts/Plugins/GeoModel/GeoModelTree.hpp" @@ -20,6 +21,8 @@ #include "Acts/Utilities/Helpers.hpp" #include "Acts/Utilities/RangeXD.hpp" +#include + #include using namespace Acts::detail; @@ -128,6 +131,14 @@ Acts::GeoModelBlueprintCreater::create(const GeometryContext& gctx, blueprint.topNode = createNode(cache, gctx, topEntry->second, blueprintTableMap, Extent()); + // Export to dot graph if configured + if (!options.dotGraph.empty()) { + std::ofstream dotFile(options.dotGraph); + Experimental::detail::BlueprintDrawer::dotStream(dotFile, + *blueprint.topNode); + dotFile.close(); + } + // Return the ready-to-use blueprint return blueprint; } diff --git a/Tests/UnitTests/Core/Geometry/CylinderVolumeStackTests.cpp b/Tests/UnitTests/Core/Geometry/CylinderVolumeStackTests.cpp index ab3c23a3dfc..19850be4bbe 100644 --- a/Tests/UnitTests/Core/Geometry/CylinderVolumeStackTests.cpp +++ b/Tests/UnitTests/Core/Geometry/CylinderVolumeStackTests.cpp @@ -813,6 +813,165 @@ 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); +} + +// original size +// <---------------> +// +---------------+ +// | | +// | | +// | Volume 1 | +// | | +// | | +// +---------------+ +// first resize +// <--------------------------> +// +---------------+----------+ +// | | | +// | | | +// | Volume 1 | Gap | +// | | | Gap is +// | | | reused!--+ +// +---------------+----------+ | +// second resize | +// <-----------------------------------> | +// +---------------+-------------------+ | +// | | | | +// | | | | +// | Volume 1 | Gap |<-----+ +// | | | +// | | | +// +---------------+-------------------+ +// +BOOST_AUTO_TEST_CASE(ResizeGapMultiple) { + Transform3 trf = Transform3::Identity(); + auto bounds = std::make_shared(70, 100, 100.0); + Volume vol{trf, bounds}; + + BOOST_TEST_CONTEXT("Positive") { + std::vector volumes = {&vol}; + CylinderVolumeStack stack(volumes, BinningValue::binZ, + CylinderVolumeStack::AttachmentStrategy::Gap, + CylinderVolumeStack::ResizeStrategy::Gap, + *logger); + + BOOST_CHECK_EQUAL(volumes.size(), 1); + BOOST_CHECK(stack.gaps().empty()); + + stack.update(std::make_shared(30.0, 100, 200), + trf * Translation3{Vector3::UnitZ() * 100}, *logger); + BOOST_CHECK_EQUAL(volumes.size(), 2); + BOOST_CHECK_EQUAL(stack.gaps().size(), 1); + + BOOST_CHECK_EQUAL(stack.gaps().front()->center()[eZ], 200.0); + const auto* cylBounds = dynamic_cast( + &stack.gaps().front()->volumeBounds()); + BOOST_REQUIRE_NE(cylBounds, nullptr); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eHalfLengthZ), + 100.0); + + stack.update(std::make_shared(30.0, 100, 300), + trf * Translation3{Vector3::UnitZ() * 200}, *logger); + + BOOST_CHECK_EQUAL(volumes.size(), 2); + // No additional gap volume was added! + BOOST_CHECK_EQUAL(stack.gaps().size(), 1); + + BOOST_CHECK_EQUAL(stack.gaps().front()->center()[eZ], 300.0); + cylBounds = dynamic_cast( + &stack.gaps().front()->volumeBounds()); + BOOST_REQUIRE_NE(cylBounds, nullptr); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eHalfLengthZ), + 200.0); + } + + BOOST_TEST_CONTEXT("Negative") { + std::vector volumes = {&vol}; + CylinderVolumeStack stack(volumes, BinningValue::binZ, + CylinderVolumeStack::AttachmentStrategy::Gap, + CylinderVolumeStack::ResizeStrategy::Gap, + *logger); + + BOOST_CHECK_EQUAL(volumes.size(), 1); + BOOST_CHECK(stack.gaps().empty()); + + stack.update(std::make_shared(30.0, 100, 200), + trf * Translation3{Vector3::UnitZ() * -100}, *logger); + BOOST_CHECK_EQUAL(volumes.size(), 2); + BOOST_CHECK_EQUAL(stack.gaps().size(), 1); + + BOOST_CHECK_EQUAL(stack.gaps().front()->center()[eZ], -200.0); + const auto* cylBounds = dynamic_cast( + &stack.gaps().front()->volumeBounds()); + BOOST_REQUIRE_NE(cylBounds, nullptr); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eHalfLengthZ), + 100.0); + + stack.update(std::make_shared(30.0, 100, 300), + trf * Translation3{Vector3::UnitZ() * -200}, *logger); + + BOOST_CHECK_EQUAL(volumes.size(), 2); + // No additional gap volume was added! + BOOST_CHECK_EQUAL(stack.gaps().size(), 1); + + BOOST_CHECK_EQUAL(stack.gaps().front()->center()[eZ], -300.0); + cylBounds = dynamic_cast( + &stack.gaps().front()->volumeBounds()); + BOOST_REQUIRE_NE(cylBounds, nullptr); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eHalfLengthZ), + 200.0); + } +} + BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE(RDirection) @@ -1503,6 +1662,82 @@ BOOST_DATA_TEST_CASE( } } +BOOST_AUTO_TEST_CASE(ResizeGapMultiple) { + Transform3 trf = Transform3::Identity(); + auto bounds = std::make_shared(100, 200, 100); + Volume vol{trf, bounds}; + + BOOST_TEST_CONTEXT("Outer") { + std::vector volumes = {&vol}; + CylinderVolumeStack stack(volumes, BinningValue::binR, + CylinderVolumeStack::AttachmentStrategy::Gap, + CylinderVolumeStack::ResizeStrategy::Gap, + *logger); + + BOOST_CHECK_EQUAL(volumes.size(), 1); + BOOST_CHECK(stack.gaps().empty()); + + stack.update(std::make_shared(100, 250, 100), trf, + *logger); + BOOST_CHECK_EQUAL(volumes.size(), 2); + BOOST_CHECK_EQUAL(stack.gaps().size(), 1); + + const auto* cylBounds = dynamic_cast( + &stack.gaps().front()->volumeBounds()); + BOOST_REQUIRE_NE(cylBounds, nullptr); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eMinR), 200); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eMaxR), 250); + + stack.update(std::make_shared(100, 300, 100), trf, + *logger); + + BOOST_CHECK_EQUAL(volumes.size(), 2); + // No additional gap volume was added! + BOOST_CHECK_EQUAL(stack.gaps().size(), 1); + + cylBounds = dynamic_cast( + &stack.gaps().front()->volumeBounds()); + BOOST_REQUIRE_NE(cylBounds, nullptr); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eMinR), 200); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eMaxR), 300); + } + + BOOST_TEST_CONTEXT("Inner") { + std::vector volumes = {&vol}; + CylinderVolumeStack stack(volumes, BinningValue::binR, + CylinderVolumeStack::AttachmentStrategy::Gap, + CylinderVolumeStack::ResizeStrategy::Gap, + *logger); + + BOOST_CHECK_EQUAL(volumes.size(), 1); + BOOST_CHECK(stack.gaps().empty()); + + stack.update(std::make_shared(50, 200, 100), trf, + *logger); + BOOST_CHECK_EQUAL(volumes.size(), 2); + BOOST_CHECK_EQUAL(stack.gaps().size(), 1); + + const auto* cylBounds = dynamic_cast( + &stack.gaps().front()->volumeBounds()); + BOOST_REQUIRE_NE(cylBounds, nullptr); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eMinR), 50); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eMaxR), 100); + + stack.update(std::make_shared(0, 200, 100), trf, + *logger); + + BOOST_CHECK_EQUAL(volumes.size(), 2); + // No additional gap volume was added! + BOOST_CHECK_EQUAL(stack.gaps().size(), 1); + + cylBounds = dynamic_cast( + &stack.gaps().front()->volumeBounds()); + BOOST_REQUIRE_NE(cylBounds, nullptr); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eMinR), 0); + BOOST_CHECK_EQUAL(cylBounds->get(CylinderVolumeBounds::eMaxR), 100); + } +} + BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE(Common) diff --git a/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp b/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp index dc032060a5b..e5345268e84 100644 --- a/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp +++ b/Tests/UnitTests/Core/Geometry/PortalLinkTests.cpp @@ -6,6 +6,7 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at https://mozilla.org/MPL/2.0/. +#include #include #include #include @@ -2135,7 +2136,161 @@ BOOST_AUTO_TEST_CASE(CompositeConstruction) { BOOST_AUTO_TEST_SUITE(PortalMerging) -BOOST_AUTO_TEST_CASE(TrivialTrivial) {} +BOOST_DATA_TEST_CASE(TrivialTrivial, + (boost::unit_test::data::make(0, -135, 180, 45) * + boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm}, + Vector3{20_mm, 0_mm, 0_mm}, + Vector3{0_mm, 20_mm, 0_mm}, + Vector3{20_mm, 20_mm, 0_mm}, + Vector3{0_mm, 0_mm, 20_mm})), + angle, offset) { + Transform3 base = + AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset); + + BOOST_TEST_CONTEXT("RDirection") { + auto vol1 = std::make_shared( + base, std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + base, std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol3 = std::make_shared( + base, std::make_shared(40_mm, 50_mm, 100_mm)); + + auto disc1 = + Surface::makeShared(Transform3::Identity(), 30_mm, 60_mm); + + auto disc2 = + Surface::makeShared(Transform3::Identity(), 60_mm, 90_mm); + + auto disc3 = + Surface::makeShared(Transform3::Identity(), 90_mm, 120_mm); + + auto trivial1 = std::make_unique(disc1, *vol1); + BOOST_REQUIRE(trivial1); + auto trivial2 = std::make_unique(disc2, *vol2); + BOOST_REQUIRE(trivial2); + auto trivial3 = std::make_unique(disc3, *vol3); + BOOST_REQUIRE(trivial3); + + auto grid1 = trivial1->makeGrid(BinningValue::binR); + auto compGridTrivial = PortalLinkBase::merge( + std::move(grid1), std::make_unique(*trivial2), + BinningValue::binR, *logger); + BOOST_REQUIRE(compGridTrivial); + BOOST_CHECK_EQUAL(dynamic_cast(*compGridTrivial) + .makeGrid(gctx, *logger), + nullptr); + + auto composite = PortalLinkBase::merge( + std::move(trivial1), std::move(trivial2), BinningValue::binR, *logger); + BOOST_REQUIRE(composite); + + auto grid12 = + dynamic_cast(*composite).makeGrid(gctx, *logger); + BOOST_REQUIRE(grid12); + + BOOST_CHECK_EQUAL( + grid12->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(), + vol1.get()); + + BOOST_CHECK_EQUAL( + grid12->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(), + vol2.get()); + + composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3), + BinningValue::binR, *logger); + BOOST_REQUIRE(composite); + + auto grid123 = + dynamic_cast(*composite).makeGrid(gctx, *logger); + BOOST_REQUIRE(grid123); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(), + vol1.get()); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(), + vol2.get()); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(), + vol3.get()); + } + + BOOST_TEST_CONTEXT("ZDirection") { + auto vol1 = std::make_shared( + base, std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + base * Translation3{Vector3::UnitZ() * 200}, + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol3 = std::make_shared( + base * Translation3{Vector3::UnitZ() * 400}, + std::make_shared(30_mm, 40_mm, 100_mm)); + + auto cyl1 = Surface::makeShared(base, 40_mm, 100_mm); + + auto cyl2 = Surface::makeShared( + base * Translation3{Vector3::UnitZ() * 200}, 40_mm, 100_mm); + + auto cyl3 = Surface::makeShared( + base * Translation3{Vector3::UnitZ() * 400}, 40_mm, 100_mm); + + auto trivial1 = std::make_unique(cyl1, *vol1); + BOOST_REQUIRE(trivial1); + auto trivial2 = std::make_unique(cyl2, *vol2); + BOOST_REQUIRE(trivial2); + auto trivial3 = std::make_unique(cyl3, *vol3); + BOOST_REQUIRE(trivial3); + + auto grid1 = trivial1->makeGrid(BinningValue::binZ); + auto compGridTrivial = PortalLinkBase::merge( + std::move(grid1), std::make_unique(*trivial2), + BinningValue::binZ, *logger); + BOOST_REQUIRE(compGridTrivial); + BOOST_CHECK_EQUAL(dynamic_cast(*compGridTrivial) + .makeGrid(gctx, *logger), + nullptr); + + auto composite = PortalLinkBase::merge( + std::move(trivial1), std::move(trivial2), BinningValue::binZ, *logger); + BOOST_REQUIRE(composite); + + auto grid12 = + dynamic_cast(*composite).makeGrid(gctx, *logger); + BOOST_REQUIRE(grid12); + + BOOST_CHECK_EQUAL( + grid12->resolveVolume(gctx, Vector2{40_mm, -40_mm}).value(), + vol1.get()); + + BOOST_CHECK_EQUAL( + grid12->resolveVolume(gctx, Vector2{40_mm, 40_mm}).value(), vol2.get()); + + composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3), + BinningValue::binZ, *logger); + BOOST_REQUIRE(composite); + + auto grid123 = + dynamic_cast(*composite).makeGrid(gctx, *logger); + BOOST_REQUIRE(grid123); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{40_mm, -110_mm}).value(), + vol1.get()); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{40_mm, -10_mm}).value(), + vol2.get()); + + BOOST_CHECK_EQUAL( + grid123->resolveVolume(gctx, Vector2{40_mm, 190_mm}).value(), + vol3.get()); + } +} BOOST_AUTO_TEST_CASE(TrivialGridR) { auto vol1 = std::make_shared( @@ -2167,32 +2322,14 @@ BOOST_AUTO_TEST_CASE(TrivialGridR) { auto merged = PortalLinkBase::merge(copy(trivial), copy(gridR), BinningValue::binR, *logger); BOOST_REQUIRE(merged); - - auto* mergedGrid = dynamic_cast(merged.get()); - BOOST_REQUIRE(mergedGrid); - - mergedGrid->printContents(std::cout); - - BOOST_CHECK_EQUAL(mergedGrid->grid().axes().size(), 1); - Axis axisExpected{AxisBound, {30_mm, 45_mm, 60_mm, 90_mm}}; - BOOST_CHECK_EQUAL(*mergedGrid->grid().axes().front(), axisExpected); + BOOST_CHECK_NE(dynamic_cast(merged.get()), nullptr); } BOOST_TEST_CONTEXT("Orthogonal") { auto merged = PortalLinkBase::merge(copy(gridPhi), copy(trivial), BinningValue::binR, *logger); BOOST_REQUIRE(merged); - - auto* mergedGrid = dynamic_cast(merged.get()); - BOOST_REQUIRE(mergedGrid); - - mergedGrid->printContents(std::cout); - - BOOST_CHECK_EQUAL(mergedGrid->grid().axes().size(), 2); - Axis axis1Expected{AxisBound, 30_mm, 90_mm, 2}; - Axis axis2Expected{AxisClosed, -M_PI, M_PI, 2}; - BOOST_CHECK_EQUAL(*mergedGrid->grid().axes().front(), axis1Expected); - BOOST_CHECK_EQUAL(*mergedGrid->grid().axes().back(), axis2Expected); + BOOST_CHECK_NE(dynamic_cast(merged.get()), nullptr); } } @@ -2227,38 +2364,14 @@ BOOST_AUTO_TEST_CASE(TrivialGridPhi) { auto merged = PortalLinkBase::merge(copy(trivial), copy(gridPhi), BinningValue::binPhi, *logger); BOOST_REQUIRE(merged); - - auto* mergedGrid = dynamic_cast(merged.get()); - BOOST_REQUIRE(mergedGrid); - - mergedGrid->printContents(std::cout); - - BOOST_CHECK_EQUAL(mergedGrid->grid().axes().size(), 1); - const auto& axis = *mergedGrid->grid().axes().front(); - BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable); - BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound); - std::vector expectedBins{-90_degree, 30_degree, 60_degree, - 90_degree}; - CHECK_CLOSE_REL(axis.getBinEdges(), expectedBins, 1e-6); + BOOST_CHECK_NE(dynamic_cast(merged.get()), nullptr); } BOOST_TEST_CONTEXT("Orthogonal") { auto merged = PortalLinkBase::merge(copy(gridR), copy(trivial), BinningValue::binPhi, *logger); BOOST_REQUIRE(merged); - - auto* mergedGrid = dynamic_cast(merged.get()); - BOOST_REQUIRE(mergedGrid); - - mergedGrid->printContents(std::cout); - - BOOST_CHECK_EQUAL(mergedGrid->grid().axes().size(), 2); - const auto& axis1 = *mergedGrid->grid().axes().front(); - const auto& axis2 = *mergedGrid->grid().axes().back(); - Axis axis1Expected{AxisBound, 30_mm, 100_mm, 2}; - BOOST_CHECK_EQUAL(axis1, axis1Expected); - std::vector expectedBins{-90_degree, 30_degree, 90_degree}; - CHECK_CLOSE_REL(axis2.getBinEdges(), expectedBins, 1e-6); + BOOST_CHECK_NE(dynamic_cast(merged.get()), nullptr); } } diff --git a/Tests/UnitTests/Core/Geometry/PortalShellTests.cpp b/Tests/UnitTests/Core/Geometry/PortalShellTests.cpp index 4c44617524a..15c5a9c7953 100644 --- a/Tests/UnitTests/Core/Geometry/PortalShellTests.cpp +++ b/Tests/UnitTests/Core/Geometry/PortalShellTests.cpp @@ -68,7 +68,7 @@ BOOST_AUTO_TEST_CASE(ConstructionFromVolume) { SingleCylinderPortalShell shell1{cyl1}; BOOST_CHECK_EQUAL(shell1.size(), 4); - using enum CylinderPortalShell::Face; + using enum CylinderVolumeBounds::Face; const auto* pDisc = shell1.portal(PositiveDisc); BOOST_REQUIRE_NE(pDisc, nullptr); @@ -299,7 +299,7 @@ BOOST_AUTO_TEST_CASE(ConstructionFromVolume) { // inner cylinder BOOST_AUTO_TEST_CASE(PortalAssignment) { - using enum CylinderPortalShell::Face; + using enum CylinderVolumeBounds::Face; TrackingVolume vol( Transform3::Identity(), std::make_shared(30_mm, 100_mm, 100_mm)); @@ -350,7 +350,7 @@ BOOST_AUTO_TEST_CASE(PortalAssignment) { BOOST_AUTO_TEST_SUITE(CylinderStack) BOOST_AUTO_TEST_CASE(ZDirection) { - using enum CylinderPortalShell::Face; + using enum CylinderVolumeBounds::Face; BOOST_TEST_CONTEXT("rMin>0") { TrackingVolume vol1( Transform3{Translation3{Vector3::UnitZ() * -100_mm}}, @@ -447,7 +447,7 @@ BOOST_AUTO_TEST_CASE(ZDirection) { } BOOST_AUTO_TEST_CASE(RDirection) { - using enum CylinderPortalShell::Face; + using enum CylinderVolumeBounds::Face; BOOST_TEST_CONTEXT("rMin>0") { TrackingVolume vol1( Transform3::Identity(), @@ -610,7 +610,7 @@ BOOST_AUTO_TEST_CASE(NestedStacks) { gctx, {&stack, &shell3}, BinningValue::binZ, *logger}; BOOST_CHECK(stack2.isValid()); - using enum CylinderPortalShell::Face; + using enum CylinderVolumeBounds::Face; auto lookup = [](auto& shell, CylinderPortalShell::Face face, Vector3 position, @@ -726,7 +726,7 @@ BOOST_AUTO_TEST_CASE(ConnectOuter) { SingleCylinderPortalShell shell{cyl1}; - using enum CylinderPortalShell::Face; + using enum CylinderVolumeBounds::Face; BOOST_CHECK_EQUAL( shell.portal(OuterCylinder)->getLink(Direction::AlongNormal), nullptr); BOOST_CHECK_EQUAL( @@ -749,7 +749,7 @@ BOOST_AUTO_TEST_CASE(ConnectOuter) { } BOOST_AUTO_TEST_CASE(RegisterInto) { - using enum CylinderPortalShell::Face; + using enum CylinderVolumeBounds::Face; TrackingVolume vol1( Transform3::Identity(), std::make_shared(0_mm, 100_mm, 100_mm)); diff --git a/Tests/UnitTests/Core/Geometry/PortalTests.cpp b/Tests/UnitTests/Core/Geometry/PortalTests.cpp index 49081748649..8997ccedda7 100644 --- a/Tests/UnitTests/Core/Geometry/PortalTests.cpp +++ b/Tests/UnitTests/Core/Geometry/PortalTests.cpp @@ -12,6 +12,7 @@ #include #include "Acts/Definitions/Units.hpp" +#include "Acts/Geometry/CompositePortalLink.hpp" #include "Acts/Geometry/CylinderVolumeBounds.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Geometry/GridPortalLink.hpp" @@ -171,10 +172,9 @@ BOOST_AUTO_TEST_CASE(Cylinder) { BOOST_CHECK_NE(merged12.getLink(Direction::AlongNormal), nullptr); BOOST_CHECK_EQUAL(merged12.getLink(Direction::OppositeNormal), nullptr); - auto grid12 = dynamic_cast( + auto composite12 = dynamic_cast( merged12.getLink(Direction::AlongNormal)); - BOOST_REQUIRE_NE(grid12, nullptr); - grid12->printContents(std::cout); + BOOST_REQUIRE_NE(composite12, nullptr); BOOST_CHECK_EQUAL( merged12 @@ -528,6 +528,67 @@ BOOST_AUTO_TEST_CASE(Material) { BOOST_CHECK_EQUAL(portal12.surface().surfaceMaterial(), material.get()); } +BOOST_AUTO_TEST_CASE(GridCreationOnFuse) { + Transform3 base = Transform3::Identity(); + + auto vol1 = std::make_shared( + base, std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol2 = std::make_shared( + base, std::make_shared(30_mm, 40_mm, 100_mm)); + + auto vol3 = std::make_shared( + base, std::make_shared(40_mm, 50_mm, 100_mm)); + + auto vol4 = std::make_shared( + base, std::make_shared(40_mm, 50_mm, 100_mm)); + + auto disc1 = + Surface::makeShared(Transform3::Identity(), 30_mm, 60_mm); + + auto disc2 = + Surface::makeShared(Transform3::Identity(), 60_mm, 90_mm); + + auto disc3 = + Surface::makeShared(Transform3::Identity(), 90_mm, 120_mm); + + auto trivial1 = std::make_unique(disc1, *vol1); + BOOST_REQUIRE(trivial1); + auto trivial2 = std::make_unique(disc2, *vol2); + BOOST_REQUIRE(trivial2); + auto trivial3 = std::make_unique(disc3, *vol3); + BOOST_REQUIRE(trivial3); + + std::vector> links; + links.push_back(std::move(trivial1)); + links.push_back(std::move(trivial2)); + links.push_back(std::move(trivial3)); + + auto composite = std::make_unique(std::move(links), + BinningValue::binR); + + auto discOpposite = + Surface::makeShared(Transform3::Identity(), 30_mm, 120_mm); + + auto trivialOpposite = + std::make_unique(discOpposite, *vol4); + + Portal aPortal{gctx, std::move(composite), nullptr}; + Portal bPortal{gctx, nullptr, std::move(trivialOpposite)}; + + Portal fused = Portal::fuse(gctx, aPortal, bPortal, *logger); + + BOOST_CHECK_NE(dynamic_cast( + fused.getLink(Direction::OppositeNormal)), + nullptr); + + const auto* grid = dynamic_cast( + fused.getLink(Direction::AlongNormal)); + BOOST_REQUIRE_NE(grid, nullptr); + + BOOST_CHECK_EQUAL(grid->grid().axes().front()->getNBins(), 3); +} + BOOST_AUTO_TEST_SUITE_END() // Fusing BOOST_AUTO_TEST_CASE(Construction) {