Merge branch 'main' into gx2f-oom

Core/include/Acts/MagneticField/MultiRangeBField.hpp

@@ -0,0 +1,67 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+#pragma once
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/MagneticField/MagneticFieldError.hpp"
+#include "Acts/MagneticField/MagneticFieldProvider.hpp"
+#include "Acts/Utilities/RangeXD.hpp"
+
+namespace Acts {
+
+/// @ingroup MagneticField
+///
+/// @brief Magnetic field provider modelling a magnetic field consisting of
+/// several (potentially overlapping) regions of constant values.
+class MultiRangeBField final : public MagneticFieldProvider {
+ private:
+  struct Cache {
+    explicit Cache(const MagneticFieldContext& /*unused*/);
+
+    std::optional<std::size_t> index = {};
+  };
+
+  using BFieldRange = std::pair<RangeXD<3, double>, Vector3>;
+
+  // The different ranges and their corresponding field vectors. Note that
+  // regions positioned _later_ in this vector take priority over earlier
+  // regions.
+  std::vector<BFieldRange> fieldRanges;
+
+ public:
+  /// @brief Construct a magnetic field from a vector of ranges.
+  ///
+  /// @warning These ranges are listed in increasing order of precedence,
+  /// i.e. ranges further along the vector have higher priority.
+  explicit MultiRangeBField(const std::vector<BFieldRange>& ranges);
+
+  explicit MultiRangeBField(std::vector<BFieldRange>&& ranges);
+
+  /// @brief Construct a cache object.
+  MagneticFieldProvider::Cache makeCache(
+      const MagneticFieldContext& mctx) const override;
+
+  /// @brief Request the value of the magnetic field at a given position.
+  ///
+  /// @param [in] position Global 3D position for the lookup.
+  /// @param [in, out] cache Cache object.
+  /// @returns A successful value containing a field vector if the given
+  /// location is contained inside any of the regions, or a failure value
+  /// otherwise.
+  Result<Vector3> getField(const Vector3& position,
+                           MagneticFieldProvider::Cache& cache) const override;
+
+  /// @brief Get the field gradient at a given position.
+  ///
+  /// @warning This is not currently implemented.
+  Result<Vector3> getFieldGradient(
+      const Vector3& position, ActsMatrix<3, 3>& /*unused*/,
+      MagneticFieldProvider::Cache& cache) const override;
+};
+}  // namespace Acts

Core/include/Acts/Utilities/AlgebraHelpers.hpp

@@ -178,7 +178,7 @@ inline ActsMatrix<A::RowsAtCompileTime, B::ColsAtCompileTime> blockedMult(
 /// Calculate the inverse of an Eigen matrix after checking if it can be
 /// numerically inverted. This allows to catch potential FPEs before they occur.
 /// For matrices up to 4x4, the inverse is computed directly. For larger
-/// matrices, the FullPivLU is used.
+/// matrices, and dynamic matrices the FullPivLU is used.
 ///
 /// @tparam Derived Eigen derived concrete type
 /// @tparam Result Eigen result type defaulted to input type

Core/src/MagneticField/CMakeLists.txt

@@ -1,4 +1,8 @@
 target_sources(
     ActsCore
-    PRIVATE BFieldMapUtils.cpp SolenoidBField.cpp MagneticFieldError.cpp
+    PRIVATE
+        BFieldMapUtils.cpp
+        SolenoidBField.cpp
+        MagneticFieldError.cpp
+        MultiRangeBField.cpp
 )

Core/src/MagneticField/MultiRangeBField.cpp

@@ -0,0 +1,78 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+#include "Acts/MagneticField/MultiRangeBField.hpp"
+
+namespace Acts {
+
+MultiRangeBField::Cache::Cache(const MagneticFieldContext& /*unused*/) {}
+
+MultiRangeBField::MultiRangeBField(const std::vector<BFieldRange>& ranges)
+    : fieldRanges(ranges) {}
+
+MultiRangeBField::MultiRangeBField(
+    std::vector<MultiRangeBField::BFieldRange>&& ranges)
+    : fieldRanges(std::move(ranges)) {}
+
+MagneticFieldProvider::Cache MultiRangeBField::makeCache(
+    const MagneticFieldContext& mctx) const {
+  return MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
+}
+
+Result<Vector3> MultiRangeBField::getField(
+    const Vector3& position, MagneticFieldProvider::Cache& cache) const {
+  // Because we assume that the regions are in increasing order of
+  // precedence, we can iterate over the array, remembering the _last_
+  // region that contained the given point. At the end of the loop, this
+  // region will then be the one we query for its field vector.
+  std::optional<std::size_t> foundRange = {};
+
+  // The cache object for this type contains an optional integer; if the
+  // integer is set, it indicates the index of the region in which the last
+  // access succeeded. This acts as a cache because if the current access is
+  // still within that region, it is guaranteed that none of the regions
+  // with lower priority -- i.e. that are earlier in the region vector --
+  // can be relevant to the current access. Thus, we request the cache index
+  // if it exists and perform a membership check on it; if that succeeds, we
+  // remember the corresponding region as a candidate.
+  if (Cache& lCache = cache.as<Cache>();
+      lCache.index.has_value() &&
+      std::get<0>(fieldRanges[*lCache.index])
+          .contains({position[0], position[1], position[2]})) {
+    foundRange = *lCache.index;
+  }
+
+  // Now, we iterate over the ranges. If we already have a range candidate,
+  // we start iteration just after the existing candidate; otherwise we
+  // start from the beginning.
+  for (std::size_t i = (foundRange.has_value() ? (*foundRange) + 1 : 0);
+       i < fieldRanges.size(); ++i) {
+    if (std::get<0>(fieldRanges[i])
+            .contains({position[0], position[1], position[2]})) {
+      foundRange = i;
+    }
+  }
+
+  // Update the cache using the result of this access.
+  cache.as<Cache>().index = foundRange;
+
+  // If we found a valid range, return the corresponding vector; otherwise
+  // return an out-of-bounds error.
+  if (foundRange.has_value()) {
+    return Result<Vector3>::success(std::get<1>(fieldRanges[*foundRange]));
+  } else {
+    return Result<Vector3>::failure(MagneticFieldError::OutOfBounds);
+  }
+}
+
+Result<Vector3> MultiRangeBField::getFieldGradient(
+    const Vector3& position, ActsMatrix<3, 3>& /*unused*/,
+    MagneticFieldProvider::Cache& cache) const {
+  return getField(position, cache);
+}
+}  // namespace Acts

Examples/Io/Obj/CMakeLists.txt

@@ -2,6 +2,7 @@ add_library(
     ActsExamplesIoObj
     SHARED
     src/ObjTrackingGeometryWriter.cpp
+    src/ObjSimHitWriter.cpp
     src/ObjPropagationStepsWriter.cpp
 )
 

Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp

@@ -0,0 +1,84 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "ActsExamples/EventData/SimHit.hpp"
+#include "ActsExamples/Framework/ProcessCode.hpp"
+#include "ActsExamples/Framework/WriterT.hpp"
+
+#include <cstdint>
+#include <mutex>
+#include <string>
+
+namespace ActsExamples {
+struct AlgorithmContext;
+
+/// Write out a simhit collection before detector digitization as wavefront obj
+/// file(s per event).
+///
+/// This writes two files per event, one for the hits one for the trajectory.
+/// The latter can be smoothed using a spline interpolation.
+///
+///     event000000001-<stem>.obj
+///     event000000001-<stem>_trajectory.obj
+///     event000000002-<stem>.obj
+///     event000000002-<stem>_trajectory.obj
+///
+class ObjSimHitWriter : public WriterT<SimHitContainer> {
+ public:
+  struct Config {
+    /// Input sim hit collection to write.
+    std::string inputSimHits;
+    /// Where to place output files
+    std::string outputDir;
+    /// Output filename stem.
+    std::string outputStem = "simhits";
+    /// Number of decimal digits for floating point precision in output.
+    std::size_t outputPrecision = std::numeric_limits<float>::max_digits10;
+    /// Draw line connections between hits
+    bool drawConnections = true;
+    /// Momentum threshold for hits
+    double momentumThreshold = 0.05 * Acts::UnitConstants::GeV;
+    /// Momentum threshold for trajectories
+    double momentumThresholdTraj = 0.05 * Acts::UnitConstants::GeV;
+    /// Number of points to interpolate between hits
+    std::size_t nInterpolatedPoints = 10;
+  };
+
+  /// Construct the particle writer.
+  ///
+  /// @param config is the configuration object
+  /// @param level is the logging level
+  ObjSimHitWriter(const Config& config, Acts::Logging::Level level);
+
+  /// Ensure underlying file is closed.
+  ~ObjSimHitWriter() override = default;
+
+  /// End-of-run hook
+  ProcessCode finalize() override;
+
+  /// Get readonly access to the config parameters
+  const Config& config() const { return m_cfg; }
+
+ protected:
+  /// Type-specific write implementation.
+  ///
+  /// @param[in] ctx is the algorithm context
+  /// @param[in] hits are the hits to be written
+  ProcessCode writeT(const AlgorithmContext& ctx,
+                     const SimHitContainer& hits) override;
+
+ private:
+  Config m_cfg;
+  std::mutex m_writeMutex;
+};
+
+}  // namespace ActsExamples

Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp

@@ -6,7 +6,7 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
-#include "ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp"
+#include "ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp"
 
 namespace ActsExamples {