algorithms/reco/TrackClusterMatch.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025 Tristan Protzman
0003
0004 #include <edm4eic/EDM4eicVersion.h> // Needs edm4eic::TrackClusterMatch
0005 #include <edm4eic/Track.h>
0006 #include <fmt/core.h>
0007 #include <podio/RelationRange.h>
0008 #include <cstdint>
0009 #include <gsl/pointers>
0010 #include <optional>
0011 #include <set>
0012 #include <stdexcept>
0013 #include <vector>
0014 #if EDM4EIC_VERSION_MAJOR >= 8
0015
0016 #include <DD4hep/Detector.h>
0017 #include <edm4eic/ClusterCollection.h>
0018 #include <edm4eic/TrackPoint.h>
0019 #include <edm4hep/Vector3f.h>
0020 #include <edm4hep/utils/vector_utils.h>
0021
0022 #include "algorithms/reco/TrackClusterMatch.h"
0023 #include "algorithms/reco/TrackClusterMatchConfig.h"
0024
0025 namespace eicrecon {
0026 void TrackClusterMatch::process(const TrackClusterMatch::Input& input,
0027                                 const TrackClusterMatch::Output& output) const {
0028   auto [tracks, clusters]  = input;
0029   auto [matched_particles] = output;
0030   trace("We have {} tracks and {} clusters", tracks->size(), clusters->size());
0031
0032   // Validate the configuration
0033   if (m_cfg.matching_distance <= 0) {
0034     throw std::runtime_error(fmt::format("Invalid matching distance: {}", m_cfg.matching_distance));
0035   }
0036   if (m_cfg.calo_id.empty()) {
0037     throw std::runtime_error("Calorimeter ID must be set in the configuration");
0038   }
0039
0040   std::set<int> used_tracks;
0041   // Loop across each cluster, and find the cloeset projected track
0042   for (auto cluster : *clusters) {
0043     const double cluster_eta = edm4hep::utils::eta(cluster.getPosition());
0044     const double cluster_phi = edm4hep::utils::angleAzimuthal(cluster.getPosition());
0045     trace("Cluster at eta={}, phi={}", cluster_eta, cluster_phi);
0046     // TODO: Get the detector ID so I can check if a projection is in the appropriate calorimeter
0047
0048     std::optional<edm4eic::TrackSegment> closest_segment;
0049     int closest_segment_id  = -1;
0050     double closest_distance = m_cfg.matching_distance;
0051     // Loop through each track segment, and its points
0052     for (int closest_id = 0; auto track : *tracks) {
0053       if (used_tracks.contains(closest_id)) {
0054         trace("Skipping track segment already used");
0055         continue;
0056       }
0057       for (auto point : track.getPoints()) {
0058         // Check if the point is at the calorimeter
0059         // int id = m_detector->volumeManager().lookupDetector(cluster.getHits()[0].getCellID()).id(); // TODO: Find programmatic way to get detector cluster is from
0060         uint32_t calo_id = m_geo.detector()->constant<int>(m_cfg.calo_id);
0061         bool is_calo     = point.system == calo_id;
0062         bool is_surface  = point.surface == 1;
0063
0064         if (!is_calo || !is_surface) {
0065           trace("Skipping track point not at the calorimeter");
0066           continue;
0067         }
0068         const double track_eta = edm4hep::utils::eta(point.position);
0069         const double track_phi = edm4hep::utils::angleAzimuthal(point.position);
0070         trace("Track point at eta={}, phi={}", track_eta, track_phi);
0071
0072         double delta = distance(cluster.getPosition(), point.position);
0073         trace("Distance between cluster and track point: {}", delta);
0074         if (delta < closest_distance) {
0075           closest_distance   = delta;
0076           closest_segment    = track;
0077           closest_segment_id = closest_id;
0078         }
0079       }
0080       closest_id++;
0081     }
0082     // If we found a point, create a new particle
0083     if (closest_segment) {
0084       trace("Found a closest point at distance {}", closest_distance);
0085       auto particle = matched_particles->create();
0086       particle.setCluster(cluster);
0087       particle.setTrack(closest_segment.value().getTrack());
0088       used_tracks.insert(closest_segment_id);
0089     }
0090   }
0091   trace("Matched {} particles", matched_particles->size());
0092 }
0093
0094 double TrackClusterMatch::distance(const edm4hep::Vector3f& v1, const edm4hep::Vector3f& v2) {
0095   double cluster_eta = edm4hep::utils::eta(v1);
0096   double cluster_phi = edm4hep::utils::angleAzimuthal(v1);
0097   double track_eta   = edm4hep::utils::eta(v2);
0098   double track_phi   = edm4hep::utils::angleAzimuthal(v2);
0099   double deta        = cluster_eta - track_eta;
0100   double dphi        = Phi_mpi_pi(cluster_phi - track_phi);
0101   return std::hypot(deta, dphi);
0102 }
0103 } // namespace eicrecon
0104
0105 #endif // EDM4EIC_VERSION_MAJOR >= 8