diff --git a/inc/CaloHitInfoMC.hh b/inc/CaloHitInfoMC.hh index ef75d558..534d414b 100644 --- a/inc/CaloHitInfoMC.hh +++ b/inc/CaloHitInfoMC.hh @@ -21,6 +21,7 @@ namespace mu2e std::vector momentumIns; // list of the momentum of the SimParticle when entering in the disk std::vector simParticleIds; // list of simparticle ids std::vector simRels; // relationship to the particle that deposited the most energy in the calo Hit + std::vector entrantSimIds; // calo-entrant (shower originator) SimParticle id per deposit, aligned with simParticleIds; filled only when calo.mc.entrantTag is configured int clusterIdx_; // Cluster index int caloHitIdx_; // index into calohits branch, -1 if unset; calohitsmc is NOT index-aligned with calohits diff --git a/inc/InfoMCStructHelper.hh b/inc/InfoMCStructHelper.hh index a9a0fcde..edbf12d9 100644 --- a/inc/InfoMCStructHelper.hh +++ b/inc/InfoMCStructHelper.hh @@ -19,6 +19,7 @@ #include "EventNtuple/inc/CaloClusterInfoMC.hh" #include "EventNtuple/inc/CaloHitInfoMC.hh" #include "EventNtuple/inc/CaloDigiMCInfo.hh" +#include "Offline/MCDataProducts/inc/CaloHitEntrant.hh" #include "EventNtuple/inc/MCStepInfo.hh" #include "EventNtuple/inc/MCStepSummaryInfo.hh" #include "EventNtuple/inc/SurfaceStepInfo.hh" @@ -75,7 +76,7 @@ namespace mu2e { void fillVDInfo(KalSeed const& kseed, const KalSeedMC& kseedmc, std::vector>& all_vdinfos); void fillHitInfoMCs(const KalSeed& kseed, const KalSeedMC& kseedmc, std::vector>& all_tshinfomcs); void fillCaloClusterInfoMC(CaloClusterMC const& ccmc, std::vector& ccimc); - void fillCaloHitInfoMC(CaloHitMC const& chmc, std::vector& chimc, int clusterIdx = -1); + void fillCaloHitInfoMC(CaloHitMC const& chmc, std::vector& chimc, int clusterIdx = -1, const CaloHitEntrant* entrant = nullptr); void fillCaloDigiMCInfo(CaloShowerSim const& shower, std::vector& calodigimc); void fillCaloDigiSimInfos(CaloShowerSim const& shower, std::vector& cdsis); void fillCaloSimInfos(CaloClusterMC const& ccmc, std::vector& csis); diff --git a/src/EventNtupleMaker_module.cc b/src/EventNtupleMaker_module.cc index 94991962..fcee2987 100644 --- a/src/EventNtupleMaker_module.cc +++ b/src/EventNtupleMaker_module.cc @@ -9,6 +9,7 @@ #include "Offline/MCDataProducts/inc/KalSeedMC.hh" #include "Offline/MCDataProducts/inc/CaloClusterMC.hh" #include "Offline/MCDataProducts/inc/CaloHitMC.hh" +#include "Offline/MCDataProducts/inc/CaloHitEntrant.hh" #include "Offline/MCDataProducts/inc/ProtonBunchTimeMC.hh" #include "Offline/RecoDataProducts/inc/KalSeed.hh" #include "Offline/RecoDataProducts/inc/KalSeedAssns.hh" @@ -220,6 +221,7 @@ namespace mu2e { fhicl::Atom fillDigis {Name("fillDigis"), Comment("Fill standalone calodigismc. branch")}; fhicl::Atom fillDigiSim {Name("fillDigiSim"), Comment("Fill calodigisim. branch")}; fhicl::Atom showerSimTag{Name("showerSimTag"), Comment("Tag for CaloShowerSim collection")}; + fhicl::Atom entrantTag {Name("entrantTag"), Comment("Tag for CaloHitEntrantCollection (calo-entrant truth in calohitsmc.entrantSimIds); empty disables"), art::InputTag()}; }; fhicl::Table mc{Name("mc"), Comment("Calorimeter MC filling options")}; }; @@ -366,6 +368,7 @@ namespace mu2e { std::map>> _allMCVDInfos; art::Handle _ccmcch; art::Handle _chmcch; + art::Handle _caloEntrants; std::map> _allMCTCHIs; // hit level info branches std::map>> _allTSHIs; @@ -982,8 +985,33 @@ namespace mu2e { } } if(fillCaloHitsMC()){ + // Optional calo-entrant truth: the CaloHitEntrantCollection is + // index-parallel to the CaloHitMCCollection by construction. + const CaloHitEntrantCollection* entrants = nullptr; + if(!_conf.calo().mc().entrantTag().empty()){ + event.getByLabel(_conf.calo().mc().entrantTag(),_caloEntrants); + if(!_caloEntrants.isValid()){ + throw cet::exception("EventNtuple") << "CaloHitEntrantCollection not found for entrantTag \"" << _conf.calo().mc().entrantTag().encode() << "\"\n"; + } + entrants = _caloEntrants.product(); + if(entrants->size() != _chmcch->size()){ + throw cet::exception("EventNtuple") << "CaloHitEntrantCollection size " << entrants->size() << " != CaloHitMCCollection size " << _chmcch->size() << "\n"; + } + // Source identity: each entrant carries a Ptr to the CaloHitMC it was + // computed for; require it to be entry i of the collection this fill + // iterates, so a product built from a different CaloHitMCCollection + // (or reordered) cannot be flattened silently. + for(size_t i = 0; i < entrants->size(); ++i){ + const auto& hmcPtr = (*entrants)[i].caloHitMC(); + if(hmcPtr.id() != _chmcch.id() || hmcPtr.key() != i){ + throw cet::exception("EventNtuple") << "CaloHitEntrant " << i << " references CaloHitMC (ProductID " << hmcPtr.id() << ", key " << hmcPtr.key() << ") but this fill iterates (ProductID " << _chmcch.id() << ", key " << i << ")\n"; + } + } + } + size_t entrantIdx = 0; for(const auto& hitmc : *_chmcch.product()){ - _infoMCStructHelper.fillCaloHitInfoMC(hitmc,_caloHIMCs); + _infoMCStructHelper.fillCaloHitInfoMC(hitmc,_caloHIMCs,-1,entrants ? &(*entrants)[entrantIdx] : nullptr); + ++entrantIdx; } } if(fillCaloClsMC()){ diff --git a/src/InfoMCStructHelper.cc b/src/InfoMCStructHelper.cc index a738373e..671d4e17 100644 --- a/src/InfoMCStructHelper.cc +++ b/src/InfoMCStructHelper.cc @@ -6,6 +6,7 @@ #include "Offline/MCDataProducts/inc/StepPointMC.hh" #include "Offline/MCDataProducts/inc/SimParticle.hh" #include "Offline/MCDataProducts/inc/MCRelationship.hh" +#include "cetlib_except/exception.h" #include "Offline/TrackerGeom/inc/Tracker.hh" #include "Offline/Mu2eUtilities/inc/TwoLinePCA.hh" @@ -381,9 +382,13 @@ namespace mu2e { ccimcs.push_back(ccimc); } - void InfoMCStructHelper::fillCaloHitInfoMC(CaloHitMC const& chmc, std::vector& chimcs, int clusterIdx) { + void InfoMCStructHelper::fillCaloHitInfoMC(CaloHitMC const& chmc, std::vector& chimcs, int clusterIdx, const CaloHitEntrant* entrant) { CaloHitInfoMC chimc; auto const& edeps = chmc.energyDeposits(); + if (entrant != nullptr && entrant->entrants().size() != edeps.size()) { + throw cet::exception("EventNtuple") << "CaloHitEntrant has " << entrant->entrants().size() + << " entrants but CaloHitMC has " << edeps.size() << " energy deposits\n"; + } chimc.crystalID_ = chmc.crystalID(); chimc.nsim = edeps.size(); chimc.eDep = chmc.totalEnergyDep(); @@ -393,6 +398,7 @@ namespace mu2e { if (chimc.nsim > 0){ chimc.eprimary = edeps.front().energyDep(); chimc.tprimary = edeps.front().time(); + size_t iedep = 0; for (auto const& edep : edeps){ auto simid = edep.sim()->id().asInt(); chimc.tDeps.push_back(edep.time()); @@ -400,6 +406,13 @@ namespace mu2e { chimc.momentumIns.push_back(edep.momentumIn()); chimc.simParticleIds.push_back(simid); chimc.simRels.push_back(MCRelationship(edep.sim(),edeps.front().sim())); + // entrantSimIds aligned with simParticleIds; -1 marks a deposit the + // upstream CaloEntrantTruthMaker could not resolve + if (entrant != nullptr) { + const auto& ep = entrant->entrants()[iedep]; + chimc.entrantSimIds.push_back(ep.isNonnull() ? ep->id().asInt() : -1); + } + ++iedep; } } chimcs.push_back(chimc);