Merge branch 'ex-digi-particle-measurement-maps' of github.com:andiwa…

…nd/acts into ex-digi-particle-measurement-maps

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 {
 

Examples/Io/Obj/src/ObjSimHitWriter.cpp

@@ -0,0 +1,190 @@
+// 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 "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 {
+
+/// @brief Helper function to interpolate points
+///
+/// @tparam input_vector_type
+/// @param inputs input vector points
+/// @param nPoints number of interpolation points
+///
+/// @return std::vector<Acts::Vector3> interpolated points
+template <typename input_vector_type>
+std::vector<Acts::Vector3> interpolatedPoints(
+    const std::vector<input_vector_type>& inputs, std::size_t nPoints) {
+  std::vector<Acts::Vector3> output;
+
+  if (nPoints < 2) {
+    // No interpolation done return simply the output vector
+    for (const auto& input : inputs) {
+      output.push_back(input.template head<3>());
+    }
+
+  } else {
+    Eigen::MatrixXd points(3, inputs.size());
+    for (std::size_t 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);
+
+    double step = 1. / (nPoints - 1);
+    for (std::size_t i = 0; i < nPoints; ++i) {
+      double t = i * step;
+      output.push_back(spline3D(t));
+    }
+  }
+  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::scoped_lock 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
+    std::size_t 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 (std::size_t 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
@@ -737,7 +758,8 @@ def addGeant4(
         particlesPostSelected,
         outputDirCsv,
         outputDirRoot,
-        logLevel,
+        outputDirObj,
+        logLevel=logLevel,
     )
 
     return s