|
17 | 17 | #include <podio/ObjectID.h> |
18 | 18 | #include <cmath> |
19 | 19 | #include <gsl/pointers> |
| 20 | +#include <iterator> |
20 | 21 | #include <map> |
| 22 | +#include <utility> |
21 | 23 |
|
22 | 24 | #include "MatchClusters.h" |
23 | 25 |
|
@@ -119,43 +121,60 @@ void MatchClusters::process(const MatchClusters::Input& input, |
119 | 121 | } |
120 | 122 |
|
121 | 123 | // get a map of mcID --> cluster |
122 | | -// input: clusters --> all clusters |
| 124 | +// For each cluster, pick the best associated MC particle by association weight. |
| 125 | +// Returns a map keyed by mcID and valued with the selected cluster. |
123 | 126 | std::map<int, edm4eic::Cluster> MatchClusters::indexedClusters( |
124 | 127 | const edm4eic::ClusterCollection* clusters, |
125 | 128 | const edm4eic::MCRecoClusterParticleAssociationCollection* associations) const { |
126 | 129 |
|
127 | | - std::map<int, edm4eic::Cluster> matched = {}; |
| 130 | + // temporary map: mcID -> (cluster, weight) so we can choose the cluster with highest weight per mcID |
| 131 | + std::map<int, std::pair<edm4eic::Cluster, float>> bestForMc; |
128 | 132 |
|
129 | | - // loop over clusters |
| 133 | + // loop over clusters and pick their best MC association by weight |
130 | 134 | for (const auto cluster : *clusters) { |
131 | 135 |
|
132 | | - int mcID = -1; |
| 136 | + int bestMcID = -1; |
| 137 | + float bestWeight = -1.F; |
133 | 138 |
|
134 | | - // find associated particle |
| 139 | + // find best associated MC particle for this cluster (largest association weight) |
135 | 140 | for (const auto assoc : *associations) { |
136 | 141 | if (assoc.getRec() == cluster) { |
137 | | - mcID = assoc.getSim().getObjectID().index; |
138 | | - break; |
| 142 | + const int candMcID = assoc.getSim().getObjectID().index; |
| 143 | + const float w = assoc.getWeight(); |
| 144 | + if (w > bestWeight) { |
| 145 | + bestWeight = w; |
| 146 | + bestMcID = candMcID; |
| 147 | + } |
139 | 148 | } |
140 | 149 | } |
141 | 150 |
|
142 | | - trace(" --> Found cluster with mcID {} and energy {}", mcID, cluster.getEnergy()); |
| 151 | + trace(" --> Found cluster with best mcID {} weight {} and energy {}", bestMcID, bestWeight, |
| 152 | + cluster.getEnergy()); |
143 | 153 |
|
144 | | - if (mcID < 0) { |
| 154 | + if (bestMcID < 0) { |
145 | 155 | trace(" --> WARNING: no valid MC truth link found, skipping cluster..."); |
146 | 156 | continue; |
147 | 157 | } |
148 | 158 |
|
149 | | - const bool duplicate = matched.contains(mcID); |
150 | | - if (duplicate) { |
151 | | - trace(" --> WARNING: this is a duplicate mcID, keeping the higher energy cluster"); |
152 | | - |
153 | | - if (cluster.getEnergy() < matched[mcID].getEnergy()) { |
154 | | - continue; |
| 159 | + // For this mcID, keep the cluster with the highest association weight (tie-break by energy). |
| 160 | + auto it = bestForMc.find(bestMcID); |
| 161 | + if (it == bestForMc.end()) { |
| 162 | + bestForMc.emplace(bestMcID, std::make_pair(cluster, bestWeight)); |
| 163 | + } else { |
| 164 | + const float existingWeight = it->second.second; |
| 165 | + if (bestWeight > existingWeight || |
| 166 | + (bestWeight == existingWeight && cluster.getEnergy() > it->second.first.getEnergy())) { |
| 167 | + it->second = std::make_pair(cluster, bestWeight); |
155 | 168 | } |
156 | 169 | } |
157 | | - matched[mcID] = cluster; |
158 | 170 | } |
| 171 | + |
| 172 | + // Convert to the old API: map<int, edm4eic::Cluster> |
| 173 | + std::map<int, edm4eic::Cluster> matched; |
| 174 | + for (const auto& kv : bestForMc) { |
| 175 | + matched.emplace(kv.first, kv.second.first); |
| 176 | + } |
| 177 | + |
159 | 178 | return matched; |
160 | 179 | } |
161 | 180 |
|
|
0 commit comments