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 <fmt/core.h>
0006 #include <podio/RelationRange.h>
0007 #include <stdint.h>
0008 #include <gsl/pointers>
0009 #include <optional>
0010 #include <set>
0011 #include <vector>
0012 #if EDM4EIC_VERSION_MAJOR >= 8
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, const TrackClusterMatch::Output& output) const {
0025         auto [tracks, clusters] = input;
0026         auto [matched_particles] = output;
0027         trace("We have {} tracks and {} clusters", tracks->size(), clusters->size());
0028
0029         std::set<int> used_tracks;
0030         // Loop across each cluster, and find the cloeset projected track
0031         for (auto cluster : *clusters) {
0032             const double cluster_eta = edm4hep::utils::eta(cluster.getPosition());
0033             const double cluster_phi = edm4hep::utils::angleAzimuthal(cluster.getPosition());
0034             trace("Cluster at eta={}, phi={}", cluster_eta, cluster_phi);
0035             // TODO: Get the detector ID so I can check if a projection is in the appropriate calorimeter
0036
0037             std::optional<edm4eic::TrackSegment> closest_segment;
0038             int closest_segment_id = -1;
0039             double closest_distance = m_cfg.matching_distance;
0040             // Loop through each track segment, and its points
0041             for (int closest_id = 0; auto track : *tracks) {
0042                 if (used_tracks.count(closest_id) > 0) {
0043                     trace("Skipping track segment already used");
0044                     continue;
0045                 }
0046                 for (auto point : track.getPoints()) {
0047                     // Check if the point is at the calorimeter
0048                     // int id = m_detector->volumeManager().lookupDetector(cluster.getHits()[0].getCellID()).id(); // TODO: Find programmatic way to get detector cluster is from
0049                     uint32_t ecal_barrel_id = m_geo.detector()->constant<int>("EcalBarrel_ID");
0050                     uint32_t hcal_barrel_id = m_geo.detector()->constant<int>("HcalBarrel_ID");
0051                     bool is_ecal = point.system == ecal_barrel_id;
0052                     bool is_hcal = point.system == hcal_barrel_id;
0053                     bool is_surface = point.surface == 1;
0054
0055                     if (!(is_ecal || is_hcal) || !is_surface) {
0056                         trace("Skipping track point not at the calorimeter");
0057                         continue;
0058                     }
0059                     const double track_eta = edm4hep::utils::eta(point.position);
0060                     const double track_phi = edm4hep::utils::angleAzimuthal(point.position);
0061                     trace("Track point at eta={}, phi={}", track_eta, track_phi);
0062
0063                     double delta = distance(cluster.getPosition(), point.position);
0064                     trace("Distance between cluster and track point: {}", delta);
0065                     if (delta < closest_distance) {
0066                         closest_distance = delta;
0067                         closest_segment = track;
0068                         closest_segment_id = closest_id;
0069                     }
0070                 }
0071                 closest_id++;
0072             }
0073             // If we found a point, create a new particle
0074             if (closest_segment) {
0075                 trace("Found a closest point at distance {}", closest_distance);
0076                 auto particle = matched_particles->create();
0077                 particle.setCluster(cluster);
0078                 particle.setTrack(closest_segment.value().getTrack());
0079                 used_tracks.insert(closest_segment_id);
0080             }
0081         }
0082         trace("Matched {} particles", matched_particles->size());
0083     }
0084
0085     double TrackClusterMatch::distance(const edm4hep::Vector3f& v1, const edm4hep::Vector3f& v2) const {
0086         double cluster_eta = edm4hep::utils::eta(v1);
0087         double cluster_phi = edm4hep::utils::angleAzimuthal(v1);
0088         double track_eta = edm4hep::utils::eta(v2);
0089         double track_phi = edm4hep::utils::angleAzimuthal(v2);
0090         double deta = cluster_eta - track_eta;
0091         double dphi = Phi_mpi_pi(cluster_phi - track_phi);
0092         return std::hypot(deta, dphi);
0093     }
0094 }
0095
0096 #endif // EDM4EIC_VERSION_MAJOR >= 8