feat: Use truth vertex EDM in Examples (acts-project#2998)

Currently truth vertices are not well defined. This is improved here by modifying the existing but unused truth vertex EDM and wiring it to the event generators and relevant performance writers.

blocked by
- acts-project#2904

CI/physmon/workflows/physmon_simulation.py

@@ -58,6 +58,7 @@
                 ]
             ],
             outputParticles="particles_input",
+            outputVertices="vertices_input",
             randomNumbers=rnd,
         )
     )

Core/include/Acts/Utilities/MultiIndex.hpp

@@ -142,6 +142,9 @@ class MultiIndex {
   friend constexpr bool operator==(MultiIndex lhs, MultiIndex rhs) {
     return lhs.m_value == rhs.m_value;
   }
+  friend constexpr bool operator!=(MultiIndex lhs, MultiIndex rhs) {
+    return lhs.m_value != rhs.m_value;
+  }
   friend inline std::ostream& operator<<(std::ostream& os, MultiIndex idx) {
     // one level is always defined
     os << idx.level(0u);

Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.cpp

@@ -1,13 +1,14 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
+// Copyright (C) 2019-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/Generators/EventGenerator.hpp"
 
+#include "ActsExamples/EventData/SimVertex.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
 #include "ActsFatras/EventData/Barcode.hpp"
 #include "ActsFatras/EventData/Particle.hpp"
@@ -22,6 +23,9 @@ ActsExamples::EventGenerator::EventGenerator(const Config& cfg,
   if (m_cfg.outputParticles.empty()) {
     throw std::invalid_argument("Missing output particles collection");
   }
+  if (m_cfg.outputVertices.empty()) {
+    throw std::invalid_argument("Missing output vertices collection");
+  }
   if (m_cfg.generators.empty()) {
     throw std::invalid_argument("No generators are configured");
   }
@@ -30,6 +34,7 @@ ActsExamples::EventGenerator::EventGenerator(const Config& cfg,
   }
 
   m_outputParticles.initialize(m_cfg.outputParticles);
+  m_outputVertices.initialize(m_cfg.outputVertices);
 }
 
 std::string ActsExamples::EventGenerator::name() const {
@@ -44,6 +49,7 @@ ActsExamples::EventGenerator::availableEvents() const {
 ActsExamples::ProcessCode ActsExamples::EventGenerator::read(
     const AlgorithmContext& ctx) {
   SimParticleContainer particles;
+  SimVertexContainer vertices;
 
   auto rng = m_cfg.randomNumbers->spawnGenerator(ctx);
 
@@ -59,32 +65,49 @@ ActsExamples::ProcessCode ActsExamples::EventGenerator::read(
       // generate primary vertex position
       auto vertexPosition = (*generate.vertex)(rng);
       // generate particles associated to this vertex
-      auto vertexParticles = (*generate.particles)(rng);
+      auto [newVertices, newParticles] = (*generate.particles)(rng);
 
       ACTS_VERBOSE("Generate vertex at " << vertexPosition.transpose());
 
-      auto updateParticleInPlace = [&](ActsFatras::Particle& particle) {
+      auto updateParticleInPlace = [&](SimParticle& particle) {
         // only set the primary vertex, leave everything else as-is
         // using the number of primary vertices as the index ensures
         // that barcode=0 is not used, since it is used elsewhere
         // to signify elements w/o an associated particle.
-        const auto pid = ActsFatras::Barcode(particle.particleId())
+        const auto pid = SimBarcode(particle.particleId())
                              .setVertexPrimary(nPrimaryVertices);
         // move particle to the vertex
         const auto pos4 = (vertexPosition + particle.fourPosition()).eval();
         ACTS_VERBOSE(" - particle at " << pos4.transpose());
         // `withParticleId` returns a copy because it changes the identity
         particle = particle.withParticleId(pid).setPosition4(pos4);
       };
-      for (auto& vertexParticle : vertexParticles) {
+      for (auto& vertexParticle : newParticles) {
         updateParticleInPlace(vertexParticle);
       }
 
+      auto updateVertexInPlace = [&](SimVertex& vertex) {
+        // only set the primary vertex, leave everything else as-is
+        // using the number of primary vertices as the index ensures
+        // that barcode=0 is not used, since it is used elsewhere
+        // to signify elements w/o an associated particle.
+        vertex.id = SimVertexBarcode(vertex.vertexId())
+                        .setVertexPrimary(nPrimaryVertices);
+        // move vertex
+        const auto pos4 = (vertexPosition + vertex.position4).eval();
+        ACTS_VERBOSE(" - vertex at " << pos4.transpose());
+        vertex.position4 = pos4;
+      };
+      for (auto& vertex : newVertices) {
+        updateVertexInPlace(vertex);
+      }
+
       ACTS_VERBOSE("event=" << ctx.eventNumber << " generator=" << iGenerate
                             << " primary_vertex=" << nPrimaryVertices
-                            << " n_particles=" << vertexParticles.size());
+                            << " n_particles=" << newParticles.size());
 
-      particles.merge(std::move(vertexParticles));
+      particles.merge(std::move(newParticles));
+      vertices.merge(std::move(newVertices));
     }
   }
 
@@ -94,5 +117,6 @@ ActsExamples::ProcessCode ActsExamples::EventGenerator::read(
 
   // move generated event to the store
   m_outputParticles(ctx, std::move(particles));
+  m_outputVertices(ctx, std::move(vertices));
   return ProcessCode::SUCCESS;
 }

Examples/Algorithms/Generators/ActsExamples/Generators/EventGenerator.hpp

@@ -15,6 +15,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "ActsExamples/EventData/SimParticle.hpp"
+#include "ActsExamples/EventData/SimVertex.hpp"
 #include "ActsExamples/Framework/DataHandle.hpp"
 #include "ActsExamples/Framework/IReader.hpp"
 #include "ActsExamples/Framework/ProcessCode.hpp"
@@ -62,39 +63,42 @@ class EventGenerator final : public ActsExamples::IReader {
   };
 
   /// @brief Generator interface for a vertex position
-  struct VertexGenerator {
+  struct PrimaryVertexPositionGenerator {
     /// @brief Virtual destructor required
-    virtual ~VertexGenerator() = default;
+    virtual ~PrimaryVertexPositionGenerator() = default;
     /// @brief Generate a vertex position
     ///
     /// @param rng Shared random number generator instance
     /// @return Acts::Vector4 The vertex position
     virtual Acts::Vector4 operator()(RandomEngine& rng) const = 0;
   };
 
-  /// @brief Generator interface particles for a vertex
+  /// @brief Generator interface for vertices and particles
   struct ParticlesGenerator {
     /// @brief Virtual destructor required
     virtual ~ParticlesGenerator() = default;
-    /// @brief Generate particles for a vertex
+    /// @brief Generate vertices and particles for one interaction
     /// @note This method cannot be `const` because the Pythia8 generator
     ///       uses the Pythia8 interfaces, which is non-const
     ///
     /// @param rng Shared random number generator instance
-    /// @return SimParticleContainer The populated particle container for the vertex
-    virtual SimParticleContainer operator()(RandomEngine& rng) = 0;
+    /// @return The vertex and particle containers
+    virtual std::pair<SimVertexContainer, SimParticleContainer> operator()(
+        RandomEngine& rng) = 0;
   };
 
   /// @brief Combined struct which contains all generator components
   struct Generator {
     std::shared_ptr<MultiplicityGenerator> multiplicity = nullptr;
-    std::shared_ptr<VertexGenerator> vertex = nullptr;
+    std::shared_ptr<PrimaryVertexPositionGenerator> vertex = nullptr;
     std::shared_ptr<ParticlesGenerator> particles = nullptr;
   };
 
   struct Config {
     /// Name of the output particles collection.
     std::string outputParticles;
+    /// Name of the vertex collection.
+    std::string outputVertices;
     /// List of generators that should be used to generate the event.
     std::vector<Generator> generators;
     /// The random number service.
@@ -121,6 +125,7 @@ class EventGenerator final : public ActsExamples::IReader {
 
   WriteDataHandle<SimParticleContainer> m_outputParticles{this,
                                                           "OutputParticles"};
+  WriteDataHandle<SimVertexContainer> m_outputVertices{this, "OutputVertices"};
 };
 
 }  // namespace ActsExamples

Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.cpp

@@ -11,6 +11,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Common.hpp"
 #include "Acts/Definitions/ParticleData.hpp"
+#include "ActsExamples/EventData/SimHit.hpp"
 #include "ActsFatras/EventData/Barcode.hpp"
 #include "ActsFatras/EventData/Particle.hpp"
 
@@ -19,8 +20,9 @@
 #include <random>
 #include <utility>
 
-ActsExamples::ParametricParticleGenerator::ParametricParticleGenerator(
-    const Config& cfg)
+namespace ActsExamples {
+
+ParametricParticleGenerator::ParametricParticleGenerator(const Config& cfg)
     : m_cfg(cfg),
       m_charge(cfg.charge.value_or(Acts::findCharge(m_cfg.pdg).value_or(0))),
       m_mass(cfg.mass.value_or(Acts::findMass(m_cfg.pdg).value_or(0))),
@@ -36,8 +38,8 @@ ActsExamples::ParametricParticleGenerator::ParametricParticleGenerator(
       m_etaMin(-std::log(std::tan(0.5 * m_cfg.thetaMin))),
       m_etaMax(-std::log(std::tan(0.5 * m_cfg.thetaMax))) {}
 
-ActsExamples::SimParticleContainer
-ActsExamples::ParametricParticleGenerator::operator()(RandomEngine& rng) {
+std::pair<SimVertexContainer, SimParticleContainer>
+ParametricParticleGenerator::operator()(RandomEngine& rng) {
   using UniformIndex = std::uniform_int_distribution<std::uint8_t>;
   using UniformReal = std::uniform_real_distribution<double>;
 
@@ -58,13 +60,18 @@ ActsExamples::ParametricParticleGenerator::operator()(RandomEngine& rng) {
   UniformReal etaDist(m_etaMin, m_etaMax);
   UniformReal pDist(m_cfg.pMin, m_cfg.pMax);
 
-  SimParticleContainer particles;
-  particles.reserve(m_cfg.numParticles);
+  SimVertexContainer::sequence_type vertices;
+  SimParticleContainer::sequence_type particles;
+
+  // create the primary vertex
+  auto& primaryVertex =
+      vertices.emplace_back(0, SimVertex::Vector4(0., 0., 0., 0.));
 
   // counter will be reused as barcode particle number which must be non-zero.
   for (std::size_t ip = 1; ip <= m_cfg.numParticles; ++ip) {
     // all particles are treated as originating from the same primary vertex
-    const auto pid = ActsFatras::Barcode(0u).setParticle(ip);
+    const auto pid = SimBarcode(0u).setParticle(ip);
+    primaryVertex.outgoing.insert(pid);
 
     // draw parameters
     const unsigned int type = particleTypeChoice(rng);
@@ -100,5 +107,10 @@ ActsExamples::ParametricParticleGenerator::operator()(RandomEngine& rng) {
     particles.insert(particles.end(), std::move(particle));
   }
 
-  return particles;
+  std::pair<SimVertexContainer, SimParticleContainer> out;
+  out.first.insert(vertices.begin(), vertices.end());
+  out.second.insert(particles.begin(), particles.end());
+  return out;
 }
+
+}  // namespace ActsExamples

Examples/Algorithms/Generators/ActsExamples/Generators/ParametricParticleGenerator.hpp

@@ -67,7 +67,8 @@ class ParametricParticleGenerator : public EventGenerator::ParticlesGenerator {
   ParametricParticleGenerator(const Config& cfg);
 
   /// Generate a single primary vertex with the given number of particles.
-  SimParticleContainer operator()(RandomEngine& rng) override;
+  std::pair<SimVertexContainer, SimParticleContainer> operator()(
+      RandomEngine& rng) override;
 
  private:
   Config m_cfg;

Examples/Algorithms/Generators/ActsExamples/Generators/VertexGenerators.hpp

@@ -20,7 +20,8 @@
 
 namespace ActsExamples {
 
-struct FixedVertexGenerator : public EventGenerator::VertexGenerator {
+struct FixedPrimaryVertexPositionGenerator
+    : public EventGenerator::PrimaryVertexPositionGenerator {
   /// The fixed vertex position and time.
   Acts::Vector4 fixed = Acts::Vector4::Zero();
 
@@ -29,7 +30,8 @@ struct FixedVertexGenerator : public EventGenerator::VertexGenerator {
   }
 };
 
-struct GaussianVertexGenerator : public EventGenerator::VertexGenerator {
+struct GaussianPrimaryVertexPositionGenerator
+    : public EventGenerator::PrimaryVertexPositionGenerator {
   // standard deviation comes first, since it is more likely to be modified
   /// Vertex position and time standard deviation.
   Acts::Vector4 stddev = {0.0, 0.0, 0.0, 0.0};