diff --git a/parameters/config.xml b/parameters/config.xml
index 3eb548b..48492b7 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,240 +1,232 @@
 <nuisance>
 <!-- # ###################################################### -->
 <!-- # NUISANCE CONFIGURATION OPTIONS -->
 <!-- # This file is read in by default at runtime -->
 <!-- # If you want to override on a case by case bases use -q at runtime -->
 <!-- # ###################################################### -->
 
 <!-- # MAIN Configs -->
 <!-- # ###################################################### -->
 
 <!-- # Logger goes from -->
 <!-- # 1 Quiet -->
 <!-- # 2 Fitter -->
 <!-- # 3 Samples -->
 <!-- # 4 Reconfigure Loops -->
 <!-- # 5 Every Event print out (SHOUT) -->
 <!-- # -1 DEBUGGING -->
 <config VERBOSITY='4'/>
 
 <!-- # ERROR goes from -->
 <!-- # 0 NONE -->
 <!-- # 1 FATAL -->
 <!-- # 2 WARN -->
 <config ERROR='2'/>
 <config TRACE='0'/>
 
 <config cores='1' />
 <config spline_test_throws='50' />
 <config spline_cores='1' />
 <config spline_chunks='20' />
 <config spline_procchunk='-1' />
 
 <config Electron_NThetaBins='4' />
 <config Electron_NEnergyBins='4' />
 <config Electron_ThetaWidth='1.0' />
 <config Electron_EnergyWidth='0.10' />
 
 <!-- Do we want to remove FSI, undefined and nuclear particles from the GENIE particle stack? -->
 <config RemoveFSIParticles='0' />
 <config RemoveUndefParticles='0' />
 <config RemoveNuclearParticles='0'/>
 <config MINERvASaveExtraCCQE='0' />
 
 <!-- # Input Configs -->
 <!-- # ###################################################### -->
 
 <!-- # Default Requirements file for the externalDataFitter Package -->
 <!-- # MAX Events : -1 is no cut. Events will be scaled automatically to give good xsec predictions. -->
 <config MAXEVENTS='-1'/>
 <!-- Include empty stacks in the THStack -->
 <config includeemptystackhists='0'/>
 
 <!-- # Turn on/off event manager -->
 <!-- # EventManager enables us to only loop number of events once for multiple projections of the same measurements -->
 <!-- # e.g. MiniBooNE CC1pi+ Q2 and MiniBooNE CC1pi+ Tmu would ordinarily require 2 reconfigures, but with this enabled it requires only one -->
 <config EventManager='1'/>
 
 <!-- # Event Directories -->
 <!-- # Can setup default directories and use @EVENT_DIR/path to link to it -->
 <config EVENT_DIR='/data2/stowell/NIWG/'/>
 <config NEUT_EVENT_DIR='/data2/stowell/NIWG/neut/fit_samples_neut5.3.3/'/>
 <config GENIE_EVENT_DIR='/data2/stowell/NIWG/genie/fit_samples_R.2.10.0/'/>
 <config NUWRO_EVENT_DIR='/data2/stowell/NIWG/nuwro/fit_samples/'/>
 <config GIBUU_EVENT_DIR='/data/GIBUU/DIR/'/>
 <config SaveNuWroExtra='0' />
 
 <!-- # In PrepareGENIE the reconstructed splines can be saved into the file -->
 <config save_genie_splines='1'/>
 
 <!-- # In InputHandler the option to regenerate NuWro flux/xsec plots is available -->
 <!-- # Going to move this to its own app soon -->
 
 <!-- # DEVEL CONFIG OPTION, don't touch! -->
 <config CacheSize='0'/>
 
 <!-- # ReWeighting Configuration Options -->
 <!-- # ###################################################### -->
 
 <!-- # Convert Dials in output statements using dial conversion card -->
 <config convert_dials='0'/>
 
 <!-- # Vetos can be used to specify RW dials NOT to be loaded in -->
 <!-- # Useful if one specific one has an issue -->
 <config FitWeight_fNIWGRW_veto=''/>
 <config FitWeight_fNuwroRW_veto=''/>
 <config FitWeight_fNeutRW_veto=''/>
 <config FitWeight_fGenieRW_veto=''/>
 
 <!-- # Output Options -->
 <!-- # ###################################################### -->
 
 <!-- # Save Nominal prediction with all rw engines at default -->
 <config savenominal='0'/>
 
 <!-- # Save prefit with values at starting values -->
 <config saveprefit='0'/>
 
 <!-- # Here's the full list of drawing options -->
 <!-- DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/WGHT/WEIGHTS/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS'/>
 <!-- #config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/RESIDUAL/MATRIX/FLUX/MASK/MAP -->
 <!-- #config drawopts DATA/MC -->
 <config drawopts='DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS/PROJ/CANVSLICEMC'/>
 
 <config nuisflat_SavePreFSI='true' />
 
 <config InterpolateSigmaQ0Histogram='1' />
 <config InterpolateSigmaQ0HistogramRes='100' />
 <config InterpolateSigmaQ0HistogramThrow='1' />
 <config InterpolateSigmaQ0HistogramNTHROWS='100000' />
 
 <!-- # Save the shape scaling applied with option SHAPE into the main MC hist -->
 <config saveshapescaling='0'/>
 
 <config CorrectGENIEMECNorm='1'/>
 
 <!-- # Set style of 1D output histograms -->
 <config linecolour='1'/>
 <config linestyle='1'/>
 <config linewidth='1'/>
 
 <!-- # For GenericFlux -->
 <config isLiteMode='0'/>
 
 <!-- # Statistical Options -->
 <!-- # ###################################################### -->
 
 <!-- # Add MC Statistical error to likelihoods -->
 <config addmcerror='0'/>
 
 <!-- # NUISMIN Configurations -->
 <!-- # ###################################################### -->
 
 <config MAXCALLS='1000000'/>
 <config MAXITERATIONS='1000000'/>
 <config TOLERANCE='0.001'/>
 
 <!-- # Number of events required in low stats routines -->
 <config LOWSTATEVENTS='25000'/>
 
 
 <!-- # Error band generator configs -->
 <!-- # ###################################################### -->
 
 <!-- # For -f ErrorBands creates error bands for given measurements -->
 <!-- # How many throws do we want (The higher the better precision) -->
 <config error_throws='500'/>
 
 <!-- # Are we throwing uniform or according to Gaussian? -->
 <!-- # Only use uniform if wanting to study the limits of a dial. -->
 <config error_uniform='0'/>
 <config WriteSeperateStacks='1'/>
 
 <!-- # Other Individual Case Configs -->
 <!-- # ###################################################### -->
 
 <!-- # Covariance throw options for fake data studies with MiniBooNE data. -->
 <config thrown_covariance='FULL'/>
 <config throw_mc_stat='0.0'/>
 <config throw_diag_syst='0'/>
 <config throw_corr_syst='0'/>
 <config throw_mc_stat='0'/>
 
 <!-- # Apply a shift to the muon momentum before calculation of Q2 -->
 <config muon_momentum_shift='0.0'/>
 <config muon_momentum_throw='0'/>
 
 <!-- # MINERvA Specific Configs -->
 <config MINERvA_CCinc_XSec_2DEavq3_nu_hadron_cut='0'/>
 <config MINERvA_CCinc_XSec_2DEavq3_nu_useq3true='0'/>
 <config Modes_split_PN_NN='0'/>
 
 <!-- Use only signal events when reconfiguring -->
 <config SignalReconfigures='false'/>
 <config FullEventOnSignalReconfigure="true"/>
 
 <!-- # SciBooNE specific -->
 <config SciBarDensity='1.04'/>
 <config SciBarRecoDist='12.0'/>
 <config PenetratingMuonEnergy='1.4'/>
 <config FlatEfficiency='1.0'/>
 <config NumRangeSteps='50'/>
 <config UseProton='true'/>
 <config UseZackEff='false'/>
 
 <config MINERvADensity='1.04'/>
 <config MINERvARecoDist='10.0'/>
 
 <!-- Different way of reweighting in GENIE -->
 <config GENIEWeightEngine_CCRESMode="kModeMaMv"/>
 <!--
 <config GENIEWeightEngine_CCRESMode="kModeMa"/>
 <<<<<<< HEAD
 -->
 <config GENIEWeightEngine_CCQEMode="kModeMa"/>
-||||||| merged common ancestors
--->
-<config GENIEWeightEngine_CCQEMode="kModeZExp"/> Using z-expansion
-=======
---> 
-<!-- Normal MAQE Llewelyn-Smith -->
-<config GENIEWeightEngine_CCQEMode="kModeMa"/>
->>>>>>> feature_prefsi
 <!--
 <config GENIEWeightEngine_CCQEMode="kModeNormAndMaShape"/> Taking shape into account (don't really use this...)
 <config GENIEWeightEngine_CCQEMode="kModeZExp"/> Using z-expansion
 -->
 
 <!-- CCQE/2p2h/1pi Gaussian enhancement method -->
 <!-- Apply tilt-shift weights or normal Gaussian parameters-->
 <config Gaussian_Enhancement="Normal" />
 <!--
 <config Gaussian_Enhancement="Tilt-Shift" />
 -->
 
 <config NToyThrows='100000' />
 
 <!-- Use NOvA Weights or not -->
 <config NOvA_Weights="false" />
 
 <!-- Tune name, for GENIE v3+ -->
 <config GENIETune="G18_02a_00_000" />
 <config GENIEEventGeneratorList="Default" />
 
 <!--
 <config GENIEXSecModelCCRES="genie::ReinSehgalRESPXSec" />
 <config GENIEXSecModelCOH="genie::ReinSehgalCOHPiPXSec" />
 <config GENIEXSecModelCCQE="genie::LwlynSmithQELCCPXSec" />
 -->
 
 <!--
 <config GENIEXSecModelCCQE="genie::NievesQELCCPXSec" />
 <config GENIEXSecModelCCRES="genie::BergerSehgalRESPXSec2014" />
 <config GENIEXSecModelCOH="genie::BergerSehgalCOHPiPXSec2015" />
 -->
 
 <config UseShapeCovar="0" />
 <config CC0piNBINS="156" />
 
 
 </nuisance>
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 5d5839d..b3b5420 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,586 +1,583 @@
 // 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"
 
 #pragma push_macro("LOG")
 #undef LOG
 #pragma push_macro("ERROR")
 #undef ERROR
 #ifdef GENIE_PRE_R3
 #include "Messenger/Messenger.h"
 #else
 #include "Framework/Messenger/Messenger.h"
 #endif
 #pragma pop_macro("LOG")
 #pragma pop_macro("ERROR")
 
 
 #include "InputUtils.h"
 
 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;
 
   genie::Messenger::Instance()->SetPriorityLevel("GHepUtils", pFATAL);
 
   // 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");
 
   // Are we running with NOvA weights
   fNOvAWeights = FitPar::Config().GetParB("NOvA_Weights");
   MAQEw = 1.0;
   NonResw = 1.0;
   RPAQEw = 1.0;
   RPARESw = 1.0;
   MECw = 1.0;
   DISw = 1.0;
   NOVAw = 1.0;
 
   // 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]);
     }
 
     // Check for precomputed weights
     TTree *weighttree = (TTree*)inp_file->Get("nova_wgts");
     if (fNOvAWeights) {
       if (!weighttree) {
         THROW("Did not find nova_wgts tree in file " << inputs[inp_it] << " but you specified it" << std::endl);
       } else {
         LOG(FIT) << "Found nova_wgts tree in file " << inputs[inp_it] << std::endl;
       }
     }
 
     // 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());
     if (weighttree != NULL) fGENIETree->AddFriend(weighttree);
   }
 
   // Registor all our file inputs
   SetupJointInputs();
 
   // Assign to tree
   fEventType = kGENIE;
   fGenieNtpl = NULL;
   fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
 
   // Set up the custom weights
   if (fNOvAWeights) {
     fGENIETree->SetBranchAddress("MAQEwgt", &MAQEw);
     fGENIETree->SetBranchAddress("nonResNormWgt", &NonResw);
     fGENIETree->SetBranchAddress("RPAQEWgt", &RPAQEw);
     fGENIETree->SetBranchAddress("RPARESWgt", &RPARESw);
     fGENIETree->SetBranchAddress("MECWgt", &MECw);
     fGENIETree->SetBranchAddress("DISWgt", &DISw);
     fGENIETree->SetBranchAddress("nova2018CVWgt", &NOVAw);
   }
 
   // Libraries should be seen but not heard...
   StopTalking();
   fGENIETree->GetEntry(0);
   StartTalking();
 
   // Create Fit Event
   fNUISANCEEvent = new FitEvent();
   fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
 
   if (fSaveExtra) {
     fGenieInfo = new GENIEGeneratorInfo();
     fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
   }
 
   fNUISANCEEvent->HardReset();
 };
 
 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;
 
   // Clear the previous event (See Note 1 in ROOT TClonesArray documentation)
   if (fGenieNtpl) {
     fGenieNtpl->Clear();
   }
 
   // Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
   fGENIETree->GetEntry(entry);
 
   // 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 = static_cast<GHepRecord*>(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 = static_cast<GHepRecord*>(fGenieNtpl->event);
   if (!fGenieGHep) return;
 
   // Convert GENIE Reaction Code
   fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
 
   // Set Event Info
   fNUISANCEEvent->fEventNo = 0.0;
   fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
   // Have a bool storing if interaction happened on free or bound nucleon
   bool IsFree = false;
   // Set the TargetPDG
   if (fGenieGHep->TargetNucleus() != NULL) {
     fNUISANCEEvent->fTargetPDG = fGenieGHep->TargetNucleus()->Pdg();
     IsFree = false;
   // Sometimes GENIE scatters off free nucleons, electrons, photons 
   // In which TargetNucleus is NULL and we need to find the initial state particle
   } else {
     // Check the particle is an initial state particle 
     // Follows GHepRecord::TargetNucleusPosition but doesn't do check on pdg::IsIon
     GHepParticle *p = fGenieGHep->Particle(1);
     // Check that particle 1 actually exists
     if (!p) {
       ERR(FTL) << "Can't find particle 1 for GHepRecord" << std::endl;
       throw;
     }
     // If not an ion but is an initial state particle
     if (!pdg::IsIon(p->Pdg()) && 
         p->Status() == kIStInitialState) {
       IsFree = true;
       fNUISANCEEvent->fTargetPDG = p->Pdg();
     // Catch if something strange happens: 
     // Here particle 1 is not an initial state particle OR
     // particle 1 is an ion OR
     // both
     } else {
       if (pdg::IsIon(p->Pdg())) {
         ERR(FTL) << "Particle 1 in GHepRecord stack is an ion but isn't an initial state particle" << std::endl;
         throw;
       } else {
         ERR(FTL) << "Particle 1 in GHepRecord stack is not an ion but is an initial state particle" << std::endl;
         throw;
       }
     }
   }
   // Set the A and Z and H from the target PDG
   // Depends on if we scattered off a free or bound nucleon
   if (!IsFree) {
     fNUISANCEEvent->fTargetA = TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
     fNUISANCEEvent->fTargetZ = TargetUtils::GetTargetZFromPDG(fNUISANCEEvent->fTargetPDG);
     fNUISANCEEvent->fTargetH = 0;
   } else {
     // If free proton scattering
     if (fNUISANCEEvent->fTargetPDG == 2212) {
       fNUISANCEEvent->fTargetA = 1;
       fNUISANCEEvent->fTargetZ = 1;
       fNUISANCEEvent->fTargetH = 1;
       // If free neutron scattering
     } else if (fNUISANCEEvent->fTargetPDG == 2112) {
       fNUISANCEEvent->fTargetA = 0;
       fNUISANCEEvent->fTargetZ = 1;
       fNUISANCEEvent->fTargetH = 0;
       // If neither
     } else {
       fNUISANCEEvent->fTargetA = 0;
       fNUISANCEEvent->fTargetZ = 0;
       fNUISANCEEvent->fTargetH = 0;
     }
   }
   fNUISANCEEvent->fBound = !IsFree;
   fNUISANCEEvent->InputWeight = 1.0;  //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
 
   // And the custom weights
   if (fNOvAWeights) {
     fNUISANCEEvent->CustomWeight = NOVAw;
     fNUISANCEEvent->CustomWeightArray[0] = MAQEw;
     fNUISANCEEvent->CustomWeightArray[1] = NonResw;
     fNUISANCEEvent->CustomWeightArray[2] = RPAQEw;
     fNUISANCEEvent->CustomWeightArray[3] = RPARESw;
     fNUISANCEEvent->CustomWeightArray[4] = MECw;
     fNUISANCEEvent->CustomWeightArray[5] = NOVAw;
   } else {
     fNUISANCEEvent->CustomWeight = 1.0;
     fNUISANCEEvent->CustomWeightArray[0] = 1.0;
     fNUISANCEEvent->CustomWeightArray[1] = 1.0;
     fNUISANCEEvent->CustomWeightArray[2] = 1.0;
     fNUISANCEEvent->CustomWeightArray[3] = 1.0;
     fNUISANCEEvent->CustomWeightArray[4] = 1.0;
     fNUISANCEEvent->CustomWeightArray[5] = 1.0;
   }
 
   // Get N Particle Stack
   unsigned int npart = fGenieGHep->GetEntries();
   unsigned int kmax = fNUISANCEEvent->kMaxParticles;
   if (npart > kmax) {
     ERR(WRN) << "GENIE has too many particles, expanding stack." << std::endl;
     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();
 
     // Set if the particle was on the fundamental vertex
     fNUISANCEEvent->fPrimaryVertex[curpart] = (p->FirstMother() < 2);
-    //std::cout << curpart << " " << p->Pdg() << " " << p->FirstMother() << std::endl;
-    //std::cout << p->LastMother() << std::endl;
-    //std::cout << "**" << std::endl;
 
     // 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) << "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/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index a8cec51..59b267a 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,383 +1,382 @@
 // 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 "GenericFlux_Vectors.h"
 
 GenericFlux_Vectors::GenericFlux_Vectors(std::string name,
                                          std::string inputfile, FitWeight *rw,
                                          std::string type,
                                          std::string fakeDataFile) {
   // Measurement Details
   fName = name;
   eventVariables = NULL;
 
   // Define our energy range for flux calcs
   EnuMin = 0.;
   EnuMax = 1E10;  // Arbritrarily high energy limit
 
   if (Config::HasPar("EnuMin")) {
     EnuMin = Config::GetParD("EnuMin");
   }
 
   if (Config::HasPar("EnuMax")) {
     EnuMax = Config::GetParD("EnuMax");
   }
 
   SavePreFSI = Config::Get().GetParB("nuisflat_SavePreFSI");
   LOG(SAM) << "Running GenericFlux_Vectors saving pre-FSI particles? " << SavePreFSI << std::endl;
 
   // Set default fitter flags
   fIsDiag = true;
   fIsShape = false;
   fIsRawEvents = false;
 
   // This function will sort out the input files automatically and parse all the
   // inputs,flags,etc.
   // There may be complex cases where you have to do this by hand, but usually
   // this will do.
   Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
 
   eventVariables = NULL;
 
   // Setup fDataHist as a placeholder
   this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
   this->SetupDefaultHist();
   fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
   covar = StatUtils::GetInvert(fFullCovar);
 
   // 1. The generator is organised in SetupMeasurement so it gives the
   // cross-section in "per nucleon" units.
   //    So some extra scaling for a specific measurement may be required. For
   //    Example to get a "per neutron" measurement on carbon
   //    which we do here, we have to multiple by the number of nucleons 12 and
   //    divide by the number of neutrons 6.
   // N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is
   // often included here in other classes that directly integrate the event
   // histogram. This method is used here as it now respects EnuMin and EnuMax
   // correctly.
   this->fScaleFactor =
       (this->PredictedEventRate("width", 0, EnuMax) / double(fNEvents)) /
       this->TotalIntegratedFlux("width");
 
   LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor
            << " [= " << (GetEventHistogram()->Integral("width") * 1E-38) << "/("
            << (fNEvents + 0.) << "*" << TotalIntegratedFlux("width") << ")]"
            << std::endl;
 
   if (fScaleFactor <= 0.0) {
     ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl;
     throw;
   }
 
   // Setup our TTrees
   this->AddEventVariablesToTree();
   this->AddSignalFlagsToTree();
 }
 
 void GenericFlux_Vectors::AddEventVariablesToTree() {
   // Setup the TTree to save everything
   if (!eventVariables) {
     Config::Get().out->cd();
     eventVariables = new TTree((this->fName + "_VARS").c_str(),
                                (this->fName + "_VARS").c_str());
   }
 
   LOG(SAM) << "Adding Event Variables" << std::endl;
 
   eventVariables->Branch("Mode", &Mode, "Mode/I");
   eventVariables->Branch("cc", &cc, "cc/B");
   eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I");
   eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F");
   eventVariables->Branch("tgt", &tgt, "tgt/I");
   eventVariables->Branch("tgta", &tgta, "tgta/I");
   eventVariables->Branch("tgtz", &tgtz, "tgtz/I");
   eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I");
   eventVariables->Branch("ELep", &ELep, "ELep/F");
   eventVariables->Branch("CosLep", &CosLep, "CosLep/F");
 
   // Basic interaction kinematics
   eventVariables->Branch("Q2", &Q2, "Q2/F");
   eventVariables->Branch("q0", &q0, "q0/F");
   eventVariables->Branch("q3", &q3, "q3/F");
   eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F");
   eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F");
   eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F");
   eventVariables->Branch("W", &W, "W/F");
   eventVariables->Branch("W_genie", &W_genie, "W_genie/F");
   eventVariables->Branch("x", &x, "x/F");
   eventVariables->Branch("y", &y, "y/F");
   eventVariables->Branch("Eav", &Eav, "Eav/F");
   eventVariables->Branch("EavAlt", &EavAlt, "EavAlt/F");
 
   eventVariables->Branch("dalphat", &dalphat, "dalphat/F");
   eventVariables->Branch("dpt", &dpt, "dpt/F");
   eventVariables->Branch("dphit", &dphit, "dphit/F");
   eventVariables->Branch("pnreco_C", &pnreco_C, "pnreco_C/F");
 
   // Save outgoing particle vectors
   eventVariables->Branch("nfsp", &nfsp, "nfsp/I");
   eventVariables->Branch("px", px, "px[nfsp]/F");
   eventVariables->Branch("py", py, "py[nfsp]/F");
   eventVariables->Branch("pz", pz, "pz[nfsp]/F");
   eventVariables->Branch("E", E, "E[nfsp]/F");
   eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I");
   eventVariables->Branch("status", status, "status[nfsp]/I");
   eventVariables->Branch("isalive", isalive, "isalive[nfsp]/O");
   eventVariables->Branch("isprimary", isprimary, "isprimary[nfsp]/O");
 
   // Event Scaling Information
   eventVariables->Branch("Weight", &Weight, "Weight/F");
   eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
   eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
   // Should be a double because may be 1E-39 and less
   eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D");
 
   // The customs
   eventVariables->Branch("CustomWeight", &CustomWeight, "CustomWeight/F");
   eventVariables->Branch("CustomWeightArray", CustomWeightArray, "CustomWeightArray[6]/F");
 
   return;
 }
 
 void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {
 
   ResetVariables();
 
   // Fill Signal Variables
   FillSignalFlags(event);
   LOG(DEB) << "Filling signal" << std::endl;
 
   // Now fill the information
   Mode = event->Mode;
   cc = (abs(event->Mode) < 30);
 
   // Get the incoming neutrino and outgoing lepton
   FitParticle *nu = event->GetNeutrinoIn();
   FitParticle *lep = event->GetHMFSAnyLepton();
 
   PDGnu = nu->fPID;
   Enu_true = nu->fP.E() / 1E3;
   tgt = event->fTargetPDG;
   tgta = event->fTargetA;
   tgtz = event->fTargetZ;
 
   if (lep != NULL) {
     PDGLep = lep->fPID;
     ELep = lep->fP.E() / 1E3;
     CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));
 
     // Basic interaction kinematics
     Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6;
     q0 = (nu->fP - lep->fP).E() / 1E3;
     q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3;
 
     // These assume C12 binding from MINERvA... not ideal
     Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
     Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);
 
     Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event)/1.E3;
     EavAlt = FitUtils::Eavailable(event)/1.E3;
 
     // Get W_true with assumption of initial state nucleon at rest
     float m_n = (float)PhysConst::mass_proton;
     // Q2 assuming nucleon at rest
     W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
     // True Q2
     W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
     x = Q2 / (2 * m_n * q0);
     y = 1 - ELep / Enu_true;
 
     dalphat = FitUtils::Get_STV_dalphat(event, PDGnu, true);
     dpt = FitUtils::Get_STV_dpt(event, PDGnu, true);
     dphit = FitUtils::Get_STV_dphit(event, PDGnu, true);
     pnreco_C = FitUtils::Get_pn_reco_C(event, PDGnu, true);
   }
 
   // Loop over the particles and store all the final state particles in a vector
   for (UInt_t i = 0; i < event->Npart(); ++i) {
 
     bool part_alive = event->PartInfo(i)->fIsAlive &&
       event->PartInfo(i)->Status() == kFinalState;
     if (!SavePreFSI) {
       if (!part_alive) continue;
     }
-    //std::cout << event->PartInfo(i)->Status() << std::endl;
 
     partList.push_back(event->PartInfo(i));
   }
 
   // Save outgoing particle vectors
   nfsp = (int)partList.size();
 
   for (int i = 0; i < nfsp; ++i) {
     px[i] = partList[i]->fP.X() / 1E3;
     py[i] = partList[i]->fP.Y() / 1E3;
     pz[i] = partList[i]->fP.Z() / 1E3;
     E[i] = partList[i]->fP.E() / 1E3;
     pdg[i] = partList[i]->fPID;
     status[i] = partList[i]->Status();
     isalive[i] = partList[i]->fIsAlive;
     isprimary[i] = event->fPrimaryVertex[i];
   }
 
 #ifdef __GENIE_ENABLED__
   if (event->fType == kGENIE) {
     EventRecord *  gevent      = static_cast<EventRecord*>(event->genie_event->event);
     const Interaction * interaction = gevent->Summary();
     const Kinematics &   kine       = interaction->Kine();
     double W_genie  = kine.W();
   }
 #endif
 
   // Fill event weights
   Weight = event->RWWeight * event->InputWeight;
   RWWeight = event->RWWeight;
   InputWeight = event->InputWeight;
   // And the Customs
   CustomWeight = event->CustomWeight;
   for (int i = 0; i < 6; ++i) {
     CustomWeightArray[i] = event->CustomWeightArray[i];
   }
 
   // Fill the eventVariables Tree
   eventVariables->Fill();
   return;
 };
 
 //********************************************************************
 void GenericFlux_Vectors::ResetVariables() {
   //********************************************************************
 
   cc = false;
 
   // Reset all Function used to extract any variables of interest to the event
   Mode = PDGnu = tgt = tgta = tgtz = PDGLep = 0;
 
   Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = Eav = EavAlt = -999.9;
 
   W_genie = -999;
   // Other fun variables
   // MINERvA-like ones
   dalphat = dpt = dphit = pnreco_C = -999.99;
 
   nfsp = 0;
   for (int i = 0; i < kMAX; ++i){
     px[i] = py[i] = pz[i] = E[i] = -999;
     pdg[i] = 0;
     status[i] = -999;
     isalive[i] = false;
     isprimary[i] = false;
   }
 
   Weight = InputWeight = RWWeight = 0.0;
 
   CustomWeight = 0.0;
   for (int i = 0; i < 6; ++i) CustomWeightArray[i] = 0.0;
 
   partList.clear();
 
   flagCCINC = flagNCINC = flagCCQE = flagCC0pi = flagCCQELike = flagNCEL = flagNC0pi = flagCCcoh = flagNCcoh = flagCC1pip = flagNC1pip = flagCC1pim = flagNC1pim = flagCC1pi0 = flagNC1pi0 = false;
 }
 
 //********************************************************************
 void GenericFlux_Vectors::FillSignalFlags(FitEvent *event) {
   //********************************************************************
 
   // Some example flags are given from SignalDef.
   // See src/Utils/SignalDef.cxx for more.
   int nuPDG = event->PartInfo(0)->fPID;
 
   // Generic signal flags
   flagCCINC = SignalDef::isCCINC(event, nuPDG);
   flagNCINC = SignalDef::isNCINC(event, nuPDG);
   flagCCQE = SignalDef::isCCQE(event, nuPDG);
   flagCCQELike = SignalDef::isCCQELike(event, nuPDG);
   flagCC0pi = SignalDef::isCC0pi(event, nuPDG);
   flagNCEL = SignalDef::isNCEL(event, nuPDG);
   flagNC0pi = SignalDef::isNC0pi(event, nuPDG);
   flagCCcoh = SignalDef::isCCCOH(event, nuPDG, 211);
   flagNCcoh = SignalDef::isNCCOH(event, nuPDG, 111);
   flagCC1pip = SignalDef::isCC1pi(event, nuPDG, 211);
   flagNC1pip = SignalDef::isNC1pi(event, nuPDG, 211);
   flagCC1pim = SignalDef::isCC1pi(event, nuPDG, -211);
   flagNC1pim = SignalDef::isNC1pi(event, nuPDG, -211);
   flagCC1pi0 = SignalDef::isCC1pi(event, nuPDG, 111);
   flagNC1pi0 = SignalDef::isNC1pi(event, nuPDG, 111);
 }
 
 void GenericFlux_Vectors::AddSignalFlagsToTree() {
   if (!eventVariables) {
     Config::Get().out->cd();
     eventVariables = new TTree((this->fName + "_VARS").c_str(),
         (this->fName + "_VARS").c_str());
   }
 
   LOG(SAM) << "Adding signal flags" << std::endl;
 
   // Signal Definitions from SignalDef.cxx
   eventVariables->Branch("flagCCINC", &flagCCINC, "flagCCINC/O");
   eventVariables->Branch("flagNCINC", &flagNCINC, "flagNCINC/O");
   eventVariables->Branch("flagCCQE", &flagCCQE, "flagCCQE/O");
   eventVariables->Branch("flagCC0pi", &flagCC0pi, "flagCC0pi/O");
   eventVariables->Branch("flagCCQELike", &flagCCQELike, "flagCCQELike/O");
   eventVariables->Branch("flagNCEL", &flagNCEL, "flagNCEL/O");
   eventVariables->Branch("flagNC0pi", &flagNC0pi, "flagNC0pi/O");
   eventVariables->Branch("flagCCcoh", &flagCCcoh, "flagCCcoh/O");
   eventVariables->Branch("flagNCcoh", &flagNCcoh, "flagNCcoh/O");
   eventVariables->Branch("flagCC1pip", &flagCC1pip, "flagCC1pip/O");
   eventVariables->Branch("flagNC1pip", &flagNC1pip, "flagNC1pip/O");
   eventVariables->Branch("flagCC1pim", &flagCC1pim, "flagCC1pim/O");
   eventVariables->Branch("flagNC1pim", &flagNC1pim, "flagNC1pim/O");
   eventVariables->Branch("flagCC1pi0", &flagCC1pi0, "flagCC1pi0/O");
   eventVariables->Branch("flagNC1pi0", &flagNC1pi0, "flagNC1pi0/O");
 };
 
 void GenericFlux_Vectors::Write(std::string drawOpt) {
 
   // First save the TTree
   eventVariables->Write();
 
   // Save Flux and Event Histograms too
   GetInput()->GetFluxHistogram()->Write();
   GetInput()->GetEventHistogram()->Write();
 
   return;
 }
 
 // Override functions which aren't really necessary
 bool GenericFlux_Vectors::isSignal(FitEvent *event) {
   (void)event;
   return true;
 };
 
 void GenericFlux_Vectors::ScaleEvents() { return; }
 
 void GenericFlux_Vectors::ApplyNormScale(float norm) {
   this->fCurrentNorm = norm;
   return;
 }
 
 void GenericFlux_Vectors::FillHistograms() { return; }
 
 void GenericFlux_Vectors::ResetAll() {
   eventVariables->Reset();
   return;
 }
 
 float GenericFlux_Vectors::GetChi2() { return 0.0; }