diff --git a/src/algorithms/calorimetry/TrackClusterMergeSplitter.cc b/src/algorithms/calorimetry/TrackClusterMergeSplitter.cc
index 4fdd64359..12cf9a739 100644
--- a/src/algorithms/calorimetry/TrackClusterMergeSplitter.cc
+++ b/src/algorithms/calorimetry/TrackClusterMergeSplitter.cc
@@ -1,7 +1,9 @@
// SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (C) 2024 Derek Anderson
-#include <edm4eic/CalorimeterHit.h>
+#include <edm4hep/MCParticle.h>
+#include <edm4hep/RawCalorimeterHit.h>
+#include <edm4hep/SimCalorimeterHit.h>
#include <edm4hep/Vector3f.h>
#include <edm4hep/utils/vector_utils.h>
#include <fmt/core.h>
@@ -37,14 +39,14 @@ namespace eicrecon {
//! Process inputs
// --------------------------------------------------------------------------
/*! Primary algorithm call: algorithm ingests a collection
- * protoclusters and a collection of track projections.
- * It then decides to merge or split protoclusters according
- * to the following algorithm:
+ * clusters and a collection of track projections. It then
+ * decides to merge or split clusters according to the
+ * following algorithm:
* 1. Identify all tracks projections pointing to the
* specified calorimeter.
- * 2. Match relevant track projections to protoclusters
+ * 2. Match relevant track projections to clusters
* based on distance between projection and the energy-
- * weighted barycenter of the protocluster;
+ * weighted barycenter of the cluster;
* 3. For each cluster-track pair:
* i. Calculate the significance of the pair's
* E/p relative to the provided mean E/p and
@@ -52,10 +54,10 @@ namespace eicrecon {
* ii. If the significance is less than the
* significance specified by `minSigCut`,
* merge all clusters within `drAdd`.
- * 4. Create a protocluster for each merged cluster
- * and copy all unused protoclusters into output.
+ * 4. Create a cluster for each merged cluster
+ * and copy all unused clusters into output.
* - If multiple tracks point to the same merged
- * cluster, create a new protocluster for each
+ * cluster, create a new cluster for each
* projection with hit weighted relative to
* the track momentum.
*/
@@ -65,30 +67,40 @@ namespace eicrecon {
) const {
// grab inputs/outputs
- const auto [in_protoclusters, in_projections] = input;
- auto [out_protoclusters] = output;
+#if EDM4EIC_VERSION_MAJOR >= 7
+ const auto [in_clusters, in_projections, in_clust_associations, in_hit_associations] = input;
+#else
+ const auto [in_clusters, in_projections, in_clust_associations, in_sim_hits] = input;
+#endif
+#if EDM4EIC_VERSION_MAJOR >= 8
+ auto [out_clusters, out_matches, out_clust_associations] = output;
+#else
+ auto [out_clusters, out_clust_associations] = output;
+#endif
// exit if no clusters in collection
- if (in_protoclusters->size() == 0) {
- debug("No proto-clusters in input collection.");
+ if (in_clusters->size() == 0) {
+ debug("No clusters in input collection.");
return;
}
// ------------------------------------------------------------------------
// 1. Identify projections to calorimeter
// ------------------------------------------------------------------------
+ VecTrk vecTrack;
VecProj vecProject;
- get_projections(in_projections, vecProject);
+ get_projections(in_projections, vecProject, vecTrack);
// ------------------------------------------------------------------------
// 2. Match relevant projections to clusters
// ------------------------------------------------------------------------
+ MapToVecTrk mapTrkToMatch;
MapToVecProj mapProjToSplit;
if (vecProject.size() == 0) {
debug("No projections to match clusters to.");
return;
} else {
- match_clusters_to_tracks(in_protoclusters, vecProject, mapProjToSplit);
+ match_clusters_to_tracks(in_clusters, vecProject, vecTrack, mapProjToSplit, mapTrkToMatch);
}
// ------------------------------------------------------------------------
@@ -112,7 +124,7 @@ namespace eicrecon {
setUsedClust.insert( clustSeed );
// grab cluster energy and projection momentum
- const float eClustSeed = get_cluster_energy(clustSeed);
+ const float eClustSeed = clustSeed.getEnergy();
const float eProjSeed = m_cfg.avgEP * edm4hep::utils::magnitude(projSeed.momentum);
// ----------------------------------------------------------------------
@@ -130,14 +142,13 @@ namespace eicrecon {
}
// get eta, phi of seed
- const auto posSeed = get_cluster_position(clustSeed);
- const float etaSeed = edm4hep::utils::eta(posSeed);
- const float phiSeed = edm4hep::utils::angleAzimuthal(posSeed);
+ const float etaSeed = edm4hep::utils::eta(clustSeed.getPosition());
+ const float phiSeed = edm4hep::utils::angleAzimuthal(clustSeed.getPosition());
// loop over other clusters
float eClustSum = eClustSeed;
float sigSum = sigSeed;
- for (auto in_cluster : *in_protoclusters) {
+ for (auto in_cluster : *in_clusters) {
// ignore used clusters
if (setUsedClust.count(in_cluster)) {
@@ -145,9 +156,8 @@ namespace eicrecon {
}
// get eta, phi of cluster
- const auto posClust = get_cluster_position(in_cluster);
- const float etaClust = edm4hep::utils::eta(posClust);
- const float phiClust = edm4hep::utils::angleAzimuthal(posClust);
+ const float etaClust = edm4hep::utils::eta(in_cluster.getPosition());
+ const float phiClust = edm4hep::utils::angleAzimuthal(in_cluster.getPosition());
// get distance to seed
const float drToSeed = std::hypot(
@@ -177,7 +187,7 @@ namespace eicrecon {
}
// increment sums and output debugging
- eClustSum += get_cluster_energy(in_cluster);
+ eClustSum += in_cluster.getEnergy();
sigSum = (eClustSum - eProjSeed) / m_cfg.sigEP;
trace(
"{} clusters to merge: current sum = {}, current significance = {}, {} track(s) pointing to merged cluster",
@@ -190,21 +200,55 @@ namespace eicrecon {
} // end matched cluster-projection loop
// ------------------------------------------------------------------------
- // 4. Create an output protocluster for each merged cluster and for
+ // 4. Create an output cluster for each merged cluster and for
// each track pointing to merged cluster
// ------------------------------------------------------------------------
for (auto& [clustSeed, vecClustToMerge] : mapClustToMerge) {
+
+ // create a cluster for each projection to merged cluster
+ std::vector<edm4eic::MutableCluster> new_clusters;
+ for (const auto& proj : mapProjToSplit[clustSeed]) {
+ new_clusters.push_back( out_clusters->create() );
+ }
+
+ // merge & split as needed
merge_and_split_clusters(
vecClustToMerge,
mapProjToSplit[clustSeed],
- out_protoclusters
+ new_clusters
);
+
+ // collect associations of merged clusters
+ for (const auto& new_clust : new_clusters) {
+ collect_associations(
+ new_clust,
+ vecClustToMerge,
+ in_clust_associations,
+#if EDM4EIC_VERSION_MAJOR >= 7
+ in_hit_associations,
+#else
+ in_sim_hits,
+#endif
+ out_clust_associations
+ );
+ }
+
+#if EDM4EIC_VERSION_MAJOR >= 8
+ // and finally create a track-cluster match for each pair
+ for (std::size_t iTrk = 0; const auto& trk : mapTrkToMatch[clustSeed]) {
+ edm4eic::MutableTrackClusterMatch match = out_matches->create();
+ match.setCluster( new_clusters[iTrk] );
+ match.setTrack( trk );
+ match.setWeight( 1.0 ); // FIXME placeholder
+ trace("Matched output cluster {} to track {}", new_clusters[iTrk].getObjectID().index, trk.getObjectID().index);
+ }
+#endif
} // end clusters to merge loop
// ------------------------------------------------------------------------
// copy unused clusters to output
// ------------------------------------------------------------------------
- for (auto in_cluster : *in_protoclusters) {
+ for (auto in_cluster : *in_clusters) {
// ignore used clusters
if (setUsedClust.count(in_cluster)) {
@@ -212,13 +256,25 @@ namespace eicrecon {
}
// copy cluster and add to output collection
- edm4eic::MutableProtoCluster out_cluster = in_cluster.clone();
- out_protoclusters->push_back(out_cluster);
+ edm4eic::MutableCluster out_cluster = in_cluster.clone();
+ out_clusters->push_back(out_cluster);
trace("Copied input cluster {} onto output cluster {}",
in_cluster.getObjectID().index,
out_cluster.getObjectID().index
);
+ // if provided, copy corresponding associations to
+ // output collection
+ for (auto in_assoc : *in_clust_associations) {
+ if (in_assoc.getRec() == in_cluster) {
+ edm4eic::MutableMCRecoClusterParticleAssociation out_assoc = in_assoc.clone();
+ out_clust_associations->push_back(out_assoc);
+ trace("Copied input association {} onto output association {}",
+ in_assoc.getObjectID().index,
+ out_assoc.getObjectID().index
+ );
+ }
+ } // end association loop
} // end cluster loop
} // end 'process(Input&, Output&)'
@@ -228,9 +284,12 @@ namespace eicrecon {
// --------------------------------------------------------------------------
//! Collect projections pointing to calorimeter
// --------------------------------------------------------------------------
+ /*! FIXME remove this once cluster-track matching has been centralized
+ */
void TrackClusterMergeSplitter::get_projections(
const edm4eic::TrackSegmentCollection* projections,
- VecProj& relevant_projects
+ VecProj& relevant_projects,
+ VecTrk& relevant_trks
) const {
// return if projections are empty
@@ -247,6 +306,8 @@ namespace eicrecon {
(point.surface == 1)
) {
relevant_projects.push_back(point);
+ relevant_trks.push_back(project.getTrack());
+ break;
}
} // end point loop
} // end projection loop
@@ -262,9 +323,11 @@ namespace eicrecon {
/*! FIXME remove this once cluster-track matching has been centralized
*/
void TrackClusterMergeSplitter::match_clusters_to_tracks(
- const edm4eic::ProtoClusterCollection* clusters,
+ const edm4eic::ClusterCollection* clusters,
const VecProj& projections,
- MapToVecProj& matches
+ const VecTrk& tracks,
+ MapToVecProj& matched_projects,
+ MapToVecTrk& matched_tracks
) const {
@@ -279,7 +342,7 @@ namespace eicrecon {
const float phiProj = edm4hep::utils::angleAzimuthal(project.position);
// to store matched cluster
- edm4eic::ProtoCluster match;
+ edm4eic::Cluster match;
// find closest cluster
bool foundMatch = false;
@@ -287,9 +350,8 @@ namespace eicrecon {
for (auto cluster : *clusters) {
// get eta, phi of cluster
- const auto posClust = get_cluster_position(cluster);
- const float etaClust = edm4hep::utils::eta(posClust);
- const float phiClust = edm4hep::utils::angleAzimuthal(posClust);
+ const float etaClust = edm4hep::utils::eta(cluster.getPosition());
+ const float phiClust = edm4hep::utils::angleAzimuthal(cluster.getPosition());
// calculate distance to centroid
const float dist = std::hypot(
@@ -307,55 +369,47 @@ namespace eicrecon {
// record match if found
if (foundMatch) {
- matches[match].push_back(project);
+ matched_projects[match].push_back(project);
+ matched_tracks[match].push_back(tracks[iProject]);
trace("Matched cluster to track projection: eta-phi distance = {}", dMatch);
}
- } // end cluster loop
- debug("Finished matching clusters to track projections: {} matches", matches.size());
+ } // end projection loop
+ debug("Finished matching clusters to track projections: {} matches", matched_projects.size());
- } // end 'match_clusters_to_tracks(edm4eic::ClusterCollection*, VecTrkPoint&, MapToVecProj&)'
+ } // end 'match_clusters_to_tracks(edm4eic::ClusterCollection*, VecProj&, VecTrk&, MapToVecProj&, MapToVecTrk&)'
// --------------------------------------------------------------------------
//! Merge identified clusters and split if needed
// --------------------------------------------------------------------------
- /*! If multiple tracks are pointing to merged cluster, a new protocluster
- * is created for each track w/ hits weighted by its distance to the track
- * and the track's momentum.
+ /*! If multiple tracks are pointing to merged cluster, a new cluster
+ * is created for each track w/ hits weighted by its distance to
+ * the track and the track's momentum.
*/
void TrackClusterMergeSplitter::merge_and_split_clusters(
const VecClust& to_merge,
const VecProj& to_split,
- edm4eic::ProtoClusterCollection* out_protoclusters
+ std::vector<edm4eic::MutableCluster>& new_clusters
) const {
// if only 1 matched track, no need to split
+ // otherwise split merged cluster for each
+ // matched track
if (to_split.size() == 1) {
- edm4eic::MutableProtoCluster new_clust = out_protoclusters->create();
- for (const auto& old_clust : to_merge) {
- for (const auto& hit : old_clust.getHits()) {
- new_clust.addToHits( hit );
- new_clust.addToWeights( 1. );
- }
- trace("Merged input cluster {} into output cluster {}", old_clust.getObjectID().index, new_clust.getObjectID().index);
- }
+ make_cluster(to_merge, new_clusters.front());
return;
+ } else {
+ trace("Splitting merged cluster across {} tracks", to_split.size());
}
- // otherwise split merged cluster for each matched track
- std::vector<edm4eic::MutableProtoCluster> new_clusters;
- for (const auto& proj : to_split) {
- new_clusters.push_back( out_protoclusters->create() );
- }
- trace("Splitting merged cluster across {} tracks", to_split.size());
-
- // loop over all hits from all clusters to merge
- std::vector<float> weights( to_split.size(), 1. );
+ // calculate weights for splitting
+ VecWeights weights(to_split.size());
for (const auto& old_clust : to_merge) {
for (const auto& hit : old_clust.getHits()) {
- // calculate hit's weight for each track
+ // calculate a weight for each projection
+ double wTotal = 0.;
for (std::size_t iProj = 0; const auto& proj : to_split) {
// get track eta, phi
@@ -373,77 +427,202 @@ namespace eicrecon {
std::remainder(phiHit - phiProj, 2. * M_PI)
);
- // set weight
- weights[iProj] = std::exp(-1. * dist / m_cfg.transverseEnergyProfileScale) * mom;
- ++iProj;
- }
+ // get weight
+ const float weight = std::exp(-1. * dist / m_cfg.transverseEnergyProfileScale) * mom ;
- // normalize weights
- float wTotal = 0.;
- for (const float weight : weights) {
+ // set weight & increment sum of weights
+ weights[iProj][hit] = weight;
wTotal += weight;
- }
- for (float& weight : weights) {
- weight /= wTotal;
+ ++iProj;
}
- // add hit to each split merged cluster w/ relevant weight
- for (std::size_t iProj = 0; auto& new_clust : new_clusters) {
- new_clust.addToHits( hit );
- new_clust.addToWeights( weights[iProj] );
+ // normalize weights over all projections
+ for (std::size_t iProj = 0; iProj < to_split.size(); ++iProj) {
+ weights[iProj][hit] /= wTotal;
}
-
} // end hits to merge loop
} // end clusters to merge loop
- } // end 'merge_and_split_clusters(VecClust&, VecProj&, edm4eic::MutableCluster&)'
+ // make new clusters with relevant splitting weights
+ for (std::size_t iProj = 0; auto& new_clust : new_clusters) {
+ make_cluster(to_merge, new_clust, weights[iProj]);
+ ++iProj;
+ }
+
+ } // end 'merge_and_split_clusters(VecClust&, VecProj&, std::vector<edm4eic::MutableCluster>&)'
// --------------------------------------------------------------------------
- //! Grab current energy of protocluster
+ //! Make new cluster out of old ones
// --------------------------------------------------------------------------
- float TrackClusterMergeSplitter::get_cluster_energy(const edm4eic::ProtoCluster& clust) const {
+ void TrackClusterMergeSplitter::make_cluster(
+ const VecClust& old_clusts,
+ edm4eic::MutableCluster& new_clust,
+ std::optional<MapToWeight> split_weights
+ ) const {
+
+ // determine total no. of hits
+ std::size_t nHits = 0;
+ for (const auto& old_clust : old_clusts) {
+ nHits += old_clust.getNhits();
+ }
+ new_clust.setNhits(nHits);
float eClust = 0.;
- for (auto hit : clust.getHits()) {
- eClust += hit.getEnergy();
+ float wClust = 0.;
+ float tClust = 0.;
+ for (const auto& old_clust : old_clusts) {
+ for (std::size_t iHit = 0; const auto& hit : old_clust.getHits()) {
+
+ // get weight and update if needed
+ float weight = old_clust.getHitContributions()[iHit] / hit.getEnergy();
+ if (split_weights.has_value()) {
+ weight *= split_weights.value().at(hit);
+ }
+
+ // update running tallies
+ eClust += hit.getEnergy() * weight;
+ tClust += (hit.getTime() - tClust) * (hit.getEnergy() / eClust);
+ wClust += weight;
+
+ // add hits and increment counter
+ new_clust.addToHits(hit);
+ new_clust.addToHitContributions(hit.getEnergy() * weight);
+ ++iHit;
+ } // end hit loop
+ trace("Merged input cluster {} into output cluster {}", old_clust.getObjectID().index, new_clust.getObjectID().index);
+ } // end cluster loop
+
+ // update cluster position by taking energy-weighted
+ // average of positions of clusters to merge
+ edm4hep::Vector3f rClust = new_clust.getPosition();
+ for (const auto& old_clust : old_clusts) {
+ rClust = rClust + ((old_clust.getEnergy() / eClust) * old_clust.getPosition());
}
- return eClust / m_cfg.sampFrac;
- } // end 'get_cluster_energy(edm4eic::ProtoCluster&)'
+ // set parameters
+ new_clust.setEnergy(eClust);
+ new_clust.setEnergyError(0.);
+ new_clust.setTime(tClust);
+ new_clust.setTimeError(0.);
+ new_clust.setPosition(rClust);
+ new_clust.setPositionError({});
+
+ } // end 'merge_cluster(VecClust&)'
// --------------------------------------------------------------------------
- //! Get current center of protocluster
+ //! Collect associations of merged clusters
// --------------------------------------------------------------------------
- edm4hep::Vector3f TrackClusterMergeSplitter::get_cluster_position(const edm4eic::ProtoCluster& clust) const {
+ /*! This collects the associations of all clusters being merged, and
+ * and updates the weights of each to reflect the new merged cluster.
+ *
+ * As in `CalorimeterClusterRecoCoG`, an association is made for each
+ * contributing primary with a weight equal to the ratio of the
+ * contributed energy energy over total sim hit energy.
+ */
+ void TrackClusterMergeSplitter::collect_associations(
+ const edm4eic::MutableCluster& new_clust,
+ const VecClust& old_clusts,
+ const edm4eic::MCRecoClusterParticleAssociationCollection* old_clust_assocs,
+#if EDM4EIC_VERSION_MAJOR >= 7
+ const edm4eic::MCRecoCalorimeterHitAssociationCollection* old_hit_assocs,
+#else
+ const edm4hep::SimCalorimeterHitCollection* old_sim_hits,
+#endif
+ edm4eic::MCRecoClusterParticleAssociationCollection* new_clust_assocs
+ ) const {
- // grab total energy
- const float eClust = get_cluster_energy(clust) * m_cfg.sampFrac;
+ // lambda to compare MCParticles
+ auto compare = [](const edm4hep::MCParticle& lhs, const edm4hep::MCParticle& rhs) {
+ if (lhs.getObjectID().collectionID == rhs.getObjectID().collectionID) {
+ return (lhs.getObjectID().index < rhs.getObjectID().index);
+ } else {
+ return (lhs.getObjectID().collectionID < rhs.getObjectID().collectionID);
+ }
+ };
- // calculate energy-weighted center
- float wTotal = 0.;
- edm4hep::Vector3f position(0., 0., 0.);
- for (auto hit : clust.getHits()) {
+ // create a map to hold associated particles and
+ // their weight
+ std::map<edm4hep::MCParticle, double, decltype(compare)> mapMCParToWeight(compare);
- // calculate weight
- float weight = hit.getEnergy() / eClust;
- wTotal += weight;
+ // ------------------------------------------------------------------------
+ // loop through clusters
+ // ------------------------------------------------------------------------
+ double eMergeSimHitSum = 0.;
+ for (const auto& clust : old_clusts) {
- // update cluster position
- position = position + (hit.getPosition() * weight);
+ // ----------------------------------------------------------------------
+ // get total sim hit energy of old clust
+ // ----------------------------------------------------------------------
+ double eOldSimHitSum = 0.;
+ for (const auto& recHit : clust.getHits()) {
+#if EDM4EIC_VERSION_MAJOR >= 7
+ for (const auto& hitAssoc : *old_hit_assocs) {
+ if (hitAssoc.getRawHit() == recHit.getRawHit()) {
+ eOldSimHitSum += hitAssoc.getSimHit().getEnergy();
+ }
+ } // end hit association loop
+#else
+ for (const auto& simHit : *old_sim_hits) {
+ if (simHit.getCellID() == recHit.getCellID()) {
+ eOldSimHitSum += simHit.getEnergy();
+ break;
+ }
+ } // end sim hit loop
+#endif
+ } // end rec hit loop
+ eMergeSimHitSum += eOldSimHitSum;
+
+ // ----------------------------------------------------------------------
+ // now collect all associations for this old cluster
+ // ----------------------------------------------------------------------
+ for (const auto& clustAssoc : *old_clust_assocs) {
+ if (clustAssoc.getRec() == clust) {
+ mapMCParToWeight[ clustAssoc.getSim() ] += (clustAssoc.getWeight() * eOldSimHitSum);
+ }
+ } // end association loop
+ } // end cluster loop
+
+ // ------------------------------------------------------------------------
+ // scale weights by sum of sim hit energies
+ // ------------------------------------------------------------------------
+ double wAssocSum = 0.;
+ for (auto& [par, weight] : mapMCParToWeight) {
+ weight /= eMergeSimHitSum;
+ wAssocSum += weight;
}
- float norm = 1.;
- if (wTotal == 0.) {
- warning("Total weight of 0 in position calculation!");
- } else {
- norm = wTotal;
+ // normalize weights if not already
+ if (wAssocSum != 1.0) {
+ trace("Sum of weight not unity, normalizing by sum ({})", wAssocSum);
+ for (auto& [par, weight] : mapMCParToWeight) {
+ weight /= wAssocSum;
+ }
+ }
+
+ // ------------------------------------------------------------------------
+ // and finally create an association for each unique primary
+ // ------------------------------------------------------------------------
+ for (const auto& [par, weight] : mapMCParToWeight) {
+ auto assoc = new_clust_assocs->create();
+ assoc.setRecID( new_clust.getObjectID().index );
+ assoc.setSimID( par.getObjectID().index );
+ assoc.setWeight( weight );
+ assoc.setRec( new_clust );
+ assoc.setSim( par );
+ trace("Associated output cluster {} to MC Particle {} (pid = {}, status = {}, energy = {}) with weight ({})",
+ new_clust.getObjectID().index,
+ par.getObjectID().index,
+ par.getPDG(),
+ par.getGeneratorStatus(),
+ par.getEnergy(),
+ weight
+ );
}
- return position / norm;
- } // end 'get_cluster_position(edm4eic::ProtoCluster&)'
+ } // end 'collect_associations(...)'
} // end eicrecon namespace
diff --git a/src/algorithms/calorimetry/TrackClusterMergeSplitter.h b/src/algorithms/calorimetry/TrackClusterMergeSplitter.h
index e38e7e608..e3a86eabd 100644
--- a/src/algorithms/calorimetry/TrackClusterMergeSplitter.h
+++ b/src/algorithms/calorimetry/TrackClusterMergeSplitter.h
@@ -5,13 +5,24 @@
#include <DD4hep/Detector.h>
#include <algorithms/algorithm.h>
-#include <edm4eic/ProtoClusterCollection.h>
+#include <edm4eic/CalorimeterHit.h>
+#include <edm4eic/ClusterCollection.h>
+#include <edm4eic/EDM4eicVersion.h>
+#include <edm4eic/MCRecoClusterParticleAssociationCollection.h>
+#if EDM4EIC_VERSION_MAJOR >= 8
+#include <edm4eic/TrackClusterMatchCollection.h>
+#endif
+#include <edm4eic/Track.h>
#include <edm4eic/TrackPoint.h>
#include <edm4eic/TrackSegmentCollection.h>
-#include <edm4hep/Vector3f.h>
-#include <podio/ObjectID.h>
+#if EDM4EIC_VERSION_MAJOR >= 7
+#include <edm4eic/MCRecoCalorimeterHitAssociationCollection.h>
+#else
+#include <edm4hep/SimCalorimeterHitCollection.h>
+#endif
#include <algorithm>
#include <map>
+#include <optional>
#include <set>
#include <string>
#include <string_view>
@@ -26,14 +37,14 @@
namespace eicrecon {
// --------------------------------------------------------------------------
- //! Comparator struct for protoclusters
+ //! Comparator struct for object IDs
// --------------------------------------------------------------------------
- /*! Organizes protoclusters by their ObjectID's in decreasing collection
+ /*! Organizes objects by their ObjectID's in decreasing collection
* ID first, and second by decreasing index second.
*/
- struct CompareProto {
+ template <typename T> struct CompareObjectID {
- bool operator() (const edm4eic::ProtoCluster& lhs, const edm4eic::ProtoCluster& rhs) const {
+ bool operator() (const T& lhs, const T& rhs) const {
if (lhs.getObjectID().collectionID == rhs.getObjectID().collectionID) {
return (lhs.getObjectID().index < rhs.getObjectID().index);
} else {
@@ -43,16 +54,24 @@ namespace eicrecon {
}; // end CompareObjectID
+ // specialization for clusters, hits
+ using CompareClust = CompareObjectID<edm4eic::Cluster>;
+ using CompareHit = CompareObjectID<edm4eic::CalorimeterHit>;
+
// --------------------------------------------------------------------------
//! Convenience types
// --------------------------------------------------------------------------
+ using VecTrk = std::vector<edm4eic::Track>;
using VecProj = std::vector<edm4eic::TrackPoint>;
- using VecClust = std::vector<edm4eic::ProtoCluster>;
- using SetClust = std::set<edm4eic::ProtoCluster, CompareProto>;
- using MapToVecProj = std::map<edm4eic::ProtoCluster, VecProj, CompareProto>;
- using MapToVecClust = std::map<edm4eic::ProtoCluster, VecClust, CompareProto>;
+ using VecClust = std::vector<edm4eic::Cluster>;
+ using SetClust = std::set<edm4eic::Cluster, CompareClust>;
+ using MapToVecTrk = std::map<edm4eic::Cluster, VecTrk, CompareClust>;
+ using MapToVecProj = std::map<edm4eic::Cluster, VecProj, CompareClust>;
+ using MapToVecClust = std::map<edm4eic::Cluster, VecClust, CompareClust>;
+ using MapToWeight = std::map<edm4eic::CalorimeterHit, double, CompareHit>;
+ using VecWeights = std::vector<MapToWeight>;
@@ -61,11 +80,21 @@ namespace eicrecon {
// --------------------------------------------------------------------------
using TrackClusterMergeSplitterAlgorithm = algorithms::Algorithm<
algorithms::Input<
- edm4eic::ProtoClusterCollection,
- edm4eic::TrackSegmentCollection
+ edm4eic::ClusterCollection,
+ edm4eic::TrackSegmentCollection,
+ std::optional<edm4eic::MCRecoClusterParticleAssociationCollection>,
+#if EDM4EIC_VERSION_MAJOR >= 7
+ std::optional<edm4eic::MCRecoCalorimeterHitAssociationCollection>
+#else
+ std::optional<edm4hep::SimCalorimeterHitCollection>
+#endif
>,
algorithms::Output<
- edm4eic::ProtoClusterCollection
+ edm4eic::ClusterCollection,
+#if EDM4EIC_VERSION_MAJOR >= 8
+ edm4eic::TrackClusterMatchCollection,
+#endif
+ std::optional<edm4eic::MCRecoClusterParticleAssociationCollection>
>
>;
@@ -74,8 +103,8 @@ namespace eicrecon {
// --------------------------------------------------------------------------
//! Track-Based Cluster Merger/Splitter
// --------------------------------------------------------------------------
- /*! An algorithm which takes a collection of proto-clusters, matches
- * track projections, and then decides to merge or split those proto-
+ /*! An algorithm which takes a collection of clusters, matches
+ * track projections, and then decides to merge or split those
* clusters based on average E/p from simulations.
*
* Heavily inspired by Eur. Phys. J. C (2017) 77:466
@@ -91,8 +120,16 @@ namespace eicrecon {
TrackClusterMergeSplitter(std::string_view name) :
TrackClusterMergeSplitterAlgorithm {
name,
- {"InputProtoClusterCollection", "InputTrackProjections"},
- {"OutputProtoClusterCollection"},
+#if EDM4EIC_VERSION_MAJOR >= 7
+ {"InputClusterCollection", "InputTrackProjections", "InputMCClusterAssociations", "InputMCHitAssociations"},
+#else
+ {"InputClusterCollection", "InputTrackProjections", "InputMCClusterAssociations", "InputMCHits"},
+#endif
+#if EDM4EIC_VERSION_MAJOR >= 8
+ {"OutputProtoClusterCollection", "OutputTrackClusterMatches", "OutputMCClusterAssociations"},
+#else
+ {"OutputProtoClusterCollection", "OutputMCClusterAssociations"},
+#endif
"Merges or splits clusters based on tracks projected to them."
} {}
@@ -103,11 +140,15 @@ namespace eicrecon {
private:
// private methods
- void get_projections(const edm4eic::TrackSegmentCollection* projections, VecProj& relevant_projects) const;
- void match_clusters_to_tracks(const edm4eic::ProtoClusterCollection* clusters, const VecProj& projections, MapToVecProj& matches) const;
- void merge_and_split_clusters(const VecClust& to_merge, const VecProj& to_split, edm4eic::ProtoClusterCollection* out_protoclusters) const;
- float get_cluster_energy(const edm4eic::ProtoCluster& clust) const;
- edm4hep::Vector3f get_cluster_position(const edm4eic::ProtoCluster& clust) const;
+ void get_projections(const edm4eic::TrackSegmentCollection* projections, VecProj& relevant_projects, VecTrk& relevant_trks) const;
+ void match_clusters_to_tracks(const edm4eic::ClusterCollection* clusters, const VecProj& projections, const VecTrk& tracks, MapToVecProj& matched_projects, MapToVecTrk& matched_tracks) const;
+ void merge_and_split_clusters(const VecClust& to_merge, const VecProj& to_split, std::vector<edm4eic::MutableCluster>& new_clusters) const;
+ void make_cluster(const VecClust& old_clusts, edm4eic::MutableCluster& new_clust, std::optional<MapToWeight> split_weights = std::nullopt) const;
+#if EDM4EIC_VERSION_MAJOR >= 7
+ void collect_associations(const edm4eic::MutableCluster& new_clust, const VecClust& old_clusts, const edm4eic::MCRecoClusterParticleAssociationCollection* old_clust_assocs, const edm4eic::MCRecoCalorimeterHitAssociationCollection* old_hit_assocs, edm4eic::MCRecoClusterParticleAssociationCollection* new_clust_assocs) const;
+#else
+ void collect_associations(const edm4eic::MutableCluster& new_clust, const VecClust& old_clusts, const edm4eic::MCRecoClusterParticleAssociationCollection* old_clust_assocs, const edm4hep::SimCalorimeterHitCollection* old_sim_hits, edm4eic::MCRecoClusterParticleAssociationCollection* new_clust_assocs) const;
+#endif
// calorimeter id
int m_idCalo {0};
diff --git a/src/detectors/BHCAL/BHCAL.cc b/src/detectors/BHCAL/BHCAL.cc
index f205f4b89..13d7f17e7 100644
--- a/src/detectors/BHCAL/BHCAL.cc
+++ b/src/detectors/BHCAL/BHCAL.cc
@@ -194,10 +194,20 @@ extern "C" {
app->Add(
new JOmniFactoryGeneratorT<TrackClusterMergeSplitter_factory>(
- "HcalBarrelSplitMergeProtoClusters",
- {"HcalBarrelIslandProtoClusters",
- "CalorimeterTrackProjections"},
- {"HcalBarrelSplitMergeProtoClusters"},
+ "HcalBarrelSplitMergeClustersWithoutShapes",
+ {"HcalBarrelClusters",
+ "CalorimeterTrackProjections",
+ "HcalBarrelClusterAssociations",
+#if EDM4EIC_VERSION_MAJOR >= 7
+ "HcalBarrelRawHitAssociations"},
+#else
+ "HcalBarrelHits"},
+#endif
+ {"HcalBarrelSplitMergeClustersWithoutShapes",
+#if EDM4EIC_VERSION_MAJOR >= 8
+ "HcalBarrelTrackSplitMergeClusterMatches",
+#endif
+ "HcalBarrelSplitMergeClusterAssociationsWithoutShapes"},
{
.idCalo = "HcalBarrel_ID",
.minSigCut = -2.0,
@@ -211,23 +221,6 @@ extern "C" {
)
);
- app->Add(
- new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
- "HcalBarrelSplitMergeClustersWithoutShapes",
- {"HcalBarrelSplitMergeProtoClusters", // edm4eic::ProtoClusterCollection
- "HcalBarrelHits"}, // edm4hep::SimCalorimeterHitCollection
- {"HcalBarrelSplitMergeClustersWithoutShapes", // edm4eic::Cluster
- "HcalBarrelSplitMergeClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
- {
- .energyWeight = "log",
- .sampFrac = 1.0,
- .logWeightBase = 6.2,
- .enableEtaBounds = false
- },
- app // TODO: Remove me once fixed
- )
- );
-
app->Add(
new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
"HcalBarrelSplitMergeClusters",
@@ -235,10 +228,12 @@ extern "C" {
"HcalBarrelSplitMergeClusterAssociationsWithoutShapes"},
{"HcalBarrelSplitMergeClusters",
"HcalBarrelSplitMergeClusterAssociations"},
- {},
+ {
+ .energyWeight = "log",
+ .logWeightBase = 6.2
+ },
app
)
);
-
}
}
diff --git a/src/detectors/EEMC/EEMC.cc b/src/detectors/EEMC/EEMC.cc
index bac5fdd00..d86210964 100644
--- a/src/detectors/EEMC/EEMC.cc
+++ b/src/detectors/EEMC/EEMC.cc
@@ -185,25 +185,6 @@ extern "C" {
)
);
- app->Add(
- new JOmniFactoryGeneratorT<TrackClusterMergeSplitter_factory>(
- "EcalEndcapNSplitMergeProtoClusters",
- {"EcalEndcapNIslandProtoClusters",
- "CalorimeterTrackProjections"},
- {"EcalEndcapNSplitMergeProtoClusters"},
- {
- .idCalo = "EcalEndcapN_ID",
- .minSigCut = -1.0,
- .avgEP = 1.0,
- .sigEP = 0.10,
- .drAdd = 0.08,
- .sampFrac = 1.0,
- .transverseEnergyProfileScale = 1.0
- },
- app // TODO: remove me once fixed
- )
- );
-
#if EDM4EIC_VERSION_MAJOR >= 8
app->Add(new JOmniFactoryGeneratorT<CalorimeterParticleIDPreML_factory>(
"EcalEndcapNParticleIDPreML",
@@ -248,19 +229,31 @@ extern "C" {
#endif
app->Add(
- new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
- "EcalEndcapNSplitMergeClustersWithoutShapes",
- {"EcalEndcapNSplitMergeProtoClusters", // edm4eic::ProtoClusterCollection
- "EcalEndcapNHits"}, // edm4hep::SimCalorimeterHitCollection
- {"EcalEndcapNSplitMergeClustersWithoutShapes", // edm4eic::Cluster
- "EcalEndcapNSplitMergeClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
+ new JOmniFactoryGeneratorT<TrackClusterMergeSplitter_factory>(
+ "EcalEndcapNSplitMergeClustersWithoutShapes",
+ {"EcalEndcapNClustersWithoutPID",
+ "CalorimeterTrackProjections",
+ "EcalEndcapNClusterAssociations",
+#if EDM4EIC_VERSION_MAJOR >= 7
+ "EcalEndcapNRawHitAssociations"},
+#else
+ "EcalEndcapNHits"},
+#endif
+ {"EcalEndcapNSplitMergeClustersWithoutShapes",
+#if EDM4EIC_VERSION_MAJOR >= 8
+ "EcalEndcapNTrackSplitMergeClusterMatches",
+#endif
+ "EcalEndcapNSplitMergeClusterAssociationsWithoutShapes"},
{
- .energyWeight = "log",
+ .idCalo = "EcalEndcapN_ID",
+ .minSigCut = -1.0,
+ .avgEP = 1.0,
+ .sigEP = 0.10,
+ .drAdd = 0.08,
.sampFrac = 1.0,
- .logWeightBase = 3.6,
- .enableEtaBounds = false
+ .transverseEnergyProfileScale = 1.0
},
- app // TODO: Remove me once fixed
+ app // TODO: remove me once fixed
)
);
diff --git a/src/detectors/EHCAL/EHCAL.cc b/src/detectors/EHCAL/EHCAL.cc
index 59c2283ea..033e81652 100644
--- a/src/detectors/EHCAL/EHCAL.cc
+++ b/src/detectors/EHCAL/EHCAL.cc
@@ -163,10 +163,20 @@ extern "C" {
);
app->Add(
new JOmniFactoryGeneratorT<TrackClusterMergeSplitter_factory>(
- "HcalEndcapNSplitMergeProtoClusters",
- {"HcalEndcapNIslandProtoClusters",
- "CalorimeterTrackProjections"},
- {"HcalEndcapNSplitMergeProtoClusters"},
+ "HcalEndcapNSplitMergeClustersWithoutShapes",
+ {"HcalEndcapNClusters",
+ "CalorimeterTrackProjections",
+ "HcalEndcapNClusterAssociations",
+#if EDM4EIC_VERSION_MAJOR >= 7
+ "HcalEndcapNRawHitAssociations"},
+#else
+ "HcalEndcapNHits"},
+#endif
+ {"HcalEndcapNSplitMergeClustersWithoutShapes",
+#if EDM4EIC_VERSION_MAJOR >= 8
+ "HcalEndcapNTrackSplitMergeClusterMatches",
+#endif
+ "HcalEndcapNSplitMergeClusterAssociationsWithoutShapes"},
{
.idCalo = "HcalEndcapN_ID",
.minSigCut = -2.0,
@@ -179,22 +189,6 @@ extern "C" {
app // TODO: remove me once fixed
)
);
- app->Add(
- new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
- "HcalEndcapNSplitMergeClustersWithoutShapes",
- {"HcalEndcapNSplitMergeProtoClusters", // edm4eic::ProtoClusterCollection
- "HcalEndcapNHits"}, // edm4hep::SimCalorimeterHitCollection
- {"HcalEndcapNSplitMergeClustersWithoutShapes", // edm4eic::Cluster
- "HcalEndcapNSplitMergeClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
- {
- .energyWeight = "log",
- .sampFrac = 1.0,
- .logWeightBase = 6.2,
- .enableEtaBounds = false
- },
- app // TODO: Remove me once fixed
- )
- );
app->Add(
new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
"HcalEndcapNSplitMergeClusters",
diff --git a/src/detectors/FEMC/FEMC.cc b/src/detectors/FEMC/FEMC.cc
index 209a5fbd4..c2ecf44b3 100644
--- a/src/detectors/FEMC/FEMC.cc
+++ b/src/detectors/FEMC/FEMC.cc
@@ -160,10 +160,20 @@ extern "C" {
app->Add(
new JOmniFactoryGeneratorT<TrackClusterMergeSplitter_factory>(
- "EcalEndcapPSplitMergeProtoClusters",
- {"EcalEndcapPIslandProtoClusters",
- "CalorimeterTrackProjections"},
- {"EcalEndcapPSplitMergeProtoClusters"},
+ "EcalEndcapPSplitMergeClustersWithoutShapes",
+ {"EcalEndcapPClusters",
+ "CalorimeterTrackProjections",
+ "EcalEndcapPClusterAssociations",
+#if EDM4EIC_VERSION_MAJOR >= 7
+ "EcalEndcapPRawHitAssociations"},
+#else
+ "EcalEndcapPHits"},
+#endif
+ {"EcalEndcapPSplitMergeClustersWithoutShapes",
+#if EDM4EIC_VERSION_MAJOR >= 8
+ "EcalEndcapPTrackSplitMergeClusterMatches",
+#endif
+ "EcalEndcapPSplitMergeClusterAssociationsWithoutShapes"},
{
.idCalo = "EcalEndcapP_ID",
.minSigCut = -2.0,
@@ -177,23 +187,6 @@ extern "C" {
)
);
- app->Add(
- new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
- "EcalEndcapPSplitMergeClustersWithoutShapes",
- {"EcalEndcapPSplitMergeProtoClusters", // edm4eic::ProtoClusterCollection
- "EcalEndcapPHits"}, // edm4hep::SimCalorimeterHitCollection
- {"EcalEndcapPSplitMergeClustersWithoutShapes", // edm4eic::Cluster
- "EcalEndcapPSplitMergeClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
- {
- .energyWeight = "log",
- .sampFrac = 1.0,
- .logWeightBase = 3.6,
- .enableEtaBounds = false
- },
- app // TODO: Remove me once fixed
- )
- );
-
app->Add(
new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
"EcalEndcapPSplitMergeClusters",
diff --git a/src/detectors/FHCAL/FHCAL.cc b/src/detectors/FHCAL/FHCAL.cc
index a7673d976..d25565548 100644
--- a/src/detectors/FHCAL/FHCAL.cc
+++ b/src/detectors/FHCAL/FHCAL.cc
@@ -363,10 +363,20 @@ extern "C" {
app->Add(
new JOmniFactoryGeneratorT<TrackClusterMergeSplitter_factory>(
- "LFHCALSplitMergeProtoClusters",
- {"LFHCALIslandProtoClusters",
- "CalorimeterTrackProjections"},
- {"LFHCALSplitMergeProtoClusters"},
+ "LFHCALSplitMergeClustersWithoutShapes",
+ {"LFHCALClusters",
+ "CalorimeterTrackProjections",
+ "LFHCALClusterAssociations",
+#if EDM4EIC_VERSION_MAJOR >= 7
+ "LFHCALRawHitAssocitions"},
+#else
+ "LFHCALHit"},
+#endif
+ {"LFHCALSplitMergeClustersWithoutShapes",
+#if EDM4EIC_VERSION_MAJOR >= 8
+ "LFHCALTrackSplitMergeClusterMatches",
+#endif
+ "LFHCALSplitMergeClusterAssociationsWithoutShapes"},
{
.idCalo = "LFHCAL_ID",
.minSigCut = -2.0,
@@ -380,23 +390,6 @@ extern "C" {
)
);
- app->Add(
- new JOmniFactoryGeneratorT<CalorimeterClusterRecoCoG_factory>(
- "LFHCALSplitMergeClustersWithoutShapes",
- {"LFHCALSplitMergeProtoClusters", // edm4eic::ProtoClusterCollection
- "LFHCALHits"}, // edm4hep::SimCalorimeterHitCollection
- {"LFHCALSplitMergeClustersWithoutShapes", // edm4eic::Cluster
- "LFHCALSplitMergeClusterAssociationsWithoutShapes"}, // edm4eic::MCRecoClusterParticleAssociation
- {
- .energyWeight = "log",
- .sampFrac = 1.0,
- .logWeightBase = 4.5,
- .enableEtaBounds = false
- },
- app // TODO: Remove me once fixed
- )
- );
-
app->Add(
new JOmniFactoryGeneratorT<CalorimeterClusterShape_factory>(
"LFHCALSplitMergeClusters",
@@ -405,7 +398,9 @@ extern "C" {
{"LFHCALSplitMergeClusters",
"LFHCALSplitMergeClusterAssociations"},
{
- .longitudinalShowerInfoAvailable = true
+ .longitudinalShowerInfoAvailable = true,
+ .energyWeight = "log",
+ .logWeightBase = 4.5
},
app
)
diff --git a/src/factories/calorimetry/TrackClusterMergeSplitter_factory.h b/src/factories/calorimetry/TrackClusterMergeSplitter_factory.h
index e5d6daec9..fd710da14 100644
--- a/src/factories/calorimetry/TrackClusterMergeSplitter_factory.h
+++ b/src/factories/calorimetry/TrackClusterMergeSplitter_factory.h
@@ -7,6 +7,8 @@
#include <string>
// dd4hep utilities
#include <DD4hep/Detector.h>
+// edm utilities
+#include <edm4eic/EDM4eicVersion.h>
// eicrecon components
#include "extensions/jana/JOmniFactory.h"
#include "services/geometry/dd4hep/DD4hep_service.h"
@@ -27,11 +29,21 @@ namespace eicrecon {
std::unique_ptr<AlgoT> m_algo;
// input collections
- PodioInput<edm4eic::ProtoCluster> m_protoclusters_input {this};
+ PodioInput<edm4eic::Cluster> m_clusters_input {this};
PodioInput<edm4eic::TrackSegment> m_track_projections_input {this};
+ PodioInput<edm4eic::MCRecoClusterParticleAssociation> m_cluster_association_input {this};
+#if EDM4EIC_VERSION_MAJOR >= 7
+ PodioInput<edm4eic::MCRecoCalorimeterHitAssociation> m_hit_association_input {this};
+#else
+ PodioInput<edm4hep::SimCalorimeterHit> m_sim_hit_input {this};
+#endif
// output collections
- PodioOutput<edm4eic::ProtoCluster> m_protoclusters_output {this};
+ PodioOutput<edm4eic::Cluster> m_clusters_output {this};
+#if EDM4EIC_VERSION_MAJOR >= 8
+ PodioOutput<edm4eic::TrackClusterMatch> m_track_cluster_match_output {this};
+#endif
+ PodioOutput<edm4eic::MCRecoClusterParticleAssociation> m_cluster_association_output {this};
// parameter bindings
ParameterRef<std::string> m_idCalo {this, "idCalo", config().idCalo};
@@ -60,8 +72,16 @@ namespace eicrecon {
void Process(int64_t run_number, uint64_t event_number) {
m_algo->process(
- {m_protoclusters_input(), m_track_projections_input()},
- {m_protoclusters_output().get()}
+#if EDM4EIC_VERSION_MAJOR >= 7
+ {m_clusters_input(), m_track_projections_input(), m_cluster_association_input(), m_hit_association_input()},
+#else
+ {m_clusters_input(), m_track_projections_input(), m_cluster_association_input(), m_sim_hit_input()},
+#endif
+#if EDM4EIC_VERSION_MAJOR >= 8
+ {m_clusters_output().get(), m_track_cluster_match_output().get(), m_cluster_association_output().get()}
+#else
+ {m_clusters_output().get(), m_cluster_association_output().get()}
+#endif
);
}