diff --git a/Examples/Framework/include/ActsExamples/EventData/AverageSimHits.hpp b/Examples/Framework/include/ActsExamples/EventData/AverageSimHits.hpp index 8191c0b6987..a88ca90aea9 100644 --- a/Examples/Framework/include/ActsExamples/EventData/AverageSimHits.hpp +++ b/Examples/Framework/include/ActsExamples/EventData/AverageSimHits.hpp @@ -9,6 +9,7 @@ #pragma once #include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Logger.hpp" @@ -51,11 +52,19 @@ inline std::tuple averageSimHits( // check their validity again. const auto& simHit = *simHits.nth(simHitIdx); + // We use the thickness of the detector element as tolerance, because Geant4 + // treats the Surfaces as volumes and thus it is not ensured, that each hit + // lies exactly on the Acts::Surface + const auto tolerance = + surface.associatedDetectorElement() != nullptr + ? surface.associatedDetectorElement()->thickness() + : Acts::s_onSurfaceTolerance; + // transforming first to local positions and average that ensures that the // averaged position is still on the surface. the averaged global position // might not be on the surface anymore. auto result = surface.globalToLocal(gCtx, simHit.position(), - simHit.direction(), 0.5_um); + simHit.direction(), tolerance); if (result.ok()) { avgLocal += result.value(); } else { diff --git a/Examples/Io/EDM4hep/CMakeLists.txt b/Examples/Io/EDM4hep/CMakeLists.txt index 26e61adf5c9..9f7ddf37057 100644 --- a/Examples/Io/EDM4hep/CMakeLists.txt +++ b/Examples/Io/EDM4hep/CMakeLists.txt @@ -3,11 +3,10 @@ add_library( src/EDM4hepMeasurementReader.cpp src/EDM4hepMeasurementWriter.cpp src/EDM4hepMultiTrajectoryWriter.cpp + src/EDM4hepReader.cpp src/EDM4hepTrackWriter.cpp src/EDM4hepTrackReader.cpp - src/EDM4hepParticleReader.cpp src/EDM4hepParticleWriter.cpp - src/EDM4hepSimHitReader.cpp src/EDM4hepSimHitWriter.cpp src/EDM4hepUtil.cpp) target_include_directories( diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp index bf46a34009d..0029740f80e 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp @@ -18,6 +18,7 @@ #include #include +#include namespace ActsExamples { @@ -68,7 +69,9 @@ class EDM4hepMeasurementReader final : public IReader { std::pair m_eventsRange; std::unique_ptr m_logger; - podio::ROOTFrameReader m_reader; + tbb::enumerable_thread_specific m_reader; + + podio::ROOTFrameReader& reader(); WriteDataHandle m_outputMeasurements{ this, "OutputMeasurements"}; diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp index 27ca66b4184..f6914a4e8a1 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp @@ -66,6 +66,8 @@ class EDM4hepMeasurementWriter final : public WriterT { podio::ROOTFrameWriter m_writer; + std::mutex m_writeMutex; + ReadDataHandle m_inputClusters{this, "InputClusters"}; }; diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMultiTrajectoryWriter.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMultiTrajectoryWriter.hpp index 5bb1384b51d..a0d4e27da87 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMultiTrajectoryWriter.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMultiTrajectoryWriter.hpp @@ -68,6 +68,7 @@ class EDM4hepMultiTrajectoryWriter : public WriterT { private: Config m_cfg; + std::mutex m_writeMutex; podio::ROOTFrameWriter m_writer; ReadDataHandle> diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepParticleReader.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepParticleReader.hpp deleted file mode 100644 index cac002e20a8..00000000000 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepParticleReader.hpp +++ /dev/null @@ -1,68 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2022 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/Utilities/Logger.hpp" -#include "ActsExamples/EventData/SimParticle.hpp" -#include "ActsExamples/Framework/DataHandle.hpp" -#include "ActsExamples/Framework/IReader.hpp" - -#include -#include - -#include -#include - -namespace ActsExamples { - -/// Read particles from EDM4hep. -/// -/// Inpersistent information: -/// - particle ID -/// - process -class EDM4hepParticleReader final : public IReader { - public: - struct Config { - /// Where to read input file from. - std::string inputPath; - /// Name of the particle collection in EDM4hep. - std::string inputParticles = "MCParticles"; - /// Which particle collection to read into. - std::string outputParticles; - }; - - /// Construct the particle reader. - /// - /// @param config is the configuration object - /// @param level is the logging level - EDM4hepParticleReader(const Config& config, Acts::Logging::Level level); - - std::string name() const final; - - /// Return the available events range. - std::pair availableEvents() const final; - - /// Read out data from the input stream. - ProcessCode read(const ActsExamples::AlgorithmContext& ctx) final; - - /// Readonly access to the config - const Config& config() const { return m_cfg; } - - private: - Config m_cfg; - std::pair m_eventsRange; - std::unique_ptr m_logger; - - podio::ROOTFrameReader m_reader; - - WriteDataHandle m_outputParticles{this, - "OutputParticles"}; -}; - -} // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepParticleWriter.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepParticleWriter.hpp index 31582497819..1fb59a067e3 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepParticleWriter.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepParticleWriter.hpp @@ -11,6 +11,7 @@ #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/Framework/WriterT.hpp" +#include #include #include @@ -55,6 +56,8 @@ class EDM4hepParticleWriter final : public WriterT { private: Config m_cfg; + std::mutex m_writeMutex; + podio::ROOTFrameWriter m_writer; }; diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepReader.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepReader.hpp new file mode 100644 index 00000000000..eadeec143ae --- /dev/null +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepReader.hpp @@ -0,0 +1,128 @@ +// 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/TrackingGeometry.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/EventData/SimParticle.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IReader.hpp" +#include "ActsFatras/EventData/Particle.hpp" + +#include +#include + +#include +#include +#include +#include + +namespace ActsExamples { + +namespace DD4hep { +struct DD4hepDetector; +} + +/// Read particles from EDM4hep. +/// +/// Inpersistent information: +/// - particle ID +/// - process +class EDM4hepReader final : public IReader { + public: + struct Config { + /// Where to read input file from. + std::string inputPath; + /// Name of the particle collection in EDM4hep. + std::string inputParticles = "MCParticles"; + /// Names of the sim hit collections + std::vector inputSimHits{}; + /// Particles at creation + std::string outputParticlesInitial; + /// Particles at their endpoints + std::string outputParticlesFinal; + /// Particles from the generator + std::string outputParticlesGenerator; + + /// Output simulated (truth) hits collection. + std::string outputSimHits; + + /// Directory into which to write graphviz files for particles + /// Empty string means no output + std::string graphvizOutput = ""; + + /// DD4hep detector for cellID resolution. + std::shared_ptr dd4hepDetector; + + /// Tracking geometry for cellID resolution. + std::shared_ptr trackingGeometry; + + /// Whether to sort sim hits in time to produce index sequence + bool sortSimHitsInTime = false; + }; + + using ParentRelationship = std::unordered_map; + + /// Construct the particle reader. + /// + /// @param config is the configuration object + /// @param level is the logging level + EDM4hepReader(const Config& config, Acts::Logging::Level level); + + std::string name() const final; + + /// Return the available events range. + std::pair availableEvents() const final; + + /// Read out data from the input stream. + ProcessCode read(const ActsExamples::AlgorithmContext& ctx) final; + + /// Readonly access to the config + const Config& config() const { return m_cfg; } + + void processChildren(const edm4hep::MCParticle& particle, SimBarcode parentId, + SimParticleContainer::sequence_type& particles, + ParentRelationship& parentRelationship, + std::unordered_map& particleMap, + std::size_t& nSecondaryVertices, + std::size_t& maxGen) const; + + static void setSubParticleIds( + const SimParticleContainer::sequence_type::iterator& begin, + const SimParticleContainer::sequence_type::iterator& end); + + private: + const Acts::Logger& logger() const { return *m_logger; } + + Config m_cfg; + std::pair m_eventsRange; + std::unique_ptr m_logger; + + std::unordered_map m_surfaceMap; + + tbb::enumerable_thread_specific m_reader; + + podio::ROOTFrameReader& reader(); + + WriteDataHandle m_outputParticlesInitial{ + this, "OutputParticlesInitial"}; + WriteDataHandle m_outputParticlesFinal{ + this, "OutputParticlesFinal"}; + WriteDataHandle m_outputParticlesGenerator{ + this, "OutputParticlesGenerator"}; + + WriteDataHandle m_outputSimHits{this, "OutputSimHits"}; + + void graphviz(std::ostream& os, + const SimParticleContainer::sequence_type& particles, + const ParentRelationship& parents) const; +}; + +} // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepSimHitReader.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepSimHitReader.hpp deleted file mode 100644 index 31d45078612..00000000000 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepSimHitReader.hpp +++ /dev/null @@ -1,81 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2022 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/Utilities/Logger.hpp" -#include "ActsExamples/DD4hepDetector/DD4hepDetector.hpp" -#include "ActsExamples/EventData/SimHit.hpp" -#include "ActsExamples/EventData/SimParticle.hpp" -#include "ActsExamples/Framework/DataHandle.hpp" -#include "ActsExamples/Framework/IReader.hpp" - -#include -#include - -#include -#include - -namespace ActsExamples { - -/// Read in a simhit collection from EDM4hep. -/// -/// Inpersistent information: -/// - particle ID -/// - after4 momentum -/// - hit index -/// - digitization channel -class EDM4hepSimHitReader final : public IReader { - public: - struct Config { - /// Where to the read input file from. - std::string inputPath; - /// Name of the particle collection in EDM4hep. - std::string inputParticles = "MCParticles"; - /// Name of the sim tracker hit collection in EDM4hep - std::string inputSimTrackerHits = "ActsSimTrackerHits"; - /// Which particle collection to read into. - std::string outputParticles; - /// Output simulated (truth) hits collection. - std::string outputSimHits; - /// DD4hep detector for cellID resolution. - std::shared_ptr dd4hepDetector; - }; - - /// Construct the simhit reader. - /// - /// @param config is the configuration object - /// @param level is the logging level - EDM4hepSimHitReader(const Config& config, Acts::Logging::Level level); - - std::string name() const final; - - /// Return the available events range. - std::pair availableEvents() const final; - - /// Read out data from the input stream. - ProcessCode read(const ActsExamples::AlgorithmContext& ctx) final; - - /// Readonly access to the config - const Config& config() const { return m_cfg; } - - private: - Config m_cfg; - std::pair m_eventsRange; - std::unique_ptr m_logger; - - podio::ROOTFrameReader m_reader; - - const Acts::Logger& logger() const { return *m_logger; } - - WriteDataHandle m_outputSimHits{this, "OutputSimHits"}; - WriteDataHandle m_outputParticles{this, - "OutputParticles"}; -}; - -} // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepSimHitWriter.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepSimHitWriter.hpp index e33fb1fb411..db2b2e92a3b 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepSimHitWriter.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepSimHitWriter.hpp @@ -67,6 +67,8 @@ class EDM4hepSimHitWriter final : public WriterT { podio::ROOTFrameWriter m_writer; + std::mutex m_writeMutex; + ReadDataHandle m_inputParticles{this, "InputParticles"}; }; diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepTrackReader.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepTrackReader.hpp index 57ce35af7a8..556f809a3bf 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepTrackReader.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepTrackReader.hpp @@ -19,6 +19,7 @@ #include #include +#include namespace ActsExamples { @@ -54,11 +55,15 @@ class EDM4hepTrackReader : public IReader { ProcessCode read(const ActsExamples::AlgorithmContext& ctx) final; private: + std::pair m_eventsRange; + Config m_cfg; WriteDataHandle m_outputTracks{this, "OutputTracks"}; - podio::ROOTFrameReader m_reader; + tbb::enumerable_thread_specific m_reader; + + podio::ROOTFrameReader& reader(); std::unique_ptr m_logger; diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepTrackWriter.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepTrackWriter.hpp index f98e13a1456..624abbea306 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepTrackWriter.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepTrackWriter.hpp @@ -56,6 +56,8 @@ class EDM4hepTrackWriter : public WriterT { private: Config m_cfg; + std::mutex m_writeMutex; + podio::ROOTFrameWriter m_writer; }; diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp index 11b731125bc..38e12be6306 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp @@ -31,10 +31,15 @@ namespace ActsExamples { namespace EDM4hepUtil { using MapParticleIdFrom = - std::function; + std::function; using MapParticleIdTo = std::function; +inline ActsFatras::Barcode zeroParticleMapper( + const edm4hep::MCParticle& /*particle*/) { + return 0; +} + using MapGeometryIdFrom = std::function; using MapGeometryIdTo = @@ -45,8 +50,9 @@ using MapGeometryIdTo = /// Inpersistent information: /// - particle ID /// - process -ActsFatras::Particle readParticle(const edm4hep::MCParticle& from, - const MapParticleIdFrom& particleMapper); +ActsFatras::Particle readParticle( + const edm4hep::MCParticle& from, + const MapParticleIdFrom& particleMapper = zeroParticleMapper); /// Write a Fatras particle into EDM4hep. /// diff --git a/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp index 17fe4164747..66fbcce36b1 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp @@ -33,9 +33,7 @@ EDM4hepMeasurementReader::EDM4hepMeasurementReader( throw std::invalid_argument("Missing measurement output collection"); } - m_reader.openFile(m_cfg.inputPath); - - m_eventsRange = std::make_pair(0, m_reader.getEntries("events")); + m_eventsRange = std::make_pair(0, reader().getEntries("events")); m_outputMeasurements.initialize(m_cfg.outputMeasurements); m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap); @@ -59,7 +57,7 @@ ProcessCode EDM4hepMeasurementReader::read(const AlgorithmContext& ctx) { IndexMultimap measurementSimHitsMap; IndexSourceLinkContainer sourceLinks; - podio::Frame frame = m_reader.readEntry("events", ctx.eventNumber); + podio::Frame frame = reader().readEntry("events", ctx.eventNumber); const auto& trackerHitPlaneCollection = frame.get("ActsTrackerHitsPlane"); @@ -87,4 +85,14 @@ ProcessCode EDM4hepMeasurementReader::read(const AlgorithmContext& ctx) { return ProcessCode::SUCCESS; } +podio::ROOTFrameReader& EDM4hepMeasurementReader::reader() { + bool exists = false; + auto& reader = m_reader.local(exists); + if (!exists) { + reader.openFile(m_cfg.inputPath); + } + + return reader; +} + } // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/src/EDM4hepMeasurementWriter.cpp b/Examples/Io/EDM4hep/src/EDM4hepMeasurementWriter.cpp index 501a416d63e..35f93cf97fe 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepMeasurementWriter.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepMeasurementWriter.cpp @@ -67,6 +67,7 @@ ActsExamples::ProcessCode EDM4hepMeasurementWriter::writeT( frame.put(std::move(hitsPlane), "ActsTrackerHitsPlane"); frame.put(std::move(hits), "ActsTrackerHitsRaw"); + std::lock_guard guard(m_writeMutex); m_writer.writeFrame(frame, "events"); return ActsExamples::ProcessCode::SUCCESS; diff --git a/Examples/Io/EDM4hep/src/EDM4hepMultiTrajectoryWriter.cpp b/Examples/Io/EDM4hep/src/EDM4hepMultiTrajectoryWriter.cpp index 0481ac4bda0..bfd48049c8c 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepMultiTrajectoryWriter.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepMultiTrajectoryWriter.cpp @@ -63,6 +63,7 @@ ProcessCode EDM4hepMultiTrajectoryWriter::writeT( frame.put(std::move(trackCollection), "ActsTracks"); + std::lock_guard guard(m_writeMutex); m_writer.writeFrame(frame, "events"); return ProcessCode::SUCCESS; diff --git a/Examples/Io/EDM4hep/src/EDM4hepParticleReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepParticleReader.cpp deleted file mode 100644 index 96c6d5572b2..00000000000 --- a/Examples/Io/EDM4hep/src/EDM4hepParticleReader.cpp +++ /dev/null @@ -1,72 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2022 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/EDM4hep/EDM4hepParticleReader.hpp" - -#include "Acts/Definitions/Units.hpp" -#include "ActsExamples/EventData/SimParticle.hpp" -#include "ActsExamples/Framework/WhiteBoard.hpp" -#include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp" -#include "ActsExamples/Utilities/Paths.hpp" - -#include - -#include - -namespace ActsExamples { - -EDM4hepParticleReader::EDM4hepParticleReader( - const EDM4hepParticleReader::Config& config, Acts::Logging::Level level) - : m_cfg(config), - m_logger(Acts::getDefaultLogger("EDM4hepParticleReader", level)) { - if (m_cfg.outputParticles.empty()) { - throw std::invalid_argument("Missing output collection"); - } - - m_reader.openFile(m_cfg.inputPath); - - m_eventsRange = std::make_pair(0, m_reader.getEntries("events")); - - m_outputParticles.initialize(m_cfg.outputParticles); -} - -std::string EDM4hepParticleReader::name() const { - return "EDM4hepParticleReader"; -} - -std::pair EDM4hepParticleReader::availableEvents() - const { - return m_eventsRange; -} - -ProcessCode EDM4hepParticleReader::read(const AlgorithmContext& ctx) { - podio::Frame frame = m_reader.readEntry("events", ctx.eventNumber); - const auto& mcParticleCollection = - frame.get(m_cfg.inputParticles); - - SimParticleContainer::sequence_type unordered; - - for (const auto& mcParticle : mcParticleCollection) { - auto particle = - EDM4hepUtil::readParticle(mcParticle, [](const edm4hep::MCParticle& p) { - ActsFatras::Barcode result; - result.setParticle(EDM4hepUtil::podioObjectIDToInteger(p.id())); - return result; - }); - unordered.push_back(particle); - } - - // Write ordered particles container to the EventStore - SimParticleContainer particles; - particles.insert(unordered.begin(), unordered.end()); - m_outputParticles(ctx, std::move(particles)); - - return ProcessCode::SUCCESS; -} - -} // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/src/EDM4hepParticleWriter.cpp b/Examples/Io/EDM4hep/src/EDM4hepParticleWriter.cpp index 9d64576b6fd..c4b5533e599 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepParticleWriter.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepParticleWriter.cpp @@ -52,6 +52,7 @@ ProcessCode EDM4hepParticleWriter::writeT( frame.put(std::move(mcParticleCollection), m_cfg.outputParticles); + std::lock_guard guard(m_writeMutex); m_writer.writeFrame(frame, "events"); return ProcessCode::SUCCESS; diff --git a/Examples/Io/EDM4hep/src/EDM4hepReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepReader.cpp new file mode 100644 index 00000000000..80aa321f90d --- /dev/null +++ b/Examples/Io/EDM4hep/src/EDM4hepReader.cpp @@ -0,0 +1,515 @@ +// 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/EDM4hep/EDM4hepReader.hpp" + +#include "Acts/Definitions/Units.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp" +#include "ActsExamples/DD4hepDetector/DD4hepDetector.hpp" +#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/EventData/SimParticle.hpp" +#include "ActsExamples/Framework/WhiteBoard.hpp" +#include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp" +#include "ActsExamples/Utilities/Paths.hpp" +#include "ActsFatras/EventData/Barcode.hpp" + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace ActsExamples { + +EDM4hepReader::EDM4hepReader(const Config& config, Acts::Logging::Level level) + : m_cfg(config), + m_logger(Acts::getDefaultLogger("EDM4hepParticleReader", level)) { + if (m_cfg.outputParticlesInitial.empty()) { + throw std::invalid_argument("Missing output collection initial particles"); + } + + if (m_cfg.outputParticlesFinal.empty()) { + throw std::invalid_argument("Missing output collection final particles"); + } + + if (m_cfg.outputParticlesGenerator.empty()) { + throw std::invalid_argument( + "Missing output collection generator particles"); + } + + if (m_cfg.outputSimHits.empty()) { + throw std::invalid_argument("Missing output collection sim hits"); + } + + m_eventsRange = std::make_pair(0, reader().getEntries("events")); + + m_outputParticlesInitial.initialize(m_cfg.outputParticlesInitial); + m_outputParticlesFinal.initialize(m_cfg.outputParticlesFinal); + m_outputParticlesGenerator.initialize(m_cfg.outputParticlesGenerator); + m_outputSimHits.initialize(m_cfg.outputSimHits); + + m_cfg.trackingGeometry->visitSurfaces([&](const auto* surface) { + const auto* detElement = dynamic_cast( + surface->associatedDetectorElement()); + + if (detElement == nullptr) { + ACTS_ERROR("Surface has no associated detector element"); + return; + } + + const auto translation = detElement->sourceElement() + .nominal() + .worldTransformation() + .GetTranslation(); + Acts::Vector3 position; + position << translation[0], translation[1], translation[2]; + position *= Acts::UnitConstants::cm; + + m_surfaceMap.insert({detElement->sourceElement().key(), surface}); + }); +} + +podio::ROOTFrameReader& EDM4hepReader::reader() { + bool exists = false; + auto& reader = m_reader.local(exists); + if (!exists) { + reader.openFile(m_cfg.inputPath); + } + + return reader; +} + +std::string EDM4hepReader::name() const { + return "EDM4hepReader"; +} + +std::pair EDM4hepReader::availableEvents() const { + return m_eventsRange; +} + +namespace { +std::string vid(unsigned int vtx) { + return "V" + std::to_string(vtx); +} + +std::string pid(const SimParticle& particle) { + return "P" + std::to_string(particle.particleId().value()); +} + +std::string plabel(const SimParticle& particle) { + using namespace Acts::UnitLiterals; + std::stringstream ss; + ss << particle.pdg() << "\\n(" << particle.particleId() << ")\\n" + << "p=" << std::setprecision(3) << particle.absoluteMomentum() / 1_GeV + << " GeV"; + return ss.str(); +} + +} // namespace + +void EDM4hepReader::graphviz( + std::ostream& os, const SimParticleContainer::sequence_type& particles, + const ParentRelationship& parents) const { + os << "digraph Event {\n"; + + std::set primaryVertices; + + for (const auto& particle : particles) { + if (particle.particleId().generation() == 0) { + primaryVertices.insert(particle.particleId().vertexPrimary()); + + os << vid(particle.particleId().vertexPrimary()) << " -> " + << pid(particle) << ";\n"; + } + + os << pid(particle) << " [label=\"" << plabel(particle) << "\"];\n"; + } + + for (const auto [childIdx, parentIdx] : parents) { + const auto& child = particles[childIdx]; + const auto& parent = particles[parentIdx]; + os << pid(parent) << " -> " << pid(child); + + if (parent.particleId().vertexSecondary() == + child.particleId().vertexSecondary()) { + os << " [style=dashed]"; + } + + os << ";\n"; + } + + for (unsigned int vtx : primaryVertices) { + os << vid(vtx) << " [label=\"PV" << vtx << "\"];\n"; + } + + os << "}"; +} + +ProcessCode EDM4hepReader::read(const AlgorithmContext& ctx) { + podio::Frame frame = reader().readEntry("events", ctx.eventNumber); + const auto& mcParticleCollection = + frame.get(m_cfg.inputParticles); + + ACTS_DEBUG("Reading EDM4hep inputs"); + + SimParticleContainer::sequence_type unordered; + + // Read particles from the input file + // Find particles without parents and group them by vtx position to find + // primary vertices + std::vector>> + primaryVertices; + for (const auto& particle : mcParticleCollection) { + if (particle.parents_size() > 0) { + // not a primary vertex + continue; + } + const auto& vtx = particle.getVertex(); + Acts::Vector3 vtxPos = {vtx[0], vtx[1], vtx[2]}; + vtxPos /= Acts::UnitConstants::mm; + + // linear search for vector + auto it = std::find_if( + primaryVertices.begin(), primaryVertices.end(), + [&vtxPos]( + const std::pair>& + pair) { return pair.first == vtxPos; }); + + if (it == primaryVertices.end()) { + ACTS_DEBUG("Found primary vertex at " << vtx.x << ", " << vtx.y << ", " + << vtx.z); + primaryVertices.push_back({vtxPos, {particle}}); + } else { + it->second.push_back(particle); + } + } + + ACTS_DEBUG("Found " << primaryVertices.size() << " primary vertices"); + + // key: child, value: parent + ParentRelationship parentRelationship; + + // key: input particle index, value: index in the unordered particle + // container + std::unordered_map edm4hepParticleMap; + + std::size_t nPrimaryVertices = 0; + // Walk the particle tree + for (const auto& [vtxPos, particles] : primaryVertices) { + nPrimaryVertices += 1; + ACTS_DEBUG("Walking particle tree for primary vertex at " + << vtxPos.x() << ", " << vtxPos.y() << ", " << vtxPos.z()); + std::size_t nParticles = 0; + std::size_t nSecondaryVertices = 0; + std::size_t maxGen = 0; + auto startSize = unordered.size(); + for (const auto& inParticle : particles) { + nParticles += 1; + SimParticle particle{EDM4hepUtil::readParticle(inParticle)}; + particle.setParticleId(SimBarcode{} + .setParticle(nParticles) + .setVertexPrimary(nPrimaryVertices)); + ACTS_VERBOSE("+ add particle " << particle); + ACTS_VERBOSE(" - at " << particle.position().transpose()); + ACTS_VERBOSE(" - createdInSim: " << inParticle.isCreatedInSimulation()); + ACTS_VERBOSE(" - vertexIsNotEndpointOfParent: " + << inParticle.vertexIsNotEndpointOfParent()); + ACTS_VERBOSE(" - isStopped: " << inParticle.isStopped()); + ACTS_VERBOSE(" - endpoint: " << inParticle.getEndpoint().x << ", " + << inParticle.getEndpoint().y << ", " + << inParticle.getEndpoint().z); + const auto pid = particle.particleId(); + unordered.push_back(std::move(particle)); + edm4hepParticleMap[inParticle.getObjectID().index] = unordered.size() - 1; + processChildren(inParticle, pid, unordered, parentRelationship, + edm4hepParticleMap, nSecondaryVertices, maxGen); + } + ACTS_VERBOSE("Primary vertex complete, produced " + << (unordered.size() - startSize) << " particles and " + << nSecondaryVertices << " secondary vertices in " << maxGen + << " generations"); + setSubParticleIds(std::next(unordered.begin(), startSize), unordered.end()); + } + + ACTS_DEBUG("Found " << unordered.size() << " particles"); + + // @TODO: Order simhits by time + + SimParticleContainer particlesFinal; + SimParticleContainer particlesGenerator; + for (const auto& inParticle : mcParticleCollection) { + const std::size_t index = + edm4hepParticleMap.find(inParticle.getObjectID().index)->second; + const auto& particleInitial = unordered.at(index); + if (!inParticle.isCreatedInSimulation()) { + particlesGenerator.insert(particleInitial); + } + SimParticle particleFinal = particleInitial; + + float time = inParticle.getTime() * Acts::UnitConstants::ns; + for (const auto& daughter : inParticle.getDaughters()) { + if (!daughter.vertexIsNotEndpointOfParent()) { + time = daughter.getTime() * Acts::UnitConstants::ns; + break; + } + } + + particleFinal.setPosition4( + inParticle.getEndpoint()[0] * Acts::UnitConstants::mm, + inParticle.getEndpoint()[1] * Acts::UnitConstants::mm, + inParticle.getEndpoint()[2] * Acts::UnitConstants::mm, time); + + Acts::Vector3 momentumFinal = {inParticle.getMomentumAtEndpoint()[0], + inParticle.getMomentumAtEndpoint()[1], + inParticle.getMomentumAtEndpoint()[2]}; + particleFinal.setDirection(momentumFinal.normalized()); + particleFinal.setAbsoluteMomentum(momentumFinal.norm()); + + ACTS_VERBOSE("- Updated particle initial -> final, position: " + << particleInitial.fourPosition().transpose() << " -> " + << particleFinal.fourPosition().transpose()); + ACTS_VERBOSE(" momentum: " + << particleInitial.fourMomentum().transpose() << " -> " + << particleFinal.fourMomentum().transpose()); + + particlesFinal.insert(particleFinal); + } + + // Write ordered particles container to the EventStore + SimParticleContainer particlesInitial; + particlesInitial.insert(unordered.begin(), unordered.end()); + + if (!m_cfg.graphvizOutput.empty()) { + std::string path = perEventFilepath(m_cfg.graphvizOutput, "particles.dot", + ctx.eventNumber); + std::ofstream dot(path); + graphviz(dot, unordered, parentRelationship); + } + + SimHitContainer simHits; + + ACTS_DEBUG("Reading sim hits from " << m_cfg.inputSimHits.size() + << " sim hit collections"); + for (const auto& name : m_cfg.inputSimHits) { + const auto& inputHits = frame.get(name); + + for (const auto& hit : inputHits) { + auto simHit = EDM4hepUtil::readSimHit( + hit, + [&](const auto& inParticle) { + ACTS_VERBOSE("SimHit has source particle: " + << hit.getMCParticle().getObjectID().index); + auto it = edm4hepParticleMap.find(inParticle.getObjectID().index); + if (it == edm4hepParticleMap.end()) { + ACTS_ERROR( + "SimHit has source particle that we did not see before"); + return SimBarcode{}; + } + const auto& particle = unordered.at(it->second); + ACTS_VERBOSE("- " << inParticle.getObjectID().index << " -> " + << particle.particleId()); + return particle.particleId(); + }, + [&](std::uint64_t cellId) { + ACTS_VERBOSE("CellID: " << cellId); + + const auto& vm = m_cfg.dd4hepDetector->geometryService->detector() + .volumeManager(); + + const auto detElement = vm.lookupDetElement(cellId); + + ACTS_VERBOSE(" -> detElement: " << detElement.name()); + ACTS_VERBOSE(" -> id: " << detElement.id()); + ACTS_VERBOSE(" -> key: " << detElement.key()); + + Acts::Vector3 position; + position << detElement.nominal() + .worldTransformation() + .GetTranslation()[0], + detElement.nominal().worldTransformation().GetTranslation()[1], + detElement.nominal().worldTransformation().GetTranslation()[2]; + position *= Acts::UnitConstants::cm; + + ACTS_VERBOSE(" -> detElement position: " << position.transpose()); + + auto it = m_surfaceMap.find(detElement.key()); + if (it == m_surfaceMap.end()) { + ACTS_ERROR("Unable to find surface for detElement " + << detElement.name() << " with cellId " << cellId); + } + const auto* surface = it->second; + ACTS_VERBOSE(" -> surface: " << surface->geometryId()); + return surface->geometryId(); + }); + + simHits.insert(std::move(simHit)); + } + } + + if (m_cfg.sortSimHitsInTime) { + ACTS_DEBUG("Sorting sim hits in time"); + std::multimap hitsByParticle; + + for (std::size_t i = 0; i < simHits.size(); ++i) { + hitsByParticle.insert({simHits.nth(i)->particleId(), i}); + } + + for (auto it = hitsByParticle.begin(), end = hitsByParticle.end(); + it != end; it = hitsByParticle.upper_bound(it->first)) { + std::cout << "Particle " << it->first << " has " + << hitsByParticle.count(it->first) << " hits" << std::endl; + + std::vector hitIndices; + hitIndices.reserve(hitsByParticle.count(it->first)); + for (auto hitIndex : makeRange(hitsByParticle.equal_range(it->first))) { + hitIndices.push_back(hitIndex.second); + } + + if (logger().doPrint(Acts::Logging::VERBOSE)) { + ACTS_VERBOSE("Before sorting:"); + for (const auto& hitIdx : hitIndices) { + ACTS_VERBOSE(" - " << hitIdx << " / " << simHits.nth(hitIdx)->index() + << " " << simHits.nth(hitIdx)->time()); + } + } + + std::sort(hitIndices.begin(), hitIndices.end(), + [&](std::size_t a, std::size_t b) { + return simHits.nth(a)->time() < simHits.nth(b)->time(); + }); + + for (std::size_t i = 0; i < hitIndices.size(); ++i) { + auto& hit = *simHits.nth(hitIndices[i]); + SimHit updatedHit{hit.geometryId(), hit.particleId(), + hit.fourPosition(), hit.momentum4Before(), + hit.momentum4After(), int32_t(i)}; + hit = updatedHit; + } + + if (logger().doPrint(Acts::Logging::VERBOSE)) { + ACTS_VERBOSE("After sorting:"); + for (const auto& hitIdx : hitIndices) { + ACTS_VERBOSE(" - " << hitIdx << " / " << simHits.nth(hitIdx)->index() + << " " << simHits.nth(hitIdx)->time()); + } + } + } + } + + m_outputParticlesInitial(ctx, std::move(particlesInitial)); + m_outputParticlesFinal(ctx, std::move(particlesFinal)); + m_outputParticlesGenerator(ctx, std::move(particlesGenerator)); + + m_outputSimHits(ctx, std::move(simHits)); + + return ProcessCode::SUCCESS; +} + +void EDM4hepReader::processChildren( + const edm4hep::MCParticle& inParticle, SimBarcode parentId, + SimParticleContainer::sequence_type& particles, + ParentRelationship& parentRelationship, + std::unordered_map& particleMap, + std::size_t& nSecondaryVertices, std::size_t& maxGen) const { + constexpr auto indent = [&](std::size_t n) { + std::string result; + for (std::size_t i = 0; i < n; ++i) { + result += " "; + } + return result; + }; + + const std::size_t gen = parentId.generation(); + maxGen = std::max(maxGen, gen); + + ACTS_VERBOSE(indent(gen) << " - processing daughters for input particle " + << inParticle.id()); + ACTS_VERBOSE(indent(gen) << " -> found " << inParticle.daughters_size() + << " daughter(s)"); + + bool parentDecayed = + std::any_of(inParticle.daughters_begin(), inParticle.daughters_end(), + [](const edm4hep::MCParticle& daughter) { + return !daughter.vertexIsNotEndpointOfParent(); + }); + std::size_t secondaryVertex = 0; + if (parentDecayed) { + ACTS_VERBOSE(indent(gen) << " -> parent decays"); + secondaryVertex = ++nSecondaryVertices; + } + + std::size_t parentIndex = particles.size() - 1; + + std::size_t nParticles = 0; + for (const auto& daughter : inParticle.getDaughters()) { + SimParticle particle = EDM4hepUtil::readParticle(daughter); + + auto pid = parentId.makeDescendant(nParticles++); + if (daughter.vertexIsNotEndpointOfParent()) { + // incoming particle survived, interaction via descendant + } else { + // incoming particle decayed + pid = pid.setVertexSecondary(secondaryVertex); + } + particle.setParticleId(pid); + + ACTS_VERBOSE(indent(particle.particleId().generation()) + << "+ add particle " << particle); + ACTS_VERBOSE(indent(particle.particleId().generation()) + << " - generation: " << particle.particleId().generation()); + ACTS_VERBOSE(indent(particle.particleId().generation()) + << " - at " << particle.position().transpose()); + ACTS_VERBOSE(indent(particle.particleId().generation()) + << " - createdInSim: " + << daughter.isCreatedInSimulation()); + ACTS_VERBOSE(indent(particle.particleId().generation()) + << " - vertexIsNotEndpointOfParent: " + << daughter.vertexIsNotEndpointOfParent()); + ACTS_VERBOSE(indent(particle.particleId().generation()) + << " - isStopped: " << daughter.isStopped()); + ACTS_VERBOSE(indent(particle.particleId().generation()) + << " - endpoint: " << daughter.getEndpoint().x << ", " + << daughter.getEndpoint().y << ", " + << daughter.getEndpoint().z); + + particles.push_back(std::move(particle)); + particleMap[daughter.getObjectID().index] = particles.size() - 1; + parentRelationship[particles.size() - 1] = parentIndex; + processChildren(daughter, pid, particles, parentRelationship, particleMap, + nSecondaryVertices, maxGen); + } +} + +void EDM4hepReader::setSubParticleIds( + const SimParticleContainer::sequence_type::iterator& begin, + const SimParticleContainer::sequence_type::iterator& end) { + std::vector numByGeneration; + numByGeneration.reserve(10); + + for (auto it = begin; it != end; ++it) { + auto& particle = *it; + const auto pid = particle.particleId(); + if (pid.generation() >= numByGeneration.size()) { + numByGeneration.resize(pid.generation() + 1, 0); + } + unsigned int nextSubParticle = numByGeneration[pid.generation()]++; + + auto newPid = particle.particleId().setSubParticle(nextSubParticle); + particle.setParticleId(newPid); + } +} + +} // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/src/EDM4hepSimHitReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepSimHitReader.cpp deleted file mode 100644 index aa036811f3f..00000000000 --- a/Examples/Io/EDM4hep/src/EDM4hepSimHitReader.cpp +++ /dev/null @@ -1,105 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2022 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/EDM4hep/EDM4hepSimHitReader.hpp" - -#include "Acts/Definitions/Units.hpp" -#include "ActsExamples/EventData/SimHit.hpp" -#include "ActsExamples/EventData/SimParticle.hpp" -#include "ActsExamples/Framework/WhiteBoard.hpp" -#include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp" - -#include -#include -#include - -namespace ActsExamples { - -EDM4hepSimHitReader::EDM4hepSimHitReader( - const EDM4hepSimHitReader::Config& config, Acts::Logging::Level level) - : m_cfg(config), - m_logger(Acts::getDefaultLogger("EDM4hepSimHitReader", level)) { - m_reader.openFile(m_cfg.inputPath); - - m_outputParticles.maybeInitialize(m_cfg.outputParticles); - m_outputSimHits.initialize(m_cfg.outputSimHits); - - m_eventsRange = std::make_pair(0, m_reader.getEntries("events")); -} - -std::string EDM4hepSimHitReader::EDM4hepSimHitReader::name() const { - return "EDM4hepSimHitReader"; -} - -std::pair EDM4hepSimHitReader::availableEvents() - const { - return m_eventsRange; -} - -ProcessCode EDM4hepSimHitReader::read(const AlgorithmContext& ctx) { - podio::Frame frame = m_reader.readEntry("events", ctx.eventNumber); - - const auto& mcParticleCollection = - frame.get(m_cfg.inputParticles); - - if (!m_cfg.outputParticles.empty()) { - SimParticleContainer::sequence_type unordered; - - for (const auto& mcParticle : mcParticleCollection) { - auto particle = EDM4hepUtil::readParticle( - mcParticle, [](const edm4hep::MCParticle& p) { - ActsFatras::Barcode result; - result.setParticle(EDM4hepUtil::podioObjectIDToInteger(p.id())); - return result; - }); - unordered.push_back(particle); - } - - // Write ordered particles container to the EventStore - SimParticleContainer particles; - particles.insert(unordered.begin(), unordered.end()); - m_outputParticles(ctx, std::move(particles)); - } - - SimHitContainer::sequence_type unordered; - - const auto& simTrackerHitCollection = - frame.get(m_cfg.inputSimTrackerHits); - - for (const auto& simTrackerHit : simTrackerHitCollection) { - try { - auto hit = EDM4hepUtil::readSimHit( - simTrackerHit, - [](const edm4hep::MCParticle& particle) { - ActsFatras::Barcode result; - result.setParticle( - EDM4hepUtil::podioObjectIDToInteger(particle.id())); - return result; - }, - [&](std::uint64_t cellId) { - auto detElement = m_cfg.dd4hepDetector->geometryService->detector() - .volumeManager() - .lookupDetElement(cellId); - Acts::GeometryIdentifier result = detElement.volumeID(); - return result; - }); - unordered.push_back(std::move(hit)); - } catch (...) { - ACTS_ERROR("EDM4hepSimHitReader: failed to convert SimTrackerHit"); - continue; - } - } - - SimHitContainer simHits; - simHits.insert(unordered.begin(), unordered.end()); - m_outputSimHits(ctx, std::move(simHits)); - - return ProcessCode::SUCCESS; -} - -} // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/src/EDM4hepSimHitWriter.cpp b/Examples/Io/EDM4hep/src/EDM4hepSimHitWriter.cpp index 4421ce8e3a0..711a87d77ca 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepSimHitWriter.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepSimHitWriter.cpp @@ -89,6 +89,7 @@ ProcessCode EDM4hepSimHitWriter::writeT(const AlgorithmContext& ctx, frame.put(std::move(mcParticles), m_cfg.outputParticles); frame.put(std::move(simTrackerHitCollection), m_cfg.outputSimTrackerHits); + std::lock_guard lock{m_writeMutex}; m_writer.writeFrame(frame, "events"); return ProcessCode::SUCCESS; diff --git a/Examples/Io/EDM4hep/src/EDM4hepTrackReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepTrackReader.cpp index a634639e489..0db917a2714 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepTrackReader.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepTrackReader.cpp @@ -28,12 +28,12 @@ EDM4hepTrackReader::EDM4hepTrackReader(const Config& config, m_outputTracks.initialize(m_cfg.outputTracks); - m_reader.openFile(m_cfg.inputPath); + m_eventsRange = {0, reader().getEntries("events")}; } std::pair EDM4hepTrackReader::availableEvents() const { - return {0, m_reader.getEntries("events")}; + return m_eventsRange; } std::string EDM4hepTrackReader::EDM4hepTrackReader::name() const { @@ -41,7 +41,7 @@ std::string EDM4hepTrackReader::EDM4hepTrackReader::name() const { } ProcessCode EDM4hepTrackReader::read(const AlgorithmContext& ctx) { - podio::Frame frame = m_reader.readEntry("events", ctx.eventNumber); + podio::Frame frame = reader().readEntry("events", ctx.eventNumber); const auto& trackCollection = frame.get(m_cfg.inputTracks); @@ -66,4 +66,14 @@ ProcessCode EDM4hepTrackReader::read(const AlgorithmContext& ctx) { return ProcessCode::SUCCESS; } +podio::ROOTFrameReader& EDM4hepTrackReader::reader() { + bool exists = false; + auto& reader = m_reader.local(exists); + if (!exists) { + reader.openFile(m_cfg.inputPath); + } + + return reader; +} + } // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/src/EDM4hepTrackWriter.cpp b/Examples/Io/EDM4hep/src/EDM4hepTrackWriter.cpp index 3911bfa3c8f..c9d50032351 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepTrackWriter.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepTrackWriter.cpp @@ -48,6 +48,7 @@ ProcessCode EDM4hepTrackWriter::writeT(const AlgorithmContext& context, frame.put(std::move(trackCollection), m_cfg.outputTracks); + std::lock_guard guard(m_writeMutex); m_writer.writeFrame(frame, "events"); return ProcessCode::SUCCESS; diff --git a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp index aef0d6cb191..8ec14043753 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp @@ -8,6 +8,7 @@ #include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp" +#include "Acts/Definitions/Common.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/Charge.hpp" #include "Acts/EventData/MultiTrajectory.hpp" @@ -22,6 +23,8 @@ #include "edm4hep/TrackState.h" +using namespace Acts::UnitLiterals; + namespace ActsExamples { ActsFatras::Particle EDM4hepUtil::readParticle( @@ -42,13 +45,11 @@ ActsFatras::Particle EDM4hepUtil::readParticle( from.getTime() * Acts::UnitConstants::ns); // Only used for direction; normalization/units do not matter - to.setDirection(from.getMomentum()[0], from.getMomentum()[1], - from.getMomentum()[2]); + Acts::Vector3 momentum = {from.getMomentum()[0], from.getMomentum()[1], + from.getMomentum()[2]}; + to.setDirection(momentum.normalized()); - to.setAbsoluteMomentum(std::hypot(from.getMomentum()[0], - from.getMomentum()[1], - from.getMomentum()[2]) * - Acts::UnitConstants::GeV); + to.setAbsoluteMomentum(momentum.norm() * 1_GeV); return to; } @@ -70,38 +71,36 @@ ActsFatras::Hit EDM4hepUtil::readSimHit( const edm4hep::SimTrackerHit& from, const MapParticleIdFrom& particleMapper, const MapGeometryIdFrom& geometryMapper) { ActsFatras::Barcode particleId = particleMapper(from.getMCParticle()); - Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID()); - const auto mass = from.getMCParticle().getMass(); - const Acts::ActsVector<3> momentum{ - from.getMomentum().x * Acts::UnitConstants::GeV, - from.getMomentum().y * Acts::UnitConstants::GeV, - from.getMomentum().z * Acts::UnitConstants::GeV, + const auto mass = from.getMCParticle().getMass() * 1_GeV; + const Acts::Vector3 momentum{ + from.getMomentum().x * 1_GeV, + from.getMomentum().y * 1_GeV, + from.getMomentum().z * 1_GeV, }; const auto energy = std::hypot(momentum.norm(), mass); - ActsFatras::Hit::Vector4 pos4{ - from.getPosition().x * Acts::UnitConstants::mm, - from.getPosition().y * Acts::UnitConstants::mm, - from.getPosition().z * Acts::UnitConstants::mm, - from.getTime() * Acts::UnitConstants::ns, + Acts::Vector4 pos4{ + from.getPosition().x * 1_mm, + from.getPosition().y * 1_mm, + from.getPosition().z * 1_mm, + from.getTime() * 1_ns, }; - ActsFatras::Hit::Vector4 mom4{ + Acts::Vector4 mom4{ momentum.x(), momentum.y(), momentum.z(), energy, }; - // TODO no EDM4hep equivalent? - ActsFatras::Hit::Vector4 delta4{ - 0 * Acts::UnitConstants::GeV, 0 * Acts::UnitConstants::GeV, - 0 * Acts::UnitConstants::GeV, - 0 * Acts::UnitConstants::GeV, // sth.getEDep() - }; + Acts::Vector4 delta4 = Acts::Vector4::Zero(); + delta4[Acts::eEnergy] = -from.getEDep() * Acts::UnitConstants::GeV; + + Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID()); - // TODO no EDM4hep equivalent? + // Can extract from time, but we need a complete picture of the trajectory + // first int32_t index = -1; return ActsFatras::Hit(geometryId, particleId, pos4, mom4, mom4 + delta4, diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp index f1864b86b60..ae179a88b0c 100644 --- a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp +++ b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp @@ -195,8 +195,8 @@ int ActsExamples::VertexPerformanceWriter::getNumberOfReconstructableVertices( // traverse the array for frequency for (const auto& p : collection) { - int secVtxId = p.particleId().vertexSecondary(); - if (secVtxId != 0) { + int generation = p.particleId().generation(); + if (generation > 0) { // truthparticle from secondary vtx continue; } @@ -221,8 +221,8 @@ int ActsExamples::VertexPerformanceWriter::getNumberOfTruePriVertices( std::set allPriVtxIds; for (const auto& p : collection) { int priVtxId = p.particleId().vertexPrimary(); - int secVtxId = p.particleId().vertexSecondary(); - if (secVtxId != 0) { + int generation = p.particleId().generation(); + if (generation > 0) { // truthparticle from secondary vtx continue; } @@ -511,9 +511,9 @@ ActsExamples::ProcessCode ActsExamples::VertexPerformanceWriter::writeT( for (std::size_t j = 0; j < associatedTruthParticles.size(); ++j) { const auto& particle = associatedTruthParticles[j]; int priVtxId = particle.particleId().vertexPrimary(); - int secVtxId = particle.particleId().vertexSecondary(); + int generation = particle.particleId().generation(); - if (secVtxId != 0) { + if (generation > 0) { // truthparticle from secondary vtx continue; } diff --git a/Examples/Python/src/EDM4hep.cpp b/Examples/Python/src/EDM4hep.cpp index 4daedbadada..a20c171603c 100644 --- a/Examples/Python/src/EDM4hep.cpp +++ b/Examples/Python/src/EDM4hep.cpp @@ -7,12 +7,12 @@ // file, You can obtain one at http://mozilla.org/MPL/2.0/. #include "Acts/Plugins/Python/Utilities.hpp" +#include "ActsExamples/DD4hepDetector/DD4hepDetector.hpp" #include "ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp" #include "ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp" #include "ActsExamples/Io/EDM4hep/EDM4hepMultiTrajectoryWriter.hpp" -#include "ActsExamples/Io/EDM4hep/EDM4hepParticleReader.hpp" #include "ActsExamples/Io/EDM4hep/EDM4hepParticleWriter.hpp" -#include "ActsExamples/Io/EDM4hep/EDM4hepSimHitReader.hpp" +#include "ActsExamples/Io/EDM4hep/EDM4hepReader.hpp" #include "ActsExamples/Io/EDM4hep/EDM4hepSimHitWriter.hpp" #include "ActsExamples/Io/EDM4hep/EDM4hepTrackReader.hpp" #include "ActsExamples/Io/EDM4hep/EDM4hepTrackWriter.hpp" @@ -33,9 +33,11 @@ void addEDM4hep(Context& ctx) { auto mex = ctx.get("examples"); auto edm4hep = mex.def_submodule("_edm4hep"); - ACTS_PYTHON_DECLARE_READER(ActsExamples::EDM4hepSimHitReader, edm4hep, - "EDM4hepSimHitReader", inputPath, inputParticles, - outputSimHits, dd4hepDetector); + ACTS_PYTHON_DECLARE_READER( + ActsExamples::EDM4hepReader, edm4hep, "EDM4hepReader", inputPath, + inputParticles, inputSimHits, outputParticlesInitial, + outputParticlesFinal, outputParticlesGenerator, outputSimHits, + graphvizOutput, dd4hepDetector, trackingGeometry, sortSimHitsInTime); ACTS_PYTHON_DECLARE_WRITER(ActsExamples::EDM4hepSimHitWriter, edm4hep, "EDM4hepSimHitWriter", inputSimHits, @@ -51,10 +53,6 @@ void addEDM4hep(Context& ctx) { "EDM4hepMeasurementWriter", inputMeasurements, inputClusters, outputPath); - ACTS_PYTHON_DECLARE_READER(ActsExamples::EDM4hepParticleReader, edm4hep, - "EDM4hepParticleReader", inputPath, inputParticles, - outputParticles); - ACTS_PYTHON_DECLARE_WRITER(ActsExamples::EDM4hepParticleWriter, edm4hep, "EDM4hepParticleWriter", inputParticles, outputPath, outputParticles); diff --git a/Examples/Python/tests/test_reader.py b/Examples/Python/tests/test_reader.py index 8100ef5489f..5850e3cd051 100644 --- a/Examples/Python/tests/test_reader.py +++ b/Examples/Python/tests/test_reader.py @@ -346,6 +346,7 @@ def generate_input_test_edm4hep_simhit_reader(input, output): ddsim.compactFile = input ddsim.enableGun = True ddsim.gun.direction = (1, 0, 0) + ddsim.gun.particle = "pi-" ddsim.gun.distribution = "eta" ddsim.numberOfEvents = 10 ddsim.outputFile = output @@ -355,8 +356,8 @@ def generate_input_test_edm4hep_simhit_reader(input, output): @pytest.mark.slow @pytest.mark.edm4hep @pytest.mark.skipif(not edm4hepEnabled, reason="EDM4hep is not set up") -def test_edm4hep_simhit_reader(tmp_path): - from acts.examples.edm4hep import EDM4hepSimHitReader +def test_edm4hep_simhit_particle_reader(tmp_path): + from acts.examples.edm4hep import EDM4hepReader tmp_file = str(tmp_path / "output_edm4hep.root") odd_xml_file = str(getOpenDataDetectorDirectory() / "xml" / "OpenDataDetector.xml") @@ -373,17 +374,34 @@ def test_edm4hep_simhit_reader(tmp_path): s = Sequencer(numThreads=1) s.addReader( - EDM4hepSimHitReader( + EDM4hepReader( level=acts.logging.INFO, inputPath=tmp_file, + inputSimHits=[ + "PixelBarrelReadout", + "PixelEndcapReadout", + "ShortStripBarrelReadout", + "ShortStripEndcapReadout", + "LongStripBarrelReadout", + "LongStripEndcapReadout", + ], + outputParticlesGenerator="particles_input", + outputParticlesInitial="particles_initial", + outputParticlesFinal="particles_final", outputSimHits="simhits", dd4hepDetector=detector, + trackingGeometry=trackingGeometry, ) ) alg = AssertCollectionExistsAlg("simhits", "check_alg", acts.logging.WARNING) s.addAlgorithm(alg) + alg = AssertCollectionExistsAlg( + "particles_input", "check_alg", acts.logging.WARNING + ) + s.addAlgorithm(alg) + s.run() assert alg.events_seen == 10 @@ -437,55 +455,6 @@ def test_edm4hep_measurement_reader(tmp_path, fatras, conf_const): assert alg.events_seen == 10 -@pytest.mark.edm4hep -@pytest.mark.skipif(not edm4hepEnabled, reason="EDM4hep is not set up") -def test_edm4hep_particle_reader(tmp_path, conf_const, ptcl_gun): - from acts.examples.edm4hep import ( - EDM4hepParticleWriter, - EDM4hepParticleReader, - ) - - s = Sequencer(numThreads=1, events=10, logLevel=acts.logging.WARNING) - evGen = ptcl_gun(s) - - out = tmp_path / "particles_edm4hep.root" - - out.mkdir() - - s.addWriter( - conf_const( - EDM4hepParticleWriter, - acts.logging.WARNING, - inputParticles=evGen.config.outputParticles, - outputPath=str(out), - ) - ) - - s.run() - - # reset the seeder - s = Sequencer(numThreads=1, logLevel=acts.logging.WARNING) - - s.addReader( - conf_const( - EDM4hepParticleReader, - acts.logging.WARNING, - inputPath=str(out), - outputParticles="input_particles", - ) - ) - - alg = AssertCollectionExistsAlg( - "input_particles", "check_alg", acts.logging.WARNING - ) - - s.addAlgorithm(alg) - - s.run() - - assert alg.events_seen == 10 - - @pytest.mark.edm4hep @pytest.mark.skipif(not edm4hepEnabled, reason="EDM4hep is not set up") def test_edm4hep_tracks_reader(tmp_path): diff --git a/Examples/Scripts/Python/full_chain_odd.py b/Examples/Scripts/Python/full_chain_odd.py index 41249f585cf..9a8040b37b7 100755 --- a/Examples/Scripts/Python/full_chain_odd.py +++ b/Examples/Scripts/Python/full_chain_odd.py @@ -11,6 +11,7 @@ addGeant4, ParticleSelectorConfig, addDigitization, + addParticleSelection, ) from acts.examples.reconstruction import ( addSeeding, @@ -28,10 +29,14 @@ ) from common import getOpenDataDetectorDirectory from acts.examples.odd import getOpenDataDetector +import acts.examples.edm4hep parser = argparse.ArgumentParser(description="Full chain with the OpenDataDetector") parser.add_argument("--events", "-n", help="Number of events", type=int, default=100) + +parser.add_argument("--skip", "-s", help="Number of events", type=int, default=0) +parser.add_argument("--edm4hep", help="Use edm4hep inputs", type=pathlib.Path) parser.add_argument( "--geant4", help="Use Geant4 instead of fatras", action="store_true" ) @@ -75,81 +80,126 @@ s = acts.examples.Sequencer( events=args["events"], - numThreads=1, + skip=args["skip"], + numThreads=1 if g4_simulation else -1, outputDir=str(outputDir), ) -if not ttbar: - addParticleGun( - s, - MomentumConfig(1.0 * u.GeV, 10.0 * u.GeV, transverse=True), - EtaConfig(-3.0, 3.0), - PhiConfig(0.0, 360.0 * u.degree), - ParticleConfig(4, acts.PdgParticle.eMuon, randomizeCharge=True), - vtxGen=acts.examples.GaussianVertexGenerator( - mean=acts.Vector4(0, 0, 0, 0), - stddev=acts.Vector4(0.0125 * u.mm, 0.0125 * u.mm, 55.5 * u.mm, 1.0 * u.ns), - ), - multiplicity=200, - rnd=rnd, +if args["edm4hep"]: + edm4hepReader = acts.examples.edm4hep.EDM4hepReader( + inputPath=str(args["edm4hep"]), + inputSimHits=[ + "PixelBarrelReadout", + "PixelEndcapReadout", + "ShortStripBarrelReadout", + "ShortStripEndcapReadout", + "LongStripBarrelReadout", + "LongStripEndcapReadout", + ], + outputParticlesGenerator="particles_input", + outputParticlesInitial="particles_initial", + outputParticlesFinal="particles_final", + outputSimHits="simhits", + graphvizOutput="graphviz", + dd4hepDetector=detector, + trackingGeometry=trackingGeometry, + sortSimHitsInTime=True, + level=acts.logging.INFO, ) -else: - addPythia8( - s, - hardProcess=["Top:qqbar2ttbar=on"], - npileup=50, - vtxGen=acts.examples.GaussianVertexGenerator( - mean=acts.Vector4(0, 0, 0, 0), - stddev=acts.Vector4(0.0125 * u.mm, 0.0125 * u.mm, 55.5 * u.mm, 5.0 * u.ns), - ), - rnd=rnd, - outputDirRoot=outputDir, - # outputDirCsv=outputDir, - ) -if g4_simulation: - if s.config.numThreads != 1: - raise ValueError("Geant 4 simulation does not support multi-threading") - - # Pythia can sometime simulate particles outside the world volume, a cut on the Z of the track help mitigate this effect - # Older version of G4 might not work, this as has been tested on version `geant4-11-00-patch-03` - # For more detail see issue #1578 - addGeant4( + s.addReader(edm4hepReader) + s.addWhiteboardAlias("particles", edm4hepReader.config.outputParticlesGenerator) + + addParticleSelection( s, - detector, - trackingGeometry, - field, - preSelectParticles=ParticleSelectorConfig( + config=ParticleSelectorConfig( rho=(0.0, 24 * u.mm), absZ=(0.0, 1.0 * u.m), eta=(-3.0, 3.0), pt=(150 * u.MeV, None), removeNeutral=True, ), - outputDirRoot=outputDir, - # outputDirCsv=outputDir, - rnd=rnd, - killVolume=trackingGeometry.worldVolume, - killAfterTime=25 * u.ns, + inputParticles="particles", + outputParticles="particles_selected", ) + + else: - addFatras( - s, - trackingGeometry, - field, - preSelectParticles=ParticleSelectorConfig( - rho=(0.0, 24 * u.mm), - absZ=(0.0, 1.0 * u.m), - eta=(-3.0, 3.0), - pt=(150 * u.MeV, None), - removeNeutral=True, + if not ttbar: + addParticleGun( + s, + MomentumConfig(1.0 * u.GeV, 10.0 * u.GeV, transverse=True), + EtaConfig(-3.0, 3.0), + PhiConfig(0.0, 360.0 * u.degree), + ParticleConfig(4, acts.PdgParticle.eMuon, randomizeCharge=True), + vtxGen=acts.examples.GaussianVertexGenerator( + mean=acts.Vector4(0, 0, 0, 0), + stddev=acts.Vector4( + 0.0125 * u.mm, 0.0125 * u.mm, 55.5 * u.mm, 1.0 * u.ns + ), + ), + multiplicity=200, + rnd=rnd, + ) + else: + addPythia8( + s, + hardProcess=["Top:qqbar2ttbar=on"], + npileup=50, + vtxGen=acts.examples.GaussianVertexGenerator( + mean=acts.Vector4(0, 0, 0, 0), + stddev=acts.Vector4( + 0.0125 * u.mm, 0.0125 * u.mm, 55.5 * u.mm, 5.0 * u.ns + ), + ), + rnd=rnd, + outputDirRoot=outputDir, + # outputDirCsv=outputDir, + ) + + if g4_simulation: + if s.config.numThreads != 1: + raise ValueError("Geant 4 simulation does not support multi-threading") + + # Pythia can sometime simulate particles outside the world volume, a cut on the Z of the track help mitigate this effect + # Older version of G4 might not work, this as has been tested on version `geant4-11-00-patch-03` + # For more detail see issue #1578 + addGeant4( + s, + detector, + trackingGeometry, + field, + preSelectParticles=ParticleSelectorConfig( + rho=(0.0, 24 * u.mm), + absZ=(0.0, 1.0 * u.m), + eta=(-3.0, 3.0), + pt=(150 * u.MeV, None), + removeNeutral=True, + ), + outputDirRoot=outputDir, + # outputDirCsv=outputDir, + rnd=rnd, + killVolume=trackingGeometry.worldVolume, + killAfterTime=25 * u.ns, + ) + else: + addFatras( + s, + trackingGeometry, + field, + preSelectParticles=ParticleSelectorConfig( + rho=(0.0, 24 * u.mm), + absZ=(0.0, 1.0 * u.m), + eta=(-3.0, 3.0), + pt=(150 * u.MeV, None), + removeNeutral=True, + ) + if ttbar + else ParticleSelectorConfig(), + enableInteractions=True, + outputDirRoot=outputDir, + # outputDirCsv=outputDir, + rnd=rnd, ) - if ttbar - else ParticleSelectorConfig(), - enableInteractions=True, - outputDirRoot=outputDir, - # outputDirCsv=outputDir, - rnd=rnd, - ) addDigitization( s, @@ -168,15 +218,6 @@ TruthSeedRanges(pt=(1.0 * u.GeV, None), eta=(-3.0, 3.0), nHits=(9, None)) if ttbar else TruthSeedRanges(), - initialSigmas=[ - 1 * u.mm, - 1 * u.mm, - 1 * u.degree, - 1 * u.degree, - 0.1 / u.GeV, - 1 * u.ns, - ], - initialVarInflation=[1.0] * 6, geoSelectionConfigFile=oddSeedingSel, outputDirRoot=outputDir, # outputDirCsv=outputDir, @@ -223,9 +264,7 @@ addAmbiguityResolution( s, AmbiguityResolutionConfig( - maximumSharedHits=3, - maximumIterations=1000000, - nMeasurementsMin=7, + maximumSharedHits=3, maximumIterations=1000000, nMeasurementsMin=7 ), outputDirRoot=outputDir, writeCovMat=True, diff --git a/Fatras/include/ActsFatras/EventData/Barcode.hpp b/Fatras/include/ActsFatras/EventData/Barcode.hpp index 84e2dd0ba78..35df07d91f1 100644 --- a/Fatras/include/ActsFatras/EventData/Barcode.hpp +++ b/Fatras/include/ActsFatras/EventData/Barcode.hpp @@ -99,7 +99,7 @@ class Barcode : public Acts::MultiIndex { using Base::Value; // Construct an invalid barcode with all levels set to zero. - Barcode() : Base(Base::Zeros()) {} + constexpr Barcode() : Base(Base::Zeros()) {} Barcode(const Barcode&) = default; Barcode(Barcode&&) = default; Barcode& operator=(const Barcode&) = default; @@ -117,15 +117,30 @@ class Barcode : public Acts::MultiIndex { constexpr Value subParticle() const { return level(4); } /// Set the primary vertex identifier. - constexpr Barcode& setVertexPrimary(Value id) { return set(0, id), *this; } + constexpr Barcode& setVertexPrimary(Value id) { + set(0, id); + return *this; + } /// Set the secondary vertex identifier. - constexpr Barcode& setVertexSecondary(Value id) { return set(1, id), *this; } + constexpr Barcode& setVertexSecondary(Value id) { + set(1, id); + return *this; + } /// Set the parent particle identifier. - constexpr Barcode& setParticle(Value id) { return set(2, id), *this; } + constexpr Barcode& setParticle(Value id) { + set(2, id); + return *this; + } /// Set the particle identifier. - constexpr Barcode& setGeneration(Value id) { return set(3, id), *this; } + constexpr Barcode& setGeneration(Value id) { + set(3, id); + return *this; + } /// Set the process identifier. - constexpr Barcode& setSubParticle(Value id) { return set(4, id), *this; } + constexpr Barcode& setSubParticle(Value id) { + set(4, id); + return *this; + } /// Construct a new barcode representing a descendant particle. ///