From 7efad56ce50ed433ac5ea14dcef2a51d8cc3dce4 Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Wed, 4 Dec 2024 19:08:28 +0100 Subject: [PATCH] feat: introduction obj simhit writer (#3180) This PR adds a wavefront `.obj` sim hit writer to the Examples. Added (optional) spline fit based on `Eigen::Spline` In addition it corrects the include path of Obj based writers to `Io/Ob` it still was legacy `Plugins/Obj`. ## Summary by CodeRabbit - **New Features** - Introduced `ObjSimHitWriter` class for writing simulation hit data to OBJ files. - Added optional parameter `outputDirObj` to several simulation functions for enhanced output management. - Enhanced argument parser in `full_chain_odd.py` for more granular control over particle gun configurations. - **Bug Fixes** - Improved error handling during file operations in the `ObjSimHitWriter`. - **Chores** - Restructured file organization for better clarity and modularity. - Removed the `sim_digi_odd.py` script as part of cleanup. ![ttbar-mu20000](https://github.com/acts-project/acts/assets/26623879/e65d33c3-2223-4584-8ec7-43a38ccbf172) ![Screenshot 2024-05-13 at 14 54 35](https://github.com/acts-project/acts/assets/26623879/4132b4e8-4d22-4bc3-856e-ffc3715a2ff3) --- Examples/Io/Obj/CMakeLists.txt | 1 + .../Obj/ObjPropagationStepsWriter.hpp | 0 .../ActsExamples/Io/Obj/ObjSimHitWriter.hpp | 84 ++++++++ .../Obj/ObjTrackingGeometryWriter.hpp | 0 .../{Plugins => Io}/Obj/ObjWriterOptions.hpp | 0 .../Io/Obj/src/ObjPropagationStepsWriter.cpp | 2 +- Examples/Io/Obj/src/ObjSimHitWriter.cpp | 190 ++++++++++++++++++ .../Io/Obj/src/ObjTrackingGeometryWriter.cpp | 2 +- .../Python/python/acts/examples/simulation.py | 24 ++- Examples/Python/src/Output.cpp | 9 +- Examples/Scripts/Python/full_chain_odd.py | 8 + Examples/Scripts/Python/geant4.py | 1 + 12 files changed, 316 insertions(+), 5 deletions(-) rename Examples/Io/Obj/include/ActsExamples/{Plugins => Io}/Obj/ObjPropagationStepsWriter.hpp (100%) create mode 100644 Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp rename Examples/Io/Obj/include/ActsExamples/{Plugins => Io}/Obj/ObjTrackingGeometryWriter.hpp (100%) rename Examples/Io/Obj/include/ActsExamples/{Plugins => Io}/Obj/ObjWriterOptions.hpp (100%) create mode 100644 Examples/Io/Obj/src/ObjSimHitWriter.cpp diff --git a/Examples/Io/Obj/CMakeLists.txt b/Examples/Io/Obj/CMakeLists.txt index 9ca888f726a..bce7bf19e6b 100644 --- a/Examples/Io/Obj/CMakeLists.txt +++ b/Examples/Io/Obj/CMakeLists.txt @@ -2,6 +2,7 @@ add_library( ActsExamplesIoObj SHARED src/ObjTrackingGeometryWriter.cpp + src/ObjSimHitWriter.cpp src/ObjPropagationStepsWriter.cpp ) diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp similarity index 100% rename from Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp rename to Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp diff --git a/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp new file mode 100644 index 00000000000..8ccd92645b9 --- /dev/null +++ b/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 +#include +#include + +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-.obj +/// event000000001-_trajectory.obj +/// event000000002-.obj +/// event000000002-_trajectory.obj +/// +class ObjSimHitWriter : public WriterT { + 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::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 diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp similarity index 100% rename from Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp rename to Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjWriterOptions.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjWriterOptions.hpp similarity index 100% rename from Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjWriterOptions.hpp rename to Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjWriterOptions.hpp diff --git a/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp b/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp index 6d0388984cc..0f196ce162c 100644 --- a/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp +++ b/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 { diff --git a/Examples/Io/Obj/src/ObjSimHitWriter.cpp b/Examples/Io/Obj/src/ObjSimHitWriter.cpp new file mode 100644 index 00000000000..564ff9c5409 --- /dev/null +++ b/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 +#include +#include +#include + +#include + +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 interpolated points +template +std::vector interpolatedPoints( + const std::vector& inputs, std::size_t nPoints) { + std::vector 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 spline3D = + Eigen::SplineFitting>::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> 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 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; +} diff --git a/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp b/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp index 4937757c85f..d2dbe52e745 100644 --- a/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp +++ b/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" diff --git a/Examples/Python/python/acts/examples/simulation.py b/Examples/Python/python/acts/examples/simulation.py index b60b05466a0..d076a79b30b 100644 --- a/Examples/Python/python/acts/examples/simulation.py +++ b/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 diff --git a/Examples/Python/src/Output.cpp b/Examples/Python/src/Output.cpp index 4281f2db26e..251ca2b6357 100644 --- a/Examples/Python/src/Output.cpp +++ b/Examples/Python/src/Output.cpp @@ -25,6 +25,9 @@ #include "ActsExamples/Io/Csv/CsvTrackParameterWriter.hpp" #include "ActsExamples/Io/Csv/CsvTrackWriter.hpp" #include "ActsExamples/Io/Csv/CsvTrackingGeometryWriter.hpp" +#include "ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp" +#include "ActsExamples/Io/Obj/ObjSimHitWriter.hpp" +#include "ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp" #include "ActsExamples/Io/Root/RootBFieldWriter.hpp" #include "ActsExamples/Io/Root/RootMaterialTrackWriter.hpp" #include "ActsExamples/Io/Root/RootMaterialWriter.hpp" @@ -46,8 +49,6 @@ #include "ActsExamples/Io/Root/TrackFitterPerformanceWriter.hpp" #include "ActsExamples/Io/Root/VertexNTupleWriter.hpp" #include "ActsExamples/MaterialMapping/IMaterialWriter.hpp" -#include "ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp" -#include "ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp" #include "ActsExamples/TrackFinding/ITrackParamsLookupReader.hpp" #include "ActsExamples/TrackFinding/ITrackParamsLookupWriter.hpp" @@ -111,6 +112,10 @@ void addOutput(Context& ctx) { "ObjPropagationStepsWriter", collection, outputDir, outputScalor, outputPrecision); + ACTS_PYTHON_DECLARE_WRITER(ActsExamples::ObjSimHitWriter, mex, + "ObjSimHitWriter", inputSimHits, outputDir, + outputStem, outputPrecision, drawConnections); + { auto c = py::class_(m, "ViewConfig").def(py::init<>()); diff --git a/Examples/Scripts/Python/full_chain_odd.py b/Examples/Scripts/Python/full_chain_odd.py index 8b86402b820..7ddaa02054f 100755 --- a/Examples/Scripts/Python/full_chain_odd.py +++ b/Examples/Scripts/Python/full_chain_odd.py @@ -136,6 +136,12 @@ default=True, action=argparse.BooleanOptionalAction, ) +parser.add_argument( + "--output-obj", + help="Switch obj output on/off", + default=True, + action=argparse.BooleanOptionalAction, +) args = parser.parse_args() @@ -277,6 +283,7 @@ ), outputDirRoot=outputDir if args.output_root else None, outputDirCsv=outputDir if args.output_csv else None, + outputDirObj=outputDir if args.output_obj else None, rnd=rnd, killVolume=trackingGeometry.highestTrackingVolume, killAfterTime=25 * u.ns, @@ -305,6 +312,7 @@ enableInteractions=True, outputDirRoot=outputDir if args.output_root else None, outputDirCsv=outputDir if args.output_csv else None, + outputDirObj=outputDir if args.output_obj else None, rnd=rnd, ) diff --git a/Examples/Scripts/Python/geant4.py b/Examples/Scripts/Python/geant4.py index 42cb5db95bb..65e154c2704 100755 --- a/Examples/Scripts/Python/geant4.py +++ b/Examples/Scripts/Python/geant4.py @@ -36,6 +36,7 @@ def runGeant4( field, outputDirCsv=outputDir / "csv", outputDirRoot=outputDir, + outputDirObj=outputDir / "obj", rnd=rnd, materialMappings=materialMappings, volumeMappings=volumeMappings,