Implement Calo Remnant Combiner (PFA2)

src/algorithms/CMakeLists.txt

@@ -3,6 +3,7 @@ cmake_minimum_required(VERSION 3.16)
 add_subdirectory(interfaces)
 add_subdirectory(calorimetry)
 add_subdirectory(tracking)
+add_subdirectory(particle)
 add_subdirectory(pid)
 add_subdirectory(pid_lut)
 add_subdirectory(digi)

src/algorithms/particle/CMakeLists.txt

@@ -0,0 +1,18 @@
+cmake_minimum_required(VERSION 3.16)
+
+set(PLUGIN_NAME "algorithms_particle")
+
+# Function creates ${PLUGIN_NAME}_plugin and ${PLUGIN_NAME}_library targets
+# Setting default includes, libraries and installation paths
+plugin_add(${PLUGIN_NAME} WITH_SHARED_LIBRARY WITHOUT_PLUGIN)
+
+# The macro grabs sources as *.cc *.cpp *.c and headers as *.h *.hh *.hpp Then
+# correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds
+# headers to the correct installation directory
+plugin_glob_all(${PLUGIN_NAME})
+
+# Find dependencies
+plugin_add_algorithms(${PLUGIN_NAME})
+plugin_add_dd4hep(${PLUGIN_NAME})
+plugin_add_event_model(${PLUGIN_NAME})
+plugin_add_eigen3(${PLUGIN_NAME})

src/algorithms/particle/CaloRemnantCombiner.cc

@@ -0,0 +1,154 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Subhadip Pal
+
+#include <Evaluator/DD4hepUnits.h>
+#include <algorithms/logger.h>
+#include <edm4eic/ClusterCollection.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+
+#include <edm4hep/Vector3f.h>
+#include <edm4hep/utils/vector_utils.h>
+#include <fmt/core.h>
+#include <podio/ObjectID.h>
+#include <cmath>
+#include <gsl/pointers>
+#include <set>
+
+#include "algorithms/particle/CaloRemnantCombinerConfig.h"
+
+#include "CaloRemnantCombiner.h"
+
+namespace eicrecon {
+
+void CaloRemnantCombiner::process(const CaloRemnantCombiner::Input& input,
+                                  const CaloRemnantCombiner::Output& output) const {
+
+  const auto [calo_clusters]    = input;
+  auto [out_neutral_candidates] = output;
+
+  std::vector<bool> visits_ecal(calo_clusters[0]->size(), false);
+  std::vector<bool> visits_hcal(calo_clusters[1]->size(), false);
+
+  std::vector<bool> visits_ecal_true(calo_clusters[0]->size(), true);
+  std::vector<bool> visits_hcal_true(calo_clusters[1]->size(), true);
+
+  while (visits_ecal != visits_ecal_true) {
+
+    edm4eic::MutableReconstructedParticle neutral_candidate_eh;
+
+    // Step 1: Find the seed Ecal cluster with highest energy
+    std::size_t seed_ecal_index = findSeedCluster_index(*calo_clusters[0], visits_ecal);
+
+    if (seed_ecal_index == -1) {
+      info("No Seed Ecal cluster found for remnant combination.");
+    }
+
+    if (seed_ecal_index != -1) {
+      // Get the cluster indices for merging
+      std::set<std::size_t> ecalcluster_indices = getcluster_indices_for_merging(
+          *calo_clusters[0], visits_ecal, seed_ecal_index, m_cfg.delta_r_add_em, *calo_clusters[0]);
+
+      for (const auto& idx : ecalcluster_indices) {
+        neutral_candidate_eh.addToClusters((*calo_clusters[0])[idx]);
+      }
+
+      std::set<std::size_t> e_hcalcluster_indices = getcluster_indices_for_merging(
+          *calo_clusters[1], visits_hcal, seed_ecal_index, m_cfg.delta_r_add_h, *calo_clusters[0]);
+
+      for (const auto& idx : e_hcalcluster_indices) {
+        neutral_candidate_eh.addToClusters((*calo_clusters[1])[idx]);
+      }
+
+      out_neutral_candidates->push_back(neutral_candidate_eh);
+
+    } // end of if (seed_ecal_index != -1)
+
+  } // end of while (visits_ecal != visits_ecal_true)
+
+  while (visits_hcal != visits_hcal_true) {
+
+    edm4eic::MutableReconstructedParticle neutral_candidate_h;
+    std::size_t seed_rem_hcal_index = findSeedCluster_index(*calo_clusters[1], visits_hcal);
+
+    if (seed_rem_hcal_index == -1) {
+      info("No Seed Hcal cluster found for remnant combination.");
+    }
+
+    if (seed_rem_hcal_index != -1) {
+
+      std::set<std::size_t> rem_hcalcluster_indices;
+
+      rem_hcalcluster_indices =
+          getcluster_indices_for_merging(*calo_clusters[1], visits_hcal, seed_rem_hcal_index,
+                                         m_cfg.delta_r_add_h, *calo_clusters[1]);
+
+      for (const auto& idx : rem_hcalcluster_indices) {
+        neutral_candidate_h.addToClusters((*calo_clusters[1])[idx]);
+      }
+      out_neutral_candidates->push_back(neutral_candidate_h);
+
+    } // end of if (seed_rem_hcal_index != -1)
+
+  } // end of while (visits_hcal != visits_hcal_true)
+
+} // end of process
+
+std::size_t CaloRemnantCombiner::findSeedCluster_index(const edm4eic::ClusterCollection& clusters,
+                                                       std::vector<bool>& visits) {
+  double max_cluster_energy      = -1.0;
+  std::size_t seed_cluster_index = -1;
+  for (std::size_t i = 0; i < clusters.size(); ++i) {
+
+    if (visits[i]) {
+      continue;
+    }
+
+    if (clusters[i].getEnergy() > max_cluster_energy) {
+      max_cluster_energy = clusters[i].getEnergy();
+      seed_cluster_index = i;
+    }
+  }
+  return seed_cluster_index;
+}
+
+std::set<std::size_t> CaloRemnantCombiner::getcluster_indices_for_merging(
+    const edm4eic::ClusterCollection& clusters, std::vector<bool>& visits,
+    std::size_t seed_cluster_index, double delta_r, const edm4eic::ClusterCollection& seed) {
+  std::set<std::size_t> cluster_indices;
+
+  for (std::size_t i = 0; i < clusters.size(); ++i) {
+    if (visits[i]) {
+      continue;
+    }
+    // get the distance between current cluster and seed cluster
+
+    // Using simple Euclidean distance in 3D space
+    ///double distance = edm4hep::utils::magnitude(clusters[i].getPosition() - seed[seed_cluster_index].getPosition());
+
+    // Using delta R in the transverse plane (x-y plane)
+    /*double dx = clusters[i].getPosition().x - seed[seed_cluster_index].getPosition().x;
+    double dy = clusters[i].getPosition().y - seed[seed_cluster_index].getPosition().y;
+    double distance = std::sqrt(dx * dx + dy * dy);*/
+
+    // Using angular distance (delta R) in eta-phi space
+    edm4hep::Vector3f seed_pos    = seed[seed_cluster_index].getPosition();
+    edm4hep::Vector3f cluster_pos = clusters[i].getPosition();
+
+    float eta_seed    = edm4hep::utils::eta(seed_pos);
+    float phi_seed    = edm4hep::utils::angleAzimuthal(seed_pos);
+    float eta_cluster = edm4hep::utils::eta(cluster_pos);
+    float phi_cluster = edm4hep::utils::angleAzimuthal(cluster_pos);
+
+    float dphi     = phi_cluster - phi_seed;
+    float deta     = eta_cluster - eta_seed;
+    float distance = std::sqrt(deta * deta + dphi * dphi);
+
+    if (distance < delta_r) { // distance threshold for merging
+      cluster_indices.insert(i);
+      visits[i] = true;
+    }
+  }
+  return cluster_indices;
+}
+
+} // namespace eicrecon

src/algorithms/particle/CaloRemnantCombiner.h

@@ -0,0 +1,41 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Subhadip Pal
+
+#pragma once
+
+#include <algorithms/algorithm.h>
+#include <edm4eic/ClusterCollection.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+#include <string>
+#include <string_view>
+
+#include "CaloRemnantCombinerConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using CaloRemnantCombinerAlgorithm =
+    algorithms::Algorithm<algorithms::Input<std::vector<edm4eic::ClusterCollection>>,
+                          algorithms::Output<edm4eic::ReconstructedParticleCollection>>;
+
+class CaloRemnantCombiner : public CaloRemnantCombinerAlgorithm,
+                            public WithPodConfig<CaloRemnantCombinerConfig> {
+
+public:
+  CaloRemnantCombiner(std::string_view name)
+      : CaloRemnantCombinerAlgorithm{name,
+                                     {"CaloClusters"},
+                                     {"NeutralParticleCandidate"},
+                                     "make neutral candidates from remnant clusters"} {}
+
+  void init() final {};
+  void process(const Input&, const Output&) const final;
+  static std::size_t findSeedCluster_index(const edm4eic::ClusterCollection& clusters,
+                                           std::vector<bool>& visits);
+  static std::set<std::size_t>
+  getcluster_indices_for_merging(const edm4eic::ClusterCollection& clusters,
+                                 std::vector<bool>& visits, std::size_t seed_cluster_index,
+                                 double delta_r, const edm4eic::ClusterCollection& seed);
+};
+
+} // namespace eicrecon

src/algorithms/particle/CaloRemnantCombinerConfig.h

@@ -0,0 +1,16 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Subhadip Pal
+
+#pragma once
+
+#include <Evaluator/DD4hepUnits.h>
+
+namespace eicrecon {
+
+struct CaloRemnantCombinerConfig {
+
+  double delta_r_add_em = 0.03;
+  double delta_r_add_h  = 0.15;
+};
+
+} // namespace eicrecon

src/factories/CMakeLists.txt

@@ -2,6 +2,7 @@ add_subdirectory(calorimetry)
 add_subdirectory(digi)
 add_subdirectory(fardetectors)
 add_subdirectory(meta)
+add_subdirectory(particle)
 add_subdirectory(pid)
 add_subdirectory(pid_lut)
 add_subdirectory(reco)

src/factories/particle/CMakeLists.txt

@@ -0,0 +1,2 @@
+set(PLUGIN_NAME "factories_particle")
+plugin_headers_only(${PLUGIN_NAME})

src/factories/particle/CaloRemnantCombiner_factory.h

@@ -0,0 +1,65 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Subhadip Pal
+
+#pragma once
+
+#include <algorithms/logger.h>
+#include <edm4eic/ClusterCollection.h>
+#include <edm4eic/ReconstructedParticleCollection.h>
+
+#include <spdlog/logger.h>
+#include <stdint.h>
+#include <memory>
+
+#include "algorithms/particle/CaloRemnantCombiner.h"
+#include "extensions/jana/JOmniFactory.h"
+#include "services/algorithms_init/AlgorithmsInit_service.h"
+
+namespace eicrecon {
+
+class CaloRemnantCombiner_factory
+    : public JOmniFactory<CaloRemnantCombiner_factory, CaloRemnantCombinerConfig> {
+private:
+  // Underlying algorithm
+  std::unique_ptr<eicrecon::CaloRemnantCombiner> m_algo;
+
+  // Declare inputs
+  VariadicPodioInput<edm4eic::Cluster> m_in_calo_clusters{this};
+
+  // Declare outputs
+  PodioOutput<edm4eic::ReconstructedParticle> m_out_neutral_candidates{this};
+
+  // Declare parameters
+  /*ParameterRef<double> m_min_energy_over_momentum{this, "minEnergyOverMomentum",
+                                                  config().min_energy_over_momentum};
+  ParameterRef<double> m_max_energy_over_momentum{this, "maxEnergyOverMomentum",
+                                                  config().max_energy_over_momentum};*/
+
+public:
+  void Configure() {
+    // This is called when the factory is instantiated.
+    // Use this callback to make sure the algorithm is configured.
+    // The logger, parameters, and services have all been fetched before this is called
+    m_algo = std::make_unique<CaloRemnantCombiner>(GetPrefix());
+    m_algo->level(static_cast<algorithms::LogLevel>(logger()->level()));
+    // Pass config object to algorithm
+    m_algo->applyConfig(config());
+
+    m_algo->init();
+  }
+
+  void Process(int32_t /* run_number */, uint64_t /* event_number */) {
+    // This is called on every event.
+    // Use this callback to call your Algorithm using all inputs and outputs
+    // The inputs will have already been fetched for you at this point.
+    auto in1 = m_in_calo_clusters();
+    std::vector<gsl::not_null<const edm4eic::ClusterCollection*>> in2;
+
+    std::copy(in1.cbegin(), in1.cend(), std::back_inserter(in2));
+    m_algo->process({in2}, {m_out_neutral_candidates().get()});
+
+    logger()->debug("Found {} reconstructed neutral candidates",
+                    m_out_neutral_candidates()->size());
+  }
+};
+} // namespace eicrecon

src/global/CMakeLists.txt

@@ -2,6 +2,7 @@ cmake_minimum_required(VERSION 3.16)
 
 add_subdirectory(tracking)
 add_subdirectory(reco)
+add_subdirectory(particle)
 add_subdirectory(pid)
 add_subdirectory(pid_lut)
 add_subdirectory(beam)

src/global/particle/CMakeLists.txt

@@ -0,0 +1,25 @@
+cmake_minimum_required(VERSION 3.16)
+
+# Automatically set plugin name the same as the directory name Don't forget
+# string(REPLACE " " "_" PLUGIN_NAME ${PLUGIN_NAME}) if this dir has spaces in
+# its name
+get_filename_component(PLUGIN_NAME ${CMAKE_CURRENT_LIST_DIR} NAME)
+
+# Function creates ${PLUGIN_NAME}_plugin and ${PLUGIN_NAME}_library targets
+# Setting default includes, libraries and installation paths
+plugin_add(${PLUGIN_NAME} PLUGIN_USE_CC_ONLY)
+
+# The macro grabs sources as *.cc *.cpp *.c and headers as *.h *.hh *.hpp Then
+# correctly sets sources for ${_name}_plugin and ${_name}_library targets Adds
+# headers to the correct installation directory
+plugin_glob_all(${PLUGIN_NAME})
+
+# Find dependencies
+plugin_add_dd4hep(${PLUGIN_NAME})
+plugin_add_event_model(${PLUGIN_NAME})
+
+# Add include directories (works same as target_include_directories)
+# plugin_include_directories(${PLUGIN_NAME} SYSTEM PUBLIC ...)
+
+# Add libraries (works same as target_include_directories)
+plugin_link_libraries(${PLUGIN_NAME} algorithms_particle_library)