125 {
128
129
133 <<
"' = " << segments->
size());
134
135 if (truthDiag) {
136 std::size_t truthSegs = 0, bkgSegs = 0;
139 }
142 }
143
144 SegmentEdgeGraph graph{};
145 std::vector<SegmentEdgeScore>
scores{};
148 << ": built graph with nodes=" << graph.nNodes
149 << ", edges=" << graph.nEdges);
150
151
152
153 std::vector<int32_t> nodeTruthId;
154 if (truthDiag) {
155 nodeTruthId.assign(graph.nNodes, -1);
156 std::size_t truthSegs = 0, bkgSegs = 0;
157 for (std::size_t node = 0; node < graph.nNodes; ++node) {
160 if (truthPart) {
161 nodeTruthId[node] = static_cast<int32_t>(truthPart->index());
162 ++truthSegs;
163 } else {
164 ++bkgSegs;
165 }
166 }
169 }
170
173 float minProb = std::numeric_limits<float>::max();
174 float maxProb = std::numeric_limits<float>::lowest();
175 for (
const SegmentEdgeScore& score :
scores) {
176 minProb = std::min(minProb,
score.probability);
177 maxProb = std::max(maxProb,
score.probability);
178 }
180 <<
": edge scores=" <<
scores.size()
181 << ", prob range=[" << minProb << ", " << maxProb << "]");
182 } else {
184 << ": no edge scores produced");
185 }
186
187
188
189 std::unordered_map<std::uint64_t, float> pairProbability;
190 pairProbability.reserve(
scores.size());
191 for (
const SegmentEdgeScore& score :
scores) {
192 if (
score.src >= graph.nNodes ||
score.dst >= graph.nNodes ||
194 continue;
195 }
196 const std::uint64_t
key = undirectedPairKey(
score.src,
score.dst);
197 auto [
it, inserted] = pairProbability.emplace(key,
score.probability);
198 if (!inserted)
it->second = std::max(
it->second,
score.probability);
199 }
200
201
202
203 if (truthDiag) {
204 std::size_t trueTotal = 0, truePassed = 0, bkgTotal = 0, bkgPassed = 0;
205 for (const auto& [key, probability] : pairProbability) {
207 if (first >= graph.nNodes || second >= graph.nNodes) continue;
208 const bool isTrueEdge =
211 if (isTrueEdge) {
212 ++trueTotal;
214 } else {
215 ++bkgTotal;
217 }
218 }
223 }
224
225
226
227 using WeightedEdge = std::pair<std::uint64_t, float>;
228 const auto betterWeightedEdge = [](
const WeightedEdge&
first,
229 const WeightedEdge&
second) {
232 if (probabilityOrder != 0) {
233 return probabilityOrder < 0;
234 }
236 };
237 std::vector<std::vector<WeightedEdge>> edgesByNode(graph.nNodes);
238 std::size_t thresholdPairs = 0;
239 for (const auto& [key, probability] : pairProbability) {
242 if (first >= graph.nNodes || second >= graph.nNodes) continue;
243 edgesByNode[
first].emplace_back(key, probability);
244 edgesByNode[
second].emplace_back(key, probability);
245 ++thresholdPairs;
246 }
247
248 std::size_t thresholdedNodes = 0;
249 for (const std::vector<WeightedEdge>& nodeEdges : edgesByNode) {
250 thresholdedNodes += !nodeEdges.empty();
251 }
252
253 std::unordered_map<std::uint64_t, unsigned char> nominations;
254 nominations.reserve(thresholdPairs);
255 for (std::vector<WeightedEdge>& nodeEdges : edgesByNode) {
256 std::sort(nodeEdges.begin(), nodeEdges.end(), betterWeightedEdge);
260 }
261 for (const WeightedEdge& edge : nodeEdges) {
262 ++nominations[edge.first];
263 }
264 }
265
266 std::size_t mutualTopKPairs = 0;
267 std::size_t oneSidedTopKPairs = 0;
268 for (
const auto& [_,
count] : nominations) {
270 ++mutualTopKPairs;
271 } else {
272 ++oneSidedTopKPairs;
273 }
274 }
275
276 std::unordered_set<std::uint64_t> selectedPairKeys;
277 selectedPairKeys.reserve(thresholdPairs);
279
280 std::vector<WeightedEdge> acceptedPairs;
281 acceptedPairs.reserve(thresholdPairs);
282 for (const auto& [key, probability] : pairProbability) {
284 acceptedPairs.emplace_back(key, probability);
285 }
286 std::sort(acceptedPairs.begin(), acceptedPairs.end(), betterWeightedEdge);
287
289 std::vector<unsigned int>
degree(graph.nNodes, 0);
290 for (const WeightedEdge& edge : acceptedPairs) {
291 const auto [
first,
second] = unpackPairKey(edge.first);
292 if (maxDegree != 0 &&
293 (degree[first] >= maxDegree || degree[second] >= maxDegree)) {
294 continue;
295 }
296 selectedPairKeys.insert(edge.first);
299 }
300 } else {
301 for (
const auto& [key,
count] : nominations) {
303 selectedPairKeys.insert(key);
304 }
305 }
306
307
308
309 std::size_t orphanRecoveryPairs = 0;
313 std::vector<unsigned char> selectedNode(graph.nNodes, 0);
314 for (const std::uint64_t key : selectedPairKeys) {
316 if (first < graph.nNodes) selectedNode[
first] = 1;
317 if (second < graph.nNodes) selectedNode[
second] = 1;
318 }
319 for (std::size_t node = 0; node < graph.nNodes; ++node) {
320 if (selectedNode[node] || edgesByNode[node].
empty())
continue;
321 const std::uint64_t
key = edgesByNode[node].front().first;
323 if (first >= graph.nNodes || second >= graph.nNodes) continue;
324 if (selectedPairKeys.insert(key).second) ++orphanRecoveryPairs;
325 selectedNode[
first] = 1;
327 }
328 }
329
330 std::vector<std::uint64_t> selectedPairs{selectedPairKeys.begin(),
331 selectedPairKeys.end()};
332 std::sort(selectedPairs.begin(), selectedPairs.end());
333
334 DisjointSet components{graph.nNodes};
335 std::vector<bool> activeNode(graph.nNodes, false);
336 for (const std::uint64_t key : selectedPairs) {
338 components.unite(first, second);
339 activeNode[
first] =
true;
340 activeNode[
second] =
true;
341 }
342
343 std::unordered_map<std::size_t, std::vector<std::size_t>> byRoot;
344 byRoot.reserve(graph.nNodes);
345 for (std::size_t node = 0; node < graph.nNodes; ++node) {
346 if (activeNode[node]) byRoot[components.find(node)].push_back(node);
347 }
348
349
350 std::vector<std::vector<std::size_t>> componentNodes;
351 componentNodes.reserve(byRoot.size());
352 for (auto& [_, nodes] : byRoot) {
354 componentNodes.push_back(std::move(nodes));
355 }
356 std::ranges::sort(componentNodes,
357 [](const auto& first, const auto& second) {
359 });
360
361 xAOD::ContainerDecorator<xAOD::MuonSegmentContainer, std::vector<unsigned>>
363
364 std::vector<std::uint8_t> keptNode(graph.nNodes, false);
365 std::vector<std::uint8_t> nodeDropReason;
366 if (truthDiag) {
367 nodeDropReason.assign(graph.nNodes, kReasonNotSelected);
368 }
369
370 std::size_t topologyNodes = 0;
371 std::size_t retainedNodes = 0;
372 std::size_t chamberSuppressedNodes = 0;
373 std::size_t rejectedComponents = 0;
374 std::size_t componentsKept = 0;
375 std::size_t anchors = 0;
376 std::size_t nodesRejectedByMinComponent = 0;
377 unsigned nextComponentId = 1;
378 const InferenceUtils::SegmentQualityOrder betterSegment{};
379
380 const auto isBetterNode = [&](std::size_t candidate, std::size_t incumbent) {
381 if (betterSegment(graph.segments[candidate], graph.segments[incumbent])) {
382 return true;
383 }
384 if (betterSegment(graph.segments[incumbent], graph.segments[candidate])) {
385 return false;
386 }
387 return candidate < incumbent;
388 };
389
390 for (const std::vector<std::size_t>& rawNodes : componentNodes) {
391 topologyNodes += rawNodes.size();
392 std::vector<std::size_t> retained = rawNodes;
393
394
395
397 std::unordered_map<int, std::size_t> bestByChamber;
398 bestByChamber.reserve(rawNodes.size());
399 for (const std::size_t node : rawNodes) {
401 static_cast<int>(graph.segments[node]->chamberIndex());
402 const auto found = bestByChamber.find(chamber);
403 if (found == bestByChamber.end() || isBetterNode(node,
found->second)) {
405 }
406 }
407 retained.clear();
408 retained.reserve(bestByChamber.size());
409 for (const auto& [_, node] : bestByChamber) retained.push_back(node);
410 std::sort(retained.begin(), retained.end());
411 chamberSuppressedNodes += rawNodes.size() - retained.size();
412 }
413
414
415 if (truthDiag) {
416 for (const std::size_t node : rawNodes) {
417 nodeDropReason[node] = kReasonChamberDedup;
418 }
419 for (const std::size_t node : retained) {
420 nodeDropReason[node] = kReasonComponentRejected;
421 }
422 }
423
425 nodesRejectedByMinComponent += retained.size();
426 ++rejectedComponents;
427 continue;
428 }
429
430
431
432 std::vector<std::size_t> rankedNodes{retained};
434
435 int bestRank = std::numeric_limits<int>::max();
436 for (const std::size_t node : rankedNodes) {
437 bestRank = std::min(
438 bestRank,
440 graph.segments[node]->chamberIndex())));
441 }
444 graph.segments[node]->chamberIndex())) != bestRank;
445 });
446 } else {
447 std::ranges::sort(rankedNodes, isBetterNode);
449 ? rankedNodes.size()
451 rankedNodes.size());
452 rankedNodes.resize(nAnchors);
453 }
454 if (rankedNodes.empty()) {
455 ++rejectedComponents;
456 continue;
457 }
458 std::ranges::sort(rankedNodes);
459
460 const unsigned componentId = nextComponentId++;
461 for (const std::size_t node : retained) {
462 const bool isAnchor = Acts::rangeContainsValue(rankedNodes, node);
463 decor(*graph.segments[node]) = {
464 componentId, static_cast<unsigned int>(isAnchor)};
465 keptNode[node] = true;
466 }
467 retainedNodes += retained.size();
468 anchors += rankedNodes.size();
469 ++componentsKept;
470 }
471
472 if (truthDiag) {
473 std::size_t truthSegs = 0, bkgSegs = 0;
474 for (std::size_t node = 0; node < graph.nNodes; ++node) {
475 if (!keptNode[node]) {
476 continue;
477 }
478 nodeTruthId[node] >= 0 ? ++truthSegs : ++bkgSegs;
479 }
482
483 std::vector<unsigned char> thresholded(graph.nNodes, 0);
484 for (std::size_t node = 0; node < graph.nNodes; ++node) {
485 thresholded[node] = !edgesByNode[node].empty();
486 }
488 keptNode, nodeDropReason);
489 }
490
492 auto connectedSegments =
493 std::make_unique<ConstDataVector<xAOD::MuonSegmentContainer>>(
495 connectedSegments->reserve(graph.nNodes);
496 for (std::size_t node = 0; node < graph.nNodes; ++node) {
497 if (!keptNode[node] || !graph.segments[node]) continue;
498 connectedSegments->push_back(graph.segments[node]);
499 }
500
501 const std::size_t nConnectedSegments = connectedSegments->size();
502 SG::WriteHandle<ConstDataVector<xAOD::MuonSegmentContainer>> connectedHandle{
506 << ": wrote " << nConnectedSegments
507 << " ML-connected segment(s) to '"
509 }
510
512 << ": ML components graphNodes=" << graph.nNodes
513 << ", thresholdedNodes=" << thresholdedNodes
514 << ", thresholdPairs=" << thresholdPairs
515 << ", mutualTopKPairs=" << mutualTopKPairs
516 << ", oneSidedTopKPairs=" << oneSidedTopKPairs
517 << ", selectedPairs=" << selectedPairKeys.size()
518 << ", components=" << componentsKept
519 << ", topologyNodes=" << topologyNodes
520 << ", retainedNodes=" << retainedNodes
521 << ", chamberSuppressedNodes=" << chamberSuppressedNodes
522 << ", nodesRejectedByMinComponent=" << nodesRejectedByMinComponent
524 << ", seedAnchors=" << anchors
525 << ", rejectedComponents=" << rejectedComponents
528 << ", orphanRecoveryPairs=" << orphanRecoveryPairs
531
532
545 }
546 return StatusCode::SUCCESS;
547}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_DEBUG(x,...)
std::vector< std::vector< float > > scores
static const Attributes_t empty
bool msgLvl(const MSG::Level lvl) const
std::atomic< std::size_t > m_sumOnnxBkgSegments
std::atomic< std::size_t > m_sumChamberSuppressedNodes
Gaudi::Property< bool > m_keepBestSegmentPerChamber
std::atomic< std::size_t > m_sumRetainedNodes
Gaudi::Property< bool > m_anchorInnermostLayer
SG::WriteDecorHandleKey< xAOD::MuonSegmentContainer > m_pairGateDecorKey
Per-segment payload consumed by MlMsTrackSeeder: [componentId, isSeedAnchor] Empty means that the seg...
Gaudi::Property< bool > m_useDegreeCappedComponents
std::atomic< std::size_t > m_sumOrphanRecoveryPairs
ToolHandle< ISegmentEdgeClassifierTool > m_edgeClassifier
Gaudi::Property< unsigned int > m_maxEdgesPerNode
std::atomic< std::size_t > m_sumThresholdedNodes
std::atomic< std::size_t > m_sumRetainedBkgSegments
Gaudi::Property< bool > m_requireMutualTopKEdges
std::atomic< std::size_t > m_sumOneSidedTopKPairs
std::atomic< std::size_t > m_sumInputBkgSegments
void classifyLostTruthSegments(const xAOD::MuonSegmentContainer &segments, const SegmentEdgeGraph &graph, const std::unordered_map< std::uint64_t, float > &pairProbability, const std::vector< unsigned char > &thresholded, const std::vector< std::uint8_t > &keptNode, const std::vector< unsigned char > &nodeDropReason) const
MC-only: classify every truth-labeled input segment that did not make it into the filtered container ...
std::atomic< std::size_t > m_sumTrueEdgesPassed
bool m_truthDiagEnabled
Whether truth-particle matching is available/enabled this job.
std::atomic< std::size_t > m_sumBkgEdgesPassed
std::atomic< std::size_t > m_sumBkgEdgesTotal
std::atomic< std::size_t > m_sumGraphNodes
Job-summed post-ONNX selection counters (see execute()).
Gaudi::Property< unsigned int > m_minSegmentsPerComponent
std::atomic< std::size_t > m_sumOnnxTruthSegments
std::atomic< std::size_t > m_sumInputTruthSegments
std::atomic< std::size_t > m_sumNodesRejectedByMinComponent
std::atomic< std::size_t > m_sumTrueEdgesTotal
std::atomic< std::size_t > m_sumRetainedTruthSegments
SG::WriteHandleKey< ConstDataVector< xAOD::MuonSegmentContainer > > m_filteredSegmentKey
Gaudi::Property< bool > m_recoverOrphanNodes
std::atomic< std::size_t > m_sumMutualTopKPairs
std::atomic< std::size_t > m_sumRejectedComponents
SG::ReadHandleKey< xAOD::MuonSegmentContainer > m_segmentKey
Gaudi::Property< unsigned int > m_seedAnchorsPerComponent
Gaudi::Property< float > m_pairGateThreshold
std::atomic< std::size_t > m_sumTopologyNodes
std::atomic< std::size_t > m_sumComponentsKept
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
int count(std::string s, const std::string ®x)
count how many occurances of a regx are in a string
int compareFloatDescending(float first, float second)
Three-way descending comparison which also orders NaN last.
@ VIEW_ELEMENTS
this data object is a view, it does not own its elmts
const T * get(const ReadCondHandleKey< T > &key, const EventContext &ctx)
Convenience function to retrieve an object given a ReadCondHandleKey.
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
std::size_t erase_if(T_container &container, T_Func pred)
MuonSegmentContainer_v1 MuonSegmentContainer
Definition of the current "MuonSegment container version".