diff --git a/src/app/nuissamples.cxx b/src/app/nuissamples.cxx
index 809f275..4be6a8f 100644
--- a/src/app/nuissamples.cxx
+++ b/src/app/nuissamples.cxx
@@ -1,155 +1,156 @@
 #include "config/GlobalConfiguration.hxx"
 
 #include "input/IInputHandler.hxx"
 
 #include "plugins/Instantiate.hxx"
 
 #include "samples/IDataComparison.hxx"
 
 #include "utility/StringUtility.hxx"
 
 #include <fstream>
 #include <iostream>
 #include <regex>
 #include <string>
 
 std::string name_search_term = ".*";
 bool strict_name_regex = false;
 
 std::string target_search_term = ".*";
 bool strict_target_regex = false;
 
 std::string year_search_term = ".*";
 
 std::string config_out_filename = "";
 
 bool NameOnly = false;
 
 void SayUsage(char const *argv[]) {
-  std::stringstream ss("");
-  ss << "[USAGE]: " << argv[0]
+  std::cout << "[USAGE]: " << argv[0]
      << "\n"
         "\t-n,-N,-t,-T,-y <search regex> : Filters known IDataComparisons by "
-        "the search term (-n: Name, -t: Target, -y: Year). Capitalized "
-        "versions filter on exact match."
+        "the search term \n"
+        "\t                                (-n: Name, -t: Target, -y: Year). Capitalized versions \n"
+        "\t                                filter on exact match.\n"
         "\t-o <output_file>              : Dump example sample configuration "
-        "file for matching samples.\n"
+        "file for matching \n"
+        "\t                                samples.\n"
         "\t--name-only                   : Only write out matching sample "
-        "names. (Still applies all search terms)\n";
+        "names. \n"
+        "\t                                (Still applies all search terms)\n" << std::endl;
 
-  std::cout << nuis::utility::indent_apply_width(ss.str()) << std::endl;
 }
 
 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]) == "-s") {
+    } else if (std::string(argv[opt]) == "-n") {
       name_search_term = argv[++opt];
-    } else if (std::string(argv[opt]) == "-S") {
+    } else if (std::string(argv[opt]) == "-N") {
       strict_name_regex = true;
       name_search_term = argv[++opt];
     } else if (std::string(argv[opt]) == "-t") {
       target_search_term = argv[++opt];
     } else if (std::string(argv[opt]) == "-T") {
       strict_target_regex = true;
       target_search_term = argv[++opt];
     } else if (std::string(argv[opt]) == "-y") {
       year_search_term = argv[++opt];
     } else if (std::string(argv[opt]) == "-o") {
       config_out_filename = argv[++opt];
     } else if (std::string(argv[opt]) == "--name-only") {
       NameOnly = true;
     } else {
       std::cout << "[ERROR]: Unknown option: " << argv[opt] << std::endl;
       SayUsage(argv);
       exit(1);
     }
     opt++;
   }
 }
 
 NEW_NUIS_EXCEPT(invalid_output_file);
 
 int main(int argc, char const *argv[]) {
   nuis::config::EnsureConfigurationRead("nuis.global.config.fcl");
   nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl");
 
   handleOpts(argc, argv);
 
   std::regex rpattern_name(strict_name_regex
                                ? name_search_term
                                : std::string(".*") + name_search_term + ".*");
   std::regex rpattern_target(
       strict_target_regex ? target_search_term
                           : std::string(".*") + target_search_term + ".*");
 
   std::regex rpattern_year(std::string(".*") + year_search_term + ".*");
 
   std::vector<fhicl::ParameterSet> example_sample_configs;
 
   for (std::string const &comparison_set_key :
        nuis::config::GetDocument()
            .get<fhicl::ParameterSet>("data_comparisons")
            .get_names()) {
     for (std::string const &sample_name :
          nuis::config::GetDocument().get<std::vector<std::string>>(
              std::string("data_comparisons.") + comparison_set_key)) {
 
       if (!std::regex_match(sample_name, rpattern_name)) {
         continue;
       }
 
       nuis::plugins::plugin_traits<IDataComparison>::unique_ptr_t sample =
           nuis::plugins::Instantiate<IDataComparison>(sample_name);
 
       if (!std::regex_match(sample->GetTargetMaterial(), rpattern_target)) {
         continue;
       }
 
       if (!std::regex_match(sample->GetYear(), rpattern_year)) {
         continue;
       }
 
       std::cout << sample->Name() << std::endl;
       if (!NameOnly) {
         std::cout << "\tJournal: " << sample->GetJournalReference()
                   << std::endl;
         std::cout << "\tDOI: " << sample->GetDOI() << std::endl;
         std::cout << "\tYear: " << sample->GetYear() << std::endl;
         std::cout << "\tTarget: " << sample->GetTargetMaterial() << std::endl;
         std::cout << "\tFlux: " << sample->GetFluxDescription() << std::endl;
         std::cout << "\tSignal: " << sample->GetSignalDescription()
                   << std::endl;
         std::cout << "\tDocs: \n"
                   << nuis::utility::indent_apply_width(
                          sample->GetDocumentation(), 10)
                   << std::endl;
         std::cout << "\tExample_Config: {\n"
                   << nuis::utility::indent_apply_width(
                          sample->GetExampleConfiguration().to_indented_string(),
                          12)
                   << "\n\t}\n"
                   << std::endl;
       }
 
       example_sample_configs.push_back(sample->GetExampleConfiguration());
     }
   }
 
   if (config_out_filename.length()) {
     std::ofstream out_file(config_out_filename);
 
     if (!out_file) {
       throw invalid_output_file() << "[ERROR]: Failed to open output file: "
                                   << std::quoted(config_out_filename);
     }
 
     fhicl::ParameterSet example_config;
     example_config.put("samples", example_sample_configs);
     out_file << example_config.to_indented_string();
   }
 }
diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx
index 52cadce..8968f3a 100644
--- a/src/generator/input/NEUTInputHandler.cxx
+++ b/src/generator/input/NEUTInputHandler.cxx
@@ -1,250 +1,252 @@
 #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");
 
   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();
+    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;
   }
 }
 
 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;
 }
 
-void NEUTInputHandler::RebuildEventRate() {
+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();
+    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 {
   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() {
   return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT");
 }
 
 DECLARE_PLUGIN(IInputHandler, NEUTInputHandler);
diff --git a/src/generator/input/NEUTInputHandler.hxx b/src/generator/input/NEUTInputHandler.hxx
index 36050cb..456b8cc 100644
--- a/src/generator/input/NEUTInputHandler.hxx
+++ b/src/generator/input/NEUTInputHandler.hxx
@@ -1,70 +1,70 @@
 // 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();
-  void RebuildEventRate();
+  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;
   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/persistency/ROOTOutput.cxx b/src/persistency/ROOTOutput.cxx
index 619a820..2de36b1 100644
--- a/src/persistency/ROOTOutput.cxx
+++ b/src/persistency/ROOTOutput.cxx
@@ -1,81 +1,81 @@
 #include "persistency/ROOTOutput.hxx"
 
 #include "utility/ROOTUtility.hxx"
 
 #include "config/GlobalConfiguration.hxx"
 
 #include "fhiclcpp/ParameterSet.h"
 
 #include "TFile.h"
 
 namespace nuis {
 namespace persistency {
 
 struct NamedTFile {
   std::string name;
   std::unique_ptr<TFile> file;
 
   NamedTFile() : name(""), file(nullptr) {}
   NamedTFile(NamedTFile &&other)
       : name(std::move(other.name)), file(std::move(other.file)) {}
   ~NamedTFile() {
     if (file) {
       file->Write();
       file->Close();
     }
   }
 };
 
 static std::vector<NamedTFile> Files;
 
 std::unique_ptr<TFile> &GetOutputFile(std::string const &name) {
   for (NamedTFile &file : Files) {
     if (file.name == name) {
       return file.file;
     }
   }
 
   fhicl::ParameterSet const &persistency =
       config::GetDocument().get<fhicl::ParameterSet>("persistency");
   std::string file_name = persistency.get<std::string>(name + ".output_file");
   std::string open_opts =
       persistency.get<std::string>(name + ".open_mode", "CREATE");
 
   NamedTFile ntf;
   ntf.name = name;
   ntf.file = std::make_unique<TFile>(file_name.c_str(), open_opts.c_str());
   if (!ntf.file || !ntf.file->IsOpen()) {
     if ((name == "default") && (open_opts == "CREATE")) {
       std::cout
           << "[WARN]: It appears the default file cannot be opened because it "
              "exists already, to stop you wasting processing time, the default "
              "output stream will be written to nuis.default.tmp.root in the "
-             "current director in RECREATE mode."
+             "current directory in RECREATE mode."
           << std::endl;
       file_name = "nuis.default.tmp.root";
       open_opts = "RECREATE";
       ntf.file = std::make_unique<TFile>(file_name.c_str(), open_opts.c_str());
       if (!ntf.file || !ntf.file->IsOpen()) {
         throw utility::failed_to_open_TFile()
             << "[ERROR]: Failed to open output file: " << std::quoted(file_name)
             << " in write mode with opts = " << std::quoted(open_opts);
       }
     } else {
       throw utility::failed_to_open_TFile()
           << "[ERROR]: Failed to open output file: " << std::quoted(file_name)
           << " in write mode with opts = " << std::quoted(open_opts);
     }
   }
   Files.push_back(std::move(ntf));
   return Files.back().file;
 }
 
 void CloseOpenTFiles() {
   for (NamedTFile &f : Files) {
     std::cout << "[INFO]: Closing open TFile: " << f.name << " "
               << f.file->GetName() << std::endl;
   }
   Files.clear();
 }
 } // namespace persistency
 } // namespace nuis
diff --git a/src/samples/SimpleDataComparison.hxx b/src/samples/SimpleDataComparison.hxx
index b3a5e17..b8f91e7 100644
--- a/src/samples/SimpleDataComparison.hxx
+++ b/src/samples/SimpleDataComparison.hxx
@@ -1,435 +1,451 @@
 // 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/IDataComparison.hxx"
 
 #include "event/FullEvent.hxx"
 
 #include "input/InputManager.hxx"
 
 #include "persistency/ROOTOutput.hxx"
 
 #include "utility/FileSystemUtility.hxx"
 #include "utility/HistogramUtility.hxx"
 #include "utility/KinematicUtility.hxx"
 #include "utility/ROOTUtility.hxx"
 #include "utility/StatsUtility.hxx"
 
 #include <array>
 #include <functional>
 #include <limits>
 #include <memory>
 #include <vector>
 
 template <size_t nd, typename NumericT = double,
           typename HT = typename nuis::utility::HType_Helper<nd, void>::type>
 class SimpleDataComparison : public IDataComparison {
 
   NEW_NUIS_EXCEPT(invalid_SimpleDataComparison_initialization);
   NEW_NUIS_EXCEPT(SimpleDataComparison_already_finalized);
 
 protected:
   using HistType = HT;
   using TH_Help = typename nuis::utility::TH_Helper<HistType>;
 
   static size_t const NDim = nd;
 
   nuis::input::InputManager::Input_id_t fIH_id;
   std::string write_directory;
   size_t NMaxSample_override;
   int fIsShapeOnly;
   int fIsFluxUnfolded;
 
   std::vector<bool> fSignalCache;
   std::vector<std::array<NumericT, NDim>> fProjectionCache;
+  bool fUseCache;
 
   std::string fDataInputDescriptor;
   std::unique_ptr<HistType> fData;
   std::string fMaskInputDescriptor;
   std::unique_ptr<HistType> fMask;
   std::string fCovarianceInputDescriptor;
   std::unique_ptr<TH2> fCovariance;
   std::unique_ptr<HistType> fPrediction;
   std::unique_ptr<HistType> fPrediction_xsec;
   std::unique_ptr<HistType> fPrediction_shape;
   std::unique_ptr<HistType> fPrediction_comparison;
   bool fComparisonFinalized;
 
   std::string fJournalReference;
   std::string fDOI;
   std::string fYear;
   std::string fTargetMaterial;
   std::string fFluxDescription;
   std::string fSignalDescription;
 
   nuis::utility::ENuRange energy_cut;
 
   std::function<bool(nuis::event::FullEvent const &)> IsSigFunc;
   std::function<std::array<NumericT, NDim>(nuis::event::FullEvent const &)>
       CompProjFunc;
 
   // If assigned by subclass will be called on for all events, bool signifies
   // whether the event was selected.
   std::function<void(nuis::event::FullEvent const &, bool, double)>
       ProcessExtraFunc;
 
 public:
   SimpleDataComparison() {
     fIH_id = std::numeric_limits<nuis::input::InputManager::Input_id_t>::max();
+    fUseCache = false;
     write_directory = "";
     NMaxSample_override = std::numeric_limits<size_t>::max();
     fDataInputDescriptor = "";
     fData = nullptr;
     fMaskInputDescriptor = "";
     fMask = nullptr;
     fCovarianceInputDescriptor = "";
     fCovariance = nullptr;
     fPrediction = nullptr;
     fPrediction_xsec = nullptr;
     fPrediction_shape = nullptr;
     fPrediction_comparison = nullptr;
     fComparisonFinalized = false;
     IsSigFunc = [](nuis::event::FullEvent const &) -> bool { return true; };
     CompProjFunc =
         [](nuis::event::FullEvent const &) -> std::array<NumericT, NDim> {
       std::array<NumericT, NDim> arr;
       for (NumericT &el : arr) {
         el = 0;
       }
       return arr;
     };
 
     ProcessExtraFunc =
         std::function<void(nuis::event::FullEvent const &, bool, double)>();
 
     fJournalReference = "";
     fDOI = "";
     fYear = "";
     fTargetMaterial = "";
     fFluxDescription = "";
     fSignalDescription = "";
     fIsShapeOnly = -1;
     fIsFluxUnfolded = -1;
 
     energy_cut = nuis::utility::ENuRange{std::numeric_limits<double>::max(),
                                          std::numeric_limits<double>::max()};
   }
 
+  void SetUseCache(bool uc = true) {
+    fUseCache = uc;
+    if (!uc) {
+      fSignalCache.resize(0);
+      fProjectionCache.resize(0);
+    }
+  }
+
   fhicl::ParameterSet fGlobalConfig;
   fhicl::ParameterSet fInstanceConfig;
 
   void ReadGlobalConfigDefaults() {
     fGlobalConfig = nuis::config::GetDocument().get<fhicl::ParameterSet>(
         std::string("global.sample_configuration.") + Name(),
         fhicl::ParameterSet());
 
     if (!fJournalReference.length()) {
       fJournalReference = fGlobalConfig.get<std::string>("journal_reference",
                                                          "Not yet published");
     }
     if (!fDOI.length()) {
       fDOI = fGlobalConfig.get<std::string>("DOI", "Unknown DOI");
     }
     if (!fYear.length()) {
       fYear = fGlobalConfig.get<std::string>("year", "Unknown year");
     }
 
     if (!fTargetMaterial.length()) {
       fTargetMaterial = fGlobalConfig.get<std::string>(
           "target_material", "Unknown target material");
     }
     if (!fFluxDescription.length()) {
       fFluxDescription =
           fGlobalConfig.get<std::string>("flux_description", "Unknown flux");
     }
     if (!fSignalDescription.length()) {
       fSignalDescription = fGlobalConfig.get<std::string>("signal_description",
                                                           "Unknown Signal");
     }
 
     if (fIsShapeOnly == -1) {
       fIsShapeOnly = fGlobalConfig.get<bool>("shape_only", false);
     }
     if (fIsFluxUnfolded == -1) {
       fIsFluxUnfolded = fGlobalConfig.get<bool>("flux_unfolded", false);
     }
 
     if ((energy_cut.first == std::numeric_limits<double>::max()) &&
         (fGlobalConfig.has_key("enu_range"))) {
       energy_cut = fGlobalConfig.get<std::pair<double, double>>("enu_range");
     }
   }
 
   virtual std::string GetJournalReference() { return fJournalReference; }
   virtual std::string GetDOI() { return fDOI; }
   virtual std::string GetYear() { return fYear; }
   virtual std::string GetTargetMaterial() { return fTargetMaterial; }
   virtual std::string GetFluxDescription() { return fFluxDescription; }
   virtual std::string GetSignalDescription() { return fSignalDescription; }
 
   void SetShapeOnly(bool iso) { fIsShapeOnly = iso; }
   void SetFluxUnfolded(bool ifo) { fIsFluxUnfolded = ifo; }
 
   void SetData(std::string const &data_descriptor) {
     fDataInputDescriptor = data_descriptor;
   }
 
   void SetMask(std::string const &mask_descriptor) {
     fMaskInputDescriptor = mask_descriptor;
   }
 
   void SetCovariance(std::string const &cov_descriptor) {
     fCovarianceInputDescriptor = cov_descriptor;
   }
 
   virtual void FillProjection(std::array<NumericT, NDim> const &proj,
                               NumericT event_weight) {
     TH_Help::Fill(fPrediction, proj, event_weight);
   }
 
   virtual void FinalizeComparison() {
     if (fComparisonFinalized) {
       throw SimpleDataComparison_already_finalized()
           << "[ERROR]: Attempted to re-finalize a comparison for "
              "SimpleDataComparison: "
           << std::quoted(Name());
     }
     fPrediction_xsec =
         nuis::utility::Clone(fPrediction, false, "Prediction_xsec");
 
     IInputHandler const &IH =
         nuis::input::InputManager::Get().GetInputHandler(fIH_id);
 
     TH_Help::Scale(fPrediction_xsec, 1.0, "width");
 
     // If we have a flux cut
     if (energy_cut.first != std::numeric_limits<double>::max()) {
       TH_Help::Scale(fPrediction_xsec, IH.GetXSecScaleFactor(energy_cut));
     }
 
     fPrediction_shape =
         nuis::utility::Clone(fPrediction_xsec, false, "Prediction_shape");
     if (fData) {
       TH_Help::Scale(fPrediction_shape,
                      TH_Help::Integral(fData, "width") /
                          TH_Help::Integral(fPrediction_shape, "width"));
 
     } else {
       IEventProcessor_WARN(
           "When Finalizing comparison, no Data histogram available.");
     }
 
     if (fIsFluxUnfolded) {
       // fPrediction_comparison
     } else if (fIsShapeOnly) {
       fPrediction_comparison = nuis::utility::Clone(fPrediction_shape, false,
                                                     "Prediction_comparison");
     } else {
       fPrediction_comparison = nuis::utility::Clone(fPrediction_xsec, false,
                                                     "Prediction_comparison");
     }
     fComparisonFinalized = true;
   }
 
   void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
 
     fInstanceConfig = instance_sample_configuration;
 
     SetSampleVerbosity(fInstanceConfig.get<std::string>(
         "verbosity", nuis::config::GetDocument().get<std::string>(
                          "global.sample.verbosity_default", "Reticent")));
 
     ReadGlobalConfigDefaults();
 
     if (fInstanceConfig.has_key("fake_data")) {
       fData = nuis::utility::GetHistogram<HistType>(
           fInstanceConfig.get<std::string>("fake_data_histogram"));
     } else if (!fGlobalConfig.get<bool>("has_data", true) ||
                !fInstanceConfig.get<bool>("has_data", true)) {
       // Explicitly not expecting data
     } else {
       if (!fDataInputDescriptor.length()) {
         if (!fGlobalConfig.has_key("data_descriptor")) {
           throw invalid_SimpleDataComparison_initialization()
               << "[ERROR]: SimpleDataComparison::Initialize for "
                  "IDataComparison: "
               << std::quoted(Name())
               << " failed as no input data was set by a call to "
                  "SimpleDataComparison::SetData and no data_descriptor for "
                  "this SimpleDataComparison could be found in the global "
                  "configuration.";
         }
         fDataInputDescriptor =
             fGlobalConfig.get<std::string>("data_descriptor");
       }
       fData = nuis::utility::GetHistogram<HistType>(fDataInputDescriptor);
     }
 
     if (!fPrediction) {
       if (!fData) {
         throw invalid_SimpleDataComparison_initialization()
             << "[ERROR]: SimpleDataComparison::Initialize for "
                "IDataComparison: "
             << std::quoted(Name())
             << " failed. As `has_data: false` was set in the configuration  "
                "(global: "
             << (!fGlobalConfig.get<bool>("has_data", true))
             << ", instance: " << !fInstanceConfig.get<bool>("has_data", true)
             << "), the instance constructor must supply the fPrediction "
                "binning, and it wasn't.";
       }
       fPrediction = nuis::utility::Clone(fData, true, "Prediction");
     }
 
     if (fCovarianceInputDescriptor.length()) {
       fCovariance =
           nuis::utility::GetHistogram<TH2D>(fCovarianceInputDescriptor);
     } else if (fGlobalConfig.has_key("covariance_descriptor")) {
       fCovariance = nuis::utility::GetHistogram<TH2D>(
           fGlobalConfig.get<std::string>("covariance_descriptor"));
     }
 
     if (fMaskInputDescriptor.length()) {
       fMask = nuis::utility::GetHistogram<HistType>(fMaskInputDescriptor);
     } else if (fGlobalConfig.has_key("mask_descriptor")) {
       fMask = nuis::utility::GetHistogram<HistType>(
           fGlobalConfig.get<std::string>("mask_descriptor"));
     }
 
     if (fInstanceConfig.has_key("verbosity")) {
       SetSampleVerbosity(fInstanceConfig.get<std::string>("verbosity"));
     }
 
     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;
 
-    bool DetermineSignalEvents = !fSignalCache.size();
+    bool DetermineSignalEvents = !fSignalCache.size() || !fUseCache;
 
     nuis::utility::Clear<HistType>(*fPrediction);
     fComparisonFinalized = false;
 
     size_t cache_ctr = 0;
 
     while (ev_idx < NEvsToProcess) {
+      bool is_sig = false;
+      std::array<NumericT, NDim> proj;
+
       if (DetermineSignalEvents) {
         nuis::event::FullEvent const &fev = IH.GetFullEvent(ev_idx);
-        bool is_sig = IsSigFunc(fev);
+        is_sig = IsSigFunc(fev);
         fSignalCache.push_back(is_sig);
         if (is_sig) {
-          fProjectionCache.push_back(CompProjFunc(fev));
+          proj = CompProjFunc(fev);
+          fProjectionCache.push_back(proj);
         }
         if (ProcessExtraFunc) {
           ProcessExtraFunc(fev, is_sig,
                            IH.GetEventWeight(ev_idx) * nmax_scaling);
         }
+      } else {
+        is_sig = fSignalCache[ev_idx];
+        proj = fProjectionCache[cache_ctr++];
       }
 
       if (NToShout && !(ev_idx % NToShout)) {
         IEventProcessor_INFO("\tProcessed " << ev_idx << "/" << NEvsToProcess
                                             << " events.");
       }
 
-      if (fSignalCache[ev_idx]) {
-        FillProjection(fProjectionCache[cache_ctr++],
-                       IH.GetEventWeight(ev_idx) * nmax_scaling);
+      if (is_sig) {
+        FillProjection(proj, IH.GetEventWeight(ev_idx) * nmax_scaling);
       }
 
       ev_idx++;
     }
     IEventProcessor_INFO("\t" << fProjectionCache.size() << "/" << NEvsToProcess
                               << " events passed selection.");
   }
   void Write() {
     if (!fComparisonFinalized) {
       FinalizeComparison();
     }
 
     nuis::persistency::WriteToOutputFile<HistType>(
         fPrediction_comparison, "Prediction", write_directory);
     nuis::persistency::WriteToOutputFile<HistType>(
         fPrediction_xsec, "Prediction_xsec", write_directory);
 
     if (fData) {
       nuis::persistency::WriteToOutputFile<HistType>(fData, "Data",
                                                      write_directory);
       nuis::persistency::WriteToOutputFile<HistType>(
           fPrediction_shape, "Prediction_shape", write_directory);
     }
   }
 
   double GetGOF() {
     if (!fComparisonFinalized) {
       FinalizeComparison();
     }
     if (fData && fPrediction_comparison) {
       return nuis::utility::GetChi2(fData, fPrediction_comparison);
     } else
       return std::numeric_limits<double>::max();
   }
   double GetNDOGuess() {
     if (fData) {
       return TH_Help::Nbins(fData);
     }
     return 0;
   }
 
   fhicl::ParameterSet GetExampleConfiguration() {
     fhicl::ParameterSet exps;
     exps.put<std::string>("name", Name());
     exps.put<std::string>("input_type", "NuWro/NEUT/GENIE");
     exps.put<std::string>("file", "input_events.root");
     exps.put<std::string>("write_directory", Name());
     return exps;
   }
 };
 
 typedef SimpleDataComparison<1, double, TH1D> SimpleDataComparison_1D;
 typedef SimpleDataComparison<2, double, TH2D> SimpleDataComparison_2D;
 typedef SimpleDataComparison<1, float, TH1F> SimpleDataComparison_1F;
 typedef SimpleDataComparison<2, float, TH2F> SimpleDataComparison_2F;
 typedef SimpleDataComparison<2, double, TH2Poly> SimpleDataComparison_2DPoly;
diff --git a/src/utility/HistogramUtility.hxx b/src/utility/HistogramUtility.hxx
index 40a8273..6010309 100644
--- a/src/utility/HistogramUtility.hxx
+++ b/src/utility/HistogramUtility.hxx
@@ -1,430 +1,461 @@
 // 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 "utility/ROOTUtility.hxx"
 
 #include "exception/exception.hxx"
 
 #include "string_parsers/from_string.hxx"
 
+#include "fhiclcpp/ParameterSet.h"
+
 #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);
 
 bool IsFlowBin(TAxis const &ax, Int_t bin_it);
 
 bool IsInHistogramRange(TAxis const &ax, double v);
 
 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) {
     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) {
     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) {
     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) {
     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) {
     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) {
     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,
                          bool ThrowIfMissing = true) {
 
   f->Get(hname.c_str());
   HT *h = dynamic_cast<HT *>(f->Get(hname.c_str()));
   if (!h) {
     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,
                          bool ThrowIfMissing = true) {
   TFile_ptr temp = CheckOpenTFile(fname, "READ");
   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).";
   }
 }
 
 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);
+
+template <typename HT>
+inline typename std::enable_if<HType_traits<HT>::NDim == 1,
+                               std::unique_ptr<TH1>>::type
+BuildHistFromFHiCL(fhicl::ParameterSet const &ps) {
+  std::array<double, 3> xaxis =
+      ps.get<std::array<double, 3>>("xaxis_descriptor");
+  std::unique_ptr<TH1> rtn = std::make_unique<TH1D>(
+      ps.get<std::string>("name", "").c_str(),
+      ps.get<std::string>("title", "").c_str(), xaxis[0], xaxis[1], xaxis[2]);
+  rtn->SetDirectory(nullptr);
+  return rtn;
+}
+
+template <typename HT>
+inline typename std::enable_if<HType_traits<HT>::NDim == 2,
+                               std::unique_ptr<TH2>>::type
+BuildHistFromFHiCL(fhicl::ParameterSet const &ps) {
+  std::array<double, 3> xaxis =
+      ps.get<std::array<double, 3>>("xaxis_descriptor");
+  std::array<double, 3> yaxis =
+      ps.get<std::array<double, 3>>("yaxis_descriptor");
+  std::unique_ptr<TH2> rtn =
+      std::make_unique<TH2D>(ps.get<std::string>("name", "").c_str(),
+                             ps.get<std::string>("title", "").c_str(), xaxis[0],
+                             xaxis[1], xaxis[2], yaxis[0], yaxis[1], yaxis[2]);
+  rtn->SetDirectory(nullptr);
+  return rtn;
+}
 } // namespace utility
 } // namespace nuis