algorithms/reco/TrackClusterMatch.cc

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