From 72e857e8ce4eebb2b15a45d070f72bbda930c381 Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Wed, 28 Jan 2026 14:00:59 -0600 Subject: [PATCH 1/8] initial commit of generic filter --- sbndcode/Filters/CMakeLists.txt | 1 + sbndcode/Filters/TrueSignalFilter_module.cc | 303 ++++++++++++++++++++ 2 files changed, 304 insertions(+) create mode 100644 sbndcode/Filters/TrueSignalFilter_module.cc diff --git a/sbndcode/Filters/CMakeLists.txt b/sbndcode/Filters/CMakeLists.txt index 211a11e45..bb221f835 100644 --- a/sbndcode/Filters/CMakeLists.txt +++ b/sbndcode/Filters/CMakeLists.txt @@ -37,6 +37,7 @@ cet_build_plugin(LArG4CRTFilter art::module SOURCE LArG4CRTFilter_module.cc LIBR cet_build_plugin(LArG4FakeTriggerFilter art::module SOURCE LArG4FakeTriggerFilter_module.cc LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(SimEnergyDepFakeTriggerFilter art::module SOURCE SimEnergyDepFakeTriggerFilter_module.cc LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(NumberOfHitsFilter art::module SOURCE NumberOfHitsFilter_module.cc LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(TrueSignalFilter art::module SOURCE TrueSignalFilter_module.cc LIBRARIES ${MODULE_LIBRARIES}) install_headers() install_fhicl() diff --git a/sbndcode/Filters/TrueSignalFilter_module.cc b/sbndcode/Filters/TrueSignalFilter_module.cc new file mode 100644 index 000000000..f41694f3b --- /dev/null +++ b/sbndcode/Filters/TrueSignalFilter_module.cc @@ -0,0 +1,303 @@ +/* + * Filter module for common true signal definition options + * - True nu flavor + * - CC/NC + * - GENIE mode + * - Fiducial vertex + * - Final state primary PDG + * - Final state primary KE + * - Exclude PDG list + * - Exclusive? Exact match to primary list, or allow others (except those listed in exclude list) + * + * Can pass multiple filter options with fcl parameters + * Example: Pass CC events with no final state proton OR NC events with final + * state proton (not sure why you would want to do this, but it's possible!) + * { + * nu: 14 + * cc: true + * NoFinalStatePDG: [ 2112 ] + * } + * { + * nu: 14 + * cc: False + * FinalStatePDG: [ 2112 ] + * } + */ + +#include +#include +#include + +#include "TGeoManager.h" + +#include "fhiclcpp/types/Sequence.h" +#include "fhiclcpp/types/Tuple.h" +#include "fhiclcpp/types/OptionalAtom.h" +#include "fhiclcpp/types/OptionalSequence.h" +#include "art/Framework/Core/EDFilter.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "nusimdata/SimulationBase/MCTruth.h" +#include "nusimdata/SimulationBase/MCNeutrino.h" +#include "sbndcode/RecoUtils/RecoUtils.h" + +namespace filt{ + +// All optional. Empty block config will pass all events +struct FilterBlockConfig { + fhicl::OptionalAtom NuPDG { + fhicl::Name("NuPDG"), + fhicl::Comment("PDG code of the neutrino") + }; + + fhicl::OptionalAtom InTPC { + fhicl::Name("InTPC"), + fhicl::Comment("Require interaction vertex in the TPC") + }; + + + fhicl::OptionalAtom CCNC { + fhicl::Name("CCNC"), + fhicl::Comment("If true, only CC events are accepted. If false, only NC events are accepted. If ommitted, no requirement") + }; + + fhicl::OptionalSequence GenieModes { + fhicl::Name("GENIEModes"), + fhicl::Comment("List of GENIE interaction modes") + }; + + fhicl::OptionalSequence RequiredPDGs { + fhicl::Name("RequiredPDGs"), + fhicl::Comment("List of PDG codes that must appear in the event as primaries. May repeat PDG codes for multiple particles of the same type") + }; + + fhicl::OptionalSequence> KEThresholds { + fhicl::Name("KEThresholds"), + fhicl::Comment("Minimum KE for particles to count towards either required or disallowed particle lists. Format is [PDG code, threshold (MeV)]") + }; + + fhicl::OptionalAtom Exclusive { + fhicl::Name("Exclusive"), + fhicl::Comment("If true, event must contain the exact required PDGs with no extra primary particles") + }; + + fhicl::OptionalSequence DisallowedPDGs { + fhicl::Name("DisallowedPDGs"), + fhicl::Comment("List of PDG codes which must not be present in the primaries (not used if \"Exclusive\" option is true)") + }; +}; + + +struct FilterBlock { + std::optional nu_pdg; + std::optional in_tpc; + std::optional ccnc; + std::optional> genie_modes; + std::optional> required_pdgs; + std::optional> disallowed_pdgs; + std::optional> ke_thresholds; + std::optional exclusive; + + // computed from options + // (pdg, nrequired) + std::map required_counts; +}; + + +class TrueSignalFilter : public art::EDFilter { +public: + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + fhicl::Atom GENIELabel { + Name("GENIELabel"), Comment("Label for MC truth (GENIE)") + }; + fhicl::Sequence> Filters { + Name("Filters"), + Comment("List of filter-condition blocks; event passes if it matches any block") + }; + }; + using Parameters = art::EDFilter::Table; + explicit TrueSignalFilter(const Parameters& config); + + virtual bool filter(art::Event& e) override; + +protected: + bool PassBlock(const art::Ptr, const FilterBlock&) const; + +private: + art::InputTag fGenieModuleLabel; + std::vector fFilterBlocks; +}; + + +TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : + EDFilter{pset}, fGenieModuleLabel(pset().GENIELabel()) +{ + // loop over filter blocks + auto const& cfg_blocks = pset().Filters(); + fFilterBlocks.reserve(cfg_blocks.size()); + + for (std::size_t i = 0; i < cfg_blocks.size(); ++i) { + auto const& cfg = cfg_blocks.at(i); + FilterBlock block; + + int nu_pdg; + if (cfg.NuPDG(nu_pdg)) { + block.nu_pdg = nu_pdg; + } + + bool in_tpc; + if (cfg.InTPC(in_tpc)) { + block.in_tpc = in_tpc; + } + + bool ccnc; + if (cfg.CCNC(ccnc)) { + block.ccnc = ccnc; + } + + std::vector genie_modes; + if (cfg.GenieModes(genie_modes)) { + block.genie_modes = genie_modes; + } + + std::vector required_pdgs; + if (cfg.RequiredPDGs(required_pdgs)) { + block.required_pdgs = required_pdgs; + + // Count multiplicities in the requirements. + for (int pdg : required_pdgs) { + block.required_counts[pdg]++; + } + } + + std::vector disallowed_pdgs; + if (cfg.RequiredPDGs(disallowed_pdgs)) { + block.disallowed_pdgs = disallowed_pdgs; + } + + std::vector> ke_thresholds; + if (cfg.KEThresholds(ke_thresholds)) { + for (auto& pdg_threshold : ke_thresholds) { + std::map ke_thresholds; + int pdg = std::get<0>(pdg_threshold); + float threshold = std::get<1>(pdg_threshold); + auto [it, inserted] = ke_thresholds.emplace(pdg, threshold); + + // error on duplicates + if (!inserted) { + throw cet::exception("TrueSignalFilter") + << "In Filters[" << i << "]: KEThreshold duplicate " + << "entry for PDG " << pdg << ".\n"; + } + + block.ke_thresholds = ke_thresholds; + } + } + + bool exclusive; + if (cfg.Exclusive(exclusive)) { + block.exclusive = exclusive; + } + + fFilterBlocks.push_back(std::move(block)); + } +} + + +bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterBlock& block) const { + if (mc->Origin() != simb::kBeamNeutrino) return false; + + const auto& nu = mc->GetNeutrino(); + + if (block.nu_pdg) { + if (nu.Nu().PdgCode() != block.nu_pdg.value()) return false; + } + + if (block.in_tpc) { + // technically user could request vertices outside the tpc? + TVector3 vtx{ nu.Nu().Vx(), nu.Nu().Vy(), nu.Nu().Vz() }; + if (RecoUtils::IsInsideTPC(vtx, 0) != block.in_tpc.value()) return false; + } + + if (block.ccnc) { + if (block.ccnc.value() && (nu.CCNC() != simb::kCC)) return false; + if (!block.ccnc.value() && (nu.CCNC() != simb::kNC)) return false; + } + + //get a vector of final state particles for the next few checks + std::vector primaries; + for (int i = 0; i < mc->NParticles(); ++i) { + const auto& part(mc->GetParticle(i)); + // only consider primaries + if (part.StatusCode() != 1) continue; + if (part.Mother() > 0) continue; + + // check against threshold map. Primaries not counted if below threshold + if (block.ke_thresholds) { + int pdg = part.PdgCode(); + if (block.ke_thresholds->find(pdg) != block.ke_thresholds->end()) { + float ke = part.E() - part.Mass(); + if (ke < block.ke_thresholds->at(pdg)) continue; + } + } + primaries.push_back(part.PdgCode()); + } + + if (block.required_pdgs) { + // check that primaries contains all the required PDGs + // count multiplicities in the event + std::map counts; + for (int pdg : primaries) { + counts[pdg]++; + } + + // Check that for every required PDG, the event has enough. + // if exclusive: check equality + bool exclusive = block.exclusive.value_or(false); + for (auto const& [required_pdg, nrequired] : block.required_counts) { + auto it = counts.find(required_pdg); + if (it == counts.end() || it->second < nrequired) return false; + if (exclusive && it->second != nrequired) return false; + } + + // TODO check if there are extra particles not in the required list if exclusive + } + + if (block.disallowed_pdgs) { + // check that the primaries list doesn't contain anything disallowed + // note: primaries still are only counted if they are above threshold + for (int disallowed_pdg : *block.disallowed_pdgs) { + if (std::find(primaries.begin(), primaries.end(), disallowed_pdg) + != primaries.end()) return false; + } + } + + return true; +} + + +bool TrueSignalFilter::filter(art::Event & e) { + auto mclists = e.getMany>(); + for (size_t i = 0; i != mclists.size(); i++) { + const auto& mclist(mclists.at(i)); + for (size_t j = 0; j != mclist->size(); j++) { + art::Ptr mc(mclist, j); + for (auto const& block : fFilterBlocks) { + if (PassBlock(mc, block)) { + // mctruth passes at least one filter, so we are done + return true; + } + } + } + } + + // did not pass any filters + return false; +} + + +DEFINE_ART_MODULE(TrueSignalFilter) + +} From 9b76bd79f91d58333c4797aed95f69afabf3750a Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Wed, 28 Jan 2026 17:25:56 -0600 Subject: [PATCH 2/8] add example fcl and some logging --- sbndcode/Filters/TrueSignalFilter_module.cc | 105 +++++++++++++----- sbndcode/Filters/fcls/signal_filters_sbnd.fcl | 60 ++++++++++ 2 files changed, 137 insertions(+), 28 deletions(-) create mode 100644 sbndcode/Filters/fcls/signal_filters_sbnd.fcl diff --git a/sbndcode/Filters/TrueSignalFilter_module.cc b/sbndcode/Filters/TrueSignalFilter_module.cc index f41694f3b..6ea2ce9f5 100644 --- a/sbndcode/Filters/TrueSignalFilter_module.cc +++ b/sbndcode/Filters/TrueSignalFilter_module.cc @@ -1,8 +1,8 @@ /* * Filter module for common true signal definition options - * - True nu flavor + * - True nu flavors * - CC/NC - * - GENIE mode + * - Modes * - Fiducial vertex * - Final state primary PDG * - Final state primary KE @@ -41,13 +41,16 @@ #include "nusimdata/SimulationBase/MCNeutrino.h" #include "sbndcode/RecoUtils/RecoUtils.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + + namespace filt{ // All optional. Empty block config will pass all events struct FilterBlockConfig { - fhicl::OptionalAtom NuPDG { - fhicl::Name("NuPDG"), - fhicl::Comment("PDG code of the neutrino") + fhicl::OptionalSequence NuPDGs { + fhicl::Name("NuPDGs"), + fhicl::Comment("PDG codes of the neutrino") }; fhicl::OptionalAtom InTPC { @@ -56,14 +59,14 @@ struct FilterBlockConfig { }; - fhicl::OptionalAtom CCNC { - fhicl::Name("CCNC"), + fhicl::OptionalAtom IsCC { + fhicl::Name("IsCC"), fhicl::Comment("If true, only CC events are accepted. If false, only NC events are accepted. If ommitted, no requirement") }; - fhicl::OptionalSequence GenieModes { - fhicl::Name("GENIEModes"), - fhicl::Comment("List of GENIE interaction modes") + fhicl::OptionalSequence Modes { + fhicl::Name("Modes"), + fhicl::Comment("List of interaction modes") }; fhicl::OptionalSequence RequiredPDGs { @@ -89,10 +92,10 @@ struct FilterBlockConfig { struct FilterBlock { - std::optional nu_pdg; + std::optional> nu_pdgs; std::optional in_tpc; - std::optional ccnc; - std::optional> genie_modes; + std::optional iscc; + std::optional> modes; std::optional> required_pdgs; std::optional> disallowed_pdgs; std::optional> ke_thresholds; @@ -124,6 +127,7 @@ class TrueSignalFilter : public art::EDFilter { protected: bool PassBlock(const art::Ptr, const FilterBlock&) const; + void PrintBlock(const FilterBlock&) const; private: art::InputTag fGenieModuleLabel; @@ -142,9 +146,9 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : auto const& cfg = cfg_blocks.at(i); FilterBlock block; - int nu_pdg; - if (cfg.NuPDG(nu_pdg)) { - block.nu_pdg = nu_pdg; + std::vector nu_pdgs; + if (cfg.NuPDGs(nu_pdgs)) { + block.nu_pdgs = nu_pdgs; } bool in_tpc; @@ -152,14 +156,14 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : block.in_tpc = in_tpc; } - bool ccnc; - if (cfg.CCNC(ccnc)) { - block.ccnc = ccnc; + bool iscc; + if (cfg.IsCC(iscc)) { + block.iscc = iscc; } - std::vector genie_modes; - if (cfg.GenieModes(genie_modes)) { - block.genie_modes = genie_modes; + std::vector modes; + if (cfg.Modes(modes)) { + block.modes = modes; } std::vector required_pdgs; @@ -173,7 +177,7 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : } std::vector disallowed_pdgs; - if (cfg.RequiredPDGs(disallowed_pdgs)) { + if (cfg.DisallowedPDGs(disallowed_pdgs)) { block.disallowed_pdgs = disallowed_pdgs; } @@ -201,6 +205,7 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : block.exclusive = exclusive; } + PrintBlock(block); fFilterBlocks.push_back(std::move(block)); } } @@ -211,8 +216,14 @@ bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterB const auto& nu = mc->GetNeutrino(); - if (block.nu_pdg) { - if (nu.Nu().PdgCode() != block.nu_pdg.value()) return false; + if (block.nu_pdgs) { + if (std::find(block.nu_pdgs->begin(), block.nu_pdgs->end(), nu.Nu().PdgCode()) + == block.nu_pdgs->end()) return false; + } + + if (block.modes) { + if (std::find(block.modes->begin(), block.modes->end(), nu.Mode()) + == block.modes->end()) return false; } if (block.in_tpc) { @@ -221,9 +232,9 @@ bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterB if (RecoUtils::IsInsideTPC(vtx, 0) != block.in_tpc.value()) return false; } - if (block.ccnc) { - if (block.ccnc.value() && (nu.CCNC() != simb::kCC)) return false; - if (!block.ccnc.value() && (nu.CCNC() != simb::kNC)) return false; + if (block.iscc) { + if (block.iscc.value() && (nu.CCNC() != simb::kCC)) return false; + if (!block.iscc.value() && (nu.CCNC() != simb::kNC)) return false; } //get a vector of final state particles for the next few checks @@ -298,6 +309,44 @@ bool TrueSignalFilter::filter(art::Event & e) { } +void TrueSignalFilter::PrintBlock(const FilterBlock& block) const { + const std::string kNotSet(""); + + // print helpers + mf::LogInfo log{"TrueSignalFilter"}; + log << "Filter configuration:\n"; + auto print_int_vec = [&](const std::string& name, std::optional> const& val) { + log << " - " << name << ": "; + if (!val) { + log << kNotSet << "\n"; + } + else { + log << "{ "; + for (int i : val.value()) { + log << i << ", "; + } + log << " }\n"; + } + }; + auto print_bool = [&](const std::string& name, std::optional const& val) { + log << " - " << name << ": "; + if (!val) { + log << kNotSet << "\n"; + } + else { + log << (val.value() ? "True" : "False") << "\n"; + } + }; + + print_int_vec("NuPDGs", block.nu_pdgs); + print_bool("InTPC", block.in_tpc); + print_bool("IsCC", block.iscc); + print_int_vec("Modes", block.modes); + print_int_vec("RequiredPDGs", block.required_pdgs); + print_int_vec("DisallowedPDGs", block.disallowed_pdgs); + print_bool("Exclusive", block.exclusive); +} + DEFINE_ART_MODULE(TrueSignalFilter) } diff --git a/sbndcode/Filters/fcls/signal_filters_sbnd.fcl b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl new file mode 100644 index 000000000..e0541b579 --- /dev/null +++ b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl @@ -0,0 +1,60 @@ +BEGIN_PROLOG + +sbnd_ccnue_filter: +{ + module_type: TrueSignalFilter + GENIELabel: "generator" + Filters: [ + { + NuPDGs: [ 12, -12 ] + IsCC: true + InTPC: true + } + ] +} + +sbnd_cohpi_filter: +{ + module_type: TrueSignalFilter + GENIELabel: "generator" + Filters: [ + { + KEThresholds: [ + [ 13, 27 ], + [ 211, 30 ], + [ 2212, 20 ] + ] + RequiredPDGs: [ 13, 211 ] + DisallowedPDGs: [ 111, 2212, 2112 ] + InTPC: true + Exclusive: true + } + ] +} + +sbnd_eta_filter: +{ + module_type: TrueSignalFilter + GENIELabel: "generator" + Filters: [ + { + RequiredPDGs: [ 221 ] + DisallowedPDGs: [ 111 ] + InTPC: true + } + ] +} + +sbnd_pi0_filter: +{ + module_type: TrueSignalFilter + GENIELabel: "generator" + Filters: [ + { + RequiredPDGs: [ 111 ] + InTPC: true + } + ] +} + +END_PROLOG From 9b2154d2217c02afedae507bc1546a78bacec88f Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Thu, 29 Jan 2026 10:59:37 -0600 Subject: [PATCH 3/8] support more signal definitions --- sbndcode/Filters/fcls/signal_filters_sbnd.fcl | 34 +++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/sbndcode/Filters/fcls/signal_filters_sbnd.fcl b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl index e0541b579..54fb206fe 100644 --- a/sbndcode/Filters/fcls/signal_filters_sbnd.fcl +++ b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl @@ -32,6 +32,40 @@ sbnd_cohpi_filter: ] } +sbnd_2p2h_filter: +{ + module_type: TrueSignalFilter + GENIELabel: "generator" + Filters: [ + { + NuPDGs: [ 14 ] + IsCC: true + KEThresholds: [ + [ 13, 27 ], + [ 2212, 50 ] + ] + RequiredPDGs: [ 13, 2212, 2212 ] + DisallowedPDGs: [ 111, 211 ] + InTPC: true + Exclusive: true + } + ] +} + +sbnd_nc1p_filter: +{ + module_type: TrueSignalFilter + GENIELabel: "generator" + Filters: [ + { + IsCC: false + RequiredPDGs: [ 2212 ] + InTPC: true + Exclusive: true + } + ] +} + sbnd_eta_filter: { module_type: TrueSignalFilter From 205706f273d0d93d212209f95088613868b68c95 Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Thu, 29 Jan 2026 11:03:19 -0600 Subject: [PATCH 4/8] module support --- sbndcode/Filters/TrueSignalFilter_module.cc | 63 ++++++++++++++++++++- 1 file changed, 62 insertions(+), 1 deletion(-) diff --git a/sbndcode/Filters/TrueSignalFilter_module.cc b/sbndcode/Filters/TrueSignalFilter_module.cc index 6ea2ce9f5..9e0b40d1d 100644 --- a/sbndcode/Filters/TrueSignalFilter_module.cc +++ b/sbndcode/Filters/TrueSignalFilter_module.cc @@ -27,6 +27,7 @@ #include #include #include +#include #include "TGeoManager.h" @@ -53,6 +54,17 @@ struct FilterBlockConfig { fhicl::Comment("PDG codes of the neutrino") }; + fhicl::OptionalSequence TargetPDGs { + fhicl::Name("TargetPDGs"), + fhicl::Comment("Allowed target PDG codes") + }; + + fhicl::OptionalSequence WRange { + fhicl::Name("WRange"), + fhicl::Comment("Range of allowed invariant hadronic mass (W) as (min, max). Use -1 for min or max to indicate one-sided bound.") + }; + + fhicl::OptionalAtom InTPC { fhicl::Name("InTPC"), fhicl::Comment("Require interaction vertex in the TPC") @@ -66,7 +78,7 @@ struct FilterBlockConfig { fhicl::OptionalSequence Modes { fhicl::Name("Modes"), - fhicl::Comment("List of interaction modes") + fhicl::Comment("List of allowed interaction modes") }; fhicl::OptionalSequence RequiredPDGs { @@ -93,6 +105,8 @@ struct FilterBlockConfig { struct FilterBlock { std::optional> nu_pdgs; + std::optional> target_pdgs; + std::optional> wrange; std::optional in_tpc; std::optional iscc; std::optional> modes; @@ -104,6 +118,8 @@ struct FilterBlock { // computed from options // (pdg, nrequired) std::map required_counts; + float wmin = -1.0; + float wmax = -1.0; }; @@ -151,6 +167,25 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : block.nu_pdgs = nu_pdgs; } + std::vector target_pdgs; + if (cfg.TargetPDGs(target_pdgs)) { + block.target_pdgs = target_pdgs; + } + + std::array wrange; + if (cfg.WRange(wrange)) { + block.wrange = wrange; + block.wmin = wrange.at(0); + block.wmax = wrange.at(1); + + if (block.wmin > block.wmax && block.wmax > 0) { + throw cet::exception("TrueSignalFilter") + << "In Filters[" << i << "]: WRange (" + << block.wmin << ", " << block.wmax << ") invalid\n"; + } + + } + bool in_tpc; if (cfg.InTPC(in_tpc)) { block.in_tpc = in_tpc; @@ -221,6 +256,17 @@ bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterB == block.nu_pdgs->end()) return false; } + if (block.target_pdgs) { + if (std::find(block.target_pdgs->begin(), block.target_pdgs->end(), nu.Target()) + == block.target_pdgs->end()) return false; + } + + if (block.wrange) { + float w = nu.W(); + if (block.wmin >= 0. && w < block.wmin) return false; + if (block.wmax >= 0. && w > block.wmax) return false; + } + if (block.modes) { if (std::find(block.modes->begin(), block.modes->end(), nu.Mode()) == block.modes->end()) return false; @@ -337,14 +383,29 @@ void TrueSignalFilter::PrintBlock(const FilterBlock& block) const { log << (val.value() ? "True" : "False") << "\n"; } }; + auto print_ke_thresholds = [&](const std::string& name, const std::optional>& val) { + log << " - " << name << ": "; + if (!val) { + log << kNotSet << "\n"; + } + else { + log << "\n"; + for (const auto& t : val.value()) { + log << " - " << std::get<0>(t) << ": " << std::get<1>(t) << " MeV\n"; + } + } + }; print_int_vec("NuPDGs", block.nu_pdgs); + print_int_vec("TargetPDGs", block.target_pdgs); print_bool("InTPC", block.in_tpc); print_bool("IsCC", block.iscc); print_int_vec("Modes", block.modes); print_int_vec("RequiredPDGs", block.required_pdgs); print_int_vec("DisallowedPDGs", block.disallowed_pdgs); print_bool("Exclusive", block.exclusive); + print_ke_thresholds("KEThresholds", block.ke_thresholds); + // TODO WRange } DEFINE_ART_MODULE(TrueSignalFilter) From bafe41adef98195a0880607d6bd19a6c52613b1b Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Sat, 31 Jan 2026 09:34:19 -0600 Subject: [PATCH 5/8] add debug logging --- sbndcode/Filters/TrueSignalFilter_module.cc | 152 ++++++++++++------ sbndcode/Filters/fcls/signal_filters_sbnd.fcl | 6 +- 2 files changed, 108 insertions(+), 50 deletions(-) diff --git a/sbndcode/Filters/TrueSignalFilter_module.cc b/sbndcode/Filters/TrueSignalFilter_module.cc index 9e0b40d1d..331cce4b3 100644 --- a/sbndcode/Filters/TrueSignalFilter_module.cc +++ b/sbndcode/Filters/TrueSignalFilter_module.cc @@ -15,12 +15,12 @@ * { * nu: 14 * cc: true - * NoFinalStatePDG: [ 2112 ] + * NoFinalStatePDG: [ 2112 ] * } * { * nu: 14 * cc: False - * FinalStatePDG: [ 2112 ] + * FinalStatePDG: [ 2112 ] * } */ @@ -33,11 +33,11 @@ #include "fhiclcpp/types/Sequence.h" #include "fhiclcpp/types/Tuple.h" -#include "fhiclcpp/types/OptionalAtom.h" -#include "fhiclcpp/types/OptionalSequence.h" -#include "art/Framework/Core/EDFilter.h" -#include "art/Framework/Core/ModuleMacros.h" -#include "art/Framework/Principal/Event.h" +#include "fhiclcpp/types/OptionalAtom.h" +#include "fhiclcpp/types/OptionalSequence.h" +#include "art/Framework/Core/EDFilter.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" #include "nusimdata/SimulationBase/MCTruth.h" #include "nusimdata/SimulationBase/MCNeutrino.h" #include "sbndcode/RecoUtils/RecoUtils.h" @@ -64,13 +64,11 @@ struct FilterBlockConfig { fhicl::Comment("Range of allowed invariant hadronic mass (W) as (min, max). Use -1 for min or max to indicate one-sided bound.") }; - fhicl::OptionalAtom InTPC { fhicl::Name("InTPC"), fhicl::Comment("Require interaction vertex in the TPC") }; - fhicl::OptionalAtom IsCC { fhicl::Name("IsCC"), fhicl::Comment("If true, only CC events are accepted. If false, only NC events are accepted. If ommitted, no requirement") @@ -103,18 +101,18 @@ struct FilterBlockConfig { }; -struct FilterBlock { +struct FilterBlock { std::optional> nu_pdgs; std::optional> target_pdgs; std::optional> wrange; std::optional in_tpc; - std::optional iscc; + std::optional iscc; std::optional> modes; std::optional> required_pdgs; std::optional> disallowed_pdgs; std::optional> ke_thresholds; std::optional exclusive; - + // computed from options // (pdg, nrequired) std::map required_counts; @@ -142,17 +140,17 @@ class TrueSignalFilter : public art::EDFilter { virtual bool filter(art::Event& e) override; protected: - bool PassBlock(const art::Ptr, const FilterBlock&) const; + bool PassBlock(const art::Ptr, const FilterBlock&, mf::LogDebug&) const; void PrintBlock(const FilterBlock&) const; private: - art::InputTag fGenieModuleLabel; + art::InputTag fGenieModuleLabel; std::vector fFilterBlocks; }; TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : - EDFilter{pset}, fGenieModuleLabel(pset().GENIELabel()) + EDFilter{pset}, fGenieModuleLabel{pset().GENIELabel()} { // loop over filter blocks auto const& cfg_blocks = pset().Filters(); @@ -183,25 +181,25 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : << "In Filters[" << i << "]: WRange (" << block.wmin << ", " << block.wmax << ") invalid\n"; } - + } - bool in_tpc; + bool in_tpc; if (cfg.InTPC(in_tpc)) { block.in_tpc = in_tpc; } - bool iscc; + bool iscc; if (cfg.IsCC(iscc)) { block.iscc = iscc; } - std::vector modes; + std::vector modes; if (cfg.Modes(modes)) { block.modes = modes; } - std::vector required_pdgs; + std::vector required_pdgs; if (cfg.RequiredPDGs(required_pdgs)) { block.required_pdgs = required_pdgs; @@ -211,18 +209,18 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : } } - std::vector disallowed_pdgs; + std::vector disallowed_pdgs; if (cfg.DisallowedPDGs(disallowed_pdgs)) { block.disallowed_pdgs = disallowed_pdgs; } - std::vector> ke_thresholds; + std::vector> ke_thresholds; if (cfg.KEThresholds(ke_thresholds)) { - for (auto& pdg_threshold : ke_thresholds) { - std::map ke_thresholds; + std::map ke_threshold_map; + for (const auto& pdg_threshold : ke_thresholds) { int pdg = std::get<0>(pdg_threshold); float threshold = std::get<1>(pdg_threshold); - auto [it, inserted] = ke_thresholds.emplace(pdg, threshold); + auto [it, inserted] = ke_threshold_map.emplace(pdg, threshold); // error on duplicates if (!inserted) { @@ -230,12 +228,11 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : << "In Filters[" << i << "]: KEThreshold duplicate " << "entry for PDG " << pdg << ".\n"; } - - block.ke_thresholds = ke_thresholds; } + block.ke_thresholds = ke_threshold_map; } - bool exclusive; + bool exclusive; if (cfg.Exclusive(exclusive)) { block.exclusive = exclusive; } @@ -246,56 +243,74 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : } -bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterBlock& block) const { +bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterBlock& block, mf::LogDebug& debug_log) const { + const std::string kSpace = " "; + debug_log << kSpace << "Checking MCTruth is neutrino... "; if (mc->Origin() != simb::kBeamNeutrino) return false; + debug_log << "PASS\n"; const auto& nu = mc->GetNeutrino(); if (block.nu_pdgs) { - if (std::find(block.nu_pdgs->begin(), block.nu_pdgs->end(), nu.Nu().PdgCode()) + int nu_pdg = nu.Nu().PdgCode(); + debug_log << kSpace << "Checking MCTruth has appropriate neutrino PDG (Nu PDG=" << nu_pdg << ")... "; + if (std::find(block.nu_pdgs->begin(), block.nu_pdgs->end(), nu_pdg) == block.nu_pdgs->end()) return false; + debug_log << "PASS\n"; } if (block.target_pdgs) { - if (std::find(block.target_pdgs->begin(), block.target_pdgs->end(), nu.Target()) + int target_pdg = nu.Target(); + debug_log << kSpace << "Checking MCTruth has appropriate target PDG (Target PDG=" << target_pdg << ")... "; + if (std::find(block.target_pdgs->begin(), block.target_pdgs->end(), target_pdg) == block.target_pdgs->end()) return false; + debug_log << "PASS\n"; } if (block.wrange) { float w = nu.W(); + debug_log << kSpace << "Checking MCTruth has W within range (W=" << w << ")... "; if (block.wmin >= 0. && w < block.wmin) return false; if (block.wmax >= 0. && w > block.wmax) return false; + debug_log << "PASS\n"; } if (block.modes) { - if (std::find(block.modes->begin(), block.modes->end(), nu.Mode()) + int mode = nu.Mode(); + debug_log << kSpace << "Checking MCTruth has appropriate mode (mode=" << mode << ")... "; + if (std::find(block.modes->begin(), block.modes->end(), mode) == block.modes->end()) return false; - } - - if (block.in_tpc) { - // technically user could request vertices outside the tpc? - TVector3 vtx{ nu.Nu().Vx(), nu.Nu().Vy(), nu.Nu().Vz() }; - if (RecoUtils::IsInsideTPC(vtx, 0) != block.in_tpc.value()) return false; + debug_log << "PASS\n"; } if (block.iscc) { + debug_log << kSpace << "Checking MCTruth CC/NC (CCNC=" << nu.CCNC() << ")... "; if (block.iscc.value() && (nu.CCNC() != simb::kCC)) return false; if (!block.iscc.value() && (nu.CCNC() != simb::kNC)) return false; + debug_log << "PASS\n"; + } + + if (block.in_tpc) { + TVector3 vtx{ nu.Nu().Vx(), nu.Nu().Vy(), nu.Nu().Vz() }; + debug_log << kSpace << "Checking MCTruth is inside the TPC (x=" << vtx(0) << ", y=" << vtx(1) << ", z=" << vtx(2) << ")... "; + // technically user could request vertices outside the tpc? + if (RecoUtils::IsInsideTPC(vtx, 0) != block.in_tpc.value()) return false; + debug_log << "PASS\n"; } //get a vector of final state particles for the next few checks std::vector primaries; + bool has_thresholds = block.ke_thresholds.has_value(); for (int i = 0; i < mc->NParticles(); ++i) { const auto& part(mc->GetParticle(i)); // only consider primaries if (part.StatusCode() != 1) continue; - if (part.Mother() > 0) continue; // check against threshold map. Primaries not counted if below threshold - if (block.ke_thresholds) { + if (has_thresholds) { int pdg = part.PdgCode(); if (block.ke_thresholds->find(pdg) != block.ke_thresholds->end()) { - float ke = part.E() - part.Mass(); + float ke = 1.0e3 * (part.E() - part.Mass()); if (ke < block.ke_thresholds->at(pdg)) continue; } } @@ -303,6 +318,7 @@ bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterB } if (block.required_pdgs) { + debug_log << kSpace << "Checking list of primaries ["; // check that primaries contains all the required PDGs // count multiplicities in the event std::map counts; @@ -310,25 +326,47 @@ bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterB counts[pdg]++; } + for (const auto& [k, v] : counts) { + float threshold = 0.; + if (has_thresholds) { + if (block.ke_thresholds->find(k) != block.ke_thresholds->end()) { + threshold = block.ke_thresholds->at(k); + } + } + debug_log << "\n" << kSpace << kSpace << "PDG=" << k + << ", count above " << threshold << " MeV threshold=" << v; + } + debug_log << "\n" << kSpace << "]... "; + // Check that for every required PDG, the event has enough. // if exclusive: check equality bool exclusive = block.exclusive.value_or(false); - for (auto const& [required_pdg, nrequired] : block.required_counts) { + for (const auto& [required_pdg, nrequired] : block.required_counts) { auto it = counts.find(required_pdg); if (it == counts.end() || it->second < nrequired) return false; if (exclusive && it->second != nrequired) return false; } - // TODO check if there are extra particles not in the required list if exclusive + // check if there are extra particles not in the required list if exclusive + if (exclusive) { + for (const auto& [primary_pdg, count] : counts) { + // extra particle not in the list + if (block.required_counts.find(primary_pdg) == block.required_counts.end()) return false; + } + } + + debug_log << "PASS\n"; } if (block.disallowed_pdgs) { + debug_log << kSpace << "Checking no primaries are disallowed... "; // check that the primaries list doesn't contain anything disallowed // note: primaries still are only counted if they are above threshold for (int disallowed_pdg : *block.disallowed_pdgs) { if (std::find(primaries.begin(), primaries.end(), disallowed_pdg) != primaries.end()) return false; } + debug_log << "PASS\n"; } return true; @@ -336,16 +374,30 @@ bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterB bool TrueSignalFilter::filter(art::Event & e) { + // pass debug log reference to PassBlock so that we can print messages to + // the same line without adding too many timestamps to the output logs, and + // we also only have to write "FAIL" once in the code. + // We get a logger message once for every event + mf::LogDebug debug_log{"TrueSignalFilter"}; + auto mclists = e.getMany>(); for (size_t i = 0; i != mclists.size(); i++) { const auto& mclist(mclists.at(i)); for (size_t j = 0; j != mclist->size(); j++) { art::Ptr mc(mclist, j); - for (auto const& block : fFilterBlocks) { - if (PassBlock(mc, block)) { + // for (auto const& block : fFilterBlocks) { + debug_log << "Filtering MCTruth " << j << "...\n"; + for (size_t ib = 0; ib != fFilterBlocks.size(); ib++) { + const auto& block = fFilterBlocks.at(ib); + if (PassBlock(mc, block, debug_log)) { // mctruth passes at least one filter, so we are done + debug_log << "MCTruth " << j << " passed filter " << ib << "!\n"; return true; } + else { + debug_log << "FAIL\n"; + debug_log << "MCTruth " << j << " did not pass filter " << ib << "\n"; + } } } } @@ -368,8 +420,10 @@ void TrueSignalFilter::PrintBlock(const FilterBlock& block) const { } else { log << "{ "; - for (int i : val.value()) { - log << i << ", "; + const auto& vec = val.value(); + for (size_t i = 0; i != vec.size(); i++) { + log << vec.at(i); + if (i < vec.size() - 1) log << ", "; } log << " }\n"; } @@ -390,8 +444,8 @@ void TrueSignalFilter::PrintBlock(const FilterBlock& block) const { } else { log << "\n"; - for (const auto& t : val.value()) { - log << " - " << std::get<0>(t) << ": " << std::get<1>(t) << " MeV\n"; + for (const auto& [k, v] : val.value()) { + log << " - " << k << ": " << v << " MeV\n"; } } }; diff --git a/sbndcode/Filters/fcls/signal_filters_sbnd.fcl b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl index 54fb206fe..f5bb9149b 100644 --- a/sbndcode/Filters/fcls/signal_filters_sbnd.fcl +++ b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl @@ -1,5 +1,9 @@ BEGIN_PROLOG +# Signal definitions for exclusive cross section analysis sample production +# Each block can contain multiple filters in case the signal definition cannot +# be expressed with a single set of filter options! + sbnd_ccnue_filter: { module_type: TrueSignalFilter @@ -47,7 +51,7 @@ sbnd_2p2h_filter: RequiredPDGs: [ 13, 2212, 2212 ] DisallowedPDGs: [ 111, 211 ] InTPC: true - Exclusive: true + Exclusive: false } ] } From cd617153db3de9bd38a73c27dd2bcbbdd83ffba1 Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Sat, 31 Jan 2026 11:52:08 -0600 Subject: [PATCH 6/8] add ignored list --- sbndcode/Filters/TrueSignalFilter_module.cc | 78 ++++++++++--------- sbndcode/Filters/fcls/signal_filters_sbnd.fcl | 13 +++- 2 files changed, 52 insertions(+), 39 deletions(-) diff --git a/sbndcode/Filters/TrueSignalFilter_module.cc b/sbndcode/Filters/TrueSignalFilter_module.cc index 331cce4b3..f16c20d48 100644 --- a/sbndcode/Filters/TrueSignalFilter_module.cc +++ b/sbndcode/Filters/TrueSignalFilter_module.cc @@ -1,27 +1,17 @@ /* * Filter module for common true signal definition options + * Each filter block can check * - True nu flavors * - CC/NC * - Modes * - Fiducial vertex * - Final state primary PDG - * - Final state primary KE - * - Exclude PDG list - * - Exclusive? Exact match to primary list, or allow others (except those listed in exclude list) + * - Final state primary KE thresholds + * - Disallowed PDG list + * - Exact counts? Exact match to primary list, or allow others (except those listed in exclude list or below threshold) * - * Can pass multiple filter options with fcl parameters - * Example: Pass CC events with no final state proton OR NC events with final - * state proton (not sure why you would want to do this, but it's possible!) - * { - * nu: 14 - * cc: true - * NoFinalStatePDG: [ 2112 ] - * } - * { - * nu: 14 - * cc: False - * FinalStatePDG: [ 2112 ] - * } + * Multiple filters blocks per event are supported. Event is kept if it passes + * any filter block */ #include @@ -84,19 +74,24 @@ struct FilterBlockConfig { fhicl::Comment("List of PDG codes that must appear in the event as primaries. May repeat PDG codes for multiple particles of the same type") }; - fhicl::OptionalSequence> KEThresholds { - fhicl::Name("KEThresholds"), - fhicl::Comment("Minimum KE for particles to count towards either required or disallowed particle lists. Format is [PDG code, threshold (MeV)]") + fhicl::OptionalAtom ExactCounts { + fhicl::Name("ExactCounts"), + fhicl::Comment("If true, event must contain the exact count of particles listed in the required PDG list with no extra primaries, except those in the ignored list or below KE threshold") }; - - fhicl::OptionalAtom Exclusive { - fhicl::Name("Exclusive"), - fhicl::Comment("If true, event must contain the exact required PDGs with no extra primary particles") + + fhicl::OptionalSequence IgnoredPDGs { + fhicl::Name("IgnoredPDGs"), + fhicl::Comment("List of PDG codes which are not considered when counting primaries, even if they are above KE thresholds") }; fhicl::OptionalSequence DisallowedPDGs { fhicl::Name("DisallowedPDGs"), - fhicl::Comment("List of PDG codes which must not be present in the primaries (not used if \"Exclusive\" option is true)") + fhicl::Comment("List of PDG codes which must not be present in the primaries") + }; + + fhicl::OptionalSequence> KEThresholds { + fhicl::Name("KEThresholds"), + fhicl::Comment("Minimum KE for particles to count towards either required or disallowed particle lists. Format is [PDG code, threshold (MeV)]") }; }; @@ -109,9 +104,10 @@ struct FilterBlock { std::optional iscc; std::optional> modes; std::optional> required_pdgs; + std::optional> ignored_pdgs; std::optional> disallowed_pdgs; std::optional> ke_thresholds; - std::optional exclusive; + std::optional exact_counts; // computed from options // (pdg, nrequired) @@ -232,9 +228,9 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : block.ke_thresholds = ke_threshold_map; } - bool exclusive; - if (cfg.Exclusive(exclusive)) { - block.exclusive = exclusive; + bool exact_counts; + if (cfg.ExactCounts(exact_counts)) { + block.exact_counts = exact_counts; } PrintBlock(block); @@ -301,14 +297,22 @@ bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterB //get a vector of final state particles for the next few checks std::vector primaries; bool has_thresholds = block.ke_thresholds.has_value(); + bool has_ignored = block.ignored_pdgs.has_value(); for (int i = 0; i < mc->NParticles(); ++i) { const auto& part(mc->GetParticle(i)); + debug_log << kSpace << "MC particle list " << i << " " << part.PdgCode() << " STATUS=" << part.StatusCode() << ", E=" << part.E() << "\n"; // only consider primaries if (part.StatusCode() != 1) continue; - // check against threshold map. Primaries not counted if below threshold + // don't count particles in the ignored list + int pdg = part.PdgCode(); + if (has_ignored) { + if (std::find(block.ignored_pdgs->begin(), block.ignored_pdgs->end(), pdg) + != block.ignored_pdgs->end()) continue; + } + + // don't count particles below threshold if (has_thresholds) { - int pdg = part.PdgCode(); if (block.ke_thresholds->find(pdg) != block.ke_thresholds->end()) { float ke = 1.0e3 * (part.E() - part.Mass()); if (ke < block.ke_thresholds->at(pdg)) continue; @@ -339,16 +343,17 @@ bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterB debug_log << "\n" << kSpace << "]... "; // Check that for every required PDG, the event has enough. - // if exclusive: check equality - bool exclusive = block.exclusive.value_or(false); + // if exact_counts: check equality + bool exact_counts = block.exact_counts.value_or(false); for (const auto& [required_pdg, nrequired] : block.required_counts) { auto it = counts.find(required_pdg); if (it == counts.end() || it->second < nrequired) return false; - if (exclusive && it->second != nrequired) return false; + if (exact_counts && it->second != nrequired) return false; } - // check if there are extra particles not in the required list if exclusive - if (exclusive) { + // check if there are extra particles not in the required list if exact_counts is set + // note: counts already excludes below threshold particles & ignored particles + if (exact_counts) { for (const auto& [primary_pdg, count] : counts) { // extra particle not in the list if (block.required_counts.find(primary_pdg) == block.required_counts.end()) return false; @@ -456,8 +461,9 @@ void TrueSignalFilter::PrintBlock(const FilterBlock& block) const { print_bool("IsCC", block.iscc); print_int_vec("Modes", block.modes); print_int_vec("RequiredPDGs", block.required_pdgs); + print_int_vec("IgnoredPDGs", block.ignored_pdgs); print_int_vec("DisallowedPDGs", block.disallowed_pdgs); - print_bool("Exclusive", block.exclusive); + print_bool("ExactCounts", block.exact_counts); print_ke_thresholds("KEThresholds", block.ke_thresholds); // TODO WRange } diff --git a/sbndcode/Filters/fcls/signal_filters_sbnd.fcl b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl index f5bb9149b..61b58d6c6 100644 --- a/sbndcode/Filters/fcls/signal_filters_sbnd.fcl +++ b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl @@ -4,6 +4,7 @@ BEGIN_PROLOG # Each block can contain multiple filters in case the signal definition cannot # be expressed with a single set of filter options! +# CC nue or nuebar sbnd_ccnue_filter: { module_type: TrueSignalFilter @@ -17,6 +18,7 @@ sbnd_ccnue_filter: ] } +# one muon > 27 MeV, one charged pion > 30 MeV, no protons > 20 MeV, no neutrons, no pi0 sbnd_cohpi_filter: { module_type: TrueSignalFilter @@ -31,11 +33,12 @@ sbnd_cohpi_filter: RequiredPDGs: [ 13, 211 ] DisallowedPDGs: [ 111, 2212, 2112 ] InTPC: true - Exclusive: true + ExactCounts: true } ] } +# one muon > 27 MeV, two protons > 50 MeV, no pions, any number of neutrons sbnd_2p2h_filter: { module_type: TrueSignalFilter @@ -49,13 +52,15 @@ sbnd_2p2h_filter: [ 2212, 50 ] ] RequiredPDGs: [ 13, 2212, 2212 ] + IgnoredPDGs: [ 2112 ] DisallowedPDGs: [ 111, 211 ] InTPC: true - Exclusive: false + ExactCounts: true } ] } +# one proton only! sbnd_nc1p_filter: { module_type: TrueSignalFilter @@ -65,11 +70,12 @@ sbnd_nc1p_filter: IsCC: false RequiredPDGs: [ 2212 ] InTPC: true - Exclusive: true + ExactCounts: true } ] } +# one eta + anything except pi0s sbnd_eta_filter: { module_type: TrueSignalFilter @@ -83,6 +89,7 @@ sbnd_eta_filter: ] } +# pi0 + anything sbnd_pi0_filter: { module_type: TrueSignalFilter From 8c730ec533bda962fa3d095961181e09901a87bf Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Mon, 2 Feb 2026 10:27:06 -0600 Subject: [PATCH 7/8] add k filter --- sbndcode/Filters/TrueSignalFilter_module.cc | 13 ++++++------ sbndcode/Filters/fcls/signal_filters_sbnd.fcl | 21 +++++++++++++++++-- 2 files changed, 25 insertions(+), 9 deletions(-) diff --git a/sbndcode/Filters/TrueSignalFilter_module.cc b/sbndcode/Filters/TrueSignalFilter_module.cc index f16c20d48..7142735b4 100644 --- a/sbndcode/Filters/TrueSignalFilter_module.cc +++ b/sbndcode/Filters/TrueSignalFilter_module.cc @@ -41,12 +41,12 @@ namespace filt{ struct FilterBlockConfig { fhicl::OptionalSequence NuPDGs { fhicl::Name("NuPDGs"), - fhicl::Comment("PDG codes of the neutrino") + fhicl::Comment("Allowed PDG codes for the neutrino") }; fhicl::OptionalSequence TargetPDGs { fhicl::Name("TargetPDGs"), - fhicl::Comment("Allowed target PDG codes") + fhicl::Comment("Allowed PDG codes for the target") }; fhicl::OptionalSequence WRange { @@ -71,7 +71,7 @@ struct FilterBlockConfig { fhicl::OptionalSequence RequiredPDGs { fhicl::Name("RequiredPDGs"), - fhicl::Comment("List of PDG codes that must appear in the event as primaries. May repeat PDG codes for multiple particles of the same type") + fhicl::Comment("List of PDG codes that must appear in the event as primaries. May repeat PDG codes to require multiple particles of the same type") }; fhicl::OptionalAtom ExactCounts { @@ -86,12 +86,12 @@ struct FilterBlockConfig { fhicl::OptionalSequence DisallowedPDGs { fhicl::Name("DisallowedPDGs"), - fhicl::Comment("List of PDG codes which must not be present in the primaries") + fhicl::Comment("List of PDG codes which must not be present in the list of primaries if they are above KE thresholds") }; fhicl::OptionalSequence> KEThresholds { fhicl::Name("KEThresholds"), - fhicl::Comment("Minimum KE for particles to count towards either required or disallowed particle lists. Format is [PDG code, threshold (MeV)]") + fhicl::Comment("Minimum KE required for particles to count towards the list of primaries. Format is [PDG code, threshold (MeV)]") }; }; @@ -239,7 +239,7 @@ TrueSignalFilter::TrueSignalFilter(const Parameters& pset) : } -bool TrueSignalFilter::PassBlock(const art::Ptr mc, const FilterBlock& block, mf::LogDebug& debug_log) const { +bool TrueSignalFilter::PassBlock(const art::Ptr& mc, const FilterBlock& block, mf::LogDebug& debug_log) const { const std::string kSpace = " "; debug_log << kSpace << "Checking MCTruth is neutrino... "; if (mc->Origin() != simb::kBeamNeutrino) return false; @@ -390,7 +390,6 @@ bool TrueSignalFilter::filter(art::Event & e) { const auto& mclist(mclists.at(i)); for (size_t j = 0; j != mclist->size(); j++) { art::Ptr mc(mclist, j); - // for (auto const& block : fFilterBlocks) { debug_log << "Filtering MCTruth " << j << "...\n"; for (size_t ib = 0; ib != fFilterBlocks.size(); ib++) { const auto& block = fFilterBlocks.at(ib); diff --git a/sbndcode/Filters/fcls/signal_filters_sbnd.fcl b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl index 61b58d6c6..f2ed7be0e 100644 --- a/sbndcode/Filters/fcls/signal_filters_sbnd.fcl +++ b/sbndcode/Filters/fcls/signal_filters_sbnd.fcl @@ -1,8 +1,8 @@ BEGIN_PROLOG # Signal definitions for exclusive cross section analysis sample production -# Each block can contain multiple filters in case the signal definition cannot -# be expressed with a single set of filter options! +# The "Filters" option accepts a list of parameter sets in case the signal +# definition cannot be expressed with a single set of filter options # CC nue or nuebar sbnd_ccnue_filter: @@ -102,4 +102,21 @@ sbnd_pi0_filter: ] } +# Charged or neutral kaon + anything +sbnd_kaon_filter: +{ + module_type: TrueSignalFilter + GENIELabel: "generator" + Filters: [ + { + RequiredPDGs: [ 321 ] + InTPC: true + }, + { + RequiredPDGs: [ 311 ] + InTPC: true + } + ] +} + END_PROLOG From ae8f468fcaecb68aade57c0741799d6d487dfefc Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Fri, 6 Feb 2026 17:04:18 -0600 Subject: [PATCH 8/8] fix ampersand --- sbndcode/Filters/TrueSignalFilter_module.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbndcode/Filters/TrueSignalFilter_module.cc b/sbndcode/Filters/TrueSignalFilter_module.cc index 7142735b4..4cbd4f074 100644 --- a/sbndcode/Filters/TrueSignalFilter_module.cc +++ b/sbndcode/Filters/TrueSignalFilter_module.cc @@ -136,7 +136,7 @@ class TrueSignalFilter : public art::EDFilter { virtual bool filter(art::Event& e) override; protected: - bool PassBlock(const art::Ptr, const FilterBlock&, mf::LogDebug&) const; + bool PassBlock(const art::Ptr&, const FilterBlock&, mf::LogDebug&) const; void PrintBlock(const FilterBlock&) const; private: