diff --git a/cmake/setup.sh.in b/cmake/setup.sh.in
index 123aaef..44d8d3f 100644
--- a/cmake/setup.sh.in
+++ b/cmake/setup.sh.in
@@ -1,155 +1,159 @@
 # 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/sh
 
 ### Adapted from https://unix.stackexchange.com/questions/4965/keep-duplicates-out-of-path-on-source
 function add_to_PATH () {
   for d; do
 
     d=$(cd -- "$d" && { pwd -P || pwd; }) 2>/dev/null  # canonicalize symbolic links
     if [ -z "$d" ]; then continue; fi  # skip nonexistent directory
 
     if [ "$d" == "/usr/bin" ] || [ "$d" == "/usr/bin64" ] || [ "$d" == "/usr/local/bin" ] || [ "$d" == "/usr/local/bin64" ]; then
       case ":$PATH:" in
         *":$d:"*) :;;
         *) export PATH=$PATH:$d;;
       esac
     else
       case ":$PATH:" in
         *":$d:"*) :;;
         *) export PATH=$d:$PATH;;
       esac
     fi
   done
 }
 
 function add_to_LD_LIBRARY_PATH () {
   for d; do
 
     d=$(cd -- "$d" && { pwd -P || pwd; }) 2>/dev/null  # canonicalize symbolic links
     if [ -z "$d" ]; then continue; fi  # skip nonexistent directory
 
     if [ "$d" == "/usr/lib" ] || [ "$d" == "/usr/lib64" ] || [ "$d" == "/usr/local/lib" ] || [ "$d" == "/usr/local/lib64" ]; then
       case ":$LD_LIBRARY_PATH:" in
         *":$d:"*) :;;
         *) export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$d;;
       esac
     else
       case ":$LD_LIBRARY_PATH:" in
         *":$d:"*) :;;
         *) export LD_LIBRARY_PATH=$d:$LD_LIBRARY_PATH;;
       esac
     fi
   done
 }
 
 if [ "@EXTRA_SETUP_SCRIPT@" ]; then
   if [ ! -e @EXTRA_SETUP_SCRIPT@ ]; then
       echo "[WARN]: Extra setup script \"@EXTRA_SETUP_SCRIPT@\" requested, but could not be found. Skipping..."
   else
     echo "[INFO]: Sourcing extra setup from \"@EXTRA_SETUP_SCRIPT@\"."
     . @EXTRA_SETUP_SCRIPT@
   fi
 fi
 
 add_to_PATH "@CMAKE_INSTALL_PREFIX@/bin"
 
 add_to_LD_LIBRARY_PATH "@CMAKE_INSTALL_PREFIX@/lib"
 
 if [ ! "${ROOTSYS}" ]; then
   echo "[INFO]: Sourcing ROOT from: @CMAKE_ROOTSYS@"
   source "@CMAKE_ROOTSYS@/bin/thisroot.sh"
 fi
 
 if [ "@USE_T2K@" != "FALSE" ]; then
   echo "[INFO]: Adding T2K paths to the environment."
   export T2KREWEIGHT=@T2KREWEIGHT@
   add_to_LD_LIBRARY_PATH "@T2KREWEIGHT@/lib"
 fi
 
 if [ "@USE_NIWG@" != "FALSE" ]; then
   echo "[INFO]: Adding NIWG paths to the environment."
   export NIWG=@NIWG_ROOT@
   export NIWGREWEIGHT_INPUTS=@NIWG_ROOT@/inputs
   add_to_LD_LIBRARY_PATH "@NIWG_ROOT@"
 fi
 
 if [ "@USE_NEUT@" != "FALSE" ]; then
   echo "[INFO]: Adding NEUT library paths to the environment."
   export NEUT_ROOT=@NEUT_ROOT@
   export CERN=@CERN@
   export CERN_LEVEL=@CERN_LEVEL@
   add_to_LD_LIBRARY_PATH "${NEUT_LIB_DIR}" "${NEUT_ROOT}/src/reweight"
 fi
 
 if [ "@USE_NuWro@" != "FALSE" ]; then
   if [ "@NUWRO_BUILT_FROM_FILE@" == "FALSE" ]; then
     echo "[INFO]: Adding NuWro library paths to the environment."
     export NUWRO="@NUWRO@"
 
     add_to_PATH "@NUWRO@/bin"
     add_to_LD_LIBRARY_PATH "@NUWRO@/build/@CMAKE_SYSTEM_NAME@/lib"
 
     if [ "@NUWRO_INC@" ]; then
       export NUWRO_INC=@NUWRO_INC@
     fi
 
   else
     echo "[INFO]: NuWro support included from input event file."
   fi
 fi
 
 if [ "@NEED_PYTHIA6@" != "FALSE" ]; then
   echo "[INFO]: Adding PYTHIA6 library paths to the environment."
   export PYTHIA6="@PYTHIA6@"
 
   add_to_LD_LIBRARY_PATH "@PYTHIA6@"
 fi
 
 if [ "@USE_GENIE@" != "FALSE" ]; then
   echo "[INFO]: Adding GENIE paths to the environment."
 
   export GENIE="@GENIE@"
 
   export LHAPDF_LIB="@LHAPDF_LIB@"
   export LHAPDF_INC="@LHAPDF_INC@"
 
   export LIBXML2_LIB="@LIBXML2_LIB@"
   export LIBXML2_INC="@LIBXML2_INC@"
 
   export LOG4CPP_LIB="@LOG4CPP_LIB@"
   export LOG4CPP_INC="@LOG4CPP_INC@"
 
   if [ "@LHAPATH@" ]; then
     export LHAPATH="@LHAPATH@"
   fi
 
   add_to_PATH "@GENIE@/bin"
 
   add_to_LD_LIBRARY_PATH "@GENIE@/lib" "@LHAPDF_LIB@" "@LIBXML2_LIB@" "@LOG4CPP_LIB@"
 
 fi
 
 if [ "@BUILD_GiBUU@" != "FALSE" ]; then
   echo "[INFO]: Sourcing GiBUU tools."
   source @CMAKE_BINARY_DIR@/GiBUUTools/src/GiBUUTools-build/Linux/setup.sh
 fi
 
+if [ "@NUSYST_ROOT@" != "" ]; then
+  source @NUSYST_ROOT@/setup.sh
+fi
+
 export NUISANCE="@CMAKE_SOURCE_DIR@"
diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx
index 7516fee..9c33652 100644
--- a/src/InputHandler/FitEvent.cxx
+++ b/src/InputHandler/FitEvent.cxx
@@ -1,466 +1,466 @@
 // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is pddrt of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 #include "FitEvent.h"
 #include "TObjArray.h"
 #include <iostream>
 
 FitEvent::FitEvent() {
   fGenInfo = NULL;
   kRemoveFSIParticles = true;
   kRemoveUndefParticles = true;
 
-  AllocateParticleStack(25);
+  AllocateParticleStack(50);
 };
 
 void FitEvent::AddGeneratorInfo(GeneratorInfoBase *gen) {
   fGenInfo = gen;
   gen->AllocateParticleStack(kMaxParticles);
 }
 
 void FitEvent::AllocateParticleStack(int stacksize) {
 
   kMaxParticles = stacksize;
 
   if (!stacksize) {
     return;
   }
 
   LOG(DEB) << "Allocating particle stack of size: " << stacksize << std::endl;
 
   fParticleList = new FitParticle *[kMaxParticles];
 
   fParticleMom = new double *[kMaxParticles];
   fParticleState = new UInt_t[kMaxParticles];
   fParticlePDG = new int[kMaxParticles];
 
   fOrigParticleMom = new double *[kMaxParticles];
   fOrigParticleState = new UInt_t[kMaxParticles];
   fOrigParticlePDG = new int[kMaxParticles];
 
   for (size_t i = 0; i < kMaxParticles; i++) {
     fParticleList[i] = NULL;
     fParticleMom[i] = new double[4];
     fOrigParticleMom[i] = new double[4];
   }
 
   if (fGenInfo)
     fGenInfo->AllocateParticleStack(kMaxParticles);
 }
 
 void FitEvent::ExpandParticleStack(int stacksize) {
   DeallocateParticleStack();
   AllocateParticleStack(stacksize);
 }
 
 void FitEvent::DeallocateParticleStack() {
 
   if (!kMaxParticles) {
     return;
   }
 
   for (size_t i = 0; i < kMaxParticles; i++) {
     if (fParticleList[i])
       delete fParticleList[i];
     delete[] fParticleMom[i];
     delete[] fOrigParticleMom[i];
   }
   delete[] fParticleMom;
   delete[] fOrigParticleMom;
 
   delete[] fParticleList;
 
   delete[] fParticleState;
   delete[] fParticlePDG;
 
   delete[] fOrigParticleState;
   delete[] fOrigParticlePDG;
 
   if (fGenInfo)
     fGenInfo->DeallocateParticleStack();
 
   kMaxParticles = 0;
 
   fNParticles = 0;
 }
 
 void FitEvent::ClearFitParticles() {
   for (size_t i = 0; i < kMaxParticles; i++) {
     fParticleList[i] = NULL;
   }
 }
 
 void FitEvent::FreeFitParticles() {
   for (size_t i = 0; i < kMaxParticles; i++) {
     FitParticle *fp = fParticleList[i];
     if (fp)
       delete fp;
     fParticleList[i] = NULL;
   }
 }
 
 void FitEvent::ResetParticleList() {
   for (unsigned int i = 0; i < kMaxParticles; i++) {
     FitParticle *fp = fParticleList[i];
     if (fp)
       delete fp;
     fParticleList[i] = NULL;
   }
 }
 
 void FitEvent::HardReset() {
   for (unsigned int i = 0; i < kMaxParticles; i++) {
     fParticleList[i] = NULL;
   }
 }
 
 void FitEvent::ResetEvent() {
   Mode = 9999;
   fEventNo = -1;
   fTotCrs = -1.0;
   fTargetA = -1;
   fTargetZ = -1;
   fTargetH = -1;
   fBound = false;
   fNParticles = 0;
 
   if (fGenInfo)
     fGenInfo->Reset();
 
   for (unsigned int i = 0; i < kMaxParticles; i++) {
     if (fParticleList[i])
       delete fParticleList[i];
     fParticleList[i] = NULL;
 
     continue;
 
     fParticlePDG[i] = 0;
     fParticleState[i] = kUndefinedState;
     fParticleMom[i][0] = 0.0;
     fParticleMom[i][1] = 0.0;
     fParticleMom[i][2] = 0.0;
     fParticleMom[i][3] = 0.0;
 
     fOrigParticlePDG[i] = 0;
     fOrigParticleState[i] = kUndefinedState;
     fOrigParticleMom[i][0] = 0.0;
     fOrigParticleMom[i][1] = 0.0;
     fOrigParticleMom[i][2] = 0.0;
     fOrigParticleMom[i][3] = 0.0;
   }
 }
 
 void FitEvent::OrderStack() {
   // Copy current stack
   int npart = fNParticles;
 
   for (int i = 0; i < npart; i++) {
     fOrigParticlePDG[i] = fParticlePDG[i];
     fOrigParticleState[i] = fParticleState[i];
     fOrigParticleMom[i][0] = fParticleMom[i][0];
     fOrigParticleMom[i][1] = fParticleMom[i][1];
     fOrigParticleMom[i][2] = fParticleMom[i][2];
     fOrigParticleMom[i][3] = fParticleMom[i][3];
   }
 
   // Now run loops for each particle
   fNParticles = 0;
   int stateorder[6] = {kInitialState,   kFinalState,     kFSIState,
                        kNuclearInitial, kNuclearRemnant, kUndefinedState};
 
   for (int s = 0; s < 6; s++) {
     for (int i = 0; i < npart; i++) {
       if ((UInt_t)fOrigParticleState[i] != (UInt_t)stateorder[s])
         continue;
 
       fParticlePDG[fNParticles] = fOrigParticlePDG[i];
       fParticleState[fNParticles] = fOrigParticleState[i];
       fParticleMom[fNParticles][0] = fOrigParticleMom[i][0];
       fParticleMom[fNParticles][1] = fOrigParticleMom[i][1];
       fParticleMom[fNParticles][2] = fOrigParticleMom[i][2];
       fParticleMom[fNParticles][3] = fOrigParticleMom[i][3];
 
       fNParticles++;
     }
   }
 
   if (LOG_LEVEL(DEB)) {
     LOG(DEB) << "Ordered stack" << std::endl;
     for (int i = 0; i < fNParticles; i++) {
       LOG(DEB) << "Particle " << i << ". " << fParticlePDG[i] << " "
                << fParticleMom[i][0] << " " << fParticleMom[i][1] << " "
                << fParticleMom[i][2] << " " << fParticleMom[i][3] << " "
                << fParticleState[i] << std::endl;
     }
   }
 
   if (fNParticles != npart) {
     ERR(FTL) << "Dropped some particles when ordering the stack!" << std::endl;
   }
 
   return;
 }
 
 void FitEvent::Print() {
   if (LOG_LEVEL(FIT)) {
     LOG(FIT) << "FITEvent print" << std::endl;
     LOG(FIT) << "Mode: " << Mode << ", Weight: " << InputWeight << std::endl;
     LOG(FIT) << "Particles: " << fNParticles << std::endl;
     LOG(FIT) << " -> Particle Stack " << std::endl;
     for (int i = 0; i < fNParticles; i++) {
       LOG(FIT) << " -> -> " << i << ". " << fParticlePDG[i] << " "
                << fParticleState[i] << " "
                << "  Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1]
                << ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3]
                << ")." << std::endl;
     }
   }
   return;
 }
 
 /* Read/Write own event class */
 void FitEvent::SetBranchAddress(TChain *tn) {
   tn->SetBranchAddress("Mode", &Mode);
 
   tn->SetBranchAddress("EventNo", &fEventNo);
   tn->SetBranchAddress("TotCrs", &fTotCrs);
   tn->SetBranchAddress("TargetA", &fTargetA);
   tn->SetBranchAddress("TargetH", &fTargetH);
   tn->SetBranchAddress("Bound", &fBound);
 
   tn->SetBranchAddress("RWWeight", &SavedRWWeight);
   tn->SetBranchAddress("InputWeight", &InputWeight);
 }
 
 void FitEvent::AddBranchesToTree(TTree *tn) {
   tn->Branch("Mode", &Mode, "Mode/I");
 
   tn->Branch("EventNo", &fEventNo, "EventNo/i");
   tn->Branch("TotCrs", &fTotCrs, "TotCrs/D");
   tn->Branch("TargetA", &fTargetA, "TargetA/I");
   tn->Branch("TargetH", &fTargetH, "TargetH/I");
   tn->Branch("Bound", &fBound, "Bound/O");
 
   tn->Branch("RWWeight", &RWWeight, "RWWeight/D");
   tn->Branch("InputWeight", &InputWeight, "InputWeight/D");
 
   tn->Branch("NParticles", &fNParticles, "NParticles/I");
   tn->Branch("ParticleState", fOrigParticleState,
              "ParticleState[NParticles]/i");
   tn->Branch("ParticlePDG", fOrigParticlePDG, "ParticlePDG[NParticles]/I");
   tn->Branch("ParticleMom", fOrigParticleMom, "ParticleMom[NParticles][4]/D");
 }
 
 // ------- EVENT ACCESS FUNCTION --------- //
 TLorentzVector FitEvent::GetParticleP4(int index) const {
   if (index == -1 or index >= fNParticles)
     return TLorentzVector();
   return TLorentzVector(fParticleMom[index][0], fParticleMom[index][1],
                         fParticleMom[index][2], fParticleMom[index][3]);
 }
 
 TVector3 FitEvent::GetParticleP3(int index) const {
   if (index == -1 or index >= fNParticles)
     return TVector3();
   return TVector3(fParticleMom[index][0], fParticleMom[index][1],
                   fParticleMom[index][2]);
 }
 
 double FitEvent::GetParticleMom(int index) const {
   if (index == -1 or index >= fNParticles)
     return 0.0;
   return sqrt(fParticleMom[index][0] * fParticleMom[index][0] +
               fParticleMom[index][1] * fParticleMom[index][1] +
               fParticleMom[index][2] * fParticleMom[index][2]);
 }
 
 double FitEvent::GetParticleMom2(int index) const {
   if (index == -1 or index >= fNParticles)
     return 0.0;
   return fabs((fParticleMom[index][0] * fParticleMom[index][0] +
                fParticleMom[index][1] * fParticleMom[index][1] +
                fParticleMom[index][2] * fParticleMom[index][2]));
 }
 
 double FitEvent::GetParticleE(int index) const {
   if (index == -1 or index >= fNParticles)
     return 0.0;
   return fParticleMom[index][3];
 }
 
 int FitEvent::GetParticleState(int index) const {
   if (index == -1 or index >= fNParticles)
     return kUndefinedState;
   return (fParticleState[index]);
 }
 
 int FitEvent::GetParticlePDG(int index) const {
   if (index == -1 or index >= fNParticles)
     return 0;
   return (fParticlePDG[index]);
 }
 
 FitParticle *FitEvent::GetParticle(int const i) {
   // Check Valid Index
   if (i == -1) {
     return NULL;
   }
 
   // Check Valid
   if (i > fNParticles) {
     ERR(FTL) << "Requesting particle beyond stack!" << std::endl
              << "i = " << i << " N = " << fNParticles << std::endl
              << "Mode = " << Mode << std::endl;
     throw;
   }
 
   if (!fParticleList[i]) {
     /*
     std::cout << "Creating particle with values i " << i << " ";
     std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] <<  " " <<
     fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
     std::cout << fParticlePDG[i] << " " << fParticleState[i] << std::endl;
     */
     fParticleList[i] = new FitParticle(fParticleMom[i][0], fParticleMom[i][1],
                                        fParticleMom[i][2], fParticleMom[i][3],
                                        fParticlePDG[i], fParticleState[i]);
   } else {
     /*
     std::cout << "Filling particle with values i " << i << " ";
     std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] <<  " " <<
     fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
     std::cout << fParticlePDG[i] << " "<< fParticleState[i] <<std::endl;
     */
     fParticleList[i]->SetValues(fParticleMom[i][0], fParticleMom[i][1],
                                 fParticleMom[i][2], fParticleMom[i][3],
                                 fParticlePDG[i], fParticleState[i]);
   }
 
   return fParticleList[i];
 }
 
 bool FitEvent::HasParticle(int const pdg, int const state) const {
   bool found = false;
   for (int i = 0; i < fNParticles; i++) {
     if (state != -1 && fParticleState[i] != (uint)state)
       continue;
     if (fParticlePDG[i] == pdg)
       found = true;
   }
   return found;
 }
 
 int FitEvent::NumParticle(int const pdg, int const state) const {
   int nfound = 0;
   for (int i = 0; i < fNParticles; i++) {
     if (state != -1 and fParticleState[i] != (uint)state)
       continue;
     if (pdg == 0 or fParticlePDG[i] == pdg)
       nfound += 1;
   }
   return nfound;
 }
 
 std::vector<int> FitEvent::GetAllParticleIndices(int const pdg,
                                                  int const state) const {
   std::vector<int> indexlist;
   for (int i = 0; i < fNParticles; i++) {
     if (state != -1 and fParticleState[i] != (uint)state)
       continue;
     if (pdg == 0 or fParticlePDG[i] == pdg) {
       indexlist.push_back(i);
     }
   }
   return indexlist;
 }
 
 std::vector<FitParticle *> FitEvent::GetAllParticle(int const pdg,
                                                     int const state) {
   std::vector<int> indexlist = GetAllParticleIndices(pdg, state);
   std::vector<FitParticle *> plist;
   for (std::vector<int>::iterator iter = indexlist.begin();
        iter != indexlist.end(); iter++) {
     plist.push_back(GetParticle((*iter)));
   }
   return plist;
 }
 
 int FitEvent::GetHMParticleIndex(int const pdg, int const state) const {
   double maxmom2 = -9999999.9;
   int maxind = -1;
   for (int i = 0; i < fNParticles; i++) {
     if (state != -1 and fParticleState[i] != (uint)state)
       continue;
     if (pdg == 0 or fParticlePDG[i] == pdg) {
       double newmom2 = GetParticleMom2(i);
       if (newmom2 > maxmom2) {
         maxind = i;
         maxmom2 = newmom2;
       }
     }
   }
 
   return maxind;
 }
 
 int FitEvent::GetBeamNeutrinoIndex(void) const {
   for (int i = 0; i < fNParticles; i++) {
     if (fParticleState[i] != kInitialState)
       continue;
     int pdg = abs(fParticlePDG[i]);
     if (pdg == 12 or pdg == 14 or pdg == 16) {
       return i;
     }
   }
   return 0;
 }
 
 int FitEvent::GetBeamElectronIndex(void) const {
   return GetHMISParticleIndex(11);
 }
 
 int FitEvent::GetBeamPionIndex(void) const {
   return GetHMISParticleIndex(PhysConst::pdg_pions);
 }
 
 int FitEvent::NumFSMesons() {
   int nMesons = 0;
 
   for (int i = 0; i < fNParticles; i++) {
     if (fParticleState[i] != kFinalState)
       continue;
     if (abs(fParticlePDG[i]) >= 111 && abs(fParticlePDG[i]) <= 557)
       nMesons += 1;
   }
 
   return nMesons;
 }
 
 int FitEvent::NumFSLeptons(void) const {
   int nLeptons = 0;
 
   for (int i = 0; i < fNParticles; i++) {
     if (fParticleState[i] != kFinalState)
       continue;
     if (abs(fParticlePDG[i]) == 11 || abs(fParticlePDG[i]) == 13 ||
         abs(fParticlePDG[i]) == 15)
       nLeptons += 1;
   }
 
   return nLeptons;
 }
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 278295d..aa9a02a 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,540 +1,543 @@
 // 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/>.
  *******************************************************************************/
 #ifdef __GENIE_ENABLED__
 #include "GENIEInputHandler.h"
 
 #include "InputUtils.h"
 
 #ifdef __DUNERWT_ENABLED__
 
 #include "systematicstools/utility/ParameterAndProviderConfigurationUtility.hh"
 
 #include "fhiclcpp/make_ParameterSet.h"
 
 #endif
 
 GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
 
 void GENIEGeneratorInfo::AddBranchesToTree(TTree *tn) {
   tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
 }
 
 void GENIEGeneratorInfo::SetBranchesFromTree(TTree *tn) {
   tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
 }
 
 void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
   fGenieParticlePDGs = new int[stacksize];
 }
 
 void GENIEGeneratorInfo::DeallocateParticleStack() {
   delete fGenieParticlePDGs;
 }
 
 void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord *ntpl) {
   Reset();
 
   // Check for GENIE Event
   if (!ntpl)
     return;
   if (!ntpl->event)
     return;
 
   // Cast Event Record
   GHepRecord *ghep = static_cast<GHepRecord *>(ntpl->event);
   if (!ghep)
     return;
 
   // Fill Particle Stack
   GHepParticle *p = 0;
   TObjArrayIter iter(ghep);
 
   // Loop over all particles
   int i = 0;
   while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
     if (!p)
       continue;
 
     // Get PDG
     fGenieParticlePDGs[i] = p->Pdg();
     i++;
   }
 }
 
 void GENIEGeneratorInfo::Reset() {
   for (int i = 0; i < kMaxParticles; i++) {
     fGenieParticlePDGs[i] = 0;
   }
 }
 
 GENIEInputHandler::GENIEInputHandler(std::string const &handle,
                                      std::string const &rawinputs) {
   LOG(SAM) << "Creating GENIEInputHandler : " << handle << std::endl;
 
   // Run a joint input handling
   fName = handle;
 
   // Setup the TChain
   fGENIETree = new TChain("gtree");
   fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
   fCacheSize = FitPar::Config().GetParI("CacheSize");
   fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
 
   // Loop over all inputs and grab flux, eventhist, and nevents
   std::vector<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
   for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
     // Open File for histogram access
     TFile *inp_file = new TFile(
         InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
     if (!inp_file or inp_file->IsZombie()) {
       THROW("GENIE File IsZombie() at : '"
             << inputs[inp_it] << "'" << std::endl
             << "Check that your file paths are correct and the file exists!"
             << std::endl
             << "$ ls -lh " << inputs[inp_it]);
     }
 
     // Get Flux/Event hist
     TH1D *fluxhist = (TH1D *)inp_file->Get("nuisance_flux");
     TH1D *eventhist = (TH1D *)inp_file->Get("nuisance_events");
     if (!fluxhist or !eventhist) {
       ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
       inp_file->ls();
       THROW("GENIE FILE doesn't contain flux/xsec info."
             << std::endl
             << "Try running the app PrepareGENIE first on :" << inputs[inp_it]
             << std::endl
             << "$ PrepareGENIE -h");
     }
 
     // Get N Events
     TTree *genietree = (TTree *)inp_file->Get("gtree");
     if (!genietree) {
       ERROR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
       THROW("Check your inputs, they may need to be completely regenerated!");
       throw;
     }
     int nevents = genietree->GetEntries();
     if (nevents <= 0) {
       THROW("Trying to a TTree with "
             << nevents << " to TChain from : " << inputs[inp_it]);
     }
 
     // Register input to form flux/event rate hists
     RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
 
     // Add To TChain
     fGENIETree->AddFile(inputs[inp_it].c_str());
   }
 
   // Registor all our file inputs
   SetupJointInputs();
 
   // Assign to tree
   fEventType = kGENIE;
   fGenieNtpl = NULL;
   fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
 
   // Libraries should be seen but not heard...
   StopTalking();
   fGENIETree->GetEntry(0);
   StartTalking();
 
 #ifndef __DUNERWT_ENABLED__
   // Create Fit Event
   fNUISANCEEvent = new FitEvent();
   fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
 
   if (fSaveExtra) {
     fGenieInfo = new GENIEGeneratorInfo();
     fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
   }
 
   fNUISANCEEvent->HardReset();
 #else
   std::vector<nuiskey> HandlerOpts = Config::QueryKeys("GENIEInputHandler");
   fUseCache = HandlerOpts.size() && HandlerOpts.front().Has("UseCache") &&
               HandlerOpts.front().GetB("UseCache");
 
   DUNERwtCachedResponseReader = nullptr;
   HaveCachedResponseReader = false;
   if (fUseCache && (inputs.size() == 1)) {
     std::vector<nuiskey> DuneRwtCacheParams =
         Config::QueryKeys("DUNERwtResponseCache");
     for (nuiskey &key : DuneRwtCacheParams) {
       if (key.Has("Input") && (key.GetS("Input") == inputs.front()) &&
           key.Has("CacheFile") && key.Has("ParameterFHiCL")) {
 
         fhicl::ParameterSet ps =
             fhicl::make_ParameterSet(key.GetS("ParameterFHiCL"));
 
         fhicl::ParameterSet syst_providers = ps.get<fhicl::ParameterSet>(
             "generated_systematic_provider_configuration");
 
         systtools::param_header_map_t configuredParameterHeaders =
             systtools::BuildParameterHeaders(syst_providers);
 
         DUNERwtCachedResponseReader =
             std::make_unique<systtools::PrecalculatedResponseReader<5>>(
                 key.GetS("CacheFile"), "resp_tree",
                 configuredParameterHeaders.size());
         HaveCachedResponseReader = true;
         break;
       }
     }
   } else {
     fNUISANCEEvent = new FitEvent();
     fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
 
     if (fSaveExtra) {
       fGenieInfo = new GENIEGeneratorInfo();
       fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
     }
 
     fNUISANCEEvent->HardReset();
   }
 #endif
 };
 
 GENIEInputHandler::~GENIEInputHandler() {
   // if (fGenieGHep) delete fGenieGHep;
   // if (fGenieNtpl) delete fGenieNtpl;
   // if (fGENIETree) delete fGENIETree;
   // if (fGenieInfo) delete fGenieInfo;
 }
 
 void GENIEInputHandler::CreateCache() {
   if (fCacheSize > 0) {
     // fGENIETree->SetCacheEntryRange(0, fNEvents);
     fGENIETree->AddBranchToCache("*", 1);
     fGENIETree->SetCacheSize(fCacheSize);
   }
 }
 
 void GENIEInputHandler::RemoveCache() {
   // fGENIETree->SetCacheEntryRange(0, fNEvents);
   fGENIETree->AddBranchToCache("*", 0);
   fGENIETree->SetCacheSize(0);
 }
 
 FitEvent *GENIEInputHandler::GetNuisanceEvent(const UInt_t entry,
                                               const bool lightweight) {
   if (entry >= (UInt_t)fNEvents)
     return NULL;
 
   // Reduce memory pressure from the cache by clearing out the last entry each
   // time.
   if (entry && rwEvs[entry - 1].NParticles()) {
     rwEvs[entry - 1].DeallocateParticleStack();
   }
 
   // Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
   fGENIETree->GetEntry(entry);
 
 #ifdef __DUNERWT_ENABLED__
   if (entry >= rwEvs.size()) {
     rwEvs.push_back(FitEvent());
     if (HaveCachedResponseReader) {
       rwEvs.back().DUNERwtPolyResponses =
           DUNERwtCachedResponseReader->GetEventResponse(entry);
       rwEvs.back().HasDUNERwtPolyResponses = true;
     }
   }
   rwEvs[entry].SetGenieEvent(fGenieNtpl);
 
   fNUISANCEEvent = &rwEvs[entry];
 #endif
 
   // Run NUISANCE Vector Filler
   if (!lightweight) {
     CalcNUISANCEKinematics();
   }
 #ifdef __PROB3PP_ENABLED__
   else {
     // Check for GENIE Event
     if (!fGenieNtpl)
       return NULL;
     if (!fGenieNtpl->event)
       return NULL;
 
     // Cast Event Record
     fGenieGHep = fGenieNtpl->event;
     if (!fGenieGHep)
       return NULL;
 
     TObjArrayIter iter(fGenieGHep);
     genie::GHepParticle *p;
     while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
       if (!p) {
         continue;
       }
 
       // Get Status
       int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
       if (state != genie::kIStInitialState) {
         continue;
       }
       fNUISANCEEvent->probe_E = p->E() * 1.E3;
       fNUISANCEEvent->probe_pdg = p->Pdg();
       break;
     }
   }
 #endif
 
   // Setup Input scaling for joint inputs
   fNUISANCEEvent->InputWeight = GetInputWeight(entry);
 
   return fNUISANCEEvent;
 }
 
 int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle *p,
                                               int mode) {
   /*
     kIStUndefined                  = -1,
     kIStInitialState               =  0,   / generator-level initial state /
     kIStStableFinalState           =  1,   / generator-level final state:
     particles to be tracked by detector-level MC /
     kIStIntermediateState          =  2,
     kIStDecayedState               =  3,
     kIStCorrelatedNucleon          = 10,
     kIStNucleonTarget              = 11,
     kIStDISPreFragmHadronicState   = 12,
     kIStPreDecayResonantState      = 13,
     kIStHadronInTheNucleus         = 14,   / hadrons inside the nucleus:
     marked for hadron transport modules to act on /
     kIStFinalStateNuclearRemnant   = 15,   / low energy nuclear fragments
     entering the record collectively as a 'hadronic blob' pseudo-particle /
     kIStNucleonClusterTarget       = 16,   // for composite nucleons before
     phase space decay
   */
 
   int state = kUndefinedState;
   switch (p->Status()) {
   case genie::kIStNucleonTarget:
   case genie::kIStInitialState:
   case genie::kIStCorrelatedNucleon:
   case genie::kIStNucleonClusterTarget:
     state = kInitialState;
     break;
 
   case genie::kIStStableFinalState:
     state = kFinalState;
     break;
 
   case genie::kIStHadronInTheNucleus:
     if (abs(mode) == 2)
       state = kInitialState;
     else
       state = kFSIState;
     break;
 
   case genie::kIStPreDecayResonantState:
   case genie::kIStDISPreFragmHadronicState:
   case genie::kIStIntermediateState:
     state = kFSIState;
     break;
 
   case genie::kIStFinalStateNuclearRemnant:
   case genie::kIStUndefined:
   case genie::kIStDecayedState:
   default:
     break;
   }
 
   // Flag to remove nuclear part in genie
   if (p->Pdg() > 1000000) {
     if (state == kInitialState)
       state = kNuclearInitial;
     else if (state == kFinalState)
       state = kNuclearRemnant;
   }
 
   return state;
 }
 #endif
 
 #ifdef __GENIE_ENABLED__
 int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord *gheprec) {
   // Electron Scattering
   if (gheprec->Summary()->ProcInfo().IsEM()) {
     if (gheprec->Summary()->InitState().ProbePdg() == 11) {
       if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
         return 1;
       else if (gheprec->Summary()->ProcInfo().IsMEC())
         return 2;
       else if (gheprec->Summary()->ProcInfo().IsResonant())
         return 13;
       else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
         return 26;
       else {
         ERROR(WRN,
               "Unknown GENIE Electron Scattering Mode!"
                   << std::endl
                   << "ScatteringTypeId = "
                   << gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
                   << "InteractionTypeId = "
                   << gheprec->Summary()->ProcInfo().InteractionTypeId()
                   << std::endl
                   << genie::ScatteringType::AsString(
                          gheprec->Summary()->ProcInfo().ScatteringTypeId())
                   << " "
                   << genie::InteractionType::AsString(
                          gheprec->Summary()->ProcInfo().InteractionTypeId())
                   << " " << gheprec->Summary()->ProcInfo().IsMEC());
         return 0;
       }
     }
 
     // Weak CC
   } else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
     // CC MEC
     if (gheprec->Summary()->ProcInfo().IsMEC()) {
       if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
         return 2;
       else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
         return -2;
 
       // CC OTHER
     } else {
       return utils::ghep::NeutReactionCode(gheprec);
     }
 
     // Weak NC
   } else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
     // NC MEC
     if (gheprec->Summary()->ProcInfo().IsMEC()) {
       if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
         return 32;
       else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
         return -32;
 
       // NC OTHER
     } else {
       return utils::ghep::NeutReactionCode(gheprec);
     }
   }
 
   return 0;
 }
 
 void GENIEInputHandler::CalcNUISANCEKinematics() {
   // Reset all variables
   fNUISANCEEvent->ResetEvent();
 
   // Check for GENIE Event
   if (!fGenieNtpl)
     return;
   if (!fGenieNtpl->event)
     return;
 
   // Cast Event Record
   fGenieGHep = fGenieNtpl->event;
   if (!fGenieGHep)
     return;
 
   // Convert GENIE Reaction Code
   fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
 
   // Set Event Info
   fNUISANCEEvent->fEventNo = 0.0;
   fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
   fNUISANCEEvent->fTargetA = 0.0;
   fNUISANCEEvent->fTargetZ = 0.0;
   fNUISANCEEvent->fTargetH = 0;
   fNUISANCEEvent->fBound = 0.0;
   fNUISANCEEvent->InputWeight =
       1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
 
   // Get N Particle Stack
   unsigned int npart = fGenieGHep->GetEntries();
   unsigned int kmax = fNUISANCEEvent->kMaxParticles;
   if (npart > kmax) {
     fNUISANCEEvent->ExpandParticleStack(npart);
   }
 
   // Fill Particle Stack
   GHepParticle *p = 0;
   TObjArrayIter iter(fGenieGHep);
   fNUISANCEEvent->fNParticles = 0;
 
   // Loop over all particles
   while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
     if (!p)
       continue;
 
     // Get Status
     int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
 
     // Remove Undefined
     if (kRemoveUndefParticles && state == kUndefinedState)
       continue;
 
     // Remove FSI
     if (kRemoveFSIParticles && state == kFSIState)
       continue;
 
     if (kRemoveNuclearParticles &&
         (state == kNuclearInitial || state == kNuclearRemnant))
       continue;
 
     // Fill Vectors
     int curpart = fNUISANCEEvent->fNParticles;
     fNUISANCEEvent->fParticleState[curpart] = state;
 
     // Mom
     fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
     fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
     fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
     fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
 
     // PDG
     fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
 
     // Add to N particle count
     fNUISANCEEvent->fNParticles++;
 
     // Extra Check incase GENIE fails.
     if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
-      ERR(WRN) << "Number of GENIE Particles exceeds maximum!" << std::endl;
+      ERR(WRN) << "Number of GENIE Particles exceeds maximum (Max: " << kmax
+               << ", GHEP: " << fGenieGHep->GetEntries()
+               << ", Added: " << fNUISANCEEvent->fNParticles << ")!"
+               << std::endl;
       ERR(WRN) << "Extend kMax, or run without including FSI particles!"
                << std::endl;
       break;
     }
   }
 
   // Fill Extra Stack
   if (fSaveExtra)
     fGenieInfo->FillGeneratorInfo(fGenieNtpl);
 
   // Run Initial, FSI, Final, Other ordering.
   fNUISANCEEvent->OrderStack();
 
   FitParticle *ISNeutralLepton =
       fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
   if (ISNeutralLepton) {
     fNUISANCEEvent->probe_E = ISNeutralLepton->E();
     fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
   }
 
   return;
 }
 
 void GENIEInputHandler::Print() {}
 
 #endif
diff --git a/src/Reweight/DUNERwtWeightEngine.cxx b/src/Reweight/DUNERwtWeightEngine.cxx
index 3477010..00f5d3d 100644
--- a/src/Reweight/DUNERwtWeightEngine.cxx
+++ b/src/Reweight/DUNERwtWeightEngine.cxx
@@ -1,167 +1,164 @@
 // 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/>.
  *******************************************************************************/
 
 #include "DUNERwtWeightEngine.h"
 
 #include <limits>
 #include <string>
 
 DUNERwtWeightEngine::DUNERwtWeightEngine() { Config(); }
 
 void DUNERwtWeightEngine::Config() {
   std::vector<nuiskey> DuneRwtParam = Config::QueryKeys("DUNERwt");
 
   if (DuneRwtParam.size() < 1) {
     ERROR(WRN, "Oscillation parameters specified but no OscParam element "
                "configuring the experimental characteristics found.\nExpect at "
                "least <OscParam baseline_km=\"XXX\" />. Pausing for "
                "10...");
     sleep(10);
     return;
   }
 
   std::string fhicl_name = DuneRwtParam.front().GetS("ConfigFHiCL");
 
   DUNErwt.LoadConfiguration(fhicl_name);
 }
 
 int DUNERwtWeightEngine::ConvDial(std::string name) {
   if (!DUNErwt.HaveHeader(name)) {
     THROW("DUNERwtWeightEngine passed dial: "
           << name << " that it does not understand.");
   }
   return DUNErwt.GetHeaderId(name);
 }
 
 void DUNERwtWeightEngine::IncludeDial(std::string name, double startval) {
   EnabledParams.push_back({systtools::paramId_t(ConvDial(name)), startval});
 }
 
 void DUNERwtWeightEngine::SetDialValue(int nuisenum, double val) {
 
   systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
   systtools::ParamValue &pval =
       GetParamElementFromContainer(EnabledParams, DuneRwtEnum);
   fHasChanged = (pval.val - val) > std::numeric_limits<double>::epsilon();
   pval.val = val;
 }
 void DUNERwtWeightEngine::SetDialValue(std::string name, double val) {
   if (!IsDialIncluded(name)) {
     THROW("DUNERwtWeightEngine passed dial: " << name
                                               << " that is not enabled.");
   }
 
   systtools::ParamValue &pval =
       GetParamElementFromContainer(EnabledParams, ConvDial(name));
   fHasChanged = (pval.val - val) > std::numeric_limits<double>::epsilon();
   pval.val = val;
 }
 
 bool DUNERwtWeightEngine::IsDialIncluded(std::string name) {
   return IsDialIncluded(ConvDial(name));
 }
 bool DUNERwtWeightEngine::IsDialIncluded(int nuisenum) {
   systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
   return systtools::ContainterHasParam(EnabledParams, DuneRwtEnum);
 }
 
 double DUNERwtWeightEngine::GetDialValue(std::string name) {
   if (!IsDialIncluded(name)) {
     THROW("DUNERwtWeightEngine passed dial: " << name
                                               << " that is not enabled.");
   }
   systtools::ParamValue &pval =
       GetParamElementFromContainer(EnabledParams, ConvDial(name));
   return pval.val;
 }
 double DUNERwtWeightEngine::GetDialValue(int nuisenum) {
   if (!IsDialIncluded(nuisenum)) {
     THROW("DUNERwtWeightEngine passed dial: " << nuisenum
                                               << " that is not enabled.");
   }
   systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
   systtools::ParamValue &pval =
       GetParamElementFromContainer(EnabledParams, DuneRwtEnum);
   return pval.val;
 }
 
 void DUNERwtWeightEngine::Reconfigure(bool silent) { fHasChanged = false; };
 
 bool DUNERwtWeightEngine::NeedsEventReWeight() {
   if (fHasChanged) {
     return true;
   }
   return false;
 }
 
 double DUNERwtWeightEngine::CalcWeight(BaseFitEvt *evt) {
 
   double weight = 1;
 
   if (evt->HasDUNERwtPolyResponses) {
     for (size_t i = 0; i < EnabledParams.size(); ++i) {
       if (DUNErwt.IsSplineParam(EnabledParams[i].pid)) {
         if (!ContainterHasParam(evt->DUNERwtPolyResponses,
                                 EnabledParams[i].pid)) {
           continue;
         }
         weight *= GetParamElementFromContainer(evt->DUNERwtPolyResponses,
                                                EnabledParams[i].pid)
                       .resp.eval(EnabledParams[i].val);
 
-        evt->DUNERwtResponses =
-            DUNErwt.GetEventResponses(*evt->genie_event->event);
-
       } else {
         if (!evt->HasDUNERwtResponses) {
           evt->DUNERwtResponses =
               DUNErwt.GetEventResponses(*evt->genie_event->event);
           evt->HasDUNERwtResponses = true;
         }
         weight *= DUNErwt.GetDiscreteResponse(EnabledParams[i].pid,
                                               int(EnabledParams[i].val),
                                               evt->DUNERwtResponses);
       }
     }
   } else {
     if (!evt->HasDUNERwtResponses) {
       evt->DUNERwtResponses =
           DUNErwt.GetEventResponses(*evt->genie_event->event);
       evt->HasDUNERwtResponses = true;
     }
 
     for (size_t i = 0; i < EnabledParams.size(); ++i) {
       if (DUNErwt.IsSplineParam(EnabledParams[i].pid)) {
         weight *= DUNErwt.GetParameterResponse(
             EnabledParams[i].pid, EnabledParams[i].val, evt->DUNERwtResponses);
 
       } else {
         weight *= DUNErwt.GetDiscreteResponse(EnabledParams[i].pid,
                                               int(EnabledParams[i].val),
                                               evt->DUNERwtResponses);
       }
     }
   }
 
   return weight;
 }
 
 void DUNERwtWeightEngine::Print() {
   std::cout << "DUNERwtWeightEngine: " << std::endl;
 }