feat: add Truthiness algorithm for event quality assessment

Author: wdconinc

src/algorithms/reco/Truthiness.cc

@@ -0,0 +1,105 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Wouter Deconinck
+
+#include "Truthiness.h"
+
+#include <edm4eic/MCRecoParticleAssociationCollection.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+#include <edm4hep/MCParticleCollection.h>
+#include <edm4hep/utils/vector_utils.h>
+#include <cmath>
+#include <map>
+#include <set>
+
+namespace eicrecon {
+
+void Truthiness::process(const Truthiness::Input& input,
+                         const Truthiness::Output& /* output */) const {
+
+  const auto [mc_particles, rc_particles, associations] = input;
+
+  double truthiness = 0.0;
+
+  // Track which MC particles and reconstructed particles are covered by associations
+  std::set<edm4hep::MCParticle> associated_mc_particles;
+  std::set<edm4eic::ReconstructedParticle> associated_rc_particles;
+
+  // Process all associations
+  for (const auto& assoc : *associations) {
+    const auto& mc_part = assoc.getSim();
+    const auto& rc_part = assoc.getRec();
+
+    // Track that these particles have associations
+    associated_mc_particles.insert(mc_part);
+    associated_rc_particles.insert(rc_part);
+
+    // Get particle properties
+    const auto mc_momentum = mc_part.getMomentum();
+    const auto rc_momentum = rc_part.getMomentum();
+
+    const double mc_p      = edm4hep::utils::magnitude(mc_momentum);
+    const double mc_energy = std::sqrt(mc_p * mc_p + mc_part.getMass() * mc_part.getMass());
+    const double rc_energy = rc_part.getEnergy();
+
+    // Calculate energy difference squared
+    const double energy_diff    = mc_energy - rc_energy;
+    const double energy_penalty = energy_diff * energy_diff;
+
+    // Calculate momentum component differences squared
+    const double px_diff          = mc_momentum.x - rc_momentum.x;
+    const double py_diff          = mc_momentum.y - rc_momentum.y;
+    const double pz_diff          = mc_momentum.z - rc_momentum.z;
+    const double momentum_penalty = px_diff * px_diff + py_diff * py_diff + pz_diff * pz_diff;
+
+    // PDG code mismatch penalty
+    const double pdg_penalty = (mc_part.getPDG() != rc_part.getPDG()) ? 1.0 : 0.0;
+
+    const double assoc_penalty = energy_penalty + momentum_penalty + pdg_penalty;
+
+    debug("Association: MC PDG={} (E={:.3f}, p=[{:.3f},{:.3f},{:.3f}]) <-> "
+          "RC PDG={} (E={:.3f}, p=[{:.3f},{:.3f},{:.3f}])",
+          mc_part.getPDG(), mc_energy, mc_momentum.x, mc_momentum.y, mc_momentum.z,
+          rc_part.getPDG(), rc_energy, rc_momentum.x, rc_momentum.y, rc_momentum.z);
+    debug("  Energy penalty: {:.6f}, Momentum penalty: {:.6f}, PDG penalty: {:.0f}", energy_penalty,
+          momentum_penalty, pdg_penalty);
+
+    truthiness += assoc_penalty;
+  }
+
+  // Penalty for unassociated charged MC particles with generator status 2
+  int unassociated_mc_count = 0;
+  for (const auto& mc_part : *mc_particles) {
+    if (mc_part.getGeneratorStatus() == 2 && mc_part.getCharge() != 0.0) {
+      if (associated_mc_particles.find(mc_part) == associated_mc_particles.end()) {
+        unassociated_mc_count++;
+        debug("Unassociated MC particle: PDG={}, charge={:.1f}, status={}", mc_part.getPDG(),
+              mc_part.getCharge(), mc_part.getGeneratorStatus());
+      }
+    }
+  }
+  const double mc_penalty = static_cast<double>(unassociated_mc_count);
+  debug("Unassociated charged MC particles (status 2): {} (penalty: {:.0f})", unassociated_mc_count,
+        mc_penalty);
+  truthiness += mc_penalty;
+
+  // Penalty for unassociated reconstructed particles
+  int unassociated_rc_count = 0;
+  for (const auto& rc_part : *rc_particles) {
+    if (associated_rc_particles.find(rc_part) == associated_rc_particles.end()) {
+      unassociated_rc_count++;
+      debug("Unassociated reconstructed particle: PDG={}, E={:.3f}, p=[{:.3f},{:.3f},{:.3f}]",
+            rc_part.getPDG(), rc_part.getEnergy(), rc_part.getMomentum().x, rc_part.getMomentum().y,
+            rc_part.getMomentum().z);
+    }
+  }
+  const double rc_penalty = static_cast<double>(unassociated_rc_count);
+  debug("Unassociated reconstructed particles: {} (penalty: {:.0f})", unassociated_rc_count,
+        rc_penalty);
+  truthiness += rc_penalty;
+
+  // Report final truthiness
+  info("Event truthiness: {:.6f} (from {} associations, {} unassociated MC, {} unassociated RC)",
+       truthiness, associations->size(), unassociated_mc_count, unassociated_rc_count);
+}
+
+} // namespace eicrecon

src/algorithms/reco/Truthiness.h

@@ -0,0 +1,39 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Wouter Deconinck
+
+#pragma once
+
+#include <algorithms/algorithm.h>
+#include <edm4eic/MCRecoParticleAssociation.h>
+#include <edm4eic/ReconstructedParticle.h>
+#include <edm4hep/MCParticle.h>
+#include <edm4hep/utils/vector_utils.h>
+#include <string>
+#include <string_view>
+
+#include "TruthinessConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using TruthinessAlgorithm = algorithms::Algorithm<
+    algorithms::Input<edm4hep::MCParticleCollection, edm4eic::ReconstructedParticleCollection,
+                      edm4eic::MCRecoParticleAssociationCollection>,
+    algorithms::Output<>>;
+
+class Truthiness : public TruthinessAlgorithm, public WithPodConfig<TruthinessConfig> {
+
+public:
+  Truthiness(std::string_view name)
+      : TruthinessAlgorithm{
+            name,
+            {"inputMCParticles", "inputReconstructedParticles", "inputAssociations"},
+            {},
+            "Calculate truthiness metric comparing reconstructed particles to MC "
+            "truth."} {}
+
+  void init() final {};
+  void process(const Input&, const Output&) const final;
+};
+
+} // namespace eicrecon

src/algorithms/reco/TruthinessConfig.h

@@ -0,0 +1,12 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Wouter Deconinck
+
+#pragma once
+
+namespace eicrecon {
+
+struct TruthinessConfig {
+  // No configuration parameters needed at this time
+};
+
+} // namespace eicrecon

src/factories/reco/Truthiness_factory.h

@@ -0,0 +1,44 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Wouter Deconinck
+
+#pragma once
+
+#include <JANA/JEvent.h>
+#include <memory>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "algorithms/reco/Truthiness.h"
+#include "extensions/jana/JOmniFactory.h"
+#include "services/algorithms_init/AlgorithmsInit_service.h"
+
+namespace eicrecon {
+
+class Truthiness_factory : public JOmniFactory<Truthiness_factory> {
+
+public:
+  using FactoryT = JOmniFactory<Truthiness_factory>;
+
+private:
+  std::unique_ptr<Truthiness> m_algo;
+
+  FactoryT::PodioInput<edm4hep::MCParticle> m_mc_particles_input{this};
+  FactoryT::PodioInput<edm4eic::ReconstructedParticle> m_rc_particles_input{this};
+  FactoryT::PodioInput<edm4eic::MCRecoParticleAssociation> m_associations_input{this};
+
+  FactoryT::Service<AlgorithmsInit_service> m_algorithmsInit{this};
+
+public:
+  void Configure() {
+    m_algo = std::make_unique<Truthiness>(this->GetPrefix());
+    m_algo->level(static_cast<algorithms::LogLevel>(this->logger()->level()));
+    m_algo->init();
+  }
+
+  void Process(int32_t /* run_number */, uint64_t /* event_number */) {
+    m_algo->process({m_mc_particles_input(), m_rc_particles_input(), m_associations_input()}, {});
+  }
+};
+
+} // namespace eicrecon

src/global/reco/reco.cc

@@ -42,6 +42,7 @@
 #include "factories/reco/ScatteredElectronsTruth_factory.h"
 #include "factories/reco/TrackClusterMatch_factory.h"
 #include "factories/reco/TransformBreitFrame_factory.h"
+#include "factories/reco/Truthiness_factory.h"
 #include "factories/reco/UndoAfterBurnerMCParticles_factory.h"
 
 extern "C" {
@@ -267,5 +268,10 @@ void InitPlugin(JApplication* app) {
 
   app->Add(new JOmniFactoryGeneratorT<PrimaryVertices_factory>(
       "PrimaryVertices", {"CentralTrackVertices"}, {"PrimaryVertices"}, {}, app));
+
+  // Truthiness metric for event quality assessment
+  app->Add(new JOmniFactoryGeneratorT<Truthiness_factory>(
+      "Truthiness", {"MCParticles", "ReconstructedParticles", "ReconstructedParticleAssociations"},
+      {}, app));
 }
 } // extern "C"