No OneTemporary
Actions

Size

70 KB

Subscribers

None

View Options

	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

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:28 PM (1 d, 16 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3804950
Default Alt Text: (70 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions