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/CMakeLists.txt b/ubsim/CMakeLists.txt index 37003e9d3..36d23bb1c 100644 --- a/ubsim/CMakeLists.txt +++ b/ubsim/CMakeLists.txt @@ -7,3 +7,4 @@ add_subdirectory(SNStreamSim) add_subdirectory(Simulation) add_subdirectory(Filters) add_subdirectory(EventGenerator) +add_subdirectory(NGEM) 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) diff --git a/ubsim/NGEM/CMakeLists.txt b/ubsim/NGEM/CMakeLists.txt new file mode 100644 index 000000000..a50374748 --- /dev/null +++ b/ubsim/NGEM/CMakeLists.txt @@ -0,0 +1,25 @@ +cet_add_compiler_flags( CXX -DHIDE_GENIE_LOG_XXX ) + +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() +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..b995e3ea3 --- /dev/null +++ b/ubsim/NGEM/GENIEGenNGEM_module.cc @@ -0,0 +1,800 @@ +//////////////////////////////////////////////////////////////////////// +// +// +// 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_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 "nurandom/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 "nugen/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 "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" +#include "DeleteOneRandomPhoton.h" +#include "GenerateIsotropicSinglePhoton.h" +#include "ModifyParticle.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 + + 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", "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 + 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. + + }; +} + +namespace evgen{ + + //____________________________________________________________________________ + GENIEGenNGEM::GENIEGenNGEM(fhicl::ParameterSet const& pset) + : art::EDProducer(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())) + , 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)) + , fManuallyDecayModifiedPi0s(pset.get("ManuallyDecayModifiedPi0s", false)) + , fDeleteRandomModifiedGamma(pset.get("DeleteRandomModifiedGamma", false)) + { + 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->GDMLFile(), + 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; + 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.); + + 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), art::fullRun()); + + 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), art::subRunFragment()); + + 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; + } + } + + 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 (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 (fManuallyDecayModifiedPi0s){ //fModifyParticleVar.at(i_v) == "mass" && fModifyParticlePdg == 111) { + 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) { + 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/ModifyParticle.h b/ubsim/NGEM/ModifyParticle.h new file mode 100644 index 000000000..98cdac535 --- /dev/null +++ b/ubsim/NGEM/ModifyParticle.h @@ -0,0 +1,193 @@ +#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" + +#include "GenerateIsotropicSinglePhoton.h" + +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) { + + 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) { // 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_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(); + 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_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; + 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 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 + 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") { + double new_pE = new_part_var; + // protect against numerical issues / invalid input + if (new_pE < 0) { + std::cout << "Warning: Invalid particle parameters (negative energy or momentum)" << std::endl; + continue; + } + // Compute new |p| from E^2 = p^2 + m^2 + double new_p2 = new_pE*new_pE - part_mass*part_mass; + // 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 + 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(); + 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); + 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); + TLorentzVector mom(new_px, new_py, new_pz, pE); + new_part_momentum = mom; + } else if (modVar == "p") { + double new_p = new_part_var; + // protect against numerical issues / invalid input + if (new_p < 0) { + std::cout << "Warning: Invalid particle parameters (negative energy or momentum)" << std::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); + + // 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); + } 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(), 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); + 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(int modPdg, std::string modVar, + const std::vector& varBinEdges, const std::vector& varBinProbs, + int seed, simb::MCTruth& mcTruth) { + simb::MCTruth newTruth; + ModifyParticle(modPdg, modVar, varBinEdges, varBinProbs, seed, 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..1aa8c404d --- /dev/null +++ b/ubsim/NGEM/prodgenie_bnb_modifiedPi0s.fcl @@ -0,0 +1,33 @@ +#include "prodgenie_common_uboone_Filtered.fcl" + +process_name: GenieGenFiltered + +outputs.out1.fileName: "prodgenie_bnb_nu_filtered_modpi0_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 # 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 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) + +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 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]