diff --git a/scripts/to_configure/nuisbuild.in b/scripts/to_configure/nuisbuild.in
index 316b790..1daaeb1 100644
--- a/scripts/to_configure/nuisbuild.in
+++ b/scripts/to_configure/nuisbuild.in
@@ -1,101 +1,109 @@
 # 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/>.
 ################################################################################
 
 #!/bin/bash
 
 EXTRA_COMPILE_FLAGS=""
 EXTRA_LINK_FLAGS=""
 INPUT_FILE_NAME=""
 OUTPUT_LIB_NAME=""
+DEBUG_FLAGS=""
 
 SCRIPTNAME=${0}
 
 while [[ ${#} -gt 0 ]]; do
 
   key="$1"
   case $key in
 
       -i|--implementation-file)
 
       if [[ ${#} -lt 2 ]]; then
         echo "[ERROR]: ${1} expected a value."
         exit 1
       fi
 
       INPUT_FILE_NAME="$2"
       echo "[OPT]: Reading event processor implementation from: ${INPUT_FILE_NAME}"
       shift # past argument
       ;;
 
       --compile-flags)
 
       if [[ ${#} -lt 2 ]]; then
         echo "[ERROR]: ${1} expected a value."
         exit 1
       fi
 
       EXTRA_COMPILE_FLAGS="$2"
       echo "[OPT]: Including extra compile flags: ${EXTRA_COMPILE_FLAGS}"
       shift # past argument
       ;;
 
+      -g|--debug)
+
+      DEBUG_FLAGS=" -g -O0 "
+      echo "[OPT]: Building with debugging symbols."
+      ;;
+
       --link-flags)
 
       if [[ ${#} -lt 2 ]]; then
         echo "[ERROR]: ${1} expected a value."
         exit 1
       fi
 
       EXTRA_LINK_FLAGS="$2"
       echo "[OPT]: Including extra linker flags: ${EXTRA_LINK_FLAGS}"
       shift # past argument
       ;;
 
       -?|--help)
 
       echo "[RUNLIKE] ${SCRIPTNAME}"
       echo -e "\t-i|--implementation-file   : Path to IEventProcessor subclass to build into dynamically loadable shared object."
       echo -e "\t--compile-flags            : Extra compile flags to pass to the compiler."
       echo -e "\t--link-flags               : Extra linker flags to pass to the linker."
+      echo -e "\t-g|--debug                 : Add default debugging compiler flags.."
       echo -e "\t-?|--help                  : Print this message."
       exit 0
       ;;
 
       *)
               # unknown option
       echo "Unknown option $1"
       exit 1
       ;;
 
   esac
   shift
 done
 
 OUTPUT_LIB_NAME=${INPUT_FILE_NAME%.*}.so
 
 #Removes the NUISANCE filename descriptor that uses a gmake trick to insert the file name into the compiler options
 NUISANCE_CXX_FLAGS_GMAKE_SCRUB=$(echo '@NUISANCE_CXX_FLAGS@' | sed s/-D__FILENAME__=\".*\"//g)
 
-echo "@CMAKE_CXX_COMPILER@ ${INPUT_FILE_NAME} -shared -fPIC -o ${OUTPUT_LIB_NAME} ${NUISANCE_CXX_FLAGS_GMAKE_SCRUB} ${EXTRA_COMPILE_FLAGS} -D__FILENAME__=\"${INPUT_FILE_NAME##*/}\" -I. -I@CMAKE_INSTALL_PREFIX@/include  @CMAKE_LINK_FLAGS@ ${EXTRA_LINK_FLAGS} -L@CMAKE_INSTALL_PREFIX@/lib @NUISANCE_LINK_DIRS@  @NUISANCE_DEPEND_LIBS@"
+echo "@CMAKE_CXX_COMPILER@ ${INPUT_FILE_NAME} -shared -fPIC -o ${OUTPUT_LIB_NAME} ${NUISANCE_CXX_FLAGS_GMAKE_SCRUB} ${EXTRA_COMPILE_FLAGS} -D__FILENAME__=\"${INPUT_FILE_NAME##*/}\" -I. -I@CMAKE_INSTALL_PREFIX@/include ${DEBUG_FLAGS} @CMAKE_LINK_FLAGS@ ${EXTRA_LINK_FLAGS} -L@CMAKE_INSTALL_PREFIX@/lib @NUISANCE_LINK_DIRS@  @NUISANCE_DEPEND_LIBS@"
 
-if ! @CMAKE_CXX_COMPILER@ ${INPUT_FILE_NAME} -shared -fPIC -o ${OUTPUT_LIB_NAME} ${NUISANCE_CXX_FLAGS_GMAKE_SCRUB} ${EXTRA_COMPILE_FLAGS} -D__FILENAME__=\"${INPUT_FILE_NAME##*/}\" -I. -I@CMAKE_INSTALL_PREFIX@/include  @CMAKE_LINK_FLAGS@ ${EXTRA_LINK_FLAGS} -L@CMAKE_INSTALL_PREFIX@/lib @NUISANCE_LINK_DIRS@  @NUISANCE_DEPEND_LIBS@; then
+if ! @CMAKE_CXX_COMPILER@ ${INPUT_FILE_NAME} -shared -fPIC -o ${OUTPUT_LIB_NAME} ${NUISANCE_CXX_FLAGS_GMAKE_SCRUB} ${EXTRA_COMPILE_FLAGS} -D__FILENAME__=\"${INPUT_FILE_NAME##*/}\" -I. -I@CMAKE_INSTALL_PREFIX@/include ${DEBUG_FLAGS} @CMAKE_LINK_FLAGS@ ${EXTRA_LINK_FLAGS} -L@CMAKE_INSTALL_PREFIX@/lib @NUISANCE_LINK_DIRS@  @NUISANCE_DEPEND_LIBS@; then
   echo "[ERROR]: Failed to compile ${INPUT_FILE_NAME} -> ${OUTPUT_LIB_NAME}"
 else
   echo "[INFO]: Successfully build: ${OUTPUT_LIB_NAME}."
 fi
diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx
index dda5a6d..dbd85bd 100644
--- a/src/generator/input/NEUTInputHandler.cxx
+++ b/src/generator/input/NEUTInputHandler.cxx
@@ -1,129 +1,198 @@
 #include "generator/input/NEUTInputHandler.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);
 
-  fFlux = GetHistogramFromROOTFile<TH1>(fInputTreeFile.file, "flux_numu");
-  fEvtrt = GetHistogramFromROOTFile<TH1>(fInputTreeFile.file, "evtrt_numu");
+  fFlux =
+      GetHistogramFromROOTFile<TH1>(fInputTreeFile.file, "flux_numu", false);
+  if (fFlux) {
+    fEvtrt = GetHistogramFromROOTFile<TH1>(fInputTreeFile.file, "evtrt_numu");
 
-  fFileWeight = fEvtrt->Integral() / (fFlux->Integral() * double(GetNEvents()));
+    fFileWeight =
+        fEvtrt->Integral() / (fFlux->Integral() * double(GetNEvents()));
+  } 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;
+  }
+}
+
+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;
+    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;
 }
 
 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 {
   if (idx >= GetNEvents()) {
     throw IInputHandler::invalid_entry()
         << "[ERROR]: Attempted to get entry " << idx
         << " from an InputHandler with only " << GetNEvents();
   }
   fInputTreeFile.tree->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);
   }
 
   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(){
+GeneratorManager::Generator_id_t NEUTInputHandler::GetGeneratorId() {
   return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT");
 }
 
 DECLARE_PLUGIN(IInputHandler, NEUTInputHandler);
diff --git a/src/generator/input/NEUTInputHandler.hxx b/src/generator/input/NEUTInputHandler.hxx
index 55cbf8c..7fc1c3f 100644
--- a/src/generator/input/NEUTInputHandler.hxx
+++ b/src/generator/input/NEUTInputHandler.hxx
@@ -1,61 +1,69 @@
 // 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;
 
   std::unique_ptr<TH1> fFlux;
   std::unique_ptr<TH1> fEvtrt;
   double fFileWeight;
 
   bool fKeepIntermediates;
 
+  double GetMonoEXSecWeight();
+
 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;
   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();
 };
diff --git a/src/generator/input/NuWroInputHandler.cxx b/src/generator/input/NuWroInputHandler.cxx
index 7769c14..0b50cb7 100644
--- a/src/generator/input/NuWroInputHandler.cxx
+++ b/src/generator/input/NuWroInputHandler.cxx
@@ -1,127 +1,132 @@
 #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()){
+    fInputTree.tree->GetEntry(0);
+    std::cout << "[INFO]: Average NuWro XSec weight: " << (fTreeEvent->weight / double(GetNEvents())) << std::endl;
+  }
 }
 MinimalEvent const &NuWroInputHandler::GetMinimalEvent(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);
 
   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);
   }
 
   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(){
   return GeneratorManager::Get().EnsureGeneratorRegistered("NuWro");
 }
 
 DECLARE_PLUGIN(IInputHandler, NuWroInputHandler);
diff --git a/src/generator/utility/NEUTUtility.cxx b/src/generator/utility/NEUTUtility.cxx
index 56cd4a9..b7f6310 100644
--- a/src/generator/utility/NEUTUtility.cxx
+++ b/src/generator/utility/NEUTUtility.cxx
@@ -1,90 +1,93 @@
 #include "generator/utility/NEUTUtility.hxx"
 
+#include "generator/input/NEUTInputHandler.hxx"
+
 #include "exception/exception.hxx"
 
 #include "utility/InteractionChannelUtility.hxx"
 #include "utility/PDGCodeUtility.hxx"
 
 #include "neutpart.h"
+#include "neutvect.h"
 
 using namespace nuis::event;
 using namespace nuis::utility;
 
 namespace nuis {
 namespace neuttools {
 NEW_NUIS_EXCEPT(unexpected_NEUT_particle_state);
 
 Particle::Status_t GetNeutParticleStatus(NeutPart const &part, Channel_t mode) {
 
 #ifdef DEBUG_NEUT_UTILITY
   std::cout << "[DEBUG]: Mode: " << mode << ", Part: { Status: " << part.fStatus
             << ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << " }."
             << std::endl;
 #endif
 
   // Remove Pauli blocked events, probably just single pion events
   if (part.fStatus == 5) {
     return Particle::Status_t::kBlocked;
 
     // fStatus == -1 means initial  state
   } else if (part.fIsAlive == false && part.fStatus == -1) {
     return Particle::Status_t::kPrimaryInitialState;
 
     // NEUT has a bit of a strange convention for fIsAlive and fStatus
     // combinations
     // for NC and neutrino particle isAlive true/false and status 2 means
     // final state particle
     // for other particles in NC status 2 means it's an FSI particle
     // for CC it means it was an FSI particle
   } else if (part.fStatus == 2) {
     // NC case is a little strange... The outgoing neutrino might be alive or
     // not alive. Remaining particles with status 2 are FSI particles that
     // reinteracted
     if (IsNC(mode) && IsNeutralLepton(part.fPID)) {
       return Particle::Status_t::kNuclearLeaving;
       // The usual CC case
     } else if (part.fIsAlive == true) {
       // return Particle::Status_t::kIntermediate;
       throw unexpected_NEUT_particle_state()
           << "[ERROR] Found unexpected NEUT particle in neutvect stack: Mode: "
           << mode << ", Part: { Status: " << part.fStatus
           << ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << " }.";
     }
 
   } else if ((part.fIsAlive == true) && (part.fStatus == 2) &&
              IsNeutralLepton(part.fPID)) {
     return Particle::Status_t::kNuclearLeaving;
 
   } else if ((part.fIsAlive == true) && (part.fStatus == 0)) {
     return Particle::Status_t::kNuclearLeaving;
   } else if (!part.fIsAlive && ((part.fStatus == 1) || (part.fStatus == 3) ||
                                 (part.fStatus == 4) || (part.fStatus == 7) ||
                                 (part.fStatus == 8))) {
     return Particle::Status_t::kIntermediate;
 
     // There's one hyper weird case where fStatus = -3. This apparently
     // corresponds to a nucleon being ejected via pion FSI when there is "data
     // available"
   } else if (!part.fIsAlive && (part.fStatus == -3)) {
     return Particle::Status_t::kUnknown;
     // NC neutrino outgoing
   } else if (!part.fIsAlive && part.fStatus == 0 &&
              IsNeutralLepton(part.fPID)) {
     return Particle::Status_t::kNuclearLeaving;
 
     // Warn if we still find alive particles without classifying them
   } else if (part.fIsAlive == true) {
     throw unexpected_NEUT_particle_state()
         << "[ERROR]: Undefined NEUT state "
         << " Alive: " << part.fIsAlive << " Status: " << part.fStatus
         << " PDG: " << part.fPID;
   }
   // Warn if we find dead particles that we haven't classified
   throw unexpected_NEUT_particle_state()
       << "[ERROR]: Undefined NEUT state "
       << " Alive: " << part.fIsAlive << " Status: " << part.fStatus
       << " PDG: " << part.fPID;
 }
 
 } // namespace neuttools
 } // namespace nuis
diff --git a/src/generator/utility/NEUTUtility.hxx b/src/generator/utility/NEUTUtility.hxx
index c191e52..fdcf697 100644
--- a/src/generator/utility/NEUTUtility.hxx
+++ b/src/generator/utility/NEUTUtility.hxx
@@ -1,15 +1,17 @@
-#ifndef GENERATOR_UTILITY_NEUTUTILITY_HXX_SEEN
-#define GENERATOR_UTILITY_NEUTUTILITY_HXX_SEEN
+#pragma once
 
 #include "event/Particle.hxx"
 #include "event/types.hxx"
 
+#include "exception/exception.hxx"
+
+#include <memory>
+
 class NeutPart;
 
 namespace nuis {
 namespace neuttools {
 nuis::event::Particle::Status_t GetNeutParticleStatus(NeutPart const &,
                                                       nuis::event::Channel_t);
 }
 } // namespace nuis
-#endif
diff --git a/src/samples/SimpleMCStudy.hxx b/src/samples/SimpleMCStudy.hxx
index c8a55d5..b6da7ae 100644
--- a/src/samples/SimpleMCStudy.hxx
+++ b/src/samples/SimpleMCStudy.hxx
@@ -1,108 +1,111 @@
 // 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 "samples/IEventProcessor.hxx"
 
 #include "event/FullEvent.hxx"
 
 #include "input/InputManager.hxx"
 
 #include <functional>
 
 class SimpleMCStudy : public IEventProcessor {
 protected:
   nuis::input::InputManager::Input_id_t fIH_id;
   std::string write_directory;
   size_t NMaxSample_override;
 
   std::function<void(nuis::event::FullEvent const &, double)>
       ProcessEventFunction;
 
 public:
   SimpleMCStudy() {
     fIH_id = std::numeric_limits<nuis::input::InputManager::Input_id_t>::max();
     write_directory = "";
     NMaxSample_override = std::numeric_limits<size_t>::max();
 
-    ProcessEventFunction = std::function<void(nuis::event::FullEvent const &, double)>();
+    ProcessEventFunction =
+        std::function<void(nuis::event::FullEvent const &, double)>();
   }
 
   fhicl::ParameterSet fGlobalConfig;
   fhicl::ParameterSet fInstanceConfig;
 
   void ReadGlobalConfigDefaults() {
     fGlobalConfig = nuis::config::GetDocument().get<fhicl::ParameterSet>(
         std::string("global.sample_configuration.") + Name(),
         fhicl::ParameterSet());
   }
 
   void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
 
     fInstanceConfig = instance_sample_configuration;
 
-    if (fInstanceConfig.has_key("verbosity")) {
-      SetSampleVerbosity(fInstanceConfig.get<std::string>("verbosity"));
-    } else {
-      SetSampleVerbosity("Reticent");
-    }
+    SetSampleVerbosity(fInstanceConfig.get<std::string>(
+        "verbosity", nuis::config::GetDocument().get<std::string>(
+                         "global.sample.verbosity_default", "Reticent")));
 
     ReadGlobalConfigDefaults();
 
     NMaxSample_override =
         fInstanceConfig.get<size_t>("nmax", std::numeric_limits<size_t>::max());
 
     write_directory =
         fInstanceConfig.get<std::string>("write_directory", Name());
 
     fIH_id =
         nuis::input::InputManager::Get().EnsureInputLoaded(fInstanceConfig);
   }
 
   void ProcessSample(size_t nmax = std::numeric_limits<size_t>::max()) {
     if (fIH_id ==
         std::numeric_limits<nuis::input::InputManager::Input_id_t>::max()) {
       throw uninitialized_IEventProcessor();
     }
     IInputHandler const &IH =
         nuis::input::InputManager::Get().GetInputHandler(fIH_id);
 
     nmax = std::min(NMaxSample_override, nmax);
 
     size_t NEvsToProcess = std::min(nmax, IH.GetNEvents());
 
     double nmax_scaling = double(IH.GetNEvents()) / double(NEvsToProcess);
 
     size_t NToShout = NEvsToProcess / 10;
     IEventProcessor_INFO("Sample " << std::quoted(Name()) << " processing "
                                    << NEvsToProcess << " events.");
 
     IInputHandler::ev_index_t ev_idx = 0;
 
     while (ev_idx < NEvsToProcess) {
+      if (NToShout && !(ev_idx % NToShout)) {
+        IEventProcessor_INFO("\t\t Done " << ev_idx << "/" << NEvsToProcess
+                                          << " events.");
+      }
       nuis::event::FullEvent const &fev = IH.GetFullEvent(ev_idx);
       ProcessEventFunction(fev, IH.GetEventWeight(ev_idx) * nmax_scaling);
       ev_idx++;
     }
 
     IEventProcessor_INFO("Sample " << std::quoted(Name()) << " processed "
                                    << NEvsToProcess << " events.");
   }
 };
diff --git a/src/utility/HistogramUtility.hxx b/src/utility/HistogramUtility.hxx
index b3f8ef4..851d84c 100644
--- a/src/utility/HistogramUtility.hxx
+++ b/src/utility/HistogramUtility.hxx
@@ -1,519 +1,522 @@
 // 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_HISTOGRAMUTILITY_HXX_SEEN
 #define UTILITY_HISTOGRAMUTILITY_HXX_SEEN
 
 #include "utility/ROOTUtility.hxx"
 
 #include "exception/exception.hxx"
 
 #include "string_parsers/from_string.hxx"
 
 #include "TAxis.h"
 #include "TH1D.h"
 #include "TH1F.h"
 #include "TH2D.h"
 #include "TH2F.h"
 #include "TH2Poly.h"
 
 #include <iomanip>
 #include <memory>
 #include <string>
 #include <vector>
 
 namespace nuis {
 namespace utility {
 
 NEW_NUIS_EXCEPT(unimplemented_GetHistogram_method);
 NEW_NUIS_EXCEPT(invalid_histogram_descriptor);
 NEW_NUIS_EXCEPT(invalid_histogram_name);
 NEW_NUIS_EXCEPT(failed_to_clone);
 
 inline bool IsFlowBin(TAxis const &ax, Int_t bin_it) {
   return ((bin_it <= 0) || (bin_it >= (ax.GetNbins() + 1)));
 }
 
 inline bool IsInHistogramRange(TAxis const &ax, double v) {
   Int_t bin_it = ax.FindFixBin(v);
   return !IsFlowBin(ax, bin_it);
 }
 
 template <typename HT> struct HType_traits {};
 template <> struct HType_traits<TH1> {
   using type = TH1;
   static size_t const NDim = 1;
   using NumericT = double;
   static std::string name() { return "TH1"; }
 };
 template <> struct HType_traits<TH1D> {
   using type = TH1D;
   static size_t const NDim = 1;
   using NumericT = double;
   static std::string name() { return "TH1D"; }
 };
 template <> struct HType_traits<TH1F> {
   using type = TH1F;
   static size_t const NDim = 1;
   using NumericT = float;
   static std::string name() { return "TH1F"; }
 };
 template <> struct HType_traits<TH2> {
   using type = TH2;
   static size_t const NDim = 2;
   using NumericT = double;
   static std::string name() { return "TH2"; }
 };
 template <> struct HType_traits<TH2D> {
   using type = TH2D;
   static size_t const NDim = 2;
   using NumericT = double;
   static std::string name() { return "TH2D"; }
 };
 template <> struct HType_traits<TH2F> {
   using type = TH2F;
   static size_t const NDim = 2;
   using NumericT = float;
   static std::string name() { return "TH2F"; }
 };
 template <> struct HType_traits<TH2Poly> {
   using type = TH2Poly;
   static size_t const NDim = 0;
   using NumericT = double;
   static std::string name() { return "TH2Poly"; }
 };
 
 template <size_t, typename T = double> struct HType_Helper {};
 template <> struct HType_Helper<1, void> {
   using type = TH1;
   static size_t const NDim = HType_traits<type>::NDim;
   using NumericT = HType_traits<type>::NumericT;
   static std::string name() { return HType_traits<type>::name(); }
 };
 template <> struct HType_Helper<1, double> {
   using type = TH1D;
   static size_t const NDim = HType_traits<type>::NDim;
   using NumericT = HType_traits<type>::NumericT;
   static std::string name() { return HType_traits<type>::name(); }
 };
 template <> struct HType_Helper<1, float> {
   using type = TH1F;
   static size_t const NDim = HType_traits<type>::NDim;
   using NumericT = HType_traits<type>::NumericT;
   static std::string name() { return HType_traits<type>::name(); }
 };
 template <> struct HType_Helper<2, void> {
   using type = TH2;
   static size_t const NDim = HType_traits<type>::NDim;
   using NumericT = HType_traits<type>::NumericT;
   static std::string name() { return HType_traits<type>::name(); }
 };
 template <> struct HType_Helper<2, double> {
   using type = TH2D;
   static size_t const NDim = HType_traits<type>::NDim;
   using NumericT = HType_traits<type>::NumericT;
   static std::string name() { return HType_traits<type>::name(); }
 };
 template <> struct HType_Helper<2, float> {
   using type = TH2F;
   static size_t const NDim = HType_traits<type>::NDim;
   using NumericT = HType_traits<type>::NumericT;
   static std::string name() { return HType_traits<type>::name(); }
 };
 
 template <typename HT, typename Enable = void> struct TH_Helper {};
 
 template <typename HT>
 struct TH_Helper<HT,
                  typename std::enable_if<HType_traits<HT>::NDim == 1>::type> {
   static size_t const NDim = HType_traits<HT>::NDim;
   using NumericT = typename HType_traits<HT>::NumericT;
   static std::string name() { return HType_traits<HT>::name(); }
 
   static Int_t NbinsIncludeFlow(HT const &h) {
     return h.GetXaxis()->GetNbins() + 2;
   }
   static Int_t Nbins(HT const &h) { return h.GetXaxis()->GetNbins(); }
 
   static bool IsFlowBin(HT const &h, Int_t bin_it) {
     return nuis::utility::IsFlowBin(*h.GetXaxis(), bin_it);
   }
 
-  static void Fill(HT &h, std::array<NumericT, NDim> const &v,
-                   double w = 1) {
+  static void Fill(HT &h, std::array<NumericT, NDim> const &v, double w = 1) {
     h.Fill(v[0], w);
   }
 
   static void Scale(HT &h, NumericT SF, char const *opt = "") {
     h.Scale(SF, opt);
   }
 
   static double Integral(HT const &h, char const *opt = "") {
     return h.Integral(opt);
   }
 
   static Int_t NbinsIncludeFlow(std::unique_ptr<HT> const &h) {
     return NbinsIncludeFlow(*h);
   }
   static Int_t Nbins(std::unique_ptr<HT> const &h) { return Nbins(*h); }
 
   static bool IsFlowBin(std::unique_ptr<HT> const &h, Int_t bin_it) {
     return IsFlowBin(*h, bin_it);
   }
 
-  static void Fill(std::unique_ptr<HT> &h,
-                   std::array<NumericT, NDim> const &v, double w = 1) {
+  static void Fill(std::unique_ptr<HT> &h, std::array<NumericT, NDim> const &v,
+                   double w = 1) {
     Fill(*h, v, w);
   }
 
   static void Scale(std::unique_ptr<HT> &h, NumericT SF, char const *opt = "") {
     Scale(*h, SF, opt);
   }
 
   static double Integral(std::unique_ptr<HT> const &h, char const *opt = "") {
     return Integral(*h, opt);
   }
 };
 
 template <typename HT>
 struct TH_Helper<HT,
                  typename std::enable_if<HType_traits<HT>::NDim == 2>::type> {
   static size_t const NDim = HType_traits<HT>::NDim;
   using NumericT = typename HType_traits<HT>::NumericT;
   static std::string name() { return HType_traits<HT>::name(); }
 
   // TH2 ***************************************************************
   static Int_t NbinsIncludeFlow(HT const &h) {
     return (h.GetXaxis()->GetNbins() + 2) * (h.GetYaxis()->GetNbins() + 2);
   }
   static Int_t Nbins(HT const &h) {
     return (h.GetXaxis()->GetNbins()) * (h.GetYaxis()->GetNbins());
   }
 
   static bool IsFlowBin(HT const &h, Int_t xbin_it, Int_t ybin_it) {
     return nuis::utility::IsFlowBin(*h.GetXaxis(), xbin_it) ||
            nuis::utility::IsFlowBin(*h.GetYaxis(), ybin_it);
   }
 
-  static void Fill(HT &h, std::array<NumericT, NDim> const &v,
-                   double w = 1) {
+  static void Fill(HT &h, std::array<NumericT, NDim> const &v, double w = 1) {
     h.Fill(v[0], v[1], w);
   }
 
   static void Scale(HT &h, NumericT SF, char const *opt = "") {
     h.Scale(SF, opt);
   }
 
   static double Integral(HT const &h, char const *opt = "") {
     return h.Integral(opt);
   }
 
   static Int_t NbinsIncludeFlow(std::unique_ptr<HT> const &h) {
     return NbinsIncludeFlow(*h);
   }
   static Int_t Nbins(std::unique_ptr<HT> const &h) { return Nbins(*h); }
 
   static bool IsFlowBin(std::unique_ptr<HT> const &h, Int_t xbin_it,
                         Int_t ybin_it) {
     return IsFlowBin(*h, xbin_it, ybin_it);
   }
 
-  static void Fill(std::unique_ptr<HT> &h,
-                   std::array<NumericT, NDim> const &v, double w = 1) {
+  static void Fill(std::unique_ptr<HT> &h, std::array<NumericT, NDim> const &v,
+                   double w = 1) {
     Fill(*h, v, w);
   }
 
   static void Scale(std::unique_ptr<HT> &h, NumericT SF, char const *opt = "") {
     Scale(*h, SF, opt);
   }
 
   static double Integral(std::unique_ptr<HT> const &h, char const *opt = "") {
     return Integral(*h, opt);
   }
 };
 
 template <typename HT>
 struct TH_Helper<
     HT, typename std::enable_if<std::is_same<HT, TH2Poly>::value>::type> {
   static size_t const NDim = 2;
   using NumericT = typename HType_traits<HT>::NumericT;
   static std::string name() { return HType_traits<HT>::name(); }
 
   // TH2Poly ***************************************************************
   static Int_t NbinsIncludeFlow(HT const &h) { return h.GetNumberOfBins() + 9; }
   static Int_t Nbins(HT const &h) { return h.GetNumberOfBins(); }
 
   static bool IsFlowBin(HT const &h, Int_t bin_it) { return (bin_it < 0); }
 
-  static void Fill(HT &h, std::array<NumericT, NDim> const &v,
-                   double w = 1) {
+  static void Fill(HT &h, std::array<NumericT, NDim> const &v, double w = 1) {
     h.Fill(v[0], v[1], w);
   }
 
   static void Scale(HT &h, NumericT SF, char const *opt = "") {
 
     bool width = (std::string(opt).find("width") != std::string::npos);
     size_t nbins = Nbins(h);
     for (size_t bin_it = 0; bin_it < nbins; ++bin_it) {
       double bin_area = 1;
 
       if (width) {
         TH2PolyBin *poly_bin =
             dynamic_cast<TH2PolyBin *>(h.GetBins()->At(bin_it));
 
         bin_area = poly_bin->GetArea();
       }
 
       h.SetBinContent(bin_it + 1,
                       h.GetBinContent(bin_it + 1) * (SF / bin_area));
       h.SetBinError(bin_it + 1, h.GetBinError(bin_it + 1) * (SF / bin_area));
     }
   }
 
   static double Integral(HT &h, char const *opt = "") {
     bool width = (std::string(opt).find("width") != std::string::npos);
     size_t nbins = Nbins(h);
     double integral = 0;
     for (size_t bin_it = 0; bin_it < nbins; ++bin_it) {
       double bin_area = 1;
 
       if (width) {
         TH2PolyBin *poly_bin =
             dynamic_cast<TH2PolyBin *>(h.GetBins()->At(bin_it));
 
         bin_area = poly_bin->GetArea();
       }
 
       integral += h.GetBinContent(bin_it + 1) * bin_area;
     }
     return integral;
   }
 
   static Int_t NbinsIncludeFlow(std::unique_ptr<HT> const &h) {
     return NbinsIncludeFlow(*h);
   }
   static Int_t Nbins(std::unique_ptr<HT> const &h) { return Nbins(*h); }
 
   static bool IsFlowBin(std::unique_ptr<HT> const &h, Int_t bin_it) {
     return IsFlowBin(*h, bin_it);
   }
 
-  static void Fill(std::unique_ptr<HT> &h,
-                   std::array<NumericT, NDim> const &v, double w = 1) {
+  static void Fill(std::unique_ptr<HT> &h, std::array<NumericT, NDim> const &v,
+                   double w = 1) {
     Fill(*h, v, w);
   }
 
   static void Scale(std::unique_ptr<HT> &h, NumericT SF, char const *opt = "") {
     Scale(*h, SF, opt);
   }
 
   static double Integral(std::unique_ptr<HT> &h, char const *opt = "") {
     return Integral(*h, opt);
   }
 };
 
 template <typename HT>
 void Clear(typename std::enable_if<HType_traits<HT>::NDim != 0, HT>::type &h) {
   for (Int_t bin_it = 0; bin_it < TH_Helper<HT>::NbinsIncludeFlow(h);
        ++bin_it) {
     h.SetBinContent(bin_it, 0);
     h.SetBinError(bin_it, 0);
   }
 }
 
 template <typename HT>
 void Clear(
     typename std::enable_if<std::is_same<HT, TH2Poly>::value, HT>::type &h) {
   h.ClearBinContents();
 }
 
 template <typename HT>
 inline std::unique_ptr<HT> Clone(HT const &source, bool clear = false,
                                  std::string const &clone_name = "") {
   std::unique_ptr<HT> target(dynamic_cast<HT *>(
       source.Clone(clone_name.size() ? clone_name.c_str() : "")));
   if (!target) {
     throw failed_to_clone()
         << "[ERROR]: Failed to clone a " << TH_Helper<HT>::name() << ".";
   }
   target->SetDirectory(nullptr);
 
   if (clear) {
     Clear<HT>(*target);
   }
 
   return target;
 }
 template <typename HT>
 inline std::unique_ptr<HT> Clone(std::unique_ptr<HT> const &source,
                                  bool clear = false,
                                  std::string const &clone_name = "") {
   return Clone(*source, clear, clone_name);
 }
 
 template <typename HT>
-inline std::unique_ptr<HT> GetHistogramFromROOTFile(TFile_ptr &f,
-                                                    std::string const &hname) {
+inline std::unique_ptr<HT>
+GetHistogramFromROOTFile(TFile_ptr &f, std::string const &hname,
+                         bool ThrowIfMissing = true) {
 
   f->Get(hname.c_str());
   HT *h = dynamic_cast<HT *>(f->Get(hname.c_str()));
   if (!h) {
-    throw invalid_histogram_name()
-        << "[ERROR]: Failed to get " << TH_Helper<HT>::name() << " named "
-        << std::quoted(hname) << " from input file "
-        << std::quoted(f->GetName());
+    if (ThrowIfMissing) {
+      throw invalid_histogram_name()
+          << "[ERROR]: Failed to get " << TH_Helper<HT>::name() << " named "
+          << std::quoted(hname) << " from input file "
+          << std::quoted(f->GetName());
+    } else {
+      return nullptr;
+    }
   }
   std::unique_ptr<HT> clone = Clone<HT>(*h);
   return clone;
 }
 
 template <typename HT>
-inline std::unique_ptr<HT> GetHistogramFromROOTFile(std::string const &fname,
-                                                    std::string const &hname) {
+inline std::unique_ptr<HT>
+GetHistogramFromROOTFile(std::string const &fname, std::string const &hname,
+                         bool ThrowIfMissing = true) {
   TFile_ptr temp = CheckOpenTFile(fname, "READ");
-  return GetHistogramFromROOTFile<HT>(temp, hname);
+  return GetHistogramFromROOTFile<HT>(temp, hname, ThrowIfMissing);
 }
 
 template <typename HT>
 std::unique_ptr<HT> GetHistogram(std::string const &input_descriptor) {
   std::vector<std::string> split_descriptor =
       fhicl::string_parsers::ParseToVect<std::string>(input_descriptor, ";",
                                                       true, true);
 
   if (split_descriptor.size() == 1) { // assume text
     throw unimplemented_GetHistogram_method();
   } else if (split_descriptor.size() == 2) {
     return GetHistogramFromROOTFile<HT>(split_descriptor[0],
                                         split_descriptor[1]);
   } else {
     throw invalid_histogram_descriptor()
         << "[ERROR]: Failed to parse histogram descriptor: "
         << std::quoted(input_descriptor)
         << " as an input histogram (Text/ROOT).";
   }
 }
 
 NEW_NUIS_EXCEPT(invalid_TH2Poly);
 NEW_NUIS_EXCEPT(invalid_PolyBinSpecifierList);
 
 // List of bins to put in a row, X/Y
 struct PolyBinSpecifier {
   double X, Y;
   bool UseXAxis;
 };
 
 constexpr PolyBinSpecifier XPolyBinSpec(double X, double Y) {
   return PolyBinSpecifier{X + std::numeric_limits<double>::epsilon() * 1E2,
                           Y + std::numeric_limits<double>::epsilon() * 1E2,
                           true};
 }
 
 constexpr PolyBinSpecifier YPolyBinSpec(double X, double Y) {
   return PolyBinSpecifier{X + std::numeric_limits<double>::epsilon() * 1E2,
                           Y + std::numeric_limits<double>::epsilon() * 1E2,
                           false};
 }
 
 std::vector<std::unique_ptr<TH1>> GetTH2PolySlices(
     std::unique_ptr<TH2Poly> &hinp,
     std::vector<std::vector<PolyBinSpecifier>> const &BinsSpecifiers) {
 
   std::vector<std::unique_ptr<TH1>> slices;
 
   size_t sl_it = 0;
   for (auto &slice_spec : BinsSpecifiers) {
     std::vector<double> Binning;
     std::vector<double> BinContent;
     std::vector<double> BinError;
     bool UseXAxis = false;
     size_t bin_ctr = 0;
     for (auto poly_bin_spec : slice_spec) {
       Int_t bin_it = hinp->FindBin(poly_bin_spec.X, poly_bin_spec.Y);
       if (bin_it < 1) {
         std::cout << "[WARN]: When searching for matching bin: { X: "
                   << poly_bin_spec.X << ", Y: " << poly_bin_spec.Y
                   << "} got flow bin: " << bin_it << std::endl;
         continue;
       }
       TH2PolyBin *poly_bin =
           dynamic_cast<TH2PolyBin *>(hinp->GetBins()->At(bin_it - 1));
 
       if (!bin_ctr) {
         UseXAxis = poly_bin_spec.UseXAxis;
       } else if (UseXAxis != poly_bin_spec.UseXAxis) {
         throw invalid_PolyBinSpecifierList()
             << "[ERROR]: For slice: " << sl_it
             << " of TH2Poly: " << std::quoted(hinp->GetName())
             << " bin specifier: " << bin_ctr << " was set to use the "
             << (poly_bin_spec.UseXAxis ? "X" : "Y")
             << " axis as the dependent axis of the slice, but previous bins "
                "were set to use the "
             << (poly_bin_spec.UseXAxis ? "X" : "Y") << " axis.";
       }
 
       if ((poly_bin_spec.X < poly_bin->GetXMin()) ||
           (poly_bin_spec.X >= poly_bin->GetXMax()) ||
           (poly_bin_spec.Y < poly_bin->GetYMin()) ||
           (poly_bin_spec.Y >= poly_bin->GetYMax())) {
         std::cout
             << "[WARN]: Found bin doesn't seem to contain expected point: { X: "
             << poly_bin_spec.X << ", Y: " << poly_bin_spec.Y
             << "}, got bin_it = " << bin_it
             << " which had x_low: " << poly_bin->GetXMin()
             << ", and x_up: " << poly_bin->GetXMax()
             << ", y_low: " << poly_bin->GetYMin()
             << ", y_up: " << poly_bin->GetYMax() << std::endl;
       }
 
       double low = UseXAxis ? poly_bin->GetXMin() : poly_bin->GetYMin();
       double up = UseXAxis ? poly_bin->GetXMax() : poly_bin->GetYMax();
 
       if (!Binning.size()) { // Add low edge
         Binning.push_back(low);
       } else if (std::abs(Binning.back() - low) >
                  (std::numeric_limits<double>::epsilon() * 1E2)) {
         BinContent.push_back(0);
         BinError.push_back(0);
         Binning.push_back(low);
       }
 
       BinContent.push_back(hinp->GetBinContent(bin_it));
       BinError.push_back(hinp->GetBinError(bin_it));
       Binning.push_back(up);
     }
     slices.emplace_back(new TH1D(
         (std::string(hinp->GetName()) + "_slice" + std::to_string(sl_it++))
             .c_str(),
         (std::string(";") +
          (UseXAxis ? hinp->GetXaxis() : hinp->GetYaxis())->GetTitle() + ";" +
          hinp->GetZaxis()->GetTitle())
             .c_str(),
         Binning.size() - 1, Binning.data()));
 
     for (size_t bin_it = 0; bin_it < BinContent.size(); ++bin_it) {
       slices.back()->SetBinContent(bin_it + 1, BinContent[bin_it]);
       slices.back()->SetBinError(bin_it + 1, BinError[bin_it]);
     }
   }
   return slices;
 }
 
 } // namespace utility
 } // namespace nuis
 #endif