From 6970e2cb9c51d227fe1e12562600bc06d8152c50 Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Wed, 28 Aug 2024 12:53:57 +0200 Subject: [PATCH 01/10] feat: detray material conversion (#3546) This PR has a brush over the the `Detray/Plugin` conversion and adds a first version of a material converter. It does: - changes `DetrayConverter` into a class that also allows screen output - splits off `DetrayGeometryConverter` for geometry conversion - introduce `DetrayMaterialConverter` for the surface material conversion - introduce `DetrayConversionUtils` for common code This PR should be followed by one by @EleniXoch that introduces a `DetraySurfaceGridConverter` and the first prototype for surfaced based geometries should be complete. --- Core/include/Acts/Utilities/BinningData.hpp | 2 +- Examples/Python/src/Detray.cpp | 25 +- Plugins/Detray/CMakeLists.txt | 9 +- .../Plugins/Detray/DetrayConversionUtils.hpp | 78 ++++ .../Acts/Plugins/Detray/DetrayConverter.hpp | 206 +++++------ .../Detray/DetrayGeometryConverter.hpp | 112 ++++++ .../Detray/DetrayMaterialConverter.hpp | 61 +++ Plugins/Detray/src/DetrayConversionUtils.cpp | 82 ++++ Plugins/Detray/src/DetrayConverter.cpp | 301 +-------------- .../Detray/src/DetrayGeometryConverter.cpp | 350 ++++++++++++++++++ .../Detray/src/DetrayMaterialConverter.cpp | 228 ++++++++++++ 11 files changed, 1008 insertions(+), 446 deletions(-) create mode 100644 Plugins/Detray/include/Acts/Plugins/Detray/DetrayConversionUtils.hpp create mode 100644 Plugins/Detray/include/Acts/Plugins/Detray/DetrayGeometryConverter.hpp create mode 100644 Plugins/Detray/include/Acts/Plugins/Detray/DetrayMaterialConverter.hpp create mode 100644 Plugins/Detray/src/DetrayConversionUtils.cpp create mode 100644 Plugins/Detray/src/DetrayGeometryConverter.cpp create mode 100644 Plugins/Detray/src/DetrayMaterialConverter.cpp 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/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/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; +} From ef092ca9fd31e795859a9cecfa49fa18b86d5d1e Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Wed, 28 Aug 2024 15:54:55 +0200 Subject: [PATCH 02/10] chore: Make license check work with single files (#3569) Co-authored-by: kodiakhq[bot] <49736102+kodiakhq[bot]@users.noreply.github.com> --- CI/check_license.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 0258d12b53940d390f019a80ec922a102771aaf9 Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Wed, 28 Aug 2024 20:52:24 +0200 Subject: [PATCH 03/10] refactor(gx2f): early exit for `addToGx2fSums` (#3568) --- .../TrackFitting/GlobalChiSquareFitter.hpp | 121 ++++++++++-------- 1 file changed, 65 insertions(+), 56 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index deed734ac92..e8d283943da 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(); + + const auto safeInvCovMeasurement = safeInverse(covarianceMeasurement); + if (!safeInvCovMeasurement) { + ACTS_WARNING("addToGx2fSums: safeInvCovMeasurement failed."); + ACTS_VERBOSE(" covarianceMeasurement:\n" << covarianceMeasurement); + return; + } - ActsVector measurement = trackState.template calibrated(); + const BoundVector predicted = trackState.predicted(); - ActsSquareMatrix covarianceMeasurement = - trackState.template calibratedCovariance(); + const ActsVector measurement = + trackState.template calibrated(); - ActsMatrix projector = + 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: From 468e9372a3d0e8910b7384c294cd95c9644a7b7b Mon Sep 17 00:00:00 2001 From: "Alexander J. Pfleger" <70842573+AJPfleger@users.noreply.github.com> Date: Thu, 29 Aug 2024 01:14:16 +0200 Subject: [PATCH 04/10] refactor(gx2f): early exit for the Actor (#3566) The design for the Actor of the GX2F is clear now. Since we have only a clear forward propagation, we can skip all navigation states without a surface. Therefore, we do not need to reserve any options afterwards and early exits are possible. --- .../TrackFitting/GlobalChiSquareFitter.hpp | 325 +++++++++--------- 1 file changed, 168 insertions(+), 157 deletions(-) diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp index e8d283943da..829cd022f60 100644 --- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp +++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp @@ -464,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; + } + + ++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, ..."); + } - // 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; + // 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."); + + return; + } - // 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 + // 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; } }; From 3b32e5d816413ec0c88ebbe8734d9f48b4d223d6 Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Thu, 29 Aug 2024 07:41:41 +0200 Subject: [PATCH 05/10] refactor: Rework G4 surface mapping to make it more robust (#3562) Should now work with ITk annulus bounds and be more robust in general. Also adds some dedicated debugging facilities to ease debugging in case of missing surfaces. --- .../Geant4/SensitiveSurfaceMapper.hpp | 70 ++-- .../Geant4/src/Geant4Simulation.cpp | 9 +- .../Geant4/src/SensitiveSurfaceMapper.cpp | 332 +++++++++++++++--- .../Python/python/acts/examples/simulation.py | 2 +- Examples/Python/src/Geant4Component.cpp | 86 +++-- Examples/Scripts/Python/geant4.py | 2 + 6 files changed, 391 insertions(+), 110 deletions(-) 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/Geant4Simulation.cpp b/Examples/Algorithms/Geant4/src/Geant4Simulation.cpp index fce1c03db1e..7564f923462 100644 --- a/Examples/Algorithms/Geant4/src/Geant4Simulation.cpp +++ b/Examples/Algorithms/Geant4/src/Geant4Simulation.cpp @@ -282,8 +282,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); 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/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/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/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 From 809b378e51249733567ae24bffe6cdffa9058503 Mon Sep 17 00:00:00 2001 From: Stephen Nicholas Swatman Date: Thu, 29 Aug 2024 09:37:24 +0200 Subject: [PATCH 06/10] feat: Add covfie magnetic field plugin (#3479) This commit adds a new Acts plugin that adds support for magnetic fields implemented using the covfie library. This commit is based on #3117. Closes #3117. Depends on #3478. Virtually all credit goes to @fredevb. Co-authored-by: Fred <92879080+fredevb@users.noreply.github.com> --- Examples/Python/CMakeLists.txt | 3 + Examples/Python/src/Covfie.cpp | 55 +++++ Examples/Python/src/CovfieStub.cpp | 13 ++ Examples/Python/src/ModuleEntry.cpp | 2 + Examples/Python/tests/helpers/__init__.py | 7 + Examples/Python/tests/test_covfie.py | 63 +++++ Plugins/CMakeLists.txt | 1 + Plugins/Covfie/CMakeLists.txt | 16 ++ .../Acts/Plugins/Covfie/FieldConversion.hpp | 64 ++++++ Plugins/Covfie/src/FieldConversion.cpp | 197 ++++++++++++++++ Tests/DownstreamProject/CMakeLists.txt | 10 +- Tests/UnitTests/Plugins/CMakeLists.txt | 1 + Tests/UnitTests/Plugins/Covfie/CMakeLists.txt | 2 + .../Covfie/CovfieFieldConversionTest.cpp | 217 ++++++++++++++++++ cmake/ActsConfig.cmake.in | 4 + 15 files changed, 653 insertions(+), 2 deletions(-) create mode 100644 Examples/Python/src/Covfie.cpp create mode 100644 Examples/Python/src/CovfieStub.cpp create mode 100644 Examples/Python/tests/test_covfie.py create mode 100644 Plugins/Covfie/CMakeLists.txt create mode 100644 Plugins/Covfie/include/Acts/Plugins/Covfie/FieldConversion.hpp create mode 100644 Plugins/Covfie/src/FieldConversion.cpp create mode 100644 Tests/UnitTests/Plugins/Covfie/CMakeLists.txt create mode 100644 Tests/UnitTests/Plugins/Covfie/CovfieFieldConversionTest.cpp 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/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/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/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/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/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) From d6b112d4ccda8e560f91e19e215f32fe4b8cee15 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Thu, 29 Aug 2024 11:10:13 +0200 Subject: [PATCH 07/10] refactor: Remove `calculateTrackQuantities` from Core CKF (#3567) Originally https://github.com/acts-project/acts/pull/3536 but reverted with https://github.com/acts-project/acts/pull/3565 due to FPEs in Athena. --- Remove `calculateTrackQuantities` which recalculates measurements and holes from the track since we already accumulate these numbers inside the CKF actor. --- .../include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) 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); } }; From 45692302448dc93af3cb53385256a9af5f92a69c Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Thu, 29 Aug 2024 13:09:18 +0200 Subject: [PATCH 08/10] feat: Gen3 geometry Portals (#3501) Part of #3502 This PR implements *Portals*. 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. Portals can be **fused** and **merged**. **Fusing** is the combination of two portal linkson the same logical surfaces. The actual surface instances can be different, as long as they are geometrically equivalent (within numerical precistion). The resulting portal will have one portal along the shared surface's normal vector, and one opposite that vector. ``` portal1 portal2 +---+ +---+ | | | | | | | | <----+ | + | +----> | | | | | | | | +---+ +---+ ``` The input portals need to have compatible link loadout, e.g. one portal needs to have the *along normal* slot filled, and the other one needs to have the *opposite normal* slot filled. If portals share a filled slot, the function throws an exception. **Merging** 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 ``` 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 --- Core/include/Acts/Definitions/Algebra.hpp | 2 + Core/include/Acts/Geometry/Portal.hpp | 234 +++++++ Core/src/Geometry/CMakeLists.txt | 1 + Core/src/Geometry/Portal.cpp | 332 ++++++++++ Tests/UnitTests/Core/Detector/CMakeLists.txt | 2 +- Tests/UnitTests/Core/Geometry/CMakeLists.txt | 1 + Tests/UnitTests/Core/Geometry/PortalTests.cpp | 602 ++++++++++++++++++ 7 files changed, 1173 insertions(+), 1 deletion(-) create mode 100644 Core/include/Acts/Geometry/Portal.hpp create mode 100644 Core/src/Geometry/Portal.cpp create mode 100644 Tests/UnitTests/Core/Geometry/PortalTests.cpp 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/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/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 From e2063f04b02fd942ee4ffd08b105c04d18e529e7 Mon Sep 17 00:00:00 2001 From: Benjamin Huth <37871400+benjaminhuth@users.noreply.github.com> Date: Thu, 29 Aug 2024 16:13:30 +0200 Subject: [PATCH 09/10] fix: Add G4 log level tweaking on algorithm initialization (#3570) This adjusts the G4 loglevel to the level provided by the algorithm at algorithm construction time. Removes some unused loglevel integers from the `Geant4Manager`. --- .../include/ActsExamples/Geant4/Geant4Manager.hpp | 8 +++----- Examples/Algorithms/Geant4/src/Geant4Manager.cpp | 14 ++++++-------- .../Algorithms/Geant4/src/Geant4Simulation.cpp | 11 ++++++----- 3 files changed, 15 insertions(+), 18 deletions(-) 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/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 7564f923462..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"); } @@ -339,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); From ba521200fce09c3a32e11b41f4f237fe27c6aaf2 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Thu, 29 Aug 2024 20:46:11 +0200 Subject: [PATCH 10/10] refactor: Remove target volume estimation from `Navigator` (#3242) This is https://github.com/acts-project/acts/pull/3217 after it was reverted with https://github.com/acts-project/acts/pull/3241. There were a couple of failures monitored in Athena which could not be fixed right away. This is an attempt to remove (IMO) unnecessary complexity from the Navigator. We have a target volume estimation which runs a couple of times to guess the target volume based on an intersection with the target surface. IMO if we know the target surface we should also know the target volume and the user is responsible for passing this information. --- Core/include/Acts/Propagator/Navigator.hpp | 155 +----------------- .../Core/Propagator/NavigatorTests.cpp | 48 +++--- 2 files changed, 26 insertions(+), 177 deletions(-) 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/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)); } }