algorithms/reco/TrackClusterMatch.cc

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