Skip to content

Implement Calo Remnant Combiner (PFA2) #2195

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

Implement Calo Remnant Combiner (PFA2) #2195

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ad853fe
PFA2 : first implementation
palsgit Nov 19, 2025
5f26b12
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Nov 19, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
113 changes: 60 additions & 53 deletions src/algorithms/particle/CaloRemnantCombiner.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,11 @@

namespace eicrecon {

void CaloRemnantCombiner::process(const CaloRemnantCombiner::Input& input, const CaloRemnantCombiner::Output& output) const {
void CaloRemnantCombiner::process(const CaloRemnantCombiner::Input& input,
Copy link
Contributor

@ruse-traveler ruse-traveler Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Definitely want a docstring here! This is a good place for a high-level summary of what the algorithm is doing!

Suggested change
void CaloRemnantCombiner::process(const CaloRemnantCombiner::Input& input,
// ----------------------------------------------------------------------------
//! Process inputs
// ----------------------------------------------------------------------------
/*! Construct a candidate neutral particle via the
* following algorithm.
* 1. Repeat the following the steps until every EMCal
* cluster has been used:
* a. Identify seed EMCal cluster
* b. Identify all EMCal clusters and HCal clusters which
* lie within a radius of deltaRAddEM and deltaRAddH
* around seed EMCal cluster respectively
* c. Combine all identified clusters into a neutral particle
* candidate
* 2. Repeat the following steps until every HCal
* cluster has been used:
* a. Identify seed HCal cluster
* b. Identify all HCal clusters which lie within a
* radius of deltaRAddH around seed HCal
* cluster
* c. Combine all identified clusters into a neutral particle
* candidate
*/
void CaloRemnantCombiner::process(const CaloRemnantCombiner::Input& input,

const CaloRemnantCombiner::Output& output) const {

const auto [calo_clusters] = input;
auto [out_neutral_candidates] = output;
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);
Expand All @@ -33,66 +34,70 @@ void CaloRemnantCombiner::process(const CaloRemnantCombiner::Input& input, const

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);
edm4eic::MutableReconstructedParticle neutral_candidate_eh;

if (seed_ecal_index == -1) {
info("No Seed Ecal cluster found for remnant combination.");
}
// 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) {
// 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]);
if (seed_ecal_index == -1) {
info("No Seed Ecal cluster found for remnant combination.");
Copy link
Member

@veprbl veprbl Nov 20, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
info("No Seed Ecal cluster found for remnant combination.");
debug("No Seed Ecal cluster found for remnant combination");

info messages are too "loud" for a common algorithm message

}

for (const auto& idx : ecalcluster_indices) {
neutral_candidate_eh.addToClusters((*calo_clusters[0])[idx]);
}
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]);

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 : ecalcluster_indices) {
neutral_candidate_eh.addToClusters((*calo_clusters[0])[idx]);
}

for (const auto& idx : e_hcalcluster_indices) {
neutral_candidate_eh.addToClusters((*calo_clusters[1])[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]);

out_neutral_candidates->push_back(neutral_candidate_eh);
for (const auto& idx : e_hcalcluster_indices) {
neutral_candidate_eh.addToClusters((*calo_clusters[1])[idx]);
}

} // end of if (seed_ecal_index != -1)
out_neutral_candidates->push_back(neutral_candidate_eh);

} // end of while (visits_ecal != visits_ecal_true)
} // 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.");
}
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) {
if (seed_rem_hcal_index == -1) {
info("No Seed Hcal cluster found for remnant combination.");
}

std::set<std::size_t> rem_hcalcluster_indices;
if (seed_rem_hcal_index != -1) {

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]);
std::set<std::size_t> rem_hcalcluster_indices;

for (const auto& idx : rem_hcalcluster_indices) {
neutral_candidate_h.addToClusters((*calo_clusters[1])[idx]);
}
out_neutral_candidates->push_back(neutral_candidate_h);
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 if (seed_rem_hcal_index != -1)

} // end of while (visits_hcal != visits_hcal_true)

} // end of process
} // 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) {
std::size_t CaloRemnantCombiner::findSeedCluster_index(const edm4eic::ClusterCollection& clusters,
Copy link
Contributor

@ruse-traveler ruse-traveler Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another good place for a docstring! Here you'll definitely want to note the criteria for a seed cluster and what it's for.

Suggested change
std::size_t CaloRemnantCombiner::findSeedCluster_index(const edm4eic::ClusterCollection& clusters,
// ----------------------------------------------------------------------------
//! Find seed cluster
// ----------------------------------------------------------------------------
/*! Identifies a seed (highest energy) cluster in a collection
* which sets the center of the cone in which clusters are
* combined.
*/
std::size_t CaloRemnantCombiner::find_seed_cluster_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;
Expand All @@ -102,12 +107,14 @@ std::size_t CaloRemnantCombiner::findSeedCluster_index(const edm4eic::ClusterCol
max_cluster_energy = clusters[i].getEnergy();
seed_cluster_index = i;
}
}
return seed_cluster_index;
}
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;
std::set<std::size_t> CaloRemnantCombiner::getcluster_indices_for_merging(
Copy link
Contributor

@ruse-traveler ruse-traveler Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is another good place for a docstring: you'll want to briefly explain what the method is doing!

Suggested change
std::set<std::size_t> CaloRemnantCombiner::getcluster_indices_for_merging(
// ----------------------------------------------------------------------------
//! Find cluster indices for merging
// ----------------------------------------------------------------------------
/*! Creates a set of indices corresponding to clusters which lie within
* a radius of `delta_r_add` around the seed cluster with index
* `seed_cluster_index`.
*/
std::set<std::size_t> CaloRemnantCombiner::get_cluster_indices_for_merging(

const edm4eic::ClusterCollection& clusters, std::vector<bool>& visits,
std::size_t seed_cluster_index, double delta_r, const edm4eic::ClusterCollection& seed) {
Comment on lines +115 to +116
Copy link
Contributor

@ruse-traveler ruse-traveler Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd suggest making the delta_r argument here parallel the config name (makes the code a little more self-describing).

Suggested change
const edm4eic::ClusterCollection& clusters, std::vector<bool>& visits,
std::size_t seed_cluster_index, double delta_r, const edm4eic::ClusterCollection& seed) {
const edm4eic::ClusterCollection& clusters, std::vector<bool>& visits,
std::size_t seed_cluster_index, double delta_r_add, const edm4eic::ClusterCollection& seed) {

std::set<std::size_t> cluster_indices;

for (std::size_t i = 0; i < clusters.size(); ++i) {
if (visits[i]) {
Expand All @@ -117,23 +124,23 @@ std::set<std::size_t> CaloRemnantCombiner::getcluster_indices_for_merging(const

// 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 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_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 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
Copy link
Contributor

@ruse-traveler ruse-traveler Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Following up on the above comment: this way I think you can drop the comment after the braces!

Suggested change
if (distance < delta_r) { // distance threshold for merging
if (distance < delta_r_add) {

Expand Down
24 changes: 14 additions & 10 deletions src/algorithms/particle/CaloRemnantCombiner.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,24 +14,28 @@

namespace eicrecon {

using CaloRemnantCombinerAlgorithm = algorithms::Algorithm<
algorithms::Input<std::vector<edm4eic::ClusterCollection>>,
algorithms::Output<edm4eic::ReconstructedParticleCollection>>;
using CaloRemnantCombinerAlgorithm =
Copy link
Contributor

@ruse-traveler ruse-traveler Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This would be another good place for a docstring, especially since the ordering of collections in the input vector matters!

Suggested change
using CaloRemnantCombinerAlgorithm =
// --------------------------------------------------------------------------
//! Algorithm input/output
// --------------------------------------------------------------------------
/*! Input is a vector of calorimeter cluster collections. For now:
* - 1st entry in the vector should be the EMCal collection, and
* - 2nd entry in the vector should be the HCal collection.
* This can be generalized in the future.
*/
using CaloRemnantCombinerAlgorithm =

algorithms::Algorithm<algorithms::Input<std::vector<edm4eic::ClusterCollection>>,
algorithms::Output<edm4eic::ReconstructedParticleCollection>>;

class CaloRemnantCombiner : public CaloRemnantCombinerAlgorithm,
Copy link
Contributor

@ruse-traveler ruse-traveler Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd suggest adding docstrings in several places. These should get propagated to our doxygen page, and help explain what's happening in the code for anyone who's looking at the source code!

Suggested change
class CaloRemnantCombiner : public CaloRemnantCombinerAlgorithm,
// ==========================================================================
//! Calorimeter Remnant Cluster Combiner
// ==========================================================================
/*! An algorithm which takes multiple calorimeter cluster collections and combines them into
* neutral-particle candidates based on distance matching.
*/
class CaloRemnantCombiner : public CaloRemnantCombinerAlgorithm,

public WithPodConfig<CaloRemnantCombinerConfig> {
public WithPodConfig<CaloRemnantCombinerConfig> {

public:
CaloRemnantCombiner(std::string_view name)
: CaloRemnantCombinerAlgorithm{name,
{"CaloClusters"},
{"NeutralParticleCandidate"},
"make neutral candidates from remnant clusters"} {}
{"CaloClusters"},
{"NeutralParticleCandidate"},
"make neutral candidates from remnant clusters"} {}

void init() final {};
Copy link
Contributor

@ruse-traveler ruse-traveler Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since you're not overwriting this method, I'd suggest dropping it from here!

Suggested change
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);
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);
Comment on lines +33 to +38
Copy link
Contributor

@ruse-traveler ruse-traveler Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd suggest flagging these as const and being consistent with snake- vs. camel-case!

Suggested change
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);
static std::size_t find_seed_cluster_index(const edm4eic::ClusterCollection& clusters,
std::vector<bool>& visits) const;
static std::set<std::size_t>
get_cluster_indices_for_merging(const edm4eic::ClusterCollection& clusters,
std::vector<bool>& visits, std::size_t seed_cluster_index,
double delta_r, const edm4eic::ClusterCollection& seed) const;

};

} // namespace eicrecon
} // namespace eicrecon
3 changes: 1 addition & 2 deletions src/algorithms/particle/CaloRemnantCombinerConfig.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,7 @@ namespace eicrecon {
struct CaloRemnantCombinerConfig {

double delta_r_add_em = 0.03;
double delta_r_add_h = 0.15;

double delta_r_add_h = 0.15;
};

} // namespace eicrecon
2 changes: 1 addition & 1 deletion src/factories/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,4 +6,4 @@ add_subdirectory(particle)
add_subdirectory(pid)
add_subdirectory(pid_lut)
add_subdirectory(reco)
add_subdirectory(tracking)
add_subdirectory(tracking)
9 changes: 5 additions & 4 deletions src/factories/particle/CaloRemnantCombiner_factory.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,10 @@
#include "extensions/jana/JOmniFactory.h"
#include "services/algorithms_init/AlgorithmsInit_service.h"


namespace eicrecon {

class CaloRemnantCombiner_factory : public JOmniFactory<CaloRemnantCombiner_factory, CaloRemnantCombinerConfig> {
class CaloRemnantCombiner_factory
: public JOmniFactory<CaloRemnantCombiner_factory, CaloRemnantCombinerConfig> {
private:
// Underlying algorithm
std::unique_ptr<eicrecon::CaloRemnantCombiner> m_algo;
Expand Down Expand Up @@ -58,7 +58,8 @@ class CaloRemnantCombiner_factory : public JOmniFactory<CaloRemnantCombiner_fact
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());
logger()->debug("Found {} reconstructed neutral candidates",
m_out_neutral_candidates()->size());
}
};
} // namespace eicrecon
} // namespace eicrecon
2 changes: 1 addition & 1 deletion src/global/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,4 @@ add_subdirectory(reco)
add_subdirectory(particle)
add_subdirectory(pid)
add_subdirectory(pid_lut)
add_subdirectory(beam)
add_subdirectory(beam)
4 changes: 2 additions & 2 deletions src/global/particle/particle.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,8 @@ void InitPlugin(JApplication* app) {
/* TODO add PFA2 EEEMCal here */
/* TODO add PFA2 EHCal here */
app->Add(new JOmniFactoryGeneratorT<CaloRemnantCombiner_factory>(
"ReconstructedNeutralCandidates", {"EcalEndcapNClusters","HcalEndcapNClusters"}, {"ReconstructedNeutralCandidates"}, {}, app));
"ReconstructedNeutralCandidates", {"EcalEndcapNClusters", "HcalEndcapNClusters"},
{"ReconstructedNeutralCandidates"}, {}, app));

// central ------------------------------------------------------------

Expand Down Expand Up @@ -114,6 +115,5 @@ void InitPlugin(JApplication* app) {
/* TODO add PFA3 FHCal insert here */

/* TODO collect reconstructed particles here */

}
} // extern "C"
Loading