feat: Add uber EDM4hep reader (acts-project#2939)

This unifies the EDM4hep particle and simhit reader with correct
particle barcodes, and working hit to particle association. I've tested
to work ~ ok on ODD with ddsim inputs. ddsim produces weird hits that I
haven't debugged yet, but given these inputs, the reco part seems to
work somewhat ok.


![hits](https://github.com/acts-project/acts/assets/1058585/6f337824-151b-4bb6-bb66-0e30201a64c3)

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<Acts::Vector2, Acts::Vector4, Acts::Vector3> 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 {

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(

Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp

@@ -18,6 +18,7 @@
 #include <string>
 
 #include <podio/ROOTFrameReader.h>
+#include <tbb/enumerable_thread_specific.h>
 
 namespace ActsExamples {
 
@@ -68,7 +69,9 @@ class EDM4hepMeasurementReader final : public IReader {
   std::pair<std::size_t, std::size_t> m_eventsRange;
   std::unique_ptr<const Acts::Logger> m_logger;
 
-  podio::ROOTFrameReader m_reader;
+  tbb::enumerable_thread_specific<podio::ROOTFrameReader> m_reader;
+
+  podio::ROOTFrameReader& reader();
 
   WriteDataHandle<MeasurementContainer> m_outputMeasurements{
       this, "OutputMeasurements"};

Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementWriter.hpp

@@ -66,6 +66,8 @@ class EDM4hepMeasurementWriter final : public WriterT<MeasurementContainer> {
 
   podio::ROOTFrameWriter m_writer;
 
+  std::mutex m_writeMutex;
+
   ReadDataHandle<ClusterContainer> m_inputClusters{this, "InputClusters"};
 };
 

Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMultiTrajectoryWriter.hpp

@@ -68,6 +68,7 @@ class EDM4hepMultiTrajectoryWriter : public WriterT<TrajectoriesContainer> {
  private:
   Config m_cfg;
 
+  std::mutex m_writeMutex;
   podio::ROOTFrameWriter m_writer;
 
   ReadDataHandle<IndexMultimap<ActsFatras::Barcode>>

Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepParticleWriter.hpp

@@ -11,6 +11,7 @@
 #include "ActsExamples/EventData/SimParticle.hpp"
 #include "ActsExamples/Framework/WriterT.hpp"
 
+#include <mutex>
 #include <string>
 
 #include <podio/ROOTFrameWriter.h>
@@ -55,6 +56,8 @@ class EDM4hepParticleWriter final : public WriterT<SimParticleContainer> {
  private:
   Config m_cfg;
 
+  std::mutex m_writeMutex;
+
   podio::ROOTFrameWriter m_writer;
 };
 

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 <memory>
+#include <string>
+
+#include <DD4hep/DetElement.h>
+#include <edm4hep/MCParticleCollection.h>
+#include <podio/ROOTFrameReader.h>
+#include <tbb/enumerable_thread_specific.h>
+
+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<std::string> 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<DD4hep::DD4hepDetector> dd4hepDetector;
+
+    /// Tracking geometry for cellID resolution.
+    std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry;
+
+    /// Whether to sort sim hits in time to produce index sequence
+    bool sortSimHitsInTime = false;
+  };
+
+  using ParentRelationship = std::unordered_map<std::size_t, std::size_t>;
+
+  /// 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<std::size_t, std::size_t> 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<int, std::size_t>& 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<std::size_t, std::size_t> m_eventsRange;
+  std::unique_ptr<const Acts::Logger> m_logger;
+
+  std::unordered_map<unsigned int, const Acts::Surface*> m_surfaceMap;
+
+  tbb::enumerable_thread_specific<podio::ROOTFrameReader> m_reader;
+
+  podio::ROOTFrameReader& reader();
+
+  WriteDataHandle<SimParticleContainer> m_outputParticlesInitial{
+      this, "OutputParticlesInitial"};
+  WriteDataHandle<SimParticleContainer> m_outputParticlesFinal{
+      this, "OutputParticlesFinal"};
+  WriteDataHandle<SimParticleContainer> m_outputParticlesGenerator{
+      this, "OutputParticlesGenerator"};
+
+  WriteDataHandle<SimHitContainer> m_outputSimHits{this, "OutputSimHits"};
+
+  void graphviz(std::ostream& os,
+                const SimParticleContainer::sequence_type& particles,
+                const ParentRelationship& parents) const;
+};
+
+}  // namespace ActsExamples

Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepSimHitWriter.hpp

@@ -67,6 +67,8 @@ class EDM4hepSimHitWriter final : public WriterT<SimHitContainer> {
 
   podio::ROOTFrameWriter m_writer;
 
+  std::mutex m_writeMutex;
+
   ReadDataHandle<SimParticleContainer> m_inputParticles{this, "InputParticles"};
 };