From 3c8f6257b48237d2357d28dfd7f9377d86fc2d0f Mon Sep 17 00:00:00 2001 From: Sam Zhou Date: Wed, 29 Jul 2026 12:23:27 -0500 Subject: [PATCH 1/3] Add CaloEntrantTruthMaker: per-hit calo-entrant truth product For each CaloHitMC energy deposit, walk SimParticle::parent() and store the calo-entrant ancestor (highest ancestor that also deposited in the same disk) in a new transient CaloHitEntrantCollection, index-parallel to the input CaloHitMCCollection. Crystal/disk resolution prefers CaloHitMC::crystalID() when filled; files produced before that member existed (reads back -1) fall back to the CaloCluster<->CaloClusterMC pairing. Cardinality or null-Ptr violations of that pairing throw rather than degrade. Analysis choices (purity cut, ambiguity, cluster IDs) deliberately stay downstream; the product records ancestry facts only. Registered in both build systems (CMakeLists cet_build_plugin + SConscript; MF_MessageLogger added to the SCons plugin link list for the mf logging call). Validated byte-exact against the Python ancestry-walk reference on 200 MDC2025 events under envset p103: 6554/6554 contribution roots identical (fingerprint-matched across ntuple formats). --- CaloMC/CMakeLists.txt | 10 + CaloMC/src/CaloEntrantTruthMaker_module.cc | 203 +++++++++++++++++++++ CaloMC/src/SConscript | 1 + MCDataProducts/inc/CaloHitEntrant.hh | 36 ++++ 4 files changed, 250 insertions(+) create mode 100644 CaloMC/src/CaloEntrantTruthMaker_module.cc create mode 100644 MCDataProducts/inc/CaloHitEntrant.hh diff --git a/CaloMC/CMakeLists.txt b/CaloMC/CMakeLists.txt index 9f9468b0a3..cd96137767 100644 --- a/CaloMC/CMakeLists.txt +++ b/CaloMC/CMakeLists.txt @@ -85,6 +85,16 @@ cet_build_plugin(CaloShowerUpdater art::module Offline::MCDataProducts ) +cet_build_plugin(CaloEntrantTruthMaker art::module + REG_SOURCE src/CaloEntrantTruthMaker_module.cc + LIBRARIES REG + + Offline::CalorimeterGeom + Offline::GeometryService + Offline::MCDataProducts + Offline::RecoDataProducts +) + configure_file(${CMAKE_CURRENT_SOURCE_DIR}/fcl/prolog.fcl ${CURRENT_BINARY_DIR} fcl/prolog.fcl) install_source(SUBDIRS src) diff --git a/CaloMC/src/CaloEntrantTruthMaker_module.cc b/CaloMC/src/CaloEntrantTruthMaker_module.cc new file mode 100644 index 0000000000..bbcadfffd2 --- /dev/null +++ b/CaloMC/src/CaloEntrantTruthMaker_module.cc @@ -0,0 +1,203 @@ +// +// Assign each calorimeter MC energy deposit to its calo-entrant ancestor: +// the highest SimParticle in the Geant4 parent chain that also deposited +// energy in the same calorimeter disk (the shower originator). +// +// For each CaloHitMC in the input collection this module emits a +// CaloHitEntrant whose entrants vector is aligned with +// CaloHitMC::energyDeposits(); the output collection is index-parallel +// to the input CaloHitMCCollection. +// +// Crystal (hence disk) resolution prefers CaloHitMC::crystalID() when it +// is filled. For files produced before that member existed (reads back +// -1) the crystal is resolved through the CaloCluster <-> CaloClusterMC +// pairing (index-parallel collections with positionally-matched hit +// lists, the invariant CaloClusterTruthMatch establishes). CaloHitMC +// entries whose crystal cannot be determined either way get null +// (unresolved) entrants. +// +// Grouping hits by entrant collapses secondary shower products +// (bremsstrahlung photons etc.) into their parent shower, recovering +// true shower membership for clustering truth definitions. Purity cuts +// and ambiguity handling are left to consumers. +// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "fhiclcpp/types/Atom.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "cetlib_except/exception.h" + +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" +#include "Offline/RecoDataProducts/inc/CaloCluster.hh" +#include "Offline/MCDataProducts/inc/CaloHitMC.hh" +#include "Offline/MCDataProducts/inc/CaloClusterMC.hh" +#include "Offline/MCDataProducts/inc/CaloHitEntrant.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" + +#include +#include +#include +#include + +namespace mu2e { + + class CaloEntrantTruthMaker : public art::EDProducer + { + public: + struct Config + { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + fhicl::Atom caloHitMCTag { Name("caloHitMCTag"), Comment("CaloHitMCCollection input tag") }; + fhicl::Atom caloClusterTag { Name("caloClusterTag"), Comment("CaloClusterCollection input tag (legacy crystal resolution)") }; + fhicl::Atom caloClusterMCTag { Name("caloClusterMCTag"), Comment("CaloClusterMCCollection input tag (legacy crystal resolution)") }; + fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diagnostic level"), 0 }; + }; + + explicit CaloEntrantTruthMaker(const art::EDProducer::Table& config) : + EDProducer{config}, + caloHitMCToken_ {consumes(config().caloHitMCTag())}, + caloClusterToken_ {consumes(config().caloClusterTag())}, + caloClusterMCToken_{consumes(config().caloClusterMCTag())}, + diagLevel_ (config().diagLevel()) + { + produces(); + } + + void produce(art::Event& event) override; + + private: + art::ProductToken caloHitMCToken_; + art::ProductToken caloClusterToken_; + art::ProductToken caloClusterMCToken_; + int diagLevel_; + }; + + //-------------------------------------------------------------------- + void CaloEntrantTruthMaker::produce(art::Event& event) + { + const auto& hitMCs = *event.getValidHandle(caloHitMCToken_); + const auto& clusters = *event.getValidHandle(caloClusterToken_); + const auto& clusterMCs = *event.getValidHandle(caloClusterMCToken_); + const Calorimeter& cal = *(GeomHandle()); + + // Legacy crystal resolution: CaloCluster <-> CaloClusterMC pairing. + // CaloClusterTruthMatch emits exactly one CaloClusterMC per + // CaloCluster in input order with one CaloHitMC per hit, so any + // cardinality mismatch is a wiring error, not a soft condition. + if (clusters.size() != clusterMCs.size()) + { + throw cet::exception("CALOENTRANT") + << "CaloCluster (" << clusters.size() << ") and CaloClusterMC (" + << clusterMCs.size() << ") collections are not parallel\n"; + } + std::map hitmcCrystal; + for (size_t ic = 0; ic < clusters.size(); ++ic) + { + const auto& hits = clusters[ic].caloHitsPtrVector(); + const auto& hitmcs = clusterMCs[ic].caloHitMCs(); + if (hits.size() != hitmcs.size()) + { + throw cet::exception("CALOENTRANT") + << "cluster " << ic << ": " << hits.size() << " CaloHits vs " + << hitmcs.size() << " CaloHitMCs\n"; + } + for (size_t ih = 0; ih < hits.size(); ++ih) + { + if (hits[ih].isNull() || hitmcs[ih].isNull()) + { + throw cet::exception("CALOENTRANT") + << "cluster " << ic << " hit " << ih << ": null Ptr in cluster hit lists\n"; + } + hitmcCrystal[hitmcs[ih].key()] = hits[ih]->crystalID(); + } + } + + // Crystal of a CaloHitMC: the member when filled (newer + // productions), else the cluster-pairing lookup (legacy files). + auto crystalOf = [&](size_t idx, const CaloHitMC& chmc) -> int + { + if (chmc.crystalID() >= 0) return chmc.crystalID(); + const auto found = hitmcCrystal.find(idx); + return (found != hitmcCrystal.end()) ? found->second : -1; + }; + + // Pass 1: record which disks each SimParticle deposited in. + // Keyed by SimParticle id to match the id convention exposed to + // ntuple consumers (SimParticle::id().asInt()). + std::map> simDisks; + for (size_t idx = 0; idx < hitMCs.size(); ++idx) + { + const int crystal = crystalOf(idx, hitMCs[idx]); + if (crystal < 0) continue; + int disk = cal.crystal(crystal).diskID(); + for (const auto& edep : hitMCs[idx].energyDeposits()) + { + if (edep.sim().isNonnull()) simDisks[edep.sim()->id().asInt()].insert(disk); + } + } + + // Pass 2: for each deposit, walk the parent chain upward; the + // entrant is the highest ancestor that also deposited in the same + // disk. The walk stops at the chain root or at an unresolvable + // parent pointer (event compression). + auto output = std::make_unique(); + output->reserve(hitMCs.size()); + std::map, art::Ptr> entrantCache; // (simId, disk) -> entrant + + for (size_t idx = 0; idx < hitMCs.size(); ++idx) + { + const auto& chmc = hitMCs[idx]; + CaloHitEntrant che; + + const int crystal = crystalOf(idx, chmc); + if (crystal < 0) + { + // not referenced by any cluster and no crystal member: + // entrants unresolvable; keep the collection index-parallel + // with null entries + che.entrants.assign(chmc.nParticles(), art::Ptr()); + output->push_back(std::move(che)); + continue; + } + const int disk = cal.crystal(crystal).diskID(); + che.entrants.reserve(chmc.nParticles()); + + for (const auto& edep : chmc.energyDeposits()) + { + const art::Ptr& sim = edep.sim(); + if (sim.isNull()) { che.entrants.emplace_back(); continue; } + + const auto key = std::make_pair(sim->id().asInt(), disk); + auto it = entrantCache.find(key); + if (it == entrantCache.end()) + { + art::Ptr entrant = sim; + for (auto p = sim->parent(); p.isNonnull(); p = p->parent()) + { + auto fnd = simDisks.find(p->id().asInt()); + if (fnd != simDisks.end() && fnd->second.count(disk)) entrant = p; + } + it = entrantCache.emplace(key, entrant).first; + } + che.entrants.push_back(it->second); + } + output->push_back(std::move(che)); + } + + if (diagLevel_ > 0) + { + mf::LogInfo("CaloEntrantTruthMaker") << "produced " << output->size() + << " CaloHitEntrant for " << hitMCs.size() << " CaloHitMC" + << " (" << hitmcCrystal.size() << " cluster-referenced)"; + } + + event.put(std::move(output)); + } +} + +DEFINE_ART_MODULE(mu2e::CaloEntrantTruthMaker) diff --git a/CaloMC/src/SConscript b/CaloMC/src/SConscript index 5e88eeca92..b7a8da22e9 100644 --- a/CaloMC/src/SConscript +++ b/CaloMC/src/SConscript @@ -69,6 +69,7 @@ helper.make_plugins( [ mainlib, 'art_Persistency_Provenance', 'art_Utilities', 'canvas', + 'MF_MessageLogger', 'fhiclcpp', 'fhiclcpp_types', 'tbb', diff --git a/MCDataProducts/inc/CaloHitEntrant.hh b/MCDataProducts/inc/CaloHitEntrant.hh new file mode 100644 index 0000000000..3505afb092 --- /dev/null +++ b/MCDataProducts/inc/CaloHitEntrant.hh @@ -0,0 +1,36 @@ +#ifndef MCDataProducts_CaloHitEntrant_hh +#define MCDataProducts_CaloHitEntrant_hh +// +// Calo-entrant truth assignment for one CaloHitMC. +// +// For each MC energy deposit in the hit (aligned with +// CaloHitMC::energyDeposits()), stores the SimParticle that originated +// the shower on this disk: the highest ancestor in the Geant4 parent +// chain that also deposited energy in the same calorimeter disk. +// A particle with no such ancestor (e.g. a cross-disk secondary or a +// pileup particle arriving from outside) is its own calo-entrant. +// +// Grouping deposits by entrant recovers true shower membership for +// clustering truth definitions. Analysis-level choices (purity cuts, +// ambiguity handling, cluster ID assignment) are deliberately left to +// consumers; this product records only the ancestry facts. +// +// The collection is index-parallel to the CaloHitMCCollection it was +// produced from. +// + +#include "canvas/Persistency/Common/Ptr.h" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include + +namespace mu2e +{ + struct CaloHitEntrant + { + std::vector> entrants; + }; + + using CaloHitEntrantCollection = std::vector; +} + +#endif From d3ba04d8de4c48d94fbdc8006843fcab3577075e Mon Sep 17 00:00:00 2001 From: Sam Zhou Date: Sat, 1 Aug 2026 05:36:45 -0500 Subject: [PATCH 2/3] Address review: cluster-level disk lookup, dictionary, optional tags Review feedback on #1911: - Legacy fallback no longer pairs cluster hit lists positionally. CaloClusterTruthMatch filters unmatched hits and re-sorts the MC list by energy, so only the collection-level zip is contracted. Each CaloHitMC referenced by a CaloClusterMC now inherits CaloCluster::diskID() directly (the disk is all the pairing was used for), tolerating unmatched hits instead of throwing on them. A ProductID cross-check rejects a CaloClusterMC built against a different CaloHitMCCollection than caloHitMCTag. - CaloHitEntrant gets a ROOT dictionary (classes.h + classes_def.xml) and is persistable like every other Mu2e data product. Each entry carries an art::Ptr back-reference, so the product is self-describing and consumers can verify source identity instead of trusting configuration. Members are private with accessors. - caloClusterTag/caloClusterMCTag are optional (mayConsume, empty default, must be set together): on files with CaloHitMC::crystalID() filled the module runs without any cluster products. A once-per-job LogWarning flags legacy files with no fallback configured (entrants null). - simDisks/entrantCache keyed by art::Ptr (ProductID + key) instead of bare SimParticle id, removing the single-collection assumption; the parent walk throws a clear message on an unreadable parent Ptr (dropped ancestor collection); conventional Parameters alias added. Re-validated against the Python ancestry-walk reference on 200 MDC2025 events: 3903/3903 MC hits matched, 6554/6554 contribution roots identical. Smoke-tested: no-cluster-tags job runs with all-null entrants plus warning; a keep-everything RootOutput job writes the product branch (write-side smoke; full read-back done in independent review). --- CaloMC/src/CaloEntrantTruthMaker_module.cc | 182 ++++++++++++--------- MCDataProducts/inc/CaloHitEntrant.hh | 28 +++- MCDataProducts/src/classes.h | 1 + MCDataProducts/src/classes_def.xml | 5 + 4 files changed, 138 insertions(+), 78 deletions(-) diff --git a/CaloMC/src/CaloEntrantTruthMaker_module.cc b/CaloMC/src/CaloEntrantTruthMaker_module.cc index bbcadfffd2..f8799faca5 100644 --- a/CaloMC/src/CaloEntrantTruthMaker_module.cc +++ b/CaloMC/src/CaloEntrantTruthMaker_module.cc @@ -5,16 +5,21 @@ // // For each CaloHitMC in the input collection this module emits a // CaloHitEntrant whose entrants vector is aligned with -// CaloHitMC::energyDeposits(); the output collection is index-parallel -// to the input CaloHitMCCollection. +// CaloHitMC::energyDeposits() and which carries an art::Ptr back to its +// CaloHitMC; the output collection is index-parallel to the input +// CaloHitMCCollection. // -// Crystal (hence disk) resolution prefers CaloHitMC::crystalID() when it -// is filled. For files produced before that member existed (reads back -// -1) the crystal is resolved through the CaloCluster <-> CaloClusterMC -// pairing (index-parallel collections with positionally-matched hit -// lists, the invariant CaloClusterTruthMatch establishes). CaloHitMC -// entries whose crystal cannot be determined either way get null -// (unresolved) entrants. +// Disk resolution prefers CaloHitMC::crystalID() when it is filled. For +// files produced before that member existed (reads back -1) the disk is +// resolved through the optional CaloCluster / CaloClusterMC pair: +// CaloClusterTruthMatch emits exactly one CaloClusterMC per CaloCluster +// in input order, so the two collections zip by index, and every +// CaloHitMC referenced by a CaloClusterMC inherits the corresponding +// CaloCluster::diskID(). The per-cluster hit lists themselves are NOT +// positionally matched (the MC list keeps only truth-matched hits and is +// re-sorted by MC energy), so no hit-by-hit pairing is attempted; only +// the disk, a cluster-level property, is taken. CaloHitMC entries whose +// disk cannot be determined either way get null (unresolved) entrants. // // Grouping hits by entrant collapses secondary shower products // (bremsstrahlung photons etc.) into their parent shower, recovering @@ -42,6 +47,7 @@ #include #include #include +#include namespace mu2e { @@ -53,18 +59,27 @@ namespace mu2e { using Name = fhicl::Name; using Comment = fhicl::Comment; fhicl::Atom caloHitMCTag { Name("caloHitMCTag"), Comment("CaloHitMCCollection input tag") }; - fhicl::Atom caloClusterTag { Name("caloClusterTag"), Comment("CaloClusterCollection input tag (legacy crystal resolution)") }; - fhicl::Atom caloClusterMCTag { Name("caloClusterMCTag"), Comment("CaloClusterMCCollection input tag (legacy crystal resolution)") }; + fhicl::Atom caloClusterTag { Name("caloClusterTag"), Comment("CaloClusterCollection input tag (legacy disk resolution; empty disables)"), art::InputTag() }; + fhicl::Atom caloClusterMCTag { Name("caloClusterMCTag"), Comment("CaloClusterMCCollection input tag (legacy disk resolution; empty disables)"), art::InputTag() }; fhicl::Atom diagLevel { Name("diagLevel"), Comment("Diagnostic level"), 0 }; }; + using Parameters = art::EDProducer::Table; - explicit CaloEntrantTruthMaker(const art::EDProducer::Table& config) : + // member initialization follows declaration order below: + // tokens, useClusterFallback_, diagLevel_ + explicit CaloEntrantTruthMaker(const Parameters& config) : EDProducer{config}, caloHitMCToken_ {consumes(config().caloHitMCTag())}, - caloClusterToken_ {consumes(config().caloClusterTag())}, - caloClusterMCToken_{consumes(config().caloClusterMCTag())}, + caloClusterToken_ {mayConsume(config().caloClusterTag())}, + caloClusterMCToken_{mayConsume(config().caloClusterMCTag())}, + useClusterFallback_(!config().caloClusterTag().empty()), diagLevel_ (config().diagLevel()) { + if (config().caloClusterTag().empty() != config().caloClusterMCTag().empty()) + { + throw cet::exception("CALOENTRANT") + << "caloClusterTag and caloClusterMCTag must be both set or both empty\n"; + } produces(); } @@ -74,126 +89,147 @@ namespace mu2e { art::ProductToken caloHitMCToken_; art::ProductToken caloClusterToken_; art::ProductToken caloClusterMCToken_; - int diagLevel_; + bool useClusterFallback_; + int diagLevel_; + bool warnedNoFallback_ = false; }; //-------------------------------------------------------------------- void CaloEntrantTruthMaker::produce(art::Event& event) { - const auto& hitMCs = *event.getValidHandle(caloHitMCToken_); - const auto& clusters = *event.getValidHandle(caloClusterToken_); - const auto& clusterMCs = *event.getValidHandle(caloClusterMCToken_); - const Calorimeter& cal = *(GeomHandle()); - - // Legacy crystal resolution: CaloCluster <-> CaloClusterMC pairing. - // CaloClusterTruthMatch emits exactly one CaloClusterMC per - // CaloCluster in input order with one CaloHitMC per hit, so any - // cardinality mismatch is a wiring error, not a soft condition. - if (clusters.size() != clusterMCs.size()) + const auto hitMCHandle = event.getValidHandle(caloHitMCToken_); + const auto& hitMCs = *hitMCHandle; + const Calorimeter& cal = *(GeomHandle()); + + // Legacy disk resolution: every CaloHitMC referenced by a + // CaloClusterMC inherits the disk of the corresponding CaloCluster. + // Only the collection-level zip is an invariant of + // CaloClusterTruthMatch (one CaloClusterMC per CaloCluster, input + // order); the per-cluster hit lists are filtered and re-sorted, so + // they are never paired positionally here. + std::map hitmcDisk; + if (useClusterFallback_) { - throw cet::exception("CALOENTRANT") - << "CaloCluster (" << clusters.size() << ") and CaloClusterMC (" - << clusterMCs.size() << ") collections are not parallel\n"; - } - std::map hitmcCrystal; - for (size_t ic = 0; ic < clusters.size(); ++ic) - { - const auto& hits = clusters[ic].caloHitsPtrVector(); - const auto& hitmcs = clusterMCs[ic].caloHitMCs(); - if (hits.size() != hitmcs.size()) + const auto& clusters = *event.getValidHandle(caloClusterToken_); + const auto& clusterMCs = *event.getValidHandle(caloClusterMCToken_); + if (clusters.size() != clusterMCs.size()) { throw cet::exception("CALOENTRANT") - << "cluster " << ic << ": " << hits.size() << " CaloHits vs " - << hitmcs.size() << " CaloHitMCs\n"; + << "CaloCluster (" << clusters.size() << ") and CaloClusterMC (" + << clusterMCs.size() << ") collections are not parallel\n"; } - for (size_t ih = 0; ih < hits.size(); ++ih) + for (size_t ic = 0; ic < clusters.size(); ++ic) { - if (hits[ih].isNull() || hitmcs[ih].isNull()) + const int disk = clusters[ic].diskID(); + for (const auto& hmc : clusterMCs[ic].caloHitMCs()) { - throw cet::exception("CALOENTRANT") - << "cluster " << ic << " hit " << ih << ": null Ptr in cluster hit lists\n"; + if (hmc.isNull()) continue; + if (hmc.id() != hitMCHandle.id()) + { + throw cet::exception("CALOENTRANT") + << "CaloClusterMC references CaloHitMCCollection " << hmc.id() + << " but caloHitMCTag resolves to " << hitMCHandle.id() + << "; caloClusterMCTag and caloHitMCTag are inconsistent\n"; + } + hitmcDisk[hmc.key()] = disk; } - hitmcCrystal[hitmcs[ih].key()] = hits[ih]->crystalID(); } } - // Crystal of a CaloHitMC: the member when filled (newer - // productions), else the cluster-pairing lookup (legacy files). - auto crystalOf = [&](size_t idx, const CaloHitMC& chmc) -> int + // Disk of a CaloHitMC: from the crystalID member when filled (newer + // productions), else the cluster lookup (legacy files), else -1. + auto diskOf = [&](size_t idx, const CaloHitMC& chmc) -> int { - if (chmc.crystalID() >= 0) return chmc.crystalID(); - const auto found = hitmcCrystal.find(idx); - return (found != hitmcCrystal.end()) ? found->second : -1; + if (chmc.crystalID() >= 0) return cal.crystal(chmc.crystalID()).diskID(); + const auto found = hitmcDisk.find(idx); + return (found != hitmcDisk.end()) ? found->second : -1; }; - // Pass 1: record which disks each SimParticle deposited in. - // Keyed by SimParticle id to match the id convention exposed to - // ntuple consumers (SimParticle::id().asInt()). - std::map> simDisks; + // Pass 1: record which disks each SimParticle deposited in. Keyed by + // art::Ptr (ProductID + key) so deposits referencing more than one + // SimParticle collection cannot collide. + std::map, std::set> simDisks; for (size_t idx = 0; idx < hitMCs.size(); ++idx) { - const int crystal = crystalOf(idx, hitMCs[idx]); - if (crystal < 0) continue; - int disk = cal.crystal(crystal).diskID(); + const int disk = diskOf(idx, hitMCs[idx]); + if (disk < 0) continue; for (const auto& edep : hitMCs[idx].energyDeposits()) { - if (edep.sim().isNonnull()) simDisks[edep.sim()->id().asInt()].insert(disk); + if (edep.sim().isNonnull()) simDisks[edep.sim()].insert(disk); } } // Pass 2: for each deposit, walk the parent chain upward; the // entrant is the highest ancestor that also deposited in the same - // disk. The walk stops at the chain root or at an unresolvable - // parent pointer (event compression). + // disk. The walk stops at the chain root; a parent Ptr that cannot + // be read back (ancestor collection dropped from the file) throws. auto output = std::make_unique(); output->reserve(hitMCs.size()); - std::map, art::Ptr> entrantCache; // (simId, disk) -> entrant + std::map, int>, art::Ptr> entrantCache; + size_t nUnresolved(0); for (size_t idx = 0; idx < hitMCs.size(); ++idx) { const auto& chmc = hitMCs[idx]; - CaloHitEntrant che; + std::vector> entrants; - const int crystal = crystalOf(idx, chmc); - if (crystal < 0) + const int disk = diskOf(idx, chmc); + if (disk < 0) { - // not referenced by any cluster and no crystal member: + // no crystalID member and not referenced by any cluster: // entrants unresolvable; keep the collection index-parallel // with null entries - che.entrants.assign(chmc.nParticles(), art::Ptr()); - output->push_back(std::move(che)); + entrants.assign(chmc.nParticles(), art::Ptr()); + ++nUnresolved; + output->emplace_back(art::Ptr(hitMCHandle, idx), std::move(entrants)); continue; } - const int disk = cal.crystal(crystal).diskID(); - che.entrants.reserve(chmc.nParticles()); + entrants.reserve(chmc.nParticles()); for (const auto& edep : chmc.energyDeposits()) { const art::Ptr& sim = edep.sim(); - if (sim.isNull()) { che.entrants.emplace_back(); continue; } + if (sim.isNull()) { entrants.emplace_back(); continue; } - const auto key = std::make_pair(sim->id().asInt(), disk); + const auto key = std::make_pair(sim, disk); auto it = entrantCache.find(key); if (it == entrantCache.end()) { art::Ptr entrant = sim; for (auto p = sim->parent(); p.isNonnull(); p = p->parent()) { - auto fnd = simDisks.find(p->id().asInt()); + if (!p.isAvailable()) + { + throw cet::exception("CALOENTRANT") + << "CaloEntrantTruthMaker: SimParticle parent (ProductID " << p.id() + << ", key " << p.key() << ") cannot be read from this file" + << " - an ancestor SimParticle collection was dropped;" + << " cannot walk the parent chain\n"; + } + const auto fnd = simDisks.find(p); if (fnd != simDisks.end() && fnd->second.count(disk)) entrant = p; } it = entrantCache.emplace(key, entrant).first; } - che.entrants.push_back(it->second); + entrants.push_back(it->second); } - output->push_back(std::move(che)); + output->emplace_back(art::Ptr(hitMCHandle, idx), std::move(entrants)); } + if (nUnresolved > 0 && !useClusterFallback_ && !warnedNoFallback_) + { + warnedNoFallback_ = true; + mf::LogWarning("CaloEntrantTruthMaker") + << nUnresolved << " of " << hitMCs.size() << " CaloHitMC have no crystalID" + << " and no cluster fallback is configured; their entrants are null." + << " Set caloClusterTag/caloClusterMCTag to resolve disks on legacy files."; + } if (diagLevel_ > 0) { mf::LogInfo("CaloEntrantTruthMaker") << "produced " << output->size() << " CaloHitEntrant for " << hitMCs.size() << " CaloHitMC" - << " (" << hitmcCrystal.size() << " cluster-referenced)"; + << " (" << hitmcDisk.size() << " cluster-referenced, " + << nUnresolved << " unresolved)"; } event.put(std::move(output)); diff --git a/MCDataProducts/inc/CaloHitEntrant.hh b/MCDataProducts/inc/CaloHitEntrant.hh index 3505afb092..845f3f71fc 100644 --- a/MCDataProducts/inc/CaloHitEntrant.hh +++ b/MCDataProducts/inc/CaloHitEntrant.hh @@ -10,24 +10,42 @@ // A particle with no such ancestor (e.g. a cross-disk secondary or a // pileup particle arriving from outside) is its own calo-entrant. // +// caloHitMC() points back to the hit this assignment was computed for, +// so the product is self-describing: consumers join through the Ptr +// (or verify it) instead of trusting that their configured +// CaloHitMCCollection is the one the producer read. The collection is +// also index-parallel to that CaloHitMCCollection by construction. +// // Grouping deposits by entrant recovers true shower membership for // clustering truth definitions. Analysis-level choices (purity cuts, // ambiguity handling, cluster ID assignment) are deliberately left to // consumers; this product records only the ancestry facts. // -// The collection is index-parallel to the CaloHitMCCollection it was -// produced from. -// #include "canvas/Persistency/Common/Ptr.h" +#include "Offline/MCDataProducts/inc/CaloHitMC.hh" #include "Offline/MCDataProducts/inc/SimParticle.hh" +#include #include namespace mu2e { - struct CaloHitEntrant + class CaloHitEntrant { - std::vector> entrants; + public: + CaloHitEntrant() = default; + CaloHitEntrant(art::Ptr caloHitMC, + std::vector> entrants) : + caloHitMC_(std::move(caloHitMC)), + entrants_ (std::move(entrants)) + {} + + const art::Ptr& caloHitMC() const {return caloHitMC_;} + const std::vector>& entrants () const {return entrants_;} + + private: + art::Ptr caloHitMC_; + std::vector> entrants_; }; using CaloHitEntrantCollection = std::vector; diff --git a/MCDataProducts/src/classes.h b/MCDataProducts/src/classes.h index d97d8ddd29..599c746f68 100644 --- a/MCDataProducts/src/classes.h +++ b/MCDataProducts/src/classes.h @@ -60,6 +60,7 @@ #include "Offline/MCDataProducts/inc/CaloShowerStep.hh" #include "Offline/MCDataProducts/inc/CaloShowerSim.hh" #include "Offline/MCDataProducts/inc/CaloHitMC.hh" +#include "Offline/MCDataProducts/inc/CaloHitEntrant.hh" #include "Offline/MCDataProducts/inc/CaloClusterMC.hh" #include "Offline/MCDataProducts/inc/CaloEDepMC.hh" // straws diff --git a/MCDataProducts/src/classes_def.xml b/MCDataProducts/src/classes_def.xml index 350fc1d0d6..e69a90d402 100644 --- a/MCDataProducts/src/classes_def.xml +++ b/MCDataProducts/src/classes_def.xml @@ -195,6 +195,11 @@ include="Math/Vector3D.h, Math/Vector4D.h, CLHEP/Vector/LorentzVector.h, CLHEP/V + + + + + From af6d54c1c3f11c56f61a729e9aba2da45dd80e66 Mon Sep 17 00:00:00 2001 From: Sam Zhou Date: Sun, 2 Aug 2026 11:11:29 -0500 Subject: [PATCH 3/3] Harden the cluster-disk fallback and guard the depositor Ptr Independent review of the previous commit found three gaps: - The cluster-wide disk fallback assumed disk-local clusters, but association strategy 2 (ClusterAssociator) can merge proto-clusters across disks while CaloCluster keeps only the seed's disk; such a cluster would have received silently wrong truth. The fallback now verifies every reco hit of a cluster sits on the cluster's disk and throws on a mixed-disk cluster (the filtered, energy-sorted MC list cannot recover per-hit disks there). Configure the cluster tags only when the fallback is actually needed. - The HMC-key -> disk map was unchecked last-write-wins. Keys are now bounds-checked against the collection, and a conflicting second disk assignment throws with both cluster indices and both disks. - The parent-chain walk guarded ancestor Ptrs but dereferenced the depositor SimParticle Ptr unguarded; sim.isAvailable() is now checked before the first dereference, with the same diagnostic style. No behavior change on disk-local, uniquely-referenced input: the same 200 MDC2025 events revalidate identically (200/200 events, 3903/3903 MC hits, 6554/6554 contribution roots, PARITY OK; no-tags and RootOutput smokes re-pass). --- CaloMC/src/CaloEntrantTruthMaker_module.cc | 57 ++++++++++++++++++++-- 1 file changed, 52 insertions(+), 5 deletions(-) diff --git a/CaloMC/src/CaloEntrantTruthMaker_module.cc b/CaloMC/src/CaloEntrantTruthMaker_module.cc index f8799faca5..ea06010f29 100644 --- a/CaloMC/src/CaloEntrantTruthMaker_module.cc +++ b/CaloMC/src/CaloEntrantTruthMaker_module.cc @@ -18,8 +18,13 @@ // CaloCluster::diskID(). The per-cluster hit lists themselves are NOT // positionally matched (the MC list keeps only truth-matched hits and is // re-sorted by MC energy), so no hit-by-hit pairing is attempted; only -// the disk, a cluster-level property, is taken. CaloHitMC entries whose -// disk cannot be determined either way get null (unresolved) entrants. +// the disk, a cluster-level property, is taken. That is only sound for +// disk-local clusters, which is verified against the reco hits' crystals: +// a mixed-disk cluster (possible under cluster-association strategy 2, +// which does not enforce same-disk membership) throws rather than +// receiving a guessed disk. Configure the cluster tags only when the +// fallback is actually needed. CaloHitMC entries whose disk cannot be +// determined either way get null (unresolved) entrants. // // Grouping hits by entrant collapses secondary shower products // (bremsstrahlung photons etc.) into their parent shower, recovering @@ -107,7 +112,7 @@ namespace mu2e { // CaloClusterTruthMatch (one CaloClusterMC per CaloCluster, input // order); the per-cluster hit lists are filtered and re-sorted, so // they are never paired positionally here. - std::map hitmcDisk; + std::map> hitmcDisk; // HMC key -> (disk, first claiming cluster) if (useClusterFallback_) { const auto& clusters = *event.getValidHandle(caloClusterToken_); @@ -121,6 +126,26 @@ namespace mu2e { for (size_t ic = 0; ic < clusters.size(); ++ic) { const int disk = clusters[ic].diskID(); + + // The cluster-wide disk is only meaningful for a disk-local + // cluster. Association strategy 2 (ClusterAssociator) can merge + // proto-clusters across disks while CaloCluster keeps only the + // seed's disk; the filtered, re-sorted MC list cannot recover + // per-hit disks for such a cluster, so fail closed instead of + // assigning a guessed disk. + for (const auto& hit : clusters[ic].caloHitsPtrVector()) + { + if (hit.isNull()) continue; + const int hitDisk = cal.crystal(hit->crystalID()).diskID(); + if (hitDisk != disk) + { + throw cet::exception("CALOENTRANT") + << "cluster " << ic << " is not disk-local: diskID " << disk + << " but hit crystal " << hit->crystalID() << " is on disk " + << hitDisk << "; the cluster-wide disk fallback cannot label it\n"; + } + } + for (const auto& hmc : clusterMCs[ic].caloHitMCs()) { if (hmc.isNull()) continue; @@ -131,7 +156,21 @@ namespace mu2e { << " but caloHitMCTag resolves to " << hitMCHandle.id() << "; caloClusterMCTag and caloHitMCTag are inconsistent\n"; } - hitmcDisk[hmc.key()] = disk; + if (hmc.key() >= hitMCs.size()) + { + throw cet::exception("CALOENTRANT") + << "CaloClusterMC " << ic << " references CaloHitMC key " << hmc.key() + << " outside collection size " << hitMCs.size() << "\n"; + } + const auto ins = hitmcDisk.emplace(hmc.key(), std::make_pair(disk, ic)); + if (!ins.second && ins.first->second.first != disk) + { + throw cet::exception("CALOENTRANT") + << "CaloHitMC key " << hmc.key() << " (ProductID " << hmc.id() + << ") is claimed by cluster " << ins.first->second.second + << " on disk " << ins.first->second.first + << " and by cluster " << ic << " on disk " << disk << "\n"; + } } } } @@ -142,7 +181,7 @@ namespace mu2e { { if (chmc.crystalID() >= 0) return cal.crystal(chmc.crystalID()).diskID(); const auto found = hitmcDisk.find(idx); - return (found != hitmcDisk.end()) ? found->second : -1; + return (found != hitmcDisk.end()) ? found->second.first : -1; }; // Pass 1: record which disks each SimParticle deposited in. Keyed by @@ -195,6 +234,14 @@ namespace mu2e { auto it = entrantCache.find(key); if (it == entrantCache.end()) { + if (!sim.isAvailable()) + { + throw cet::exception("CALOENTRANT") + << "CaloEntrantTruthMaker: depositor SimParticle (ProductID " << sim.id() + << ", key " << sim.key() << ") cannot be read from this file" + << " - its SimParticle collection was dropped;" + << " cannot walk the parent chain\n"; + } art::Ptr entrant = sim; for (auto p = sim->parent(); p.isNonnull(); p = p->parent()) {