diff --git a/CI/check_license.py b/CI/check_license.py index 46a6087a922..73bc343f763 100755 --- a/CI/check_license.py +++ b/CI/check_license.py @@ -189,7 +189,7 @@ def year_print(*pargs): exit = 0 srcs = list(srcs) nsrcs = len(srcs) - step = int(nsrcs / 20) + step = max(int(nsrcs / 20), 1) for i, src in enumerate(srcs): if any([fnmatch(src, e) for e in args.exclude]): continue diff --git a/Core/include/Acts/Definitions/Algebra.hpp b/Core/include/Acts/Definitions/Algebra.hpp index fa5edcc2252..c489f5ecd94 100644 --- a/Core/include/Acts/Definitions/Algebra.hpp +++ b/Core/include/Acts/Definitions/Algebra.hpp @@ -100,4 +100,6 @@ using AngleAxis3 = Eigen::AngleAxis; using Transform2 = Eigen::Transform; using Transform3 = Eigen::Transform; +constexpr ActsScalar s_transformEquivalentTolerance = 1e-9; + } // namespace Acts diff --git a/Core/include/Acts/Geometry/Portal.hpp b/Core/include/Acts/Geometry/Portal.hpp new file mode 100644 index 00000000000..c6ec526a42d --- /dev/null +++ b/Core/include/Acts/Geometry/Portal.hpp @@ -0,0 +1,234 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Direction.hpp" +#include "Acts/Utilities/BinningType.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "Acts/Utilities/Result.hpp" + +#include + +namespace Acts { + +class RegularSurface; +class GeometryContext; +class TrackingVolume; +class CylinderSurface; +class PlaneSurface; +class DiscSurface; +class Surface; + +class PortalLinkBase; + +/// Exception thrown when portals cannot be merged +class PortalMergingException : public std::exception { + const char* what() const noexcept override; +}; + +/// Exception thrown when portals cannot be fused +class PortalFusingException : public std::exception { + const char* what() const noexcept override; +}; + +/// A portal connects two or more neighboring volumes. Each volume has a set of +/// portals that describes which volumes lie behind the portal in that +/// direction. Portals use associated portal links to perform lookups of target +/// volumes. +/// Each portal has two links, and a corresponding surface. One link is +/// associated with the direction along the surface's normal vector, and one +/// with the opposite direction. +class Portal { + public: + /// Constructor for a portal from a single link + /// @param direction The direction of the link + /// @param link The portal link + Portal(Direction direction, std::unique_ptr link); + + /// Constructor for a portal from a surface and volume, where a trivial portal + /// link is automatically constructed. + /// @param direction The direction of the link + /// @param surface The surface from which to create the portal link + /// @param volume The volume this portal connects to in the @p direction + /// relative to the normal of @p surface. + Portal(Direction direction, std::shared_ptr surface, + TrackingVolume& volume); + + /// Constructor for a portal from two links. One of the links can be + /// `nullptr`, but at least one of them needs to be set. If both are set, they + /// need to be valid compatible links that can be fused. + /// @param gctx The geometry context + /// @param alongNormal The link along the normal of the surface + /// @param oppositeNormal The link opposite to the normal of the + Portal(const GeometryContext& gctx, + std::unique_ptr alongNormal, + std::unique_ptr oppositeNormal); + + /// Helper struct for the arguments to the portal constructor below using + /// designated initializers. + struct Arguments { + /// Aggregate over a surface and a volume with optional semantics + struct Link { + Link() = default; + /// Constructor from a surface and a volume + Link(std::shared_ptr surfaceIn, TrackingVolume& volumeIn) + : surface(std::move(surfaceIn)), volume(&volumeIn) {} + + /// The associated surface + std::shared_ptr surface = nullptr; + /// The associated volume + TrackingVolume* volume = nullptr; + }; + + /// Entry for the link along normal + /// Entry for the link opposite normal + Link alongNormal{}; + Link oppositeNormal{}; + }; + + /// Constructor that takes a geometry context and an rvalue reference to a + /// helper struct from above. This pattern allows you to use designated + /// initializers to construct this object like: + /// ```cpp + /// Portal{gctx, {.oppositeNormal = {cyl1, *vol1}}}; + /// Portal{gctx, {.alongNormal = {cyl2, *vol2}}}; + /// ``` + /// @param gctx The geometry context + /// @param args The struct containing the arguments + Portal(const GeometryContext& gctx, Arguments&& args); + + /// Fuse two portals together. Fusing is the combination of two portal links + /// on the same logical surfaces. The actual surface instances can be + /// different, as long as they are geometrically equivalent (within numerical + /// precision). The resulting portal will have one portal along the shared + /// surface's normal vector, and one opposite that vector. + /// + /// portal1 portal2 + /// +---+ +---+ + /// | | | | + /// | | | | + /// <----+ | + | +----> + /// | | | | + /// | | | | + /// +---+ +---+ + /// + /// @note The input portals need to have compatible link loadaout, e.g. one + /// portal needs to have the *along normal* slot filled, and the + /// otherone one needs to have the *opposite normal* slot filled. If + /// portals share a filled slot, the function throws an exception. + /// @note This is a destructive operation on the portals involved + /// @param gctx The geometry context + /// @param aPortal The first portal + /// @param bPortal The second portal + /// @param logger The logger to push output to + static Portal fuse(const GeometryContext& gctx, Portal& aPortal, + Portal& bPortal, const Logger& logger = getDummyLogger()); + + /// Merge two adjacent portals with each other to produce a new portal that + /// encompasses both inputs. It is the complementary operation to the fusing + /// of portals. To be able to merge portals, the surfaces of their associated + /// links need to be *mergeable*, and the portal links need to be compatible. + /// This means that both portals need to have a link along the portal surface + /// normal, opposite the normal, or both. If the equipped links are opposite + /// relative to one another (e.g. one along one opposite), the function will + /// throw an exception. + /// + /// ^ ^ + /// | | + /// portal1| portal2| + /// +-------+-------+ +-------+-------+ + /// | | + | | + /// +-------+-------+ +-------+-------+ + /// | | + /// | | + /// v v + /// + /// @note This is a destructive operation on both portals, their + /// links will be moved to produce merged links, which can fail + /// if the portal links are not compatible + /// @param gctx The geometry context + /// @param aPortal The first portal + /// @param bPortal The second portal + /// @param direction The direction of the merge (e.g. along z) + /// @param logger The logger to push output to + static Portal merge(const GeometryContext& gctx, Portal& aPortal, + Portal& bPortal, BinningValue direction, + const Logger& logger = getDummyLogger()); + + /// Resolve the volume for a 3D position and a direction + /// The @p direction is used to select the right portal link, if it is set. + /// In case no link is found in the specified direction, a `nullptr` is + /// returned. + /// @param gctx The geometry context + /// @param position The 3D position + /// @param direction The direction + /// @return The target volume (can be `nullptr`) + Result resolveVolume(const GeometryContext& gctx, + const Vector3& position, + const Vector3& direction) const; + + /// Set a link on the portal into the slot specified by the direction. + /// @note The surface associated with @p link must be logically equivalent + /// to the one of the link that's already set on the portal. + /// @param gctx The geometry context + /// @param direction The direction + /// @param link The link to set + void setLink(const GeometryContext& gctx, Direction direction, + std::unique_ptr link); + + /// Helper function create a trivial portal link based on a surface. + /// @param gctx The geometry context + /// @param direction The direction of the link to create + /// @param surface The surface + /// @note The @p surface must be logically equivalent + /// to the one of the link that's already set on the portal. + /// @param volume The target volume + void setLink(const GeometryContext& gctx, Direction direction, + std::shared_ptr surface, TrackingVolume& volume); + + /// Get the link associated with the @p direction. Can be null if the associated link is unset. + /// @param direction The direction + /// @return The link (can be null) + const PortalLinkBase* getLink(Direction direction) const; + + /// Returns true if the portal is valid, that means it has at least one + /// non-null link associated.Portals can be in an invalid state after they get + /// merged or fused with other portals. + /// @return True if the portal is valid + bool isValid() const; + + /// Create and attach a trivial portal link to the empty slot of this portal + /// @param volume The target volume to connect to + void fill(TrackingVolume& volume); + + /// Access the portal surface that is shared between the two links + /// @return The portal surface + const RegularSurface& surface() const; + + private: + /// Helper to check surface equivalence without checking material status. This + /// is needed because we allow fusing portals with surfaces that are + /// equivalent but one of them has material while the other does not. The + /// normal surface comparison would determine these surfaces as not + /// equivalent. + /// @param gctx The geometry context + /// @param a The first surface + /// @param b The second surface + /// @return True if the surfaces are equivalent + static bool isSameSurface(const GeometryContext& gctx, const Surface& a, + const Surface& b); + + std::shared_ptr m_surface; + + std::unique_ptr m_alongNormal; + std::unique_ptr m_oppositeNormal; +}; + +} // namespace Acts diff --git a/Core/include/Acts/Propagator/Navigator.hpp b/Core/include/Acts/Propagator/Navigator.hpp index 007e0b7d19a..99f6a436341 100644 --- a/Core/include/Acts/Propagator/Navigator.hpp +++ b/Core/include/Acts/Propagator/Navigator.hpp @@ -168,10 +168,6 @@ class Navigator { const Layer* currentLayer = nullptr; /// Navigation state - external state: the current surface const Surface* currentSurface = nullptr; - /// Navigation state: the target volume - const TrackingVolume* targetVolume = nullptr; - /// Navigation state: the target layer - const Layer* targetLayer = nullptr; /// Navigation state: the target surface const Surface* targetSurface = nullptr; @@ -348,16 +344,7 @@ class Navigator { if (state.navigation.startLayer != nullptr) { ACTS_VERBOSE(volInfo(state) << "Start layer to be resolved."); // We provide the layer to the resolve surface method in this case - bool startResolved = resolveSurfaces(state, stepper); - if (!startResolved && - state.navigation.startLayer == state.navigation.targetLayer) { - ACTS_VERBOSE(volInfo(state) - << "Start is target layer and we have no surface " - "candidates. Nothing left to do."); - // set the navigation break - state.navigation.navigationBreak = true; - stepper.releaseStepSize(state.stepping, ConstrainedStep::actor); - } + resolveSurfaces(state, stepper); } } @@ -386,13 +373,6 @@ class Navigator { // Navigator pre step always resets the current surface state.navigation.currentSurface = nullptr; - // Initialize the target and target volume - if (state.navigation.targetSurface != nullptr && - state.navigation.targetVolume == nullptr) { - // Find out about the target as much as you can - initializeTarget(state, stepper); - } - auto tryTargetNextSurface = [&]() { // Try targeting the surfaces - then layers - then boundaries @@ -576,19 +556,6 @@ class Navigator { state.navigation.navigationStage = Stage::boundaryTarget; ACTS_VERBOSE(volInfo(state) << "Staying focussed on boundary."); } - } else if (state.navigation.currentVolume == - state.navigation.targetVolume) { - if (state.navigation.targetSurface == nullptr) { - ACTS_WARNING(volInfo(state) - << "No further navigation action, proceed to " - "target. This is very likely an error"); - } else { - ACTS_VERBOSE(volInfo(state) - << "No further navigation action, proceed to target."); - } - // Set navigation break and release the navigation step size - state.navigation.navigationBreak = true; - stepper.releaseStepSize(state.stepping, ConstrainedStep::actor); } else { ACTS_VERBOSE(volInfo(state) << "Status could not be determined - good luck."); @@ -752,12 +719,6 @@ class Navigator { } else { ACTS_VERBOSE(volInfo(state) << "Last surface on layer reached, and no layer."); - if (state.navigation.currentVolume == state.navigation.targetVolume) { - ACTS_VERBOSE(volInfo(state) - << "This is the target volume, stop navigation."); - state.navigation.navigationBreak = true; - stepper.releaseStepSize(state.stepping, ConstrainedStep::actor); - } } } @@ -835,30 +796,8 @@ class Navigator { ++state.navigation.navLayerIndex; } - // Re-initialize target at last layer, only in case it is the target volume - // This avoids a wrong target volume estimation - if (state.navigation.currentVolume == state.navigation.targetVolume) { - initializeTarget(state, stepper); - } - // Screen output - if (logger().doPrint(Logging::VERBOSE)) { - std::ostringstream os; - os << "Last layer"; - if (state.navigation.currentVolume == state.navigation.targetVolume) { - os << " (final volume) done, proceed to target."; - } else { - os << " done, target volume boundary."; - } - logger().log(Logging::VERBOSE, os.str()); - } + ACTS_VERBOSE(volInfo(state) << "Last layer done, target volume boundary."); - // Set the navigation break if necessary - if (state.navigation.currentVolume == state.navigation.targetVolume) { - ACTS_VERBOSE(volInfo(state) - << "This is the target volume, stop navigation."); - state.navigation.navigationBreak = true; - stepper.releaseStepSize(state.stepping, ConstrainedStep::actor); - } return false; } @@ -900,19 +839,8 @@ class Navigator { "stopping navigation."); stepper.releaseStepSize(state.stepping, ConstrainedStep::actor); return false; - } else if (state.navigation.currentVolume == - state.navigation.targetVolume) { - ACTS_VERBOSE(volInfo(state) - << "In target volume: no need to resolve boundary, " - "stopping navigation."); - state.navigation.navigationBreak = true; - stepper.releaseStepSize(state.stepping, ConstrainedStep::actor); - return true; } - // Re-initialize target at volume boundary - initializeTarget(state, stepper); - // Helper function to find boundaries auto findBoundaries = [&]() -> bool { // The navigation options @@ -1005,85 +933,6 @@ class Navigator { return false; } - /// @brief Navigation (re-)initialisation for the target - /// - /// @note This is only called a few times every propagation/extrapolation - /// - /// As a straight line estimate can lead you to the wrong destination - /// Volume, this will be called at: - /// - initialization - /// - attempted volume switch - /// Target finding by association will not be done again - /// - /// @tparam propagator_state_t The state type of the propagator - /// @tparam stepper_t The type of stepper used for the propagation - /// - /// @param [in,out] state is the propagation state object - /// @param [in] stepper Stepper in use - /// - /// boolean return triggers exit to stepper - template - void initializeTarget(propagator_state_t& state, - const stepper_t& stepper) const { - if (state.navigation.targetVolume != nullptr && - state.stepping.pathAccumulated == 0.) { - ACTS_VERBOSE(volInfo(state) - << "Re-initialzing cancelled as it is the first step."); - return; - } - // Fast Navigation initialization for target: - if (state.navigation.targetSurface != nullptr && - state.navigation.targetSurface->associatedLayer() && - state.navigation.targetVolume == nullptr) { - ACTS_VERBOSE(volInfo(state) - << "Fast target initialization through association."); - ACTS_VERBOSE(volInfo(state) - << "Target surface set to " - << state.navigation.targetSurface->geometryId()); - state.navigation.targetLayer = - state.navigation.targetSurface->associatedLayer(); - state.navigation.targetVolume = - state.navigation.targetLayer->trackingVolume(); - } else if (state.navigation.targetSurface != nullptr) { - // screen output that you do a re-initialization - if (state.navigation.targetVolume != nullptr) { - ACTS_VERBOSE(volInfo(state) - << "Re-initialization of target volume triggered."); - } - // Slow navigation initialization for target: - // target volume and layer search through global search - auto targetIntersection = - state.navigation.targetSurface - ->intersect( - state.geoContext, stepper.position(state.stepping), - state.options.direction * stepper.direction(state.stepping), - BoundaryTolerance::Infinite(), state.options.surfaceTolerance) - .closest(); - if (targetIntersection.isValid()) { - ACTS_VERBOSE(volInfo(state) - << "Target estimate position (" - << targetIntersection.position().x() << ", " - << targetIntersection.position().y() << ", " - << targetIntersection.position().z() << ")"); - /// get the target volume from the intersection - auto tPosition = targetIntersection.position(); - state.navigation.targetVolume = - m_cfg.trackingGeometry->lowestTrackingVolume(state.geoContext, - tPosition); - state.navigation.targetLayer = - state.navigation.targetVolume != nullptr - ? state.navigation.targetVolume->associatedLayer( - state.geoContext, tPosition) - : nullptr; - if (state.navigation.targetVolume != nullptr) { - ACTS_VERBOSE(volInfo(state) - << "Target volume estimated : " - << state.navigation.targetVolume->volumeName()); - } - } - } - } - /// @brief Resolve the surfaces of this layer /// /// @tparam propagator_state_t The state type of the propagator diff --git a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp index 7156160a606..34dabcd343d 100644 --- a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp +++ b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp @@ -982,6 +982,8 @@ class CombinatorialKalmanFilter { typeFlags.set(TrackStateFlag::MeasurementFlag); // Increment number of measurements newBranch.nMeasurements()++; + newBranch.nDoF() += trackState.calibratedSize(); + newBranch.chi2() += trackState.chi2(); } else { ACTS_WARNING("Cannot handle this track state flags"); continue; @@ -1330,10 +1332,6 @@ class CombinatorialKalmanFilter { return error.error(); } - for (const auto& track : combKalmanResult.collectedTracks) { - calculateTrackQuantities(track); - } - return std::move(combKalmanResult.collectedTracks); } }; diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index deed734ac92..829cd022f60 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -228,7 +228,8 @@ struct Gx2FitterResult { const TrackingVolume* startVolume = nullptr; }; -/// addToGx2fSums Function +/// @brief Process measurements and fill the aMatrix and bVector +/// /// The function processes each measurement for the GX2F Actor fitting process. /// It extracts the information from the track state and adds it to aMatrix, /// bVector, and chi2sum. @@ -246,64 +247,72 @@ template void addToGx2fSums(BoundMatrix& aMatrix, BoundVector& bVector, double& chi2sum, const BoundMatrix& jacobianFromStart, const track_state_t& trackState, const Logger& logger) { - BoundVector predicted = trackState.predicted(); + // First we get back the covariance and try to invert it. If the inversion + // fails, we can already abort. + const ActsSquareMatrix covarianceMeasurement = + trackState.template calibratedCovariance(); - ActsVector measurement = trackState.template calibrated(); + const auto safeInvCovMeasurement = safeInverse(covarianceMeasurement); + if (!safeInvCovMeasurement) { + ACTS_WARNING("addToGx2fSums: safeInvCovMeasurement failed."); + ACTS_VERBOSE(" covarianceMeasurement:\n" << covarianceMeasurement); + return; + } - ActsSquareMatrix covarianceMeasurement = - trackState.template calibratedCovariance(); + const BoundVector predicted = trackState.predicted(); - ActsMatrix projector = + const ActsVector measurement = + trackState.template calibrated(); + + const ActsMatrix projector = trackState.projector().template topLeftCorner(); - ActsMatrix projJacobian = projector * jacobianFromStart; - - ActsMatrix projPredicted = projector * predicted; - - ActsVector residual = measurement - projPredicted; - - ACTS_VERBOSE("Contributions in addToGx2fSums:\n" - << "kMeasDim: " << kMeasDim << "\n" - << "predicted" << predicted.transpose() << "\n" - << "measurement: " << measurement.transpose() << "\n" - << "covarianceMeasurement:\n" - << covarianceMeasurement << "\n" - << "projector:\n" - << projector.eval() << "\n" - << "projJacobian:\n" - << projJacobian.eval() << "\n" - << "projPredicted: " << (projPredicted.transpose()).eval() - << "\n" - << "residual: " << (residual.transpose()).eval()); - - auto safeInvCovMeasurement = safeInverse(covarianceMeasurement); - - if (safeInvCovMeasurement) { - chi2sum += - (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); - aMatrix += - (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) - .eval(); - bVector += - (residual.transpose() * (*safeInvCovMeasurement) * projJacobian).eval(); - - ACTS_VERBOSE( - "aMatrixMeas:\n" - << (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) - .eval() - << "\n" - << "bVectorMeas: " - << (residual.transpose() * (*safeInvCovMeasurement) * projJacobian) - .eval() - << "\n" - << "chi2sumMeas: " - << (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0) - << "\n" - << "safeInvCovMeasurement:\n" - << (*safeInvCovMeasurement)); - } else { - ACTS_WARNING("safeInvCovMeasurement failed"); - } + const ActsMatrix projJacobian = + projector * jacobianFromStart; + + const ActsMatrix projPredicted = projector * predicted; + + const ActsVector residual = measurement - projPredicted; + + // Finally contribute to chi2sum, aMatrix, and bVector + chi2sum += (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0); + + aMatrix += + (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) + .eval(); + + bVector += (residual.transpose() * (*safeInvCovMeasurement) * projJacobian) + .eval() + .transpose(); + + ACTS_VERBOSE( + "Contributions in addToGx2fSums:\n" + << "kMeasDim: " << kMeasDim << "\n" + << "predicted" << predicted.transpose() << "\n" + << "measurement: " << measurement.transpose() << "\n" + << "covarianceMeasurement:\n" + << covarianceMeasurement << "\n" + << "projector:\n" + << projector.eval() << "\n" + << "projJacobian:\n" + << projJacobian.eval() << "\n" + << "projPredicted: " << (projPredicted.transpose()).eval() << "\n" + << "residual: " << (residual.transpose()).eval() << "\n" + << "jacobianFromStart:\n" + << jacobianFromStart << "aMatrixMeas:\n" + << (projJacobian.transpose() * (*safeInvCovMeasurement) * projJacobian) + .eval() + << "\n" + << "bVectorMeas: " + << (residual.transpose() * (*safeInvCovMeasurement) * projJacobian).eval() + << "\n" + << "chi2sumMeas: " + << (residual.transpose() * (*safeInvCovMeasurement) * residual)(0, 0) + << "\n" + << "safeInvCovMeasurement:\n" + << (*safeInvCovMeasurement)); + + return; } /// calculateDeltaParams Function @@ -361,8 +370,8 @@ class Gx2Fitter { /// @tparam calibrator_t The type of calibrator /// @tparam outlier_finder_t Type of the outlier finder class /// - /// The GX2FnActor does not rely on the measurements to be - /// sorted along the track. /// TODO is this true? + /// The GX2F tor does not rely on the measurements to be sorted along the + /// track. template class Actor { public: @@ -455,179 +464,190 @@ class Gx2Fitter { } result.startVolume = state.navigation.startVolume; - // Update: - // - Waiting for a current surface + // We are only interested in surfaces. If we are not on a surface, we + // continue the navigation auto surface = navigator.currentSurface(state.navigation); - if (surface != nullptr) { - ++result.surfaceCount; - const GeometryIdentifier geoId = surface->geometryId(); - ACTS_DEBUG("Surface " << geoId << " detected."); - const bool surfaceIsSensitive = - (surface->associatedDetectorElement() != nullptr); - const bool surfaceHasMaterial = (surface->surfaceMaterial() != nullptr); - const bool doMaterial = multipleScattering && surfaceHasMaterial; - - // Found material - if (doMaterial) { - ACTS_DEBUG(" The surface contains material, ..."); - } + if (surface == nullptr) { + return; + } - // Check if we have a measurement surface - if (auto sourcelink_it = inputMeasurements->find(geoId); - sourcelink_it != inputMeasurements->end()) { - ACTS_DEBUG(" The surface contains a measurement."); - - // Transport the covariance to the surface - stepper.transportCovarianceToBound(state.stepping, *surface, - freeToBoundCorrection); - - // TODO generalize the update of the currentTrackIndex - auto& fittedStates = *result.fittedStates; - - // Add a TrackState entry multi trajectory. This - // allocates storage for all components, which we will set later. - typename traj_t::TrackStateProxy trackStateProxy = - fittedStates.makeTrackState(Gx2fConstants::trackStateMask, - result.lastTrackIndex); - const std::size_t currentTrackIndex = trackStateProxy.index(); - - // Set the trackStateProxy components with the state from the ongoing - // propagation - { - trackStateProxy.setReferenceSurface(surface->getSharedPtr()); - // Bind the transported state to the current surface - auto res = stepper.boundState(state.stepping, *surface, false, - freeToBoundCorrection); - if (!res.ok()) { - result.result = res.error(); - return; - } - const auto& [boundParams, jacobian, pathLength] = *res; - - // Fill the track state - trackStateProxy.predicted() = boundParams.parameters(); - trackStateProxy.predictedCovariance() = state.stepping.cov; - - trackStateProxy.jacobian() = jacobian; - trackStateProxy.pathLength() = pathLength; + ++result.surfaceCount; + const GeometryIdentifier geoId = surface->geometryId(); + ACTS_DEBUG("Surface " << geoId << " detected."); + const bool surfaceIsSensitive = + (surface->associatedDetectorElement() != nullptr); + const bool surfaceHasMaterial = (surface->surfaceMaterial() != nullptr); + const bool doMaterial = multipleScattering && surfaceHasMaterial; + + // Found material + if (doMaterial) { + ACTS_DEBUG(" The surface contains material, ..."); + } + + // Here we handle all measurements + if (auto sourcelink_it = inputMeasurements->find(geoId); + sourcelink_it != inputMeasurements->end()) { + ACTS_DEBUG(" The surface contains a measurement."); + + // Transport the covariance to the surface + stepper.transportCovarianceToBound(state.stepping, *surface, + freeToBoundCorrection); + + // TODO generalize the update of the currentTrackIndex + auto& fittedStates = *result.fittedStates; + + // Add a TrackState entry multi trajectory. This + // allocates storage for all components, which we will set later. + typename traj_t::TrackStateProxy trackStateProxy = + fittedStates.makeTrackState(Gx2fConstants::trackStateMask, + result.lastTrackIndex); + const std::size_t currentTrackIndex = trackStateProxy.index(); + + // Set the trackStateProxy components with the state from the ongoing + // propagation + { + trackStateProxy.setReferenceSurface(surface->getSharedPtr()); + // Bind the transported state to the current surface + auto res = stepper.boundState(state.stepping, *surface, false, + freeToBoundCorrection); + if (!res.ok()) { + result.result = res.error(); + return; } + const auto& [boundParams, jacobian, pathLength] = *res; - // We have predicted parameters, so calibrate the uncalibrated input - // measurement - extensions.calibrator(state.geoContext, *calibrationContext, - sourcelink_it->second, trackStateProxy); + // Fill the track state + trackStateProxy.predicted() = boundParams.parameters(); + trackStateProxy.predictedCovariance() = state.stepping.cov; - // Get and set the type flags - auto typeFlags = trackStateProxy.typeFlags(); - typeFlags.set(TrackStateFlag::ParameterFlag); - if (surfaceHasMaterial) { - typeFlags.set(TrackStateFlag::MaterialFlag); - } + trackStateProxy.jacobian() = jacobian; + trackStateProxy.pathLength() = pathLength; + } + + // We have predicted parameters, so calibrate the uncalibrated input + // measurement + extensions.calibrator(state.geoContext, *calibrationContext, + sourcelink_it->second, trackStateProxy); + + // Get and set the type flags + auto typeFlags = trackStateProxy.typeFlags(); + typeFlags.set(TrackStateFlag::ParameterFlag); + if (surfaceHasMaterial) { + typeFlags.set(TrackStateFlag::MaterialFlag); + } + + // Set the measurement type flag + typeFlags.set(TrackStateFlag::MeasurementFlag); + // We count the processed measurement + ++result.processedMeasurements; + + result.lastMeasurementIndex = currentTrackIndex; + result.lastTrackIndex = currentTrackIndex; + + // TODO check for outlier first + // We count the state with measurement + ++result.measurementStates; + + // We count the processed state + ++result.processedStates; + + // Update the number of holes count only when encountering a + // measurement + result.measurementHoles = result.missedActiveSurfaces.size(); + + return; + } + + // Here we handle material for multipleScattering. If holes exist, we also + // handle them already. We create a full trackstate (unlike for simple + // holes), since we need to evaluate the material later + if (doMaterial) { + // TODO add material handling + ACTS_DEBUG( + " The surface contains no measurement, but material and maybe " + "a hole."); - // Set the measurement type flag - typeFlags.set(TrackStateFlag::MeasurementFlag); - // We count the processed measurement - ++result.processedMeasurements; - - result.lastMeasurementIndex = currentTrackIndex; - result.lastTrackIndex = currentTrackIndex; - - // TODO check for outlier first - // We count the state with measurement - ++result.measurementStates; - - // We count the processed state - ++result.processedStates; - - // Update the number of holes count only when encountering a - // measurement - result.measurementHoles = result.missedActiveSurfaces.size(); - } else if (doMaterial) { - // Here we handle material for multipleScattering. If holes exist, we - // also handle them already. We create a full trackstate (unlike for - // simple holes), since we need to evaluate the material later - // TODO add material handling + return; + } + + // Here we handle holes. If material hasn't been handled before (because + // multipleScattering is turned off), we will also handle it here + if (surfaceIsSensitive || surfaceHasMaterial) { + if (multipleScattering) { ACTS_DEBUG( - " The surface contains no measurement, but material and maybe " - "a hole."); - } else if (surfaceIsSensitive || surfaceHasMaterial) { - // Here we handle holes. If material hasn't been handled before - // (because multipleScattering is turned off), we will also handle it - // here - if (multipleScattering) { - ACTS_DEBUG( - " The surface contains no measurement, but maybe a hole."); - } else { - ACTS_DEBUG( - " The surface contains no measurement, but maybe a hole " - "and/or material."); - } + " The surface contains no measurement, but maybe a hole."); + } else { + ACTS_DEBUG( + " The surface contains no measurement, but maybe a hole " + "and/or material."); + } - // We only create track states here if there is already a measurement - // detected (no holes before the first measurement) or if we encounter - // material - const bool precedingMeasurementExists = - (result.measurementStates > 0); - if (!precedingMeasurementExists && !surfaceHasMaterial) { - ACTS_DEBUG( - " Ignoring hole, because there are no preceding " - "measurements."); - return; - } + // We only create track states here if there is already a measurement + // detected (no holes before the first measurement) or if we encounter + // material + const bool precedingMeasurementExists = (result.measurementStates > 0); + if (!precedingMeasurementExists && !surfaceHasMaterial) { + ACTS_DEBUG( + " Ignoring hole, because there are no preceding " + "measurements."); + return; + } - auto& fittedStates = *result.fittedStates; - - // Add a TrackState entry multi trajectory. This - // allocates storage for all components, which we will set later. - typename traj_t::TrackStateProxy trackStateProxy = - fittedStates.makeTrackState(Gx2fConstants::trackStateMask, - result.lastTrackIndex); - const std::size_t currentTrackIndex = trackStateProxy.index(); - - // Set the trackStateProxy components with the state from the - // ongoing propagation - { - trackStateProxy.setReferenceSurface(surface->getSharedPtr()); - // Bind the transported state to the current surface - auto res = stepper.boundState(state.stepping, *surface, false, - freeToBoundCorrection); - if (!res.ok()) { - result.result = res.error(); - return; - } - const auto& [boundParams, jacobian, pathLength] = *res; - - // Fill the track state - trackStateProxy.predicted() = boundParams.parameters(); - trackStateProxy.predictedCovariance() = state.stepping.cov; - - trackStateProxy.jacobian() = jacobian; - trackStateProxy.pathLength() = pathLength; + auto& fittedStates = *result.fittedStates; + + // Add a TrackState entry multi trajectory. This + // allocates storage for all components, which we will set later. + typename traj_t::TrackStateProxy trackStateProxy = + fittedStates.makeTrackState(Gx2fConstants::trackStateMask, + result.lastTrackIndex); + const std::size_t currentTrackIndex = trackStateProxy.index(); + + // Set the trackStateProxy components with the state from the + // ongoing propagation + { + trackStateProxy.setReferenceSurface(surface->getSharedPtr()); + // Bind the transported state to the current surface + auto res = stepper.boundState(state.stepping, *surface, false, + freeToBoundCorrection); + if (!res.ok()) { + result.result = res.error(); + return; } + const auto& [boundParams, jacobian, pathLength] = *res; - // Get and set the type flags - auto typeFlags = trackStateProxy.typeFlags(); - typeFlags.set(TrackStateFlag::ParameterFlag); - if (surfaceHasMaterial) { - ACTS_DEBUG(" It is material."); - typeFlags.set(TrackStateFlag::MaterialFlag); - } + // Fill the track state + trackStateProxy.predicted() = boundParams.parameters(); + trackStateProxy.predictedCovariance() = state.stepping.cov; - // Set hole only, if we are on a sensitive surface - if (surfaceIsSensitive && precedingMeasurementExists) { - ACTS_DEBUG(" It is a hole."); - typeFlags.set(TrackStateFlag::HoleFlag); - // Count the missed surface - result.missedActiveSurfaces.push_back(surface); - } + trackStateProxy.jacobian() = jacobian; + trackStateProxy.pathLength() = pathLength; + } - result.lastTrackIndex = currentTrackIndex; + // Get and set the type flags + auto typeFlags = trackStateProxy.typeFlags(); + typeFlags.set(TrackStateFlag::ParameterFlag); + if (surfaceHasMaterial) { + ACTS_DEBUG(" It is material."); + typeFlags.set(TrackStateFlag::MaterialFlag); + } - ++result.processedStates; - } else { - ACTS_DEBUG(" The surface contains no measurement/material/hole."); + // Set hole only, if we are on a sensitive surface + if (surfaceIsSensitive && precedingMeasurementExists) { + ACTS_DEBUG(" It is a hole."); + typeFlags.set(TrackStateFlag::HoleFlag); + // Count the missed surface + result.missedActiveSurfaces.push_back(surface); } + + result.lastTrackIndex = currentTrackIndex; + + ++result.processedStates; + + return; } + + ACTS_DEBUG(" The surface contains no measurement/material/hole."); + return; } }; diff --git a/Core/include/Acts/Utilities/BinningData.hpp b/Core/include/Acts/Utilities/BinningData.hpp index 9b674109853..9ad848ba152 100644 --- a/Core/include/Acts/Utilities/BinningData.hpp +++ b/Core/include/Acts/Utilities/BinningData.hpp @@ -531,7 +531,7 @@ class BinningData { /// String screen output method /// @param indent the current indentation /// @return a string containing the screen information - std::string toString(const std::string& indent) const { + std::string toString(const std::string& indent = "") const { std::stringstream sl; sl << indent << "BinngingData object:" << '\n'; sl << indent << " - type : " << static_cast(type) diff --git a/Core/src/Geometry/CMakeLists.txt b/Core/src/Geometry/CMakeLists.txt index 00d9618f062..fa5d70c1c02 100644 --- a/Core/src/Geometry/CMakeLists.txt +++ b/Core/src/Geometry/CMakeLists.txt @@ -35,6 +35,7 @@ target_sources( Volume.cpp VolumeBounds.cpp CylinderVolumeStack.cpp + Portal.cpp GridPortalLink.cpp GridPortalLinkMerging.cpp TrivialPortalLink.cpp diff --git a/Core/src/Geometry/Portal.cpp b/Core/src/Geometry/Portal.cpp new file mode 100644 index 00000000000..4b1e65356af --- /dev/null +++ b/Core/src/Geometry/Portal.cpp @@ -0,0 +1,332 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Geometry/Portal.hpp" + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/PortalLinkBase.hpp" +#include "Acts/Geometry/TrivialPortalLink.hpp" +#include "Acts/Surfaces/RegularSurface.hpp" +#include "Acts/Utilities/BinningType.hpp" + +#include +#include +#include + +namespace Acts { + +const char* PortalMergingException::what() const noexcept { + return "Failure to merge portals"; +} + +const char* PortalFusingException::what() const noexcept { + return "Failure to fuse portals"; +} + +Portal::Portal(Direction direction, std::unique_ptr link) { + if (link == nullptr) { + throw std::invalid_argument("Link must not be null"); + } + + m_surface = link->surfacePtr(); + + if (direction == Direction::AlongNormal) { + m_alongNormal = std::move(link); + } else { + m_oppositeNormal = std::move(link); + } +} + +Portal::Portal(Direction direction, std::shared_ptr surface, + TrackingVolume& volume) + : Portal(direction, + std::make_unique(std::move(surface), volume)) {} + +Portal::Portal(const GeometryContext& gctx, + std::unique_ptr alongNormal, + std::unique_ptr oppositeNormal) { + if (alongNormal == nullptr && oppositeNormal == nullptr) { + throw std::invalid_argument("At least one link must be provided"); + } + + if (alongNormal != nullptr) { + setLink(gctx, Direction::AlongNormal, std::move(alongNormal)); + } + if (oppositeNormal != nullptr) { + setLink(gctx, Direction::OppositeNormal, std::move(oppositeNormal)); + } +} + +Portal::Portal(const GeometryContext& gctx, Arguments&& args) { + if (!args.alongNormal.surface && !args.oppositeNormal.surface) { + throw std::invalid_argument("At least one link must be provided"); + } + + if (args.alongNormal.surface) { + setLink(gctx, Direction::AlongNormal, + std::make_unique( + std::move(args.alongNormal.surface), *args.alongNormal.volume)); + } + if (args.oppositeNormal.surface) { + setLink(gctx, Direction::OppositeNormal, + std::make_unique( + std::move(args.oppositeNormal.surface), + *args.oppositeNormal.volume)); + } +} + +void Portal::setLink(const GeometryContext& gctx, Direction direction, + std::unique_ptr link) { + if (link == nullptr) { + throw std::invalid_argument("Link must not be null"); + } + + auto& target = + direction == Direction::AlongNormal ? m_alongNormal : m_oppositeNormal; + const auto& other = + direction == Direction::AlongNormal ? m_oppositeNormal : m_alongNormal; + + // check if surfaces are identical + if (m_surface != nullptr && + !isSameSurface(gctx, link->surface(), *m_surface)) { + throw PortalFusingException(); + } + + // check if they both have material but are not the same surface + if (m_surface != nullptr && (m_surface.get() != &link->surface()) && + link->surface().surfaceMaterial() != nullptr && + m_surface->surfaceMaterial() != nullptr) { + throw PortalFusingException(); + } + + target = std::move(link); + + if (other == nullptr) { + // We don't have an existing surface, take the one we just got + m_surface = target->surfacePtr(); + return; + } + + if (target->surface().surfaceMaterial() != nullptr) { + // new link has material: assign that to existing link + m_surface = target->surfacePtr(); + other->setSurface(m_surface); + } else { + // none have material, or the existing surface had material: assign the + // existing surface by convention + target->setSurface(m_surface); + } +} + +void Portal::setLink(const GeometryContext& gctx, Direction direction, + std::shared_ptr surface, + TrackingVolume& volume) { + setLink(gctx, direction, + std::make_unique(std::move(surface), volume)); +} + +const PortalLinkBase* Portal::getLink(Direction direction) const { + if (direction == Direction::AlongNormal) { + return m_alongNormal.get(); + } else { + return m_oppositeNormal.get(); + } +} + +Result Portal::resolveVolume( + const GeometryContext& gctx, const Vector3& position, + const Vector3& direction) const { + assert(m_surface != nullptr); + const Vector3 normal = m_surface->normal(gctx, position); + Direction side = Direction::fromScalarZeroAsPositive(normal.dot(direction)); + + const PortalLinkBase* link = side == Direction::AlongNormal + ? m_alongNormal.get() + : m_oppositeNormal.get(); + + if (link == nullptr) { + // no link is attached in this direction => this is the end of the world as + // we know it. (i feel fine) + return nullptr; + } else { + auto res = link->resolveVolume(gctx, position); + if (!res.ok()) { + return res.error(); + } + return *res; + } +} + +bool Portal::isValid() const { + return m_alongNormal != nullptr || m_oppositeNormal != nullptr; +} + +const RegularSurface& Portal::surface() const { + assert(m_surface != nullptr); + return *m_surface; +} + +Portal Portal::merge(const GeometryContext& gctx, Portal& aPortal, + Portal& bPortal, BinningValue direction, + const Logger& logger) { + ACTS_DEBUG("Merging to portals along " << direction); + + if (&aPortal == &bPortal) { + ACTS_ERROR("Cannot merge a portal with itself"); + throw PortalMergingException{}; + } + + if (aPortal.m_surface->surfaceMaterial() != nullptr || + bPortal.m_surface->surfaceMaterial() != nullptr) { + ACTS_ERROR("Cannot merge portals with material"); + throw PortalMergingException{}; + } + + std::unique_ptr mergedAlongNormal = nullptr; + std::unique_ptr mergedOppositeNormal = nullptr; + + bool aHasAlongNormal = aPortal.m_alongNormal != nullptr; + bool aHasOppositeNormal = aPortal.m_oppositeNormal != nullptr; + bool bHasAlongNormal = bPortal.m_alongNormal != nullptr; + bool bHasOppositeNormal = bPortal.m_oppositeNormal != nullptr; + + if (aHasAlongNormal != bHasAlongNormal || + aHasOppositeNormal != bHasOppositeNormal) { + ACTS_ERROR("Portals do not have the same links attached"); + throw PortalMergingException(); + } + + if (aPortal.m_alongNormal != nullptr) { + if (bPortal.m_alongNormal == nullptr) { + ACTS_ERROR( + "Portal A has link along normal, while b does not. This is not " + "supported"); + throw PortalMergingException(); + } + + ACTS_VERBOSE("Portals have links along normal, merging"); + mergedAlongNormal = PortalLinkBase::merge(std::move(aPortal.m_alongNormal), + std::move(bPortal.m_alongNormal), + direction, logger); + } + + if (aPortal.m_oppositeNormal != nullptr) { + if (bPortal.m_oppositeNormal == nullptr) { + ACTS_ERROR( + "Portal A has link opposite normal, while b does not. This is not " + "supported"); + throw PortalMergingException(); + } + + ACTS_VERBOSE("Portals have links opposite normal, merging"); + mergedOppositeNormal = PortalLinkBase::merge( + std::move(aPortal.m_oppositeNormal), + std::move(bPortal.m_oppositeNormal), direction, logger); + } + + aPortal.m_surface.reset(); + bPortal.m_surface.reset(); + return Portal{gctx, std::move(mergedAlongNormal), + std::move(mergedOppositeNormal)}; +} + +Portal Portal::fuse(const GeometryContext& gctx, Portal& aPortal, + Portal& bPortal, const Logger& logger) { + ACTS_DEBUG("Fusing two portals"); + if (&aPortal == &bPortal) { + ACTS_ERROR("Cannot merge a portal with itself"); + throw PortalMergingException{}; + } + + bool aHasAlongNormal = aPortal.m_alongNormal != nullptr; + bool aHasOppositeNormal = aPortal.m_oppositeNormal != nullptr; + bool bHasAlongNormal = bPortal.m_alongNormal != nullptr; + bool bHasOppositeNormal = bPortal.m_oppositeNormal != nullptr; + + if (aPortal.m_surface == nullptr || bPortal.m_surface == nullptr) { + ACTS_ERROR("Portals have no surface"); + throw PortalFusingException(); + } + + if (aPortal.m_surface->associatedDetectorElement() != nullptr || + bPortal.m_surface->associatedDetectorElement() != nullptr) { + ACTS_ERROR("Cannot fuse portals with detector elements"); + throw PortalFusingException(); + } + + if (!isSameSurface(gctx, *aPortal.m_surface, *bPortal.m_surface)) { + ACTS_ERROR("Portals have different surfaces"); + throw PortalFusingException(); + } + + if (aPortal.m_surface->surfaceMaterial() != nullptr && + bPortal.m_surface->surfaceMaterial() != nullptr) { + ACTS_ERROR("Cannot fuse portals if both have material"); + throw PortalFusingException(); + } + + if (aHasAlongNormal == bHasAlongNormal || + aHasOppositeNormal == bHasOppositeNormal) { + ACTS_ERROR("Portals have the same links attached"); + throw PortalFusingException(); + } + + 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)}; + } 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)}; + } +} + +bool Portal::isSameSurface(const GeometryContext& gctx, const Surface& a, + const Surface& b) { + if (&a == &b) { + return true; + } + + if (a.type() != b.type()) { + return false; + } + + if (a.bounds() != b.bounds()) { + return false; + } + + if (!a.transform(gctx).isApprox(b.transform(gctx), + s_transformEquivalentTolerance)) { + return false; + } + + return true; +}; + +void Portal::fill(TrackingVolume& volume) { + if (m_alongNormal != nullptr && m_oppositeNormal != nullptr) { + throw std::logic_error{"Portal is already filled"}; + } + + if (m_surface == nullptr) { + throw std::logic_error{"Portal has no existing link set, can't fill"}; + } + + if (m_alongNormal == nullptr) { + m_alongNormal = std::make_unique(m_surface, volume); + } else { + assert(m_oppositeNormal == nullptr); + m_oppositeNormal = std::make_unique(m_surface, volume); + } +} + +} // namespace Acts diff --git a/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/Geant4Manager.hpp b/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/Geant4Manager.hpp index b9d6f427547..726b6378ce1 100644 --- a/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/Geant4Manager.hpp +++ b/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/Geant4Manager.hpp @@ -37,12 +37,11 @@ class Geant4Manager; /// and loading the Geant4 library which should reset it to its original state. struct Geant4Handle { std::mutex mutex; - int logLevel{}; std::unique_ptr runManager; G4VUserPhysicsList *physicsList; std::string physicsListName; - Geant4Handle(int logLevel, std::unique_ptr runManager, + Geant4Handle(std::unique_ptr runManager, std::unique_ptr physicsList, std::string physicsListName); Geant4Handle(const Geant4Handle &) = delete; @@ -66,12 +65,11 @@ class Geant4Manager { std::shared_ptr currentHandle() const; /// This can only be called once due to Geant4 limitations - std::shared_ptr createHandle(int logLevel, - const std::string &physicsList); + std::shared_ptr createHandle(const std::string &physicsList); /// This can only be called once due to Geant4 limitations std::shared_ptr createHandle( - int logLevel, std::unique_ptr physicsList, + std::unique_ptr physicsList, std::string physicsListName); /// Registers a named physics list factory to the manager for easy diff --git a/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/SensitiveSurfaceMapper.hpp b/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/SensitiveSurfaceMapper.hpp index bdb1f3753a0..3a4bc3f2f3f 100644 --- a/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/SensitiveSurfaceMapper.hpp +++ b/Examples/Algorithms/Geant4/include/ActsExamples/Geant4/SensitiveSurfaceMapper.hpp @@ -28,34 +28,37 @@ class Surface; namespace ActsExamples { -// Helper struct to find the sensitive surface candidates -struct SensitiveCandidates { - std::shared_ptr trackingGeometry = nullptr; - /// Find the sensitive surfaces for a given position +struct SensitiveCandidatesBase { + /// Get the sensitive surfaces for a given position /// - /// This fulfills the concept of a SensitiveCandidates + /// @param gctx the geometry context + /// @param position the position to look for sensitive surfaces + /// + /// @return a vector of sensitive surfaces + virtual std::vector queryPosition( + const Acts::GeometryContext& gctx, + const Acts::Vector3& position) const = 0; + + /// Get all sensitive surfaces /// /// @param gctx the geometry context /// @param position the position to look for sensitive surfaces /// /// @return a vector of sensitive surfaces - std::vector operator()( - const Acts::GeometryContext& gctx, const Acts::Vector3& position) const { - std::vector surfaces; - - if (trackingGeometry != nullptr) { - auto layer = trackingGeometry->associatedLayer(gctx, position); - - if (layer->surfaceArray() != nullptr) { - for (const auto& surface : layer->surfaceArray()->surfaces()) { - if (surface->associatedDetectorElement() != nullptr) { - surfaces.push_back(surface); - } - } - } - } - return surfaces; - } + virtual std::vector queryAll() const = 0; + + virtual ~SensitiveCandidatesBase() = default; +}; + +/// Implementation of the SensitiveCandidates for Gen1 geometry +struct SensitiveCandidates : public SensitiveCandidatesBase { + std::shared_ptr trackingGeometry = nullptr; + + std::vector queryPosition( + const Acts::GeometryContext& gctx, + const Acts::Vector3& position) const override; + + std::vector queryAll() const override; }; /// This Mapper takes a (non-const) Geant4 geometry and maps @@ -72,9 +75,6 @@ class SensitiveSurfaceMapper { /// This prefix is used to indicate a sensitive volume that is matched constexpr static std::string_view mappingPrefix = "ActsSensitive#"; - using SensitiveCandidates = std::function( - const Acts::GeometryContext&, const Acts::Vector3&)>; - /// Configuration struct for the surface mapper struct Config { /// For which G4 material names we try to find a mapping @@ -84,7 +84,7 @@ class SensitiveSurfaceMapper { std::vector volumeMappings; /// The sensitive surfaces that are being mapped to - SensitiveCandidates candidateSurfaces; + std::shared_ptr candidateSurfaces; }; /// State object that coutns the assignments and makes @@ -94,6 +94,10 @@ class SensitiveSurfaceMapper { /// there can be replicas) std::multimap g4VolumeToSurfaces; + + /// Record of the missing volumes + std::vector> + missingVolumes; }; /// Constructor with: @@ -118,6 +122,20 @@ class SensitiveSurfaceMapper { G4VPhysicalVolume* g4PhysicalVolume, const Acts::Transform3& motherTransform) const; + /// Function that checks the success of the mapping, and exposes + /// some additional information for debugging + /// + /// @param state state object after a call to remapSensitiveNames + /// @param gctx the geometry context + /// @param writeMissingG4VolsAsObj write the Geant4 volumes that are + /// not mapped to 'missing_g4_volumes.obj' in the working directory + /// @param writeMissingSurfacesAsObj write the sensitive surfaces that + /// where not mapped to 'missing_acts_surfaces.obj' in the working directory + /// @return Returns true only if all sensitive surfaces where mapped + bool checkMapping(const State& state, const Acts::GeometryContext& gctx, + bool writeMissingG4VolsAsObj = false, + bool writeMissingSurfacesAsObj = false) const; + protected: /// Configuration object Config m_cfg; diff --git a/Examples/Algorithms/Geant4/src/Geant4Manager.cpp b/Examples/Algorithms/Geant4/src/Geant4Manager.cpp index ffadf6e2b18..46e2dbe867c 100644 --- a/Examples/Algorithms/Geant4/src/Geant4Manager.cpp +++ b/Examples/Algorithms/Geant4/src/Geant4Manager.cpp @@ -31,12 +31,10 @@ namespace ActsExamples { -Geant4Handle::Geant4Handle(int _logLevel, - std::unique_ptr _runManager, +Geant4Handle::Geant4Handle(std::unique_ptr _runManager, std::unique_ptr _physicsList, std::string _physicsListName) - : logLevel(_logLevel), - runManager(std::move(_runManager)), + : runManager(std::move(_runManager)), physicsList(_physicsList.release()), physicsListName(std::move(_physicsListName)) { if (runManager == nullptr) { @@ -81,12 +79,12 @@ std::shared_ptr Geant4Manager::currentHandle() const { } std::shared_ptr Geant4Manager::createHandle( - int logLevel, const std::string& physicsList) { - return createHandle(logLevel, createPhysicsList(physicsList), physicsList); + const std::string& physicsList) { + return createHandle(createPhysicsList(physicsList), physicsList); } std::shared_ptr Geant4Manager::createHandle( - int logLevel, std::unique_ptr physicsList, + std::unique_ptr physicsList, std::string physicsListName) { if (!m_handle.expired()) { throw std::runtime_error("creating a second handle is prohibited"); @@ -100,7 +98,7 @@ std::shared_ptr Geant4Manager::createHandle( auto runManager = std::unique_ptr( G4RunManagerFactory::CreateRunManager(G4RunManagerType::SerialOnly)); - auto handle = std::make_shared(logLevel, std::move(runManager), + auto handle = std::make_shared(std::move(runManager), std::move(physicsList), std::move(physicsListName)); diff --git a/Examples/Algorithms/Geant4/src/Geant4Simulation.cpp b/Examples/Algorithms/Geant4/src/Geant4Simulation.cpp index fce1c03db1e..dbfc14242ae 100644 --- a/Examples/Algorithms/Geant4/src/Geant4Simulation.cpp +++ b/Examples/Algorithms/Geant4/src/Geant4Simulation.cpp @@ -90,6 +90,8 @@ void ActsExamples::Geant4SimulationBase::commonInitialization() { runManager().SetUserInitialization(m_detectorConstruction); runManager().InitializeGeometry(); } + + m_geant4Instance->tweakLogging(m_geant4Level); } G4RunManager& ActsExamples::Geant4SimulationBase::runManager() const { @@ -175,10 +177,10 @@ ActsExamples::Geant4SimulationBase::geant4Handle() const { ActsExamples::Geant4Simulation::Geant4Simulation(const Config& cfg, Acts::Logging::Level level) : Geant4SimulationBase(cfg, "Geant4Simulation", level), m_cfg(cfg) { - m_geant4Instance = m_cfg.geant4Handle - ? m_cfg.geant4Handle - : Geant4Manager::instance().createHandle( - m_geant4Level, m_cfg.physicsList); + m_geant4Instance = + m_cfg.geant4Handle + ? m_cfg.geant4Handle + : Geant4Manager::instance().createHandle(m_cfg.physicsList); if (m_geant4Instance->physicsListName != m_cfg.physicsList) { throw std::runtime_error("inconsistent physics list"); } @@ -282,8 +284,13 @@ ActsExamples::Geant4Simulation::Geant4Simulation(const Config& cfg, "Remapping selected volumes from Geant4 to Acts::Surface::GeometryID"); cfg.sensitiveSurfaceMapper->remapSensitiveNames( sState, Acts::GeometryContext{}, g4World, Acts::Transform3::Identity()); - ACTS_INFO("Remapping successful for " << sState.g4VolumeToSurfaces.size() - << " selected volumes."); + + auto allSurfacesMapped = cfg.sensitiveSurfaceMapper->checkMapping( + sState, Acts::GeometryContext{}, false, false); + if (!allSurfacesMapped) { + ACTS_WARNING( + "Not all sensitive surfaces have been mapped to Geant4 volumes!"); + } sensitiveSteppingActionAccess->assignSurfaceMapping( sState.g4VolumeToSurfaces); @@ -334,7 +341,6 @@ ActsExamples::Geant4MaterialRecording::Geant4MaterialRecording( m_cfg.geant4Handle ? m_cfg.geant4Handle : Geant4Manager::instance().createHandle( - m_geant4Level, std::make_unique( m_logger->cloneWithSuffix("MaterialPhysicsList")), physicsListName); diff --git a/Examples/Algorithms/Geant4/src/SensitiveSurfaceMapper.cpp b/Examples/Algorithms/Geant4/src/SensitiveSurfaceMapper.cpp index 731aa79ac1e..341e0531b07 100644 --- a/Examples/Algorithms/Geant4/src/SensitiveSurfaceMapper.cpp +++ b/Examples/Algorithms/Geant4/src/SensitiveSurfaceMapper.cpp @@ -9,15 +9,125 @@ #include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp" #include "Acts/Definitions/Units.hpp" +#include "Acts/Surfaces/AnnulusBounds.hpp" +#include "Acts/Surfaces/ConvexPolygonBounds.hpp" +#include "Acts/Visualization/GeometryView3D.hpp" +#include "Acts/Visualization/ObjVisualization3D.hpp" #include #include +#include #include +#include #include #include #include +#include #include +#include +#include +#include + +// Add some type traits for boost::geometry, so we can use the machinery +// directly with Acts::Vector2 / Eigen::Matrix +namespace boost::geometry::traits { + +template +struct tag> { + using type = point_tag; +}; +template +struct dimension> : std::integral_constant {}; +template +struct coordinate_type> { + using type = T; +}; +template +struct coordinate_system> { + using type = boost::geometry::cs::cartesian; +}; + +template +struct access, Index> { + static_assert(Index < D, "Out of range"); + using Point = Eigen::Matrix; + using CoordinateType = typename coordinate_type::type; + static inline CoordinateType get(Point const& p) { return p[Index]; } + static inline void set(Point& p, CoordinateType const& value) { + p[Index] = value; + } +}; + +} // namespace boost::geometry::traits + +namespace { +void writeG4Polyhedron( + Acts::IVisualization3D& visualizer, const G4Polyhedron& polyhedron, + const Acts::Transform3& trafo = Acts::Transform3::Identity(), + Acts::ColorRGB color = {0, 0, 0}) { + constexpr double convertLength = CLHEP::mm / Acts::UnitConstants::mm; + + for (int i = 1; i <= polyhedron.GetNoFacets(); ++i) { + // This is a bit ugly but I didn't find an easy way to compute this + // beforehand. + constexpr std::size_t maxPoints = 1000; + G4Point3D points[maxPoints]; + int nPoints = 0; + polyhedron.GetFacet(i, nPoints, points); + assert(static_cast(nPoints) < maxPoints); + + std::vector faces; + for (int j = 0; j < nPoints; ++j) { + faces.emplace_back(points[j][0] * convertLength, + points[j][1] * convertLength, + points[j][2] * convertLength); + faces.back() = trafo * faces.back(); + } + + visualizer.face(faces, color); + } +} +} // namespace + +std::vector +ActsExamples::SensitiveCandidates::queryPosition( + const Acts::GeometryContext& gctx, const Acts::Vector3& position) const { + std::vector surfaces; + + if (trackingGeometry == nullptr) { + return surfaces; + } + + // In case we do not find a layer at this position for whatever reason + const auto layer = trackingGeometry->associatedLayer(gctx, position); + if (layer == nullptr) { + return surfaces; + } + + const auto surfaceArray = layer->surfaceArray(); + if (surfaceArray == nullptr) { + return surfaces; + } + + for (const auto& surface : surfaceArray->surfaces()) { + if (surface->associatedDetectorElement() != nullptr) { + surfaces.push_back(surface); + } + } + return surfaces; +} + +std::vector ActsExamples::SensitiveCandidates::queryAll() + const { + std::vector surfaces; + + const bool restrictToSensitives = true; + trackingGeometry->visitSurfaces( + [&](auto surface) { surfaces.push_back(surface); }, restrictToSensitives); + + return surfaces; +} ActsExamples::SensitiveSurfaceMapper::SensitiveSurfaceMapper( const Config& cfg, std::unique_ptr logger) @@ -30,33 +140,44 @@ void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames( // Make sure the unit conversion is correct constexpr double convertLength = CLHEP::mm / Acts::UnitConstants::mm; + auto g4ToActsVector = [](const G4ThreeVector& g4vec) { + return Acts::Vector3(g4vec[0] * convertLength, g4vec[1] * convertLength, + g4vec[2] * convertLength); + }; + + auto actsToG4Vec = [](const Acts::Vector3& actsVec) { + return G4ThreeVector(actsVec[0] / convertLength, actsVec[1] / convertLength, + actsVec[2] / convertLength); + }; + auto g4LogicalVolume = g4PhysicalVolume->GetLogicalVolume(); auto g4SensitiveDetector = g4LogicalVolume->GetSensitiveDetector(); // Get the transform of the G4 object - auto g4Translation = g4PhysicalVolume->GetTranslation(); - auto g4Rotation = g4PhysicalVolume->GetRotation(); - Acts::Vector3 g4RelPosition(g4Translation[0] * convertLength, - g4Translation[1] * convertLength, - g4Translation[2] * convertLength); - Acts::Translation3 translation(g4RelPosition); - Acts::Transform3 transform; - if (g4Rotation == nullptr) { - transform = motherTransform * translation; - } else { - Acts::RotationMatrix3 rotation; - rotation << g4Rotation->xx(), g4Rotation->yx(), g4Rotation->zx(), - g4Rotation->xy(), g4Rotation->yy(), g4Rotation->zy(), g4Rotation->xz(), - g4Rotation->yz(), g4Rotation->zz(); - transform = motherTransform * (translation * rotation); - } - Acts::Vector3 g4AbsPosition = transform * Acts::Vector3::Zero(); + Acts::Transform3 localG4ToGlobal; + { + auto g4Translation = g4PhysicalVolume->GetTranslation(); + auto g4Rotation = g4PhysicalVolume->GetRotation(); + Acts::Vector3 g4RelPosition = g4ToActsVector(g4Translation); + Acts::Translation3 translation(g4RelPosition); + if (g4Rotation == nullptr) { + localG4ToGlobal = motherTransform * translation; + } else { + Acts::RotationMatrix3 rotation; + rotation << g4Rotation->xx(), g4Rotation->yx(), g4Rotation->zx(), + g4Rotation->xy(), g4Rotation->yy(), g4Rotation->zy(), + g4Rotation->xz(), g4Rotation->yz(), g4Rotation->zz(); + localG4ToGlobal = motherTransform * (translation * rotation); + } + } + + Acts::Vector3 g4AbsPosition = localG4ToGlobal * Acts::Vector3::Zero(); if (G4int nDaughters = g4LogicalVolume->GetNoDaughters(); nDaughters > 0) { // Step down to all daughters for (G4int id = 0; id < nDaughters; ++id) { remapSensitiveNames(state, gctx, g4LogicalVolume->GetDaughter(id), - transform); + localG4ToGlobal); } return; } @@ -72,45 +193,154 @@ void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames( std::find(m_cfg.volumeMappings.begin(), m_cfg.volumeMappings.end(), volumeName) != m_cfg.volumeMappings.end(); - if (isSensitive || isMappedMaterial || isMappedVolume) { - // Prepare the mapped surface - const Acts::Surface* mappedSurface = nullptr; + if (!(isSensitive || isMappedMaterial || isMappedVolume)) { + ACTS_VERBOSE("Did not try mapping '" + << g4PhysicalVolume->GetName() << "' at " + << g4AbsPosition.transpose() + << " because g4SensitiveDetector (=" << g4SensitiveDetector + << ") is null and volume name (=" << volumeName + << ") and material name (=" << volumeMaterialName + << ") were not found"); + return; + } + + // Prepare the mapped surface + const Acts::Surface* mappedSurface = nullptr; + + std::vector candidateSurfaces; + const auto g4Polyhedron = g4LogicalVolume->GetSolid()->GetPolyhedron(); + for (int i = 1; i < g4Polyhedron->GetNoVertices(); ++i) { + auto vtx = g4ToActsVector(g4Polyhedron->GetVertex(i)); + auto vtxGlobal = localG4ToGlobal * vtx; + + candidateSurfaces = m_cfg.candidateSurfaces->queryPosition(gctx, vtxGlobal); + + if (!candidateSurfaces.empty()) { + break; + } + } + + // Fall back to query all surfaces + if (candidateSurfaces.empty()) { + ACTS_DEBUG("No candidate surfaces for volume '" << volumeName << "' at " + << g4AbsPosition.transpose() + << ", query all surfaces"); + candidateSurfaces = m_cfg.candidateSurfaces->queryAll(); + } + + ACTS_VERBOSE("Found " << candidateSurfaces.size() + << " candidate surfaces for " << volumeName); + + for (const auto& candidateSurface : candidateSurfaces) { + if (candidateSurface->center(gctx).isApprox(g4AbsPosition)) { + ACTS_VERBOSE("Successful match with center matching"); + mappedSurface = candidateSurface; + break; + } else if (candidateSurface->bounds().type() == + Acts::SurfaceBounds::eAnnulus) { + const auto& bounds = + *static_cast(&candidateSurface->bounds()); + + const auto vertices = bounds.vertices(0); - // Get the candidate surfaces - auto candidateSurfaces = m_cfg.candidateSurfaces(gctx, g4AbsPosition); - for (const auto& candidateSurface : candidateSurfaces) { - // Check for center matching at the moment (needs to be extended) - if (candidateSurface->center(gctx).isApprox(g4AbsPosition)) { + constexpr bool clockwise = false; + constexpr bool closed = false; + using Polygon = + boost::geometry::model::polygon; + + Polygon poly; + boost::geometry::assign_points(poly, vertices); + + Acts::Vector2 boundsCentroidSurfaceFrame = Acts::Vector2::Zero(); + boost::geometry::centroid(poly, boundsCentroidSurfaceFrame); + + Acts::Vector3 boundsCentroidGlobal{boundsCentroidSurfaceFrame[0], + boundsCentroidSurfaceFrame[1], 0.0}; + boundsCentroidGlobal = + candidateSurface->transform(gctx) * boundsCentroidGlobal; + + const auto boundsCentroidG4Frame = + localG4ToGlobal.inverse() * boundsCentroidGlobal; + + if (g4LogicalVolume->GetSolid()->Inside( + actsToG4Vec(boundsCentroidG4Frame)) != EInside::kOutside) { + ACTS_VERBOSE("Successful match with centroid matching"); mappedSurface = candidateSurface; break; } } + } - // A mapped surface was found, a new name will be set that G4PhysVolume - if (mappedSurface != nullptr) { - ACTS_VERBOSE("Found matching surface " << mappedSurface->geometryId() - << " at position " - << g4RelPosition.transpose()); - // Check if the prefix is not yet assigned - if (volumeName.find(mappingPrefix) == std::string::npos) { - // Set the new name - std::string mappedName = std::string(mappingPrefix) + volumeName; - g4PhysicalVolume->SetName(mappedName); - } - // Insert into the multi-map - state.g4VolumeToSurfaces.insert({g4PhysicalVolume, mappedSurface}); - } else { - ACTS_VERBOSE("No mapping found for '" - << volumeName << "' with material '" << volumeMaterialName - << "' at position " << g4RelPosition.transpose()); + if (mappedSurface == nullptr) { + ACTS_DEBUG("No mapping found for '" + << volumeName << "' with material '" << volumeMaterialName + << "' at position " << g4AbsPosition.transpose()); + state.missingVolumes.emplace_back(g4PhysicalVolume, localG4ToGlobal); + return; + } + + // A mapped surface was found, a new name will be set that G4PhysVolume + ACTS_DEBUG("Matched " << volumeName << " to " << mappedSurface->geometryId() + << " at position " << g4AbsPosition.transpose()); + // Check if the prefix is not yet assigned + if (volumeName.find(mappingPrefix) == std::string::npos) { + // Set the new name + std::string mappedName = std::string(mappingPrefix) + volumeName; + g4PhysicalVolume->SetName(mappedName); + } + // Insert into the multi-map + state.g4VolumeToSurfaces.insert({g4PhysicalVolume, mappedSurface}); +} + +bool ActsExamples::SensitiveSurfaceMapper::checkMapping( + const State& state, const Acts::GeometryContext& gctx, + bool writeMissingG4VolsAsObj, bool writeMissingSurfacesAsObj) const { + auto allSurfaces = m_cfg.candidateSurfaces->queryAll(); + std::sort(allSurfaces.begin(), allSurfaces.end()); + + std::vector found; + for (const auto [_, surfacePtr] : state.g4VolumeToSurfaces) { + found.push_back(surfacePtr); + } + std::sort(found.begin(), found.end()); + auto newEnd = std::unique(found.begin(), found.end()); + found.erase(newEnd, found.end()); + + std::vector missing; + std::set_difference(allSurfaces.begin(), allSurfaces.end(), found.begin(), + found.end(), std::back_inserter(missing)); + + ACTS_INFO("Number of overall sensitive surfaces: " << allSurfaces.size()); + ACTS_INFO("Number of mapped volume->surface mappings: " + << state.g4VolumeToSurfaces.size()); + ACTS_INFO( + "Number of sensitive surfaces that are not mapped: " << missing.size()); + ACTS_INFO("Number of G4 volumes without a matching Surface: " + << state.missingVolumes.size()); + + if (writeMissingG4VolsAsObj) { + Acts::ObjVisualization3D visualizer; + for (const auto& [g4vol, trafo] : state.missingVolumes) { + auto polyhedron = g4vol->GetLogicalVolume()->GetSolid()->GetPolyhedron(); + writeG4Polyhedron(visualizer, *polyhedron, trafo); } - } else { - ACTS_VERBOSE("Did not try mapping '" - << g4PhysicalVolume->GetName() << "' at " - << g4RelPosition.transpose() - << " because g4SensitiveDetector (=" << g4SensitiveDetector - << ") is null and volume name (=" << volumeName - << ") and material name (=" << volumeMaterialName - << ") were not found"); + + std::ofstream os("missing_g4_volumes.obj"); + visualizer.write(os); } + + if (writeMissingSurfacesAsObj) { + Acts::ObjVisualization3D visualizer; + Acts::ViewConfig vcfg; + vcfg.nSegments = 720; + for (auto srf : missing) { + Acts::GeometryView3D::drawSurface(visualizer, *srf, gctx, + Acts::Transform3::Identity(), vcfg); + } + + std::ofstream os("missing_acts_surfaces.obj"); + visualizer.write(os); + } + + return missing.empty(); } diff --git a/Examples/Python/CMakeLists.txt b/Examples/Python/CMakeLists.txt index 1ab76a804a2..5ed407bd889 100644 --- a/Examples/Python/CMakeLists.txt +++ b/Examples/Python/CMakeLists.txt @@ -97,8 +97,11 @@ endif() if(ACTS_BUILD_PLUGIN_TRACCC) target_link_libraries(ActsPythonBindings PUBLIC ActsPluginDetray) target_sources(ActsPythonBindings PRIVATE src/Detray.cpp) + target_link_libraries(ActsPythonBindings PUBLIC ActsPluginCovfie) + target_sources(ActsPythonBindings PRIVATE src/Covfie.cpp) else() target_sources(ActsPythonBindings PRIVATE src/DetrayStub.cpp) + target_sources(ActsPythonBindings PRIVATE src/CovfieStub.cpp) endif() if(ACTS_BUILD_PLUGIN_ACTSVG) diff --git a/Examples/Python/python/acts/examples/simulation.py b/Examples/Python/python/acts/examples/simulation.py index 82d7144006f..125ce422c16 100644 --- a/Examples/Python/python/acts/examples/simulation.py +++ b/Examples/Python/python/acts/examples/simulation.py @@ -708,7 +708,7 @@ def addGeant4( smmConfig.volumeMappings = volumeMappings smmConfig.materialMappings = materialMappings sensitiveMapper = SensitiveSurfaceMapper.create( - smmConfig, acts.logging.INFO, trackingGeometry + smmConfig, customLogLevel(), trackingGeometry ) # Simulation diff --git a/Examples/Python/src/Covfie.cpp b/Examples/Python/src/Covfie.cpp new file mode 100644 index 00000000000..d4ad5d30b00 --- /dev/null +++ b/Examples/Python/src/Covfie.cpp @@ -0,0 +1,55 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Plugins/Covfie/FieldConversion.hpp" +#include "Acts/Plugins/Python/Utilities.hpp" + +#include +#include + +namespace py = pybind11; +using namespace pybind11::literals; + +namespace Acts::Python { + +namespace { +template +void declareCovfieField(py::module& m, const std::string& fieldName) { + using view_t = typename field_t::view_t; + m.def("toView", + [](const field_t& field) { return typename field_t::view_t(field); }); + py::class_>(m, fieldName.c_str()); + py::class_>( + m, (fieldName + std::string("View")).c_str()) + .def("at", &view_t::template at); +} +} // namespace + +void addCovfie(Context& ctx) { + auto main = ctx.get("main"); + auto m = main.def_submodule("covfie", "Submodule for covfie conversion"); + + declareCovfieField(m, + "CovfieConstantField"); + declareCovfieField( + m, "CovfieAffineLinearStridedField"); + + m.def("makeCovfieField", + py::overload_cast( + &Acts::CovfiePlugin::covfieField)); + m.def("makeCovfieField", py::overload_cast( + &Acts::CovfiePlugin::covfieField)); + m.def("makeCovfieField", + py::overload_cast&, + const Acts::Vector3&, const Acts::Vector3&>( + &Acts::CovfiePlugin::covfieField)); +} + +} // namespace Acts::Python diff --git a/Examples/Python/src/CovfieStub.cpp b/Examples/Python/src/CovfieStub.cpp new file mode 100644 index 00000000000..bb2ea09bfe1 --- /dev/null +++ b/Examples/Python/src/CovfieStub.cpp @@ -0,0 +1,13 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2021 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Plugins/Python/Utilities.hpp" + +namespace Acts::Python { +void addCovfie(Context& /* ctx */) {} +} // namespace Acts::Python diff --git a/Examples/Python/src/Detray.cpp b/Examples/Python/src/Detray.cpp index 34d23aad495..5239ad3ce36 100644 --- a/Examples/Python/src/Detray.cpp +++ b/Examples/Python/src/Detray.cpp @@ -13,13 +13,11 @@ #include #include +#include #include #include #include -#include "detray/builders/detector_builder.hpp" -#include "detray/io/frontend/detector_reader_config.hpp" - namespace py = pybind11; using namespace pybind11::literals; @@ -39,26 +37,5 @@ void addDetray(Context& ctx) { } { mex.def("writeToJson", &DetrayConverter::writeToJson); } - /** - { - /// @brief Converts an Acts::Detector to a detray::detector - mex.def("convertDetectorToDetray", - [](const Acts::GeometryContext& gctx, - const Acts::Experimental::Detector& acts_detector, - const std::string& name) -> auto - {//detector - - // Create a host memory resource - vecmem::host_memory_resource host_mr; - // Convert Acts detector to detray detector using the - detray converter function auto det_tuple = - DetrayConverter::detrayConvert(acts_detector, gctx, host_mr); - - return true; //TO DO:: cannot return tuple - - }); - } - - **/ } } // namespace Acts::Python diff --git a/Examples/Python/src/Geant4Component.cpp b/Examples/Python/src/Geant4Component.cpp index 991908c29cb..edcf28c535f 100644 --- a/Examples/Python/src/Geant4Component.cpp +++ b/Examples/Python/src/Geant4Component.cpp @@ -66,6 +66,31 @@ namespace Acts::Python { void addGeant4HepMC3(Context& ctx); } +struct ExperimentalSensitiveCandidates : public SensitiveCandidatesBase { + std::shared_ptr detector; + + /// Find the sensitive surfaces for a given position + std::vector queryPosition( + const Acts::GeometryContext& gctx, + const Acts::Vector3& position) const override { + std::vector surfaces; + // Here's the detector volume + auto volume = detector->findDetectorVolume(gctx, position); + if (volume != nullptr) { + for (const auto& surface : volume->surfaces()) { + if (surface->associatedDetectorElement() != nullptr) { + surfaces.push_back(surface); + } + } + } + return surfaces; + } + + std::vector queryAll() const override { + throw std::runtime_error("not implemented"); + } +}; + PYBIND11_MODULE(ActsPythonBindingsGeant4, mod) { py::class_>( @@ -100,6 +125,7 @@ PYBIND11_MODULE(ActsPythonBindingsGeant4, mod) { { using Config = SensitiveSurfaceMapper::Config; + using State = SensitiveSurfaceMapper::State; auto sm = py::class_>( @@ -109,6 +135,8 @@ PYBIND11_MODULE(ActsPythonBindingsGeant4, mod) { cfg, getDefaultLogger("SensitiveSurfaceMapper", level)); })); + py::class_(sm, "State").def(py::init<>()); + auto c = py::class_(sm, "Config").def(py::init<>()); ACTS_PYTHON_STRUCT_BEGIN(c, Config); ACTS_PYTHON_MEMBER(materialMappings); @@ -116,46 +144,44 @@ PYBIND11_MODULE(ActsPythonBindingsGeant4, mod) { ACTS_PYTHON_MEMBER(candidateSurfaces); ACTS_PYTHON_STRUCT_END(); - sm.def( - "create", [](const Config& cfg, Acts::Logging::Level level, - const std::shared_ptr tGeometry) { - // Set a new surface finder - Config ccfg = cfg; - ccfg.candidateSurfaces = ActsExamples::SensitiveCandidates{tGeometry}; - return std::make_shared( - ccfg, getDefaultLogger("SensitiveSurfaceMapper", level)); - }); + sm.def("create", + [](const Config& cfg, Acts::Logging::Level level, + const std::shared_ptr tGeometry) { + // Set a new surface finder + Config ccfg = cfg; + auto candidateSurfaces = + std::make_shared(); + candidateSurfaces->trackingGeometry = tGeometry; + ccfg.candidateSurfaces = candidateSurfaces; + return std::make_shared( + ccfg, getDefaultLogger("SensitiveSurfaceMapper", level)); + }); sm.def("create", [](const Config& cfg, Acts::Logging::Level level, const std::shared_ptr& detector) { // Helper struct to find the sensitive surface candidates - struct SensitiveCandidates { - std::shared_ptr detector; - - /// Find the sensitive surfaces for a given position - std::vector operator()( - const Acts::GeometryContext& gctx, - const Acts::Vector3& position) const { - std::vector surfaces; - // Here's the detector volume - auto volume = detector->findDetectorVolume(gctx, position); - if (volume != nullptr) { - for (const auto& surface : volume->surfaces()) { - if (surface->associatedDetectorElement() != nullptr) { - surfaces.push_back(surface); - } - } - } - return surfaces; - } - }; + // Set a new surface finder Config ccfg = cfg; - ccfg.candidateSurfaces = SensitiveCandidates{detector}; + auto candidateSurfaces = + std::make_shared(); + candidateSurfaces->detector = detector; + ccfg.candidateSurfaces = candidateSurfaces; return std::make_shared( ccfg, getDefaultLogger("SensitiveSurfaceMapper", level)); }); + + sm.def( + "remapSensitiveNames", + [](SensitiveSurfaceMapper& self, State& state, GeometryContext& gctx, + DetectorConstructionFactory& factory, Transform3& transform) { + return self.remapSensitiveNames( + state, gctx, factory.factorize()->Construct(), transform); + }, + "state"_a, "gctx"_a, "g4physicalVolume"_a, "motherTransform"_a); + sm.def("checkMapping", &SensitiveSurfaceMapper::checkMapping, "state"_a, + "gctx"_a, "writeMappedAsObj"_a, "writeMissingAsObj"_a); } { diff --git a/Examples/Python/src/ModuleEntry.cpp b/Examples/Python/src/ModuleEntry.cpp index 53e5d0dd3af..3cb00f50bd3 100644 --- a/Examples/Python/src/ModuleEntry.cpp +++ b/Examples/Python/src/ModuleEntry.cpp @@ -82,6 +82,7 @@ void addSvg(Context& ctx); void addObj(Context& ctx); void addOnnx(Context& ctx); void addOnnxNeuralCalibrator(Context& ctx); +void addCovfie(Context& ctx); } // namespace Acts::Python @@ -148,4 +149,5 @@ PYBIND11_MODULE(ActsPythonBindings, m) { addSvg(ctx); addOnnx(ctx); addOnnxNeuralCalibrator(ctx); + addCovfie(ctx); } diff --git a/Examples/Python/tests/helpers/__init__.py b/Examples/Python/tests/helpers/__init__.py index 67f337609c7..386c370f487 100644 --- a/Examples/Python/tests/helpers/__init__.py +++ b/Examples/Python/tests/helpers/__init__.py @@ -58,6 +58,13 @@ except ImportError: onnxEnabled = False +try: + from acts import covfie + + covfieEnabled = True +except ImportError: + covfieEnabled = False + try: import acts.examples diff --git a/Examples/Python/tests/test_covfie.py b/Examples/Python/tests/test_covfie.py new file mode 100644 index 00000000000..3c1d7ca0d3d --- /dev/null +++ b/Examples/Python/tests/test_covfie.py @@ -0,0 +1,63 @@ +import pathlib, acts, acts.examples +import pytest + +from helpers import covfieEnabled + + +@pytest.mark.skipif(not covfieEnabled, reason="Covfie plugin not available") +def test_constant_field_conversion(): + from acts import covfie + + v = acts.Vector3(1, 2, 3) + af = acts.ConstantBField(v) + cf = covfie.makeCovfieField(af) + view = covfie.toView(cf) + points = [(0, 0, 1), (1, 1, 1), (1, 0, 2)] + for x, y, z in points: + assert view.at(x, y, z) == [1, 2, 3] + + +@pytest.mark.skipif(not covfieEnabled, reason="Covfie plugin not available") +def test_root_field_conversion(): + from acts import covfie + + current_file_path = pathlib.Path(__file__).resolve().parent + p = ( + current_file_path.parent.parent.parent + / "thirdparty" + / "OpenDataDetector" + / "data" + / "odd-bfield.root" + ) + + af = acts.examples.MagneticFieldMapXyz(str(p)) + bc = acts.MagneticFieldContext() + fc = af.makeCache(bc) + + cf = covfie.makeCovfieField(af) + view = covfie.toView(cf) + points = [ + (9300.0, 4700.0, 11200.0), + (9999.0, 9999.0, 14300.0), + (-2900.0, -4700.0, 5200.0), + (-2900.0, -4800.0, 9100.0), + (-2900.0, -5200.0, -8800.0), + (-4400.0, 4800.0, -12700.0), + (-6600.0, 1900.0, 7700.0), + (-9700.0, -900.0, 12700.0), + (-9999.0, -9999.0, -13000.0), + (9999.0, 0, 14900.0), + ] + + error_margin_half_width = 0.0001 + for x, y, z in points: + val = af.getField(acts.Vector3(x, y, z), fc) + Bx1, By1, Bz1 = val[0], val[1], val[2] + + Bx2, By2, Bz2 = tuple(view.at(x, y, z)) + + assert ( + abs(Bx1 - Bx2) < error_margin_half_width + and abs(By1 - By2) < error_margin_half_width + and abs(Bz1 - Bz2) < error_margin_half_width + ) diff --git a/Examples/Scripts/Python/geant4.py b/Examples/Scripts/Python/geant4.py index 85ccbbfa25e..42cb5db95bb 100755 --- a/Examples/Scripts/Python/geant4.py +++ b/Examples/Scripts/Python/geant4.py @@ -17,6 +17,7 @@ def runGeant4( field, outputDir, materialMappings=["Silicon"], + volumeMappings=[], s: acts.examples.Sequencer = None, ): s = s or acts.examples.Sequencer(events=100, numThreads=1) @@ -37,6 +38,7 @@ def runGeant4( outputDirRoot=outputDir, rnd=rnd, materialMappings=materialMappings, + volumeMappings=volumeMappings, ) return s diff --git a/Plugins/CMakeLists.txt b/Plugins/CMakeLists.txt index b0d006a61c1..ab6d375a27b 100644 --- a/Plugins/CMakeLists.txt +++ b/Plugins/CMakeLists.txt @@ -11,6 +11,7 @@ add_component_if(Legacy PluginLegacy ACTS_BUILD_PLUGIN_LEGACY) add_component_if(Onnx PluginOnnx ACTS_BUILD_PLUGIN_ONNX) add_component_if(ExaTrkX PluginExaTrkX ACTS_BUILD_PLUGIN_EXATRKX) add_component_if(Detray PluginDetray ACTS_BUILD_PLUGIN_TRACCC) +add_component_if(Covfie PluginCovfie ACTS_BUILD_PLUGIN_TRACCC) # dependent plugins. depend either on a independent plugins or on one another add_component_if(TGeo PluginTGeo ACTS_BUILD_PLUGIN_TGEO) diff --git a/Plugins/Covfie/CMakeLists.txt b/Plugins/Covfie/CMakeLists.txt new file mode 100644 index 00000000000..0245904b035 --- /dev/null +++ b/Plugins/Covfie/CMakeLists.txt @@ -0,0 +1,16 @@ +add_library(ActsPluginCovfie SHARED src/FieldConversion.cpp) + +target_include_directories( + ActsPluginCovfie + PUBLIC + $ + $ +) +target_link_libraries(ActsPluginCovfie PUBLIC ActsCore covfie::core) + +install( + TARGETS ActsPluginCovfie + EXPORT ActsPluginCovfieTargets + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} +) +install(DIRECTORY include/Acts DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) diff --git a/Plugins/Covfie/include/Acts/Plugins/Covfie/FieldConversion.hpp b/Plugins/Covfie/include/Acts/Plugins/Covfie/FieldConversion.hpp new file mode 100644 index 00000000000..cfb1d6b516c --- /dev/null +++ b/Plugins/Covfie/include/Acts/Plugins/Covfie/FieldConversion.hpp @@ -0,0 +1,64 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// acts includes +#include "Acts/MagneticField/BFieldMapUtils.hpp" +#include "Acts/MagneticField/ConstantBField.hpp" +#include "Acts/MagneticField/MagneticFieldProvider.hpp" + +namespace Acts::CovfiePlugin { + +using BuilderBackend = + covfie::backend::strided>; + +using InterpolatedField = covfie::field>>>; + +using ConstantField = covfie::field< + covfie::backend::constant>; + +/// @brief Creates a covfie field from an interpolated magnetic field. +/// @param magneticField The acts interpolated magnetic field. +/// @return An affine linear strided covfie field. +InterpolatedField covfieField( + const Acts::InterpolatedMagneticField& magneticField); + +/// @brief Creates a covfie field from a constant B field. +/// @param magneticField The acts constant magnetic field. +/// @return A constant covfie field. +ConstantField covfieField(const Acts::ConstantBField& magneticField); + +/// @brief Creates a covfie field from a magnetic field provider by sampling it. +/// The field must be defined within min (inclusive) and max (inclusive). +/// @param magneticField The acts magnetic field provider. +/// @param cache The acts cache. +/// @param nPoints 3D array of containing the number of bins for each axis. +/// @param min (min_x, min_y, min_z) +/// @param max (max_x, max_y, max_z) +/// @return An affine linear strided covfie field. +InterpolatedField covfieField(const Acts::MagneticFieldProvider& magneticField, + Acts::MagneticFieldProvider::Cache& cache, + const std::array& nPoints, + const Acts::Vector3& min, + const Acts::Vector3& max); + +} // namespace Acts::CovfiePlugin diff --git a/Plugins/Covfie/src/FieldConversion.cpp b/Plugins/Covfie/src/FieldConversion.cpp new file mode 100644 index 00000000000..6f182f25885 --- /dev/null +++ b/Plugins/Covfie/src/FieldConversion.cpp @@ -0,0 +1,197 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Plugins/Covfie/FieldConversion.hpp" + +#include "Acts/Utilities/Axis.hpp" +#include "Acts/Utilities/AxisFwd.hpp" +#include "Acts/Utilities/Grid.hpp" + +#include +#include +#include + +namespace Acts::CovfiePlugin { + +/// @brief Creates a strided covfie field that stores the values of the +/// magnetic field in the volume given by min and max using a fixed sample +/// spacing (determined by nPoints). +/// +/// @param magneticField The acts magnetic field. +/// @param cache The acts cache. +/// @param nPoints 3D array of containing the number of bins for each axis. +/// @param min (min_x, min_y, min_z) +/// @param max (max_x, max_y, max_z) +/// @return A strided covfie field. +template +auto newBuilder(const magnetic_field_t& magneticField, + typename magnetic_field_t::Cache& cache, + const std::array& nPoints, + const Acts::Vector3& min, const Acts::Vector3& max) { + using Field = covfie::field; + + // Hack to avoid the fact that the domain of ACTS magnetic fields is defined + // as a half-open interval. Has the potential to introduce very minor + // floating point errors, but no easy way to fix this right now. + // TODO: Fix the aforementioned problem. + std::vector maxima = { + std::nexttoward(max[0], -std::numeric_limits::infinity()), + std::nexttoward(max[1], -std::numeric_limits::infinity()), + std::nexttoward(max[1], -std::numeric_limits::infinity()), + }; + + Field field(covfie::make_parameter_pack( + Field::backend_t::configuration_t{nPoints[0], nPoints[1], nPoints[2]})); + + Field::view_t view(field); + + std::array sampleSpacing = { + (max.x() - min.x()) / (nPoints[0] - 1), + (max.y() - min.y()) / (nPoints[1] - 1), + (max.z() - min.z()) / (nPoints[2] - 1)}; + + for (std::size_t x = 0; x < nPoints[0]; x++) { + for (std::size_t y = 0; y < nPoints[1]; y++) { + for (std::size_t z = 0; z < nPoints[2]; z++) { + Acts::Vector3 position{ + std::min(x * sampleSpacing[0] + min[0], maxima[0]), + std::min(y * sampleSpacing[1] + min[1], maxima[1]), + std::min(z * sampleSpacing[2] + min[2], maxima[2])}; + + Field::view_t::output_t& p = view.at(x, y, z); + Result result = magneticField.getField(position, cache); + + if (!result.ok()) { + throw std::runtime_error("Field lookup failed!"); + } + + Acts::Vector3 rv = *result; + p[0] = static_cast(rv[0]); + p[1] = static_cast(rv[1]); + p[2] = static_cast(rv[2]); + } + } + } + + return field; +} + +/// @brief Generate the affine covfie configuration (scaling and rotation) +/// given the size of the field (min and max) +/// +/// @param nPoints 3D array of containing the number of bins for each axis. +/// @param min (min_x, min_y, min_z) +/// @param max (max_x, max_y, max_z) +/// @return The affine field configuration. +template +typename backend_t::configuration_t affineConfiguration( + const std::array& nPoints, const Acts::Vector3& min, + const Acts::Vector3& max) { + auto scaling = covfie::algebra::affine<3>::scaling( + static_cast((nPoints[0] - 1) / (max[0] - min[0])), + static_cast((nPoints[1] - 1) / (max[1] - min[1])), + static_cast((nPoints[2] - 1) / (max[2] - min[2]))); + + auto translation = covfie::algebra::affine<3>::translation( + static_cast(-min[0]), static_cast(-min[1]), + static_cast(-min[2])); + + return {scaling * translation}; +} + +/// @brief Uses std::nextafter to generates a clamp backend +/// configuration where arguments min and max hold floating point values. +/// @param min (min_x, min_y, min_z) +/// @param max (max_x, max_y, max_z) +/// @return The clamp field configuration. +template +typename backend_t::configuration_t clampConfigurationFloat( + const Acts::Vector3& min, const Acts::Vector3& max) { + return {{std::nextafter(static_cast(min[0]), + std::numeric_limits::infinity()), + std::nextafter(static_cast(min[1]), + std::numeric_limits::infinity()), + std::nextafter(static_cast(min[2]), + std::numeric_limits::infinity())}, + {std::nextafter(static_cast(max[0]), + -std::numeric_limits::infinity()), + std::nextafter(static_cast(max[1]), + -std::numeric_limits::infinity()), + std::nextafter(static_cast(max[2]), + -std::numeric_limits::infinity())}}; +} + +/// @brief Creates a covfie field from a generic magnetic field. +/// @param magneticField The generic magnetic field. +/// @param cache The cache. +/// @param nPoints 3D array of containing the number of bins for each axis. +/// @param min (min_x, min_y, min_z) +/// @param max (max_x, max_y, max_z) +/// @return A clamp affine linear strided covfie field. +template +InterpolatedField covfieFieldLinear(const magnetic_field_t& magneticField, + typename magnetic_field_t::Cache& cache, + const std::array& nPoints, + const Acts::Vector3& min, + const Acts::Vector3& max) { + auto builder = newBuilder(magneticField, cache, nPoints, min, max); + InterpolatedField field(covfie::make_parameter_pack( + clampConfigurationFloat(min, max), + affineConfiguration(nPoints, min, + max), + InterpolatedField::backend_t::backend_t::backend_t::configuration_t{}, + builder.backend())); + + return field; +} + +/// @brief Creates a covfie field from a magnetic field provider by sampling it. +/// @param magneticField The acts magnetic field provider. +/// @param cache The acts cache. +/// @param nPoints 3D array of containing the number of bins for each axis. +/// @param min (min_x, min_y, min_z) +/// @param max (max_x, max_y, max_z) +/// @return A clamp affine linear strided covfie field. +InterpolatedField covfieField(const Acts::MagneticFieldProvider& magneticField, + Acts::MagneticFieldProvider::Cache& cache, + const std::array& nPoints, + const Acts::Vector3& min, + const Acts::Vector3& max) { + return covfieFieldLinear(magneticField, cache, nPoints, min, max); +} + +/// @brief Creates a covfie field from an interpolated magnetic field. +/// @param magneticField The acts interpolated magnetic field. +/// @return A clamp affine linear strided covfie field. +InterpolatedField covfieField( + const Acts::InterpolatedMagneticField& magneticField) { + Acts::MagneticFieldContext ctx; + auto cache = magneticField.makeCache(ctx); + const std::vector& old_min = magneticField.getMin(); + const std::vector& old_max = magneticField.getMax(); + const std::vector& old_nbins = magneticField.getNBins(); + Acts::Vector3 min{old_min.at(0), old_min.at(1), old_min.at(2)}; + Acts::Vector3 max{old_max.at(0), old_max.at(1), old_max.at(2)}; + std::array nPoints{old_nbins.at(0), old_nbins.at(1), + old_nbins.at(2)}; + return covfieFieldLinear(magneticField, cache, nPoints, min, max); +} + +/// @brief Creates a covfie field from a constant B field. +/// @param magneticField The acts constant magnetic field. +/// @return A constant covfie field. +ConstantField covfieField(const Acts::ConstantBField& magneticField) { + auto B = magneticField.getField(); + ConstantField field( + covfie::make_parameter_pack(ConstantField::backend_t::configuration_t{ + static_cast(B[0]), static_cast(B[1]), + static_cast(B[2])})); + return field; +} + +} // namespace Acts::CovfiePlugin diff --git a/Plugins/Detray/CMakeLists.txt b/Plugins/Detray/CMakeLists.txt index 4dba2bf7c5d..79d56b37cf0 100644 --- a/Plugins/Detray/CMakeLists.txt +++ b/Plugins/Detray/CMakeLists.txt @@ -1,4 +1,11 @@ -add_library(ActsPluginDetray SHARED src/DetrayConverter.cpp) +add_library( + ActsPluginDetray + SHARED + src/DetrayConversionUtils.cpp + src/DetrayConverter.cpp + src/DetrayGeometryConverter.cpp + src/DetrayMaterialConverter.cpp +) add_dependencies(ActsPluginDetray detray::core covfie::core vecmem::core) diff --git a/Plugins/Detray/include/Acts/Plugins/Detray/DetrayConversionUtils.hpp b/Plugins/Detray/include/Acts/Plugins/Detray/DetrayConversionUtils.hpp new file mode 100644 index 00000000000..6e5a6788a63 --- /dev/null +++ b/Plugins/Detray/include/Acts/Plugins/Detray/DetrayConversionUtils.hpp @@ -0,0 +1,78 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Utilities/BinningData.hpp" +#include "Acts/Utilities/BinningType.hpp" + +#include + +#include +#include +#include + +namespace Acts { + +using DetrayDetector = detray::detector; + +namespace DetrayConversionUtils { + +/// Detray conversion options +struct Options { + /// Option to switch on/off the material conversion + bool convertMaterial = true; + /// Option to switch on/off the surface grid conversin + bool convertSurfaceGrids = true; + /// Option to switch on/off the export to json + bool writeToJson = false; +}; + +/// Detray conversion cache object +/// +/// This object is used to synchronize link information between the +/// different converters (geometry, material, surface grids) +struct GeometryIdCache { + /// This is a multimap to pass volume local surface link information + /// The portal splitting requires a multimap implementation here + std::multimap localSurfaceLinks; + /// This is a map to pass on volume link information + std::map volumeLinks; +}; + +/// Convert the binning option +/// +/// @param bOption the binning option +/// +/// @return a detray binning option +detray::axis::bounds convertBinningOption(BinningOption bOption); + +/// Convert the binning value +/// +/// @param bValue the binning value +/// +/// @return a detray binning value +detray::axis::label convertBinningValue(BinningValue bValue); + +/// Convert the binning type +/// +/// @param bType the binning type +/// +/// @return a detray binning type +detray::axis::binning convertBinningType(BinningType bType); + +/// Convert the binning data to an axis +/// +/// @param bData the binning data to be converted +/// +/// @return a detray axis payload +detray::io::axis_payload convertBinningData(const BinningData& bData); + +} // namespace DetrayConversionUtils +} // namespace Acts diff --git a/Plugins/Detray/include/Acts/Plugins/Detray/DetrayConverter.hpp b/Plugins/Detray/include/Acts/Plugins/Detray/DetrayConverter.hpp index ff332d6d2a6..ffb799cc29b 100644 --- a/Plugins/Detray/include/Acts/Plugins/Detray/DetrayConverter.hpp +++ b/Plugins/Detray/include/Acts/Plugins/Detray/DetrayConverter.hpp @@ -8,137 +8,99 @@ #pragma once -#include "Acts/Definitions/Algebra.hpp" #include "Acts/Detector/Detector.hpp" -#include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Geometry/VolumeBounds.hpp" +#include "Acts/Plugins/Detray/DetrayConversionUtils.hpp" +#include "Acts/Plugins/Detray/DetrayGeometryConverter.hpp" +#include "Acts/Plugins/Detray/DetrayMaterialConverter.hpp" +#include "Acts/Utilities/Logger.hpp" -#include +#include -#include "detray/builders/detector_builder.hpp" -#include "detray/core/detector.hpp" -#include "detray/definitions/geometry.hpp" -#include "detray/io/common/geometry_reader.hpp" -#include "detray/io/frontend/detector_writer.hpp" -#include "detray/io/frontend/payloads.hpp" -#include "detray/utils/consistency_checker.hpp" +#include +#include +#include namespace Acts { -class Surface; -class SurfaceBounds; - -namespace Experimental { -class DetectorVolume; -class Portal; -} // namespace Experimental - -namespace DetrayConverter { - -using DetrayDetector = detray::detector; - -/// Write the detector to json output -/// -/// @param dDetector is the detray detector (converted) -/// @param names a name map for the detector volumes -/// @param writer_cfg the writer configuration -void writeToJson(const DetrayDetector& dDetector, - const typename DetrayDetector::name_map& names = {}, - detray::io::detector_writer_config writer_cfg = {}); - -/// Conversion method for transform objects to detray::transform payloads -/// -/// @param t the transform to be converted -/// -/// @return the transform_payload(translation, rotation) -detray::io::transform_payload convertTransform(const Transform3& t); - -/// Conversion method for surface bounds to detray::mask payloads -/// -/// @param bounds the bounds object to be converted -/// @param portal the flag for conversion into detray portal format -/// -/// @return the mask_payload representing the bounds -detray::io::mask_payload convertMask(const SurfaceBounds& bounds, - bool portal = false); - -/// Conversion method for surface objects to detray::surface payloads -/// -/// @param gctx the geometry context -/// @param surface the surface to be converted -/// @param portal the flag for conversion into detray portal format -/// -/// @return the surface_payload for portals and volumes by @param Surface acts object -detray::io::surface_payload convertSurface(const GeometryContext& gctx, - const Surface& surface, - bool portal = false); -/// Conversion method for Portal object to detray::portal payloads -/// -/// @param gctx the geometry context -/// @param portal the portal to be converted -/// @param ip the portal index -/// @param volume the volume to which the portal belongs -/// @param orientedSurfaces the oriented surfaces of the portal -/// @param detectorVolumes the detector volumes for the link lookup -/// -/// @note due to portal splitting this can add up in N portals for one initial one -/// -/// @brief convert the acts portal to detray surface payload and populate the payload -std::vector convertPortal( - const GeometryContext& gctx, const Experimental::Portal& portal, - std::size_t ip, const Experimental::DetectorVolume& volume, - const std::vector& orientedSurfaces, - const std::vector& detectorVolumes); - -/// Conversion method for volume objects to detray::volume payloads -/// -/// @param gctx the geometry context -/// @param volume the volume to be converted -/// @param detectorVolumes the detector volumes for the link lookup -/// -/// @return the volume_payload for portals and volumes by @param volume acts object -detray::io::volume_payload convertVolume( - const GeometryContext& gctx, const Experimental::DetectorVolume& volume, - const std::vector& detectorVolumes); - -/// Conversion method for (common) header payload -/// -/// @param detector is the detector to be converted -/// -/// @return a geometry header payload -detray::io::geo_header_payload convertHead( - const Acts::Experimental::Detector& detector); - -/// Convert an Acts::Experimental::Detector to a detray::detector object -/// -/// @param gctx the geometry context -/// @param detector the detector to be converted -/// @param mr the memory resource to be used -/// -/// @returns a detector of requested return type -template -std::tuple convertDetector( - const Acts::GeometryContext& gctx, - const Acts::Experimental::Detector& detector, vecmem::memory_resource& mr) { - detray::io::detector_payload detectorPayload; - for (const auto volume : detector.volumes()) { - detectorPayload.volumes.push_back( - convertVolume(gctx, *volume, detector.volumes())); +using namespace Experimental; + +class DetrayConverter { + public: + /// Constructor with logger + DetrayConverter(std::unique_ptr logger = + getDefaultLogger("DetrayConverter", Logging::INFO)); + + /// Convert an Acts::Experimental::Detector to a detray::detector object + /// + /// @param gctx the geometry context + /// @param detector the detector to be converted + /// @param mr the memory resource to be used + /// @param options the conversion options + /// + /// @returns a detector of requested return type + template + detector_t convert( + const GeometryContext& gctx, const Detector& detector, + vecmem::memory_resource& mr, + [[maybe_unused]] const DetrayConversionUtils::Options& options = {}) { + // The building cache object + DetrayConversionUtils::GeometryIdCache geoIdCache; + + typename detector_t::name_map names = {{0u, detector.name()}}; + + // build detector + detray::detector_builder detectorBuilder{}; + // (1) geometry + detray::io::detector_payload detectorPayload = + DetrayGeometryConverter::convertDetector(geoIdCache, gctx, detector, + logger()); + detray::io::geometry_reader::convert(detectorBuilder, names, + detectorPayload); + // (2) material + if constexpr (detray::detail::has_material_grids_v) { + if (options.convertMaterial) { + detray::io::detector_grids_payload + materialPayload = + DetrayMaterialConverter::convertSurfaceMaterialGrids( + geoIdCache, detector, logger()); + detray::io::material_map_reader<>::convert( + detectorBuilder, names, materialPayload); + } + } + + detector_t detrayDetector(detectorBuilder.build(mr)); + + // Checks and print + detray::detail::check_consistency(detrayDetector); + + // If configured, write the detector to json + if (options.writeToJson) { + // Create a writer configuration and write it out + detray::io::detector_writer_config writerConfig{}; + writerConfig.m_write_material = options.convertMaterial; + writerConfig.m_write_grids = options.convertSurfaceGrids; + writeToJson(detrayDetector, names, writerConfig); + } + + return detrayDetector; } - typename detector_t::name_map names = {{0u, detector.name()}}; - // build detector - detray::detector_builder detectorBuilder{}; - detray::io::geometry_reader::convert(detectorBuilder, names, - detectorPayload); - detector_t detrayDetector(detectorBuilder.build(mr)); + /// Write the detector to json output + /// + /// @param dDetector is the detray detector (converted) + /// @param names a name map for the detector volumes + /// @param writer_cfg the writer configuration + static void writeToJson(const DetrayDetector& dDetector, + const typename DetrayDetector::name_map& names = {}, + detray::io::detector_writer_config writer_cfg = {}); - // checks and print - detray::detail::check_consistency(detrayDetector); - converterPrint(detrayDetector, names); + private: + /// The logger instance + std::unique_ptr m_logger = nullptr; - return {std::move(detrayDetector), mr}; -} + // Return the logging instance + const Acts::Logger& logger() const { return *m_logger; } +}; -} // namespace DetrayConverter } // namespace Acts diff --git a/Plugins/Detray/include/Acts/Plugins/Detray/DetrayGeometryConverter.hpp b/Plugins/Detray/include/Acts/Plugins/Detray/DetrayGeometryConverter.hpp new file mode 100644 index 00000000000..7268b30be7d --- /dev/null +++ b/Plugins/Detray/include/Acts/Plugins/Detray/DetrayGeometryConverter.hpp @@ -0,0 +1,112 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Detector/Detector.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/VolumeBounds.hpp" +#include "Acts/Plugins/Detray/DetrayConversionUtils.hpp" +#include "Acts/Utilities/Logger.hpp" + +#include + +#include +#include +#include +#include +#include + +namespace Acts { + +class Surface; +class SurfaceBounds; + +namespace Experimental { +class DetectorVolume; +class Detector; +class Portal; +} // namespace Experimental + +namespace DetrayGeometryConverter { + +/// Conversion method for transform objects to detray::transform payloads +/// +/// @param t the transform to be converted +/// +/// @return the transform_payload(translation, rotation) +detray::io::transform_payload convertTransform(const Transform3& t); + +/// Conversion method for surface bounds to detray::mask payloads +/// +/// @param bounds the bounds object to be converted +/// @param portal the flag for conversion into detray portal format +/// +/// @return the mask_payload representing the bounds +detray::io::mask_payload convertMask(const SurfaceBounds& bounds, + bool portal = false); + +/// Conversion method for surface objects to detray::surface payloads +/// +/// @param gctx the geometry context +/// @param surface the surface to be converted +/// @param portal the flag for conversion into detray portal format +/// +/// @return the surface_payload for portals and volumes by @param Surface acts object +detray::io::surface_payload convertSurface(const GeometryContext& gctx, + const Surface& surface, + bool portal = false); +/// Conversion method for Portal object to detray::portal payloads +/// +/// @param gctx the geometry context +/// @param portal the portal to be converted +/// @param ip the portal index +/// @param volume the volume to which the portal belongs +/// @param orientedSurfaces the oriented surfaces of the portal +/// @param detectorVolumes the detector volumes for the link lookup +/// +/// @note due to portal splitting this can add up in N portals for one initial one +/// +/// @brief convert the acts portal to detray surface payload and populate the payload +std::vector convertPortal( + const GeometryContext& gctx, const Experimental::Portal& portal, + std::size_t ip, const Experimental::DetectorVolume& volume, + const std::vector& orientedSurfaces, + const std::vector& detectorVolumes); + +/// Conversion method for volume objects to detray::volume payloads +/// +/// @param geoIdCache [in, out] object +/// @param gctx the geometry context +/// @param volume the volume to be converted +/// @param detectorVolumes the detector volumes for the link lookup +/// @param logger the logger object for screen output +/// +/// @return the volume_payload for portals and volumes by @param volume acts object +detray::io::volume_payload convertVolume( + DetrayConversionUtils::GeometryIdCache& geoIdCache, + const GeometryContext& gctx, const Experimental::DetectorVolume& volume, + const std::vector& detectorVolumes, + const Acts::Logger& logger); + +/// Conversion method for detector objects to detray::detector payload +/// +/// @param geoIdCache [in, out] object +/// @param gctx the geometry context +/// @param detector the detector to be converted +/// @param logger the logger object for screen output +/// +/// @return the detector_payload for portals and volumes by @param detector acts object +detray::io::detector_payload convertDetector( + DetrayConversionUtils::GeometryIdCache& geoIdCache, + const GeometryContext& gctx, const Experimental::Detector& detector, + const Acts::Logger& logger); + +} // namespace DetrayGeometryConverter +} // namespace Acts diff --git a/Plugins/Detray/include/Acts/Plugins/Detray/DetrayMaterialConverter.hpp b/Plugins/Detray/include/Acts/Plugins/Detray/DetrayMaterialConverter.hpp new file mode 100644 index 00000000000..d7edf20b607 --- /dev/null +++ b/Plugins/Detray/include/Acts/Plugins/Detray/DetrayMaterialConverter.hpp @@ -0,0 +1,61 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Material/ISurfaceMaterial.hpp" +#include "Acts/Material/MaterialSlab.hpp" +#include "Acts/Plugins/Detray/DetrayConversionUtils.hpp" +#include "Acts/Utilities/BinningData.hpp" +#include "Acts/Utilities/Logger.hpp" + +#include "detray/io/frontend/payloads.hpp" + +namespace Acts { + +namespace Experimental { +class Detector; +} + +namespace DetrayMaterialConverter { + +/// Conversion method for material slab +/// +/// @param materialSlab the material slab object +/// +/// @return a material slab payload +detray::io::material_slab_payload convertMaterialSlab( + const MaterialSlab& materialSlab); + +/// Conversion method for surface material objects +/// +/// @param detector the detector object +/// @param logger the logger object for screen output +/// +/// @return a surface material +detray::io::grid_payload +convertSurfaceMaterial(const ISurfaceMaterial& material, + const Acts::Logger& logger); + +/// Conversion method for material grids +/// +/// @param geoIdCache object to have the link association from the geometry building +/// @param detector the detector object +/// @param logger the logger object for screen output +/// +/// @return the volume_payload for portals and volumes by @param volume acts object +detray::io::detector_grids_payload +convertSurfaceMaterialGrids( + const DetrayConversionUtils::GeometryIdCache& geoIdCache, + const Experimental::Detector& detector, const Logger& logger); + +} // namespace DetrayMaterialConverter + +} // namespace Acts diff --git a/Plugins/Detray/src/DetrayConversionUtils.cpp b/Plugins/Detray/src/DetrayConversionUtils.cpp new file mode 100644 index 00000000000..6c26781d8bd --- /dev/null +++ b/Plugins/Detray/src/DetrayConversionUtils.cpp @@ -0,0 +1,82 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Plugins/Detray/DetrayConversionUtils.hpp" + +detray::axis::label Acts::DetrayConversionUtils::convertBinningValue( + BinningValue bValue) { + switch (bValue) { + case BinningValue::binX: + return detray::axis::label::e_x; + case BinningValue::binY: + return detray::axis::label::e_y; + case BinningValue::binZ: + return detray::axis::label::e_z; + case BinningValue::binR: + return detray::axis::label::e_r; + case BinningValue::binPhi: + return detray::axis::label::e_phi; + case BinningValue::binRPhi: + return detray::axis::label::e_rphi; + default: + throw std::invalid_argument( + "DetrayMaterialConverter: unknown binning value detected."); + } +} + +detray::axis::bounds Acts::DetrayConversionUtils::convertBinningOption( + BinningOption bOption) { + // That's a bit of a mind bender, but the conversion is correct + // closed -> axis are closed, i.e. circular + // open -> axis are not closed, but the range is closed (no overflow bin) -> + // closed + switch (bOption) { + case BinningOption::closed: + return detray::axis::bounds::e_circular; + case BinningOption::open: + return detray::axis::bounds::e_closed; + default: + throw std::invalid_argument( + "DetrayMaterialConverter: unknown binning option detected."); + } +} + +detray::axis::binning Acts::DetrayConversionUtils::convertBinningType( + BinningType bType) { + switch (bType) { + case BinningType::equidistant: + return detray::axis::binning::e_regular; + case BinningType::arbitrary: + return detray::axis::binning::e_irregular; + default: + throw std::invalid_argument( + "DetrayMaterialConverter: unknown binning type detected."); + } +} + +detray::io::axis_payload Acts::DetrayConversionUtils::convertBinningData( + const BinningData& bData) { + detray::io::axis_payload axis; + + axis.bins = bData.bins(); + // Set the binning type + axis.binning = convertBinningType(bData.type); + // Set the binning option + axis.bounds = convertBinningOption(bData.option); + // Set the binning value + axis.label = convertBinningValue(bData.binvalue); + // Set the binning range + axis.edges = {}; + if (bData.type == BinningType::equidistant) { + axis.edges = {bData.min, bData.max}; + } else { + axis.edges.insert(axis.edges.end(), bData.boundaries().begin(), + bData.boundaries().end()); + } + return axis; +} diff --git a/Plugins/Detray/src/DetrayConverter.cpp b/Plugins/Detray/src/DetrayConverter.cpp index 4f2ac983ad0..3b72b998ab5 100644 --- a/Plugins/Detray/src/DetrayConverter.cpp +++ b/Plugins/Detray/src/DetrayConverter.cpp @@ -8,43 +8,10 @@ #include "Acts/Plugins/Detray/DetrayConverter.hpp" -#include "Acts/Detector/Detector.hpp" -#include "Acts/Detector/DetectorVolume.hpp" -#include "Acts/Detector/Portal.hpp" -#include "Acts/Navigation/PortalNavigation.hpp" -#include "Acts/Plugins/Json/DetrayJsonHelper.hpp" -#include "Acts/Surfaces/CylinderBounds.hpp" -#include "Acts/Surfaces/CylinderSurface.hpp" -#include "Acts/Surfaces/DiscSurface.hpp" -#include "Acts/Surfaces/RadialBounds.hpp" -#include "Acts/Surfaces/RegularSurface.hpp" -#include "Acts/Surfaces/Surface.hpp" -#include "Acts/Surfaces/SurfaceBounds.hpp" +Acts::DetrayConverter::DetrayConverter( + std::unique_ptr logger) + : m_logger(std::move(logger)) {} -#include "detray/io/frontend/detector_writer.hpp" - -using namespace detray; - -namespace { - -/// Find the position of the volume to point to -/// -/// @param volume the volume to find -/// @param the collection of volumes -/// -/// @note return -1 if not found, to be interpreted by the caller -int findVolume( - const Acts::Experimental::DetectorVolume* volume, - const std::vector& volumes) { - auto candidate = std::find(volumes.begin(), volumes.end(), volume); - if (candidate != volumes.end()) { - return std::distance(volumes.begin(), candidate); - } - return -1; -} -} // namespace - -/// detray geometry writer function, debug purposes void Acts::DetrayConverter::writeToJson( const DetrayDetector& dDetector, const typename DetrayDetector::name_map& names, @@ -52,265 +19,3 @@ void Acts::DetrayConverter::writeToJson( writer_cfg.format(detray::io::format::json); detray::io::write_detector(dDetector, names, writer_cfg); } - -detray::io::transform_payload Acts::DetrayConverter::convertTransform( - const Transform3& t) { - detray::io::transform_payload tfPayload; - auto translation = t.translation(); - tfPayload.tr = {translation.x(), translation.y(), translation.z()}; - - const auto rotation = t.rotation(); - tfPayload.rot = {rotation(0, 0), rotation(0, 1), rotation(0, 2), - rotation(1, 0), rotation(1, 1), rotation(1, 2), - rotation(2, 0), rotation(2, 1), rotation(2, 2)}; - return tfPayload; -} - -detray::io::mask_payload Acts::DetrayConverter::convertMask( - const Acts::SurfaceBounds& bounds, bool portal) { - detray::io::mask_payload maskPayload; - auto [shape, boundaries] = DetrayJsonHelper::maskFromBounds(bounds, portal); - maskPayload.shape = static_cast(shape); - maskPayload.boundaries = static_cast>(boundaries); - // default maskPayload.volume_link - - return maskPayload; -} - -detray::io::surface_payload Acts::DetrayConverter::convertSurface( - const Acts::GeometryContext& gctx, const Surface& surface, bool portal) { - detray::io::surface_payload surfacePayload; - - surfacePayload.transform = convertTransform(surface.transform(gctx)); - surfacePayload.source = surface.geometryId().value(); - surfacePayload.barcode = std::nullopt; - surfacePayload.type = static_cast( - portal ? surface_id::e_portal - : (surface.geometryId().sensitive() > 0 - ? detray::surface_id::e_sensitive - : detray::surface_id::e_passive)); - surfacePayload.mask = convertMask(surface.bounds()); - return surfacePayload; -} - -std::vector Acts::DetrayConverter::convertPortal( - const GeometryContext& gctx, const Experimental::Portal& portal, - std::size_t ip, const Experimental::DetectorVolume& volume, - const std::vector& orientedSurfaces, - const std::vector& detectorVolumes) { - std::vector portals{}; - - const RegularSurface& surface = portal.surface(); - const auto& volumeLinks = portal.portalNavigation(); - - // First assumption for outside link (along direction) - std::size_t outside = 1u; - - // Find out if you need to take the outside or inside volume - // for planar surfaces that's easy - if (surface.type() != Acts::Surface::SurfaceType::Cylinder) { - // Get the two volume center - const auto volumeCenter = volume.transform(gctx).translation(); - const auto surfaceCenter = surface.center(gctx); - const auto surfaceNormal = surface.normal(gctx, surfaceCenter); - // Get the direction from the volume to the surface, correct link - const auto volumeToSurface = surfaceCenter - volumeCenter; - if (volumeToSurface.dot(surfaceNormal) < 0.) { - outside = 0u; - } - } else { - // This is a cylinder portal, inner cover reverses the normal - if (ip == 3u) { - outside = 0u; - } - } - - const auto& outsideLink = volumeLinks[outside]; - // Grab the corresponding volume link - // If it is a single link, we are done - const auto* instance = outsideLink.instance(); - // Single link cast - auto singleLink = - dynamic_cast( - instance); - - auto [surfaceAdjusted, insidePointer] = orientedSurfaces[ip]; - - // Single link detected - just write it out, we use the oriented surface - // in order to make sure the size is adjusted - if (singleLink != nullptr) { - // Single link can be written out - std::size_t vLink = findVolume(singleLink->object(), detectorVolumes); - auto portalPayload = convertSurface(gctx, *surfaceAdjusted, true); - portalPayload.mask.volume_link.link = vLink; - portals.push_back(portalPayload); - } else { - // Multi link detected - 1D - auto multiLink1D = - dynamic_cast( - instance); - if (multiLink1D != nullptr) { - // Resolve the multi link 1D - auto boundaries = - multiLink1D->indexedUpdater.grid.axes()[0u]->getBinEdges(); - const auto& cast = multiLink1D->indexedUpdater.casts[0u]; - const auto& transform = multiLink1D->indexedUpdater.transform; - const auto& volumes = multiLink1D->indexedUpdater.extractor.dVolumes; - - // Apply the correction from local to global boundaries - ActsScalar gCorr = VectorHelpers::cast(transform.translation(), cast); - std::for_each(boundaries.begin(), boundaries.end(), - [&gCorr](ActsScalar& b) { b -= gCorr; }); - - // Get the volume indices - auto surfaceType = surfaceAdjusted->type(); - std::vector vIndices = {}; - for (const auto& v : volumes) { - vIndices.push_back(findVolume(v, detectorVolumes)); - } - - // Pick the surface dimension - via poly - std::array clipRange = {0., 0.}; - std::vector boundValues = surfaceAdjusted->bounds().values(); - if (surfaceType == Surface::SurfaceType::Cylinder && - cast == BinningValue::binZ) { - ActsScalar zPosition = surfaceAdjusted->center(gctx).z(); - clipRange = { - zPosition - boundValues[CylinderBounds::BoundValues::eHalfLengthZ], - zPosition + boundValues[CylinderBounds::BoundValues::eHalfLengthZ]}; - } else if (surfaceType == Surface::SurfaceType::Disc && - cast == BinningValue::binR) { - clipRange = {boundValues[RadialBounds::BoundValues::eMinR], - boundValues[RadialBounds::BoundValues::eMaxR]}; - } else { - throw std::runtime_error( - "PortalDetrayConverter: surface type not (yet) supported for " - "detray " - "conversion, only cylinder and disc are currently supported."); - } - - // Need to clip the parameter space to the surface dimension - std::vector clippedIndices = {}; - std::vector clippedBoundaries = {}; - clippedBoundaries.push_back(clipRange[0u]); - for (const auto [ib, b] : enumerate(boundaries)) { - if (ib > 0) { - unsigned int vI = vIndices[ib - 1u]; - ActsScalar highEdge = boundaries[ib]; - if (boundaries[ib - 1] >= clipRange[1u]) { - break; - } - if (highEdge <= clipRange[0u] || - std::abs(highEdge - clipRange[0u]) < 1e-5) { - continue; - } - if (highEdge > clipRange[1u]) { - highEdge = clipRange[1u]; - } - clippedIndices.push_back(vI); - clippedBoundaries.push_back(highEdge); - } - } - // Interpret the clipped information - // - // Clipped cylinder case - if (surfaceType == Surface::SurfaceType::Cylinder) { - for (auto [ib, b] : enumerate(clippedBoundaries)) { - if (ib > 0) { - // Create sub surfaces - std::array - subBoundValues = {}; - for (auto [ibv, bv] : enumerate(boundValues)) { - subBoundValues[ibv] = bv; - } - subBoundValues[CylinderBounds::BoundValues::eHalfLengthZ] = - 0.5 * (b - clippedBoundaries[ib - 1u]); - auto subBounds = std::make_shared(subBoundValues); - auto subTransform = Transform3::Identity(); - subTransform.pretranslate(Vector3( - 0., 0., - clippedBoundaries[ib - 1u] + - subBoundValues[CylinderBounds::BoundValues::eHalfLengthZ])); - auto subSurface = - Surface::makeShared(subTransform, subBounds); - - auto portalPayload = convertSurface(gctx, *subSurface, true); - portalPayload.mask.volume_link.link = clippedIndices[ib - 1u]; - portals.push_back(portalPayload); - } - } - } else { - for (auto [ib, b] : enumerate(clippedBoundaries)) { - if (ib > 0) { - // Create sub surfaces - std::array - subBoundValues = {}; - for (auto [ibv, bv] : enumerate(boundValues)) { - subBoundValues[ibv] = bv; - } - subBoundValues[RadialBounds::BoundValues::eMinR] = - clippedBoundaries[ib - 1u]; - subBoundValues[RadialBounds::BoundValues::eMaxR] = b; - auto subBounds = std::make_shared(subBoundValues); - auto subSurface = Surface::makeShared( - portal.surface().transform(gctx), subBounds); - - auto portalPayload = convertSurface(gctx, *subSurface, true); - portalPayload.mask.volume_link.link = clippedIndices[ib - 1u]; - portals.push_back(portalPayload); - } - } - } - - } else { - // End of world portal - // Write surface with invalid link - auto portalPayload = convertSurface(gctx, *surfaceAdjusted, true); - using NavigationLink = - typename DetrayDetector::surface_type::navigation_link; - portalPayload.mask.volume_link.link = - std::numeric_limits::max(); - - portals.push_back(portalPayload); - } - } - - return portals; -} - -detray::io::volume_payload Acts::DetrayConverter::convertVolume( - const Acts::GeometryContext& gctx, - const Acts::Experimental::DetectorVolume& volume, - const std::vector& detectorVolumes) { - detray::io::volume_payload volumePayload; - volumePayload.name = volume.name(); - volumePayload.index.link = findVolume(&volume, detectorVolumes); - volumePayload.transform = convertTransform(volume.transform(gctx)); - - // iterate over surfaces and portals keeping the same surf_pd.index_in_coll - std::size_t sIndex = 0; - for (const auto surface : volume.surfaces()) { - io::surface_payload surfacePayload = convertSurface(gctx, *surface, false); - - surfacePayload.index_in_coll = sIndex++; - surfacePayload.mask.volume_link.link = - volumePayload.index.link; // link surface' mask to volume - volumePayload.surfaces.push_back(surfacePayload); - } - - auto orientedSurfaces = - volume.volumeBounds().orientedSurfaces(volume.transform(gctx)); - - int portalCounter = 0; - for (const auto& [ip, p] : enumerate(volume.portals())) { - auto portals = - convertPortal(gctx, *p, ip, volume, orientedSurfaces, detectorVolumes); - std::for_each(portals.begin(), portals.end(), [&](auto& portalPayload) { - portalPayload.index_in_coll = sIndex++; - volumePayload.surfaces.push_back(portalPayload); - portalCounter++; - }); - } - - return volumePayload; -} diff --git a/Plugins/Detray/src/DetrayGeometryConverter.cpp b/Plugins/Detray/src/DetrayGeometryConverter.cpp new file mode 100644 index 00000000000..9e531e0be21 --- /dev/null +++ b/Plugins/Detray/src/DetrayGeometryConverter.cpp @@ -0,0 +1,350 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Plugins/Detray/DetrayGeometryConverter.hpp" + +#include "Acts/Detector/Detector.hpp" +#include "Acts/Detector/DetectorVolume.hpp" +#include "Acts/Detector/Portal.hpp" +#include "Acts/Navigation/PortalNavigation.hpp" +#include "Acts/Plugins/Json/DetrayJsonHelper.hpp" +#include "Acts/Surfaces/CylinderBounds.hpp" +#include "Acts/Surfaces/CylinderSurface.hpp" +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" +#include "Acts/Surfaces/RegularSurface.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Surfaces/SurfaceBounds.hpp" + +#include "detray/io/frontend/detector_writer.hpp" + +using namespace detray; + +namespace { + +/// Find the position of the volume to point to +/// +/// @param volume the volume to find +/// @param the collection of volumes +/// +/// @note return -1 if not found, to be interpreted by the caller +int findVolume( + const Acts::Experimental::DetectorVolume* volume, + const std::vector& volumes) { + auto candidate = std::find(volumes.begin(), volumes.end(), volume); + if (candidate != volumes.end()) { + return std::distance(volumes.begin(), candidate); + } + return -1; +} +} // namespace + +detray::io::transform_payload Acts::DetrayGeometryConverter::convertTransform( + const Transform3& t) { + detray::io::transform_payload tfPayload; + auto translation = t.translation(); + tfPayload.tr = {translation.x(), translation.y(), translation.z()}; + + const auto rotation = t.rotation(); + tfPayload.rot = {rotation(0, 0), rotation(0, 1), rotation(0, 2), + rotation(1, 0), rotation(1, 1), rotation(1, 2), + rotation(2, 0), rotation(2, 1), rotation(2, 2)}; + return tfPayload; +} + +detray::io::mask_payload Acts::DetrayGeometryConverter::convertMask( + const Acts::SurfaceBounds& bounds, bool portal) { + detray::io::mask_payload maskPayload; + auto [shape, boundaries] = DetrayJsonHelper::maskFromBounds(bounds, portal); + maskPayload.shape = static_cast(shape); + maskPayload.boundaries = static_cast>(boundaries); + // default maskPayload.volume_link + + return maskPayload; +} + +detray::io::surface_payload Acts::DetrayGeometryConverter::convertSurface( + const GeometryContext& gctx, const Surface& surface, bool portal) { + detray::io::surface_payload surfacePayload; + + surfacePayload.transform = convertTransform(surface.transform(gctx)); + surfacePayload.source = surface.geometryId().value(); + surfacePayload.barcode = std::nullopt; + surfacePayload.type = static_cast( + portal ? surface_id::e_portal + : (surface.geometryId().sensitive() > 0 + ? detray::surface_id::e_sensitive + : detray::surface_id::e_passive)); + surfacePayload.mask = convertMask(surface.bounds()); + return surfacePayload; +} + +std::vector +Acts::DetrayGeometryConverter::convertPortal( + const GeometryContext& gctx, const Experimental::Portal& portal, + std::size_t ip, const Experimental::DetectorVolume& volume, + const std::vector& orientedSurfaces, + const std::vector& detectorVolumes) { + std::vector portals{}; + + const RegularSurface& surface = portal.surface(); + const auto& volumeLinks = portal.portalNavigation(); + + // First assumption for outside link (along direction) + std::size_t outside = 1u; + + // Find out if you need to take the outside or inside volume + // for planar surfaces that's easy + if (surface.type() != Acts::Surface::SurfaceType::Cylinder) { + // Get the two volume center + const auto volumeCenter = volume.transform(gctx).translation(); + const auto surfaceCenter = surface.center(gctx); + const auto surfaceNormal = surface.normal(gctx, surfaceCenter); + // Get the direction from the volume to the surface, correct link + const auto volumeToSurface = surfaceCenter - volumeCenter; + if (volumeToSurface.dot(surfaceNormal) < 0.) { + outside = 0u; + } + } else { + // This is a cylinder portal, inner cover reverses the normal + if (ip == 3u) { + outside = 0u; + } + } + + const auto& outsideLink = volumeLinks[outside]; + // Grab the corresponding volume link + // If it is a single link, we are done + const auto* instance = outsideLink.instance(); + // Single link cast + auto singleLink = + dynamic_cast( + instance); + + auto [surfaceAdjusted, insidePointer] = orientedSurfaces[ip]; + + // Single link detected - just write it out, we use the oriented surface + // in order to make sure the size is adjusted + if (singleLink != nullptr) { + // Single link can be written out + std::size_t vLink = findVolume(singleLink->object(), detectorVolumes); + auto portalPayload = convertSurface(gctx, *surfaceAdjusted, true); + portalPayload.mask.volume_link.link = vLink; + portals.push_back(portalPayload); + } else { + // Multi link detected - 1D + auto multiLink1D = + dynamic_cast( + instance); + if (multiLink1D != nullptr) { + // Resolve the multi link 1D + auto boundaries = + multiLink1D->indexedUpdater.grid.axes()[0u]->getBinEdges(); + const auto& cast = multiLink1D->indexedUpdater.casts[0u]; + const auto& transform = multiLink1D->indexedUpdater.transform; + const auto& volumes = multiLink1D->indexedUpdater.extractor.dVolumes; + + // Apply the correction from local to global boundaries + ActsScalar gCorr = VectorHelpers::cast(transform.translation(), cast); + std::for_each(boundaries.begin(), boundaries.end(), + [&gCorr](ActsScalar& b) { b -= gCorr; }); + + // Get the volume indices + auto surfaceType = surfaceAdjusted->type(); + std::vector vIndices = {}; + for (const auto& v : volumes) { + vIndices.push_back(findVolume(v, detectorVolumes)); + } + + // Pick the surface dimension - via poly + std::array clipRange = {0., 0.}; + std::vector boundValues = surfaceAdjusted->bounds().values(); + if (surfaceType == Surface::SurfaceType::Cylinder && + cast == BinningValue::binZ) { + ActsScalar zPosition = surfaceAdjusted->center(gctx).z(); + clipRange = { + zPosition - boundValues[CylinderBounds::BoundValues::eHalfLengthZ], + zPosition + boundValues[CylinderBounds::BoundValues::eHalfLengthZ]}; + } else if (surfaceType == Surface::SurfaceType::Disc && + cast == BinningValue::binR) { + clipRange = {boundValues[RadialBounds::BoundValues::eMinR], + boundValues[RadialBounds::BoundValues::eMaxR]}; + } else { + throw std::runtime_error( + "PortalDetrayGeometryConverter: surface type not (yet) supported " + "for " + "detray " + "conversion, only cylinder and disc are currently supported."); + } + + // Need to clip the parameter space to the surface dimension + std::vector clippedIndices = {}; + std::vector clippedBoundaries = {}; + clippedBoundaries.push_back(clipRange[0u]); + for (const auto [ib, b] : enumerate(boundaries)) { + if (ib > 0) { + unsigned int vI = vIndices[ib - 1u]; + ActsScalar highEdge = boundaries[ib]; + if (boundaries[ib - 1] >= clipRange[1u]) { + break; + } + if (highEdge <= clipRange[0u] || + std::abs(highEdge - clipRange[0u]) < 1e-5) { + continue; + } + if (highEdge > clipRange[1u]) { + highEdge = clipRange[1u]; + } + clippedIndices.push_back(vI); + clippedBoundaries.push_back(highEdge); + } + } + // Interpret the clipped information + // + // Clipped cylinder case + if (surfaceType == Surface::SurfaceType::Cylinder) { + for (auto [ib, b] : enumerate(clippedBoundaries)) { + if (ib > 0) { + // Create sub surfaces + std::array + subBoundValues = {}; + for (auto [ibv, bv] : enumerate(boundValues)) { + subBoundValues[ibv] = bv; + } + subBoundValues[CylinderBounds::BoundValues::eHalfLengthZ] = + 0.5 * (b - clippedBoundaries[ib - 1u]); + auto subBounds = std::make_shared(subBoundValues); + auto subTransform = Transform3::Identity(); + subTransform.pretranslate(Vector3( + 0., 0., + clippedBoundaries[ib - 1u] + + subBoundValues[CylinderBounds::BoundValues::eHalfLengthZ])); + auto subSurface = + Surface::makeShared(subTransform, subBounds); + + auto portalPayload = convertSurface(gctx, *subSurface, true); + portalPayload.mask.volume_link.link = clippedIndices[ib - 1u]; + portals.push_back(portalPayload); + } + } + } else { + for (auto [ib, b] : enumerate(clippedBoundaries)) { + if (ib > 0) { + // Create sub surfaces + std::array + subBoundValues = {}; + for (auto [ibv, bv] : enumerate(boundValues)) { + subBoundValues[ibv] = bv; + } + subBoundValues[RadialBounds::BoundValues::eMinR] = + clippedBoundaries[ib - 1u]; + subBoundValues[RadialBounds::BoundValues::eMaxR] = b; + auto subBounds = std::make_shared(subBoundValues); + auto subSurface = Surface::makeShared( + portal.surface().transform(gctx), subBounds); + + auto portalPayload = convertSurface(gctx, *subSurface, true); + portalPayload.mask.volume_link.link = clippedIndices[ib - 1u]; + portals.push_back(portalPayload); + } + } + } + + } else { + // End of world portal + // Write surface with invalid link + auto portalPayload = convertSurface(gctx, *surfaceAdjusted, true); + using NavigationLink = + typename DetrayDetector::surface_type::navigation_link; + portalPayload.mask.volume_link.link = + std::numeric_limits::max(); + + portals.push_back(portalPayload); + } + } + + return portals; +} + +detray::io::volume_payload Acts::DetrayGeometryConverter::convertVolume( + DetrayConversionUtils::GeometryIdCache& geoIdCache, + const GeometryContext& gctx, + const Acts::Experimental::DetectorVolume& volume, + const std::vector& detectorVolumes, + const Acts::Logger& logger) { + ACTS_DEBUG("DetrayGeometryConverter: converting volume " + << volume.name() << " with " << volume.surfaces().size() + << " surfaces and " << volume.portals().size() << " portals"); + + detray::io::volume_payload volumePayload; + volumePayload.name = volume.name(); + volumePayload.index.link = findVolume(&volume, detectorVolumes); + volumePayload.transform = convertTransform(volume.transform(gctx)); + + // Remember the link + geoIdCache.volumeLinks[volume.geometryId()] = volumePayload.index.link; + + // iterate over surfaces and portals keeping the same surf_pd.index_in_coll + std::size_t sIndex = 0; + for (const auto surface : volume.surfaces()) { + io::surface_payload surfacePayload = convertSurface(gctx, *surface, false); + // Set the index in the collection & remember it in the cache + surfacePayload.index_in_coll = sIndex++; + geoIdCache.localSurfaceLinks.insert( + {surface->geometryId(), surfacePayload.index_in_coll.value()}); + // Set mask to volume link + surfacePayload.mask.volume_link.link = + volumePayload.index.link; // link surface' mask to volume + volumePayload.surfaces.push_back(surfacePayload); + } + + auto orientedSurfaces = + volume.volumeBounds().orientedSurfaces(volume.transform(gctx)); + + // Iterate over portals + int portalCounter = 0; + for (const auto& [ip, p] : enumerate(volume.portals())) { + auto portals = + convertPortal(gctx, *p, ip, volume, orientedSurfaces, detectorVolumes); + + GeometryIdentifier geoID = p->surface().geometryId(); + std::for_each(portals.begin(), portals.end(), [&](auto& portalPayload) { + // Set the index in the collection & remember it in the cache + portalPayload.index_in_coll = sIndex++; + geoIdCache.localSurfaceLinks.insert( + {geoID, portalPayload.index_in_coll.value()}); + // Add it to the surfaces + volumePayload.surfaces.push_back(portalPayload); + portalCounter++; + }); + } + ACTS_VERBOSE(" > " << volume.portals().size() + << " initial ACTS portals split into final " + << portalCounter << " detray portals"); + + return volumePayload; +} + +detray::io::detector_payload Acts::DetrayGeometryConverter::convertDetector( + DetrayConversionUtils::GeometryIdCache& geoIdCache, + const GeometryContext& gctx, const Acts::Experimental::Detector& detector, + const Acts::Logger& logger) { + ACTS_DEBUG("DetrayGeometryConverter: converting detector" + << detector.name() << " with " << detector.volumes().size() + << " volumes."); + + // The detector payload which will be handed back + detray::io::detector_payload detectorPayload; + + for (const auto volume : detector.volumes()) { + detectorPayload.volumes.push_back( + convertVolume(geoIdCache, gctx, *volume, detector.volumes(), logger)); + } + + return detectorPayload; +} diff --git a/Plugins/Detray/src/DetrayMaterialConverter.cpp b/Plugins/Detray/src/DetrayMaterialConverter.cpp new file mode 100644 index 00000000000..c777b331270 --- /dev/null +++ b/Plugins/Detray/src/DetrayMaterialConverter.cpp @@ -0,0 +1,228 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Plugins/Detray/DetrayMaterialConverter.hpp" + +#include "Acts/Detector/Detector.hpp" +#include "Acts/Material/BinnedSurfaceMaterial.hpp" +#include "Acts/Material/HomogeneousSurfaceMaterial.hpp" +#include "Acts/Plugins/Detray/DetrayConversionUtils.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Utilities/BinUtility.hpp" +#include "Acts/Utilities/BinningType.hpp" + +namespace { + +struct MaterialSurfaceSelector { + std::vector surfaces = {}; + + /// @param surface is the test surface + void operator()(const Acts::Surface* surface) { + if (surface->surfaceMaterial() != nullptr) { + if (std::find(surfaces.begin(), surfaces.end(), surface) == + surfaces.end()) { + surfaces.push_back(surface); + } + } + } +}; + +/// This creates dummy axes to allow homogeneous material for the moment +/// to be represented as grid surface material +std::vector homogeneousAxesPayloads() { + Acts::BinningData bDataX(Acts::BinningValue::binX, -1, 1); + bDataX.option = Acts::BinningOption::closed; + Acts::BinningData bDataY(Acts::BinningValue::binY, -1, 1); + bDataY.option = Acts::BinningOption::closed; + auto axisPayloadX = Acts::DetrayConversionUtils::convertBinningData(bDataX); + auto axisPayloadY = Acts::DetrayConversionUtils::convertBinningData(bDataY); + + return {axisPayloadX, axisPayloadY}; +} + +} // namespace + +detray::io::material_slab_payload +Acts::DetrayMaterialConverter::convertMaterialSlab( + const MaterialSlab& materialSlab) { + detray::io::material_slab_payload slab; + // Fill the material parameters and the thickness + const auto& material = materialSlab.material(); + slab.thickness = materialSlab.thickness(); + slab.mat = detray::io::material_payload{ + {material.X0(), material.L0(), material.Ar(), material.Z(), + material.massDensity(), material.molarDensity(), 0.}}; + slab.type = detray::io::material_id::slab; + return slab; +} + +detray::io::grid_payload +Acts::DetrayMaterialConverter::convertSurfaceMaterial( + const ISurfaceMaterial& material, const Logger& logger) { + detray::io::grid_payload + materialGrid; + + // Check the material types + // (1) homogeneous -> 1 x 1 bin grid with closed axes + auto homogeneousMaterial = + dynamic_cast(&material); + if (homogeneousMaterial != nullptr) { + ACTS_VERBOSE( + "DetrayMaterialConverter: found homogeneous surface material, this " + "will be modelled as a 1x1 bin grid"); + // A single bin entry: convert it and fill it + detray::io::material_slab_payload slab = convertMaterialSlab( + homogeneousMaterial->materialSlab(Vector3{0., 0., 0.})); + detray::io::grid_bin_payload slabBin{ + {0, 0}, {slab}}; + // Filling axes and bins + materialGrid.axes = homogeneousAxesPayloads(); + materialGrid.bins = {slabBin}; + return materialGrid; + } + // (2) - binned material -> convert into grid structure + auto binnedMaterial = dynamic_cast(&material); + if (binnedMaterial != nullptr) { + ACTS_VERBOSE("DetrayMaterialConverter: found binned surface material"); + + // BinUtility modifications + bool swapped = false; + // Get the bin utility (make a copy as we may modify it) + // Detray expects 2-dimensional grid, currently supported are + // x-y, r-phi, phi-z + BinUtility bUtility = binnedMaterial->binUtility(); + // Turn the bin value into a 2D grid + if (bUtility.dimensions() == 1u) { + if (bUtility.binningData()[0u].binvalue == BinningValue::binR) { + // Turn to R-Phi + bUtility += BinUtility(1u, -M_PI, M_PI, closed, BinningValue::binR); + } else if (bUtility.binningData()[0u].binvalue == BinningValue::binZ) { + // Turn to Phi-Z - swap needed + BinUtility nbUtility(1u, -M_PI, M_PI, closed, BinningValue::binPhi); + nbUtility += bUtility; + bUtility = std::move(nbUtility); + swapped = true; + } else { + std::runtime_error("Unsupported binning for Detray"); + } + } else if (bUtility.dimensions() == 2u && + bUtility.binningData()[0u].binvalue == BinningValue::binZ && + bUtility.binningData()[1u].binvalue == BinningValue::binPhi) { + BinUtility nbUtility(bUtility.binningData()[1u]); + nbUtility += bUtility.binningData()[0u]; + bUtility = std::move(nbUtility); + swapped = true; + } + + BinningValue bVal0 = bUtility.binningData()[0u].binvalue; + BinningValue bVal1 = bUtility.binningData()[1u].binvalue; + + // Translate into grid index type + detray::io::material_id gridIndexType = detray::io::material_id::unknown; + if (bVal0 == BinningValue::binR && bVal1 == BinningValue::binPhi) { + gridIndexType = detray::io::material_id::ring2_map; + } else if (bVal0 == BinningValue::binPhi && bVal1 == BinningValue::binZ) { + gridIndexType = detray::io::material_id::concentric_cylinder2_map; + } else if (bVal0 == BinningValue::binX && bVal1 == BinningValue::binY) { + gridIndexType = detray::io::material_id::rectangle2_map; + } else { + std::runtime_error("Unsupported binning for Detray"); + } + + detray::io::typed_link_payload linkPayload{ + gridIndexType, 0u}; + materialGrid.grid_link = linkPayload; + + // Now convert the (modified) bin utility + for (const auto& bData : bUtility.binningData()) { + auto axisPayload = DetrayConversionUtils::convertBinningData(bData); + materialGrid.axes.push_back(axisPayload); + } + + // Convert the material slabs from the material matrix + auto materialMatrix = binnedMaterial->fullMaterial(); + for (std::size_t ib1 = 0; ib1 < materialMatrix.size(); ++ib1) { + for (std::size_t ib0 = 0; ib0 < materialMatrix[0u].size(); ++ib0) { + // Translate into a local bin + std::size_t lb0 = swapped ? ib1 : ib0; + std::size_t lb1 = swapped ? ib0 : ib1; + detray::io::material_slab_payload slab = + convertMaterialSlab(materialMatrix[ib1][ib0]); + detray::io::grid_bin_payload slabBin{ + {static_cast(lb0), static_cast(lb1)}, + {slab}}; + // Fill into the grid + materialGrid.bins.push_back(slabBin); + } + } + return materialGrid; + } + + throw std::invalid_argument( + "DetrayMaterialConverter: unknown surface material type detected."); +} + +detray::io::detector_grids_payload +Acts::DetrayMaterialConverter::convertSurfaceMaterialGrids( + const DetrayConversionUtils::GeometryIdCache& geoIdCache, + const Experimental::Detector& detector, const Logger& logger) { + // The material grid payload + detray::io::detector_grids_payload + materialGrids; + + using DetrayMaterialGrid = + detray::io::grid_payload; + + // Loop over the volumes in order to assign the right volume links + for (const auto& volume : detector.volumes()) { + // Per volume surface selector + MaterialSurfaceSelector selector; + volume->visitSurfaces(selector); + ACTS_DEBUG("DetrayMaterialConverter: found " + << selector.surfaces.size() + << " surfaces/portals with material in volume " + << volume->name()); + // Find the voluem index first + auto volumeIndex = geoIdCache.volumeLinks.find(volume->geometryId()); + if (volumeIndex != geoIdCache.volumeLinks.end()) { + std::vector volumeMaterialGrids = {}; + // Now convert the surfaces + for (const auto& surface : selector.surfaces) { + // Find the surface index + auto surfaceIndices = + geoIdCache.localSurfaceLinks.equal_range(surface->geometryId()); + DetrayMaterialGrid materialGrid = + convertSurfaceMaterial(*surface->surfaceMaterial(), logger); + // Loop over the equal range and fill one grid each, this is needed + // as the initial portal could be split into multiple surfaces + for (auto itr = surfaceIndices.first; itr != surfaceIndices.second; + ++itr) { + // Fill the surface index + materialGrid.owner_link = + detray::io::single_link_payload{itr->second}; + // Fill the grid + volumeMaterialGrids.push_back(materialGrid); + } + } + // Register the grids of this volume + materialGrids.grids.insert({volumeIndex->second, volumeMaterialGrids}); + + } else { + ACTS_WARNING( + "DetrayMaterialConverter: volume not found in cache, should not " + "happen."); + } + } + // Return the material grids payload + return materialGrids; +} diff --git a/Tests/DownstreamProject/CMakeLists.txt b/Tests/DownstreamProject/CMakeLists.txt index 27f2f10583b..d41277017a5 100644 --- a/Tests/DownstreamProject/CMakeLists.txt +++ b/Tests/DownstreamProject/CMakeLists.txt @@ -7,7 +7,7 @@ find_package( Acts CONFIG REQUIRED - COMPONENTS Core Fatras PluginJson PluginLegacy PluginTGeo + COMPONENTS Core Fatras PluginJson PluginLegacy PluginTGeo PluginCovfie ) # place artifacts in GNU-like paths, e.g. binaries in `/bin` @@ -24,7 +24,13 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY add_executable(ShowActsVersion ShowActsVersion.cpp) target_link_libraries( ShowActsVersion - PRIVATE ActsCore ActsFatras ActsPluginJson ActsPluginLegacy ActsPluginTGeo + PRIVATE + ActsCore + ActsFatras + ActsPluginJson + ActsPluginLegacy + ActsPluginTGeo + ActsPluginCovfie ) option(DD4HEP "Build with DD4hep" ON) diff --git a/Tests/UnitTests/Core/Detector/CMakeLists.txt b/Tests/UnitTests/Core/Detector/CMakeLists.txt index 5ad3c568db0..eff298adb84 100644 --- a/Tests/UnitTests/Core/Detector/CMakeLists.txt +++ b/Tests/UnitTests/Core/Detector/CMakeLists.txt @@ -21,7 +21,7 @@ add_unittest(ReferenceGenerators ReferenceGeneratorsTests.cpp) add_unittest(SupportSurfacesHelper SupportSurfacesHelperTests.cpp) add_unittest(ProtoDetector ProtoDetectorTests.cpp) add_unittest(ProtoBinning ProtoBinningTests.cpp) -add_unittest(Portal PortalTests.cpp) +add_unittest(DetectorPortal PortalTests.cpp) add_unittest(PortalGenerators PortalGeneratorsTests.cpp) add_unittest(VolumeStructureBuilder VolumeStructureBuilderTests.cpp) add_unittest(MultiWireStructureBuilder MultiWireStructureBuilderTests.cpp) diff --git a/Tests/UnitTests/Core/Geometry/CMakeLists.txt b/Tests/UnitTests/Core/Geometry/CMakeLists.txt index 3b058a882b4..1fa422bb0a6 100644 --- a/Tests/UnitTests/Core/Geometry/CMakeLists.txt +++ b/Tests/UnitTests/Core/Geometry/CMakeLists.txt @@ -33,3 +33,4 @@ add_unittest(VolumeBounds VolumeBoundsTests.cpp) add_unittest(Volume VolumeTests.cpp) add_unittest(CylinderVolumeStack CylinderVolumeStackTests.cpp) add_unittest(PortalLink PortalLinkTests.cpp) +add_unittest(Portal PortalTests.cpp) diff --git a/Tests/UnitTests/Core/Geometry/PortalTests.cpp b/Tests/UnitTests/Core/Geometry/PortalTests.cpp new file mode 100644 index 00000000000..b36086dc8c9 --- /dev/null +++ b/Tests/UnitTests/Core/Geometry/PortalTests.cpp @@ -0,0 +1,602 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2024 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include +#include +#include +#include + +#include "Acts/Definitions/Units.hpp" +#include "Acts/Geometry/CylinderVolumeBounds.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/GridPortalLink.hpp" +#include "Acts/Geometry/Portal.hpp" +#include "Acts/Geometry/TrackingVolume.hpp" +#include "Acts/Geometry/TrivialPortalLink.hpp" +#include "Acts/Material/HomogeneousSurfaceMaterial.hpp" +#include "Acts/Surfaces/CylinderSurface.hpp" +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Surfaces/SurfaceMergingException.hpp" +#include "Acts/Utilities/BinningType.hpp" +#include "Acts/Utilities/ThrowAssert.hpp" + +#include + +using namespace Acts::UnitLiterals; + +namespace Acts::Test { + +auto logger = Acts::getDefaultLogger("UnitTests", Acts::Logging::VERBOSE); + +struct Fixture { + Logging::Level m_level; + Fixture() { + m_level = Acts::Logging::getFailureThreshold(); + Acts::Logging::setFailureThreshold(Acts::Logging::FATAL); + } + + ~Fixture() { Acts::Logging::setFailureThreshold(m_level); } +}; + +std::shared_ptr makeDummyVolume() { + return std::make_shared( + Transform3::Identity(), + std::make_shared(30_mm, 40_mm, 100_mm)); +} + +GeometryContext gctx; + +BOOST_FIXTURE_TEST_SUITE(Geometry, Fixture) + +BOOST_AUTO_TEST_SUITE(Portals) +BOOST_AUTO_TEST_SUITE(Merging) + +BOOST_AUTO_TEST_CASE(Cylinder) { + auto vol1 = makeDummyVolume(); + vol1->setVolumeName("vol1"); + auto vol2 = makeDummyVolume(); + vol2->setVolumeName("vol2"); + + auto cyl1 = Surface::makeShared( + Transform3{Translation3{Vector3::UnitZ() * -100_mm}}, 50_mm, 100_mm); + + auto cyl2 = Surface::makeShared( + Transform3{Translation3{Vector3::UnitZ() * 100_mm}}, 50_mm, 100_mm); + + Portal portal1{Direction::AlongNormal, + std::make_unique(cyl1, *vol1)}; + BOOST_CHECK(portal1.isValid()); + + BOOST_CHECK_EQUAL( + portal1 + .resolveVolume(gctx, Vector3{50_mm, 0_mm, -100_mm}, Vector3::UnitX()) + .value(), + vol1.get()); + + BOOST_CHECK_EQUAL( + portal1 + .resolveVolume(gctx, Vector3{50_mm, 0_mm, -100_mm}, -Vector3::UnitX()) + .value(), + nullptr); + + Portal portal2{Direction::AlongNormal, cyl2, *vol2}; + BOOST_CHECK(portal2.isValid()); + + BOOST_CHECK_EQUAL( + portal2 + .resolveVolume(gctx, Vector3{50_mm, 0_mm, 100_mm}, -Vector3::UnitX()) + .value(), + nullptr); + + BOOST_CHECK_EQUAL( + portal2 + .resolveVolume(gctx, Vector3{50_mm, 0_mm, 100_mm}, Vector3::UnitX()) + .value(), + vol2.get()); + + Portal portal3{gctx, std::make_unique(cyl2, *vol2), + nullptr}; + BOOST_CHECK(portal3.isValid()); + + BOOST_CHECK_NE(portal3.getLink(Direction::AlongNormal), nullptr); + BOOST_CHECK_EQUAL(portal3.getLink(Direction::OppositeNormal), nullptr); + + Portal portal4{gctx, nullptr, + std::make_unique(cyl2, *vol2)}; + BOOST_CHECK(portal4.isValid()); + + BOOST_CHECK_EQUAL(portal4.getLink(Direction::AlongNormal), nullptr); + BOOST_CHECK_NE(portal4.getLink(Direction::OppositeNormal), nullptr); + + // Not mergeable because 1 has portal along but 4 has portal oppsite + // ^ + // | + // portal1| portal2 + // +-------+-------+ + +---------------+ + // | | | | + // +---------------+ +-------+-------+ + // | + // | + // v + BOOST_CHECK_THROW( + Portal::merge(gctx, portal1, portal4, BinningValue::binZ, *logger), + PortalMergingException); + + // This call leaves both valid because the exception is thrown before the + // pointers are moved + BOOST_CHECK(portal1.isValid()); + BOOST_CHECK(portal2.isValid()); + + BOOST_CHECK_EQUAL( + portal2.resolveVolume(gctx, Vector3{50_mm, 0_mm, 50_mm}, Vector3::UnitX()) + .value(), + vol2.get()); + + BOOST_CHECK_EQUAL( + portal2 + .resolveVolume(gctx, Vector3{50_mm, 0_mm, 50_mm}, -Vector3::UnitX()) + .value(), + nullptr); + + // Cannot merge in binRPhi + BOOST_CHECK_THROW( + Portal::merge(gctx, portal1, portal2, BinningValue::binRPhi, *logger), + SurfaceMergingException); + + // The call above leaves both portals invalid because the exception is thrown + // after the pointers are moved (during durface merging) + BOOST_CHECK(!portal1.isValid()); + BOOST_CHECK(!portal2.isValid()); + + // ^ ^ + // | | + // portal1| portal2| + // +-------+-------+ + +-------+-------+ + // | | | | + // +---------------+ +---------------+ + + // Reset portals to valid to continue + portal1 = Portal{gctx, {.alongNormal = {cyl1, *vol1}}}; + portal2 = Portal{gctx, {.alongNormal = {cyl2, *vol2}}}; + + Portal merged12 = + Portal::merge(gctx, portal1, portal2, BinningValue::binZ, *logger); + BOOST_CHECK_NE(merged12.getLink(Direction::AlongNormal), nullptr); + BOOST_CHECK_EQUAL(merged12.getLink(Direction::OppositeNormal), nullptr); + + auto grid12 = dynamic_cast( + merged12.getLink(Direction::AlongNormal)); + BOOST_REQUIRE_NE(grid12, nullptr); + grid12->printContents(std::cout); + + BOOST_CHECK_EQUAL( + merged12 + .resolveVolume(gctx, Vector3{50_mm, 0_mm, -50_mm}, Vector3::UnitX()) + .value(), + vol1.get()); + + BOOST_CHECK_EQUAL( + merged12 + .resolveVolume(gctx, Vector3{50_mm, 0_mm, 50_mm}, Vector3::UnitX()) + .value(), + vol2.get()); + + BOOST_CHECK_EQUAL( + merged12 + .resolveVolume(gctx, Vector3{50_mm, 0_mm, -50_mm}, -Vector3::UnitX()) + .value(), + nullptr); + + BOOST_CHECK_EQUAL( + merged12 + .resolveVolume(gctx, Vector3{50_mm, 0_mm, 50_mm}, -Vector3::UnitX()) + .value(), + nullptr); + + portal1 = Portal{gctx, {.alongNormal = {cyl1, *vol1}}}; + + // Can't merge with self + BOOST_CHECK_THROW( + Portal::merge(gctx, portal1, portal1, BinningValue::binZ, *logger), + PortalMergingException); + + // Can't merge because the surfaces are the same + portal1 = Portal{gctx, {.alongNormal = {cyl1, *vol1}}}; + portal2 = Portal{gctx, {.alongNormal = {cyl1, *vol2}}}; + BOOST_CHECK_THROW( + Portal::merge(gctx, portal1, portal2, BinningValue::binZ, *logger), + AssertionFailureException); + + // Can't merge because surface has material + auto material = + std::make_shared(MaterialSlab{}); // vacuum + cyl2->assignSurfaceMaterial(material); + portal1 = Portal{gctx, {.alongNormal = {cyl1, *vol1}}}; + portal2 = Portal{gctx, {.alongNormal = {cyl2, *vol2}}}; + BOOST_CHECK_THROW( + Portal::merge(gctx, portal1, portal2, BinningValue::binZ, *logger), + PortalMergingException); +} + +BOOST_AUTO_TEST_CASE(Disc) { + auto vol1 = makeDummyVolume(); + vol1->setVolumeName("vol1"); + auto vol2 = makeDummyVolume(); + vol2->setVolumeName("vol2"); + auto vol3 = makeDummyVolume(); + vol3->setVolumeName("vol3"); + auto vol4 = makeDummyVolume(); + vol4->setVolumeName("vol4"); + + auto disc1 = Surface::makeShared( + Transform3::Identity(), std::make_shared(50_mm, 100_mm)); + + auto disc2 = Surface::makeShared( + Transform3::Identity(), std::make_shared(100_mm, 150_mm)); + + Portal portal1{ + gctx, {.alongNormal = {disc1, *vol1}, .oppositeNormal = {disc1, *vol2}}}; + + Portal portal2{ + gctx, {.alongNormal = {disc2, *vol3}, .oppositeNormal = {disc2, *vol4}}}; + + BOOST_CHECK(portal1.isValid()); + BOOST_CHECK(portal2.isValid()); + + BOOST_CHECK_EQUAL( + portal1.resolveVolume(gctx, Vector3{55_mm, 0_mm, 0_mm}, Vector3::UnitZ()) + .value(), + vol1.get()); + BOOST_CHECK_EQUAL( + portal1.resolveVolume(gctx, Vector3{55_mm, 0_mm, 0_mm}, -Vector3::UnitZ()) + .value(), + vol2.get()); + + BOOST_CHECK_EQUAL( + portal2.resolveVolume(gctx, Vector3{105_mm, 0_mm, 0_mm}, Vector3::UnitZ()) + .value(), + vol3.get()); + BOOST_CHECK_EQUAL( + portal2 + .resolveVolume(gctx, Vector3{105_mm, 0_mm, 0_mm}, -Vector3::UnitZ()) + .value(), + vol4.get()); + + BOOST_CHECK_THROW( + Portal::merge(gctx, portal1, portal2, BinningValue::binZ, *logger), + AssertionFailureException); + + BOOST_CHECK(portal1.isValid()); + BOOST_CHECK(portal2.isValid()); + + BOOST_CHECK_THROW( + Portal::merge(gctx, portal1, portal2, BinningValue::binPhi, *logger), + SurfaceMergingException); + + // Portals not valid anymore because they were moved before the exception was + // thrown + BOOST_CHECK(!portal1.isValid()); + BOOST_CHECK(!portal2.isValid()); + + // recreate them + portal1 = Portal{ + gctx, {.alongNormal = {disc1, *vol1}, .oppositeNormal = {disc1, *vol2}}}; + + portal2 = Portal{ + gctx, {.alongNormal = {disc2, *vol3}, .oppositeNormal = {disc2, *vol4}}}; + + // ^ ^ + // | | + // portal1| portal2| + // +-------+-------+ +-------+-------+ + // | | + | | + // +-------+-------+ +-------+-------+ + // | | + // | | + // v v + Portal merged12 = + Portal::merge(gctx, portal1, portal2, BinningValue::binR, *logger); + + BOOST_CHECK_EQUAL( + merged12.resolveVolume(gctx, Vector3{55_mm, 0_mm, 0_mm}, Vector3::UnitZ()) + .value(), + vol1.get()); + BOOST_CHECK_EQUAL( + merged12 + .resolveVolume(gctx, Vector3{55_mm, 0_mm, 0_mm}, -Vector3::UnitZ()) + .value(), + vol2.get()); + + BOOST_CHECK_EQUAL( + merged12 + .resolveVolume(gctx, Vector3{105_mm, 0_mm, 0_mm}, Vector3::UnitZ()) + .value(), + vol3.get()); + BOOST_CHECK_EQUAL( + merged12 + .resolveVolume(gctx, Vector3{105_mm, 0_mm, 0_mm}, -Vector3::UnitZ()) + .value(), + vol4.get()); + + // Can't merge because surface has material + auto material = + std::make_shared(MaterialSlab{}); // vacuum + disc2->assignSurfaceMaterial(material); + portal1 = Portal{ + gctx, {.alongNormal = {disc1, *vol1}, .oppositeNormal = {disc1, *vol2}}}; + portal2 = Portal{ + gctx, {.alongNormal = {disc2, *vol3}, .oppositeNormal = {disc2, *vol4}}}; + BOOST_CHECK_THROW( + Portal::merge(gctx, portal1, portal2, BinningValue::binR, *logger), + PortalMergingException); +} + +BOOST_AUTO_TEST_SUITE_END() // Merging + +BOOST_AUTO_TEST_SUITE(Fusing) + +BOOST_AUTO_TEST_CASE(Separated) { + auto vol1 = makeDummyVolume(); + vol1->setVolumeName("vol1"); + auto vol2 = makeDummyVolume(); + vol2->setVolumeName("vol2"); + + auto cyl1 = Surface::makeShared(Transform3::Identity(), + 50_mm, 100_mm); + + auto cyl2 = Surface::makeShared(Transform3::Identity(), + 60_mm, 100_mm); + + Portal portal1{gctx, {.oppositeNormal = {cyl1, *vol1}}}; + + Portal portal2{gctx, {.alongNormal = {cyl2, *vol2}}}; + + BOOST_CHECK(portal1.isValid()); + BOOST_CHECK(portal2.isValid()); + + // Surfaces have a 10mm gap in r + BOOST_CHECK_THROW(Portal::fuse(gctx, portal1, portal2, *logger), + PortalFusingException); + + BOOST_CHECK(portal1.isValid()); + BOOST_CHECK(portal2.isValid()); + + // Same way can't set cyl2 as other link + BOOST_CHECK_THROW(portal1.setLink(gctx, Direction::AlongNormal, cyl2, *vol2), + PortalFusingException); + BOOST_CHECK_EQUAL(portal1.getLink(Direction::AlongNormal), nullptr); + + Portal portal1b{gctx, {.oppositeNormal = {cyl1, *vol1}}}; + BOOST_CHECK(portal1b.isValid()); + + // portal1 portal1b + // +---+ +---+ + // | | | | + // | | | | + // <----+ | + <----+ | + // | | | | + // | | | | + // +---+ +---+ + BOOST_CHECK_THROW(Portal::fuse(gctx, portal1, portal1b, *logger), + PortalFusingException); + BOOST_CHECK(portal1.isValid()); + BOOST_CHECK(portal1b.isValid()); + + auto disc1 = Surface::makeShared( + Transform3::Identity(), std::make_shared(50_mm, 100_mm)); + + auto disc2 = Surface::makeShared( + Transform3{Translation3{Vector3{0, 0, 5_mm}}}, + std::make_shared(50_mm, 100_mm)); + + // portal2 portal2b + // +---+ +---+ + // | | | | + // | | | | + // | +----> + | +----> + // | | | | + // | | | | + // +---+ +---+ + Portal portal2b{gctx, {.alongNormal = {disc2, *vol2}}}; + + BOOST_CHECK_THROW(Portal::fuse(gctx, portal2, portal2b, *logger), + PortalFusingException); + BOOST_CHECK(portal1.isValid()); + BOOST_CHECK(portal2.isValid()); + + // portal2 portal2c + // +---+ +---+ + // | | | | + // | | | | + // <----+ | + <----+ +----> + // | | | | + // | | | | + // +---+ +---+ + Portal portal2c{ + gctx, {.alongNormal = {disc2, *vol1}, .oppositeNormal = {disc2, *vol2}}}; + BOOST_CHECK(portal2c.isValid()); + + BOOST_CHECK_THROW(Portal::fuse(gctx, portal2, portal2c, *logger), + PortalFusingException); + BOOST_CHECK(portal2.isValid()); + BOOST_CHECK(portal2c.isValid()); +} + +BOOST_AUTO_TEST_CASE(Success) { + auto vol1 = makeDummyVolume(); + vol1->setVolumeName("vol1"); + auto vol2 = makeDummyVolume(); + vol2->setVolumeName("vol2"); + + auto cyl1 = Surface::makeShared(Transform3::Identity(), + 50_mm, 100_mm); + + auto cyl2 = Surface::makeShared(Transform3::Identity(), + 50_mm, 100_mm); + + BOOST_CHECK(*cyl1 == *cyl2); + + // portal1 portal2 + // +---+ +---+ + // | | | | + // | | | | + // <----+ | + | +----> + // | | | | + // | | | | + // +---+ +---+ + Portal portal1{gctx, {.oppositeNormal = {cyl1, *vol1}}}; + BOOST_CHECK_EQUAL(&portal1.getLink(Direction::OppositeNormal)->surface(), + cyl1.get()); + + Portal portal2{gctx, {.alongNormal = {cyl2, *vol2}}}; + + BOOST_CHECK(portal1.isValid()); + BOOST_CHECK(portal2.isValid()); + + Portal portal3 = Portal::fuse(gctx, portal1, portal2, *logger); + // Input portals get invalidated by the fuse + BOOST_CHECK(!portal1.isValid()); + BOOST_CHECK(!portal2.isValid()); + BOOST_CHECK(portal3.isValid()); + + BOOST_CHECK_EQUAL(portal3.surface().surfaceMaterial(), nullptr); + + // Portal surface is set to the one from "along", because it gets set first + BOOST_CHECK_EQUAL(&portal3.surface(), cyl2.get()); + // "Opposite" gets the already-set surface set as well + BOOST_CHECK_EQUAL(&portal3.getLink(Direction::OppositeNormal)->surface(), + cyl2.get()); +} + +BOOST_AUTO_TEST_CASE(Material) { + auto vol1 = makeDummyVolume(); + auto vol2 = makeDummyVolume(); + + auto cyl1 = Surface::makeShared(Transform3::Identity(), + 50_mm, 100_mm); + + auto cyl2 = Surface::makeShared(Transform3::Identity(), + 50_mm, 100_mm); + + // portal1 portal2 + // +---+ +---+ + // | | | | + // | | | | + // <----+ | + | +----> + // | | | | + // | | | | + // +---+ +---+ + Portal portal1{gctx, {.oppositeNormal = {cyl1, *vol1}}}; + Portal portal2{gctx, {.alongNormal = {cyl2, *vol2}}}; + + auto material = + std::make_shared(MaterialSlab{}); // vacuum + + cyl1->assignSurfaceMaterial(material); + + Portal portal12 = Portal::fuse(gctx, portal1, portal2, *logger); + + // cyl1 had material, so this surface needs to be retained + BOOST_CHECK_EQUAL(&portal12.surface(), cyl1.get()); + BOOST_CHECK_EQUAL(portal12.surface().surfaceMaterial(), material.get()); + + // Reset portals + portal1 = Portal{gctx, {.oppositeNormal = {cyl1, *vol1}}}; + portal2 = Portal{gctx, {.alongNormal = {cyl2, *vol2}}}; + cyl2->assignSurfaceMaterial(material); + + // Both have material, this should fail + BOOST_CHECK_THROW(Portal::fuse(gctx, portal1, portal2, *logger), + PortalFusingException); + // Portals should stay valid + BOOST_CHECK(portal1.isValid()); + BOOST_CHECK(portal2.isValid()); + + cyl1->assignSurfaceMaterial(nullptr); + + portal12 = Portal::fuse(gctx, portal1, portal2, *logger); + + // cyl2 had material, so this surface needs to be retained + BOOST_CHECK_EQUAL(&portal12.surface(), cyl2.get()); + BOOST_CHECK_EQUAL(portal12.surface().surfaceMaterial(), material.get()); +} + +BOOST_AUTO_TEST_SUITE_END() // Fusing + +BOOST_AUTO_TEST_CASE(Construction) { + auto vol1 = makeDummyVolume(); + + // Displaced surfaces fail + auto disc1 = Surface::makeShared( + Transform3::Identity(), std::make_shared(50_mm, 100_mm)); + + auto disc2 = Surface::makeShared( + Transform3{Translation3{Vector3{0, 0, 5_mm}}}, + std::make_shared(50_mm, 100_mm)); + + BOOST_CHECK_THROW(std::make_unique( + gctx, std::make_unique(disc1, *vol1), + std::make_unique(disc2, *vol1)), + PortalFusingException); + + BOOST_CHECK_THROW((Portal{gctx, nullptr, nullptr}), std::invalid_argument); + BOOST_CHECK_THROW(Portal(gctx, {}), std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE(InvalidConstruction) { + BOOST_CHECK_THROW(Portal(Direction::AlongNormal, nullptr), + std::invalid_argument); + + auto vol1 = makeDummyVolume(); + + BOOST_CHECK_THROW(Portal(Direction::AlongNormal, nullptr, *vol1), + std::invalid_argument); + + auto disc1 = Surface::makeShared( + Transform3::Identity(), std::make_shared(50_mm, 100_mm)); + Portal portal(Direction::AlongNormal, disc1, *vol1); + + BOOST_CHECK_THROW(portal.setLink(gctx, Direction::AlongNormal, nullptr), + std::invalid_argument); +} + +BOOST_AUTO_TEST_CASE(PortalFill) { + auto vol1 = makeDummyVolume(); + auto vol2 = makeDummyVolume(); + + auto cyl1 = Surface::makeShared(Transform3::Identity(), + 50_mm, 100_mm); + + Portal portal1{gctx, {.oppositeNormal = {cyl1, *vol1}}}; + Portal portal2{gctx, {.alongNormal = {cyl1, *vol2}}}; + + // Fuse these to make portal 1 and 2 empty + Portal::fuse(gctx, portal1, portal2, *logger); + + BOOST_CHECK_THROW(portal1.fill(*vol2), std::logic_error); + + portal1 = Portal{gctx, {.oppositeNormal = {cyl1, *vol1}}}; + portal2 = Portal{gctx, {.alongNormal = {cyl1, *vol2}}}; + + BOOST_CHECK_EQUAL(portal1.getLink(Direction::AlongNormal), nullptr); + BOOST_CHECK_NE(portal1.getLink(Direction::OppositeNormal), nullptr); + + portal1.fill(*vol2); + BOOST_CHECK_NE(portal1.getLink(Direction::AlongNormal), nullptr); + BOOST_CHECK_NE(portal1.getLink(Direction::OppositeNormal), nullptr); + + BOOST_CHECK_THROW(portal1.fill(*vol2), std::logic_error); +} + +BOOST_AUTO_TEST_SUITE_END() // Portals + +BOOST_AUTO_TEST_SUITE_END() // Geometry + +} // namespace Acts::Test diff --git a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp index 8bf56b6855b..1e4afdf3b6c 100644 --- a/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp +++ b/Tests/UnitTests/Core/Propagator/NavigatorTests.cpp @@ -291,18 +291,18 @@ bool testNavigatorStateVectors(Navigator::State& state, std::size_t navSurf, /// @param [in] targetVol Target volume /// @param [in] targetLay Target layer /// @param [in] targetSurf Target surface -bool testNavigatorStatePointers( - Navigator::State& state, const TrackingVolume* worldVol, - const TrackingVolume* startVol, const Layer* startLay, - const Surface* startSurf, const Surface* currSurf, - const TrackingVolume* currVol, const TrackingVolume* targetVol, - const Layer* targetLay, const Surface* targetSurf) { - return ( - (state.worldVolume == worldVol) && (state.startVolume == startVol) && - (state.startLayer == startLay) && (state.startSurface == startSurf) && - (state.currentSurface == currSurf) && (state.currentVolume == currVol) && - (state.targetVolume == targetVol) && (state.targetLayer == targetLay) && - (state.targetSurface == targetSurf)); +bool testNavigatorStatePointers(Navigator::State& state, + const TrackingVolume* worldVol, + const TrackingVolume* startVol, + const Layer* startLay, const Surface* startSurf, + const Surface* currSurf, + const TrackingVolume* currVol, + const Surface* targetSurf) { + return ((state.worldVolume == worldVol) && (state.startVolume == startVol) && + (state.startLayer == startLay) && (state.startSurface == startSurf) && + (state.currentSurface == currSurf) && + (state.currentVolume == currVol) && + (state.targetSurface == targetSurf)); } // the surface cache & the creation of the geometry @@ -344,7 +344,7 @@ BOOST_AUTO_TEST_CASE(Navigator_status_methods) { BOOST_CHECK(testNavigatorStateVectors(state.navigation, 0u, 0u, 0u)); BOOST_CHECK(testNavigatorStatePointers(state.navigation, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, - nullptr, nullptr, nullptr)); + nullptr)); } ACTS_INFO(" b) Run with geometry but without resolving"); @@ -360,7 +360,7 @@ BOOST_AUTO_TEST_CASE(Navigator_status_methods) { BOOST_CHECK(testNavigatorStateVectors(state.navigation, 0u, 0u, 0u)); BOOST_CHECK(testNavigatorStatePointers(state.navigation, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, - nullptr, nullptr, nullptr)); + nullptr)); } ACTS_INFO( @@ -381,7 +381,7 @@ BOOST_AUTO_TEST_CASE(Navigator_status_methods) { BOOST_CHECK(testNavigatorStateVectors(state.navigation, 0u, 0u, 0u)); BOOST_CHECK(testNavigatorStatePointers(state.navigation, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, - nullptr, nullptr, nullptr)); + nullptr)); ACTS_INFO(" ii) Because of no target surface"); state.navigation.targetReached = false; @@ -390,7 +390,7 @@ BOOST_AUTO_TEST_CASE(Navigator_status_methods) { BOOST_CHECK(testNavigatorStateVectors(state.navigation, 0u, 0u, 0u)); BOOST_CHECK(testNavigatorStatePointers(state.navigation, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, - nullptr, nullptr, nullptr)); + nullptr)); ACTS_INFO(" iii) Because the target surface is reached"); const Surface* startSurf = tGeometry->getBeamline(); state.stepping.pos4.segment<3>(Acts::ePos0) = @@ -399,9 +399,9 @@ BOOST_AUTO_TEST_CASE(Navigator_status_methods) { state.navigation.targetSurface = targetSurf; navigator.postStep(state, stepper); BOOST_CHECK(testNavigatorStateVectors(state.navigation, 0u, 0u, 0u)); - BOOST_CHECK(testNavigatorStatePointers( - state.navigation, nullptr, nullptr, nullptr, nullptr, targetSurf, - nullptr, nullptr, nullptr, targetSurf)); + BOOST_CHECK(testNavigatorStatePointers(state.navigation, nullptr, nullptr, + nullptr, nullptr, targetSurf, + nullptr, targetSurf)); ACTS_INFO("(2) Test the initialisation"); ACTS_INFO(" a) Initialise without additional information"); @@ -416,16 +416,16 @@ BOOST_AUTO_TEST_CASE(Navigator_status_methods) { BOOST_CHECK(testNavigatorStateVectors(state.navigation, 0u, 0u, 0u)); BOOST_CHECK(testNavigatorStatePointers(state.navigation, worldVol, startVol, startLay, nullptr, nullptr, startVol, - nullptr, nullptr, nullptr)); + nullptr)); ACTS_INFO(" b) Initialise having a start surface"); state.navigation = Navigator::State(); state.navigation.startSurface = startSurf; navigator.initialize(state, stepper); BOOST_CHECK(testNavigatorStateVectors(state.navigation, 0u, 0u, 0u)); - BOOST_CHECK(testNavigatorStatePointers( - state.navigation, worldVol, startVol, startLay, startSurf, startSurf, - startVol, nullptr, nullptr, nullptr)); + BOOST_CHECK(testNavigatorStatePointers(state.navigation, worldVol, startVol, + startLay, startSurf, startSurf, + startVol, nullptr)); ACTS_INFO(" c) Initialise having a start volume"); state.navigation = Navigator::State(); @@ -434,7 +434,7 @@ BOOST_AUTO_TEST_CASE(Navigator_status_methods) { BOOST_CHECK(testNavigatorStateVectors(state.navigation, 0u, 0u, 0u)); BOOST_CHECK(testNavigatorStatePointers(state.navigation, worldVol, startVol, startLay, nullptr, nullptr, startVol, - nullptr, nullptr, nullptr)); + nullptr)); } } diff --git a/Tests/UnitTests/Plugins/CMakeLists.txt b/Tests/UnitTests/Plugins/CMakeLists.txt index a3ce9f373ef..0b48b9782f8 100644 --- a/Tests/UnitTests/Plugins/CMakeLists.txt +++ b/Tests/UnitTests/Plugins/CMakeLists.txt @@ -9,3 +9,4 @@ add_subdirectory_if(TGeo ACTS_BUILD_PLUGIN_TGEO) add_subdirectory_if(EDM4hep ACTS_BUILD_PLUGIN_EDM4HEP) add_subdirectory_if(FpeMonitoring ACTS_BUILD_PLUGIN_FPEMON) add_subdirectory_if(Podio ACTS_BUILD_PLUGIN_PODIO) +add_subdirectory_if(Covfie ACTS_BUILD_PLUGIN_TRACCC) diff --git a/Tests/UnitTests/Plugins/Covfie/CMakeLists.txt b/Tests/UnitTests/Plugins/Covfie/CMakeLists.txt new file mode 100644 index 00000000000..a95e6b65e25 --- /dev/null +++ b/Tests/UnitTests/Plugins/Covfie/CMakeLists.txt @@ -0,0 +1,2 @@ +set(unittest_extra_libraries ActsPluginCovfie) +add_unittest(CovfieFieldConversion CovfieFieldConversionTest.cpp) diff --git a/Tests/UnitTests/Plugins/Covfie/CovfieFieldConversionTest.cpp b/Tests/UnitTests/Plugins/Covfie/CovfieFieldConversionTest.cpp new file mode 100644 index 00000000000..3bf85cccf22 --- /dev/null +++ b/Tests/UnitTests/Plugins/Covfie/CovfieFieldConversionTest.cpp @@ -0,0 +1,217 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +// Acts include(s) +#include "Acts/Definitions/Units.hpp" +#include "Acts/MagneticField/ConstantBField.hpp" +#include "Acts/MagneticField/MagneticFieldContext.hpp" +#include "Acts/MagneticField/MagneticFieldProvider.hpp" +#include "Acts/MagneticField/SolenoidBField.hpp" + +// Covfie Plugin include(s) +#include "Acts/Plugins/Covfie/FieldConversion.hpp" + +// System include(s) +#include +#include +#include +#include + +// Boost include(s) +#include + +using namespace Acts::UnitLiterals; + +template +void checkMagneticFieldEqual(const Acts::MagneticFieldProvider& fieldProvider, + Acts::MagneticFieldProvider::Cache& cache, + view_t view, iterator_t points, + float error_margin_half_width) { + for (auto point : points) { + auto x = point[0], y = point[1], z = point[2]; + + auto lookupResult = fieldProvider.getField(Acts::Vector3{x, y, z}, cache); + if (!lookupResult.ok()) { + throw std::runtime_error{"Field lookup failure"}; + } + auto actsValueX = (*lookupResult)[0], actsValueY = (*lookupResult)[1], + actsValueZ = (*lookupResult)[2]; + + auto covfieValues = view.at(x, y, z); + auto covfieValueX = covfieValues[0], covfieValueY = covfieValues[1], + covfieValueZ = covfieValues[2]; + + auto isEqual = + std::abs(covfieValueX - actsValueX) <= error_margin_half_width && + std::abs(covfieValueY - actsValueY) <= error_margin_half_width && + std::abs(covfieValueZ - actsValueZ) <= error_margin_half_width; + + std::stringstream ss; + ss << "Fields are not equal at position (" << x << ", " << y << ", " << z + << "). Acts: (" << actsValueX << ", " << actsValueY << ", " << actsValueZ + << "), Covfie: (" << covfieValueX << ", " << covfieValueY << ", " + << covfieValueZ << ")" << std::endl; + + BOOST_CHECK_MESSAGE(isEqual, ss.str()); + } +} + +BOOST_AUTO_TEST_SUITE(CovfiePlugin) + +BOOST_AUTO_TEST_CASE(InterpolatedMagneticField1) { + auto localToGlobalBin_xyz = [](std::array binsXYZ, + std::array nBinsXYZ) { + return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) + + binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2)); + }; + + std::vector xPos = {0., 1., 2., 3.}; + std::vector yPos = {0., 1., 2., 3.}; + std::vector zPos = {0., 1., 2., 3.}; + + std::vector bField_xyz; + for (int i = 0; i < 64; i++) { + bField_xyz.push_back(Acts::Vector3(i, i, i)); + } + + Acts::MagneticFieldContext fieldContext; + auto actsField = Acts::fieldMapXYZ(localToGlobalBin_xyz, xPos, yPos, zPos, + bField_xyz, 1, 1, false); + Acts::MagneticFieldProvider::Cache cache = actsField.makeCache(fieldContext); + + Acts::CovfiePlugin::InterpolatedField field = + Acts::CovfiePlugin::covfieField(actsField); + typename Acts::CovfiePlugin::InterpolatedField::view_t view(field); + + std::array, 14> points = {{ + {0.f, 0.f, 0.f}, + {1.f, 1.f, 1.f}, + {2.f, 2.f, 2.f}, + {2.9f, 2.9f, 2.9f}, + {1.2f, 2.5f, 0.8f}, + {0.7f, 1.9f, 2.3f}, + {2.1f, 0.3f, 1.5f}, + {0.4f, 2.8f, 2.9f}, + {1.6f, 1.2f, 0.5f}, + {2.3f, 0.6f, 2.2f}, + {1.1f, 2.7f, 1.3f}, + {0.9f, 1.4f, 2.7f}, + {2.4f, 1.8f, 0.9f}, + {0.6f, 2.2f, 2.1f}, + }}; + + checkMagneticFieldEqual(actsField, cache, view, points, 0.0001); +} + +BOOST_AUTO_TEST_CASE(InterpolatedMagneticField2) { + auto localToGlobalBin_xyz = [](std::array binsXYZ, + std::array nBinsXYZ) { + return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) + + binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2)); + }; + + std::vector xPos = {8., 12., 16., 20.}; + std::vector yPos = {8., 12., 16., 20.}; + std::vector zPos = {8., 12., 16., 20.}; + + std::vector bField_xyz; + for (int i = 0; i < 64; i++) { + bField_xyz.push_back(Acts::Vector3(i, i * i * 0.01, i)); + } + + Acts::MagneticFieldContext fieldContext; + auto actsField = Acts::fieldMapXYZ(localToGlobalBin_xyz, xPos, yPos, zPos, + bField_xyz, 1, 1, false); + Acts::MagneticFieldProvider::Cache cache = actsField.makeCache(fieldContext); + + Acts::CovfiePlugin::InterpolatedField field = + Acts::CovfiePlugin::covfieField(actsField); + typename Acts::CovfiePlugin::InterpolatedField::view_t view(field); + + std::array, 14> points = {{ + {8.f, 8.f, 8.f}, + {12.f, 12.f, 12.f}, + {16.f, 16.f, 16.f}, + {19.9, 19.9, 19.9}, + {8.1f, 10.2f, 12.3f}, + {9.4f, 11.5f, 13.6f}, + {10.7f, 12.8f, 14.9f}, + {11.0f, 13.1f, 15.2f}, + {12.3f, 14.4f, 16.5f}, + {13.6f, 15.7f, 17.8f}, + {14.9f, 16.0f, 18.1f}, + {16.2f, 17.3f, 19.4f}, + {17.5f, 18.6f, 19.7f}, + {18.8f, 19.9f, 14.0f}, + }}; + + checkMagneticFieldEqual(actsField, cache, view, points, 0.0001f); +} + +BOOST_AUTO_TEST_CASE(ConstantMagneticField1) { + Acts::ConstantBField actsField(Acts::Vector3{1.3f, 2.5f, 2.f}); + Acts::MagneticFieldContext ctx; + Acts::MagneticFieldProvider::Cache cache = actsField.makeCache(ctx); + + Acts::CovfiePlugin::ConstantField field = + Acts::CovfiePlugin::covfieField(actsField); + typename Acts::CovfiePlugin::ConstantField::view_t view(field); + + std::array, 13> points = {{ + {8.f, 8.f, 8.f}, + {12.f, 12.f, 12.f}, + {16.f, 16.f, 16.f}, + {8.1f, 10.2f, 12.3f}, + {9.4f, 11.5f, 13.6f}, + {10.7f, 12.8f, 14.9f}, + {11.0f, 13.1f, 15.2f}, + {12.3f, 14.4f, 16.5f}, + {13.6f, 15.7f, 17.8f}, + {14.9f, 16.0f, 18.1f}, + {16.2f, 17.3f, 19.4f}, + {17.5f, 18.6f, 19.7f}, + {18.8f, 19.9f, 14.0f}, + }}; + + checkMagneticFieldEqual(actsField, cache, view, points, 0.0001f); +} + +BOOST_AUTO_TEST_CASE(SolenoidBField1) { + Acts::SolenoidBField::Config cfg{}; + cfg.length = 5.8_m; + cfg.radius = (2.56 + 2.46) * 0.5 * 0.5_m; + cfg.nCoils = 1154; + cfg.bMagCenter = 2_T; + Acts::SolenoidBField actsField(cfg); + Acts::MagneticFieldContext ctx; + Acts::MagneticFieldProvider::Cache cache = actsField.makeCache(ctx); + + Acts::CovfiePlugin::InterpolatedField field = Acts::CovfiePlugin::covfieField( + actsField, cache, {21UL, 21UL, 21UL}, {0., 0., 0.}, {20., 20., 20.}); + typename Acts::CovfiePlugin::InterpolatedField::view_t view(field); + + std::array, 13> points = {{ + {8.f, 8.f, 8.f}, + {12.f, 12.f, 12.f}, + {16.f, 16.f, 16.f}, + {8.1f, 10.2f, 12.3f}, + {9.4f, 11.5f, 13.6f}, + {10.7f, 12.8f, 14.9f}, + {11.0f, 13.1f, 15.2f}, + {12.3f, 14.4f, 16.5f}, + {13.6f, 15.7f, 17.8f}, + {14.9f, 16.0f, 18.1f}, + {16.2f, 17.3f, 19.4f}, + {17.5f, 18.6f, 19.7f}, + {18.8f, 19.9f, 14.0f}, + }}; + + checkMagneticFieldEqual(actsField, cache, view, points, 0.0001); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/cmake/ActsConfig.cmake.in b/cmake/ActsConfig.cmake.in index 127d8b14786..907e0469cc2 100644 --- a/cmake/ActsConfig.cmake.in +++ b/cmake/ActsConfig.cmake.in @@ -103,6 +103,10 @@ if(PluginDetray IN_LIST Acts_COMPONENTS) find_dependency(detray @detray_VERSION@ CONFIG EXACT) endif() +if (PluginCovfie IN_LIST Acts_COMPONENTS) + find_dependency(covfie @covfie_VERSION@ CONFIG EXACT) +endif() + # load **all** available components. we can not just include the requested # components since there can be interdependencies between them. if(NOT Acts_FIND_QUIETLY)