Eyandel genie ngem mcc10

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)

ubsim/CMakeLists.txt

@@ -7,3 +7,4 @@ add_subdirectory(SNStreamSim)
 add_subdirectory(Simulation)
 add_subdirectory(Filters)
 add_subdirectory(EventGenerator)
+add_subdirectory(NGEM)

ubsim/EventWeight/Calculators/G4RWManagerService.cc

@@ -0,0 +1,18 @@
+#include "G4RWManagerService.h"
+
+evwgh::G4RWManagerService::G4RWManagerService(
+    fhicl::ParameterSet const& pset) :
+    fMaterial(pset.get<fhicl::ParameterSet>("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;
+}

ubsim/EventWeight/Calculators/G4RWManagerService_service.cc

@@ -1,20 +1,3 @@
 #include "G4RWManagerService.h"
 
-evwgh::G4RWManagerService::G4RWManagerService(
-    fhicl::ParameterSet const& pset) :
-    fMaterial(pset.get<fhicl::ParameterSet>("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)
+DEFINE_ART_SERVICE(evwgh::G4RWManagerService)

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()
+

ubsim/NGEM/DeleteOneRandomPhoton.h

@@ -0,0 +1,72 @@
+#ifndef EVGEN_DELETE_ONE_RANDOM_PHOTON_H
+#define EVGEN_DELETE_ONE_RANDOM_PHOTON_H
+
+#include <vector>
+#include <iostream>
+
+// 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<int> 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