diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx
index eabd34b..fe038b5 100644
--- a/src/InputHandler/FitEvent.cxx
+++ b/src/InputHandler/FitEvent.cxx
@@ -1,493 +1,500 @@
 // 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(400);
 };
 
 void FitEvent::AddGeneratorInfo(GeneratorInfoBase *gen) {
   fGenInfo = gen;
   gen->AllocateParticleStack(kMaxParticles);
 }
 
 void FitEvent::AllocateParticleStack(int stacksize) {
   NUIS_LOG(DEB, "Allocating particle stack of size: " << stacksize);
   kMaxParticles = stacksize;
 
   fParticleList = new FitParticle *[kMaxParticles];
 
   fParticleMom = new double *[kMaxParticles];
   fParticleState = new UInt_t[kMaxParticles];
   fParticlePDG = new int[kMaxParticles];
   fPrimaryVertex = new bool[kMaxParticles];
 
   fOrigParticleMom = new double *[kMaxParticles];
   fOrigParticleState = new UInt_t[kMaxParticles];
   fOrigParticlePDG = new int[kMaxParticles];
   fOrigPrimaryVertex = new bool[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() {
   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 fPrimaryVertex;
 
   delete fOrigParticleState;
   delete fOrigParticlePDG;
   delete fOrigPrimaryVertex;
 
   if (fGenInfo)
     fGenInfo->DeallocateParticleStack();
 
   kMaxParticles = 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;
     fPrimaryVertex[i] = false;
 
     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;
     fOrigPrimaryVertex[i] = false;
   }
 }
 
 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];
     fOrigPrimaryVertex[i] = fPrimaryVertex[i];
   }
 
   // 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];
       fPrimaryVertex[fNParticles] = fOrigPrimaryVertex[i];
 
       fNParticles++;
     }
   }
 
   if (LOG_LEVEL(DEB)) {
     NUIS_LOG(DEB, "Ordered stack");
     for (int i = 0; i < fNParticles; i++) {
       NUIS_LOG(DEB, "Particle " << i << ". " << fParticlePDG[i] << " "
                             << fParticleMom[i][0] << " " << fParticleMom[i][1]
                             << " " << fParticleMom[i][2] << " "
                             << fParticleMom[i][3] << " " << fParticleState[i]);
     }
   }
 
   if (fNParticles != npart) {
     NUIS_ABORT("Dropped some particles when ordering the stack!");
   }
 
   return;
 }
 
 void FitEvent::Print() {
   if (LOG_LEVEL(FIT)) {
     NUIS_LOG(FIT, "FITEvent print");
     NUIS_LOG(FIT, "Mode: " << Mode << ", Weight: " << InputWeight);
     NUIS_LOG(FIT, "Particles: " << fNParticles);
     NUIS_LOG(FIT, " -> Particle Stack ");
     for (int i = 0; i < fNParticles; i++) {
       NUIS_LOG(FIT, " -> -> " << i << ". " << fParticlePDG[i] << " "
                           << fParticleState[i] << " "
                           << "  Mom(" << fParticleMom[i][0] << ", "
                           << fParticleMom[i][1] << ", " << fParticleMom[i][2]
                           << ", " << fParticleMom[i][3] << ").");
     }
   }
   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) {
     NUIS_ABORT("Requesting particle beyond stack!"
            << std::endl
            << "i = " << i << " N = " << fNParticles << std::endl
            << "Mode = " << Mode);
   }
 
   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::GetBeamPartIndex(void) const {
   return GetHMISParticleIndex(this->probe_pdg);
 }
 
 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;
 }
 
 // Get the outgoing lepton PDG depending on if it's a CC or NC event
 int FitEvent::GetLeptonOutPDG() {
   // Make sure the outgoing lepton has the correct PDG
   int pdgnu = PDGnu();
   int PDGout;
   if (IsCC()) {
     if (pdgnu > 0) PDGout = pdgnu-1;
     else PDGout = pdgnu+1;
   } else {
     PDGout = pdgnu;
   }
   return PDGout;
 }
 
 //********************************************************************
 // Returns the true Q2 of an event
 double FitEvent::GetQ2() {
   FitParticle *neutrino = GetNeutrinoIn();
   FitParticle *lepton = GetLeptonOut();
-  // Figure out why this happens
-  if (!neutrino || !lepton) return 0.0;
+  // Sometimes NEUT won't have an outgoing lepton because the event is Pauli blocked
+  if (!neutrino || !lepton) {
+    //if (!neutrino) std::cout << "no incoming neutrino!" << std::endl;
+    //if (!lepton) std::cout << "no outgoing lepton!" << std::endl;
+#ifdef __NEUT_ENABLED__
+    //fNeutVect->Dump();
+#endif
+    return -999;
+  }
   double Q2 =
     -1.0 * (lepton->P4() - neutrino->P4()) *
     (lepton->P4() - neutrino->P4()) / 1.E6;
   return Q2;
 }
 //********************************************************************
diff --git a/src/Reweight/MINERvAWeightCalcs.h b/src/Reweight/MINERvAWeightCalcs.h
index 9c6efe0..1a8a1ae 100644
--- a/src/Reweight/MINERvAWeightCalcs.h
+++ b/src/Reweight/MINERvAWeightCalcs.h
@@ -1,209 +1,209 @@
 #ifndef MINERVA_WEIGHT_CALCS
 #define MINERVA_WEIGHT_CALCS
 
 #include <string>
 
 #ifdef __MINERVA_RW_ENABLED__
 #ifdef __GENIE_ENABLED__
 #ifdef GENIE_PRE_R3
 #include "Conventions/Units.h"
 #include "EVGCore/EventRecord.h"
 #include "GHEP/GHepParticle.h"
 #include "GHEP/GHepRecord.h"
 #include "GHEP/GHepUtils.h"
 #include "Ntuple/NtpMCEventRecord.h"
 #include "PDG/PDGUtils.h"
 #else
 #include "Framework/Conventions/Units.h"
 #include "Framework/EventGen/EventRecord.h"
 #include "Framework/GHEP/GHepParticle.h"
 #include "Framework/GHEP/GHepRecord.h"
 #include "Framework/GHEP/GHepUtils.h"
 #include "Framework/Ntuple/NtpMCEventRecord.h"
 #include "Framework/ParticleData/PDGUtils.h"
+using namespace genie;
 #endif
 
 #include "NUISANCEWeightCalcs.h"
 #include "GeneralUtils.h"
 #include "NUISANCESyst.h"
 #include "FitEvent.h"
 #include "WeightUtils.h"
 #include "weightRPA.h"
-using namespace genie;
 class BaseFitEvt;
 
 
 namespace nuisance {
   namespace reweight {
 
     // QE Dials
     class MINERvAReWeight_QE : public NUISANCEWeightCalc {
       public:
         MINERvAReWeight_QE();
         virtual ~MINERvAReWeight_QE();
 
         double CalcWeight(BaseFitEvt* evt);
         void SetDialValue(std::string name, double val);
         void SetDialValue(int rwenum, double val);
         bool IsHandled(int rwenum);
 
         double fTwk_NormCCQE;
         double fCur_NormCCQE;
         double fDef_NormCCQE;
         bool fTweaked;
 
     };
 
 
     // MEC Dials
     class MINERvAReWeight_MEC : public NUISANCEWeightCalc {
       public:
         MINERvAReWeight_MEC();
         virtual ~MINERvAReWeight_MEC();
 
         double CalcWeight(BaseFitEvt* evt);
         void SetDialValue(std::string name, double val);
         void SetDialValue(int rwenum, double val);
         bool IsHandled(int rwenum);
 
         double fTwk_NormCCMEC;
         double fCur_NormCCMEC;
         double fDef_NormCCMEC;
         bool fTweaked;
 
     };
 
     // RES Dials
     class MINERvAReWeight_RES : public NUISANCEWeightCalc {
       public:
         MINERvAReWeight_RES();
         virtual ~MINERvAReWeight_RES();
 
         double CalcWeight(BaseFitEvt* evt);
         void SetDialValue(std::string name, double val);
         void SetDialValue(int rwenum, double val);
         bool IsHandled(int rwenum);
 
         double fTwk_NormCCRES;
         double fCur_NormCCRES;
         double fDef_NormCCRES;
         bool fTweaked;
 
     };
 
     /// RPA Weight Calculator that applies RPA systematics
     /// to GENIE events. GENIE EVENTS ONLY!
     class RikRPA : public NUISANCEWeightCalc {
       public:
         RikRPA();
         ~RikRPA();
 
         double CalcWeight(BaseFitEvt* evt);
         void SetDialValue(std::string name, double val);
         void SetDialValue(int rwenum, double val);
         bool IsHandled(int rwenum);
 
         void SetupRPACalculator(int calcenum);
         int GetRPACalcEnum(int bpdg, int tpdg);
 
         bool   fApplyDial_RPACorrection;
 
         double fTwkDial_RPALowQ2;
         double fDefDial_RPALowQ2;
         double fCurDial_RPALowQ2;
         double fErrDial_RPALowQ2;
 
         double fTwkDial_RPAHighQ2;
         double fDefDial_RPAHighQ2;
         double fCurDial_RPAHighQ2;
         double fErrDial_RPAHighQ2;
 
         bool   fApplyDial_RESRPACorrection;
 
         double fTwkDial_RESRPALowQ2;
         double fDefDial_RESRPALowQ2;
         double fCurDial_RESRPALowQ2;
         double fErrDial_RESRPALowQ2;
 
         double fTwkDial_RESRPAHighQ2;
         double fDefDial_RESRPAHighQ2;
         double fCurDial_RESRPAHighQ2;
         double fErrDial_RESRPAHighQ2;
 
         double* fEventWeights;
         bool fTweaked;
 
         const static int kMaxCalculators = 10;
         enum rpacalcenums {
           kNuMuC12,
           kNuMuO16,
           kNuMuAr40,
           kNuMuCa40,
           kNuMuFe56,
           kNuMuBarC12,
           kNuMuBarO16,
           kNuMuBarAr40,
           kNuMuBarCa40,
           kNuMuBarFe56
         };
         weightRPA* fRPACalculators[kMaxCalculators];
     };
 
     // Custom coherent tune from MINERvA
     class COHBrandon : public NUISANCEWeightCalc {
       public:
         COHBrandon();
         ~COHBrandon();
 
         double CalcWeight(BaseFitEvt* evt);
         void SetDialValue(std::string name, double val);
         void SetDialValue(int rwenum, double val);
         bool IsHandled(int rwenum);
 
         bool fApply_COHNorm;
 
         double fDef_COHNorm;
         double fCur_COHNorm;
         double fTwk_COHNorm;
 
         double fDef_COHCut;
         double fCur_COHCut;
         double fTwk_COHCut;
 
         bool fTweaked;
     };
 
     class WEnhancement : public NUISANCEWeightCalc {
       public:
         WEnhancement();
         ~WEnhancement();
 
         double CalcWeight(BaseFitEvt* evt);
         void SetDialValue(std::string name, double val);
         void SetDialValue(int rwenum, double val);
         bool IsHandled(int rwenum);
 
         bool fTweaked;
 
         bool fApply_Enhancement;
 
         double fDef_WNorm;
         double fCur_WNorm;
         double fTwk_WNorm;
 
         double fDef_WMean;
         double fCur_WMean;
         double fTwk_WMean;
 
         double fDef_WSigma;
         double fCur_WSigma;
         double fTwk_WSigma;
     };
 
   };  // namespace reweight
 };  // namespace nuisance
 
 #endif // __GENIE_ENABLED__
 #endif //__MINERVA_RW_ENABLED__
 
 #endif
diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx
index e523b9b..426f582 100644
--- a/src/Reweight/NUISANCEWeightCalcs.cxx
+++ b/src/Reweight/NUISANCEWeightCalcs.cxx
@@ -1,794 +1,792 @@
 #include "NUISANCEWeightCalcs.h"
 
 #include "FitEvent.h"
 #include "GeneralUtils.h"
 #include "NUISANCESyst.h"
 #include "WeightUtils.h"
 
 using namespace Reweight;
 
 ModeNormCalc::ModeNormCalc() { fNormRES = 1.0; }
 
 double ModeNormCalc::CalcWeight(BaseFitEvt *evt) {
   int mode = abs(evt->Mode);
   double w = 1.0;
 
   if (mode == 11 or mode == 12 or mode == 13) {
     w *= fNormRES;
   }
 
   return w;
 }
 
 void ModeNormCalc::SetDialValue(std::string name, double val) {
   SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
 }
 
 void ModeNormCalc::SetDialValue(int rwenum, double val) {
   int curenum = rwenum % 1000;
 
   // Check Handled
   if (!IsHandled(curenum)) 
     return;
   if (curenum == kModeNorm_NormRES) 
     fNormRES = val;
 }
 
 bool ModeNormCalc::IsHandled(int rwenum) {
   int curenum = rwenum % 1000;
   switch (curenum) {
   case kModeNorm_NormRES:
     return true;
   default:
     return false;
   }
 }
 
 //*****************************************************************************
 MINOSRPA::MINOSRPA() {
 
   fApply_MINOSRPA = false;
 
   fDef_MINOSRPA_A = 1.010;
   fTwk_MINOSRPA_A = 0.000;
   fDef_MINOSRPA_B = 0.156;
   fTwk_MINOSRPA_B = 0.000;
 
 }
 //
 double MINOSRPA::CalcWeight(BaseFitEvt* evt) {
   if (!fTweaked || !fApply_MINOSRPA) return 1.0;
 
   double w = 1.0;
 
   // If GENIE is enabled, use old code
 #ifdef __GENIE_ENABLED__
   // Extract the GENIE Record
   GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
   const Interaction* interaction = ghep->Summary();
   const InitialState& init_state = interaction->InitState();
   const ProcessInfo& proc_info = interaction->ProcInfo();
   const Target& tgt = init_state.Tgt();
 
-  // If not QE return 1.0
-  if (!tgt.IsNucleus()) return 1.0;
-  if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0;
+  // If not on nucleus, not resonant, or NC
+  if (!tgt.IsNucleus() || !proc_info.IsResonant() || proc_info.IsWeakNC()) return 1.0;
 
   // Extract Beam and Target PDG
   GHepParticle* neutrino = ghep->Probe();
   //   int bpdg = neutrino->Pdg();
 
   GHepParticle* target = ghep->Particle(1);
   assert(target);
   //    int tpdg = target->Pdg();
 
   // Extract Q0-Q3
   GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
   const TLorentzVector& k1 = *(neutrino->P4());
   const TLorentzVector& k2 = *(fsl->P4());
   //double q0 = fabs((k1 - k2).E());
   //double q3 = fabs((k1 - k2).Vect().Mag());
   double Q2 = fabs((k1 - k2).Mag2());
 
   w *= GetRPAWeight(Q2);
 #else
   // Get the Q2 from NUISANCE if not GENIE
   FitEvent *fevt = static_cast<FitEvent*>(evt);
   // Check the event is resonant
-  if (!fevt->IsResonant()) return 1.0;
+  if (!fevt->IsResonant() || fevt->IsNC()) return 1.0;
   int targeta = fevt->GetTargetA();
   int targetz = fevt->GetTargetZ();
   // Apply only to nuclear targets, ignore free protons
   if (targeta == 1 || targetz == 1) return 1.0;
   // Q2 in GeV2
   double Q2 = fevt->GetQ2();
   w *= GetRPAWeight(Q2);
 #endif
 
   return w;
 }
 
 // Do the actual weight calculation
 double MINOSRPA::GetRPAWeight(double Q2) {
   if (Q2 > 0.7) return 1.0;
   double w = fCur_MINOSRPA_A / (1.0 + TMath::Exp(1.0 - TMath::Sqrt(Q2) / fCur_MINOSRPA_B));
   return w;
 }
 
 bool MINOSRPA::IsHandled(int rwenum) {
   int curenum = rwenum % 1000;
   switch (curenum) {
     case Reweight::kMINERvARW_MINOSRPA_Apply:
     case Reweight::kMINERvARW_MINOSRPA_A:
     case Reweight::kMINERvARW_MINOSRPA_B:
       return true;
     default:
       return false;
   }
 }
 //
 void MINOSRPA::SetDialValue(std::string name, double val) {
   SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
 }
 //
 void MINOSRPA::SetDialValue(int rwenum, double val) {
   int curenum = rwenum % 1000;
 
   // Check Handled
   if (!IsHandled(curenum)) return;
   if (curenum == Reweight::kMINERvARW_MINOSRPA_Apply) fApply_MINOSRPA = (val > 0.5);
   if (curenum == Reweight::kMINERvARW_MINOSRPA_A) fTwk_MINOSRPA_A = val;
   if (curenum == Reweight::kMINERvARW_MINOSRPA_B) fTwk_MINOSRPA_B = val;
 
   // Check for changes
   fTweaked = (fApply_MINOSRPA ||
       fabs(fTwk_MINOSRPA_A) > 0.0 ||
       fabs(fTwk_MINOSRPA_B) > 0.0);
 
   // Update Values
   fCur_MINOSRPA_A = fDef_MINOSRPA_A * (1.0 + 0.1 * fTwk_MINOSRPA_A);
   fCur_MINOSRPA_B = fDef_MINOSRPA_B * (1.0 + 0.1 * fTwk_MINOSRPA_B);
 }
 
 //
 //
 //*****************************************************************************
 LagrangeRPA::LagrangeRPA() {
 
   fApplyRPA = false;
 
   /*
      fI1_Def = 4.18962e-01;
      fI1     = fI1_Def;
 
      fI2_Def = 7.39927e-01;
      fI2     = fI2_Def;
      */
 
   // Table VIII https://arxiv.org/pdf/1903.01558.pdf
   fR1_Def = 0.37;
   fR1     = fR1_Def;
 
   fR2_Def = 0.60;
   fR2     = fR2_Def;
 
 }
 //
 double LagrangeRPA::CalcWeight(BaseFitEvt* evt) {
 
   if (!fTweaked || !fApplyRPA) return 1.0;
   double w = 1.0;
 
   // If GENIE is enabled, use old code
 #ifdef __GENIE_ENABLED__
   // Extract the GENIE Record
   GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
   const Interaction* interaction = ghep->Summary();
   const InitialState& init_state = interaction->InitState();
   const ProcessInfo& proc_info = interaction->ProcInfo();
   const Target& tgt = init_state.Tgt();
 
-  // If not QE return 1.0
-  if (!tgt.IsNucleus()) return 1.0;
-  if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0;
+  // If not on nucleus, not resonant, or NC
+  if (!tgt.IsNucleus() || !proc_info.IsResonant() || proc_info.IsWeakNC()) return 1.0;
 
   // Extract Beam and Target PDG
   GHepParticle* neutrino = ghep->Probe();
   //   int bpdg = neutrino->Pdg();
 
   GHepParticle* target = ghep->Particle(1);
   assert(target);
   //    int tpdg = target->Pdg();
 
   // Extract Q0-Q3
   GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
   const TLorentzVector& k1 = *(neutrino->P4());
   const TLorentzVector& k2 = *(fsl->P4());
   //double q0 = fabs((k1 - k2).E());
   //double q3 = fabs((k1 - k2).Vect().Mag());
   double Q2 = fabs((k1 - k2).Mag2());
   w *= GetRPAWeight(Q2);
 #else
   // Get the Q2 from NUISANCE if not GENIE
   FitEvent *fevt = static_cast<FitEvent*>(evt);
   // Check the event is resonant
-  if (!fevt->IsResonant()) return 1.0;
+  if (!fevt->IsResonant() || fevt->IsNC()) return 1.0;
   int targeta = fevt->GetTargetA();
   int targetz = fevt->GetTargetZ();
   // Apply only to nuclear targets, ignore free protons
   if (targeta == 1 || targetz == 1) return 1.0;
   // Q2 in GeV2
   double Q2 = fevt->GetQ2();
   w *= GetRPAWeight(Q2);
 #endif
 
   return w;
 }
 //
 //
 double LagrangeRPA::GetRPAWeight(double Q2) {
   //std::cout << "Getting RPA Weight : " << Q2 << std::endl;
   if (Q2 > 0.7) return 1.0;
 
   // Keep original Lagrange RPA for documentation
   /*
      double x1 = 0.00;
      double x2 = 0.30;
      double x3 = 0.70;
 
      double y1 = 0.00;
      double y2 = fI2;
      double y3 = 1.00;
 
      double xv = Q2;
 
   // Automatically 0 because y1 choice
   double W1 = y1 * (xv-x2)*(xv-x3)/((x1-x2)*(x1-x3));
   double W2 = y2 * (xv-x1)*(xv-x3)/((x2-x1)*(x2-x3));
   double W3 = y3 * (xv-x1)*(xv-x2)/((x3-x1)*(x3-x2));
 
   double P = W1 + W2 + W3;
   double A1 = (1.0 - sqrt(1.0 - fI1));
   double R = P * (1.0 - A1) + A1;
 
   return 1.0 - (1.0-R)*(1.0-R);
   */
 
   // Equation 7 https://arxiv.org/pdf/1903.01558.pdf
   const double x1 = 0.00;
   const double x2 = 0.35;
   const double x3 = 0.70;
 
   // Equation 6 https://arxiv.org/pdf/1903.01558.pdf
   double RQ2 = fR2 *( (Q2-x1)*(Q2-x3)/((x2-x1)*(x2-x3)) )
     + (Q2-x1)*(Q2-x2)/((x3-x1)*(x3-x2));
   double weight = 1-(1-fR1)*(1-RQ2)*(1-RQ2);
 
   // Check range of R1 and R2
   // Commented out because this is implementation dependent: user may want strange double peaks
   /*
      if (fR1 > 1) return 1;
      if (fR2 > 1 || fR2 < 0.5) return 1;
      */
 
   return weight;
 }
 //
 bool LagrangeRPA::IsHandled(int rwenum) {
   int curenum = rwenum % 1000;
   switch (curenum) {
     case Reweight::kMINERvARW_LagrangeRPA_Apply:
     case Reweight::kMINERvARW_LagrangeRPA_R1:
     case Reweight::kMINERvARW_LagrangeRPA_R2:
       return true;
     default:
       return false;
   }
 }
 //
 void LagrangeRPA::SetDialValue(std::string name, double val) {
   SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
 }
 //
 void LagrangeRPA::SetDialValue(int rwenum, double val) {
   int curenum = rwenum % 1000;
 
   // Check Handled
   if (!IsHandled(curenum)) return;
   if (curenum == Reweight::kMINERvARW_LagrangeRPA_Apply) fApplyRPA = (val > 0.5);
   if (curenum == Reweight::kMINERvARW_LagrangeRPA_R1) fR1 = val;
   if (curenum == Reweight::kMINERvARW_LagrangeRPA_R2) fR2 = val;
 
   // Check for changes
   fTweaked = (fApplyRPA);
 }
 //
 
 
 BeRPACalc::BeRPACalc()
   : fBeRPA_A(0.59), fBeRPA_B(1.05), fBeRPA_D(1.13), fBeRPA_E(0.88),
   fBeRPA_U(1.2), nParams(0) {
     // A = 0.59 +/- 20%
     // B = 1.05 +/- 20%
     // D = 1.13 +/- 15%
     // E = 0.88 +/- 40%
     // U = 1.2
 }
 
 double BeRPACalc::CalcWeight(BaseFitEvt *evt) {
   int mode = abs(evt->Mode);
   double w = 1.0;
   if (nParams == 0) {
     return w;
   }
 
   // Get Q2
   // Get final state lepton
   if (mode == 1) {
     FitEvent *fevt = static_cast<FitEvent*>(evt);
     double Q2 = fevt->GetQ2();
     // Only CCQE events
     w *= calcRPA(Q2, fBeRPA_A, fBeRPA_B, fBeRPA_D, fBeRPA_E, fBeRPA_U);
   }
 
   return w;
 }
 
 void BeRPACalc::SetDialValue(std::string name, double val) {
   SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
 }
 
 void BeRPACalc::SetDialValue(int rwenum, double val) {
   int curenum = rwenum % 1000;
 
   // Check Handled
   if (!IsHandled(curenum))
     return;
   // Need 4 or 5 reconfigures
   if (curenum == kBeRPA_A)
     fBeRPA_A = val;
   else if (curenum == kBeRPA_B)
     fBeRPA_B = val;
   else if (curenum == kBeRPA_D)
     fBeRPA_D = val;
   else if (curenum == kBeRPA_E)
     fBeRPA_E = val;
   else if (curenum == kBeRPA_U)
     fBeRPA_U = val;
   nParams++;
 }
 
 bool BeRPACalc::IsHandled(int rwenum) {
   int curenum = rwenum % 1000;
   switch (curenum) {
     case kBeRPA_A:
     case kBeRPA_B:
     case kBeRPA_D:
     case kBeRPA_E:
     case kBeRPA_U:
       return true;
     default:
       return false;
   }
 }
 
 SBLOscWeightCalc::SBLOscWeightCalc() {
   fDistance = 0.0;
   fMassSplitting = 0.0;
   fSin2Theta = 0.0;
 }
 
 double SBLOscWeightCalc::CalcWeight(BaseFitEvt *evt) {
   FitEvent *fevt = static_cast<FitEvent *>(evt);
 
   FitParticle *pnu = fevt->PartInfo(0);
   double E = pnu->fP.E() / 1.E3;
 
   // Extract energy
   return GetSBLOscWeight(E);
 }
 
 void SBLOscWeightCalc::SetDialValue(std::string name, double val) {
   SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
 }
 
 void SBLOscWeightCalc::SetDialValue(int rwenum, double val) {
   int curenum = rwenum % 1000;
   if (!IsHandled(curenum))
     return;
   if (curenum == kSBLOsc_Distance)
     fDistance = val;
   if (curenum == kSBLOsc_MassSplitting)
     fMassSplitting = val;
   if (curenum == kSBLOsc_Sin2Theta)
     fSin2Theta = val;
 }
 
 bool SBLOscWeightCalc::IsHandled(int rwenum) {
   int curenum = rwenum % 1000;
   switch (curenum) {
     case kSBLOsc_Distance:
       return true;
     case kSBLOsc_MassSplitting:
       return true;
     case kSBLOsc_Sin2Theta:
       return true;
     default:
       return false;
   }
 }
 
 double SBLOscWeightCalc::GetSBLOscWeight(double E) {
   if (E <= 0.0)
     return 1.0 - 0.5 * fSin2Theta;
   return 1.0 - fSin2Theta * pow(sin(1.267 * fMassSplitting * fDistance / E), 2);
 }
 
 GaussianModeCorr::GaussianModeCorr() {
 
   // Apply the tilt-shift Gauss by Patrick
   // Alternatively set in config
   fMethod = true;
 
   // Init
   fApply_CCQE = false;
   fGausVal_CCQE[kPosNorm] = 0.0;
   fGausVal_CCQE[kPosTilt] = 0.0;
   fGausVal_CCQE[kPosPq0] = 1.0;
   fGausVal_CCQE[kPosWq0] = 1.0;
   fGausVal_CCQE[kPosPq3] = 1.0;
   fGausVal_CCQE[kPosWq3] = 1.0;
 
   fApply_2p2h = false;
   fGausVal_2p2h[kPosNorm] = 0.0;
   fGausVal_2p2h[kPosTilt] = 0.0;
   fGausVal_2p2h[kPosPq0] = 1.0;
   fGausVal_2p2h[kPosWq0] = 1.0;
   fGausVal_2p2h[kPosPq3] = 1.0;
   fGausVal_2p2h[kPosWq3] = 1.0;
 
   fApply_2p2h_PPandNN = false;
   fGausVal_2p2h_PPandNN[kPosNorm] = 0.0;
   fGausVal_2p2h_PPandNN[kPosTilt] = 0.0;
   fGausVal_2p2h_PPandNN[kPosPq0] = 1.0;
   fGausVal_2p2h_PPandNN[kPosWq0] = 1.0;
   fGausVal_2p2h_PPandNN[kPosPq3] = 1.0;
   fGausVal_2p2h_PPandNN[kPosWq3] = 1.0;
 
   fApply_2p2h_NP = false;
   fGausVal_2p2h_NP[kPosNorm] = 0.0;
   fGausVal_2p2h_NP[kPosTilt] = 0.0;
   fGausVal_2p2h_NP[kPosPq0] = 1.0;
   fGausVal_2p2h_NP[kPosWq0] = 1.0;
   fGausVal_2p2h_NP[kPosPq3] = 1.0;
   fGausVal_2p2h_NP[kPosWq3] = 1.0;
 
   fApply_CC1pi = false;
   fGausVal_CC1pi[kPosNorm] = 0.0;
   fGausVal_CC1pi[kPosTilt] = 0.0;
   fGausVal_CC1pi[kPosPq0] = 1.0;
   fGausVal_CC1pi[kPosWq0] = 1.0;
   fGausVal_CC1pi[kPosPq3] = 1.0;
   fGausVal_CC1pi[kPosWq3] = 1.0;
 
   fAllowSuppression = false;
 
   fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG");
   // fDebugStatements = true;
 }
 
 double GaussianModeCorr::CalcWeight(BaseFitEvt *evt) {
   FitEvent *fevt = static_cast<FitEvent *>(evt);
   double rw_weight = 1.0;
 
   // Get Neutrino
   if (!fevt->Npart()) {
     NUIS_ABORT("NO particles found in stack!");
   }
   FitParticle *pnu = fevt->GetNeutrinoIn();
   if (!pnu) {
     NUIS_ABORT("NO Starting particle found in stack!");
   }
 
   FitParticle *plep = fevt->GetLeptonOut();
   if (!plep) return 1.0;
 
   TLorentzVector q = pnu->fP - plep->fP;
 
   // Extra q0,q3
   double q0 = fabs(q.E()) / 1.E3;
   double q3 = fabs(q.Vect().Mag()) / 1.E3;
 
   int initialstate = -1; // Undef
   if (abs(fevt->Mode) == 2) {
     int npr = 0;
     int nne = 0;
 
     for (UInt_t j = 0; j < fevt->Npart(); j++) {
       if ((fevt->PartInfo(j))->fIsAlive)
         continue;
 
       if (fevt->PartInfo(j)->fPID == 2212) 
         npr++;
       else if (fevt->PartInfo(j)->fPID == 2112) 
         nne++;
     }
 
     if (fevt->Mode == 2 && npr == 1 && nne == 1) {
       initialstate = 2;
 
     } else if (fevt->Mode == 2 && 
         ((npr == 0 && nne == 2) || (npr == 2 && nne == 0))) {
       initialstate = 1;
     }
   }
 
   // Apply weighting
   if (fApply_CCQE && abs(fevt->Mode) == 1) {
     if (fDebugStatements)
       std::cout << "Getting CCQE Weight" << std::endl;
     double g = GetGausWeight(q0, q3, fGausVal_CCQE);
     if (g < 1.0) 
       g = 1.0;
     rw_weight *= g;
   }
 
   if (fApply_2p2h && abs(fevt->Mode) == 2) {
     if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl;
     if (fDebugStatements) std::cout << "Got q0 q3 = " << q0 << " " << q3 << " mode = " << fevt->Mode << std::endl;
     rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h);
     if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl;
   }
 
   if (fApply_2p2h_PPandNN && abs(fevt->Mode) == 2 && initialstate == 1) {
     if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl;
     rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN);
   }
 
   if (fApply_2p2h_NP && abs(fevt->Mode) == 2 && initialstate == 2) {
     if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl;
     rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP);
   }
 
   if (fApply_CC1pi && abs(fevt->Mode) >= 11 && abs(fevt->Mode) <= 13) {
     if (fDebugStatements)
       std::cout << "Getting CC1pi Weight" << std::endl;
     rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi);
   }
 
   return rw_weight;
 }
 
 void GaussianModeCorr::SetMethod(bool method) {
   fMethod = method;
   if (fMethod == true) {
     NUIS_LOG(FIT,
         " Using tilt-shift Gaussian parameters for Gaussian enhancement...");
   } else {
     NUIS_LOG(FIT, " Using Normal Gaussian parameters for Gaussian enhancement...");
   }
 };
 
 double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) {
   // The weight
   double w = 1.0;
 
   // Use tilt-shift method by Patrick
   if (fMethod) {
     if (fDebugStatements) {
       std::cout << "Using Patrick gaussian" << std::endl;
     }
     // // CCQE Without Suppression
     // double Norm = 4.82788679036;
     // double Tilt = 2.3501416116;
     // double Pq0  = 0.363964889702;
     // double Wq0  = 0.133976806938;
     // double Pq3  = 0.431769740224;
     // double Wq3  = 0.207666663434;
 
     // // Also add for CCQE at the end
     // return (w > 1.0) ? w : 1.0;
 
     // // 2p2h with suppression
     // double Norm = 15.967;
     // double Tilt = -0.455655;
     // double Pq0  = 0.214598;
     // double Wq0  = 0.0291061;
     // double Pq3  = 0.480194;
     // double Wq3  = 0.134588;
 
     double Norm = vals[kPosNorm];
     double Tilt = vals[kPosTilt];
     double Pq0 = vals[kPosPq0];
     double Wq0 = vals[kPosWq0];
     double Pq3 = vals[kPosPq3];
     double Wq3 = vals[kPosWq3];
 
     double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0);
     a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3);
 
     double b = -sin(2 * Tilt) / (4 * Wq0 * Wq0);
     b += sin(2 * Tilt) / (4 * Wq3 * Wq3);
 
     double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0);
     c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3);
 
     w = Norm;
     w *= exp(-a * (q0 - Pq0) * (q0 - Pq0));
     w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3));
     w *= exp(-c * (q3 - Pq3) * (q3 - Pq3));
 
     if (fDebugStatements) {
       std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " "
         << sin(Tilt) << std::endl;
       std::cout << "abc = " << a << " " << b << " " << c << std::endl;
       std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " "
         << Pq3 << " " << Wq3 << " " << w << std::endl;
     }
 
     if (w != w || std::isnan(w) || w < 0.0) {
       w = 0.0;
     }
 
     if (w < 1.0 and !fAllowSuppression) {
       w = 1.0;
     }
 
     return w;
 
     // Use the MINERvA Gaussian method
   } else {
     /*
      * From MINERvA and Daniel Ruterbories:
      * Old notes here: *
      * http://cdcvs.fnal.gov/cgi-bin/public-cvs/cvsweb-public.cgi/AnalysisFramework/Ana/CCQENu/ana_common/data/?cvsroot=mnvsoft
      * These parameters are slightly altered
      *
      * FRESH:
      * 10.5798
      * 0.254032
      * 0.50834
      * 0.0571035
      * 0.129051
      * 0.875287
      */
     if (fDebugStatements) {
       std::cout << "Using MINERvA Gaussian" << std::endl;
     }
 
     double norm = vals[kPosNorm];
     double meanq0 = vals[kPosTilt];
     double meanq3 = vals[kPosPq0];
     double sigmaq0 = vals[kPosWq0];
     double sigmaq3 = vals[kPosPq3];
     double corr = vals[kPosWq3];
 
     double z = (q0 - meanq0) * (q0 - meanq0) / (sigmaq0 * sigmaq0) +
       (q3 - meanq3) * (q3 - meanq3) / (sigmaq3 * sigmaq3) -
       2 * corr * (q0 - meanq0) * (q3 - meanq3) / (sigmaq0 * sigmaq3);
     double ret = 1.0;
 
     if ( fabs(1 - corr*corr) < 1.E-5 ) {
       return 1.0;
     }
 
     if ( (-0.5 * z / (1 - corr*corr)) > 200 or (-0.5 * z / (1 - corr*corr)) < -200 ) {
       return 1.0;
     } else {
       ret = norm * exp( -0.5 * z / (1 - corr*corr) );
     }
 
     if (ret != ret or ret < 0.0 or isnan(ret)) {
       return 1.0;
     }
 
     if (fAllowSuppression) return ret;
     return ret + 1.0;
     // Need to add 1 to the results
   }
   return w;
 }
 
 void GaussianModeCorr::SetDialValue(std::string name, double val) {
   SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
 }
 
 void GaussianModeCorr::SetDialValue(int rwenum, double val) {
   int curenum = rwenum % 1000;
 
   // Check Handled
   if (!IsHandled(curenum)) return;
 
   // CCQE Setting
   for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) {
     if (i == curenum) {
       int index = i - kGaussianCorr_CCQE_norm;
       fGausVal_CCQE[index] = val;
       fApply_CCQE = true;
     }
   }
 
   // 2p2h Setting
   for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) {
     if (i == curenum) {
       int index = i - kGaussianCorr_2p2h_norm;
       fGausVal_2p2h[index] = val;
       fApply_2p2h = true;
     }
   }
 
   // 2p2h_PPandNN Setting
   for (int i = kGaussianCorr_2p2h_PPandNN_norm; i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) {
     if (i == curenum) {
       int index = i - kGaussianCorr_2p2h_PPandNN_norm;
       fGausVal_2p2h_PPandNN[index] = val;
       fApply_2p2h_PPandNN = true;
     }
   }
 
   // 2p2h_NP Setting
   for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; i++) {
     if (i == curenum) {
       int index = i - kGaussianCorr_2p2h_NP_norm;
       fGausVal_2p2h_NP[index] = val;
       fApply_2p2h_NP = true;
     }
   }
 
   // CC1pi Setting
   for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) {
     if (i == curenum) {
       int index = i - kGaussianCorr_CC1pi_norm;
       fGausVal_CC1pi[index] = val;
       fApply_CC1pi = true;
     }
   }
 
   if (curenum == kGaussianCorr_AllowSuppression) {
     fAllowSuppression = (val > 0.5);
   }
 }
 
 bool GaussianModeCorr::IsHandled(int rwenum) {
   int curenum = rwenum % 1000;
   switch (curenum) {
     case kGaussianCorr_CCQE_norm:
     case kGaussianCorr_CCQE_tilt:
     case kGaussianCorr_CCQE_Pq0:
     case kGaussianCorr_CCQE_Wq0:
     case kGaussianCorr_CCQE_Pq3:
     case kGaussianCorr_CCQE_Wq3:
 
     case kGaussianCorr_2p2h_norm:
     case kGaussianCorr_2p2h_tilt:
     case kGaussianCorr_2p2h_Pq0:
     case kGaussianCorr_2p2h_Wq0:
     case kGaussianCorr_2p2h_Pq3:
     case kGaussianCorr_2p2h_Wq3:
 
     case kGaussianCorr_2p2h_PPandNN_norm:
     case kGaussianCorr_2p2h_PPandNN_tilt:
     case kGaussianCorr_2p2h_PPandNN_Pq0:
     case kGaussianCorr_2p2h_PPandNN_Wq0:
     case kGaussianCorr_2p2h_PPandNN_Pq3:
     case kGaussianCorr_2p2h_PPandNN_Wq3:
 
     case kGaussianCorr_2p2h_NP_norm:
     case kGaussianCorr_2p2h_NP_tilt:
     case kGaussianCorr_2p2h_NP_Pq0:
     case kGaussianCorr_2p2h_NP_Wq0:
     case kGaussianCorr_2p2h_NP_Pq3:
     case kGaussianCorr_2p2h_NP_Wq3:
 
     case kGaussianCorr_CC1pi_norm:
     case kGaussianCorr_CC1pi_tilt:
     case kGaussianCorr_CC1pi_Pq0:
     case kGaussianCorr_CC1pi_Wq0:
     case kGaussianCorr_CC1pi_Pq3:
     case kGaussianCorr_CC1pi_Wq3:
     case kGaussianCorr_AllowSuppression:
       return true;
     default:
       return false;
   }
 }
diff --git a/src/Reweight/NUISANCEWeightCalcs.h b/src/Reweight/NUISANCEWeightCalcs.h
index 520acfa..cd21a88 100644
--- a/src/Reweight/NUISANCEWeightCalcs.h
+++ b/src/Reweight/NUISANCEWeightCalcs.h
@@ -1,179 +1,199 @@
 #ifndef NUISANCE_WEIGHT_CALCS
 #define NUISANCE_WEIGHT_CALCS
 
 #include "BaseFitEvt.h"
 #include "BeRPA.h"
+#ifdef __GENIE_ENABLED__
+#ifdef GENIE_PRE_R3
+#include "Conventions/Units.h"
+#include "EVGCore/EventRecord.h"
+#include "GHEP/GHepParticle.h"
+#include "GHEP/GHepRecord.h"
+#include "GHEP/GHepUtils.h"
+#include "Ntuple/NtpMCEventRecord.h"
+#include "PDG/PDGUtils.h"
+#else
+#include "Framework/Conventions/Units.h"
+#include "Framework/EventGen/EventRecord.h"
+#include "Framework/GHEP/GHepParticle.h"
+#include "Framework/GHEP/GHepRecord.h"
+#include "Framework/GHEP/GHepUtils.h"
+#include "Framework/Ntuple/NtpMCEventRecord.h"
+#include "Framework/ParticleData/PDGUtils.h"
+using namespace genie;
+#endif // End GENIE pre v3
+#endif // End GENIE enabled
 
 class NUISANCEWeightCalc {
   public:
     NUISANCEWeightCalc() {};
     virtual ~NUISANCEWeightCalc() {};
 
     virtual double CalcWeight(BaseFitEvt* evt){return 1.0;};
     virtual void SetDialValue(std::string name, double val){};
     virtual void SetDialValue(int rwenum, double val){};
     virtual bool IsHandled(int rwenum){return false;};
 
     virtual void Print(){};
 
     std::map<std::string, int> fDialNameIndex;
     std::map<int, int> fDialEnumIndex;
     std::vector<double> fDialValues;
 
     std::string fName;
 };
 
 class ModeNormCalc : public NUISANCEWeightCalc {
   public:
     ModeNormCalc();
     ~ModeNormCalc(){};
 
     double CalcWeight(BaseFitEvt* evt);
     void SetDialValue(std::string name, double val);
     void SetDialValue(int rwenum, double val);
     bool IsHandled(int rwenum);
 
     double fNormRES;
 };
 
 // MINOS pion tuning, https://arxiv.org/pdf/1903.01558.pdf
 class MINOSRPA : public NUISANCEWeightCalc {
   public:
     MINOSRPA();
     ~MINOSRPA(){};
 
     double CalcWeight(BaseFitEvt* evt);
     void SetDialValue(std::string name, double val);
     void SetDialValue(int rwenum, double val);
     bool IsHandled(int rwenum);
 
     double GetRPAWeight(double Q2);
 
     bool fTweaked;
 
     bool fApply_MINOSRPA;
 
     double fTwk_MINOSRPA_A;
     double fDef_MINOSRPA_A;
     double fCur_MINOSRPA_A;
 
     double fTwk_MINOSRPA_B;
     double fDef_MINOSRPA_B;
     double fCur_MINOSRPA_B;
 
 };
 
 // MINERvA pion tuning, https://arxiv.org/pdf/1903.01558.pdf
 class LagrangeRPA : public NUISANCEWeightCalc {
   public:
     LagrangeRPA();
     ~LagrangeRPA(){};
 
     double CalcWeight(BaseFitEvt* evt);
     void SetDialValue(std::string name, double val);
     void SetDialValue(int rwenum, double val);
     bool IsHandled(int rwenum);
 
     double GetRPAWeight(double Q2);
 
     bool fTweaked;
 
     bool fApplyRPA;
 
     double fR1;
     double fR2;
     double fR1_Def;
     double fR2_Def;
 };
 
 class BeRPACalc : public NUISANCEWeightCalc {
   public:
     BeRPACalc();
     ~BeRPACalc(){};
 
     double CalcWeight(BaseFitEvt* evt);
     void SetDialValue(std::string name, double val);
     void SetDialValue(int rwenum, double val);
     bool IsHandled(int rwenum);
 
   private:
     // Parameter values
     double fBeRPA_A;
     double fBeRPA_B;
     double fBeRPA_D;
     double fBeRPA_E;
     double fBeRPA_U;
     // Counts of enabled parameters
     int nParams;
 
 };
 
 class SBLOscWeightCalc : public NUISANCEWeightCalc {
   public:
     SBLOscWeightCalc();
     ~SBLOscWeightCalc(){};
 
     double CalcWeight(BaseFitEvt* evt);
     void SetDialValue(std::string name, double val);
     void SetDialValue(int rwenum, double val);
     bool IsHandled(int rwenum);
 
     double GetSBLOscWeight(double E);
 
     double fDistance;
     double fMassSplitting;
     double fSin2Theta;
 
 };
 
 class GaussianModeCorr : public NUISANCEWeightCalc {
   public:
 
     GaussianModeCorr();
     ~GaussianModeCorr(){};
 
     double CalcWeight(BaseFitEvt* evt);
     void SetDialValue(std::string name, double val);
     void SetDialValue(int rwenum, double val);
     bool IsHandled(int rwenum);
     double GetGausWeight(double q0, double q3, double vals[]);
     // Set the Gaussian method (tilt-shift or normal Gaussian parameters)
     void SetMethod(bool method);
 
     // 5 pars describe the Gaussain
     // 0 norm.
     // 1 q0 pos
     // 2 q0 width
     // 3 q3 pos
     // 4 q3 width
     static const int kPosNorm = 0;
     static const int kPosTilt = 1;
     static const int kPosPq0  = 2;
     static const int kPosWq0  = 3;
     static const int kPosPq3  = 4;
     static const int kPosWq3  = 5;
 
     bool fApply_CCQE;
     double fGausVal_CCQE[6];
 
     bool fApply_2p2h;
     double fGausVal_2p2h[6];
 
     bool fApply_2p2h_PPandNN;
     double fGausVal_2p2h_PPandNN[6];
 
     bool fApply_2p2h_NP;
     double fGausVal_2p2h_NP[6];
 
     bool fApply_CC1pi;
     double fGausVal_CC1pi[6];
 
     bool fAllowSuppression;
 
     bool fDebugStatements;
 
     bool fMethod;
 };
 
 
 #endif