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..ea06010f29 --- /dev/null +++ b/CaloMC/src/CaloEntrantTruthMaker_module.cc @@ -0,0 +1,286 @@ +// +// 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() and which carries an art::Ptr back to its +// CaloHitMC; the output collection is index-parallel to the input +// CaloHitMCCollection. +// +// 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. 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 +// 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 +#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 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; + + // member initialization follows declaration order below: + // tokens, useClusterFallback_, diagLevel_ + explicit CaloEntrantTruthMaker(const Parameters& config) : + EDProducer{config}, + caloHitMCToken_ {consumes(config().caloHitMCTag())}, + 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(); + } + + void produce(art::Event& event) override; + + private: + art::ProductToken caloHitMCToken_; + art::ProductToken caloClusterToken_; + art::ProductToken caloClusterMCToken_; + bool useClusterFallback_; + int diagLevel_; + bool warnedNoFallback_ = false; + }; + + //-------------------------------------------------------------------- + void CaloEntrantTruthMaker::produce(art::Event& event) + { + 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; // HMC key -> (disk, first claiming cluster) + if (useClusterFallback_) + { + const auto& clusters = *event.getValidHandle(caloClusterToken_); + const auto& clusterMCs = *event.getValidHandle(caloClusterMCToken_); + if (clusters.size() != clusterMCs.size()) + { + throw cet::exception("CALOENTRANT") + << "CaloCluster (" << clusters.size() << ") and CaloClusterMC (" + << clusterMCs.size() << ") collections are not parallel\n"; + } + 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; + 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"; + } + 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"; + } + } + } + } + + // 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 cal.crystal(chmc.crystalID()).diskID(); + const auto found = hitmcDisk.find(idx); + return (found != hitmcDisk.end()) ? found->second.first : -1; + }; + + // 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 disk = diskOf(idx, hitMCs[idx]); + if (disk < 0) continue; + for (const auto& edep : hitMCs[idx].energyDeposits()) + { + 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; 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, int>, art::Ptr> entrantCache; + size_t nUnresolved(0); + + for (size_t idx = 0; idx < hitMCs.size(); ++idx) + { + const auto& chmc = hitMCs[idx]; + std::vector> entrants; + + const int disk = diskOf(idx, chmc); + if (disk < 0) + { + // no crystalID member and not referenced by any cluster: + // entrants unresolvable; keep the collection index-parallel + // with null entries + entrants.assign(chmc.nParticles(), art::Ptr()); + ++nUnresolved; + output->emplace_back(art::Ptr(hitMCHandle, idx), std::move(entrants)); + continue; + } + entrants.reserve(chmc.nParticles()); + + for (const auto& edep : chmc.energyDeposits()) + { + const art::Ptr& sim = edep.sim(); + if (sim.isNull()) { entrants.emplace_back(); continue; } + + const auto key = std::make_pair(sim, disk); + 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()) + { + 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; + } + entrants.push_back(it->second); + } + 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" + << " (" << hitmcDisk.size() << " cluster-referenced, " + << nUnresolved << " unresolved)"; + } + + 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..845f3f71fc --- /dev/null +++ b/MCDataProducts/inc/CaloHitEntrant.hh @@ -0,0 +1,54 @@ +#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. +// +// 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. +// + +#include "canvas/Persistency/Common/Ptr.h" +#include "Offline/MCDataProducts/inc/CaloHitMC.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include +#include + +namespace mu2e +{ + class CaloHitEntrant + { + 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; +} + +#endif 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 + + + + +