diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2b088f9..6b75c51 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,78 +1,79 @@
 # Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 ################################################################################
 #    This file is part of NUISANCE.
 #
 #    NUISANCE is free software: you can redistribute it and/or modify
 #    it under the terms of the GNU General Public License as published by
 #    the Free Software Foundation, either version 3 of the License, or
 #    (at your option) any later version.
 #
 #    NUISANCE is distributed in the hope that it will be useful,
 #    but WITHOUT ANY WARRANTY; without even the implied warranty of
 #    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 #    GNU General Public License for more details.
 #
 #    You should have received a copy of the GNU General Public License
 #    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
 ################################################################################
 cmake_minimum_required (VERSION 2.8 FATAL_ERROR)
 
 #Use the compilers found in the path
 find_program(CMAKE_C_COMPILER NAMES $ENV{CC} gcc PATHS ENV PATH NO_DEFAULT_PATH)
 find_program(CMAKE_CXX_COMPILER NAMES $ENV{CXX} g++ PATHS ENV PATH NO_DEFAULT_PATH)
 
 project(NUISANCE)
 
 include(ExternalProject)
 
 set (NUISANCE_VERSION_MAJOR 3)
 set (NUISANCE_VERSION_MINOR 0)
 set (NUISANCE_VERSION_REVISION 0)
 
 set (NUISANCE_VERSION_STRING "v${NUISANCE_VERSION_MAJOR}r${NUISANCE_VERSION_MINOR}")
 if(${NUISANCE_VERSION_REVISION} STRGREATER "0")
   set (NUISANCE_VERSION_STRING "${NUISANCE_VERSION_STRING}p${NUISANCE_VERSION_REVISION}")
 endif()
 
 #Set this to TRUE to enable build debugging messages
 set(BUILD_DEBUG_MSGS TRUE)
 
 include(${CMAKE_SOURCE_DIR}/cmake/cmessage.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/cacheVariables.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/parseConfigApp.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/StringAndListUtils.cmake)
 
 cmessage(STATUS "CMAKE_INSTALL_PREFIX: \"${CMAKE_INSTALL_PREFIX}\"")
 cmessage(STATUS "CMAKE_BUILD_TYPE: \"${CMAKE_BUILD_TYPE}\"")
 
 ################################################################################
 #                            Check Dependencies
 ################################################################################
 
 ##################################  ROOT  ######################################
 include(${CMAKE_SOURCE_DIR}/cmake/ROOTSetup.cmake)
 
 ################################# InputHandler #################################
 include(${CMAKE_SOURCE_DIR}/cmake/InputHandlerSetup.cmake)
 
 #################################  Pythia6/8 ###################################
 include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake)
 
 #Want this before fhiclcpp which will add the install directory
 include_directories(${CMAKE_SOURCE_DIR}/src)
 
 ################################## FHICLCPP ####################################
 include(${CMAKE_SOURCE_DIR}/cmake/fhiclcppSetup.cmake)
 
 ################################## COMPILER ####################################
+include(${CMAKE_SOURCE_DIR}/cmake/gperfSetup.cmake)
 include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake)
 
 ################################################################################
 
 add_subdirectory(src)
 
 add_subdirectory(scripts)
 
 add_subdirectory(config)
 add_subdirectory(data)
diff --git a/cmake/gperfSetup.cmake b/cmake/gperfSetup.cmake
new file mode 100644
index 0000000..3d54519
--- /dev/null
+++ b/cmake/gperfSetup.cmake
@@ -0,0 +1,60 @@
+# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
+
+################################################################################
+#    This file is part of NUISANCE.
+#
+#    NUISANCE is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation, either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    NUISANCE is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License
+#    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
+################################################################################
+
+if(DEFINED USE_GPERFTOOLS AND USE_GPERFTOOLS)
+
+  ExternalProject_Add(libunwind
+    PREFIX "${CMAKE_BINARY_DIR}/Ext"
+    GIT_REPOSITORY "https://github.com/libunwind/libunwind.git"
+    GIT_TAG v1.3.1
+    CONFIGURE_COMMAND ${CMAKE_BINARY_DIR}/Ext/src/libunwind/autogen.sh --prefix=${CMAKE_INSTALL_PREFIX}
+    UPDATE_DISCONNECTED 1
+    UPDATE_COMMAND ""
+    BUILD_COMMAND make
+    INSTALL_COMMAND make install
+    )
+
+  ExternalProject_Add(gperftools
+    PREFIX "${CMAKE_BINARY_DIR}/Ext"
+    GIT_REPOSITORY "https://github.com/gperftools/gperftools.git"
+    GIT_TAG "gperftools-2.7"
+    CONFIGURE_COMMAND ./autogen.sh && ./configure --prefix=${CMAKE_INSTALL_PREFIX} CPPFLAGS=-I${CMAKE_INSTALL_PREFIX}/include LDFLAGS=-L${CMAKE_INSTALL_PREFIX}/lib
+    BUILD_IN_SOURCE 1
+    UPDATE_DISCONNECTED 1
+    UPDATE_COMMAND ""
+    BUILD_COMMAND make
+    INSTALL_COMMAND make install
+    )
+
+  add_dependencies(gperftools libunwind)
+
+  LIST(APPEND EXTRA_CXX_FLAGS
+    -fno-builtin-malloc
+    -fno-builtin-calloc
+    -fno-builtin-realloc
+    -fno-builtin-free)
+
+
+  LIST(APPEND EXTRA_LINK_DIRS ${CMAKE_INSTALL_PREFIX}/lib)
+
+  ##Want to prepend them
+  LIST(APPEND EXTRA_LIBS tcmalloc_and_profiler)
+
+  cmessage(STATUS "Using google performance libraries")
+endif()
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 63767ba..68b52bd 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,18 +1,19 @@
 add_subdirectory(config)
 add_subdirectory(event)
 add_subdirectory(exception)
 add_subdirectory(generator)
 add_subdirectory(input)
 add_subdirectory(parameters)
 add_subdirectory(persistency)
 add_subdirectory(plugins)
 add_subdirectory(samples)
 add_subdirectory(utility)
 add_subdirectory(variation)
 
 SET(INuADataComparisons_List ${INuADataComparisons_List} PARENT_SCOPE)
 SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE)
 
 SET(GENERATOR_DEPENDENT_TARGETS ${GENERATOR_DEPENDENT_TARGETS} PARENT_SCOPE)
 
+add_subdirectory(tests)
 add_subdirectory(app)
diff --git a/src/app/CMakeLists.txt b/src/app/CMakeLists.txt
index c84c5fe..34feffb 100644
--- a/src/app/CMakeLists.txt
+++ b/src/app/CMakeLists.txt
@@ -1,15 +1,15 @@
 SET(APPS nuissamples nuiscomp nuisevsum nuisstudy)
 
 foreach(a ${APPS})
   add_executable(${a} ${a}.cxx)
-  target_link_libraries(${a} -Wl,--no-as-needed nuis_event nuis_input nuis_utility_experimental nuis_utility nuis_config nuis_persistency nuis_plugins nuis_params nuis_variation)
+  target_link_libraries(${a} -Wl,--no-as-needed nuis_event nuis_input nuis_generator_utility nuis_utility_experimental nuis_utility nuis_config nuis_persistency nuis_plugins nuis_params nuis_variation)
   target_link_libraries(${a} ${GENERATOR_DEPENDENT_TARGETS} -Wl,--as-needed)
   target_link_libraries(${a} ${NUISANCE_LINK_DIRS})
   target_link_libraries(${a} ${NUISANCE_DEPEND_LIBS})
 
   if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
     set_target_properties(${a} PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
   endif()
 
   install(TARGETS ${a} DESTINATION bin)
 endforeach()
diff --git a/src/generator/input/GENIEInputHandler.cxx b/src/generator/input/GENIEInputHandler.cxx
index cc5c146..1bfb67d 100644
--- a/src/generator/input/GENIEInputHandler.cxx
+++ b/src/generator/input/GENIEInputHandler.cxx
@@ -1,165 +1,215 @@
 #include "generator/input/GENIEInputHandler.hxx"
 
 #include "generator/utility/GENIEUtility.hxx"
 
 #include "utility/PDGCodeUtility.hxx"
 
 #ifdef GENIE_V3_INTERFACE
 #include "Framework/EventGen/EventRecord.h"
 #include "Framework/GHEP/GHepParticle.h"
 #include "Framework/GHEP/GHepRecord.h"
 #else
 #include "EVGCore/EventRecord.h"
 #include "GHEP/GHepParticle.h"
 #include "GHEP/GHepRecord.h"
 #endif
 
 #include "fhiclcpp/ParameterSet.h"
 
 using namespace nuis;
 using namespace nuis::event;
 using namespace nuis::utility;
 using namespace nuis::genietools;
 
 GENIEInputHandler::GENIEInputHandler()
     : fInputTree(), fGenieNtpl(nullptr), fFileWeight(1) {}
 GENIEInputHandler::GENIEInputHandler(GENIEInputHandler &&other)
     : fInputTree(std::move(other.fInputTree)),
       fReaderEvent(std::move(other.fReaderEvent)), fGenieNtpl(nullptr),
       fFileWeight(1) {}
 
+std::set<std::string> GENIEInputHandler::GetSplineList() {
+  std::set<std::string> splinenames;
+  size_t NEvents = 1000; // GetNEvents();
+  size_t ShoutEvery = NEvents / 100;
+  std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " GENIE events."
+            << std::flush;
+  for (size_t ev_it = 0; ev_it < NEvents; ++ev_it) {
+    if (fGenieNtpl) {
+      fGenieNtpl->Clear();
+    }
+    fInputTree.tree->GetEntry(ev_it);
+
+    if (ShoutEvery && !(ev_it % ShoutEvery)) {
+      std::cout << "\r[INFO]: Read " << ev_it << "/" << NEvents
+                << " GENIE events." << std::flush;
+    }
+
+    genie::GHepRecord *GHep =
+        static_cast<genie::GHepRecord *>(fGenieNtpl->event);
+    if (!GHep) {
+      throw invalid_GENIE_event() << "[ERROR]: GENIE event " << ev_it
+                                  << " failed to contain a GHepRecord";
+    }
+    splinenames.insert(GHep->Summary()->AsString());
+  }
+  std::cout << "[INFO]: Having read GENIE events, found " << splinenames.size()
+            << " splines: " << std::endl;
+  for (auto const &sn : splinenames) {
+    std::cout << "\t" << sn << std::endl;
+  }
+  return splinenames;
+}
+
 void GENIEInputHandler::Initialize(fhicl::ParameterSet const &ps) {
 
   fInputTree = CheckGetTTree(ps.get<std::string>("file"), "gtree");
 
   fInputTree.tree->SetBranchAddress("gmcrec", &fGenieNtpl);
 
   fKeepIntermediates = ps.get<bool>("keep_intermediates", false);
   fKeepNuclearParticles = ps.get<bool>("keep_nuclear_particles", false);
 
   fhicl::ParameterSet XSecInfo = ps.get<fhicl::ParameterSet>("xsec_info", {});
   if (XSecInfo.has_key("weight")) {
     fFileWeight = XSecInfo.get<double>("weight");
     std::cout << "[INFO]: Average GENIE XSecWeight = " << fFileWeight << ""
               << std::endl;
   } else if (XSecInfo.has_key("flux") || XSecInfo.has_key("target") ||
-             XSecInfo.has_key("spline_file")) {
+             XSecInfo.has_key("spline_file") ||
+             XSecInfo.has_key("filter_splines_by_event_tree")) {
     std::cout
         << "[INFO]: Attempting to build GENIE file weight with input splines."
         << std::endl;
-    fFileWeight = genietools::GetFileWeight(XSecInfo) / double(GetNEvents());
+
+    if (XSecInfo.get<bool>("filter_splines_by_event_tree", false)) {
+      std::cout << "[INFO]: Ensuring we only read relevant splines... this "
+                   "will take a while."
+                << std::endl;
+      fFileWeight = genietools::GetFileWeight(XSecInfo, GetSplineList()) /
+                    double(GetNEvents());
+    } else {
+      fFileWeight = genietools::GetFileWeight(XSecInfo) / double(GetNEvents());
+    }
     std::cout << "[INFO]: Average GENIE XSecWeight = " << fFileWeight << ""
               << std::endl;
   }
 }
 
-MinimalEvent const &GENIEInputHandler::GetMinimalEvent(ev_index_t idx) const {
+void GENIEInputHandler::GetEntry(ev_index_t idx) const {
   if (idx >= GetNEvents()) {
     throw IInputHandler::invalid_entry()
         << "[ERROR]: Attempted to get entry " << idx
         << " from an InputHandler with only " << GetNEvents();
   }
 
+  if (fGenieNtpl) {
+    fGenieNtpl->Clear();
+  }
   fInputTree.tree->GetEntry(idx);
+}
+
+MinimalEvent const &GENIEInputHandler::GetMinimalEvent(ev_index_t idx) const {
+  GetEntry(idx);
 
   genie::GHepRecord *GHep = static_cast<genie::GHepRecord *>(fGenieNtpl->event);
   if (!GHep) {
     throw invalid_GENIE_event()
         << "[ERROR]: GENIE event " << idx << " failed to contain a GHepRecord";
   }
 
   fReaderEvent.mode = GetEventChannel(*GHep);
 
   TObjArrayIter iter(GHep);
   genie::GHepParticle *p;
   while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
     if (!p) {
       continue;
     }
 
     Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode);
     if (state != Particle::Status_t::kPrimaryInitialState) {
       continue;
     }
     if (!IsNeutralLepton(p->Pdg(), pdgcodes::kMatterAntimatter) &&
         !IsChargedLepton(p->Pdg(), pdgcodes::kMatterAntimatter)) {
       continue;
     }
     fReaderEvent.probe_E = p->E() * 1.E3;
     fReaderEvent.probe_pdg = p->Pdg();
     break;
   }
 
   fReaderEvent.XSecWeight = fFileWeight;
 
   if (fWeightCache.size() <= idx) {
     fWeightCache.push_back(fReaderEvent.XSecWeight);
   }
 
   return fReaderEvent;
 }
 
 FullEvent const &GENIEInputHandler::GetFullEvent(ev_index_t idx) const {
   (void)GetMinimalEvent(idx);
 
   fReaderEvent.ClearParticleStack();
 
   genie::GHepRecord *GHep = static_cast<genie::GHepRecord *>(fGenieNtpl->event);
 
   unsigned int npart = GHep->GetEntries();
 
   // Fill Particle Stack
   genie::GHepParticle *p = 0;
   TObjArrayIter iter(GHep);
 
   // Loop over all particles
   while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
     if (!p) {
       continue;
     }
 
     // Get Status
     Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode);
 
     if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) {
       continue;
     }
 
     if (!fKeepNuclearParticles && IsNuclearPDG(p->Pdg())) {
       continue;
     }
 
     Particle nuis_part;
     nuis_part.pdg = p->Pdg();
     nuis_part.P4 = TLorentzVector(p->Px(), p->Py(), p->Pz(), p->E()) * 1E3;
 
     fReaderEvent.ParticleStack[static_cast<size_t>(state)].particles.push_back(
         nuis_part);
   }
 
   return fReaderEvent;
 }
 
 double GENIEInputHandler::GetEventWeight(ev_index_t idx) const {
   if (idx > fWeightCache.size()) {
     throw weight_cache_miss()
         << "[ERROR]: Failed to get cached weight for event index: " << idx;
   }
   return fWeightCache[idx];
 }
 
 double GENIEInputHandler::GetXSecScaleFactor(
     std::pair<double, double> const &EnuRange) const {
   throw input_handler_feature_unimplemented()
       << "[ERROR]: Flux cuts not yet implemented for GENIE input handler.";
 }
 
 size_t GENIEInputHandler::GetNEvents() const {
   return fInputTree.tree->GetEntries();
 }
 
 GeneratorManager::Generator_id_t GENIEInputHandler::GetGeneratorId() const {
   return GeneratorManager::Get().EnsureGeneratorRegistered("GENIE");
 }
 
 DECLARE_PLUGIN(IInputHandler, GENIEInputHandler);
diff --git a/src/generator/input/GENIEInputHandler.hxx b/src/generator/input/GENIEInputHandler.hxx
index 312b6b7..a1f3456 100644
--- a/src/generator/input/GENIEInputHandler.hxx
+++ b/src/generator/input/GENIEInputHandler.hxx
@@ -1,69 +1,72 @@
 // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is part of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 
 #pragma once
 
 #include "event/FullEvent.hxx"
 
 #include "input/IInputHandler.hxx"
 
 #include "exception/exception.hxx"
 
 #include "utility/ROOTUtility.hxx"
 
 #ifdef GENIE_V3_INTERFACE
 #include "Framework/Ntuple/NtpMCEventRecord.h"
 #else
 #include "Ntuple/NtpMCEventRecord.h"
 #endif
 
 namespace fhicl {
   class ParameterSet;
 }
 
 #include <memory>
+#include <set>
 
 class GENIEInputHandler : public IInputHandler {
   mutable nuis::utility::TreeFile fInputTree;
   mutable nuis::event::FullEvent fReaderEvent;
   mutable std::vector<double> fWeightCache;
   mutable genie::NtpMCEventRecord *fGenieNtpl;
 
   bool fKeepIntermediates;
   bool fKeepNuclearParticles;
 
   double fFileWeight;
 
 public:
   NEW_NUIS_EXCEPT(weight_cache_miss);
 
   GENIEInputHandler();
   GENIEInputHandler(GENIEInputHandler const &) = delete;
   GENIEInputHandler(GENIEInputHandler &&);
 
   void Initialize(fhicl::ParameterSet const &);
+  void GetEntry(ev_index_t) const;
   nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const;
   nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const;
   double GetEventWeight(ev_index_t idx) const;
   size_t GetNEvents() const;
+  std::set<std::string> GetSplineList();
 
   double GetXSecScaleFactor(std::pair<double, double> const &EnuRange) const;
 
   nuis::GeneratorManager::Generator_id_t GetGeneratorId() const;
 };
diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx
index 9fc688d..476b383 100644
--- a/src/generator/input/NEUTInputHandler.cxx
+++ b/src/generator/input/NEUTInputHandler.cxx
@@ -1,260 +1,264 @@
 #include "generator/input/NEUTInputHandler.hxx"
 
 #include "persistency/ROOTOutput.hxx"
 
 #include "utility/HistogramUtility.hxx"
 #include "utility/InteractionChannelUtility.hxx"
 #include "utility/PDGCodeUtility.hxx"
 
 #include "generator/utility/NEUTUtility.hxx"
 
 #include "neutpart.h"
 #include "neutvect.h"
 
 #include "fhiclcpp/ParameterSet.h"
 
 using namespace nuis;
 using namespace nuis::event;
 using namespace nuis::utility;
 using namespace nuis::neuttools;
 
 NEUTInputHandler::NEUTInputHandler() : fInputTreeFile() {}
 NEUTInputHandler::NEUTInputHandler(NEUTInputHandler &&other)
     : fInputTreeFile(std::move(other.fInputTreeFile)),
       fReaderEvent(std::move(other.fReaderEvent)) {}
 
 void NEUTInputHandler::Initialize(fhicl::ParameterSet const &ps) {
 
   fInputTreeFile = CheckGetTTree(ps.get<std::string>("file"), "neuttree");
 
   fNeutVect = nullptr;
   fInputTreeFile.tree->SetBranchAddress("vectorbranch", &fNeutVect);
 
   fKeepIntermediates = ps.get<bool>("keep_intermediates", false);
 
   std::string flux_name = ps.get<std::string>("flux_hist_name", "flux_numu");
   std::string evtrt_name = ps.get<std::string>("evtrt_hist_name", "evtrt_numu");
   fXSecRescaleFactor = ps.get<double>("xsec_rescale_factor",1);
 
   std::string override_flux_file =
       ps.get<std::string>("override_flux_file", "");
 
   if (override_flux_file.size()) {
     fFlux = GetHistogramFromROOTFile<TH1>(override_flux_file, flux_name);
     RebuildEventRate(!ps.get<bool>("flux_hist_in_MeV", false));
   } else {
     fFlux =
         GetHistogramFromROOTFile<TH1>(fInputTreeFile.file, flux_name, false);
   }
 
   if (fFlux) {
     if (!fEvtrt) {
       fEvtrt = GetHistogramFromROOTFile<TH1>(fInputTreeFile.file, evtrt_name);
     }
 
     fFileWeight = fEvtrt->Integral() * 1.0E-38 /
                   (fFlux->Integral() * double(GetNEvents()));
     std::cout << "[INFO]: Average NEUT XSecWeight = " << fFileWeight << ""
               << std::endl;
   } else {
     std::cout
         << "[INFO]: Input NEUT file doesn't have flux information, assuming "
            "mono energetic and calculating average cross-section..."
         << std::endl;
     fFileWeight = GetMonoEXSecWeight() / double(GetNEvents());
     std::cout << "[INFO]: Done (Average NEUT XSecWeight = " << fFileWeight
               << ")!" << std::endl;
   }
 
   fFileWeight *= fXSecRescaleFactor;
 }
 
 NEW_NUIS_EXCEPT(non_mono_energetic_input_file);
 double NEUTInputHandler::GetMonoEXSecWeight() {
   double xsec = 0;
   double count = 0;
   double E_first = std::numeric_limits<double>::max();
   size_t NEvents = GetNEvents();
   size_t ShoutEvery = NEvents / 100;
   std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " NEUT events."
             << std::flush;
 
   for (size_t nev_it = 0; nev_it < NEvents; ++nev_it) {
     if (ShoutEvery && !(nev_it % ShoutEvery)) {
       std::cout << "\r[INFO]: Read " << nev_it << "/" << NEvents
                 << " NEUT events." << std::flush;
     }
 
     fInputTreeFile.tree->GetEntry(nev_it);
     xsec += fNeutVect->Totcrs;
     //Assume monoE.
     return xsec * 1E-38;
 
     count++;
     if (E_first == std::numeric_limits<double>::max()) {
       NeutPart *part = fNeutVect->PartInfo(0);
       E_first = part->fP.E();
     } else {
       NeutPart *part = fNeutVect->PartInfo(0);
       if (E_first != part->fP.E()) {
         throw non_mono_energetic_input_file()
             << "[ERROR]: When calculating NEUT xsec weight for a mono "
                "energetic beam, found different energies: "
             << E_first << " and " << part->fP.E()
             << ". Cannot continue with this input file.";
       }
     }
   }
   std::cout << "\r[INFO]: Read " << NEvents << "/" << NEvents << " NEUT events."
             << std::endl;
 
   return (xsec / count) * 1E-38;
 }
 
 void NEUTInputHandler::RebuildEventRate(bool FluxInGeV) {
   auto XSec = Clone(fFlux, true, "xsec");
   auto Entry = Clone(fFlux, true, "entry");
 
   std::cout << "[INFO]: Rebuilding total cross-section prediction..."
             << std::endl;
   size_t NEvents = GetNEvents();
   size_t ShoutEvery = NEvents / 100;
   std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " NEUT events."
             << std::flush;
 
   for (size_t nev_it = 0; nev_it < NEvents; ++nev_it) {
     if (ShoutEvery && !(nev_it % ShoutEvery)) {
       std::cout << "\r[INFO]: Read " << nev_it << "/" << NEvents
                 << " NEUT events." << std::flush;
     }
 
     fInputTreeFile.tree->GetEntry(nev_it);
 
     NeutPart *part = fNeutVect->PartInfo(0);
     double E = part->fP.E() * (FluxInGeV ? 1E-3 : 1);
 
     XSec->Fill(E, fNeutVect->Totcrs);
     Entry->Fill(E);
   }
   std::cout << std::endl << "[INFO]: Done!" << std::endl;
 
   // Binned total xsec
   XSec->Divide(Entry.get());
 
   fEvtrt = Clone(XSec, false, "evtrt");
   // Event rate prediction
   fEvtrt->Multiply(fFlux.get());
 }
 
 double
 NEUTInputHandler::GetXSecScaleFactor(std::pair<double, double> const &) const {
 
   throw input_handler_feature_unimplemented()
       << "[ERROR]: Flux cuts not yet implemented for NEUT input handler.";
 }
 
 NeutVect *NEUTInputHandler::GetNeutEvent(ev_index_t idx) const {
   if (idx >= GetNEvents()) {
     throw IInputHandler::invalid_entry()
         << "[ERROR]: Attempted to get entry " << idx
         << " from an InputHandler with only " << GetNEvents();
   }
   fInputTreeFile.tree->GetEntry(idx);
 
   return fNeutVect;
 }
 
-MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const {
+void NEUTInputHandler::GetEntry(ev_index_t idx) const {
   if (idx >= GetNEvents()) {
     throw IInputHandler::invalid_entry()
         << "[ERROR]: Attempted to get entry " << idx
         << " from an InputHandler with only " << GetNEvents();
   }
   fInputTreeFile.tree->GetEntry(idx);
+}
+
+MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const {
+  GetEntry(idx);
 
   fReaderEvent.mode = IntToChannel(fNeutVect->Mode);
 
   size_t NPart = fNeutVect->Npart();
   for (size_t part_it = 0; part_it < NPart; part_it++) {
     NeutPart const &part = (*fNeutVect->PartInfo(part_it));
     if ((part.fIsAlive == false) && (part.fStatus == -1) &&
         IsNeutralLepton(part.fPID)) {
       fReaderEvent.probe_E = part.fP.T();
       fReaderEvent.probe_pdg = part.fPID;
       break;
     }
   }
 
   fReaderEvent.XSecWeight = fFileWeight;
 
   if (fWeightCache.size() <= idx) {
     fWeightCache.push_back(fReaderEvent.XSecWeight);
   }
 
   fReaderEvent.fGenEvent = fNeutVect;
 
   return fReaderEvent;
 }
 
 FullEvent const &NEUTInputHandler::GetFullEvent(ev_index_t idx) const {
   (void)GetMinimalEvent(idx);
 
   fReaderEvent.ClearParticleStack();
 
   if (fNeutVect->Ibound) {
     fReaderEvent.target_pdg =
         MakeNuclearPDG(fNeutVect->TargetA, fNeutVect->TargetZ);
   } else {
     fReaderEvent.target_pdg = MakeNuclearPDG(1, 1);
   }
 
   size_t NPart = fNeutVect->Npart();
   size_t NPrimary = fNeutVect->Nprimary();
   for (size_t part_it = 0; part_it < NPart; part_it++) {
     NeutPart const &part = (*fNeutVect->PartInfo(part_it));
 
     Particle nuis_part;
     nuis_part.pdg = part.fPID;
     nuis_part.P4 = part.fP;
 
     Particle::Status_t state = GetNeutParticleStatus(part, fReaderEvent.mode);
 
     if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) {
       continue;
     }
 
     size_t state_int = static_cast<size_t>(state);
 
     // Add status == 0 particles as pre-FSI particles until we find an
     // intermediate state particle
     if ((part_it < NPrimary) &&
         (state != Particle::Status_t::kPrimaryInitialState)) {
       fReaderEvent
           .ParticleStack[static_cast<size_t>(
               Particle::Status_t::kPrimaryFinalState)]
           .particles.push_back(nuis_part);
     }
 
     fReaderEvent.ParticleStack[state_int].particles.push_back(nuis_part);
   }
 
   return fReaderEvent;
 }
 
 double NEUTInputHandler::GetEventWeight(ev_index_t idx) const {
   if (idx > fWeightCache.size()) {
     throw weight_cache_miss()
         << "[ERROR]: Failed to get cached weight for event index: " << idx;
   }
   return fWeightCache[idx];
 }
 
 size_t NEUTInputHandler::GetNEvents() const {
   return fInputTreeFile.tree->GetEntries();
 }
 
 GeneratorManager::Generator_id_t NEUTInputHandler::GetGeneratorId() const {
   return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT");
 }
 
 DECLARE_PLUGIN(IInputHandler, NEUTInputHandler);
diff --git a/src/generator/input/NEUTInputHandler.hxx b/src/generator/input/NEUTInputHandler.hxx
index a2cc915..574e2fa 100644
--- a/src/generator/input/NEUTInputHandler.hxx
+++ b/src/generator/input/NEUTInputHandler.hxx
@@ -1,71 +1,72 @@
 // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is part of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 
 #pragma once
 
 #include "event/FullEvent.hxx"
 
 #include "input/IInputHandler.hxx"
 
 #include "utility/ROOTUtility.hxx"
 
 #include <memory>
 
 class NeutVect;
 class TH1;
 
 namespace fhicl {
 class ParameterSet;
 }
 
 class NEUTInputHandler : public IInputHandler {
   mutable nuis::utility::TreeFile fInputTreeFile;
   mutable nuis::event::FullEvent fReaderEvent;
   mutable NeutVect *fNeutVect;
   mutable std::vector<double> fWeightCache;
   double fXSecRescaleFactor;
 
   std::unique_ptr<TH1> fFlux;
   std::unique_ptr<TH1> fEvtrt;
   double fFileWeight;
 
   bool fKeepIntermediates;
 
   double GetMonoEXSecWeight();
   void RebuildEventRate(bool FluxInGeV=true);
 
 public:
   NEW_NUIS_EXCEPT(weight_cache_miss);
 
   NEUTInputHandler();
   NEUTInputHandler(NEUTInputHandler const &) = delete;
   NEUTInputHandler(NEUTInputHandler &&);
 
   void Initialize(fhicl::ParameterSet const &);
   NeutVect *GetNeutEvent(ev_index_t) const;
+  void GetEntry(ev_index_t) const;
   nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const;
   nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const;
 
   double GetEventWeight(ev_index_t idx) const;
   size_t GetNEvents() const;
 
   double GetXSecScaleFactor(std::pair<double, double> const &EnuRange) const;
 
   nuis::GeneratorManager::Generator_id_t GetGeneratorId() const;
 };
diff --git a/src/generator/input/NuWroInputHandler.cxx b/src/generator/input/NuWroInputHandler.cxx
index 15e3814..ae0d48c 100644
--- a/src/generator/input/NuWroInputHandler.cxx
+++ b/src/generator/input/NuWroInputHandler.cxx
@@ -1,134 +1,139 @@
 #include "generator/input/NuWroInputHandler.hxx"
 
 #include "generator/utility/NuWroUtility.hxx"
 
 #include "fhiclcpp/ParameterSet.h"
 
 #include "particle.h"
 using NuWroParticle = ::particle;
 
 using namespace nuis;
 using namespace nuis::event;
 using namespace nuis::utility;
 using namespace nuis::nuwrotools;
 
 NuWroInputHandler::NuWroInputHandler() : fInputTree(), fTreeEvent(nullptr) {}
 NuWroInputHandler::NuWroInputHandler(NuWroInputHandler &&other)
     : fInputTree(std::move(other.fInputTree)),
       fReaderEvent(std::move(other.fReaderEvent)),
       fTreeEvent(other.fTreeEvent) {}
 
 void NuWroInputHandler::Initialize(fhicl::ParameterSet const &ps) {
 
   fInputTree = CheckGetTTree(ps.get<std::string>("file"), "treeout");
 
   fTreeEvent = nullptr;
   fInputTree.tree->SetBranchAddress("e", &fTreeEvent);
 
   fKeepIntermediates = ps.get<bool>("keep_intermediates", false);
 
-  if(GetNEvents()){
+  if (GetNEvents()) {
     fInputTree.tree->GetEntry(0);
-    std::cout << "[INFO]: Average NuWro XSec weight: " << (fTreeEvent->weight / double(GetNEvents())) << std::endl;
+    std::cout << "[INFO]: Average NuWro XSec weight: "
+              << (fTreeEvent->weight / double(GetNEvents())) << std::endl;
   }
 }
-MinimalEvent const &NuWroInputHandler::GetMinimalEvent(ev_index_t idx) const {
+
+void NuWroInputHandler::GetEntry(ev_index_t idx) const {
   if (idx >= GetNEvents()) {
     throw IInputHandler::invalid_entry()
         << "[ERROR]: Attempted to get entry " << idx
         << " from an InputHandler with only " << GetNEvents();
   }
 
   fInputTree.tree->GetEntry(idx);
+}
+
+MinimalEvent const &NuWroInputHandler::GetMinimalEvent(ev_index_t idx) const {
 
   fReaderEvent.mode = NuWroEventChannel(*fTreeEvent);
   fReaderEvent.probe_E = fTreeEvent->in[0].E();
   fReaderEvent.probe_pdg = fTreeEvent->in[0].pdg;
   fReaderEvent.XSecWeight = fTreeEvent->weight / double(GetNEvents());
 
   if (fWeightCache.size() <= idx) {
     fWeightCache.push_back(fReaderEvent.XSecWeight);
   }
 
   return fReaderEvent;
 }
 
 FullEvent const &NuWroInputHandler::GetFullEvent(ev_index_t idx) const {
   (void)GetMinimalEvent(idx);
 
   fReaderEvent.ClearParticleStack();
 
   for (size_t p_it = 0; p_it < fTreeEvent->in.size(); ++p_it) {
     NuWroParticle &part = fTreeEvent->in[p_it];
 
     Particle nuis_part;
 
     nuis_part.pdg = part.pdg;
     nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]);
 
     fReaderEvent
         .ParticleStack[static_cast<size_t>(
             Particle::Status_t::kPrimaryInitialState)]
         .particles.push_back(nuis_part);
   }
 
   for (size_t p_it = 0; p_it < fKeepIntermediates && fTreeEvent->out.size();
        ++p_it) {
     NuWroParticle &part = fTreeEvent->out[p_it];
 
     Particle nuis_part;
 
     nuis_part.pdg = part.pdg;
     nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]);
 
     fReaderEvent
         .ParticleStack[static_cast<size_t>(
             Particle::Status_t::kPrimaryFinalState)]
         .particles.push_back(nuis_part);
   }
   for (size_t p_it = 0; (p_it < fTreeEvent->post.size()); ++p_it) {
     NuWroParticle &part = fTreeEvent->post[p_it];
 
     Particle nuis_part;
 
     nuis_part.pdg = part.pdg;
     nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]);
 
     fReaderEvent
         .ParticleStack[static_cast<size_t>(Particle::Status_t::kNuclearLeaving)]
         .particles.push_back(nuis_part);
   }
 
   fReaderEvent.fGenEvent = fTreeEvent;
 
   return fReaderEvent;
 }
 
 NuWroEvent const &NuWroInputHandler::GetNuWroEvent(ev_index_t idx) const {
   (void)GetMinimalEvent(idx);
   return *fTreeEvent;
 }
 
 double NuWroInputHandler::GetEventWeight(ev_index_t idx) const {
   if (idx > fWeightCache.size()) {
     throw weight_cache_miss()
         << "[ERROR]: Failed to get cached weight for event index: " << idx;
   }
   return fWeightCache[idx];
 }
 
 double NuWroInputHandler::GetXSecScaleFactor(
     std::pair<double, double> const &EnuRange) const {
   throw input_handler_feature_unimplemented()
       << "[ERROR]: Flux cuts not yet implemented for NuWro input handler.";
 }
 
 size_t NuWroInputHandler::GetNEvents() const {
   return fInputTree.tree->GetEntries();
 }
 
 GeneratorManager::Generator_id_t NuWroInputHandler::GetGeneratorId() const {
   return GeneratorManager::Get().EnsureGeneratorRegistered("NuWro");
 }
 
 DECLARE_PLUGIN(IInputHandler, NuWroInputHandler);
diff --git a/src/generator/input/NuWroInputHandler.hxx b/src/generator/input/NuWroInputHandler.hxx
index ec36b0c..143a39c 100644
--- a/src/generator/input/NuWroInputHandler.hxx
+++ b/src/generator/input/NuWroInputHandler.hxx
@@ -1,66 +1,67 @@
 // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is part of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 
 #pragma once
 
 #include "event/FullEvent.hxx"
 
 #include "input/IInputHandler.hxx"
 
 #include "exception/exception.hxx"
 
 #include "utility/ROOTUtility.hxx"
 
 #include "event1.h"
 using NuWroEvent = ::event;
 
 #include <memory>
 
 namespace fhicl {
 class ParameterSet;
 }
 
 class NuWroInputHandler : public IInputHandler {
   mutable nuis::utility::TreeFile fInputTree;
   mutable nuis::event::FullEvent fReaderEvent;
   mutable std::vector<double> fWeightCache;
   mutable NuWroEvent *fTreeEvent;
 
   bool fKeepIntermediates;
 
 public:
   NEW_NUIS_EXCEPT(weight_cache_miss);
 
   NuWroInputHandler();
   NuWroInputHandler(NuWroInputHandler const &) = delete;
   NuWroInputHandler(NuWroInputHandler &&);
 
   void Initialize(fhicl::ParameterSet const &);
+  void GetEntry(ev_index_t) const;
   nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const;
   nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const;
 
   NuWroEvent const &GetNuWroEvent(ev_index_t idx) const;
 
   double GetEventWeight(ev_index_t idx) const;
   size_t GetNEvents() const;
 
   double GetXSecScaleFactor(std::pair<double, double> const &EnuRange) const;
 
   nuis::GeneratorManager::Generator_id_t GetGeneratorId() const;
 };
diff --git a/src/generator/utility/GENIEUtility.cxx b/src/generator/utility/GENIEUtility.cxx
index b3be697..6011309 100644
--- a/src/generator/utility/GENIEUtility.cxx
+++ b/src/generator/utility/GENIEUtility.cxx
@@ -1,465 +1,479 @@
 // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is part of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 
 #include "generator/utility/GENIEUtility.hxx"
 #include "generator/utility/GENIESplineReader.hxx"
 #include "utility/HistogramUtility.hxx"
-#include "utility/StringUtility.hxx"
-
 #include "utility/InteractionChannelUtility.hxx"
 #include "utility/PDGCodeUtility.hxx"
+#include "utility/ROOTUtility.hxx"
+#include "utility/StringUtility.hxx"
 
 #ifdef GENIE_V3_INTERFACE
 #include "Framework/GHEP/GHepParticle.h"
 #include "Framework/GHEP/GHepRecord.h"
 #include "Framework/GHEP/GHepUtils.h"
 #include "Framework/ParticleData/PDGCodes.h"
 #else
 #include "GHEP/GHepParticle.h"
 #include "GHEP/GHepRecord.h"
 #include "GHEP/GHepUtils.h"
 #include "PDG/PDGCodes.h"
 #endif
 
 #include "fhiclcpp/ParameterSet.h"
 
 #include "string_parsers/from_string.hxx"
 
 #include "TGraph.h"
 
 namespace nuis {
 
 using namespace event;
 
 namespace genietools {
 
-static std::map<std::string, TGraph> SplineCache;
-
-TGraph dum;
-
-TGraph const &GetGENIESpline(std::string const &SplineFile,
-                             std::string const &SplineIdentifier) {
-  if (SplineCache.find(SplineFile + SplineIdentifier) != SplineCache.end()) {
-    return SplineCache.find(SplineFile + SplineIdentifier)->second;
-  }
-  return dum;
-}
-
 struct NFSParticleCount {
   size_t NProton;
   size_t NNeutron;
   size_t NPip;
   size_t NPi0;
   size_t NPim;
   size_t NOther;
 };
 
 NFSParticleCount CountPreFSIParticles(genie::GHepRecord const &ev) {
   // This code in this method is adapted from the GENIE source code found in
   // GHep/GHepUtils.cxx This method therefore carries the GENIE copyright
   // licence as copied below:
   //
   /// Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration
   /// For the full text of the license visit http://copyright.genie-mc.org
   /// or see $GENIE/LICENSE
   //
 
   genie::Target const &tgt = ev.Summary()->InitState().Tgt();
   if (!tgt.HitNucIsSet()) {
     throw invalid_GENIE_event()
         << "[ERROR]: Failed to get hit nucleon kinematics as it was not "
            "included in this GHep event. This is a fatal error.";
   }
 
   genie::GHepParticle *FSLep = ev.FinalStatePrimaryLepton();
   genie::GHepParticle *ISLep = ev.Probe();
 
   if (!FSLep || !ISLep) {
     throw invalid_GENIE_event()
         << "[ERROR]: Failed to find IS and FS lepton in event: "
         << ev.Summary()->AsString();
   }
 
   size_t NPi0 = 0, NPip = 0, NPim = 0, NProton = 0, NNeutron = 0, NOther = 0;
 
   bool nuclear_target = tgt.IsNucleus();
 
   TIter event_iter(&ev);
   genie::GHepParticle *p = 0;
 
   while ((p = dynamic_cast<genie::GHepParticle *>(event_iter.Next()))) {
     genie::GHepStatus_t ghep_ist = (genie::GHepStatus_t)p->Status();
     int ghep_pdgc = p->Pdg();
     int ghep_fm = p->FirstMother();
     int ghep_fmpdgc = (ghep_fm == -1) ? 0 : ev.Particle(ghep_fm)->Pdg();
 
     // For nuclear targets use hadrons marked as 'hadron in the nucleus'
     // which are the ones passed in the intranuclear rescattering
     // For free nucleon targets use particles marked as 'final state'
     // but make an exception for decayed pi0's,eta's (count them and not their
     // daughters)
 
     bool decayed =
         (ghep_ist == genie::kIStDecayedState &&
          (ghep_pdgc == genie::kPdgPi0 || ghep_pdgc == genie::kPdgEta));
     bool parent_included =
         (ghep_fmpdgc == genie::kPdgPi0 || ghep_fmpdgc == genie::kPdgEta);
 
     bool count_it =
         (nuclear_target && ghep_ist == genie::kIStHadronInTheNucleus) ||
         (!nuclear_target && decayed) ||
         (!nuclear_target && ghep_ist == genie::kIStStableFinalState &&
          !parent_included);
 
     if (!count_it) {
       continue;
     }
 
     if (ghep_pdgc == genie::kPdgPiP) {
       NPip++;
     } else if (ghep_pdgc == genie::kPdgPiM) {
       NPim++;
     } else if (ghep_pdgc == genie::kPdgPi0) {
       NPi0++;
     } else if (ghep_pdgc == genie::kPdgProton) {
       NProton++;
     } else if (ghep_pdgc == genie::kPdgNeutron) {
       NNeutron++;
     } else if (!utility::IsNeutralLepton(
                    ghep_pdgc, utility::pdgcodes::kMatterAntimatter) &&
                !utility::IsChargedLepton(
                    ghep_pdgc, utility::pdgcodes::kMatterAntimatter)) {
       NOther++;
     }
   }
 
   return NFSParticleCount{NProton, NNeutron, NPip, NPi0, NPim, NOther};
 }
 
 Channel_t GetEventChannel(genie::GHepRecord const &gev) {
   // Electron Scattering
   if (gev.Summary()->ProcInfo().IsEM()) {
     if (gev.Summary()->InitState().ProbePdg() == utility::pdgcodes::kElectron) {
       if (gev.Summary()->ProcInfo().IsQuasiElastic()) {
         NFSParticleCount fsparts = CountPreFSIParticles(gev);
         if (fsparts.NProton) {
           return Channel_t::kNCELP;
         } else {
           return Channel_t::kNCELN;
         }
       } else if (gev.Summary()->ProcInfo().IsMEC()) {
         return Channel_t::kNC2p2h;
       } else if (gev.Summary()->ProcInfo().IsResonant()) {
 
         NFSParticleCount fsparts = CountPreFSIParticles(gev);
         if (fsparts.NOther ||
             ((fsparts.NPip + fsparts.NPi0 + fsparts.NPim) > 1)) {
           return Channel_t::kNCTransitionMPi;
         } else if (fsparts.NPip) {
           return Channel_t::kNCSPP_NPip;
         } else if (fsparts.NPi0) {
           return fsparts.NProton ? Channel_t::kNCSPP_PPi0
                                  : Channel_t::kNCSPP_NPi0;
         } else if (fsparts.NPim) {
           return Channel_t::kNCSPP_PPim;
         }
 
         return Channel_t::kNCTransitionMPi;
 
       } else if (gev.Summary()->ProcInfo().IsDeepInelastic()) {
         return Channel_t::kNCDIS;
       } else {
         std::cout << "Unknown GENIE Electron Scattering Mode!" << std::endl
                   << "ScatteringTypeId = "
                   << gev.Summary()->ProcInfo().ScatteringTypeId() << " "
                   << "InteractionTypeId = "
                   << gev.Summary()->ProcInfo().InteractionTypeId() << std::endl
                   << genie::ScatteringType::AsString(
                          gev.Summary()->ProcInfo().ScatteringTypeId())
                   << " "
                   << genie::InteractionType::AsString(
                          gev.Summary()->ProcInfo().InteractionTypeId())
                   << " " << gev.Summary()->ProcInfo().IsMEC() << std::endl;
         return Channel_t::kUndefined;
       }
     }
 
     // Weak CC
   } else if (gev.Summary()->ProcInfo().IsWeakCC()) {
     // CC MEC
     if (gev.Summary()->ProcInfo().IsMEC()) {
       if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(),
                                    utility::pdgcodes::kMatter)) {
         return Channel_t::kCC2p2h;
       } else if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(),
                                           utility::pdgcodes::kAntimatter)) {
         return Channel_t::kCC2p2h_nub;
       }
       // CC OTHER
     } else {
       return FromNEUTCode(genie::utils::ghep::NeutReactionCode(&gev));
     }
 
     // Weak NC
   } else if (gev.Summary()->ProcInfo().IsWeakNC()) {
     // NC MEC
     if (gev.Summary()->ProcInfo().IsMEC()) {
       if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(),
                                    utility::pdgcodes::kMatter)) {
         return Channel_t::kNC2p2h;
       } else if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(),
                                           utility::pdgcodes::kAntimatter)) {
         return Channel_t::kNC2p2h_nub;
       }
 
       // NC OTHER
     } else {
       return FromNEUTCode(genie::utils::ghep::NeutReactionCode(&gev));
     }
   }
 
   return Channel_t::kUndefined;
 }
 
 Particle::Status_t GetParticleStatus(genie::GHepParticle const &p,
                                      Channel_t chan) {
   /*
      kIStUndefined                  = -1,
      kIStInitialState               =  0,   / generator-level initial state /
      kIStStableFinalState           =  1,   / generator-level final state:
      particles to be tracked by detector-level MC /
      kIStIntermediateState          =  2,
      kIStDecayedState               =  3,
      kIStCorrelatedNucleon          = 10,
      kIStNucleonTarget              = 11,
      kIStDISPreFragmHadronicState   = 12,
      kIStPreDecayResonantState      = 13,
      kIStHadronInTheNucleus         = 14,   / hadrons inside the nucleus: marked
      for hadron transport modules to act on /
      kIStFinalStateNuclearRemnant   = 15,   / low energy nuclear fragments
      entering the record collectively as a 'hadronic blob' pseudo-particle /
      kIStNucleonClusterTarget       = 16,   // for composite nucleons before
      phase space decay
      */
 
   Particle::Status_t state = Particle::Status_t::kUnknown;
   switch (p.Status()) {
   case genie::kIStNucleonTarget:
   case genie::kIStInitialState:
   case genie::kIStCorrelatedNucleon:
   case genie::kIStNucleonClusterTarget: {
     state = Particle::Status_t::kPrimaryInitialState;
     break;
   }
 
   case genie::kIStStableFinalState: {
     state = Particle::Status_t::kNuclearLeaving;
     break;
   }
 
   case genie::kIStHadronInTheNucleus: {
     state = utility::Is2p2h(chan) ? Particle::Status_t::kPrimaryInitialState
                                   : Particle::Status_t::kIntermediate;
     break;
   }
 
   case genie::kIStPreDecayResonantState:
   case genie::kIStDISPreFragmHadronicState:
   case genie::kIStIntermediateState: {
     state = Particle::Status_t::kIntermediate;
     break;
   }
 
   case genie::kIStFinalStateNuclearRemnant:
   case genie::kIStUndefined:
   case genie::kIStDecayedState:
   default: { state = Particle::Status_t::kUnknown; }
   }
 
   if (utility::IsNuclearPDG(p.Pdg())) {
     if (state == Particle::Status_t::kPrimaryInitialState) {
       state = Particle::Status_t::kPrimaryInitialState;
     } else if (state == Particle::Status_t::kNuclearLeaving) {
       state = Particle::Status_t::kPrimaryFinalState;
     }
   }
   return state;
 }
 
 struct TargetSplineBlob {
   double MolecularWeight;
   size_t NNucleons;
   std::string search_term;
   TGraph TotSpline;
 };
 
 NEW_NUIS_EXCEPT(Invalid_target_specifier);
 
-double GetFileWeight(fhicl::ParameterSet const &xsecinfo) {
+double GetFileWeight(fhicl::ParameterSet const &xsecinfo,
+                     std::set<std::string> const &spline_list) {
   double mec_scale = xsecinfo.get<double>("mec_scale", 1);
   double EMin = xsecinfo.get<double>("EMin", 0);
   double EMax = xsecinfo.get<double>("EMax", 10);
   size_t NKnots = xsecinfo.get<size_t>("NKnots", 100);
 
   std::vector<TargetSplineBlob> TargetSplines;
   int probe = xsecinfo.get<int>("nu_pdg");
 
   for (std::string const &target :
        xsecinfo.get<std::vector<std::string>>("target")) {
     std::vector<std::string> splits = nuis::utility::split(target, ";");
     if (splits.size() != 1 && splits.size() != 2) {
       throw Invalid_target_specifier()
           << "Expected to find 100ZZZAAA0;<MolecularWeight>, e.g. "
              "1000060120;1, "
              "but found: \""
           << target << "\"";
     }
     std::stringstream ss;
     ss << ".*nu:" << probe << ";tgt:" << splits[0] << ".*";
 
     size_t nuc_pdg = fhicl::string_parsers::str2T<size_t>(splits[0]);
 
     double molwght = splits.size() == 2
                          ? fhicl::string_parsers::str2T<double>(splits[1])
                          : 1;
 
     TargetSplines.push_back(
         {molwght, (nuc_pdg % 1000) / 10, ss.str(), TGraph()});
 
     std::cout << "[INFO]:\tGetting splines that match: \"" << ss.str()
               << "\" with A = " << TargetSplines.back().NNucleons
               << " and multiplied by "
               << molwght / double(TargetSplines.back().NNucleons) << std::endl;
   }
 
   size_t NTotNucleons = 0;
   std::vector<std::string> spline_patterns;
   for (auto const &spb : TargetSplines) {
     spline_patterns.push_back(spb.search_term);
     NTotNucleons += spb.MolecularWeight * spb.NNucleons;
   }
   double WeightToPerNucleon = 1.0 / double(NTotNucleons);
 
   // Read the spline file
   std::cout << "[INFO]: Reading input file: "
             << xsecinfo.get<std::string>("spline_file")
             << ", this may take a while..." << std::endl;
   GENIESplineGetter spg(spline_patterns);
   spg.ParseFile(xsecinfo.get<std::string>("spline_file").c_str());
   std::vector<std::vector<TGraph>> all_splines = spg.GetTGraphs();
 
   double step = (EMax - EMin) / double(NKnots);
 
   // Sum each targets-worth of interaction splines together, weighting for
   // molecular fraction and possibly fixing the MEC spline oddity
   for (size_t ts_it = 0; ts_it < TargetSplines.size(); ++ts_it) {
 
     TargetSplines[ts_it].TotSpline = TGraph(NKnots);
 
-    for (Int_t p_it = 0; p_it < NKnots; ++p_it) {
+    for (size_t p_it = 0; p_it < NKnots; ++p_it) {
       TargetSplines[ts_it].TotSpline.SetPoint(p_it, EMin + p_it * step, 0);
     }
 
     double target_weight = (TargetSplines[ts_it].MolecularWeight /
                             double(TargetSplines[ts_it].NNucleons));
     for (size_t c_it = 0; c_it < all_splines[ts_it].size(); ++c_it) {
       double weight = target_weight;
       if (std::string(all_splines[ts_it][c_it].GetName()).find("MEC") !=
           std::string::npos) {
         if (mec_scale != 1) {
           std::cout << "[INFO]: Weighting MEC splines by " << mec_scale
                     << std::endl;
           weight *= mec_scale;
         }
       }
-      for (Int_t p_it = 0; p_it < NKnots; ++p_it) {
+      for (size_t p_it = 0; p_it < NKnots; ++p_it) {
         double E, XSec;
         TargetSplines[ts_it].TotSpline.GetPoint(p_it, E, XSec);
 
         // Copied from GENIE/Convents/Units.h
         static const double gigaelectronvolt = 1.;
         static const double GeV = gigaelectronvolt;
         static const double meter = 5.07e+15 / GeV;
         static const double centimeter = 0.01 * meter;
         static const double centimeter2 = centimeter * centimeter;
 
         XSec += all_splines[ts_it][c_it].Eval(E) * weight / centimeter2;
 
         TargetSplines[ts_it].TotSpline.SetPoint(p_it, E, XSec);
       }
     }
   }
 
   // Sum all the correctly weighted per-target splines
   TGraph MasterSpline(NKnots);
 
-  for (Int_t p_it = 0; p_it < NKnots; ++p_it) {
+  for (size_t p_it = 0; p_it < NKnots; ++p_it) {
     MasterSpline.SetPoint(p_it, EMin + p_it * step, 0);
   }
   for (size_t c_it = 0; c_it < TargetSplines.size(); ++c_it) {
-    for (Int_t p_it = 0; p_it < NKnots; ++p_it) {
+    for (size_t p_it = 0; p_it < NKnots; ++p_it) {
       double E, XSec;
       MasterSpline.GetPoint(p_it, E, XSec);
       XSec += TargetSplines[c_it].TotSpline.Eval(E) * WeightToPerNucleon;
 
       MasterSpline.SetPoint(p_it, E, XSec);
     }
   }
 
+  nuis::utility::TFile_ptr outfile(nullptr, [](TFile *) {});
+  if (xsecinfo.has_key("spline_output_file")) {
+    outfile = nuis::utility::CheckOpenTFile(
+        xsecinfo.get<std::string>("spline_output_file"), "RECREATE");
+
+    for (auto &sp : TargetSplines) {
+      std::string spname = nuis::utility::SanitizeROOTObjectName(
+          sp.search_term.substr(2, sp.search_term.length() - 4));
+      std::cout << "[INFO]: Dumping spline: " << spname << std::endl;
+      nuis::utility::WriteToTFile(outfile, &sp.TotSpline, spname.c_str());
+    }
+    nuis::utility::WriteToTFile(outfile, &MasterSpline, "TotalXSec");
+  }
+
   fhicl::ParameterSet fluxps = xsecinfo.get<fhicl::ParameterSet>("flux");
   // If MonoE
   if (fluxps.has_key("energy_GeV")) {
     double monoE = fluxps.get<double>("energy_GeV");
 
     return MasterSpline.Eval(monoE);
   } else {
     // Use the master spline to get flux*xsec (which the events are thrown
     // according to) and return the flux-averaged total xsec which is used for
     // weighting
     std::unique_ptr<TH1> Flux =
         nuis::utility::GetHistogram<TH1>(fluxps.get<std::string>("histogram"));
+
+    if (outfile) {
+      std::unique_ptr<TH1> Fluxcp = nuis::utility::Clone<TH1>(Flux);
+      nuis::utility::WriteToTFile(outfile, Fluxcp.get(), "Flux");
+    }
+
     bool per_width = fluxps.get<bool>("is_divided_by_bin_width", true);
     double FluxIntegral = Flux->Integral(per_width ? "width" : "");
 
     for (int bi_it = 0; bi_it < Flux->GetXaxis()->GetNbins(); ++bi_it) {
       double bc = Flux->GetBinContent(bi_it + 1);
       double low_e = Flux->GetXaxis()->GetBinLowEdge(bi_it + 1);
       double up_e = Flux->GetXaxis()->GetBinUpEdge(bi_it + 1);
 
       size_t NIntSteps = 100;
       double step = (up_e - low_e) / double(NIntSteps);
 
       double avg_xsec = 0;
       for (size_t i = 0; i < NIntSteps; ++i) {
         double e = low_e + (double(i) + 0.5) * step;
         avg_xsec += MasterSpline.Eval(e);
       }
       avg_xsec /= double(NIntSteps);
 
       Flux->SetBinContent(bi_it + 1, bc * avg_xsec);
     }
 
+    if (outfile) {
+      std::unique_ptr<TH1> Fluxcp = nuis::utility::Clone<TH1>(Flux);
+      nuis::utility::WriteToTFile(outfile, Fluxcp.get(), "EvRate");
+    }
+
     double EvRateIntegral = Flux->Integral(per_width ? "width" : "");
 
     return EvRateIntegral / FluxIntegral;
   }
 }
 
 } // namespace genietools
 } // namespace nuis
diff --git a/src/generator/utility/GENIEUtility.hxx b/src/generator/utility/GENIEUtility.hxx
index 3bc447a..bb25aef 100644
--- a/src/generator/utility/GENIEUtility.hxx
+++ b/src/generator/utility/GENIEUtility.hxx
@@ -1,53 +1,53 @@
 // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is part of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 
 #pragma once
 
 #include "event/Particle.hxx"
 
 #include "event/types.hxx"
 
+#include <set>
+
 class TGraph;
 
 namespace genie {
 class GHepRecord;
 class GHepParticle;
 } // namespace genie
 
 namespace fhicl {
   class ParameterSet;
 }
 
 namespace nuis {
 namespace genietools {
 
 NEW_NUIS_EXCEPT(invalid_GENIE_event);
 
-TGraph const &GetGENIESpline(std::string const &SplineFile,
-                             std::string const &SplineIdentifier);
-
 event::Channel_t GetEventChannel(genie::GHepRecord const &);
 
 event::Particle::Status_t GetParticleStatus(genie::GHepParticle const &p,
                                             event::Channel_t chan);
 
-double GetFileWeight(fhicl::ParameterSet const &xsecinfo);
+double GetFileWeight(fhicl::ParameterSet const &xsecinfo,
+                     std::set<std::string> const &spline_list = {});
 
 } // namespace genietools
 } // namespace nuis
diff --git a/src/input/IInputHandler.hxx b/src/input/IInputHandler.hxx
index fcb467e..205ea65 100644
--- a/src/input/IInputHandler.hxx
+++ b/src/input/IInputHandler.hxx
@@ -1,122 +1,123 @@
 // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is part of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 
 #pragma once
 
 #include "plugins/traits.hxx"
 
 #include "exception/exception.hxx"
 
 #include "generator/GeneratorManager.hxx"
 
 #include <iterator>
 #include <limits>
 
 namespace fhicl {
 class ParameterSet;
 }
 
 namespace nuis {
 namespace event {
 class MinimalEvent;
 class FullEvent;
 } // namespace event
 } // namespace nuis
 
 class IInputHandler {
 public:
   struct FullEvent_const_iterator
       : public std::iterator<
             std::input_iterator_tag, nuis::event::FullEvent const, size_t,
             nuis::event::FullEvent const *, nuis::event::FullEvent const &> {
 
     FullEvent_const_iterator(size_t _idx, IInputHandler const *_ih) {
       idx = _idx;
       ih = _ih;
     }
     FullEvent_const_iterator(FullEvent_const_iterator const &other) {
       idx = other.idx;
       ih = other.ih;
     }
     FullEvent_const_iterator operator=(FullEvent_const_iterator const &other) {
       idx = other.idx;
       ih = other.ih;
       return (*this);
     }
 
     bool operator==(FullEvent_const_iterator const &other) {
       return (idx == other.idx);
     }
     bool operator!=(FullEvent_const_iterator const &other) {
       return !(*this == other);
     }
     nuis::event::FullEvent const &operator*() { return ih->GetFullEvent(idx); }
     nuis::event::FullEvent const *operator->() {
       return &ih->GetFullEvent(idx);
     }
 
     FullEvent_const_iterator operator++() {
       idx++;
       return *this;
     }
     FullEvent_const_iterator operator++(int) {
       FullEvent_const_iterator tmp(*this);
       idx++;
       return tmp;
     }
 
   private:
     size_t idx;
     IInputHandler const *ih;
   };
 
   NEW_NUIS_EXCEPT(invalid_input_file);
   NEW_NUIS_EXCEPT(invalid_entry);
   NEW_NUIS_EXCEPT(input_handler_feature_unimplemented);
 
   typedef size_t ev_index_t;
 
   virtual void Initialize(fhicl::ParameterSet const &) = 0;
+  virtual void GetEntry(ev_index_t) const = 0;
   virtual nuis::event::MinimalEvent const &
   GetMinimalEvent(ev_index_t idx) const = 0;
   virtual nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const = 0;
   virtual void RecalculateEventWeights(){};
   virtual double GetEventWeight(ev_index_t) const { return 1; };
 
   /// Allows samples that implement flux cuts to request an appropriate
   virtual double GetXSecScaleFactor(
       std::pair<double, double> const &EnuRange = std::pair<double, double>{
           std::numeric_limits<double>::max(),
           std::numeric_limits<double>::max()}) const = 0;
 
   virtual size_t GetNEvents() const = 0;
 
   FullEvent_const_iterator begin() const {
     return FullEvent_const_iterator(0, this);
   }
   FullEvent_const_iterator end() const {
     return FullEvent_const_iterator(GetNEvents(), this);
   }
 
   virtual ~IInputHandler() {}
 
   virtual nuis::GeneratorManager::Generator_id_t GetGeneratorId() const = 0;
 };
 
 DECLARE_PLUGIN_INTERFACE(IInputHandler);
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
new file mode 100644
index 0000000..c0f8bcb
--- /dev/null
+++ b/src/tests/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(memcheck)
diff --git a/src/tests/memcheck/CMakeLists.txt b/src/tests/memcheck/CMakeLists.txt
new file mode 100644
index 0000000..c3fb4e8
--- /dev/null
+++ b/src/tests/memcheck/CMakeLists.txt
@@ -0,0 +1,30 @@
+SET(TESTS InputHandlerCheck)
+
+foreach(a ${TESTS})
+  add_executable(${a} ${a}.cxx)
+  target_link_libraries(${a} -Wl,--no-as-needed nuis_event nuis_input nuis_utility_experimental nuis_utility nuis_config nuis_persistency nuis_plugins nuis_params nuis_variation InputHandlers)
+  target_link_libraries(${a} ${GENERATOR_DEPENDENT_TARGETS} -Wl,--as-needed)
+  target_link_libraries(${a} ${NUISANCE_LINK_DIRS})
+  target_link_libraries(${a} ${NUISANCE_DEPEND_LIBS})
+
+  if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
+    set_target_properties(${a} PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
+  endif()
+
+    if(USE_GENIE)
+      target_compile_options(${a} PUBLIC ${GENIE_CXX_FLAGS})
+      include_directories(${GENIE_INCLUDE_DIRS})
+    endif()
+
+    if(USE_NEUT)
+      target_compile_options(${a} PUBLIC ${NEUT_CXX_FLAGS})
+      include_directories(${NEUT_INCLUDE_DIRS})
+    endif()
+
+    if(USE_NUWRO)
+      target_compile_options(${a} PUBLIC ${NUWRO_CXX_FLAGS})
+      include_directories(${NUWRO_INCLUDE_DIRS})
+    endif()
+
+  install(TARGETS ${a} DESTINATION test)
+endforeach()
diff --git a/src/tests/memcheck/InputHandlerCheck.cxx b/src/tests/memcheck/InputHandlerCheck.cxx
new file mode 100644
index 0000000..9a795ad
--- /dev/null
+++ b/src/tests/memcheck/InputHandlerCheck.cxx
@@ -0,0 +1,156 @@
+#include "config/GlobalConfiguration.hxx"
+
+#ifdef USE_GENIE
+#include "generator/input/GENIEInputHandler.hxx"
+#endif
+#ifdef USE_NEUT
+#include "generator/input/NEUTInputHandler.hxx"
+#endif
+#include "input/InputManager.hxx"
+
+#include "event/MinimalEvent.hxx"
+
+#include "exception/exception.hxx"
+
+#include "utility/StringUtility.hxx"
+
+#include "fhiclcpp/make_ParameterSet.h"
+#include "string_parsers/from_string.hxx"
+
+#include <string>
+
+NEW_NUIS_EXCEPT(invalid_cli_arguments);
+
+void SayUsage(char const *argv[]) {
+  std::cout
+      << "[USAGE]: " << argv[0]
+      << "\n"
+         "\t-c <config.fcl>             : FHiCL file containing study "
+         "configuration. \n"
+         "\t-s <sample name>            : FHiCL key of a single sample "
+         "to run from the -c argument. \n"
+         "\t-n <NTries>                 : Number of events toload. \n"
+#ifdef USE_GENIE
+         "\t--genie                     : Just use GENIEInputHandler directly "
+#endif
+#ifdef USE_NEUT
+         "\t--neut                      : Just use NEUTInputHandler directly "
+#endif
+         "\t-G                          : Just use GetEntry directly "
+         "rather than the GetMinimalEvent interface."
+      << std::endl;
+}
+
+std::string fhicl_file = "";
+std::string named_sample = "";
+size_t NTries = 0;
+bool GetEntry = false;
+
+enum EGen {
+  kManager,
+#ifdef USE_NEUT
+  kNEUT,
+#endif
+#ifdef USE_GENIE
+  kGENIE,
+#endif
+};
+EGen WhichGen = kManager;
+
+void handleOpts(int argc, char const *argv[]) {
+  int opt = 1;
+  while (opt < argc) {
+    if ((std::string(argv[opt]) == "-?") ||
+        (std::string(argv[opt]) == "--help")) {
+      SayUsage(argv);
+      exit(0);
+    } else if (std::string(argv[opt]) == "-c") {
+      fhicl_file = argv[++opt];
+    } else if (std::string(argv[opt]) == "-s") {
+      named_sample = argv[++opt];
+    } else if (std::string(argv[opt]) == "-n") {
+      NTries = fhicl::string_parsers::str2T<size_t>(argv[++opt]);
+    } else if (std::string(argv[opt]) == "-G") {
+      GetEntry = true;
+    }
+#ifdef USE_GENIE
+    else if (std::string(argv[opt]) == "--genie") {
+      WhichGen = kGENIE;
+    }
+#endif
+#ifdef USE_NEUT
+    else if (std::string(argv[opt]) == "--neut") {
+      WhichGen = kNEUT;
+    }
+#endif
+    else {
+      std::cout << "[ERROR]: Unknown option: " << argv[opt] << std::endl;
+      SayUsage(argv);
+      exit(1);
+    }
+    opt++;
+  }
+}
+
+int main(int argc, char const *argv[]) {
+  nuis::config::EnsureConfigurationRead("nuis.global.config.fcl");
+  nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl");
+
+  handleOpts(argc, argv);
+
+  if (!fhicl_file.size()) {
+    SayUsage(argv);
+    throw invalid_cli_arguments();
+  }
+
+  nuis::config::EnsureConfigurationRead(fhicl_file);
+
+  fhicl::ParameterSet inps =
+      nuis::config::GetDocument().get<fhicl::ParameterSet>(named_sample);
+  IInputHandler const *IH;
+
+  switch (WhichGen) {
+  case kManager: {
+    auto id = nuis::input::InputManager::Get().EnsureInputLoaded(inps);
+    IH = &nuis::input::InputManager::Get().GetInputHandler(id);
+    break;
+  }
+#ifdef USE_NEUT
+  case kNEUT: {
+    IInputHandler *IIH = new NEUTInputHandler();
+    IIH->Initialize(inps);
+    IH = IIH;
+    break;
+  }
+#endif
+#ifdef USE_GENIE
+  case kGENIE: {
+    IInputHandler *IIH = new GENIEInputHandler();
+    IIH->Initialize(inps);
+    IH = IIH;
+    break;
+  }
+#endif
+  default: {
+    std::cout << "[ERROR]: Invalid input handler type." << std::endl;
+    return 1;
+  }
+  }
+
+  size_t ShoutEvery = NTries / 100;
+
+  std::cout << "[INFO]: Read " << 0 << "/" << NTries << " GENIE events."
+            << std::flush;
+  for (size_t t_it = 0; t_it < NTries; ++t_it) {
+    if (ShoutEvery && !(t_it % ShoutEvery)) {
+      std::cout << "\r[INFO]: Read " << t_it << "/" << NTries
+                << " GENIE events." << std::flush;
+    }
+    if (GetEntry) {
+      IH->GetEntry(0);
+    } else {
+      IH->GetMinimalEvent(0);
+    }
+  }
+  std::cout << std::endl;
+}
diff --git a/src/utility/ROOTUtility.hxx b/src/utility/ROOTUtility.hxx
index 75b0d47..2099789 100644
--- a/src/utility/ROOTUtility.hxx
+++ b/src/utility/ROOTUtility.hxx
@@ -1,164 +1,181 @@
 // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is part of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 
 #ifndef UTILITY_ROOTUTILITY_HXX_SEEN
 #define UTILITY_ROOTUTILITY_HXX_SEEN
 
 #include "exception/exception.hxx"
 
 #include "TFile.h"
 #include "TTree.h"
 
 #include <iomanip>
 #include <iostream>
 #include <memory>
 #include <string>
 
 namespace nuis {
 namespace utility {
 
 NEW_NUIS_EXCEPT(failed_to_open_TFile);
 NEW_NUIS_EXCEPT(failed_to_get_TTree);
+NEW_NUIS_EXCEPT(WriteToTFile_nullptr);
 
 inline void CloseTFile(TFile *f = nullptr) {
   if (f) {
     std::cout << "[INFO]: Shutting TFile: " << f->GetName() << ", " << f
               << std::endl;
     f->Close();
   }
   delete f;
 }
 
 using TFile_ptr = std::unique_ptr<TFile, decltype(&CloseTFile)>;
 inline TFile_ptr make_unique_TFile(TFile *f) {
   return TFile_ptr(f, &CloseTFile);
 }
 
 inline TFile_ptr CheckOpenTFile(std::string const &fname,
                                 char const *opts = "") {
 
   TDirectory *ogDir = gDirectory;
 
   TFile *inpF = new TFile(fname.c_str(), opts);
   if (!inpF || !inpF->IsOpen()) {
     throw failed_to_open_TFile()
         << "[ERROR]: Couldn't open input file: " << std::quoted(fname);
   }
 
   if (ogDir) {
     ogDir->cd();
   }
 
   return make_unique_TFile(inpF);
 }
 
 struct TreeFile {
   TFile_ptr file;
   TTree *tree;
   bool file_owned;
   TreeFile()
       : file(make_unique_TFile(nullptr)), tree(nullptr), file_owned(false) {}
   TreeFile(TreeFile const &other) = delete;
   TreeFile &operator=(TreeFile &&other) {
     file = std::move(other.file);
     tree = other.tree;
     file_owned = other.file_owned;
 
     other.file = nullptr;
     other.tree = nullptr;
     other.file_owned = false;
 
     return *this;
   }
 
   TreeFile(TreeFile &&other)
       : file(std::move(other.file)), tree(other.tree),
         file_owned(other.file_owned) {
     other.file = nullptr;
     other.tree = nullptr;
     other.file_owned = false;
   }
 
   TTree *operator->() { return tree; }
 
   ~TreeFile() {
     if (!file_owned) {
       file.release();
     }
   }
 };
 
 inline TreeFile CheckGetTTree(TFile *file, std::string const &tname) {
   TreeFile tf;
   tf.file = make_unique_TFile(file);
   tf.tree = dynamic_cast<TTree *>(tf.file->Get(tname.c_str()));
   tf.file_owned = false;
 
   if (!tf.tree) {
     throw failed_to_get_TTree()
         << "[ERROR]: Failed to get TTree named: " << std::quoted(tname)
         << " from TFile: " << std::quoted(file->GetName());
   }
   std::cout << "[INFO]: Opened TFile: " << file->GetName() << ", " << file
             << std::endl;
   return tf;
 }
 
 inline TreeFile CheckGetTTree(std::string const &fname,
                               std::string const &tname, char const *opts = "") {
   TreeFile tf = CheckGetTTree(CheckOpenTFile(fname, opts).release(), tname);
   tf.file_owned = true;
   return tf;
 }
 
 inline TreeFile MakeNewTTree(std::string const &fname, std::string const &tname,
                              char const *opts = "") {
   TreeFile tf;
   tf.file = CheckOpenTFile(fname, opts);
   tf.tree = new TTree(tname.c_str(), "");
   tf.tree->SetDirectory(tf.file.get());
   tf.file_owned = true;
   return tf;
 }
 
 inline TreeFile AddNewTTreeToFile(TFile *f, std::string const &tname,
                                   std::string const &dirname = "") {
   TreeFile tf;
   // This is pretty dumb...
   tf.file = make_unique_TFile(f);
   tf.tree = new TTree(tname.c_str(), "");
   if (dirname.size()) {
     TDirectory *dir = tf.file->GetDirectory(dirname.c_str());
     if (!dir) {
       dir = tf.file->mkdir(dirname.c_str());
     }
     tf.tree->SetDirectory(dir);
   } else {
     tf.tree->SetDirectory(tf.file.get());
   }
   tf.file_owned = false;
   return tf;
 }
 
 inline std::string SanitizeROOTObjectName(std::string name) { return name; }
 
+inline void WriteToTFile(TFile_ptr &tf, TObject *object,
+                         std::string const &object_name) {
+  if (!object) {
+    throw WriteToTFile_nullptr();
+  }
+
+  TDirectory *ogdir = gDirectory;
+
+  tf->cd();
+  object->Write(object_name.c_str(), TObject::kOverwrite);
+
+  if (ogdir) {
+    ogdir->cd();
+  }
+}
+
 } // namespace utility
 } // namespace nuis
 
 #endif