algorithms/reco/ClustersToParticles.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2026 ePIC Collaboration
0003
0004 #include <edm4eic/EDM4eicVersion.h>
0005 #include <edm4hep/MCParticleCollection.h>
0006 #include <edm4hep/Vector3f.h>
0007 #include <podio/ObjectID.h>
0008 #include <podio/detail/Link.h>
0009 #include <podio/detail/LinkCollectionImpl.h>
0010 #include <cmath>
0011 #include <memory>
0012 #include <tuple>
0013
0014 #include "ClustersToParticles.h"
0015
0016 namespace eicrecon {
0017
0018 void ClustersToParticles::init() {
0019   const auto& particle = m_particleSvc.particle(m_cfg.pdgCode);
0020   m_mass               = particle.mass;
0021   m_charge             = particle.charge;
0022 }
0023
0024 void ClustersToParticles::process(const ClustersToParticles::Input& input,
0025                                   const ClustersToParticles::Output& output) const {
0026   const auto [clusters, cluster_assocs] = input;
0027 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0028   auto [parts, part_links, part_assocs] = output;
0029 #else
0030   auto [parts, part_assocs] = output;
0031 #endif
0032
0033   for (const auto& cluster : *clusters) {
0034     const auto energy = cluster.getEnergy();
0035     const auto pos    = cluster.getPosition();
0036
0037     const auto momentum_mag =
0038         (energy > m_mass) ? std::sqrt(energy * energy - m_mass * m_mass) : 0.0;
0039     const auto pos_mag = std::sqrt(pos.x * pos.x + pos.y * pos.y + pos.z * pos.z);
0040
0041     edm4hep::Vector3f momentum{0, 0, 0};
0042     if (pos_mag > 0) {
0043       momentum.x = momentum_mag * pos.x / pos_mag;
0044       momentum.y = momentum_mag * pos.y / pos_mag;
0045       momentum.z = momentum_mag * pos.z / pos_mag;
0046     }
0047
0048     debug("Converting cluster: index={:<4} energy={:<8.3f} position=({:<8.3f}, {:<8.3f}, {:<8.3f})",
0049           cluster.getObjectID().index, energy, pos.x, pos.y, pos.z);
0050
0051     auto rec_part = parts->create();
0052     rec_part.addToClusters(cluster);
0053     rec_part.setPDG(m_cfg.pdgCode);
0054     rec_part.setType(0);
0055     rec_part.setEnergy(energy);
0056     rec_part.setMomentum(momentum);
0057     rec_part.setCharge(static_cast<float>(m_charge));
0058     rec_part.setMass(m_mass);
0059     rec_part.setGoodnessOfPID(0); // assume no PID until proven otherwise
0060     rec_part.setReferencePoint({0.0f, 0.0f, 0.0f});
0061
0062     // Handle associations if provided
0063     for (auto cluster_assoc : *cluster_assocs) {
0064       if (cluster_assoc.getRec() == cluster) {
0065         trace("Found cluster association: index={} -> index={}, weight={}",
0066               cluster_assoc.getRec().getObjectID().index,
0067               cluster_assoc.getSim().getObjectID().index, cluster_assoc.getWeight());
0068 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0069         auto part_link = part_links->create();
0070         part_link.setFrom(rec_part);
0071         part_link.setTo(cluster_assoc.getSim());
0072         part_link.setWeight(cluster_assoc.getWeight());
0073 #endif
0074         auto part_assoc = part_assocs->create();
0075         part_assoc.setRec(rec_part);
0076         part_assoc.setSim(cluster_assoc.getSim());
0077         part_assoc.setWeight(cluster_assoc.getWeight());
0078       }
0079     }
0080   }
0081 }
0082
0083 } // namespace eicrecon