From 851e176590514bd2a0e306d0b5896f4598fff85d Mon Sep 17 00:00:00 2001 From: Seokju Chung Date: Tue, 21 Apr 2026 09:07:23 -0500 Subject: [PATCH 1/3] Return false instead of exit --- sbncode/SinglePhotonAnalysis/NCDeltaRadiative_module.cc | 4 ++-- sbncode/SinglePhotonAnalysis/NCRadiativeResonant_module.cc | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/sbncode/SinglePhotonAnalysis/NCDeltaRadiative_module.cc b/sbncode/SinglePhotonAnalysis/NCDeltaRadiative_module.cc index aff90eb39..892d31f06 100644 --- a/sbncode/SinglePhotonAnalysis/NCDeltaRadiative_module.cc +++ b/sbncode/SinglePhotonAnalysis/NCDeltaRadiative_module.cc @@ -161,7 +161,7 @@ bool NCDeltaRadiative::filter(art::Event & e) { simb::MCParticle const & mcp = mct.GetParticle(i); if(mcp.TrackId() != i) { std::cout << "ERROR: " << __LINE__ << " " << __PRETTY_FUNCTION__ << "\nTrackId does not match index\n"; - exit(1); + return false; } if(!(mcp.StatusCode() == 1 && mcp.PdgCode() == 22)) continue; exiting_photon_parents.push_back(mcp.Mother()); @@ -173,7 +173,7 @@ bool NCDeltaRadiative::filter(art::Event & e) { //CHEKC hardcode, TPC filter: if(abs(mcp.Vx())>210 || abs(mcp.Vy())>210||mcp.Vz()>510 || mcp.Vz()<-1){ std::cout<<"OUTSIDE TPC x y z ="<210 || abs(mcp.Vy())>210||mcp.Vz()>510 || mcp.Vz()<-1){ std::cout<<"OUTSIDE TPC x y z ="< Date: Tue, 5 May 2026 09:59:58 -0500 Subject: [PATCH 2/3] Add CCDeltaRadiative filter and fcl --- sbncode/SinglePhotonAnalysis/CMakeLists.txt | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/sbncode/SinglePhotonAnalysis/CMakeLists.txt b/sbncode/SinglePhotonAnalysis/CMakeLists.txt index 35f24c2bb..4d144d72c 100644 --- a/sbncode/SinglePhotonAnalysis/CMakeLists.txt +++ b/sbncode/SinglePhotonAnalysis/CMakeLists.txt @@ -23,6 +23,13 @@ art_make_library( BASENAME_ONLY ROOT::Core ) +cet_build_plugin(CCDeltaRadiative art::module BASENAME_ONLY LIBRARIES + sbncode_SinglePhotonAnalysis_Libraries + sbncode_SinglePhotonAnalysis_SEAview + art_root_io::TFileService_service + ROOT::Tree + ) + cet_build_plugin(NCDeltaRadiative art::module BASENAME_ONLY LIBRARIES sbncode_SinglePhotonAnalysis_Libraries sbncode_SinglePhotonAnalysis_SEAview From 41d2c59e24211a5672b8a8ac0a05bb9a52b46a3c Mon Sep 17 00:00:00 2001 From: Seokju Chung Date: Tue, 5 May 2026 11:50:52 -0500 Subject: [PATCH 3/3] CC Delta Rad filter --- .../CCDeltaRadiative_module.cc | 212 ++++++++++++++++++ ...ion_tpc_CCDeltaRadiative_filtered_sbnd.fcl | 106 +++++++++ 2 files changed, 318 insertions(+) create mode 100644 sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc create mode 100644 sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl diff --git a/sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc b/sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc new file mode 100644 index 000000000..54ee42f59 --- /dev/null +++ b/sbncode/SinglePhotonAnalysis/CCDeltaRadiative_module.cc @@ -0,0 +1,212 @@ +//////////////////////////////////////////////////////////////////////// +// Class: CCDeltaRadiative +// Plugin Type: filter (art v2_05_00) +// File: CCDeltaRadiative_module.cc +// +// Selects CC events with a photon from a Delta radiative decay. +// Mirrors NCDeltaRadiative_module.cc but requires CCNC == 0 (charged +// current) and checks all four Delta charge states (1114, 2114, 2214, 2224). +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDFilter.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "art_root_io/TFileService.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art_root_io/TFileDirectory.h" + +#include + +#include "nusimdata/SimulationBase/MCTruth.h" + +#include "TTree.h" + +class CCDeltaRadiative : public art::EDFilter { + + TTree * ftree; + + int frun; + int fsubrun; + int fevent; + int fnu_pdg; + int fccnc; + int fmode; + int finteraction_type; + int fis_cc_delta_radiative; + int fparent_status_code; + int fparent_pdg; + +public: + explicit CCDeltaRadiative(fhicl::ParameterSet const & p); + + CCDeltaRadiative(CCDeltaRadiative const &) = delete; + CCDeltaRadiative(CCDeltaRadiative &&) = delete; + CCDeltaRadiative & operator = (CCDeltaRadiative const &) = delete; + CCDeltaRadiative & operator = (CCDeltaRadiative &&) = delete; + + void cout_stuff(art::Event & e, bool passed); + void FillTree(art::Event & e, + size_t const mct_index, + size_t const parent_index, + bool const is_cc_delta_radiative); + void Reset(); + bool filter(art::Event & e) override; + +}; + + +CCDeltaRadiative::CCDeltaRadiative(fhicl::ParameterSet const & p) : + art::EDFilter(p), + ftree(nullptr) { + + art::ServiceHandle tfs; + ftree = tfs->make("CCDeltaRadiativeFilter", ""); + + ftree->Branch("run", &frun, "run/I"); + ftree->Branch("subrun", &fsubrun, "subrun/I"); + ftree->Branch("event", &fevent, "event/I"); + ftree->Branch("nu_pdg", &fnu_pdg, "nu_pdg/I"); + ftree->Branch("ccnc", &fccnc, "ccnc/I"); + ftree->Branch("mode", &fmode, "mode/I"); + ftree->Branch("is_cc_delta_radiative", &fis_cc_delta_radiative, "is_cc_delta_radiative/I"); + ftree->Branch("parent_status_code", &fparent_status_code, "parent_status_code/I"); + ftree->Branch("parent_pdg", &fparent_pdg, "parent_pdg/I"); + +} + + +void CCDeltaRadiative::cout_stuff(art::Event & e, bool passed = false) { + + art::ValidHandle> const & ev_mct = + e.getValidHandle>("generator"); + + std::cout << passed << "\n" + << "==========================\n"; + for(simb::MCTruth const & mct : *ev_mct) { + std::cout << "----------------------------\n"; + for(int i = 0; i < mct.NParticles(); ++i) { + simb::MCParticle const & mcp = mct.GetParticle(i); + std::cout << mcp.TrackId() << " " << mcp.PdgCode() << " " << mcp.Mother() << " " << mcp.StatusCode() << "\n"; + } + } + +} + + +void CCDeltaRadiative::Reset() { + + frun = -1; + fsubrun = -1; + fevent = -1; + fnu_pdg = 0; + fccnc = -1; + fmode = -2; + finteraction_type = -2; + fis_cc_delta_radiative = -1; + fparent_status_code = -1; + fparent_pdg = 0; + +} + + +void CCDeltaRadiative::FillTree(art::Event & e, + size_t const mct_index, + size_t const parent_index, + bool const is_cc_delta_radiative) { + + Reset(); + + frun = e.id().run(); + fsubrun = e.id().subRun(); + fevent = e.id().event(); + + art::ValidHandle> const & ev_mct = + e.getValidHandle>("generator"); + simb::MCNeutrino const & mcn = ev_mct->at(mct_index).GetNeutrino(); + + fnu_pdg = mcn.Nu().PdgCode(); + fccnc = mcn.CCNC(); + fmode = mcn.Mode(); + finteraction_type = mcn.InteractionType(); + + fis_cc_delta_radiative = is_cc_delta_radiative; + if(parent_index != SIZE_MAX) { + fparent_status_code = ev_mct->at(mct_index).GetParticle(parent_index).StatusCode(); + fparent_pdg = ev_mct->at(mct_index).GetParticle(parent_index).PdgCode(); + } + + ftree->Fill(); + +} + + +// Returns true if |pdg| matches any Delta charge state (Δ-, Δ0, Δ+, Δ++). +static bool isDelta(int pdg) { + int apdg = std::abs(pdg); + return apdg == 1114 || apdg == 2114 || apdg == 2214 || apdg == 2224; +} + + +bool CCDeltaRadiative::filter(art::Event & e) { + + art::ValidHandle> const & ev_mct = + e.getValidHandle>("generator"); + + for(size_t i = 0; i < ev_mct->size(); ++i) { + + simb::MCTruth const & mct = ev_mct->at(i); + // Require charged-current interaction (CCNC == 0). + if(mct.GetNeutrino().CCNC() != 0) continue; + + std::vector exiting_photon_parents; + for(int j = 0; j < mct.NParticles(); ++j) { + simb::MCParticle const & mcp = mct.GetParticle(j); + if(mcp.TrackId() != j) { + std::cout << "ERROR: " << __LINE__ << " " << __PRETTY_FUNCTION__ << "\nTrackId does not match index\n"; + return false; + } + if(!(mcp.StatusCode() == 1 && mcp.PdgCode() == 22)) continue; + exiting_photon_parents.push_back(mcp.Mother()); + } + + std::vector in_nucleus_photons; + for(size_t const s : exiting_photon_parents) { + simb::MCParticle const & mcp = mct.GetParticle(s); + // Require parent vertex inside TPC. + if(std::abs(mcp.Vx()) > 210 || std::abs(mcp.Vy()) > 210 || mcp.Vz() > 510 || mcp.Vz() < -1) { + std::cout << "OUTSIDE TPC x y z =" << mcp.Vx() << " " << mcp.Vy() << " " << mcp.Vz() << std::endl; + return false; + } + if(isDelta(mcp.PdgCode())) { + if(ftree) FillTree(e, i, s, true); + return true; + } + else if(mcp.PdgCode() == 22) { + in_nucleus_photons.push_back(mcp.Mother()); + } + } + + for(size_t const s : in_nucleus_photons) { + simb::MCParticle const & mcp = mct.GetParticle(s); + if(isDelta(mcp.PdgCode())) { + if(ftree) FillTree(e, i, s, true); + return true; + } + } + + } + + if(ftree) FillTree(e, 0, SIZE_MAX, false); + return false; + +} + + +DEFINE_ART_MODULE(CCDeltaRadiative) diff --git a/sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl b/sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl new file mode 100644 index 000000000..71178e589 --- /dev/null +++ b/sbncode/SinglePhotonAnalysis/jobs/prodgenie_nu_singleinteraction_tpc_CCDeltaRadiative_filtered_sbnd.fcl @@ -0,0 +1,106 @@ +# Simulates GENIE neutrino interactions from the BNB beam +# forcing one interaction per event, inside the TPC volume +# (with 10 cm padding on each side), +# selecting only CC events with photons coming out from the Delta + + +#include "prodgenie_nu_singleinteraction_tpc_sbnd.fcl" + +#------ this could be a separated file +# +# services +# + +# +# modules +# + + + +process_name: GenieGenFiltered + + +# +# services +# +services: { + + TFileService: { fileName: "hists_genie.root" } + IFDH: {} # required by GENIEGen + @table::sbnd_basic_services # from simulationservices_sbnd.fcl + @table::sbnd_random_services # from simulationservices_sbnd.fcl + FileCatalogMetadata: @local::sbnd_file_catalog_mc # from sam_sbnd.fcl + + # since this is a configuration expected to be run for production, + # we set up message configuration accordingly: + message: @local::sbnd_message_services_prod + +} # services + + +# +# input +# + + +# +# processing +# +physics: +{ + + producers: + { + generator: @local::sbnd_genie_simple + rns: { module_type: "RandomNumberSaver" } + } + + filters: + { + CCDeltaRadFilter: + { + module_type: "CCDeltaRadiative" + } + + } + + + #define the producer and filter modules for this path, order matters, + simulate: [ rns, generator, CCDeltaRadFilter] + + #define the output stream, there could be more than one if using filters + stream1: [ out1 ] + + #trigger_paths is a keyword and contains the paths that modify the art::event, + #ie filters and producers + trigger_paths: [simulate] + + #end_paths is a keyword and contains the paths that do not modify the art::Event, + #ie analyzers and output streams. these all run simultaneously + end_paths: [stream1] + +} # physics + + +# +# output +# +outputs: +{ + out1: + { + module_type: RootOutput + fileName: "prodgenie_bnb_nu_filtered_CCDeltaRad_sbnd_%p-%tc.root" # default file name, can override from command line with -o or --output + SelectEvents: [simulate] + dataTier: "generated" + compressionLevel: 1 + } +} # outputs + + +# +# override +# THIS DOES NOT WORK, CHECK! +physics.producers.generator.TopVolume: "volTPCActive" +physics.producers.generator.BeamName: "booster" +physics.producers.generator.EventGeneratorList: "CCRES"