From 3f723b1f49ed55028af6700c4a4edec077784bce Mon Sep 17 00:00:00 2001 From: eyandel Date: Wed, 4 Mar 2026 10:17:36 -0600 Subject: [PATCH 01/19] move Lee's del 1g and isotropic 1g code into new NGEM directory and modify CMaKELists so it compiles --- CMakeLists.txt | 107 +++ ubsim/CMakeLists.txt | 1 + ubsim/NGEM/CMakeLists.txt | 54 ++ ubsim/NGEM/CMake_bad2.txt | 38 + ubsim/NGEM/CMake_noroot.txt | 54 ++ ubsim/NGEM/DeleteOneRandomPhoton.h | 72 ++ ubsim/NGEM/GENIEGenNGEM_module.cc | 746 ++++++++++++++++++ ubsim/NGEM/GenerateIsotropicSinglePhoton.h | 135 ++++ ubsim/NGEM/ManuallyDecayPi0sToTwoPhotons.h | 112 +++ ubsim/NGEM/mcc10_cmake.txt | 42 + ubsim/NGEM/prodgenie_bnb_nu_del_1g_uboone.fcl | 28 + ubsim/NGEM/prodgenie_isotropic_1g_uboone.fcl | 26 + 12 files changed, 1415 insertions(+) create mode 100644 ubsim/NGEM/CMakeLists.txt create mode 100644 ubsim/NGEM/CMake_bad2.txt create mode 100644 ubsim/NGEM/CMake_noroot.txt create mode 100644 ubsim/NGEM/DeleteOneRandomPhoton.h create mode 100644 ubsim/NGEM/GENIEGenNGEM_module.cc create mode 100644 ubsim/NGEM/GenerateIsotropicSinglePhoton.h create mode 100644 ubsim/NGEM/ManuallyDecayPi0sToTwoPhotons.h create mode 100644 ubsim/NGEM/mcc10_cmake.txt create mode 100644 ubsim/NGEM/prodgenie_bnb_nu_del_1g_uboone.fcl create mode 100644 ubsim/NGEM/prodgenie_isotropic_1g_uboone.fcl diff --git a/CMakeLists.txt b/CMakeLists.txt index 934e1d093..d866c7274 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -58,6 +58,7 @@ find_ups_product( art v1_00_00 ) find_ups_product( larcore ) find_ups_product( nutools ) find_ups_product( ppfx ) +find_ups_product( genie ) # WireCell required find_ups_product( larwirecell v1_00_00 ) @@ -82,6 +83,112 @@ include(ArtDictionary) include(ArtMake) include(BuildPlugins) +# Determine the major version number of GENIE based on the GENIE_VERSION +# variable set by the genie ups product. Define the GENIE_PRE_R3 preprocessor +# macro (needed to set the right GENIE header locations, etc.) if needed. +# +# The top-level nutools CMakeLists.txt does this as well, but I don't think +# that we can access the result of that version check in the larsim build. +parse_ups_version( ${GENIE_VERSION} ) +include_directories ( $ENV{GENIE_INC}/GENIE ) +if( ${VMAJ} LESS 3 ) + + add_definitions(-DGENIE_PRE_R3) + + # GENIE v2 library list + set(GENIE_LIB_LIST ${GBARYONRESONANCE} + ${GBASE} + ${GBODEKYANG} + ${GCHARM} + ${GCOH} + ${GDFRC} + ${GDIS} + ${GCROSSSECTIONS} + ${GDECAY} + ${GELAS} + ${GELFF} + ${GEVGCORE} + ${GEVGMODULES} + ${GGIBUU} + ${GHADRONTRANSP} + ${GFRAGMENTATION} + ${GLLEWELLYNSMITH} + ${GMEC} + ${GNUGAMMA} + ${GNUE} + ${GNUCLEAR} + ${GNUMERICAL} + ${GQPM} + ${GPDG} + ${GPDF} + ${GQEL} + ${GRES} + ${GREGISTRY} + ${GREINSEHGAL} + ${GGEO} + ${GMUELOSS} + ${GREWEIGHT} + ${GNUCLEONDECAY} ) + + #MESSAGE("--LARSIM-- GENIE_VERSION=${VMAJ} ${VMIN} -DGENIE_PRE_R3") +else() + + # GENIE v3 library list + set (GENIE_LIB_LIST ${GFWMSG} + ${GFWREG} + ${GFWALG} + ${GFWINT} + ${GFWGHEP} + ${GFWNUM} + ${GSL} # FWNUM relies on GSL + ${GFWUTL} + ${GFWPARDAT} + ${GFWEG} + ${GFWNTP} + ${GPHXSIG} + ${GPHPDF} + ${GPHNUCLST} + ${GPHCMN} + ${GPHDCY} + ${GPHHADTRANSP} + ${GPHHADNZ} + ${GPHDEEX} + ${GPHAMNGXS} + ${GPHAMNGEG} + ${GPHCHMXS} + ${GPHCOHXS} + ${GPHCOHEG} + ${GPHDISXS} + ${GPHDISEG} + ${GPHDFRCXS} + ${GPHDFRCEG} + ${GPHGLWRESXS} + ${GPHGLWRESEG} + ${GPHIBDXS} + ${GPHIBDEG} + ${GPHMNUCXS} + ${GPHMNUCEG} + ${GPHMEL} + ${GPHNUELXS} + ${GPHNUELEG} + ${GPHQELXS} + ${GPHQELEG} + ${GPHRESXS} + ${GPHRESEG} + ${GPHSTRXS} + ${GPHSTREG} + ${GPHNDCY} + ${GTLGEO} + ${GTLFLX} + ${GRWFWK} + ${GRWIO} + ${GRWCLC}) + + #MESSAGE("--LARSIM-- GENIE_VERSION=${VMAJ} ${VMIN} not defining -DGENIE_PRE_R3") +endif() + +#MESSAGE("--LARSIM-- GENIE ${VMAJ} GENIE_LIB_LIST=${GENIE_LIB_LIST}") + # ADD SOURCE CODE SUBDIRECTORIES HERE add_subdirectory(ubsim) diff --git a/ubsim/CMakeLists.txt b/ubsim/CMakeLists.txt index feb833eb8..24df4cdc6 100644 --- a/ubsim/CMakeLists.txt +++ b/ubsim/CMakeLists.txt @@ -6,3 +6,4 @@ add_subdirectory(PhotonPropagation) add_subdirectory(SNStreamSim) add_subdirectory(Simulation) add_subdirectory(Filters) +add_subdirectory(NGEM) \ No newline at end of file diff --git a/ubsim/NGEM/CMakeLists.txt b/ubsim/NGEM/CMakeLists.txt new file mode 100644 index 000000000..fa5527f53 --- /dev/null +++ b/ubsim/NGEM/CMakeLists.txt @@ -0,0 +1,54 @@ + +cet_add_compiler_flags( CXX -DHIDE_GENIE_LOG_XXX ) + +include_directories ( $ENV{LOG4CPP_INC} ) +include_directories ( $ENV{DK2NUDATA_INC} ) +include_directories ( $ENV{DK2NUGENIE_INC} ) + +# genie +cet_find_library( DK2NUDATA NAMES dk2nuTree PATHS ENV DK2NUDATA_LIB ) +cet_find_library( DK2NUGENIE NAMES dk2nuGenie PATHS ENV DK2NUGENIE_LIB ) + +art_make( + MODULE_LIBRARIES + lardataalg_MCDumpers + larsim_Simulation lardataobj_Simulation + nusimdata_SimulationBase + larcoreobj_SummaryData + larcorealg_Geometry + larcore_Geometry_Geometry_service + nutools_RandomUtils_NuRandomService_service + nutools_EventGeneratorBase_GENIE + #nurandom_RandomUtils_NuRandomService_service + #nugen_EventGeneratorBase_GENIE + ${DK2NUDATA} + ${DK2NUGENIE} + #dk2nu_Tree + #dk2nu_Genie + ${ART_FRAMEWORK_CORE} + ${ART_FRAMEWORK_PRINCIPAL} + ${ART_FRAMEWORK_SERVICES_REGISTRY} + ${ART_FRAMEWORK_SERVICES_OPTIONAL} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_TFILESERVICE_SERVICE} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_RANDOMNUMBERGENERATOR_SERVICE} + art_Persistency_Common + art_Persistency_Provenance + art_Utilities + #art_root_io_TFileService_service + canvas + ${MF_MESSAGELOGGER} + ${FHICLCPP} + cetlib cetlib_except + ${CLHEP} + ${GENIE_LIB_LIST} + ${ROOT_EGPYTHIA6} # FIXME!!! - resolving genie run time reference + ${ROOT_BASIC_LIB_LIST} + ${ROOT_EG} + ${LOG4CPP} + #Generator_Framework_GHEP + ) + +install_headers() +install_fhicl() +install_source() + diff --git a/ubsim/NGEM/CMake_bad2.txt b/ubsim/NGEM/CMake_bad2.txt new file mode 100644 index 000000000..9b811a96f --- /dev/null +++ b/ubsim/NGEM/CMake_bad2.txt @@ -0,0 +1,38 @@ +cet_build_plugin( + GENIEGenNGEM art::EDProducer + LIBRARIES + PRIVATE + ubsim::GENIEGenNGEM + fhiclcpp::fhiclcpp + lardata::Utilities + art_root_io::TFileService_service + nurandom::RandomUtils_NuRandomService_service +) + +cet_make_library( + SOURCE + LIBRARIES + PUBLIC + larsim::EventWeight_Base + larcore::Geometry_Geometry_service + nusimdata::SimulationBase + nugen::EventGeneratorBase_GENIE + Geant4::G4event + ppfx::ppfx + dk2nu::Tree + art_root_io::TFileService_service + fhiclcpp::fhiclcpp + geant4reweight::ReweightBaseLib + GENIE::GRwClc + GENIE::GFwInt + GENIE::GRwFwk + GENIE::GFwMsg + GENIE::GFwUtl + log4cpp::log4cpp + ROOT::Hist + ROOT::Tree +) + +install_headers() +install_fhicl() +install_source() \ No newline at end of file diff --git a/ubsim/NGEM/CMake_noroot.txt b/ubsim/NGEM/CMake_noroot.txt new file mode 100644 index 000000000..1226c29e0 --- /dev/null +++ b/ubsim/NGEM/CMake_noroot.txt @@ -0,0 +1,54 @@ + +cet_add_compiler_flags( CXX -DHIDE_GENIE_LOG_XXX ) + +include_directories ( $ENV{LOG4CPP_INC} ) +include_directories ( $ENV{DK2NUDATA_INC} ) +include_directories ( $ENV{DK2NUGENIE_INC} ) + +# genie +#cet_find_library( DK2NUDATA NAMES dk2nuTree PATHS ENV DK2NUDATA_LIB ) +#cet_find_library( DK2NUGENIE NAMES dk2nuGenie PATHS ENV DK2NUGENIE_LIB ) +find_library(DK2NUDATA dk2nuTree PATHS ENV DK2NUDATA_LIB REQUIRED) +find_library(DK2NUGENIE dk2nuGenie PATHS ENV DK2NUGENIE_LIB REQUIRED) + + +art_make( + MODULE_LIBRARIES + lardataalg_MCDumpers + larsim_Simulation lardataobj_Simulation + nusimdata_SimulationBase + larcoreobj_SummaryData + larcorealg_Geometry + larcore_Geometry_Geometry_service + nutools_RandomUtils_NuRandomService_service + nutools_EventGeneratorBase_GENIE + ${DK2NUDATA} + ${DK2NUGENIE} + ${ART_FRAMEWORK_CORE} + ${ART_FRAMEWORK_PRINCIPAL} + ${ART_FRAMEWORK_SERVICES_REGISTRY} + ${ART_FRAMEWORK_SERVICES_OPTIONAL} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_TFILESERVICE_SERVICE} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_RANDOMNUMBERGENERATOR_SERVICE} + art_Persistency_Common + art_Persistency_Provenance + art_Utilities + canvas + ${MF_MESSAGELOGGER} + ${FHICLCPP} + cetlib cetlib_except + ${CLHEP} + ${GENIE_LIB_LIST} + ${ROOT_EGPYTHIA6} # FIXME!!! - resolving genie run time reference + ${ROOT_BASIC_LIB_LIST} + ${ROOT_EG} + ${ROOT_GEOM} + ${ROOT_XMLIO} + ${ROOT_GDML} + ${LOG4CPP} + ) + +install_headers() +install_fhicl() +install_source() + diff --git a/ubsim/NGEM/DeleteOneRandomPhoton.h b/ubsim/NGEM/DeleteOneRandomPhoton.h new file mode 100644 index 000000000..12112bb68 --- /dev/null +++ b/ubsim/NGEM/DeleteOneRandomPhoton.h @@ -0,0 +1,72 @@ +#ifndef EVGEN_DELETE_ONE_RANDOM_PHOTON_H +#define EVGEN_DELETE_ONE_RANDOM_PHOTON_H + +#include +#include + +// ROOT +#include "TRandom3.h" + +// LArSoft +#include "nusimdata/SimulationBase/MCTruth.h" +#include "nusimdata/SimulationBase/MCParticle.h" + +namespace evgen { + + inline void DeleteOneRandomPhoton(simb::MCTruth& originalMCTruth, simb::MCTruth& newMCTruth) { + TRandom3 randomGen; + randomGen.SetSeed(0); + std::vector photon_indices; + for (int i = 0; i < originalMCTruth.NParticles(); ++i) { + int pdgCode = originalMCTruth.GetParticle(i).PdgCode(); + int statusCode = originalMCTruth.GetParticle(i).StatusCode(); + if (pdgCode == 22 && statusCode == 1) { + photon_indices.push_back(i); + } + } + int num_photons = photon_indices.size(); + std::cout << "Number of primary photons: " << num_photons << std::endl; + if (num_photons == 0) { + std::cout << "No photons to delete" << std::endl; + for (int i = 0; i < originalMCTruth.NParticles(); ++i) { + const simb::MCParticle& particle = originalMCTruth.GetParticle(i); + simb::MCParticle non_const_particle = simb::MCParticle(particle); + newMCTruth.Add(non_const_particle); + } + } else { + int index_to_delete = photon_indices[randomGen.Integer(num_photons)]; + for (int i = 0; i < originalMCTruth.NParticles(); ++i) { + if (i == index_to_delete) { + std::cout << "Deleted photon at index: " << index_to_delete << std::endl; + } else { // copying over particles that aren't getting deleted + const simb::MCParticle& particle = originalMCTruth.GetParticle(i); + simb::MCParticle non_const_particle = simb::MCParticle(particle); + newMCTruth.Add(non_const_particle); + } + } + } + simb::MCNeutrino neutrino = originalMCTruth.GetNeutrino(); + int CCNC = neutrino.CCNC(); + int mode = neutrino.Mode(); + int interactionType = neutrino.InteractionType(); + int target = neutrino.Target(); + int nucleon = neutrino.HitNuc(); + int quark = neutrino.HitQuark(); + double w = neutrino.W(); + double x = neutrino.X(); + double y = neutrino.Y(); + double qsqr = neutrino.QSqr(); + newMCTruth.SetNeutrino(CCNC, mode, interactionType, target, nucleon, quark, w, x, y, qsqr); + newMCTruth.SetOrigin(originalMCTruth.Origin()); + } + + // In-place convenience wrapper + inline void DeleteOneRandomPhoton(simb::MCTruth& mcTruth) { + simb::MCTruth newTruth; + DeleteOneRandomPhoton(mcTruth, newTruth); + mcTruth = newTruth; + } + +} // namespace evgen + +#endif // EVGEN_DELETE_ONE_RANDOM_PHOTON_H diff --git a/ubsim/NGEM/GENIEGenNGEM_module.cc b/ubsim/NGEM/GENIEGenNGEM_module.cc new file mode 100644 index 000000000..c7869888e --- /dev/null +++ b/ubsim/NGEM/GENIEGenNGEM_module.cc @@ -0,0 +1,746 @@ +//////////////////////////////////////////////////////////////////////// +// +// +// GENIE neutrino event generator +// +// brebel@fnal.gov +// +//////////////////////////////////////////////////////////////////////// +#ifndef EVGEN_GENIEGEN_H +#define EVGEN_GENIEGEN_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// ROOT includes +#include "TH1.h" +#include "TH2.h" +#include "TDatabasePDG.h" +#include "TSystem.h" +#include "TStopwatch.h" + +// Framework includes +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/SubRun.h" +#include "fhiclcpp/ParameterSet.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art/Framework/Services/Optional/TFileService.h" +#include "art/Framework/Services/Optional/TFileDirectory.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "canvas/Persistency/Common/Assns.h" +#include "art/Framework/Core/EDProducer.h" + +// art extensions +#include "nutools/RandomUtils/NuRandomService.h" + +// LArSoft includes +#include "lardataalg/MCDumpers/MCDumpers.h" // sim::dump namespace +#include "nusimdata/SimulationBase/MCTruth.h" +#include "nusimdata/SimulationBase/MCFlux.h" +#include "nusimdata/SimulationBase/GTruth.h" +#include "lardataobj/Simulation/BeamGateInfo.h" +#include "larcore/Geometry/Geometry.h" +#include "larcoreobj/SummaryData/RunData.h" +#include "larcoreobj/SummaryData/POTSummary.h" +#include "nutools/EventGeneratorBase/GENIE/GENIEHelper.h" +#include "lardata/Utilities/AssociationUtil.h" + +// dk2nu extensions +#include "dk2nu/tree/dk2nu.h" +#include "dk2nu/tree/NuChoice.h" +#include "dk2nu/genie/GDk2NuFlux.h" + +#include "GENIE/Framework/EventGen/EventRecord.h" +#include "nutools/EventGeneratorBase/GENIE/EVGBAssociationUtil.h" +#include "nutools/EventGeneratorBase/evgenbase.h" + +// Optional Modifications, if you don't want to modify normal GENIE neutrino generation events +#include "ManuallyDecayPi0sToTwoPhotons.h" +#include "DeleteOneRandomPhoton.h" +#include "GenerateIsotropicSinglePhoton.h" + +///Event Generation using GENIE, cosmics or single particles +namespace evgen { + /** + * @brief A module to check the results from the Monte Carlo generator + * + * Note on random number generator + * -------------------------------- + * + * GENIE uses a TRandom generator for its purposes. + * Since art's RandomNumberGenerator service only provides + * `CLHEP::HepRandomEngine`, the standard LArSoft/art mechanism for handling + * the random stream can't be used. + * GENIEHelper, interface to GENIE provided by nutools, creates a TRandom + * that GENIE can use. It initializes it with a random seed read from + * *RandomSeed* configuration parameter. This and all the other parameters + * are inherited from the art module (that is, `GENIEGen`) configuration. + * LArSoft meddles with this mechanism to provide support for the standard + * "Seed" parameter and NuRandomService service. + * + * Configuration parameters + * ------------------------- + * + * - *RandomSeed* (integer, optional): if specified, this value is used as + * seed for GENIE random number generator engine + * - *Seed* (unsigned integer, optional): if specified, this value is used as + * seed for GENIE random number generator engine; if *RandomSeed* is also + * specified, this value is ignored (but in the future this could turn into + * a configuration error) + * + * As custom, if the random seed is not provided by the configuration, one is + * fetched from `NuRandomService` (if available), with the behaviour in + * lar::util::FetchRandomSeed(). + */ + class GENIEGenNGEM : public art::EDProducer { + public: + explicit GENIEGenNGEM(fhicl::ParameterSet const &pset); + virtual ~GENIEGenNGEM(); + + void produce(art::Event& evt); + void beginJob(); + void beginRun(art::Run& run); + void beginSubRun(art::SubRun& sr); + void endSubRun(art::SubRun& sr); + + private: + + std::string ParticleStatus(int StatusCode); + std::string ReactionChannel(int ccnc,int mode); + + void FillHistograms(simb::MCTruth mc); + + evgb::GENIEHelper *fGENIEHelp; ///< GENIEHelper object + bool fDefinedVtxHistRange;///use defined hist range; it is useful to have for asymmetric ranges like in DP FD. + std::vector< double > fVtxPosHistRange; + + int fPassEmptySpills; ///< whether or not to kill evnets with no interactions + TStopwatch fStopwatch; ///keep track of how long it takes to run the job + + double fGlobalTimeOffset; /// The start of a simulated "beam gate". + double fRandomTimeOffset; /// The width of a simulated "beam gate". + ::sim::BeamType_t fBeamType; /// The type of beam + + double fPrevTotPOT; ///< Total POT from subruns previous to current subrun + double fPrevTotGoodPOT; ///< Total good POT from subruns previous to current subrun + + TH1F* fGenerated[6]; ///< Spectra as generated + + TH1F* fVertexX; ///< vertex location of generated events in x + TH1F* fVertexY; ///< vertex location of generated events in y + TH1F* fVertexZ; ///< vertex location of generated events in z + + TH2F* fVertexXY; ///< vertex location in xy + TH2F* fVertexXZ; ///< vertex location in xz + TH2F* fVertexYZ; ///< vertex location in yz + + TH1F* fDCosX; ///< direction cosine in x + TH1F* fDCosY; ///< direction cosine in y + TH1F* fDCosZ; ///< direction cosine in z + + TH1F* fMuMomentum; ///< momentum of outgoing muons + TH1F* fMuDCosX; ///< direction cosine of outgoing mu in x + TH1F* fMuDCosY; ///< direction cosine of outgoing mu in y + TH1F* fMuDCosZ; ///< direction cosine of outgoing mu in z + + TH1F* fEMomentum; ///< momentum of outgoing electrons + TH1F* fEDCosX; ///< direction cosine of outgoing e in x + TH1F* fEDCosY; ///< direction cosine of outgoing e in y + TH1F* fEDCosZ; ///< direction cosine of outgoing e in z + + TH1F* fCCMode; ///< CC interaction mode + TH1F* fNCMode; ///< CC interaction mode + + TH1F* fDeltaE; ///< difference in neutrino energy from MCTruth::Enu() vs TParticle + TH1F* fECons; ///< histogram to determine if energy is conserved in the event + + bool ManuallyDecayPi0s; ///< whether to manually decay pi0s to photons + bool DeleteRandomGamma; ///< whether to delete a random gamma + bool GenerateIsotropicSinglePhoton; ///< replace all particles with single isotropic gamma + std::vector SinglePhotonEnergyBinEdges; ///< configurable via FHiCL + std::vector SinglePhotonEnergyBinProbs; ///< configurable via FHiCL + std::vector SinglePhotonCosThetaBinEdges; ///< configurable via FHiCL + std::vector SinglePhotonCosThetaBinProbs; ///< configurable via FHiCL + + }; +} + +namespace evgen{ + + //____________________________________________________________________________ + GENIEGenNGEM::GENIEGenNGEM(fhicl::ParameterSet const& pset) + : fGENIEHelp(0) + , fDefinedVtxHistRange (pset.get< bool >("DefinedVtxHistRange")) + , fVtxPosHistRange (pset.get< std::vector >("VtxPosHistRange")) + , fPassEmptySpills (pset.get< bool >("PassEmptySpills")) + , fGlobalTimeOffset(pset.get< double >("GlobalTimeOffset",0)) + , fRandomTimeOffset(pset.get< double >("RandomTimeOffset",1600.)) // BNB default value + , fBeamType(::sim::kBNB) + , ManuallyDecayPi0s(pset.get("ManuallyDecayPi0s", false)) + , DeleteRandomGamma(pset.get("DeleteRandomGamma", false)) + , GenerateIsotropicSinglePhoton(pset.get("GenerateIsotropicSinglePhoton", false)) + , SinglePhotonEnergyBinEdges(pset.get>("SinglePhotonEnergyBinEdges", std::vector())) + , SinglePhotonEnergyBinProbs(pset.get>("SinglePhotonEnergyBinProbs", std::vector())) + , SinglePhotonCosThetaBinEdges(pset.get>("SinglePhotonCosThetaBinEdges", std::vector())) + , SinglePhotonCosThetaBinProbs(pset.get>("SinglePhotonCosThetaBinProbs", std::vector())) + { + fStopwatch.Start(); + + produces< std::vector >(); + produces< std::vector >(); + produces< std::vector >(); + produces< sumdata::RunData, art::InRun >(); + produces< sumdata::POTSummary, art::InSubRun >(); + produces< art::Assns >(); + produces< art::Assns >(); + produces< std::vector >(); + + // dk2nu additions + produces< std::vector >(); + produces< std::vector >(); + produces< art::Assns >(); + produces< art::Assns >(); + + std::string beam_type_name = pset.get("BeamName"); + + if(beam_type_name == "numi") + + fBeamType = ::sim::kNuMI; + + else if(beam_type_name == "booster") + + fBeamType = ::sim::kBNB; + + else + + fBeamType = ::sim::kUnknown; + + art::ServiceHandle geo; + + signed int temp_seed; // the seed read by GENIEHelper is a signed integer... + fhicl::ParameterSet GENIEconfig(pset); + if (!GENIEconfig.get_if_present("RandomSeed", temp_seed)) { // TODO use has_key() when it becomes available + // no RandomSeed specified; check for the LArSoft-style "Seed" instead: + // obtain the random seed from a service, + // unless overridden in configuration with key "Seed" + unsigned int seed; + if (!GENIEconfig.get_if_present("Seed", seed)) + seed = art::ServiceHandle()->getSeed(); + + // The seed is not passed to RandomNumberGenerator, + // since GENIE uses a TRandom generator that is owned by the GENIEHelper. + // Instead, we explicitly configure the random seed for GENIEHelper: + GENIEconfig.put("RandomSeed", seed); + } // if no RandomSeed present + + fGENIEHelp = new evgb::GENIEHelper(GENIEconfig, + geo->ROOTGeoManager(), + geo->ROOTFile(), + geo->TotalMass(pset.get< std::string>("TopVolume").c_str())); + + } + + //____________________________________________________________________________ + GENIEGenNGEM::~GENIEGenNGEM() + { + if(fGENIEHelp) delete fGENIEHelp; + fStopwatch.Stop(); + mf::LogInfo("GENIEProductionTime") << "real time to produce file: " << fStopwatch.RealTime(); + } + + //____________________________________________________________________________ + void GENIEGenNGEM::beginJob(){ + fGENIEHelp->Initialize(); + + fPrevTotPOT = 0.; + fPrevTotGoodPOT = 0.; + + // Get access to the TFile service. + art::ServiceHandle tfs; + + fGenerated[0] = tfs->make("fGenerated_necc","", 100, 0.0, 20.0); + fGenerated[1] = tfs->make("fGenerated_nebcc","", 100, 0.0, 20.0); + fGenerated[2] = tfs->make("fGenerated_nmcc","", 100, 0.0, 20.0); + fGenerated[3] = tfs->make("fGenerated_nmbcc","", 100, 0.0, 20.0); + fGenerated[4] = tfs->make("fGenerated_nnc","", 100, 0.0, 20.0); + fGenerated[5] = tfs->make("fGenerated_nbnc","", 100, 0.0, 20.0); + + fDCosX = tfs->make("fDCosX", ";dx/ds", 200, -1., 1.); + fDCosY = tfs->make("fDCosY", ";dy/ds", 200, -1., 1.); + fDCosZ = tfs->make("fDCosZ", ";dz/ds", 200, -1., 1.); + + fMuMomentum = tfs->make("fMuMomentum", ";p_{#mu} (GeV/c)", 500, 0., 50.); + fMuDCosX = tfs->make("fMuDCosX", ";dx/ds;", 200, -1., 1.); + fMuDCosY = tfs->make("fMuDCosY", ";dy/ds;", 200, -1., 1.); + fMuDCosZ = tfs->make("fMuDCosZ", ";dz/ds;", 200, -1., 1.); + + fEMomentum = tfs->make("fEMomentum", ";p_{e} (GeV/c)", 500, 0., 50.); + fEDCosX = tfs->make("fEDCosX", ";dx/ds;", 200, -1., 1.); + fEDCosY = tfs->make("fEDCosY", ";dy/ds;", 200, -1., 1.); + fEDCosZ = tfs->make("fEDCosZ", ";dz/ds;", 200, -1., 1.); + + fCCMode = tfs->make("fCCMode", ";CC Interaction Mode;", 4, 0., 4.); + fCCMode->GetXaxis()->SetBinLabel(1, "QE"); + fCCMode->GetXaxis()->SetBinLabel(2, "Res"); + fCCMode->GetXaxis()->SetBinLabel(3, "DIS"); + fCCMode->GetXaxis()->SetBinLabel(4, "Coh"); + fCCMode->GetXaxis()->CenterLabels(); + + fNCMode = tfs->make("fNCMode", ";NC Interaction Mode;", 4, 0., 4.); + fNCMode->GetXaxis()->SetBinLabel(1, "QE"); + fNCMode->GetXaxis()->SetBinLabel(2, "Res"); + fNCMode->GetXaxis()->SetBinLabel(3, "DIS"); + fNCMode->GetXaxis()->SetBinLabel(4, "Coh"); + fNCMode->GetXaxis()->CenterLabels(); + + fDeltaE = tfs->make("fDeltaE", ";#Delta E_{#nu} (GeV);", 200, -1., 1.); + fECons = tfs->make("fECons", ";#Delta E(#nu,lepton);", 500, -5., 5.); + + art::ServiceHandle geo; + double x = 2.1*geo->DetHalfWidth(); + double y = 2.1*geo->DetHalfHeight(); + double z = 2.*geo->DetLength(); + int xdiv = TMath::Nint(2*x/5.); + int ydiv = TMath::Nint(2*y/5.); + int zdiv = TMath::Nint(2*z/5.); + + if (fDefinedVtxHistRange == false) + { + fVertexX = tfs->make("fVertexX", ";x (cm)", xdiv, -0.1*x, x); + fVertexY = tfs->make("fVertexY", ";y (cm)", ydiv, -y, y); + fVertexZ = tfs->make("fVertexZ", ";z (cm)", zdiv, -0.1*z, z); + + fVertexXY = tfs->make("fVertexXY", ";x (cm);y (cm)", xdiv, -0.1*x, x, ydiv, -y, y); + fVertexXZ = tfs->make("fVertexXZ", ";z (cm);x (cm)", zdiv, -0.2*z, z, xdiv, -0.1*x, x); + fVertexYZ = tfs->make("fVertexYZ", ";z (cm);y (cm)", zdiv, -0.2*z, z, ydiv, -y, y); + } + else + { + fVertexX = tfs->make("fVertexX", ";x (cm)", xdiv, fVtxPosHistRange[0], fVtxPosHistRange[1]); + fVertexY = tfs->make("fVertexY", ";y (cm)", ydiv, fVtxPosHistRange[2], fVtxPosHistRange[3]); + fVertexZ = tfs->make("fVertexZ", ";z (cm)", zdiv, fVtxPosHistRange[4], fVtxPosHistRange[5]); + + fVertexXY = tfs->make("fVertexXY", ";x (cm);y (cm)", xdiv, fVtxPosHistRange[0], fVtxPosHistRange[1], ydiv, fVtxPosHistRange[2], fVtxPosHistRange[3]); + fVertexXZ = tfs->make("fVertexXZ", ";z (cm);x (cm)", zdiv, fVtxPosHistRange[4], fVtxPosHistRange[5], xdiv, fVtxPosHistRange[0], fVtxPosHistRange[1]); + fVertexYZ = tfs->make("fVertexYZ", ";z (cm);y (cm)", zdiv, fVtxPosHistRange[4], fVtxPosHistRange[5], ydiv, fVtxPosHistRange[2], fVtxPosHistRange[3]); + } + + } + + //____________________________________________________________________________ + void GENIEGenNGEM::beginRun(art::Run& run) + { + + // grab the geometry object to see what geometry we are using + art::ServiceHandle geo; + std::unique_ptr runcol(new sumdata::RunData(geo->DetectorName())); + + run.put(std::move(runcol)); + + return; + } + + //____________________________________________________________________________ + void GENIEGenNGEM::beginSubRun(art::SubRun& sr) + { + + fPrevTotPOT = fGENIEHelp->TotalExposure(); + fPrevTotGoodPOT = fGENIEHelp->TotalExposure(); + + return; + } + + //____________________________________________________________________________ + void GENIEGenNGEM::endSubRun(art::SubRun& sr) + { + + std::unique_ptr p(new sumdata::POTSummary()); + + p->totpot = fGENIEHelp->TotalExposure() - fPrevTotPOT; + p->totgoodpot = fGENIEHelp->TotalExposure() - fPrevTotGoodPOT; + + sr.put(std::move(p)); + + return; + } + + //____________________________________________________________________________ + void GENIEGenNGEM::produce(art::Event& evt) + { + std::unique_ptr< std::vector > truthcol (new std::vector); + std::unique_ptr< std::vector > fluxcol (new std::vector); + std::unique_ptr< std::vector > gtruthcol (new std::vector); + std::unique_ptr< art::Assns > tfassn(new art::Assns); + std::unique_ptr< art::Assns > tgtassn(new art::Assns); + std::unique_ptr< std::vector > gateCollection(new std::vector); + + std::unique_ptr< std::vector > + dk2nucol(new std::vector); + std::unique_ptr< std::vector > + nuchoicecol(new std::vector); + std::unique_ptr< art::Assns > + dk2nuassn(new art::Assns); + std::unique_ptr< art::Assns > + nuchoiceassn(new art::Assns); + + while(truthcol->size() < 1){ + while(!fGENIEHelp->Stop()){ + + simb::MCTruth truth; + simb::MCFlux flux; + simb::GTruth gTruth; + + // GENIEHelper returns a false in the sample method if + // either no neutrino was generated, or the interaction + // occurred beyond the detector's z extent - ie something we + // would never see anyway. + if(fGENIEHelp->Sample(truth, flux, gTruth)){ + + std::cout << "truth.NParticles(): " << truth.NParticles() << std::endl; + for (int i = 0; i < truth.NParticles(); ++i) { + std::cout << " pdg: " << truth.GetParticle(i).PdgCode(); + std::cout << ", track_id: " << truth.GetParticle(i).TrackId(); + std::cout << ", mother: " << truth.GetParticle(i).Mother(); + std::cout << ", mass: " << truth.GetParticle(i).Mass(); + std::cout << ", position: (" << truth.GetParticle(i).Vx() << ", " << truth.GetParticle(i).Vy() << ", " << truth.GetParticle(i).Vz() << ")"; + std::cout << ", momentum: (" << truth.GetParticle(i).Px() << ", " << truth.GetParticle(i).Py() << ", " << truth.GetParticle(i).Pz() << ")"; + std::cout << ", Gvtx: (" << truth.GetParticle(i).Gvx() << ", " << truth.GetParticle(i).Gvy() << ", " << truth.GetParticle(i).Gvz() << ", " << truth.GetParticle(i).Gvt() << ")"; + std::cout << ", status code: " << truth.GetParticle(i).StatusCode(); + std::cout << ", process: " << truth.GetParticle(i).Process(); + std::cout << std::endl; + } + + if (GenerateIsotropicSinglePhoton) { + std::cout << "Generating isotropic single photon" << std::endl; + GenerateIsotropicSinglePhotonWithBins( + SinglePhotonEnergyBinEdges, SinglePhotonEnergyBinProbs, + SinglePhotonCosThetaBinEdges, SinglePhotonCosThetaBinProbs, + truth + ); + std::cout << "truth.NParticles(): " << truth.NParticles() << std::endl; + for (int i = 0; i < truth.NParticles(); ++i) { + std::cout << " pdg: " << truth.GetParticle(i).PdgCode(); + std::cout << ", track_id: " << truth.GetParticle(i).TrackId(); + std::cout << ", mother: " << truth.GetParticle(i).Mother(); + std::cout << ", mass: " << truth.GetParticle(i).Mass(); + std::cout << ", position: (" << truth.GetParticle(i).Vx() << ", " << truth.GetParticle(i).Vy() << ", " << truth.GetParticle(i).Vz() << ")"; + std::cout << ", momentum: (" << truth.GetParticle(i).Px() << ", " << truth.GetParticle(i).Py() << ", " << truth.GetParticle(i).Pz() << ")"; + std::cout << ", Gvtx: (" << truth.GetParticle(i).Gvx() << ", " << truth.GetParticle(i).Gvy() << ", " << truth.GetParticle(i).Gvz() << ", " << truth.GetParticle(i).Gvt() << ")"; + std::cout << ", status code: " << truth.GetParticle(i).StatusCode(); + std::cout << ", process: " << truth.GetParticle(i).Process(); + std::cout << std::endl; + } + } + + if (ManuallyDecayPi0s) { + std::cout << "Manually decaying pi0s" << std::endl; + ManuallyDecayPi0sToTwoPhotons(truth); + std::cout << "truth.NParticles(): " << truth.NParticles() << std::endl; + for (int i = 0; i < truth.NParticles(); ++i) { + std::cout << " pdg: " << truth.GetParticle(i).PdgCode(); + std::cout << ", track_id: " << truth.GetParticle(i).TrackId(); + std::cout << ", mother: " << truth.GetParticle(i).Mother(); + std::cout << ", mass: " << truth.GetParticle(i).Mass(); + std::cout << ", position: (" << truth.GetParticle(i).Vx() << ", " << truth.GetParticle(i).Vy() << ", " << truth.GetParticle(i).Vz() << ")"; + std::cout << ", momentum: (" << truth.GetParticle(i).Px() << ", " << truth.GetParticle(i).Py() << ", " << truth.GetParticle(i).Pz() << ")"; + std::cout << ", Gvtx: (" << truth.GetParticle(i).Gvx() << ", " << truth.GetParticle(i).Gvy() << ", " << truth.GetParticle(i).Gvz() << ", " << truth.GetParticle(i).Gvt() << ")"; + std::cout << ", status code: " << truth.GetParticle(i).StatusCode(); + std::cout << ", process: " << truth.GetParticle(i).Process(); + std::cout << std::endl; + } + } + + if (DeleteRandomGamma) { + std::cout << "Deleting random gamma" << std::endl; + DeleteOneRandomPhoton(truth); + std::cout << "truth.NParticles(): " << truth.NParticles() << std::endl; + for (int i = 0; i < truth.NParticles(); ++i) { + std::cout << " pdg: " << truth.GetParticle(i).PdgCode(); + std::cout << ", track_id: " << truth.GetParticle(i).TrackId(); + std::cout << ", mother: " << truth.GetParticle(i).Mother(); + std::cout << ", mass: " << truth.GetParticle(i).Mass(); + std::cout << ", position: (" << truth.GetParticle(i).Vx() << ", " << truth.GetParticle(i).Vy() << ", " << truth.GetParticle(i).Vz() << ")"; + std::cout << ", momentum: (" << truth.GetParticle(i).Px() << ", " << truth.GetParticle(i).Py() << ", " << truth.GetParticle(i).Pz() << ")"; + std::cout << ", Gvtx: (" << truth.GetParticle(i).Gvx() << ", " << truth.GetParticle(i).Gvy() << ", " << truth.GetParticle(i).Gvz() << ", " << truth.GetParticle(i).Gvt() << ")"; + std::cout << ", status code: " << truth.GetParticle(i).StatusCode(); + std::cout << ", process: " << truth.GetParticle(i).Process(); + std::cout << std::endl; + } + } + + truthcol ->push_back(truth); + fluxcol ->push_back(flux); + gtruthcol->push_back(gTruth); + util::CreateAssn(*this, evt, *truthcol, *fluxcol, *tfassn, fluxcol->size()-1, fluxcol->size()); + util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *tgtassn, gtruthcol->size()-1, gtruthcol->size()); + + FillHistograms(truth); + + genie::GFluxI* fdriver = fGENIEHelp->GetFluxDriver(true); + genie::flux::GDk2NuFlux* dk2nuDriver = + dynamic_cast(fdriver); + if ( dk2nuDriver ) { + const bsim::Dk2Nu& dk2nuObj = dk2nuDriver->GetDk2Nu(); + dk2nucol ->push_back(dk2nuObj); + const bsim::NuChoice& nuchoiceObj = dk2nuDriver->GetNuChoice(); + nuchoicecol->push_back(nuchoiceObj); + util::CreateAssn(*this, evt, *truthcol, *dk2nucol, *dk2nuassn, + dk2nucol->size()-1, dk2nucol->size()); + util::CreateAssn(*this, evt, *truthcol, *nuchoicecol, *nuchoiceassn, + nuchoicecol->size()-1, nuchoicecol->size()); + } + + // check that the process code is not unsupported by GENIE + // (see issue #18025 for reference); + // if it is, print all the information we can about this truth record + if (truth.NeutrinoSet() && (truth.GetNeutrino().InteractionType() == simb::kNuanceOffset)) { + mf::LogWarning log("GENIEmissingProcessMapping"); + log << "Found an interaction that is not represented by the interaction type code in GENIE:" + "\nMCTruth record:" + "\n" + ; + sim::dump::DumpMCTruth(log, truth, 2U); // 2 trajectory points per line + log << + "\nGENIE truth record:" + "\n" + ; + sim::dump::DumpGTruth(log, gTruth); + } // if + + }// end if genie was able to make an event + + }// end event generation loop + + // check to see if we are to pass empty spills + if(truthcol->size() < 1 && fPassEmptySpills){ + MF_LOG_DEBUG("GENIEGenNGEM") << "no events made for this spill but " + << "passing it on and ending the event anyway"; + break; + } + + }// end loop while no interactions are made + + // Create a simulated "beam gate" for these neutrino events. + // We're creating a vector of these because, in a + // distant-but-possible future, we may be generating more than one + // beam gate within a simulated time window. + gateCollection->push_back(sim::BeamGateInfo( fGlobalTimeOffset, fRandomTimeOffset, fBeamType )); + + // put the collections in the event + evt.put(std::move(truthcol)); + evt.put(std::move(fluxcol)); + evt.put(std::move(gtruthcol)); + evt.put(std::move(tfassn)); + evt.put(std::move(tgtassn)); + evt.put(std::move(gateCollection)); + + evt.put(std::move(dk2nucol)); + evt.put(std::move(nuchoicecol)); + evt.put(std::move(dk2nuassn)); + evt.put(std::move(nuchoiceassn)); + + return; + } + + //...................................................................... + std::string GENIEGenNGEM::ParticleStatus(int StatusCode) + { + int code = StatusCode; + std::string ParticleStatusName; + + switch(code) + { + case -1: + ParticleStatusName = "kIStUndefined"; + break; + case 0: + ParticleStatusName = "kIStInitialState"; + break; + case 1: + ParticleStatusName = "kIStStableFinalState"; + break; + case 2: + ParticleStatusName = "kIStIntermediateState"; + break; + case 3: + ParticleStatusName = "kIStDecayedState"; + break; + case 11: + ParticleStatusName = "kIStNucleonTarget"; + break; + case 12: + ParticleStatusName = "kIStDISPreFragmHadronicState"; + break; + case 13: + ParticleStatusName = "kIStPreDecayResonantState"; + break; + case 14: + ParticleStatusName = "kIStHadronInTheNucleus"; + break; + case 15: + ParticleStatusName = "kIStFinalStateNuclearRemnant"; + break; + case 16: + ParticleStatusName = "kIStNucleonClusterTarget"; + break; + default: + ParticleStatusName = "Status Unknown"; + } + return ParticleStatusName; + } + + //...................................................................... + std::string GENIEGenNGEM::ReactionChannel(int ccnc,int mode) + { + std::string ReactionChannelName=" "; + + if(ccnc==0) + ReactionChannelName = "kCC"; + else if(ccnc==1) + ReactionChannelName = "kNC"; + else std::cout<<"Current mode unknown!! "<Fill(mc.GetNeutrino().Mode()); + if (mc.GetNeutrino().Nu().PdgCode() == 12) id = 0; + else if (mc.GetNeutrino().Nu().PdgCode() == -12) id = 1; + else if (mc.GetNeutrino().Nu().PdgCode() == 14) id = 2; + else if (mc.GetNeutrino().Nu().PdgCode() == -14) id = 3; + else return; + } + else { + fNCMode->Fill(mc.GetNeutrino().Mode()); + if (mc.GetNeutrino().Nu().PdgCode() > 0) id = 4; + else id = 5; + } + if (id==-1) abort(); + + // Fill the specta histograms + fGenerated[id]->Fill(mc.GetNeutrino().Nu().E() ); + + ///< fill the vertex histograms from the neutrino - that is always + ///< particle 0 in the list + simb::MCNeutrino mcnu = mc.GetNeutrino(); + const simb::MCParticle nu = mcnu.Nu(); + + fVertexX->Fill(nu.Vx()); + fVertexY->Fill(nu.Vy()); + fVertexZ->Fill(nu.Vz()); + + fVertexXY->Fill(nu.Vx(), nu.Vy()); + fVertexXZ->Fill(nu.Vz(), nu.Vx()); + fVertexYZ->Fill(nu.Vz(), nu.Vy()); + + double mom = nu.P(); + if(std::abs(mom) > 0.){ + fDCosX->Fill(nu.Px()/mom); + fDCosY->Fill(nu.Py()/mom); + fDCosZ->Fill(nu.Pz()/mom); + } + + + MF_LOG_DEBUG("GENIEInteractionInformation") + << std::endl + << "REACTION: " << ReactionChannel(mc.GetNeutrino().CCNC(),mc.GetNeutrino().Mode()) + << std::endl + << "-----------> Particles in the Stack = " << mc.NParticles() << std::endl + << std::setiosflags(std::ios::left) + << std::setw(20) << "PARTICLE" + << std::setiosflags(std::ios::left) + << std::setw(32) << "STATUS" + << std::setw(18) << "E (GeV)" + << std::setw(18) << "m (GeV/c2)" + << std::setw(18) << "Ek (GeV)" + << std::endl << std::endl; + + const TDatabasePDG* databasePDG = TDatabasePDG::Instance(); + + // Loop over the particle stack for this event + for(int i = 0; i < mc.NParticles(); ++i){ + simb::MCParticle part(mc.GetParticle(i)); + std::string name = databasePDG->GetParticle(part.PdgCode())->GetName(); + int code = part.StatusCode(); + std::string status = ParticleStatus(code); + double mass = part.Mass(); + double energy = part.E(); + double Ek = (energy-mass); // Kinetic Energy (GeV) + if(status=="kIStStableFinalState"||status=="kIStHadronInTheNucleus") + MF_LOG_DEBUG("GENIEFinalState") + << std::setiosflags(std::ios::left) << std::setw(20) << name + << std::setiosflags(std::ios::left) << std::setw(32) <Fill(part.P()); + fEDCosX->Fill(part.Px()/part.P()); + fEDCosY->Fill(part.Py()/part.P()); + fEDCosZ->Fill(part.Pz()/part.P()); + fECons->Fill(nu.E() - part.E()); + break; + } + else if(abs(part.PdgCode()) == 13){ + fMuMomentum->Fill(part.P()); + fMuDCosX->Fill(part.Px()/part.P()); + fMuDCosY->Fill(part.Py()/part.P()); + fMuDCosZ->Fill(part.Pz()/part.P()); + fECons->Fill(nu.E() - part.E()); + break; + } + }// end loop over particles + }//end if CC interaction + + return; + } + +} + +namespace evgen{ + + DEFINE_ART_MODULE(GENIEGenNGEM) + +} + +#endif diff --git a/ubsim/NGEM/GenerateIsotropicSinglePhoton.h b/ubsim/NGEM/GenerateIsotropicSinglePhoton.h new file mode 100644 index 000000000..b88739384 --- /dev/null +++ b/ubsim/NGEM/GenerateIsotropicSinglePhoton.h @@ -0,0 +1,135 @@ +#ifndef EVGEN_GENERATE_ISOTROPIC_SINGLE_PHOTON_H +#define EVGEN_GENERATE_ISOTROPIC_SINGLE_PHOTON_H + +#include +#include +#include +#include + +// ROOT +#include "TRandom3.h" +#include "TLorentzVector.h" +#include "TMath.h" + +// LArSoft +#include "nusimdata/SimulationBase/MCTruth.h" +#include "nusimdata/SimulationBase/MCParticle.h" + +namespace evgen { + + // Utility: sample from PDF defined by bin edges and per-bin probabilities (weights) + inline double sampleFromBinned(const std::vector& binEdges, + const std::vector& binProbs, + TRandom3& randomGen) { + const std::size_t numBins = binProbs.size(); + + if (binEdges.size() != numBins + 1) { + throw std::runtime_error("sampleFromBinned: Number of bin edges must be equal to number of weights + 1"); + } + if (numBins == 0) { + throw std::runtime_error("sampleFromBinned: Number of bins must be greater than 0"); + } + + // check that the bin edges are in ascending order + for (std::size_t i = 0; i < numBins; ++i) { + if (binEdges[i] >= binEdges[i + 1]) { + throw std::runtime_error("sampleFromBinned: Bin edges must be in ascending order"); + } + } + + // check that the bin weights are non-negative + for (std::size_t i = 0; i < numBins; ++i) { + if (binProbs[i] < 0.0) { + throw std::runtime_error("sampleFromBinned: Bin weight must be greater than or equal to 0"); + } + } + + // Build cumulative probabilities (weights are per-bin probabilities; uniform within each bin) + std::vector cumulative; + cumulative.reserve(numBins + 1); + cumulative.push_back(0.0); + for (std::size_t i = 0; i < numBins; ++i) { + const double prob = binProbs[i]; + cumulative.push_back(cumulative.back() + prob); + } + const double totalProb = cumulative.back(); + const double tol = 1e-6; + if (std::abs(totalProb - 1.0) > tol) { + throw std::runtime_error("sampleFromBinned: Sum of bin probabilities must be 1 within tolerance, but is " + std::to_string(totalProb)); + } + + // Sample probability then locate bin + const double target = randomGen.Uniform(0.0, totalProb); + std::size_t binIndex = 0; + for (std::size_t i = 0; i < numBins; ++i) { + if (target < cumulative[i + 1]) { binIndex = i; break; } + } + + // Uniform within the selected bin + const double a = binEdges[binIndex]; + const double b = binEdges[binIndex + 1]; + return a + randomGen.Uniform(0.0, 1.0) * (b - a); + } + + // Replace the MCTruth particle list with a single photon using + // binned distributions for energy and polar angle (cosTheta). + inline void GenerateIsotropicSinglePhotonWithBins( + const std::vector& energyBinEdges, + const std::vector& energyBinProbs, + const std::vector& cosThetaBinEdges, + const std::vector& cosThetaBinProbs, + simb::MCTruth& originalMCTruth, + simb::MCTruth& newMCTruth) { + TRandom3 randomGen; + randomGen.SetSeed(0); + + const double E = sampleFromBinned(energyBinEdges, energyBinProbs, randomGen); + + const double cosTheta = sampleFromBinned(cosThetaBinEdges, cosThetaBinProbs, randomGen); + const double sinTheta = std::sqrt(std::max(0.0, 1.0 - cosTheta * cosTheta)); + const double phi = randomGen.Uniform(0.0, TMath::TwoPi()); + const double px = E * sinTheta * std::cos(phi); + const double py = E * sinTheta * std::sin(phi); + const double pz = E * cosTheta; + TLorentzVector momentum(px, py, pz, E); + + const simb::MCNeutrino neutrino = originalMCTruth.GetNeutrino(); + const simb::MCParticle nu = neutrino.Nu(); + TLorentzVector position(nu.Vx(), nu.Vy(), nu.Vz(), nu.T()); + + simb::MCParticle gamma(1, 22, "primary", 0, 0.0, 1); + gamma.AddTrajectoryPoint(position, momentum); + gamma.SetGvtx(nu.Gvx(), nu.Gvy(), nu.Gvz(), nu.Gvt()); + + newMCTruth.Add(gamma); + + // Preserve neutrino/meta information + if (originalMCTruth.NeutrinoSet()) { + simb::MCNeutrino neutrino = originalMCTruth.GetNeutrino(); + newMCTruth.SetNeutrino( + neutrino.CCNC(), neutrino.Mode(), neutrino.InteractionType(), + neutrino.Target(), neutrino.HitNuc(), neutrino.HitQuark(), + neutrino.W(), neutrino.X(), neutrino.Y(), neutrino.QSqr() + ); + } + newMCTruth.SetOrigin(originalMCTruth.Origin()); + } + + // In-place wrapper requiring explicit bins + inline void GenerateIsotropicSinglePhotonWithBins( + const std::vector& energyBinEdges, + const std::vector& energybinProbs, + const std::vector& cosThetaBinEdges, + const std::vector& cosThetabinProbs, + simb::MCTruth& mcTruth) { + simb::MCTruth newTruth; + GenerateIsotropicSinglePhotonWithBins( + energyBinEdges, energybinProbs, cosThetaBinEdges, cosThetabinProbs, + mcTruth, newTruth + ); + mcTruth = newTruth; + } + +} // namespace evgen + +#endif // EVGEN_GENERATE_ISOTROPIC_SINGLE_PHOTON_H diff --git a/ubsim/NGEM/ManuallyDecayPi0sToTwoPhotons.h b/ubsim/NGEM/ManuallyDecayPi0sToTwoPhotons.h new file mode 100644 index 000000000..83ea0d0cd --- /dev/null +++ b/ubsim/NGEM/ManuallyDecayPi0sToTwoPhotons.h @@ -0,0 +1,112 @@ +#ifndef EVGEN_MANUALLYDECAYPI0S_TO_TWO_PHOTONS_H +#define EVGEN_MANUALLYDECAYPI0S_TO_TWO_PHOTONS_H + +#include +#include +#include + +// ROOT +#include "TRandom3.h" +#include "TLorentzVector.h" +#include "TVector3.h" + +// LArSoft +#include "nusimdata/SimulationBase/MCTruth.h" +#include "nusimdata/SimulationBase/MCParticle.h" + +namespace evgen { + + // note that this does not consider rarer decays of pi0s + inline void ManuallyDecayPi0sToTwoPhotons(simb::MCTruth& originalMCTruth, simb::MCTruth& newMCTruth) { + for (int i = 0; i < originalMCTruth.NParticles(); ++i) { + TRandom3 randomGen; + randomGen.SetSeed(0); + + int pdgCode = originalMCTruth.GetParticle(i).PdgCode(); + int statusCode = originalMCTruth.GetParticle(i).StatusCode(); + if (pdgCode == 111 && statusCode == 1) { // pi0, GENIE status code 1 = kIstStableFinalState + std::cout << "Decaying pi0 at index " << i << std::endl; + const simb::MCParticle& pi0 = originalMCTruth.GetParticle(i); + + TLorentzVector pi0_position(pi0.Vx(), pi0.Vy(), pi0.Vz(), pi0.T()); + TLorentzVector pi0_momentum = pi0.Momentum(); + + // 3 = kIStDecayedState from GENIE + simb::MCParticle newPi0(pi0.TrackId(), pi0.PdgCode(), pi0.Process(), pi0.Mother(), pi0.Mass(), 3); + newPi0.AddTrajectoryPoint(pi0_position, pi0_momentum); + + // track_id, pdg, process, mother, mass, status_code + simb::MCParticle gamma1(979797971, 22, "primary", pi0.TrackId(), 0.0, 1); + simb::MCParticle gamma2(979797972, 22, "primary", pi0.TrackId(), 0.0, 1); + + // GENIE vertices relative to nucleus + gamma1.SetGvtx(pi0.Gvx(), pi0.Gvy(), pi0.Gvz(), pi0.Gvt()); + gamma2.SetGvtx(pi0.Gvx(), pi0.Gvy(), pi0.Gvz(), pi0.Gvt()); + + TVector3 pi0_boost_vector = pi0_momentum.BoostVector(); + + // Calculate decay momenta + double pi0_mass = pi0.Mass(); + double theta = randomGen.Uniform(0, 2 * M_PI); + double phi = acos(1 - 2 * randomGen.Uniform(0, 1)); + double pE = pi0_mass / 2; + double px = pE * sin(phi) * cos(theta); + double py = pE * sin(phi) * sin(theta); + double pz = pE * cos(phi); + TLorentzVector rest_frame_momentum_1(px, py, pz, pE); + TLorentzVector rest_frame_momentum_2(-px, -py, -pz, pE); + + // transform the momenta to the lab frame using a Lorentz boost + TLorentzVector lab_frame_momentum_1 = rest_frame_momentum_1; + TLorentzVector lab_frame_momentum_2 = rest_frame_momentum_2; + lab_frame_momentum_1.Boost(pi0_boost_vector); + lab_frame_momentum_2.Boost(pi0_boost_vector); + + gamma1.AddTrajectoryPoint(pi0_position, lab_frame_momentum_1); + gamma2.AddTrajectoryPoint(pi0_position, lab_frame_momentum_2); + + std::cout << "Conservation check:" << std::endl; + TLorentzVector sum = lab_frame_momentum_1 + lab_frame_momentum_2; + std::cout << "4-momentum difference (should be ~0): " + << (pi0_momentum - sum).Px() << ", " + << (pi0_momentum - sum).Py() << ", " + << (pi0_momentum - sum).Pz() << ", " + << (pi0_momentum - sum).E() << std::endl; + + newMCTruth.Add(newPi0); + newMCTruth.Add(gamma1); + newMCTruth.Add(gamma2); + } else { + // Copy over non-pi0 particles + const simb::MCParticle& particle = originalMCTruth.GetParticle(i); + simb::MCParticle non_const_particle = simb::MCParticle(particle); + newMCTruth.Add(non_const_particle); + } + } + + // Copy over the neutrino information + simb::MCNeutrino neutrino = originalMCTruth.GetNeutrino(); + int CCNC = neutrino.CCNC(); + int mode = neutrino.Mode(); + int interactionType = neutrino.InteractionType(); + int target = neutrino.Target(); + int nucleon = neutrino.HitNuc(); + int quark = neutrino.HitQuark(); + double w = neutrino.W(); + double x = neutrino.X(); + double y = neutrino.Y(); + double qsqr = neutrino.QSqr(); + newMCTruth.SetNeutrino(CCNC, mode, interactionType, target, nucleon, quark, w, x, y, qsqr); + newMCTruth.SetOrigin(originalMCTruth.Origin()); + } + + // In-place convenience wrapper + inline void ManuallyDecayPi0sToTwoPhotons(simb::MCTruth& mcTruth) { + simb::MCTruth newTruth; + ManuallyDecayPi0sToTwoPhotons(mcTruth, newTruth); + mcTruth = newTruth; + } + +} // namespace evgen + +#endif // EVGEN_MANUALLYDECAYPI0S_TO_TWO_PHOTONS_H diff --git a/ubsim/NGEM/mcc10_cmake.txt b/ubsim/NGEM/mcc10_cmake.txt new file mode 100644 index 000000000..7417acf7c --- /dev/null +++ b/ubsim/NGEM/mcc10_cmake.txt @@ -0,0 +1,42 @@ +# Optional: hide Genie logging +cet_add_compiler_flags(CXX -DHIDE_GENIE_LOG_XXX) + +# Plugin 1 +cet_build_plugin( + GENIEGenNGEM art::EDProducer + LIBRARIES + PRIVATE + lardataalg::MCDumpers + lardata::Utilities + larsim::Simulation + lardataobj::Simulation + nusimdata::SimulationBase + larcoreobj::SummaryData + larcorealg::Geometry + larcore::Geometry_Geometry_service + nurandom::RandomUtils_NuRandomService_service + nugen::EventGeneratorBase_GENIE + ${DK2NUDATA} + ${DK2NUGENIE} + dk2nu::Tree + dk2nu::Genie + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art::Framework_Services_Optional_RandomNumberGenerator_service + art::Persistency_Common + art::Persistency_Provenance + art::Utilities + art_root_io::TFileService_service + canvas::canvas + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + cetlib::cetlib + cetlib_except::cetlib_except + CLHEP::CLHEP + ${GENIE_LIB_LIST} + #ROOT::EGPYTHIA6 + ROOT::Tree + ROOT::EG + ${LOG4CPP} +) diff --git a/ubsim/NGEM/prodgenie_bnb_nu_del_1g_uboone.fcl b/ubsim/NGEM/prodgenie_bnb_nu_del_1g_uboone.fcl new file mode 100644 index 000000000..29f2e32eb --- /dev/null +++ b/ubsim/NGEM/prodgenie_bnb_nu_del_1g_uboone.fcl @@ -0,0 +1,28 @@ +#include "prodgenie_common_uboone_Filtered.fcl" + +process_name: GenieGenFiltered + +outputs.out1.fileName: "prodgenie_bnb_nu_filtered_del1g_uboone_%tc_gen.root" + +source.maxEvents: 750 + +physics.producers.generator: @local::microboone_genie_simple +physics.producers.generator.GlobalTimeOffset: 3125. +physics.producers.generator.RandomTimeOffset: 1600. +physics.producers.generator.TopVolume: "volCryostat" +physics.producers.generator.BeamName: "booster" + +physics.producers.generator.module_type: "GENIEGenNGEM" +physics.producers.generator.ManuallyDecayPi0s: true +physics.producers.generator.DeleteRandomGamma: true + +physics.filters.finalstatefilter.IsVerbose: false +physics.filters.finalstatefilter.isInclusive: true +physics.filters.finalstatefilter.PDG: [22] +physics.filters.finalstatefilter.PDGCount: [1] +physics.filters.finalstatefilter.PDGCountExclusivity: [true] + +# exactly one photon in the final state +# needs to have pi0s manually decayed in order to consider those photons +# needs to have photon deletion code in order to get many events + diff --git a/ubsim/NGEM/prodgenie_isotropic_1g_uboone.fcl b/ubsim/NGEM/prodgenie_isotropic_1g_uboone.fcl new file mode 100644 index 000000000..46af4d171 --- /dev/null +++ b/ubsim/NGEM/prodgenie_isotropic_1g_uboone.fcl @@ -0,0 +1,26 @@ +#include "standard_overlay_gen_griddriver.fcl" + +// this fhicl runs GENIE, but then throws away all the GENIE particles +// and replaces them with a new randomly generated single photon + +process_name: GenieGen + +outputs.out1.fileName: "prodgenie_bnb_nu_filtered_iso1g_uboone_%tc_gen.root" + +physics.producers.generator: @local::microboone_genie_simple +physics.producers.generator.GlobalTimeOffset: 3125. +physics.producers.generator.RandomTimeOffset: 1600. +physics.producers.generator.TopVolume: "volCryostat" +physics.producers.generator.BeamName: "booster" + +physics.producers.generator.module_type: "GENIEGenNGEM" +physics.producers.generator.GenerateIsotropicSinglePhoton: true + +// these bin edges and weights form the probability distributions for random sampling of the energy and angle +// bin edges must be increasing, bin probabalities must sum to 1 + +physics.producers.generator.SinglePhotonEnergyBinEdges: [0, 0.050, 0.100, 0.300, 0.800, 1.200, 2.000, 5.000] +physics.producers.generator.SinglePhotonEnergyBinProbs: [0.01, 0.02, 0.4, 0.3, 0.2, 0.06, 0.01] + +physics.producers.generator.SinglePhotonCosThetaBinEdges: [-1.0, 0.0, 0.5, 0.75, 0.9, 1.0] +physics.producers.generator.SinglePhotonCosThetaBinProbs: [0.05, 0.1, 0.2, 0.3, 0.35] From a97ef19432686b837c69e12fede4d02be982abe0 Mon Sep 17 00:00:00 2001 From: eyandel Date: Wed, 4 Mar 2026 10:30:05 -0600 Subject: [PATCH 02/19] clean up some bad CMake files --- ubsim/NGEM/CMake_bad2.txt | 38 -------------------------- ubsim/NGEM/CMake_noroot.txt | 54 ------------------------------------- ubsim/NGEM/mcc10_cmake.txt | 42 ----------------------------- 3 files changed, 134 deletions(-) delete mode 100644 ubsim/NGEM/CMake_bad2.txt delete mode 100644 ubsim/NGEM/CMake_noroot.txt delete mode 100644 ubsim/NGEM/mcc10_cmake.txt diff --git a/ubsim/NGEM/CMake_bad2.txt b/ubsim/NGEM/CMake_bad2.txt deleted file mode 100644 index 9b811a96f..000000000 --- a/ubsim/NGEM/CMake_bad2.txt +++ /dev/null @@ -1,38 +0,0 @@ -cet_build_plugin( - GENIEGenNGEM art::EDProducer - LIBRARIES - PRIVATE - ubsim::GENIEGenNGEM - fhiclcpp::fhiclcpp - lardata::Utilities - art_root_io::TFileService_service - nurandom::RandomUtils_NuRandomService_service -) - -cet_make_library( - SOURCE - LIBRARIES - PUBLIC - larsim::EventWeight_Base - larcore::Geometry_Geometry_service - nusimdata::SimulationBase - nugen::EventGeneratorBase_GENIE - Geant4::G4event - ppfx::ppfx - dk2nu::Tree - art_root_io::TFileService_service - fhiclcpp::fhiclcpp - geant4reweight::ReweightBaseLib - GENIE::GRwClc - GENIE::GFwInt - GENIE::GRwFwk - GENIE::GFwMsg - GENIE::GFwUtl - log4cpp::log4cpp - ROOT::Hist - ROOT::Tree -) - -install_headers() -install_fhicl() -install_source() \ No newline at end of file diff --git a/ubsim/NGEM/CMake_noroot.txt b/ubsim/NGEM/CMake_noroot.txt deleted file mode 100644 index 1226c29e0..000000000 --- a/ubsim/NGEM/CMake_noroot.txt +++ /dev/null @@ -1,54 +0,0 @@ - -cet_add_compiler_flags( CXX -DHIDE_GENIE_LOG_XXX ) - -include_directories ( $ENV{LOG4CPP_INC} ) -include_directories ( $ENV{DK2NUDATA_INC} ) -include_directories ( $ENV{DK2NUGENIE_INC} ) - -# genie -#cet_find_library( DK2NUDATA NAMES dk2nuTree PATHS ENV DK2NUDATA_LIB ) -#cet_find_library( DK2NUGENIE NAMES dk2nuGenie PATHS ENV DK2NUGENIE_LIB ) -find_library(DK2NUDATA dk2nuTree PATHS ENV DK2NUDATA_LIB REQUIRED) -find_library(DK2NUGENIE dk2nuGenie PATHS ENV DK2NUGENIE_LIB REQUIRED) - - -art_make( - MODULE_LIBRARIES - lardataalg_MCDumpers - larsim_Simulation lardataobj_Simulation - nusimdata_SimulationBase - larcoreobj_SummaryData - larcorealg_Geometry - larcore_Geometry_Geometry_service - nutools_RandomUtils_NuRandomService_service - nutools_EventGeneratorBase_GENIE - ${DK2NUDATA} - ${DK2NUGENIE} - ${ART_FRAMEWORK_CORE} - ${ART_FRAMEWORK_PRINCIPAL} - ${ART_FRAMEWORK_SERVICES_REGISTRY} - ${ART_FRAMEWORK_SERVICES_OPTIONAL} - ${ART_FRAMEWORK_SERVICES_OPTIONAL_TFILESERVICE_SERVICE} - ${ART_FRAMEWORK_SERVICES_OPTIONAL_RANDOMNUMBERGENERATOR_SERVICE} - art_Persistency_Common - art_Persistency_Provenance - art_Utilities - canvas - ${MF_MESSAGELOGGER} - ${FHICLCPP} - cetlib cetlib_except - ${CLHEP} - ${GENIE_LIB_LIST} - ${ROOT_EGPYTHIA6} # FIXME!!! - resolving genie run time reference - ${ROOT_BASIC_LIB_LIST} - ${ROOT_EG} - ${ROOT_GEOM} - ${ROOT_XMLIO} - ${ROOT_GDML} - ${LOG4CPP} - ) - -install_headers() -install_fhicl() -install_source() - diff --git a/ubsim/NGEM/mcc10_cmake.txt b/ubsim/NGEM/mcc10_cmake.txt deleted file mode 100644 index 7417acf7c..000000000 --- a/ubsim/NGEM/mcc10_cmake.txt +++ /dev/null @@ -1,42 +0,0 @@ -# Optional: hide Genie logging -cet_add_compiler_flags(CXX -DHIDE_GENIE_LOG_XXX) - -# Plugin 1 -cet_build_plugin( - GENIEGenNGEM art::EDProducer - LIBRARIES - PRIVATE - lardataalg::MCDumpers - lardata::Utilities - larsim::Simulation - lardataobj::Simulation - nusimdata::SimulationBase - larcoreobj::SummaryData - larcorealg::Geometry - larcore::Geometry_Geometry_service - nurandom::RandomUtils_NuRandomService_service - nugen::EventGeneratorBase_GENIE - ${DK2NUDATA} - ${DK2NUGENIE} - dk2nu::Tree - dk2nu::Genie - art::Framework_Core - art::Framework_Principal - art::Framework_Services_Registry - art::Framework_Services_Optional_RandomNumberGenerator_service - art::Persistency_Common - art::Persistency_Provenance - art::Utilities - art_root_io::TFileService_service - canvas::canvas - messagefacility::MF_MessageLogger - fhiclcpp::fhiclcpp - cetlib::cetlib - cetlib_except::cetlib_except - CLHEP::CLHEP - ${GENIE_LIB_LIST} - #ROOT::EGPYTHIA6 - ROOT::Tree - ROOT::EG - ${LOG4CPP} -) From 7bb3366f23614109b80d34744acaa54ca511de30 Mon Sep 17 00:00:00 2001 From: eyandel Date: Wed, 4 Mar 2026 12:00:42 -0600 Subject: [PATCH 03/19] start making particle modification code --- ubsim/NGEM/ModifyParticle.h | 188 ++++++++++++++++++++++ ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl | 25 +++ 2 files changed, 213 insertions(+) create mode 100644 ubsim/NGEM/ModifyParticle.h create mode 100644 ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl diff --git a/ubsim/NGEM/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h new file mode 100644 index 000000000..35ff9d6e5 --- /dev/null +++ b/ubsim/NGEM/ModifyParticle.h @@ -0,0 +1,188 @@ +#ifndef EVGEN_MODIFYPARTICLE_H +#define EVGEN_MODIFYPARTICLE_H + +#include +#include +#include + +// ROOT +#include "TRandom3.h" +#include "TLorentzVector.h" +#include "TVector3.h" + +// LArSoft +#include "nusimdata/SimulationBase/MCTruth.h" +#include "nusimdata/SimulationBase/MCParticle.h" + +namespace evgen { + + // Utility: sample from PDF defined by bin edges and per-bin probabilities (weights) + inline double sampleFromBinned(const std::vector& binEdges, + const std::vector& binProbs, + TRandom3& randomGen) { + const std::size_t numBins = binProbs.size(); + + if (binEdges.size() != numBins + 1) { + throw std::runtime_error("sampleFromBinned: Number of bin edges must be equal to number of weights + 1"); + } + if (numBins == 0) { + throw std::runtime_error("sampleFromBinned: Number of bins must be greater than 0"); + } + + // check that the bin edges are in ascending order + for (std::size_t i = 0; i < numBins; ++i) { + if (binEdges[i] >= binEdges[i + 1]) { + throw std::runtime_error("sampleFromBinned: Bin edges must be in ascending order"); + } + } + + // check that the bin weights are non-negative + for (std::size_t i = 0; i < numBins; ++i) { + if (binProbs[i] < 0.0) { + throw std::runtime_error("sampleFromBinned: Bin weight must be greater than or equal to 0"); + } + } + + // Build cumulative probabilities (weights are per-bin probabilities; uniform within each bin) + std::vector cumulative; + cumulative.reserve(numBins + 1); + cumulative.push_back(0.0); + for (std::size_t i = 0; i < numBins; ++i) { + const double prob = binProbs[i]; + cumulative.push_back(cumulative.back() + prob); + } + const double totalProb = cumulative.back(); + const double tol = 1e-6; + if (std::abs(totalProb - 1.0) > tol) { + throw std::runtime_error("sampleFromBinned: Sum of bin probabilities must be 1 within tolerance, but is " + std::to_string(totalProb)); + } + + // Sample probability then locate bin + const double target = randomGen.Uniform(0.0, totalProb); + std::size_t binIndex = 0; + for (std::size_t i = 0; i < numBins; ++i) { + if (target < cumulative[i + 1]) { binIndex = i; break; } + } + + // Uniform within the selected bin + const double a = binEdges[binIndex]; + const double b = binEdges[binIndex + 1]; + return a + randomGen.Uniform(0.0, 1.0) * (b - a); + } + + // modify particle kinematics + inline void ModifyParticle(int modPdg, std::string modVar, + const std::vector& varBinEdges, const std::vector& varBinProbs, + simb::MCTruth& originalMCTruth, simb::MCTruth& newMCTruth, + int seed = 0) { + for (int i = 0; i < originalMCTruth.NParticles(); ++i) { + TRandom3 randomGen; + randomGen.SetSeed(seed); + + int pdgCode = originalMCTruth.GetParticle(i).PdgCode(); + int statusCode = originalMCTruth.GetParticle(i).StatusCode(); + if (pdgCode == modPdg && statusCode == 1) { // pi0, GENIE status code 1 = kIstStableFinalState + std::cout << "Modifying " << modVar << " of particle with pdg code " << modPdg << " at index " << i << std::endl; + const simb::MCParticle& part = originalMCTruth.GetParticle(i); + + TLorentzVector part_position(part.Vx(), part.Vy(), part.Vz(), part.T()); + TLorentzVector part_momentum = part.Momentum(); + + // GENIE vertices relative to nucleus + //Gvtx(part.Gvx(), part.Gvy(), part.Gvz(), part.Gvt()); + + TVector3 part_boost_vector = part_momentum.BoostVector(); + + // Calculate decay momenta + double part_mass = part.Mass(); + //double theta = randomGen.Uniform(0, 2 * M_PI); + //double phi = acos(1 - 2 * randomGen.Uniform(0, 1)); + double pE = part_mass / 2; + double px = pE * sin(phi) * cos(theta); + double py = pE * sin(phi) * sin(theta); + double pz = pE * cos(phi); + TLorentzVector rest_frame_momentum(px, py, pz, pE); + + // transform the momenta to the lab frame using a Lorentz boost + TLorentzVector lab_frame_momentum = rest_frame_momentum; + lab_frame_momentum.Boost(part_boost_vector); + + auto new_part_process = part.Process(); + auto new_part_mother = part.Mother(); + auto new_part_mass = part.Mass(); + auto new_part_position = part_position; + auto new_part_momentum = part_momentum; + + auto new_part_var = sampleFromBinned(varBinEdges, varBinProbs, randomGen); + + if (modVar == "mass"){ + new_part_mass = new_part_var; + } else if (modVar == "E") { + new_part_momentum.SetE(new_part_var); + } else if (modVar == "px") { + new_part_momentum.SetPx(new_part_var); + } else if (modVar == "py") { + new_part_momentum.SetPy(new_part_var); + } else if (modVar == "pz") { + new_part_momentum.SetPz(new_part_var); + } else if (modVar == "theta") { + double p = new_part_momentum.P(); + double phi = new_part_momentum.Phi(); + double px = p * sin(new_part_var) * cos(phi); + double py = p * sin(new_part_var) * sin(phi); + double pz = p * cos(new_part_var); + new_part_momentum.SetPx(px); + new_part_momentum.SetPy(py); + new_part_momentum.SetPz(pz); + } else if (modVar == "phi") { + double p = new_part_momentum.P(); + double theta = new_part_momentum.Theta(); + double px = p * sin(theta) * cos(new_part_var); + double py = p * sin(theta) * sin(new_part_var); + double pz = p * cos(theta); + new_part_momentum.SetPx(px); + new_part_momentum.SetPy(py); + new_part_momentum.SetPz(pz); + } else { + throw std::runtime_error("ModifyParticle: Invalid modVar " + modVar); + } + + // make new Part particle with modifications + simb::MCParticle newPart(part.TrackId(), part.PdgCode(), part.Process(), part.Mother(), part.Mass(), part.StatusCode()); + newPart.AddTrajectoryPoint(part_position, part_momentum); + + newMCTruth.Add(newPart); + } else { + // Copy over other particles + const simb::MCParticle& particle = originalMCTruth.GetParticle(i); + simb::MCParticle non_const_particle = simb::MCParticle(particle); + newMCTruth.Add(non_const_particle); + } + } + + // Copy over the neutrino information + simb::MCNeutrino neutrino = originalMCTruth.GetNeutrino(); + int CCNC = neutrino.CCNC(); + int mode = neutrino.Mode(); + int interactionType = neutrino.InteractionType(); + int target = neutrino.Target(); + int nucleon = neutrino.HitNuc(); + int quark = neutrino.HitQuark(); + double w = neutrino.W(); + double x = neutrino.X(); + double y = neutrino.Y(); + double qsqr = neutrino.QSqr(); + newMCTruth.SetNeutrino(CCNC, mode, interactionType, target, nucleon, quark, w, x, y, qsqr); + newMCTruth.SetOrigin(originalMCTruth.Origin()); + } + + // In-place convenience wrapper + inline void ModifyParticle(simb::MCTruth& mcTruth) { + simb::MCTruth newTruth; + ModifyParticle(mcTruth, newTruth); + mcTruth = newTruth; + } + +} // namespace evgen + +#endif // EVGEN_MODIFYPARTICLE_H diff --git a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl new file mode 100644 index 000000000..f7223b454 --- /dev/null +++ b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl @@ -0,0 +1,25 @@ +#include "prodgenie_common_uboone_Filtered.fcl" + +process_name: GenieGenFiltered + +outputs.out1.fileName: "prodgenie_bnb_nu_filtered_del1g_uboone_%tc_gen.root" + +source.maxEvents: 750 + +physics.producers.generator: @local::microboone_genie_simple +physics.producers.generator.GlobalTimeOffset: 3125. +physics.producers.generator.RandomTimeOffset: 1600. +physics.producers.generator.TopVolume: "volCryostat" +physics.producers.generator.BeamName: "booster" + +physics.producers.generator.module_type: "GENIEGenNGEM" +physics.producers.generator.ManuallyDecayPi0s: false +physics.producers.generator.DeleteRandomGamma: false +physics.producers.generator.GenerateIsotropicSinglePhoton: false +physics.producers.generator.ModifyPi0s: true + +physics.filters.finalstatefilter.IsVerbose: false +physics.filters.finalstatefilter.isInclusive: true +physics.filters.finalstatefilter.PDG: [111] #CC and NC PiZero +physics.filters.finalstatefilter.PDGCount: [1] +physics.filters.finalstatefilter.PDGCountExclusivity: [true] \ No newline at end of file From 7d81f3f57589be3949b46ab0fa8be210aa924939 Mon Sep 17 00:00:00 2001 From: eyandel Date: Thu, 5 Mar 2026 10:42:03 -0600 Subject: [PATCH 04/19] implement modify particle code --- ubsim/NGEM/GENIEGenNGEM_module.cc | 43 +++++++++++++ ubsim/NGEM/ModifyParticle.h | 77 +++++++++++++---------- ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl | 10 ++- 3 files changed, 96 insertions(+), 34 deletions(-) diff --git a/ubsim/NGEM/GENIEGenNGEM_module.cc b/ubsim/NGEM/GENIEGenNGEM_module.cc index c7869888e..a27735c20 100644 --- a/ubsim/NGEM/GENIEGenNGEM_module.cc +++ b/ubsim/NGEM/GENIEGenNGEM_module.cc @@ -171,6 +171,14 @@ namespace evgen { std::vector SinglePhotonCosThetaBinEdges; ///< configurable via FHiCL std::vector SinglePhotonCosThetaBinProbs; ///< configurable via FHiCL + bool ModifyParticle; ///< whether to modify particles according to configuration parameters + int ModifyParticlePdg; ///< pdg code of particle to modify + std::vector ModifyParticleVar; ///< variables to modify + //(options: "E", "px", "py", "pz", "mass", "x", "y", "z", "T", "theta", "phi", "p"); configurable via FHiCL + std::vector> ModifyParticleVarBinEdges; ///< configurable via FHiCL + std::vector> ModifyParticleVarBinProbs; ///< configurable via FHiCL + int ModifyParticleRandomSeed; ///< random seed for particle modification; configurable via FHiCL + }; } @@ -192,6 +200,12 @@ namespace evgen{ , SinglePhotonEnergyBinProbs(pset.get>("SinglePhotonEnergyBinProbs", std::vector())) , SinglePhotonCosThetaBinEdges(pset.get>("SinglePhotonCosThetaBinEdges", std::vector())) , SinglePhotonCosThetaBinProbs(pset.get>("SinglePhotonCosThetaBinProbs", std::vector())) + , ModifyParticle(pset.get("ModifyParticle", false)) + , ModifyParticlePdg(pset.get("ModifyParticlePdg", 0)) + , ModifyParticleVar(pset.get>("ModifyParticleVar", std::vector())) + , ModifyParticleVarBinEdges(pset.get>>("ModifyParticleVarBinEdges", std::vector>())) + , ModifyParticleVarBinProbs(pset.get>>("ModifyParticleVarBinProbs", std::vector>())) + , ModifyParticleRandomSeed(pset.get("ModifyParticleRandomSeed", 0)) { fStopwatch.Start(); @@ -477,6 +491,35 @@ namespace evgen{ } } + if (ModifyParticle) { + std::cout << "Modifying particles with pdg code " << ModifyParticlePdg << std::endl; + if (ModifyParticleVar.size() != ModifyParticleVarBinEdges.size() || ModifyParticleVar.size() != ModifyParticleVarBinProbs.size()) { + std::cerr << "Error: ModifyParticleVar, ModifyParticleVarBinEdges, and ModifyParticleVarBinProbs must have the same size" << std::endl; + return; + } + for (int i_v = 0; i_v < ModifyParticleVar.size(); ++i_v) { + std::cout << " Modifying variable " << ModifyParticleVar.at(i_v) << std::endl; + ModifyParticle( + ModifyParticlePdg, ModifyParticleVar.at(i_v), + ModifyParticleVarBinEdges.at(i_v), ModifyParticleVarBinProbs.at(i_v), + ModifyParticleRandomSeed, truth + ); + } + std::cout << "truth.NParticles(): " << truth.NParticles() << std::endl; + for (int i = 0; i < truth.NParticles(); ++i) { + std::cout << " pdg: " << truth.GetParticle(i).PdgCode(); + std::cout << ", track_id: " << truth.GetParticle(i).TrackId(); + std::cout << ", mother: " << truth.GetParticle(i).Mother(); + std::cout << ", mass: " << truth.GetParticle(i).Mass(); + std::cout << ", position: (" << truth.GetParticle(i).Vx() << ", " << truth.GetParticle(i).Vy() << ", " << truth.GetParticle(i).Vz() << ")"; + std::cout << ", momentum: (" << truth.GetParticle(i).Px() << ", " << truth.GetParticle(i).Py() << ", " << truth.GetParticle(i).Pz() << ")"; + std::cout << ", Gvtx: (" << truth.GetParticle(i).Gvx() << ", " << truth.GetParticle(i).Gvy() << ", " << truth.GetParticle(i).Gvz() << ", " << truth.GetParticle(i).Gvt() << ")"; + std::cout << ", status code: " << truth.GetParticle(i).StatusCode(); + std::cout << ", process: " << truth.GetParticle(i).Process(); + std::cout << std::endl; + } + } + truthcol ->push_back(truth); fluxcol ->push_back(flux); gtruthcol->push_back(gTruth); diff --git a/ubsim/NGEM/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h index 35ff9d6e5..74316b699 100644 --- a/ubsim/NGEM/ModifyParticle.h +++ b/ubsim/NGEM/ModifyParticle.h @@ -72,16 +72,15 @@ namespace evgen { // modify particle kinematics inline void ModifyParticle(int modPdg, std::string modVar, - const std::vector& varBinEdges, const std::vector& varBinProbs, - simb::MCTruth& originalMCTruth, simb::MCTruth& newMCTruth, - int seed = 0) { + const std::vector& varBinEdges, const std::vector& varBinProbs, int seed + simb::MCTruth& originalMCTruth, simb::MCTruth& newMCTruth) { for (int i = 0; i < originalMCTruth.NParticles(); ++i) { TRandom3 randomGen; randomGen.SetSeed(seed); int pdgCode = originalMCTruth.GetParticle(i).PdgCode(); int statusCode = originalMCTruth.GetParticle(i).StatusCode(); - if (pdgCode == modPdg && statusCode == 1) { // pi0, GENIE status code 1 = kIstStableFinalState + if (pdgCode == modPdg && statusCode == 1) { // GENIE status code 1 = kIstStableFinalState std::cout << "Modifying " << modVar << " of particle with pdg code " << modPdg << " at index " << i << std::endl; const simb::MCParticle& part = originalMCTruth.GetParticle(i); @@ -97,18 +96,21 @@ namespace evgen { double part_mass = part.Mass(); //double theta = randomGen.Uniform(0, 2 * M_PI); //double phi = acos(1 - 2 * randomGen.Uniform(0, 1)); - double pE = part_mass / 2; - double px = pE * sin(phi) * cos(theta); - double py = pE * sin(phi) * sin(theta); - double pz = pE * cos(phi); - TLorentzVector rest_frame_momentum(px, py, pz, pE); + double pE = part_momentum.E(); //part_mass / 2; + double px = part_momentum.Px(); //pE * sin(phi) * cos(theta); + double py = part_momentum.Py(); //pE * sin(phi) * sin(theta); + double pz = part_momentum.Pz(); //pE * cos(phi); + //TLorentzVector rest_frame_momentum(px, py, pz, pE); + + double T = part_position.T(); + double x = part_position.Vx(); + double y = part_position.Vy(); + double z = part_position.Vz(); // transform the momenta to the lab frame using a Lorentz boost - TLorentzVector lab_frame_momentum = rest_frame_momentum; - lab_frame_momentum.Boost(part_boost_vector); + //TLorentzVector lab_frame_momentum = rest_frame_momentum; + //lab_frame_momentum.Boost(part_boost_vector); - auto new_part_process = part.Process(); - auto new_part_mother = part.Mother(); auto new_part_mass = part.Mass(); auto new_part_position = part_position; auto new_part_momentum = part_momentum; @@ -118,44 +120,55 @@ namespace evgen { if (modVar == "mass"){ new_part_mass = new_part_var; } else if (modVar == "E") { - new_part_momentum.SetE(new_part_var); + new_part_momentum = TLorentzVector mom(px, py, pz, new_part_var); } else if (modVar == "px") { - new_part_momentum.SetPx(new_part_var); + new_part_momentum = TLorentzVector mom(new_part_var, py, pz, pE); } else if (modVar == "py") { - new_part_momentum.SetPy(new_part_var); + new_part_momentum = TLorentzVector mom(px, new_part_var, pz, pE); } else if (modVar == "pz") { - new_part_momentum.SetPz(new_part_var); + new_part_momentum = TLorentzVector mom(px, py, new_part_var, pE); } else if (modVar == "theta") { double p = new_part_momentum.P(); double phi = new_part_momentum.Phi(); - double px = p * sin(new_part_var) * cos(phi); - double py = p * sin(new_part_var) * sin(phi); - double pz = p * cos(new_part_var); - new_part_momentum.SetPx(px); - new_part_momentum.SetPy(py); - new_part_momentum.SetPz(pz); + double new_px = p * sin(new_part_var) * cos(phi); + double new_py = p * sin(new_part_var) * sin(phi); + double new_pz = p * cos(new_part_var); + new_part_momentum = TLorentzVector mom(new_px, new_py, new_pz, pE); } else if (modVar == "phi") { double p = new_part_momentum.P(); double theta = new_part_momentum.Theta(); - double px = p * sin(theta) * cos(new_part_var); - double py = p * sin(theta) * sin(new_part_var); - double pz = p * cos(theta); - new_part_momentum.SetPx(px); - new_part_momentum.SetPy(py); - new_part_momentum.SetPz(pz); + double new_px = p * sin(theta) * cos(new_part_var); + double new_py = p * sin(theta) * sin(new_part_var); + double new_pz = p * cos(theta); + new_part_momentum = TLorentzVector mom(new_px, new_py, new_pz, pE); + } else if (modVar == "p") { + double theta = new_part_momentum.Theta(); + double phi = new_part_momentum.Phi(); + double new_px = new_part_var * sin(theta) * cos(phi); + double new_py = new_part_var * sin(theta) * sin(phi); + double new_pz = new_part_var * cos(theta); + new_part_momentum = TLorentzVector mom(new_px, new_py, new_pz, pE); + } else if ( modVar == "x") { + new_part_position = TLorentzVector(new_part_var, y, z, T); + } else if ( modVar == "y") { + new_part_position = TLorentzVector(x, new_part_var, z, T); + } else if ( modVar == "z") { + new_part_position = TLorentzVector(x, y, new_part_var, T); + } else if ( modVar == "T") { + new_part_position = TLorentzVector(x, y, z, new_part_var); } else { throw std::runtime_error("ModifyParticle: Invalid modVar " + modVar); } // make new Part particle with modifications - simb::MCParticle newPart(part.TrackId(), part.PdgCode(), part.Process(), part.Mother(), part.Mass(), part.StatusCode()); - newPart.AddTrajectoryPoint(part_position, part_momentum); + simb::MCParticle newPart(part.TrackId(), part.PdgCode(), part.Process(), part.Mother(), new_part_mass, part.StatusCode()); + newPart.AddTrajectoryPoint(new_part_position, new_part_momentum); newMCTruth.Add(newPart); } else { // Copy over other particles const simb::MCParticle& particle = originalMCTruth.GetParticle(i); - simb::MCParticle non_const_particle = simb::MCParticle(particle); + simb::MCParticle non_const_particle = simb::MCParticle(particle); newMCTruth.Add(non_const_particle); } } diff --git a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl index f7223b454..678e2d242 100644 --- a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl +++ b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl @@ -2,7 +2,7 @@ process_name: GenieGenFiltered -outputs.out1.fileName: "prodgenie_bnb_nu_filtered_del1g_uboone_%tc_gen.root" +outputs.out1.fileName: "prodgenie_bnb_nu_filtered_modpi0_uboone_%tc_gen.root" source.maxEvents: 750 @@ -16,7 +16,13 @@ physics.producers.generator.module_type: "GENIEGenNGEM" physics.producers.generator.ManuallyDecayPi0s: false physics.producers.generator.DeleteRandomGamma: false physics.producers.generator.GenerateIsotropicSinglePhoton: false -physics.producers.generator.ModifyPi0s: true + +physics.producers.generator.ModifyParticle: true +physics.producers.generator.ModifyParticlePdg: 111 # pi0 +physics.producers.generator.ModifyParticleVar: ["mass"] # modify mass of pi0s +physics.producers.generator.ModifyParticleVarBinEdges: [[0.001, 2.0, 5.0]] # bin edges for mass modification (in GeV) +physics.producers.generator.ModifyParticleVarBinProbs: [[0.85, 0.15]] # probabilities for each mass bin +physics.producers.generator.ModifyParticleRandomSeed: 5 # random seed for particle modification, default is 0 (no seed) physics.filters.finalstatefilter.IsVerbose: false physics.filters.finalstatefilter.isInclusive: true From b80d7e4179c29e7b790c7e95d54e03bf91524139 Mon Sep 17 00:00:00 2001 From: eyandel Date: Thu, 5 Mar 2026 10:43:40 -0600 Subject: [PATCH 05/19] change mass range for testing --- ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl index 678e2d242..402c437a8 100644 --- a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl +++ b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl @@ -20,8 +20,8 @@ physics.producers.generator.GenerateIsotropicSinglePhoton: false physics.producers.generator.ModifyParticle: true physics.producers.generator.ModifyParticlePdg: 111 # pi0 physics.producers.generator.ModifyParticleVar: ["mass"] # modify mass of pi0s -physics.producers.generator.ModifyParticleVarBinEdges: [[0.001, 2.0, 5.0]] # bin edges for mass modification (in GeV) -physics.producers.generator.ModifyParticleVarBinProbs: [[0.85, 0.15]] # probabilities for each mass bin +physics.producers.generator.ModifyParticleVarBinEdges: [[0.0, 1000.0]] # bin edges for mass modification (in GeV) +physics.producers.generator.ModifyParticleVarBinProbs: [[1.00]] # probabilities for each mass bin physics.producers.generator.ModifyParticleRandomSeed: 5 # random seed for particle modification, default is 0 (no seed) physics.filters.finalstatefilter.IsVerbose: false From 08a2f5b029192672f7241078bdc3edf5a33cd7f8 Mon Sep 17 00:00:00 2001 From: eyandel Date: Tue, 17 Mar 2026 12:26:03 -0500 Subject: [PATCH 06/19] compiles now --- ubsim/NGEM/GENIEGenNGEM_module.cc | 43 ++++++------- ubsim/NGEM/ModifyParticle.h | 100 +++++++++--------------------- 2 files changed, 49 insertions(+), 94 deletions(-) diff --git a/ubsim/NGEM/GENIEGenNGEM_module.cc b/ubsim/NGEM/GENIEGenNGEM_module.cc index a27735c20..b2834ecb4 100644 --- a/ubsim/NGEM/GENIEGenNGEM_module.cc +++ b/ubsim/NGEM/GENIEGenNGEM_module.cc @@ -67,6 +67,7 @@ #include "ManuallyDecayPi0sToTwoPhotons.h" #include "DeleteOneRandomPhoton.h" #include "GenerateIsotropicSinglePhoton.h" +#include "ModifyParticle.h" ///Event Generation using GENIE, cosmics or single particles namespace evgen { @@ -171,13 +172,13 @@ namespace evgen { std::vector SinglePhotonCosThetaBinEdges; ///< configurable via FHiCL std::vector SinglePhotonCosThetaBinProbs; ///< configurable via FHiCL - bool ModifyParticle; ///< whether to modify particles according to configuration parameters - int ModifyParticlePdg; ///< pdg code of particle to modify - std::vector ModifyParticleVar; ///< variables to modify + bool fModifyParticle; ///< whether to modify particles according to configuration parameters + int fModifyParticlePdg; ///< pdg code of particle to modify + std::vector fModifyParticleVar; ///< variables to modify //(options: "E", "px", "py", "pz", "mass", "x", "y", "z", "T", "theta", "phi", "p"); configurable via FHiCL - std::vector> ModifyParticleVarBinEdges; ///< configurable via FHiCL - std::vector> ModifyParticleVarBinProbs; ///< configurable via FHiCL - int ModifyParticleRandomSeed; ///< random seed for particle modification; configurable via FHiCL + std::vector> fModifyParticleVarBinEdges; ///< configurable via FHiCL + std::vector> fModifyParticleVarBinProbs; ///< configurable via FHiCL + int fModifyParticleRandomSeed; ///< random seed for particle modification; configurable via FHiCL }; } @@ -200,12 +201,12 @@ namespace evgen{ , SinglePhotonEnergyBinProbs(pset.get>("SinglePhotonEnergyBinProbs", std::vector())) , SinglePhotonCosThetaBinEdges(pset.get>("SinglePhotonCosThetaBinEdges", std::vector())) , SinglePhotonCosThetaBinProbs(pset.get>("SinglePhotonCosThetaBinProbs", std::vector())) - , ModifyParticle(pset.get("ModifyParticle", false)) - , ModifyParticlePdg(pset.get("ModifyParticlePdg", 0)) - , ModifyParticleVar(pset.get>("ModifyParticleVar", std::vector())) - , ModifyParticleVarBinEdges(pset.get>>("ModifyParticleVarBinEdges", std::vector>())) - , ModifyParticleVarBinProbs(pset.get>>("ModifyParticleVarBinProbs", std::vector>())) - , ModifyParticleRandomSeed(pset.get("ModifyParticleRandomSeed", 0)) + , fModifyParticle(pset.get("ModifyParticle", false)) + , fModifyParticlePdg(pset.get("ModifyParticlePdg", 0)) + , fModifyParticleVar(pset.get>("ModifyParticleVar", std::vector())) + , fModifyParticleVarBinEdges(pset.get>>("ModifyParticleVarBinEdges", std::vector>())) + , fModifyParticleVarBinProbs(pset.get>>("ModifyParticleVarBinProbs", std::vector>())) + , fModifyParticleRandomSeed(pset.get("ModifyParticleRandomSeed", 0)) { fStopwatch.Start(); @@ -491,19 +492,15 @@ namespace evgen{ } } - if (ModifyParticle) { - std::cout << "Modifying particles with pdg code " << ModifyParticlePdg << std::endl; - if (ModifyParticleVar.size() != ModifyParticleVarBinEdges.size() || ModifyParticleVar.size() != ModifyParticleVarBinProbs.size()) { - std::cerr << "Error: ModifyParticleVar, ModifyParticleVarBinEdges, and ModifyParticleVarBinProbs must have the same size" << std::endl; + if (fModifyParticle) { + std::cout << "Modifying particles with pdg code " << fModifyParticlePdg << std::endl; + if (fModifyParticleVar.size() != fModifyParticleVarBinEdges.size() || fModifyParticleVar.size() != fModifyParticleVarBinProbs.size()) { + std::cerr << "Error: fModifyParticleVar, fModifyParticleVarBinEdges, and fModifyParticleVarBinProbs must have the same size" << std::endl; return; } - for (int i_v = 0; i_v < ModifyParticleVar.size(); ++i_v) { - std::cout << " Modifying variable " << ModifyParticleVar.at(i_v) << std::endl; - ModifyParticle( - ModifyParticlePdg, ModifyParticleVar.at(i_v), - ModifyParticleVarBinEdges.at(i_v), ModifyParticleVarBinProbs.at(i_v), - ModifyParticleRandomSeed, truth - ); + for (uint i_v = 0; i_v < fModifyParticleVar.size(); ++i_v) { + std::cout << " Modifying variable " << fModifyParticleVar.at(i_v) << std::endl; + ModifyParticle(fModifyParticlePdg, fModifyParticleVar.at(i_v), fModifyParticleVarBinEdges.at(i_v), fModifyParticleVarBinProbs.at(i_v), fModifyParticleRandomSeed, truth); } std::cout << "truth.NParticles(): " << truth.NParticles() << std::endl; for (int i = 0; i < truth.NParticles(); ++i) { diff --git a/ubsim/NGEM/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h index 74316b699..e3d03d44c 100644 --- a/ubsim/NGEM/ModifyParticle.h +++ b/ubsim/NGEM/ModifyParticle.h @@ -14,66 +14,15 @@ #include "nusimdata/SimulationBase/MCTruth.h" #include "nusimdata/SimulationBase/MCParticle.h" -namespace evgen { - - // Utility: sample from PDF defined by bin edges and per-bin probabilities (weights) - inline double sampleFromBinned(const std::vector& binEdges, - const std::vector& binProbs, - TRandom3& randomGen) { - const std::size_t numBins = binProbs.size(); - - if (binEdges.size() != numBins + 1) { - throw std::runtime_error("sampleFromBinned: Number of bin edges must be equal to number of weights + 1"); - } - if (numBins == 0) { - throw std::runtime_error("sampleFromBinned: Number of bins must be greater than 0"); - } - - // check that the bin edges are in ascending order - for (std::size_t i = 0; i < numBins; ++i) { - if (binEdges[i] >= binEdges[i + 1]) { - throw std::runtime_error("sampleFromBinned: Bin edges must be in ascending order"); - } - } - - // check that the bin weights are non-negative - for (std::size_t i = 0; i < numBins; ++i) { - if (binProbs[i] < 0.0) { - throw std::runtime_error("sampleFromBinned: Bin weight must be greater than or equal to 0"); - } - } - - // Build cumulative probabilities (weights are per-bin probabilities; uniform within each bin) - std::vector cumulative; - cumulative.reserve(numBins + 1); - cumulative.push_back(0.0); - for (std::size_t i = 0; i < numBins; ++i) { - const double prob = binProbs[i]; - cumulative.push_back(cumulative.back() + prob); - } - const double totalProb = cumulative.back(); - const double tol = 1e-6; - if (std::abs(totalProb - 1.0) > tol) { - throw std::runtime_error("sampleFromBinned: Sum of bin probabilities must be 1 within tolerance, but is " + std::to_string(totalProb)); - } - - // Sample probability then locate bin - const double target = randomGen.Uniform(0.0, totalProb); - std::size_t binIndex = 0; - for (std::size_t i = 0; i < numBins; ++i) { - if (target < cumulative[i + 1]) { binIndex = i; break; } - } +#include "GenerateIsotropicSinglePhoton.h" - // Uniform within the selected bin - const double a = binEdges[binIndex]; - const double b = binEdges[binIndex + 1]; - return a + randomGen.Uniform(0.0, 1.0) * (b - a); - } +namespace evgen { // modify particle kinematics inline void ModifyParticle(int modPdg, std::string modVar, - const std::vector& varBinEdges, const std::vector& varBinProbs, int seed - simb::MCTruth& originalMCTruth, simb::MCTruth& newMCTruth) { + const std::vector& varBinEdges, const std::vector& varBinProbs, + int seed, simb::MCTruth& originalMCTruth, simb::MCTruth& newMCTruth) { + for (int i = 0; i < originalMCTruth.NParticles(); ++i) { TRandom3 randomGen; randomGen.SetSeed(seed); @@ -102,16 +51,16 @@ namespace evgen { double pz = part_momentum.Pz(); //pE * cos(phi); //TLorentzVector rest_frame_momentum(px, py, pz, pE); - double T = part_position.T(); - double x = part_position.Vx(); - double y = part_position.Vy(); - double z = part_position.Vz(); + double T = part.T(); + double x = part.Vx(); + double y = part.Vy(); + double z = part.Vz(); // transform the momenta to the lab frame using a Lorentz boost //TLorentzVector lab_frame_momentum = rest_frame_momentum; //lab_frame_momentum.Boost(part_boost_vector); - auto new_part_mass = part.Mass(); + auto new_part_mass = part_mass; auto new_part_position = part_position; auto new_part_momentum = part_momentum; @@ -119,35 +68,42 @@ namespace evgen { if (modVar == "mass"){ new_part_mass = new_part_var; - } else if (modVar == "E") { - new_part_momentum = TLorentzVector mom(px, py, pz, new_part_var); + } else if (modVar == "E") { + TLorentzVector mom(px, py, pz, new_part_var); + new_part_momentum = mom; } else if (modVar == "px") { - new_part_momentum = TLorentzVector mom(new_part_var, py, pz, pE); + TLorentzVector mom(new_part_var, py, pz, pE); + new_part_momentum = mom; } else if (modVar == "py") { - new_part_momentum = TLorentzVector mom(px, new_part_var, pz, pE); + TLorentzVector mom(px, new_part_var, pz, pE); + new_part_momentum = mom; } else if (modVar == "pz") { - new_part_momentum = TLorentzVector mom(px, py, new_part_var, pE); + TLorentzVector mom(px, py, new_part_var, pE); + new_part_momentum = mom; } else if (modVar == "theta") { double p = new_part_momentum.P(); double phi = new_part_momentum.Phi(); double new_px = p * sin(new_part_var) * cos(phi); double new_py = p * sin(new_part_var) * sin(phi); double new_pz = p * cos(new_part_var); - new_part_momentum = TLorentzVector mom(new_px, new_py, new_pz, pE); + TLorentzVector mom(new_px, new_py, new_pz, pE); + new_part_momentum = mom; } else if (modVar == "phi") { double p = new_part_momentum.P(); double theta = new_part_momentum.Theta(); double new_px = p * sin(theta) * cos(new_part_var); double new_py = p * sin(theta) * sin(new_part_var); double new_pz = p * cos(theta); - new_part_momentum = TLorentzVector mom(new_px, new_py, new_pz, pE); + TLorentzVector mom(new_px, new_py, new_pz, pE); + new_part_momentum = mom; } else if (modVar == "p") { double theta = new_part_momentum.Theta(); double phi = new_part_momentum.Phi(); double new_px = new_part_var * sin(theta) * cos(phi); double new_py = new_part_var * sin(theta) * sin(phi); double new_pz = new_part_var * cos(theta); - new_part_momentum = TLorentzVector mom(new_px, new_py, new_pz, pE); + TLorentzVector mom(new_px, new_py, new_pz, pE); + new_part_momentum = mom; } else if ( modVar == "x") { new_part_position = TLorentzVector(new_part_var, y, z, T); } else if ( modVar == "y") { @@ -190,9 +146,11 @@ namespace evgen { } // In-place convenience wrapper - inline void ModifyParticle(simb::MCTruth& mcTruth) { + inline void ModifyParticle(int modPdg, std::string modVar, + const std::vector& varBinEdges, const std::vector& varBinProbs, + int seed, simb::MCTruth& mcTruth) { simb::MCTruth newTruth; - ModifyParticle(mcTruth, newTruth); + ModifyParticle(modPdg, modVar, varBinEdges, varBinProbs, seed, mcTruth, newTruth); mcTruth = newTruth; } From 53a195bceb1ed8dc47b6cb75637d661874bd4f84 Mon Sep 17 00:00:00 2001 From: eyandel Date: Tue, 17 Mar 2026 12:48:03 -0500 Subject: [PATCH 07/19] mass should be in GeV --- ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl index 402c437a8..91eb58cea 100644 --- a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl +++ b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl @@ -20,7 +20,7 @@ physics.producers.generator.GenerateIsotropicSinglePhoton: false physics.producers.generator.ModifyParticle: true physics.producers.generator.ModifyParticlePdg: 111 # pi0 physics.producers.generator.ModifyParticleVar: ["mass"] # modify mass of pi0s -physics.producers.generator.ModifyParticleVarBinEdges: [[0.0, 1000.0]] # bin edges for mass modification (in GeV) +physics.producers.generator.ModifyParticleVarBinEdges: [[0.0, 1.000]] # bin edges for mass modification (in GeV) physics.producers.generator.ModifyParticleVarBinProbs: [[1.00]] # probabilities for each mass bin physics.producers.generator.ModifyParticleRandomSeed: 5 # random seed for particle modification, default is 0 (no seed) From a33a33daae1853233f922e700c43ea878f0623fa Mon Sep 17 00:00:00 2001 From: eyandel Date: Tue, 24 Mar 2026 14:40:28 -0500 Subject: [PATCH 08/19] manualy decay pi0s when changing mass --- ubsim/NGEM/GENIEGenNGEM_module.cc | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/ubsim/NGEM/GENIEGenNGEM_module.cc b/ubsim/NGEM/GENIEGenNGEM_module.cc index b2834ecb4..f4f3bdfd1 100644 --- a/ubsim/NGEM/GENIEGenNGEM_module.cc +++ b/ubsim/NGEM/GENIEGenNGEM_module.cc @@ -501,6 +501,10 @@ namespace evgen{ for (uint i_v = 0; i_v < fModifyParticleVar.size(); ++i_v) { std::cout << " Modifying variable " << fModifyParticleVar.at(i_v) << std::endl; ModifyParticle(fModifyParticlePdg, fModifyParticleVar.at(i_v), fModifyParticleVarBinEdges.at(i_v), fModifyParticleVarBinProbs.at(i_v), fModifyParticleRandomSeed, truth); + if (fModifyParticleVar.at(i_v) == "mass" && fModifyParticlePdg == 111) { + std::cout << "Manually decaying pi0s" << std::endl; + ManuallyDecayPi0sToTwoPhotons(truth); + } } std::cout << "truth.NParticles(): " << truth.NParticles() << std::endl; for (int i = 0; i < truth.NParticles(); ++i) { From 8b98a4446d568f70549c4c1d4402a4982e1b8ed1 Mon Sep 17 00:00:00 2001 From: eyandel Date: Tue, 24 Mar 2026 16:55:48 -0500 Subject: [PATCH 09/19] fix new mass energy and momentum --- ubsim/NGEM/ModifyParticle.h | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/ubsim/NGEM/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h index e3d03d44c..cd016ce52 100644 --- a/ubsim/NGEM/ModifyParticle.h +++ b/ubsim/NGEM/ModifyParticle.h @@ -68,6 +68,14 @@ namespace evgen { if (modVar == "mass"){ new_part_mass = new_part_var; + double new_pE_rest = new_part_mass; + double new_px_rest = 0; //new_pE_rest * sin(phi) * cos(theta); + double new_py_rest = 0; //new_pE_rest * sin(phi) * sin(theta); + double new_pz_rest = 0; //new_pE_rest * cos(phi); + TLorentzVector new_rest_frame_momentum(new_px_rest, new_py_rest, new_pz_rest, new_pE_rest); + TLorentzVector new_lab_frame_momentum = new_rest_frame_momentum; + new_lab_frame_momentum.Boost(part_boost_vector); + new_part_momentum = new_lab_frame_momentum; } else if (modVar == "E") { TLorentzVector mom(px, py, pz, new_part_var); new_part_momentum = mom; From b396616c513e0b07f733c0e27f7016062a2ec971 Mon Sep 17 00:00:00 2001 From: eyandel Date: Tue, 24 Mar 2026 17:34:20 -0500 Subject: [PATCH 10/19] conserve energy when modifying mass --- ubsim/NGEM/ModifyParticle.h | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/ubsim/NGEM/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h index cd016ce52..129d0b512 100644 --- a/ubsim/NGEM/ModifyParticle.h +++ b/ubsim/NGEM/ModifyParticle.h @@ -68,14 +68,20 @@ namespace evgen { if (modVar == "mass"){ new_part_mass = new_part_var; - double new_pE_rest = new_part_mass; - double new_px_rest = 0; //new_pE_rest * sin(phi) * cos(theta); - double new_py_rest = 0; //new_pE_rest * sin(phi) * sin(theta); - double new_pz_rest = 0; //new_pE_rest * cos(phi); - TLorentzVector new_rest_frame_momentum(new_px_rest, new_py_rest, new_pz_rest, new_pE_rest); - TLorentzVector new_lab_frame_momentum = new_rest_frame_momentum; - new_lab_frame_momentum.Boost(part_boost_vector); - new_part_momentum = new_lab_frame_momentum; + double new_pE = pE; + // Compute new |p| from E^2 = p^2 + m^2 + double new_p2 = new_pE*new_pE - new_part_mass*new_part_mass; + if (new_p2 < 0) new_p2 = 0; // protect against numerical issues / invalid input + double new_p = std::sqrt(new_p2); + + // Keep original direction + TVector3 dir = part_momentum.Vect().Unit(); + + double new_px = new_p * dir.X(); + double new_py = new_p * dir.Y(); + double new_pz = new_p * dir.Z(); + TLorentzVector mom(new_px, new_py, new_pz, new_pE); + new_part_momentum = mom; } else if (modVar == "E") { TLorentzVector mom(px, py, pz, new_part_var); new_part_momentum = mom; From a43b7fe9e9dad77a8177aa2a91295defd076c5b6 Mon Sep 17 00:00:00 2001 From: eyandel Date: Tue, 7 Apr 2026 17:38:10 -0500 Subject: [PATCH 11/19] try to fix energy and momentum conservation for all options, remove px py and pz mods --- ubsim/NGEM/ModifyParticle.h | 41 ++++++++++++++++++++++--------------- 1 file changed, 25 insertions(+), 16 deletions(-) diff --git a/ubsim/NGEM/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h index 129d0b512..23a5f29e6 100644 --- a/ubsim/NGEM/ModifyParticle.h +++ b/ubsim/NGEM/ModifyParticle.h @@ -83,16 +83,19 @@ namespace evgen { TLorentzVector mom(new_px, new_py, new_pz, new_pE); new_part_momentum = mom; } else if (modVar == "E") { - TLorentzVector mom(px, py, pz, new_part_var); - new_part_momentum = mom; - } else if (modVar == "px") { - TLorentzVector mom(new_part_var, py, pz, pE); - new_part_momentum = mom; - } else if (modVar == "py") { - TLorentzVector mom(px, new_part_var, pz, pE); - new_part_momentum = mom; - } else if (modVar == "pz") { - TLorentzVector mom(px, py, new_part_var, pE); + double new_pE = new_part_var; + // Compute new |p| from E^2 = p^2 + m^2 + double new_p2 = new_pE*new_pE - part_mass*part_mass; + if (new_p2 < 0) new_p2 = 0; // protect against numerical issues / invalid input + double new_p = std::sqrt(new_p2); + + // Keep original direction + TVector3 dir = part_momentum.Vect().Unit(); + + double new_px = new_p * dir.X(); + double new_py = new_p * dir.Y(); + double new_pz = new_p * dir.Z(); + TLorentzVector mom(new_px, new_py, new_pz, new_pE); new_part_momentum = mom; } else if (modVar == "theta") { double p = new_part_momentum.P(); @@ -111,12 +114,18 @@ namespace evgen { TLorentzVector mom(new_px, new_py, new_pz, pE); new_part_momentum = mom; } else if (modVar == "p") { - double theta = new_part_momentum.Theta(); - double phi = new_part_momentum.Phi(); - double new_px = new_part_var * sin(theta) * cos(phi); - double new_py = new_part_var * sin(theta) * sin(phi); - double new_pz = new_part_var * cos(theta); - TLorentzVector mom(new_px, new_py, new_pz, pE); + double new_p = new_part_var; + // Compute new |p| from E^2 = p^2 + m^2 + double new_pE = std::sqrt(new_p*new_p + part_mass*part_mass); + if (new_pE < 0) new_pE = 0; // protect against numerical issues / invalid input + + // Keep original direction + TVector3 dir = part_momentum.Vect().Unit(); + + double new_px = new_p * dir.X(); + double new_py = new_p * dir.Y(); + double new_pz = new_p * dir.Z(); + TLorentzVector mom(new_px, new_py, new_pz, new_pE); new_part_momentum = mom; } else if ( modVar == "x") { new_part_position = TLorentzVector(new_part_var, y, z, T); From e8043ae3d1f5aa6cf70b43b1b3fbd0617a5613d5 Mon Sep 17 00:00:00 2001 From: eyandel Date: Tue, 7 Apr 2026 17:43:07 -0500 Subject: [PATCH 12/19] fixed comments --- ubsim/NGEM/GENIEGenNGEM_module.cc | 2 +- ubsim/NGEM/ModifyParticle.h | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/ubsim/NGEM/GENIEGenNGEM_module.cc b/ubsim/NGEM/GENIEGenNGEM_module.cc index f4f3bdfd1..619fc1167 100644 --- a/ubsim/NGEM/GENIEGenNGEM_module.cc +++ b/ubsim/NGEM/GENIEGenNGEM_module.cc @@ -175,7 +175,7 @@ namespace evgen { bool fModifyParticle; ///< whether to modify particles according to configuration parameters int fModifyParticlePdg; ///< pdg code of particle to modify std::vector fModifyParticleVar; ///< variables to modify - //(options: "E", "px", "py", "pz", "mass", "x", "y", "z", "T", "theta", "phi", "p"); configurable via FHiCL + //(options: "E", "p", "mass", "x", "y", "z", "T", "theta", "phi"); configurable via FHiCL std::vector> fModifyParticleVarBinEdges; ///< configurable via FHiCL std::vector> fModifyParticleVarBinProbs; ///< configurable via FHiCL int fModifyParticleRandomSeed; ///< random seed for particle modification; configurable via FHiCL diff --git a/ubsim/NGEM/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h index 23a5f29e6..e7fd45b00 100644 --- a/ubsim/NGEM/ModifyParticle.h +++ b/ubsim/NGEM/ModifyParticle.h @@ -45,10 +45,10 @@ namespace evgen { double part_mass = part.Mass(); //double theta = randomGen.Uniform(0, 2 * M_PI); //double phi = acos(1 - 2 * randomGen.Uniform(0, 1)); - double pE = part_momentum.E(); //part_mass / 2; - double px = part_momentum.Px(); //pE * sin(phi) * cos(theta); - double py = part_momentum.Py(); //pE * sin(phi) * sin(theta); - double pz = part_momentum.Pz(); //pE * cos(phi); + double pE = part_momentum.E(); + //double px = part_momentum.Px(); + //double py = part_momentum.Py(); + //double pz = part_momentum.Pz(); //TLorentzVector rest_frame_momentum(px, py, pz, pE); double T = part.T(); From f66bd3ad7e21d4406b25e920de5d75388812cdb7 Mon Sep 17 00:00:00 2001 From: eyandel Date: Wed, 8 Apr 2026 12:18:46 -0500 Subject: [PATCH 13/19] make decaying the modified pi0 a fcl parameter --- ubsim/NGEM/GENIEGenNGEM_module.cc | 4 +++- ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl | 5 +++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/ubsim/NGEM/GENIEGenNGEM_module.cc b/ubsim/NGEM/GENIEGenNGEM_module.cc index 619fc1167..92b8701d8 100644 --- a/ubsim/NGEM/GENIEGenNGEM_module.cc +++ b/ubsim/NGEM/GENIEGenNGEM_module.cc @@ -179,6 +179,7 @@ namespace evgen { std::vector> fModifyParticleVarBinEdges; ///< configurable via FHiCL std::vector> fModifyParticleVarBinProbs; ///< configurable via FHiCL int fModifyParticleRandomSeed; ///< random seed for particle modification; configurable via FHiCL + bool fManuallyDecayModifiedPi0s; ///< whether to manually decay pi0s after modification. must use if mass is modified. }; } @@ -207,6 +208,7 @@ namespace evgen{ , fModifyParticleVarBinEdges(pset.get>>("ModifyParticleVarBinEdges", std::vector>())) , fModifyParticleVarBinProbs(pset.get>>("ModifyParticleVarBinProbs", std::vector>())) , fModifyParticleRandomSeed(pset.get("ModifyParticleRandomSeed", 0)) + , fManuallyDecayModifiedPi0s(pset.get("ManuallyDecayModifiedPi0s", false)) { fStopwatch.Start(); @@ -501,7 +503,7 @@ namespace evgen{ for (uint i_v = 0; i_v < fModifyParticleVar.size(); ++i_v) { std::cout << " Modifying variable " << fModifyParticleVar.at(i_v) << std::endl; ModifyParticle(fModifyParticlePdg, fModifyParticleVar.at(i_v), fModifyParticleVarBinEdges.at(i_v), fModifyParticleVarBinProbs.at(i_v), fModifyParticleRandomSeed, truth); - if (fModifyParticleVar.at(i_v) == "mass" && fModifyParticlePdg == 111) { + if (fManuallyDecayModifiedPi0s){ //fModifyParticleVar.at(i_v) == "mass" && fModifyParticlePdg == 111) { std::cout << "Manually decaying pi0s" << std::endl; ManuallyDecayPi0sToTwoPhotons(truth); } diff --git a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl index 91eb58cea..a559fb45d 100644 --- a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl +++ b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl @@ -13,13 +13,14 @@ physics.producers.generator.TopVolume: "volCryostat" physics.producers.generator.BeamName: "booster" physics.producers.generator.module_type: "GENIEGenNGEM" -physics.producers.generator.ManuallyDecayPi0s: false -physics.producers.generator.DeleteRandomGamma: false +physics.producers.generator.ManuallyDecayPi0s: false # manually decay pi0s (before modification) +physics.producers.generator.DeleteRandomGamma: false # whether to delete one of the two photons from pi0 decay. happens before particle modification. must have ManuallyDecayPi0s is true. physics.producers.generator.GenerateIsotropicSinglePhoton: false physics.producers.generator.ModifyParticle: true physics.producers.generator.ModifyParticlePdg: 111 # pi0 physics.producers.generator.ModifyParticleVar: ["mass"] # modify mass of pi0s +physics.producers.generator.ManuallyDecayModifiedPi0s: true # whether to manually decay pi0s after modification. must use if mass is modified. physics.producers.generator.ModifyParticleVarBinEdges: [[0.0, 1.000]] # bin edges for mass modification (in GeV) physics.producers.generator.ModifyParticleVarBinProbs: [[1.00]] # probabilities for each mass bin physics.producers.generator.ModifyParticleRandomSeed: 5 # random seed for particle modification, default is 0 (no seed) From 77fa1ad700bc16a4b58f3eaa1ab19133df2c13ca Mon Sep 17 00:00:00 2001 From: eyandel Date: Wed, 8 Apr 2026 12:35:05 -0500 Subject: [PATCH 14/19] fix some checks for invalid negatve inputs and add ability to delete one modified gamma --- ubsim/NGEM/GENIEGenNGEM_module.cc | 8 +++++++- ubsim/NGEM/ModifyParticle.h | 17 ++++++++++++++--- ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl | 3 ++- 3 files changed, 23 insertions(+), 5 deletions(-) diff --git a/ubsim/NGEM/GENIEGenNGEM_module.cc b/ubsim/NGEM/GENIEGenNGEM_module.cc index 92b8701d8..109971740 100644 --- a/ubsim/NGEM/GENIEGenNGEM_module.cc +++ b/ubsim/NGEM/GENIEGenNGEM_module.cc @@ -180,6 +180,7 @@ namespace evgen { std::vector> fModifyParticleVarBinProbs; ///< configurable via FHiCL int fModifyParticleRandomSeed; ///< random seed for particle modification; configurable via FHiCL bool fManuallyDecayModifiedPi0s; ///< whether to manually decay pi0s after modification. must use if mass is modified. + bool fDeleteRandomModifiedGamma; ///< whether to delete a random gamma after modification. }; } @@ -209,6 +210,7 @@ namespace evgen{ , fModifyParticleVarBinProbs(pset.get>>("ModifyParticleVarBinProbs", std::vector>())) , fModifyParticleRandomSeed(pset.get("ModifyParticleRandomSeed", 0)) , fManuallyDecayModifiedPi0s(pset.get("ManuallyDecayModifiedPi0s", false)) + , fDeleteRandomModifiedGamma(pset.get("DeleteRandomModifiedGamma", false)) { fStopwatch.Start(); @@ -504,9 +506,13 @@ namespace evgen{ std::cout << " Modifying variable " << fModifyParticleVar.at(i_v) << std::endl; ModifyParticle(fModifyParticlePdg, fModifyParticleVar.at(i_v), fModifyParticleVarBinEdges.at(i_v), fModifyParticleVarBinProbs.at(i_v), fModifyParticleRandomSeed, truth); if (fManuallyDecayModifiedPi0s){ //fModifyParticleVar.at(i_v) == "mass" && fModifyParticlePdg == 111) { - std::cout << "Manually decaying pi0s" << std::endl; + std::cout << "Manually decaying modified pi0s" << std::endl; ManuallyDecayPi0sToTwoPhotons(truth); } + if (fDeleteRandomModifiedGamma) { + std::cout << "Deleting random modified gamma" << std::endl; + DeleteOneRandomPhoton(truth); + } } std::cout << "truth.NParticles(): " << truth.NParticles() << std::endl; for (int i = 0; i < truth.NParticles(); ++i) { diff --git a/ubsim/NGEM/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h index e7fd45b00..b7e41d190 100644 --- a/ubsim/NGEM/ModifyParticle.h +++ b/ubsim/NGEM/ModifyParticle.h @@ -71,7 +71,8 @@ namespace evgen { double new_pE = pE; // Compute new |p| from E^2 = p^2 + m^2 double new_p2 = new_pE*new_pE - new_part_mass*new_part_mass; - if (new_p2 < 0) new_p2 = 0; // protect against numerical issues / invalid input + //if the new particle mass is greater than the original particle total energy, create the new particle at rest rather than conserving energy + if (new_p2 < 0) new_p2 = 0; double new_p = std::sqrt(new_p2); // Keep original direction @@ -84,9 +85,15 @@ namespace evgen { new_part_momentum = mom; } else if (modVar == "E") { double new_pE = new_part_var; + // protect against numerical issues / invalid input + if (new_pE < 0) { + cout << "Warning: Invalid particle parameters (negative energy or momentum)" << endl; + continue; + } // Compute new |p| from E^2 = p^2 + m^2 double new_p2 = new_pE*new_pE - part_mass*part_mass; - if (new_p2 < 0) new_p2 = 0; // protect against numerical issues / invalid input + // if the new energy is set smaller than the rest mass, we will create the particle at rest + if (new_pE < part_mass) new_pE = part_mass; double new_p = std::sqrt(new_p2); // Keep original direction @@ -115,9 +122,13 @@ namespace evgen { new_part_momentum = mom; } else if (modVar == "p") { double new_p = new_part_var; + // protect against numerical issues / invalid input + if (new_p < 0) { + cout << "Warning: Invalid particle parameters (negative energy or momentum)" << endl; + continue; + } // Compute new |p| from E^2 = p^2 + m^2 double new_pE = std::sqrt(new_p*new_p + part_mass*part_mass); - if (new_pE < 0) new_pE = 0; // protect against numerical issues / invalid input // Keep original direction TVector3 dir = part_momentum.Vect().Unit(); diff --git a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl index a559fb45d..0d8a45a32 100644 --- a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl +++ b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl @@ -14,13 +14,14 @@ physics.producers.generator.BeamName: "booster" physics.producers.generator.module_type: "GENIEGenNGEM" physics.producers.generator.ManuallyDecayPi0s: false # manually decay pi0s (before modification) -physics.producers.generator.DeleteRandomGamma: false # whether to delete one of the two photons from pi0 decay. happens before particle modification. must have ManuallyDecayPi0s is true. +physics.producers.generator.DeleteRandomGamma: false # whether to delete one of the two photons from pi0 decay. happens before particle modification. must have ManuallyDecayPi0s true. physics.producers.generator.GenerateIsotropicSinglePhoton: false physics.producers.generator.ModifyParticle: true physics.producers.generator.ModifyParticlePdg: 111 # pi0 physics.producers.generator.ModifyParticleVar: ["mass"] # modify mass of pi0s physics.producers.generator.ManuallyDecayModifiedPi0s: true # whether to manually decay pi0s after modification. must use if mass is modified. +physics.producers.generator.DeleteRandomModifiedGamma: false # whether to delete a random gamma after modification. must have ManuallyDecayModifiedPi0s true. physics.producers.generator.ModifyParticleVarBinEdges: [[0.0, 1.000]] # bin edges for mass modification (in GeV) physics.producers.generator.ModifyParticleVarBinProbs: [[1.00]] # probabilities for each mass bin physics.producers.generator.ModifyParticleRandomSeed: 5 # random seed for particle modification, default is 0 (no seed) From 9b77c24d1b7c116c8429b087230e29a6299169f3 Mon Sep 17 00:00:00 2001 From: eyandel Date: Wed, 8 Apr 2026 12:36:49 -0500 Subject: [PATCH 15/19] fix missing std:: on cout --- ubsim/NGEM/ModifyParticle.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ubsim/NGEM/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h index b7e41d190..98cdac535 100644 --- a/ubsim/NGEM/ModifyParticle.h +++ b/ubsim/NGEM/ModifyParticle.h @@ -87,7 +87,7 @@ namespace evgen { double new_pE = new_part_var; // protect against numerical issues / invalid input if (new_pE < 0) { - cout << "Warning: Invalid particle parameters (negative energy or momentum)" << endl; + std::cout << "Warning: Invalid particle parameters (negative energy or momentum)" << std::endl; continue; } // Compute new |p| from E^2 = p^2 + m^2 @@ -124,7 +124,7 @@ namespace evgen { double new_p = new_part_var; // protect against numerical issues / invalid input if (new_p < 0) { - cout << "Warning: Invalid particle parameters (negative energy or momentum)" << endl; + std::cout << "Warning: Invalid particle parameters (negative energy or momentum)" << std::endl; continue; } // Compute new |p| from E^2 = p^2 + m^2 From aa40c226d4c78c3522580a5db018b5e62a712995 Mon Sep 17 00:00:00 2001 From: Herbert Greenlee Date: Tue, 28 Apr 2026 12:57:29 -0500 Subject: [PATCH 16/19] Fix one definition rule violation. --- .../Calculators/G4RWManagerService.cc | 18 ++++++++++++++++++ .../Calculators/G4RWManagerService_service.cc | 19 +------------------ 2 files changed, 19 insertions(+), 18 deletions(-) create mode 100644 ubsim/EventWeight/Calculators/G4RWManagerService.cc diff --git a/ubsim/EventWeight/Calculators/G4RWManagerService.cc b/ubsim/EventWeight/Calculators/G4RWManagerService.cc new file mode 100644 index 000000000..bb0c1c45c --- /dev/null +++ b/ubsim/EventWeight/Calculators/G4RWManagerService.cc @@ -0,0 +1,18 @@ +#include "G4RWManagerService.h" + +evwgh::G4RWManagerService::G4RWManagerService( + fhicl::ParameterSet const& pset) : + fMaterial(pset.get("Material")), + fRWManager(new G4ReweightManager({fMaterial})) {} + +evwgh::G4RWManagerService::G4RWManagerService( + fhicl::ParameterSet const& pset, + art::ActivityRegistry&) : G4RWManagerService(pset) {} + +G4ReweightManager * evwgh::G4RWManagerService::GetManager() const { + return fRWManager; +} + +const fhicl::ParameterSet & evwgh::G4RWManagerService::GetMaterial() const { + return fMaterial; +} diff --git a/ubsim/EventWeight/Calculators/G4RWManagerService_service.cc b/ubsim/EventWeight/Calculators/G4RWManagerService_service.cc index d0da237ab..321566e54 100644 --- a/ubsim/EventWeight/Calculators/G4RWManagerService_service.cc +++ b/ubsim/EventWeight/Calculators/G4RWManagerService_service.cc @@ -1,20 +1,3 @@ #include "G4RWManagerService.h" -evwgh::G4RWManagerService::G4RWManagerService( - fhicl::ParameterSet const& pset) : - fMaterial(pset.get("Material")), - fRWManager(new G4ReweightManager({fMaterial})) {} - -evwgh::G4RWManagerService::G4RWManagerService( - fhicl::ParameterSet const& pset, - art::ActivityRegistry&) : G4RWManagerService(pset) {} - -G4ReweightManager * evwgh::G4RWManagerService::GetManager() const { - return fRWManager; -} - -const fhicl::ParameterSet & evwgh::G4RWManagerService::GetMaterial() const { - return fMaterial; -} - -DEFINE_ART_SERVICE(evwgh::G4RWManagerService) \ No newline at end of file +DEFINE_ART_SERVICE(evwgh::G4RWManagerService) From 225dab1f91b421fa000b696e49bf2579b3e9e6b8 Mon Sep 17 00:00:00 2001 From: eyandel Date: Wed, 29 Apr 2026 18:36:16 -0500 Subject: [PATCH 17/19] fix example fcl to not set seed --- ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl index 0d8a45a32..1aa8c404d 100644 --- a/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl +++ b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl @@ -24,7 +24,7 @@ physics.producers.generator.ManuallyDecayModifiedPi0s: true # whether to manuall physics.producers.generator.DeleteRandomModifiedGamma: false # whether to delete a random gamma after modification. must have ManuallyDecayModifiedPi0s true. physics.producers.generator.ModifyParticleVarBinEdges: [[0.0, 1.000]] # bin edges for mass modification (in GeV) physics.producers.generator.ModifyParticleVarBinProbs: [[1.00]] # probabilities for each mass bin -physics.producers.generator.ModifyParticleRandomSeed: 5 # random seed for particle modification, default is 0 (no seed) +#physics.producers.generator.ModifyParticleRandomSeed: 5 # random seed for particle modification, default is 0 (no seed) physics.filters.finalstatefilter.IsVerbose: false physics.filters.finalstatefilter.isInclusive: true From bce09d68a35ae7974b73f22be78c3ec99d9099e0 Mon Sep 17 00:00:00 2001 From: Herbert Greenlee Date: Fri, 1 May 2026 09:54:11 -0500 Subject: [PATCH 18/19] v08_00_00_96 --- CMakeLists.txt | 106 -------------------------------------- ubsim/CMakeLists.txt | 2 +- ubsim/NGEM/CMakeLists.txt | 106 ++++++++++++++++++++++++++++++++++++++ ups/product_deps | 2 +- 4 files changed, 108 insertions(+), 108 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d866c7274..a9720e9f6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -83,112 +83,6 @@ include(ArtDictionary) include(ArtMake) include(BuildPlugins) -# Determine the major version number of GENIE based on the GENIE_VERSION -# variable set by the genie ups product. Define the GENIE_PRE_R3 preprocessor -# macro (needed to set the right GENIE header locations, etc.) if needed. -# -# The top-level nutools CMakeLists.txt does this as well, but I don't think -# that we can access the result of that version check in the larsim build. -parse_ups_version( ${GENIE_VERSION} ) -include_directories ( $ENV{GENIE_INC}/GENIE ) -if( ${VMAJ} LESS 3 ) - - add_definitions(-DGENIE_PRE_R3) - - # GENIE v2 library list - set(GENIE_LIB_LIST ${GBARYONRESONANCE} - ${GBASE} - ${GBODEKYANG} - ${GCHARM} - ${GCOH} - ${GDFRC} - ${GDIS} - ${GCROSSSECTIONS} - ${GDECAY} - ${GELAS} - ${GELFF} - ${GEVGCORE} - ${GEVGMODULES} - ${GGIBUU} - ${GHADRONTRANSP} - ${GFRAGMENTATION} - ${GLLEWELLYNSMITH} - ${GMEC} - ${GNUGAMMA} - ${GNUE} - ${GNUCLEAR} - ${GNUMERICAL} - ${GQPM} - ${GPDG} - ${GPDF} - ${GQEL} - ${GRES} - ${GREGISTRY} - ${GREINSEHGAL} - ${GGEO} - ${GMUELOSS} - ${GREWEIGHT} - ${GNUCLEONDECAY} ) - - #MESSAGE("--LARSIM-- GENIE_VERSION=${VMAJ} ${VMIN} -DGENIE_PRE_R3") -else() - - # GENIE v3 library list - set (GENIE_LIB_LIST ${GFWMSG} - ${GFWREG} - ${GFWALG} - ${GFWINT} - ${GFWGHEP} - ${GFWNUM} - ${GSL} # FWNUM relies on GSL - ${GFWUTL} - ${GFWPARDAT} - ${GFWEG} - ${GFWNTP} - ${GPHXSIG} - ${GPHPDF} - ${GPHNUCLST} - ${GPHCMN} - ${GPHDCY} - ${GPHHADTRANSP} - ${GPHHADNZ} - ${GPHDEEX} - ${GPHAMNGXS} - ${GPHAMNGEG} - ${GPHCHMXS} - ${GPHCOHXS} - ${GPHCOHEG} - ${GPHDISXS} - ${GPHDISEG} - ${GPHDFRCXS} - ${GPHDFRCEG} - ${GPHGLWRESXS} - ${GPHGLWRESEG} - ${GPHIBDXS} - ${GPHIBDEG} - ${GPHMNUCXS} - ${GPHMNUCEG} - ${GPHMEL} - ${GPHNUELXS} - ${GPHNUELEG} - ${GPHQELXS} - ${GPHQELEG} - ${GPHRESXS} - ${GPHRESEG} - ${GPHSTRXS} - ${GPHSTREG} - ${GPHNDCY} - ${GTLGEO} - ${GTLFLX} - ${GRWFWK} - ${GRWIO} - ${GRWCLC}) - - #MESSAGE("--LARSIM-- GENIE_VERSION=${VMAJ} ${VMIN} not defining -DGENIE_PRE_R3") -endif() - -#MESSAGE("--LARSIM-- GENIE ${VMAJ} GENIE_LIB_LIST=${GENIE_LIB_LIST}") - # ADD SOURCE CODE SUBDIRECTORIES HERE add_subdirectory(ubsim) diff --git a/ubsim/CMakeLists.txt b/ubsim/CMakeLists.txt index 24df4cdc6..98de4a685 100644 --- a/ubsim/CMakeLists.txt +++ b/ubsim/CMakeLists.txt @@ -6,4 +6,4 @@ add_subdirectory(PhotonPropagation) add_subdirectory(SNStreamSim) add_subdirectory(Simulation) add_subdirectory(Filters) -add_subdirectory(NGEM) \ No newline at end of file +add_subdirectory(NGEM) diff --git a/ubsim/NGEM/CMakeLists.txt b/ubsim/NGEM/CMakeLists.txt index fa5527f53..608854be8 100644 --- a/ubsim/NGEM/CMakeLists.txt +++ b/ubsim/NGEM/CMakeLists.txt @@ -9,6 +9,112 @@ include_directories ( $ENV{DK2NUGENIE_INC} ) cet_find_library( DK2NUDATA NAMES dk2nuTree PATHS ENV DK2NUDATA_LIB ) cet_find_library( DK2NUGENIE NAMES dk2nuGenie PATHS ENV DK2NUGENIE_LIB ) +# Determine the major version number of GENIE based on the GENIE_VERSION +# variable set by the genie ups product. Define the GENIE_PRE_R3 preprocessor +# macro (needed to set the right GENIE header locations, etc.) if needed. +# +# The top-level nutools CMakeLists.txt does this as well, but I don't think +# that we can access the result of that version check in the larsim build. +parse_ups_version( ${GENIE_VERSION} ) +include_directories ( $ENV{GENIE_INC}/GENIE ) +if( ${VMAJ} LESS 3 ) + + add_definitions(-DGENIE_PRE_R3) + + # GENIE v2 library list + set(GENIE_LIB_LIST ${GBARYONRESONANCE} + ${GBASE} + ${GBODEKYANG} + ${GCHARM} + ${GCOH} + ${GDFRC} + ${GDIS} + ${GCROSSSECTIONS} + ${GDECAY} + ${GELAS} + ${GELFF} + ${GEVGCORE} + ${GEVGMODULES} + ${GGIBUU} + ${GHADRONTRANSP} + ${GFRAGMENTATION} + ${GLLEWELLYNSMITH} + ${GMEC} + ${GNUGAMMA} + ${GNUE} + ${GNUCLEAR} + ${GNUMERICAL} + ${GQPM} + ${GPDG} + ${GPDF} + ${GQEL} + ${GRES} + ${GREGISTRY} + ${GREINSEHGAL} + ${GGEO} + ${GMUELOSS} + ${GREWEIGHT} + ${GNUCLEONDECAY} ) + + #MESSAGE("--LARSIM-- GENIE_VERSION=${VMAJ} ${VMIN} -DGENIE_PRE_R3") +else() + + # GENIE v3 library list + set (GENIE_LIB_LIST ${GFWMSG} + ${GFWREG} + ${GFWALG} + ${GFWINT} + ${GFWGHEP} + ${GFWNUM} + ${GSL} # FWNUM relies on GSL + ${GFWUTL} + ${GFWPARDAT} + ${GFWEG} + ${GFWNTP} + ${GPHXSIG} + ${GPHPDF} + ${GPHNUCLST} + ${GPHCMN} + ${GPHDCY} + ${GPHHADTRANSP} + ${GPHHADNZ} + ${GPHDEEX} + ${GPHAMNGXS} + ${GPHAMNGEG} + ${GPHCHMXS} + ${GPHCOHXS} + ${GPHCOHEG} + ${GPHDISXS} + ${GPHDISEG} + ${GPHDFRCXS} + ${GPHDFRCEG} + ${GPHGLWRESXS} + ${GPHGLWRESEG} + ${GPHIBDXS} + ${GPHIBDEG} + ${GPHMNUCXS} + ${GPHMNUCEG} + ${GPHMEL} + ${GPHNUELXS} + ${GPHNUELEG} + ${GPHQELXS} + ${GPHQELEG} + ${GPHRESXS} + ${GPHRESEG} + ${GPHSTRXS} + ${GPHSTREG} + ${GPHNDCY} + ${GTLGEO} + ${GTLFLX} + ${GRWFWK} + ${GRWIO} + ${GRWCLC}) + + #MESSAGE("--LARSIM-- GENIE_VERSION=${VMAJ} ${VMIN} not defining -DGENIE_PRE_R3") +endif() + +#MESSAGE("--LARSIM-- GENIE ${VMAJ} GENIE_LIB_LIST=${GENIE_LIB_LIST}") + art_make( MODULE_LIBRARIES lardataalg_MCDumpers diff --git a/ups/product_deps b/ups/product_deps index 4f9d39b40..40a7082a9 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -2,7 +2,7 @@ # The *parent* line must the first non-commented line and defines this product and version # The version must be of the form vxx_yy_zz (e.g. v01_02_03) -parent ubsim v08_00_00_94 +parent ubsim v08_00_00_96 defaultqual e17 From 5c9d860f8c3fa747702c4e946d5aee2fd73c10c0 Mon Sep 17 00:00:00 2001 From: Herbert Greenlee Date: Fri, 1 May 2026 16:23:01 -0500 Subject: [PATCH 19/19] MCC10 adaptations. --- CMakeLists.txt | 1 + ubsim/NGEM/CMakeLists.txt | 171 ++++-------------------------- ubsim/NGEM/GENIEGenNGEM_module.cc | 28 ++--- 3 files changed, 34 insertions(+), 166 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d31a18023..bdaddbc81 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,6 +38,7 @@ find_package( larcorealg REQUIRED EXPORT ) find_package( geant4reweight REQUIRED EXPORT ) find_package(LibXml2 REQUIRED EXPORT) find_package(log4cpp REQUIRED EXPORT) +find_package(GENIE REQUIRED EXPORT) # macros for dictionary and simple_plugin include(ArtDictionary) diff --git a/ubsim/NGEM/CMakeLists.txt b/ubsim/NGEM/CMakeLists.txt index 608854be8..a50374748 100644 --- a/ubsim/NGEM/CMakeLists.txt +++ b/ubsim/NGEM/CMakeLists.txt @@ -1,158 +1,23 @@ - cet_add_compiler_flags( CXX -DHIDE_GENIE_LOG_XXX ) -include_directories ( $ENV{LOG4CPP_INC} ) -include_directories ( $ENV{DK2NUDATA_INC} ) -include_directories ( $ENV{DK2NUGENIE_INC} ) - -# genie -cet_find_library( DK2NUDATA NAMES dk2nuTree PATHS ENV DK2NUDATA_LIB ) -cet_find_library( DK2NUGENIE NAMES dk2nuGenie PATHS ENV DK2NUGENIE_LIB ) - -# Determine the major version number of GENIE based on the GENIE_VERSION -# variable set by the genie ups product. Define the GENIE_PRE_R3 preprocessor -# macro (needed to set the right GENIE header locations, etc.) if needed. -# -# The top-level nutools CMakeLists.txt does this as well, but I don't think -# that we can access the result of that version check in the larsim build. -parse_ups_version( ${GENIE_VERSION} ) -include_directories ( $ENV{GENIE_INC}/GENIE ) -if( ${VMAJ} LESS 3 ) - - add_definitions(-DGENIE_PRE_R3) - - # GENIE v2 library list - set(GENIE_LIB_LIST ${GBARYONRESONANCE} - ${GBASE} - ${GBODEKYANG} - ${GCHARM} - ${GCOH} - ${GDFRC} - ${GDIS} - ${GCROSSSECTIONS} - ${GDECAY} - ${GELAS} - ${GELFF} - ${GEVGCORE} - ${GEVGMODULES} - ${GGIBUU} - ${GHADRONTRANSP} - ${GFRAGMENTATION} - ${GLLEWELLYNSMITH} - ${GMEC} - ${GNUGAMMA} - ${GNUE} - ${GNUCLEAR} - ${GNUMERICAL} - ${GQPM} - ${GPDG} - ${GPDF} - ${GQEL} - ${GRES} - ${GREGISTRY} - ${GREINSEHGAL} - ${GGEO} - ${GMUELOSS} - ${GREWEIGHT} - ${GNUCLEONDECAY} ) - - #MESSAGE("--LARSIM-- GENIE_VERSION=${VMAJ} ${VMIN} -DGENIE_PRE_R3") -else() - - # GENIE v3 library list - set (GENIE_LIB_LIST ${GFWMSG} - ${GFWREG} - ${GFWALG} - ${GFWINT} - ${GFWGHEP} - ${GFWNUM} - ${GSL} # FWNUM relies on GSL - ${GFWUTL} - ${GFWPARDAT} - ${GFWEG} - ${GFWNTP} - ${GPHXSIG} - ${GPHPDF} - ${GPHNUCLST} - ${GPHCMN} - ${GPHDCY} - ${GPHHADTRANSP} - ${GPHHADNZ} - ${GPHDEEX} - ${GPHAMNGXS} - ${GPHAMNGEG} - ${GPHCHMXS} - ${GPHCOHXS} - ${GPHCOHEG} - ${GPHDISXS} - ${GPHDISEG} - ${GPHDFRCXS} - ${GPHDFRCEG} - ${GPHGLWRESXS} - ${GPHGLWRESEG} - ${GPHIBDXS} - ${GPHIBDEG} - ${GPHMNUCXS} - ${GPHMNUCEG} - ${GPHMEL} - ${GPHNUELXS} - ${GPHNUELEG} - ${GPHQELXS} - ${GPHQELEG} - ${GPHRESXS} - ${GPHRESEG} - ${GPHSTRXS} - ${GPHSTREG} - ${GPHNDCY} - ${GTLGEO} - ${GTLFLX} - ${GRWFWK} - ${GRWIO} - ${GRWCLC}) - - #MESSAGE("--LARSIM-- GENIE_VERSION=${VMAJ} ${VMIN} not defining -DGENIE_PRE_R3") -endif() - -#MESSAGE("--LARSIM-- GENIE ${VMAJ} GENIE_LIB_LIST=${GENIE_LIB_LIST}") - -art_make( - MODULE_LIBRARIES - lardataalg_MCDumpers - larsim_Simulation lardataobj_Simulation - nusimdata_SimulationBase - larcoreobj_SummaryData - larcorealg_Geometry - larcore_Geometry_Geometry_service - nutools_RandomUtils_NuRandomService_service - nutools_EventGeneratorBase_GENIE - #nurandom_RandomUtils_NuRandomService_service - #nugen_EventGeneratorBase_GENIE - ${DK2NUDATA} - ${DK2NUGENIE} - #dk2nu_Tree - #dk2nu_Genie - ${ART_FRAMEWORK_CORE} - ${ART_FRAMEWORK_PRINCIPAL} - ${ART_FRAMEWORK_SERVICES_REGISTRY} - ${ART_FRAMEWORK_SERVICES_OPTIONAL} - ${ART_FRAMEWORK_SERVICES_OPTIONAL_TFILESERVICE_SERVICE} - ${ART_FRAMEWORK_SERVICES_OPTIONAL_RANDOMNUMBERGENERATOR_SERVICE} - art_Persistency_Common - art_Persistency_Provenance - art_Utilities - #art_root_io_TFileService_service - canvas - ${MF_MESSAGELOGGER} - ${FHICLCPP} - cetlib cetlib_except - ${CLHEP} - ${GENIE_LIB_LIST} - ${ROOT_EGPYTHIA6} # FIXME!!! - resolving genie run time reference - ${ROOT_BASIC_LIB_LIST} - ${ROOT_EG} - ${LOG4CPP} - #Generator_Framework_GHEP - ) +cet_build_plugin( + GENIEGenNGEM art::EDProducer + LIBRARIES + PRIVATE + lardata::AssociationUtil + lardataalg::MCDumpers + larcore::Geometry_Geometry_service + lardataobj::Simulation + nurandom::RandomUtils_NuRandomService_service + nugen::EventGeneratorBase + nugen::EventGeneratorBase_GENIE + art_root_io::TFileService_service + dk2nu::Genie + dk2nu::Tree + GENIE::GFwEG + ROOT::Tree + ROOT::EG +) install_headers() install_fhicl() diff --git a/ubsim/NGEM/GENIEGenNGEM_module.cc b/ubsim/NGEM/GENIEGenNGEM_module.cc index 109971740..b995e3ea3 100644 --- a/ubsim/NGEM/GENIEGenNGEM_module.cc +++ b/ubsim/NGEM/GENIEGenNGEM_module.cc @@ -33,14 +33,14 @@ #include "fhiclcpp/ParameterSet.h" #include "art/Framework/Principal/Handle.h" #include "art/Framework/Services/Registry/ServiceHandle.h" -#include "art/Framework/Services/Optional/TFileService.h" -#include "art/Framework/Services/Optional/TFileDirectory.h" +#include "art_root_io/TFileService.h" +#include "art_root_io/TFileDirectory.h" #include "messagefacility/MessageLogger/MessageLogger.h" #include "canvas/Persistency/Common/Assns.h" #include "art/Framework/Core/EDProducer.h" // art extensions -#include "nutools/RandomUtils/NuRandomService.h" +#include "nurandom/RandomUtils/NuRandomService.h" // LArSoft includes #include "lardataalg/MCDumpers/MCDumpers.h" // sim::dump namespace @@ -51,7 +51,7 @@ #include "larcore/Geometry/Geometry.h" #include "larcoreobj/SummaryData/RunData.h" #include "larcoreobj/SummaryData/POTSummary.h" -#include "nutools/EventGeneratorBase/GENIE/GENIEHelper.h" +#include "nugen/EventGeneratorBase/GENIE/GENIEHelper.h" #include "lardata/Utilities/AssociationUtil.h" // dk2nu extensions @@ -60,8 +60,8 @@ #include "dk2nu/genie/GDk2NuFlux.h" #include "GENIE/Framework/EventGen/EventRecord.h" -#include "nutools/EventGeneratorBase/GENIE/EVGBAssociationUtil.h" -#include "nutools/EventGeneratorBase/evgenbase.h" +#include "nugen/EventGeneratorBase/GENIE/EVGBAssociationUtil.h" +#include "nugen/EventGeneratorBase/evgenbase.h" // Optional Modifications, if you don't want to modify normal GENIE neutrino generation events #include "ManuallyDecayPi0sToTwoPhotons.h" @@ -189,7 +189,8 @@ namespace evgen{ //____________________________________________________________________________ GENIEGenNGEM::GENIEGenNGEM(fhicl::ParameterSet const& pset) - : fGENIEHelp(0) + : art::EDProducer(pset) + , fGENIEHelp(0) , fDefinedVtxHistRange (pset.get< bool >("DefinedVtxHistRange")) , fVtxPosHistRange (pset.get< std::vector >("VtxPosHistRange")) , fPassEmptySpills (pset.get< bool >("PassEmptySpills")) @@ -263,7 +264,7 @@ namespace evgen{ fGENIEHelp = new evgb::GENIEHelper(GENIEconfig, geo->ROOTGeoManager(), - geo->ROOTFile(), + geo->GDMLFile(), geo->TotalMass(pset.get< std::string>("TopVolume").c_str())); } @@ -325,9 +326,10 @@ namespace evgen{ fECons = tfs->make("fECons", ";#Delta E(#nu,lepton);", 500, -5., 5.); art::ServiceHandle geo; - double x = 2.1*geo->DetHalfWidth(); - double y = 2.1*geo->DetHalfHeight(); - double z = 2.*geo->DetLength(); + geo::TPCGeo const& tpc = geo->TPC(); + double x = 2.1*tpc.ActiveHalfWidth(); + double y = 2.1*tpc.ActiveHalfHeight(); + double z = 2.*tpc.ActiveLength(); int xdiv = TMath::Nint(2*x/5.); int ydiv = TMath::Nint(2*y/5.); int zdiv = TMath::Nint(2*z/5.); @@ -363,7 +365,7 @@ namespace evgen{ art::ServiceHandle geo; std::unique_ptr runcol(new sumdata::RunData(geo->DetectorName())); - run.put(std::move(runcol)); + run.put(std::move(runcol), art::fullRun()); return; } @@ -387,7 +389,7 @@ namespace evgen{ p->totpot = fGENIEHelp->TotalExposure() - fPrevTotPOT; p->totgoodpot = fGENIEHelp->TotalExposure() - fPrevTotGoodPOT; - sr.put(std::move(p)); + sr.put(std::move(p), art::subRunFragment()); return; }