File indexing completed on 2025-02-24 10:15:50
0001
0002
0003
0004 #pragma once
0005
0006 #include <limits>
0007
0008 #include <algorithms/algorithm.h>
0009 #include <fmt/format.h>
0010 #include <spdlog/spdlog.h>
0011 #include <edm4eic/ClusterCollection.h>
0012 #include <edm4eic/MCRecoClusterParticleAssociationCollection.h>
0013 #include <edm4hep/utils/vector_utils.h>
0014
0015 #include "algorithms/interfaces/WithPodConfig.h"
0016 #include "EnergyPositionClusterMergerConfig.h"
0017
0018 namespace eicrecon {
0019
0020 using EnergyPositionClusterMergerAlgorithm = algorithms::Algorithm<
0021 algorithms::Input<
0022 edm4eic::ClusterCollection,
0023 edm4eic::MCRecoClusterParticleAssociationCollection,
0024 edm4eic::ClusterCollection,
0025 edm4eic::MCRecoClusterParticleAssociationCollection
0026 >,
0027 algorithms::Output<
0028 edm4eic::ClusterCollection,
0029 edm4eic::MCRecoClusterParticleAssociationCollection
0030 >
0031 >;
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045 class EnergyPositionClusterMerger
0046 : public EnergyPositionClusterMergerAlgorithm,
0047 public WithPodConfig<EnergyPositionClusterMergerConfig> {
0048
0049 public:
0050 EnergyPositionClusterMerger(std::string_view name)
0051 : EnergyPositionClusterMergerAlgorithm{name,
0052 {"energyClusterCollection", "energyClusterAssociations",
0053 "positionClusterCollection", "positionClusterAssociations"},
0054 {"outputClusterCollection", "outputClusterAssociations"},
0055 "Merge energy and position clusters if matching."} {}
0056
0057 public:
0058
0059 void init() { }
0060
0061 void process(const Input& input, const Output& output) const final {
0062
0063 const auto [energy_clus, energy_assoc, pos_clus, pos_assoc] = input;
0064 auto [merged_clus, merged_assoc] = output;
0065
0066 debug( "Merging energy and position clusters for new event" );
0067
0068 if (energy_clus->size() == 0 && pos_clus->size() == 0) {
0069 debug( "Nothing to do for this event, returning..." );
0070 return;
0071 }
0072
0073 std::vector<bool> consumed(energy_clus->size(), false);
0074
0075
0076 for (const auto& pc : *pos_clus) {
0077
0078 trace(" --> Processing position cluster {}, energy: {}", pc.getObjectID().index, pc.getEnergy());
0079
0080
0081 int best_match = -1;
0082 double best_delta = std::numeric_limits<double>::max();
0083 for (size_t ie = 0; ie < energy_clus->size(); ++ie) {
0084 if (consumed[ie]) {
0085 continue;
0086 }
0087
0088 const auto& ec = (*energy_clus)[ie];
0089
0090 trace(" --> Evaluating energy cluster {}, energy: {}", ec.getObjectID().index, ec.getEnergy());
0091
0092
0093
0094 const double de_rel = std::abs((pc.getEnergy() - ec.getEnergy()) / ec.getEnergy());
0095 const double deta = std::abs(edm4hep::utils::eta(pc.getPosition())
0096 - edm4hep::utils::eta(ec.getPosition()));
0097
0098
0099 const double dphi = edm4hep::utils::angleAzimuthal(pc.getPosition())
0100 - edm4hep::utils::angleAzimuthal(ec.getPosition());
0101 const double dsphi = std::abs(sin(0.5 * dphi));
0102 if ((m_cfg.energyRelTolerance > 0 && de_rel > m_cfg.energyRelTolerance) ||
0103 (m_cfg.etaTolerance > 0 && deta > m_cfg.etaTolerance) ||
0104 (m_cfg.phiTolerance > 0 && dsphi > sin(0.5 * m_cfg.phiTolerance))) {
0105 continue;
0106 }
0107
0108
0109
0110
0111 const double delta = fabs(pc.getEnergy() - ec.getEnergy());
0112 if (delta < best_delta) {
0113 best_delta = delta;
0114 best_match = ie;
0115 }
0116 }
0117
0118
0119 if (best_match >= 0) {
0120
0121 const auto& ec = (*energy_clus)[best_match];
0122
0123 auto new_clus = merged_clus->create();
0124 new_clus.setEnergy(ec.getEnergy());
0125 new_clus.setEnergyError(ec.getEnergyError());
0126 new_clus.setTime(pc.getTime());
0127 new_clus.setNhits(pc.getNhits() + ec.getNhits());
0128 new_clus.setPosition(pc.getPosition());
0129 new_clus.setPositionError(pc.getPositionError());
0130 new_clus.addToClusters(pc);
0131 new_clus.addToClusters(ec);
0132
0133 trace(" --> Found matching energy cluster {}, energy: {}", ec.getObjectID().index, ec.getEnergy() );
0134 trace(" --> Created a new combined cluster {}, energy: {}", new_clus.getObjectID().index, new_clus.getEnergy() );
0135
0136
0137 auto ea = energy_assoc->begin();
0138 for (; ea != energy_assoc->end(); ++ea) {
0139 if (ea->getRec() == ec) {
0140 break;
0141 }
0142 }
0143
0144 auto pa = pos_assoc->begin();
0145 for (; pa != pos_assoc->end(); ++pa) {
0146 if (pa->getRec() == pc) {
0147 break;
0148 }
0149 }
0150 if (ea != energy_assoc->end() || pa != pos_assoc->end()) {
0151
0152 if (ea != energy_assoc->end() && pa != pos_assoc->end()) {
0153
0154 if (pa->getSimID() == ea->getSimID()) {
0155
0156 auto clusterassoc = merged_assoc->create();
0157 clusterassoc.setRecID(new_clus.getObjectID().index);
0158 clusterassoc.setSimID(ea->getSimID());
0159 clusterassoc.setWeight(1.0);
0160 clusterassoc.setRec(new_clus);
0161 clusterassoc.setSim(ea->getSim());
0162 } else {
0163
0164 debug(" --> Two associations added to {} and {}", ea->getSimID(), pa->getSimID());
0165 auto clusterassoc1 = merged_assoc->create();
0166 clusterassoc1.setRecID(new_clus.getObjectID().index);
0167 clusterassoc1.setSimID(ea->getSimID());
0168 clusterassoc1.setWeight(0.5);
0169 clusterassoc1.setRec(new_clus);
0170 clusterassoc1.setSim(ea->getSim());
0171 auto clusterassoc2 = merged_assoc->create();
0172 clusterassoc2.setRecID(new_clus.getObjectID().index);
0173 clusterassoc2.setSimID(pa->getSimID());
0174 clusterassoc2.setWeight(0.5);
0175 clusterassoc2.setRec(new_clus);
0176 clusterassoc2.setSim(pa->getSim());
0177 }
0178 } else if (ea != energy_assoc->end()) {
0179
0180 debug(" --> Only added energy cluster association to {}", ea->getSimID());
0181 auto clusterassoc = merged_assoc->create();
0182 clusterassoc.setRecID(new_clus.getObjectID().index);
0183 clusterassoc.setSimID(ea->getSimID());
0184 clusterassoc.setWeight(1.0);
0185 clusterassoc.setRec(new_clus);
0186 clusterassoc.setSim(ea->getSim());
0187 } else if (pa != pos_assoc->end()) {
0188
0189 debug(" --> Only added position cluster association to {}", pa->getSimID());
0190 auto clusterassoc = merged_assoc->create();
0191 clusterassoc.setRecID(new_clus.getObjectID().index);
0192 clusterassoc.setSimID(pa->getSimID());
0193 clusterassoc.setWeight(1.0);
0194 clusterassoc.setRec(new_clus);
0195 clusterassoc.setSim(pa->getSim());
0196 }
0197 }
0198
0199
0200 consumed[best_match] = true;
0201
0202 debug(" Matched position cluster {} with energy cluster {}", pc.getObjectID().index, ec.getObjectID().index);
0203 debug(" - Position cluster: (E: {}, phi: {}, z: {})", pc.getEnergy(),
0204 edm4hep::utils::angleAzimuthal(pc.getPosition()), pc.getPosition().z);
0205 debug(" - Energy cluster: (E: {}, phi: {}, z: {})", ec.getEnergy(),
0206 edm4hep::utils::angleAzimuthal(ec.getPosition()), ec.getPosition().z);
0207 debug(" ---> Merged cluster: (E: {}, phi: {}, z: {})", new_clus.getEnergy(),
0208 edm4hep::utils::angleAzimuthal(new_clus.getPosition()), new_clus.getPosition().z);
0209
0210
0211 } else {
0212
0213 debug(" Unmatched position cluster {}", pc.getObjectID().index);
0214
0215 }
0216
0217 }
0218 }
0219 };
0220
0221 }