feat: introduction obj simhit writer

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,94 @@
+// 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/.
+
+// 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/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>
+
+using namespace Acts::UnitLiterals;
+
+namespace ActsExamples {
+struct AlgorithmContext;
+
+/// Write out a simhit collection before detector digitization
+/// as wavefront obj files.
+///
+/// This writes one file per event containing information about the
+/// global space points, momenta (before and after interaction) and hit index
+/// into the configured output directory. By default it writes to the
+/// current working directory. Files are named using the following schema
+///
+///     event000000001-<stem>.obj
+///     event000000002-<stem>.obj
+///     ...class ObjSimHitWriter final : public WriterT<SimHitContainer> {
+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_GeV;
+    /// Momentum threshold for trajectories
+    double momentumThresholdTraj = 0.05_GeV;
+    /// Number of points to interpolate between hits
+    unsigned int 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 {
 

Examples/Io/Obj/src/ObjSimHitWriter.cpp

@@ -0,0 +1,183 @@
+// 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/.
+
+// 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 "ActsExamples/Io/Obj/ObjSimHitWriter.hpp"
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Common.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "ActsExamples/EventData/SimHit.hpp"
+#include "ActsExamples/Framework/AlgorithmContext.hpp"
+#include "ActsExamples/Utilities/Paths.hpp"
+#include "ActsFatras/EventData/Barcode.hpp"
+#include "ActsFatras/EventData/Hit.hpp"
+
+#include <fstream>
+#include <stdexcept>
+#include <unordered_map>
+#include <vector>
+
+#include <unsupported/Eigen/Splines>
+
+namespace {
+
+template <typename input_vector_type>
+std::vector<Acts::Vector3> interpolatedPoints(
+    const std::vector<input_vector_type>& inputs, unsigned int nPoints) {
+  Eigen::MatrixXd points(3, inputs.size());
+  for (unsigned int i = 0; i < inputs.size(); ++i) {
+    points.col(i) = inputs[i].template head<3>().transpose();
+  }
+  Eigen::Spline<double, 3> spline3D =
+      Eigen::SplineFitting<Eigen::Spline<double, 3>>::Interpolate(points, 2);
+
+  std::vector<Acts::Vector3> output;
+  double step = 1. / (nPoints - 1);
+  for (unsigned int i = 0; i < nPoints; ++i) {
+    double t = i * step;
+    Eigen::Vector3d point = spline3D(t);
+    output.push_back(Acts::Vector3(point[0], point[1], point[2]));
+  }
+  return output;
+}
+
+}  // namespace
+
+ActsExamples::ObjSimHitWriter::ObjSimHitWriter(
+    const ActsExamples::ObjSimHitWriter::Config& config,
+    Acts::Logging::Level level)
+    : WriterT(config.inputSimHits, "ObjSimHitWriter", level), m_cfg(config) {
+  // inputSimHits is already checked by base constructor
+  if (m_cfg.outputStem.empty()) {
+    throw std::invalid_argument("Missing output filename stem");
+  }
+}
+
+ActsExamples::ProcessCode ActsExamples::ObjSimHitWriter::writeT(
+    const AlgorithmContext& ctx, const ActsExamples::SimHitContainer& simHits) {
+  // ensure exclusive access to tree/file while writing
+  std::lock_guard<std::mutex> lock(m_writeMutex);
+
+  // open per-event file for all simhit components
+  std::string pathSimHit = perEventFilepath(
+      m_cfg.outputDir, m_cfg.outputStem + ".obj", ctx.eventNumber);
+  std::string pathSimTrajectory = perEventFilepath(
+      m_cfg.outputDir, m_cfg.outputStem + "_trajectory.obj", ctx.eventNumber);
+
+  std::ofstream osHits(pathSimHit, std::ofstream::out | std::ofstream::trunc);
+  std::ofstream osTrajectory(pathSimTrajectory,
+                             std::ofstream::out | std::ofstream::trunc);
+
+  if (!osHits || !osTrajectory) {
+    throw std::ios_base::failure("Could not open '" + pathSimHit + "' or '" +
+                                 pathSimTrajectory + "' to write");
+  }
+
+  // Only hit plotting
+  if (!m_cfg.drawConnections) {
+    // Write data from internal immplementation
+    for (const auto& simHit : simHits) {
+      // local simhit information in global coord.
+      const Acts::Vector4& globalPos4 = simHit.fourPosition();
+
+      osHits << "v " << globalPos4[Acts::ePos0] / Acts::UnitConstants::mm << " "
+             << globalPos4[Acts::ePos1] / Acts::UnitConstants::mm << " "
+             << globalPos4[Acts::ePos2] / Acts::UnitConstants::mm << std::endl;
+
+    }  // end simHit loop
+  } else {
+    // We need to associate first
+    std::unordered_map<std::size_t, std::vector<Acts::Vector4>> particleHits;
+    // Pre-loop over hits ... write those below threshold
+    for (const auto& simHit : simHits) {
+      double momentum = simHit.momentum4Before().head<3>().norm();
+      if (momentum < m_cfg.momentumThreshold) {
+        ACTS_VERBOSE("Skipping : Hit below threshold: " << momentum);
+        continue;
+      } else if (momentum < m_cfg.momentumThresholdTraj) {
+        ACTS_VERBOSE(
+            "Skipping (trajectory): Hit below threshold: " << momentum);
+        osHits << "v "
+               << simHit.fourPosition()[Acts::ePos0] / Acts::UnitConstants::mm
+               << " "
+               << simHit.fourPosition()[Acts::ePos1] / Acts::UnitConstants::mm
+               << " "
+               << simHit.fourPosition()[Acts::ePos2] / Acts::UnitConstants::mm
+               << std::endl;
+        continue;
+      }
+      ACTS_VERBOSE("Accepting : Hit above threshold: " << momentum);
+
+      if (particleHits.find(simHit.particleId().value()) ==
+          particleHits.end()) {
+        particleHits[simHit.particleId().value()] = {};
+      }
+      particleHits[simHit.particleId().value()].push_back(
+          simHit.fourPosition());
+    }
+    // Draw loop
+    unsigned int lOffset = 1;
+    for (auto& [pId, pHits] : particleHits) {
+      // Draw the particle hits
+      std::sort(pHits.begin(), pHits.end(),
+                [](const Acts::Vector4& a, const Acts::Vector4& b) {
+                  return a[Acts::eTime] < b[Acts::eTime];
+                });
+
+      osHits << "o particle_" << pId << std::endl;
+      for (const auto& hit : pHits) {
+        osHits << "v " << hit[Acts::ePos0] / Acts::UnitConstants::mm << " "
+               << hit[Acts::ePos1] / Acts::UnitConstants::mm << " "
+               << hit[Acts::ePos2] / Acts::UnitConstants::mm << std::endl;
+      }
+      osHits << '\n';
+
+      // Interpolate the points
+      std::vector<Acts::Vector3> trajectory;
+      if (pHits.size() < 3) {
+        for (const auto& hit : pHits) {
+          trajectory.push_back(hit.template head<3>());
+        }
+      } else {
+        trajectory =
+            interpolatedPoints(pHits, pHits.size() * m_cfg.nInterpolatedPoints);
+      }
+
+      osTrajectory << "o particle_trajectory_" << pId << std::endl;
+      for (const auto& hit : trajectory) {
+        osTrajectory << "v " << hit[Acts::ePos0] / Acts::UnitConstants::mm
+                     << " " << hit[Acts::ePos1] / Acts::UnitConstants::mm << " "
+                     << hit[Acts::ePos2] / Acts::UnitConstants::mm << std::endl;
+      }
+      // Draw the line
+      for (unsigned int iv = lOffset + 1; iv < lOffset + trajectory.size();
+           ++iv) {
+        osTrajectory << "l " << iv - 1 << " " << iv << '\n';
+      }
+      osTrajectory << '\n';
+      // Increase the offset count
+      lOffset += trajectory.size();
+    }
+  }
+  osHits.close();
+  osTrajectory.close();
+
+  return ActsExamples::ProcessCode::SUCCESS;
+}
+
+ActsExamples::ProcessCode ActsExamples::ObjSimHitWriter::finalize() {
+  return ActsExamples::ProcessCode::SUCCESS;
+}

Examples/Io/Obj/src/ObjTrackingGeometryWriter.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/ObjTrackingGeometryWriter.hpp"
+#include "ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp"
 
 #include "Acts/Utilities/Logger.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"

Examples/Python/python/acts/examples/simulation.py

@@ -420,6 +420,7 @@ def addFatras(
     outputSimHits: str = "simhits",
     outputDirCsv: Optional[Union[Path, str]] = None,
     outputDirRoot: Optional[Union[Path, str]] = None,
+    outputDirObj: Optional[Union[Path, str]] = None,
     logLevel: Optional[acts.logging.Level] = None,
 ) -> None:
     """This function steers the detector simulation using Fatras
@@ -444,6 +445,8 @@ def addFatras(
         the output folder for the Csv output, None triggers no output
     outputDirRoot : Path|str, path, None
         the output folder for the Root output, None triggers no output
+    outputDirObj : Path|str, path, None
+        the output folder for the Obj output, None triggers no output
     """
 
     customLogLevel = acts.examples.defaultLogging(s, logLevel)
@@ -507,6 +510,7 @@ def addFatras(
         particlesPostSelected,
         outputDirCsv,
         outputDirRoot,
+        outputDirObj,
         logLevel,
     )
 
@@ -519,6 +523,7 @@ def addSimWriters(
     particlesSimulated: str = "particles_simulated",
     outputDirCsv: Optional[Union[Path, str]] = None,
     outputDirRoot: Optional[Union[Path, str]] = None,
+    outputDirObj: Optional[Union[Path, str]] = None,
     logLevel: Optional[acts.logging.Level] = None,
 ) -> None:
     customLogLevel = acts.examples.defaultLogging(s, logLevel)
@@ -563,6 +568,19 @@ def addSimWriters(
             )
         )
 
+    if outputDirObj is not None:
+        outputDirObj = Path(outputDirObj)
+        if not outputDirObj.exists():
+            outputDirObj.mkdir()
+        s.addWriter(
+            acts.examples.ObjSimHitWriter(
+                level=customLogLevel(),
+                inputSimHits=simHits,
+                outputDir=str(outputDirObj),
+                outputStem="hits",
+            )
+        )
+
 
 def getG4DetectorConstructionFactory(
     detector: Any,
@@ -620,6 +638,7 @@ def addGeant4(
     keepParticlesWithoutHits=True,
     outputDirCsv: Optional[Union[Path, str]] = None,
     outputDirRoot: Optional[Union[Path, str]] = None,
+    outputDirObj: Optional[Union[Path, str]] = None,
     logLevel: Optional[acts.logging.Level] = None,
     killVolume: Optional[acts.Volume] = None,
     killAfterTime: float = float("inf"),
@@ -647,6 +666,8 @@ def addGeant4(
         the output folder for the Csv output, None triggers no output
     outputDirRoot : Path|str, path, None
         the output folder for the Root output, None triggers no output
+    outputDirObj : Path|str, path, None
+        the output folder for the Obj output, None triggers no output
     killVolume: acts.Volume, None
         if given, particles are killed when going outside this volume.
     killAfterTime: float
@@ -736,7 +757,8 @@ def addGeant4(
         particlesPostSelected,
         outputDirCsv,
         outputDirRoot,
-        logLevel,
+        outputDirObj,
+        logLevel=logLevel,
     )
 
     return s