diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index b3b5420..f779ec5 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,583 +1,581 @@
 // 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);
 
     // 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();
+  FitParticle* ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
+  if (ISAnyLepton) {
+    fNUISANCEEvent->probe_E = ISAnyLepton->E();
+    fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
   }
-
   return;
 }
 
 void GENIEInputHandler::Print() {}
 
 #endif
diff --git a/src/InputHandler/GIBUUInputHandler.cxx b/src/InputHandler/GIBUUInputHandler.cxx
index f61559f..1e31f41 100644
--- a/src/InputHandler/GIBUUInputHandler.cxx
+++ b/src/InputHandler/GIBUUInputHandler.cxx
@@ -1,302 +1,301 @@
 #ifdef __GiBUU_ENABLED__
 #include "GIBUUInputHandler.h"
 #include "InputUtils.h"
 
 
 GIBUUGeneratorInfo::~GIBUUGeneratorInfo() { DeallocateParticleStack(); }
 
 void GIBUUGeneratorInfo::AddBranchesToTree(TTree* tn) {
   // tn->Branch("NEUTParticleN",          fNEUTParticleN, "NEUTParticleN/I");
   // tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
   // "NEUTParticleStatusCode[NEUTParticleN]/I");
   // tn->Branch("NEUTParticleAliveCode",  fNEUTParticleAliveCode,
   // "NEUTParticleAliveCode[NEUTParticleN]/I");
 }
 
 void GIBUUGeneratorInfo::SetBranchesFromTree(TTree* tn) {
   // tn->SetBranchAddress("NEUTParticleN",          &fNEUTParticleN );
   // tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode );
   // tn->SetBranchAddress("NEUTParticleAliveCode",  &fNEUTParticleAliveCode  );
 }
 
 void GIBUUGeneratorInfo::AllocateParticleStack(int stacksize) {
   // fNEUTParticleN = 0;
   // fNEUTParticleStatusCode = new int[stacksize];
   // fNEUTParticleStatusCode = new int[stacksize];
 }
 
 void GIBUUGeneratorInfo::DeallocateParticleStack() {
   // delete fNEUTParticleStatusCode;
   // delete fNEUTParticleAliveCode;
 }
 
 void GIBUUGeneratorInfo::FillGeneratorInfo(GiBUUStdHepReader* nevent) {
   Reset();
   // for (int i = 0; i < nevent->Npart(); i++) {
   // fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
   // fNEUTParticleAliveCode[i]  = nevent->PartInfo(i)->fIsAlive;
   // fNEUTParticleN++;
   // }
 }
 
 void GIBUUGeneratorInfo::Reset() {
   // for (int i = 0; i < fNEUTParticleN; i++) {
   // fNEUTParticleStatusCode[i] = -1;
   // fNEUTParticleAliveCode[i]  = 9;
   // }
   // fNEUTParticleN = 0;
 }
 
 GIBUUInputHandler::GIBUUInputHandler(std::string const& handle,
                                      std::string const& rawinputs) {
   LOG(SAM) << "Creating GiBUUInputHandler : " << handle << std::endl;
 
   // Run a joint input handling
   fName = handle;
   fEventType = kGiBUU;
   fGIBUUTree = new TChain("giRooTracker");
 
   // 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
     LOG(SAM) << "Opening event file " << inputs[inp_it] << std::endl;
     TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ");
     if ((!inp_file) || (!inp_file->IsOpen())) {
       THROW("GiBUU file !IsOpen() at : '"
             << inputs[inp_it] << "'" << std::endl
             << "Check that your file paths are correct and the file exists!");
     }
 
     int NFluxes = bool(dynamic_cast<TH1D*>(inp_file->Get("numu_flux"))) +
                   bool(dynamic_cast<TH1D*>(inp_file->Get("numub_flux"))) +
                   bool(dynamic_cast<TH1D*>(inp_file->Get("nue_flux"))) +
                   bool(dynamic_cast<TH1D*>(inp_file->Get("nueb_flux"))) +
                   bool(dynamic_cast<TH1D*>(inp_file->Get("e_flux")));
 
     if (NFluxes != 1) {
       THROW("Found " << NFluxes << " input fluxes in " << inputs[inp_it]
                      << ". The NUISANCE GiBUU interface expects to be "
                         "passed multiple species vectors as separate "
                         "input files like: "
                         "\"GiBUU:(MINERVA_FHC_numu_evts.root,MINERVA_FHC_"
                         "numubar_evts.root,[...])\"");
     }
 
     // Get Flux/Event hist
     TH1D* fluxhist = dynamic_cast<TH1D*>(inp_file->Get("flux"));
     TH1D* eventhist = dynamic_cast<TH1D*>(inp_file->Get("evt"));
     if (!fluxhist || !eventhist) {
       ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
       inp_file->ls();
       THROW(
           "GiBUU FILE doesn't contain flux/xsec info. You may have to "
           "regenerate your MC!");
     }
 
     // Get N Events
     TTree* giRooTracker = dynamic_cast<TTree*>(inp_file->Get("giRooTracker"));
     if (!giRooTracker) {
       ERROR(FTL,
             "giRooTracker Tree not located in NEUT file: " << inputs[inp_it]);
       THROW("Check your inputs, they may need to be completely regenerated!");
       throw;
     }
     int nevents = giRooTracker->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
     fGIBUUTree->AddFile(inputs[inp_it].c_str());
   }
 
   // Registor all our file inputs
   SetupJointInputs();
 
   // Create Fit Event
   fNUISANCEEvent = new FitEvent();
 
   fGiReader = new GiBUUStdHepReader();
   fGiReader->SetBranchAddresses(fGIBUUTree);
 
   fNUISANCEEvent->HardReset();
 };
 
 FitEvent* GIBUUInputHandler::GetNuisanceEvent(const UInt_t entry,
                                               const bool lightweight) {
   // Check out of bounds
   if (entry >= (UInt_t)fNEvents) return NULL;
 
   // Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
   fGIBUUTree->GetEntry(entry);
 
   // Run NUISANCE Vector Filler
   if (!lightweight) {
     CalcNUISANCEKinematics();
   }
 #ifdef __PROB3PP_ENABLED__
   else {
     for (int i = 0; i < fGiReader->StdHepN; i++) {
       int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i],
                                          fGiReader->StdHepPdg[i]);
       if (state != kInitialState) {
         continue;
       }
       if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
                      fGiReader->StdHepPdg[i])) {
         fNUISANCEEvent->probe_E = fGiReader->StdHepP4[i][3] * 1.E3;
         fNUISANCEEvent->probe_pdg = fGiReader->StdHepPdg[i];
         break;
       }
     }
   }
 #endif
 
   fNUISANCEEvent->InputWeight *= GetInputWeight(entry);
   fNUISANCEEvent->GiRead = fGiReader;
 
   return fNUISANCEEvent;
 }
 
 int GetGIBUUParticleStatus(int status, int pdg) {
   int state = kUndefinedState;
   switch (status) {
     case 0:   // Incoming
     case 11:  // Struck nucleon
       state = kInitialState;
       break;
 
     case 1:  // Good Final State
       state = kFinalState;
       break;
 
     default:  // Other
       break;
   }
 
   // Set Nuclear States Flag
   if (pdg > 1000000) {
     if (state == kInitialState)
       state = kNuclearInitial;
     else if (state == kFinalState)
       state = kNuclearRemnant;
     else
       state = kUndefinedState;
   }
 
   return state;
 }
 
 void GIBUUInputHandler::CalcNUISANCEKinematics() {
   // Reset all variables
   fNUISANCEEvent->ResetEvent();
   FitEvent* evt = fNUISANCEEvent;
   evt->Mode = fGiReader->GiBUU2NeutCode;
   evt->fEventNo = 0.0;
   evt->fTotCrs = 0;
   evt->fTargetA = 0.0;  // Change to get these from nuclear remnant.
   evt->fTargetZ = 0.0;
   evt->fTargetH = 0;
   evt->fBound = 0.0;
 
   // Extra GiBUU Input Weight
   evt->InputWeight = fGiReader->EvtWght;
 
   // Check Stack N
   int npart = fGiReader->StdHepN;
   int kmax = evt->kMaxParticles;
   if ((UInt_t)npart > (UInt_t)kmax) {
     ERROR(WRN, "GiBUU has too many particles. Expanding Stack.");
     fNUISANCEEvent->ExpandParticleStack(npart);
   }
 
   // Create Stack
   evt->fNParticles = 0;
   for (int i = 0; i < npart; i++) {
     // State
     int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i],
                                        fGiReader->StdHepPdg[i]);
     int curpart = evt->fNParticles;
 
     // Set State
     evt->fParticleState[evt->fNParticles] = state;
 
     // Mom
     evt->fParticleMom[curpart][0] = fGiReader->StdHepP4[i][0] * 1.E3;
     evt->fParticleMom[curpart][1] = fGiReader->StdHepP4[i][1] * 1.E3;
     evt->fParticleMom[curpart][2] = fGiReader->StdHepP4[i][2] * 1.E3;
     evt->fParticleMom[curpart][3] = fGiReader->StdHepP4[i][3] * 1.E3;
 
     // PDG
     evt->fParticlePDG[curpart] = fGiReader->StdHepPdg[i];
 
     // Add to total particles
     evt->fNParticles++;
   }
 
   // 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();
+  FitParticle* ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
+  if (ISAnyLepton) {
+    fNUISANCEEvent->probe_E = ISAnyLepton->E();
+    fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
   }
 
   return;
 }
 
 void GIBUUInputHandler::Print() {}
 
 void GIBUUInputHandler::SetupJointInputs() {
   if (jointeventinputs.size() <= 1) {
     jointinput = false;
   } else if (jointeventinputs.size() > 1) {
     jointinput = true;
     jointindexswitch = 0;
   }
   fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
   if (fMaxEvents != -1 and jointeventinputs.size() > 1) {
     THROW("Can only handle joint inputs when config MAXEVENTS = -1!");
   }
 
   for (size_t i = 0; i < jointeventinputs.size(); i++) {
     double scale = double(fNEvents) / fEventHist->Integral("width");
     scale *= jointfluxinputs.at(i)->Integral("width");
 
     jointindexscale.push_back(scale);
   }
 
   fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str());
   fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str());
 
   // Setup Max Events
   if (fMaxEvents > 1 && fMaxEvents < fNEvents) {
     if (LOG_LEVEL(SAM)) {
       std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl;
     }
     fNEvents = fMaxEvents;
   }
 
   // Print out Status
   if (LOG_LEVEL(SAM)) {
     std::cout << "\t\t|-> Total Entries    : " << fNEvents << std::endl
               << "\t\t|-> Event Integral   : "
               << fEventHist->Integral("width") * 1.E-38 << " events/nucleon"
               << std::endl
               << "\t\t|-> Flux Integral    : " << fFluxHist->Integral("width")
               << " /cm2" << std::endl
               << "\t\t|-> Event/Flux       : "
               << fEventHist->Integral("width") * 1.E-38 /
                      fFluxHist->Integral("width")
               << " cm2/nucleon" << std::endl;
   }
 }
 
 #endif
diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
index b3fec8f..e4ce5aa 100644
--- a/src/InputHandler/NEUTInputHandler.cxx
+++ b/src/InputHandler/NEUTInputHandler.cxx
@@ -1,510 +1,509 @@
 #ifdef __NEUT_ENABLED__
 #include "NEUTInputHandler.h"
 #include "InputUtils.h"
 
 
 NEUTGeneratorInfo::~NEUTGeneratorInfo() { DeallocateParticleStack(); }
 
 void NEUTGeneratorInfo::AddBranchesToTree(TTree* tn) {
   tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
   tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
              "NEUTParticleStatusCode[NEUTParticleN]/I");
   tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode,
              "NEUTParticleAliveCode[NEUTParticleN]/I");
 }
 
 void NEUTGeneratorInfo::SetBranchesFromTree(TTree* tn) {
   tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN);
   tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode);
   tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode);
 }
 
 void NEUTGeneratorInfo::AllocateParticleStack(int stacksize) {
   fNEUTParticleN = 0;
   fNEUTParticleStatusCode = new int[stacksize];
   fNEUTParticleStatusCode = new int[stacksize];
 }
 
 void NEUTGeneratorInfo::DeallocateParticleStack() {
   delete fNEUTParticleStatusCode;
   delete fNEUTParticleAliveCode;
 }
 
 void NEUTGeneratorInfo::FillGeneratorInfo(NeutVect* nevent) {
   Reset();
   for (int i = 0; i < nevent->Npart(); i++) {
     fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
     fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
     fNEUTParticleN++;
   }
 }
 
 void NEUTGeneratorInfo::Reset() {
   for (int i = 0; i < fNEUTParticleN; i++) {
     fNEUTParticleStatusCode[i] = -1;
     fNEUTParticleAliveCode[i] = 9;
   }
   fNEUTParticleN = 0;
 }
 
 NEUTInputHandler::NEUTInputHandler(std::string const& handle,
                                    std::string const& rawinputs) {
   LOG(SAM) << "Creating NEUTInputHandler : " << handle << std::endl;
 
   // Run a joint input handling
   fName = handle;
 
   // Setup the TChain
   fNEUTTree = new TChain("neuttree");
   fSaveExtra = FitPar::Config().GetParB("SaveExtraNEUT");
   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(inputs[inp_it].c_str(), "READ");
     if (!inp_file or inp_file->IsZombie()) {
       THROW("NEUT 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(
         (PlotUtils::GetObjectWithName(inp_file, "flux")).c_str());
     TH1D* eventhist = (TH1D*)inp_file->Get(
         (PlotUtils::GetObjectWithName(inp_file, "evt")).c_str());
     if (!fluxhist or !eventhist) {
       ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
       inp_file->ls();
       THROW(
           "NEUT FILE doesn't contain flux/xsec info. You may have to "
           "regenerate your MC!");
     }
 
     // Get N Events
     TTree* neuttree = (TTree*)inp_file->Get("neuttree");
     if (!neuttree) {
       ERROR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]);
       THROW("Check your inputs, they may need to be completely regenerated!");
       throw;
     }
     int nevents = neuttree->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
     fNEUTTree->AddFile(inputs[inp_it].c_str());
   }
 
   // Registor all our file inputs
   SetupJointInputs();
 
   // Assign to tree
   fEventType = kNEUT;
   fNeutVect = NULL;
   fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect);
   fNEUTTree->GetEntry(0);
 
   // Create Fit Event
   fNUISANCEEvent = new FitEvent();
   fNUISANCEEvent->SetNeutVect(fNeutVect);
 
   if (fSaveExtra) {
     fNeutInfo = new NEUTGeneratorInfo();
     fNUISANCEEvent->AddGeneratorInfo(fNeutInfo);
   }
 
   fNUISANCEEvent->HardReset();
 };
 
 NEUTInputHandler::~NEUTInputHandler(){
     //  if (fNEUTTree) delete fNEUTTree;
     //  if (fNeutVect) delete fNeutVect;
     //  if (fNeutInfo) delete fNeutInfo;
 };
 
 void NEUTInputHandler::CreateCache() {
   if (fCacheSize > 0) {
     // fNEUTTree->SetCacheEntryRange(0, fNEvents);
     fNEUTTree->AddBranchToCache("vectorbranch", 1);
     fNEUTTree->SetCacheSize(fCacheSize);
   }
 }
 
 void NEUTInputHandler::RemoveCache() {
   // fNEUTTree->SetCacheEntryRange(0, fNEvents);
   fNEUTTree->AddBranchToCache("vectorbranch", 0);
   fNEUTTree->SetCacheSize(0);
 }
 
 FitEvent* NEUTInputHandler::GetNuisanceEvent(const UInt_t entry,
                                              const bool lightweight) {
   // Catch too large entries
   if (entry >= (UInt_t)fNEvents) return NULL;
 
   // Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
   fNEUTTree->GetEntry(entry);
 
   // Run NUISANCE Vector Filler
   if (!lightweight) {
     CalcNUISANCEKinematics();
   }
 #ifdef __PROB3PP_ENABLED__
   else {
 
     UInt_t npart = fNeutVect->Npart();
     for (size_t i = 0; i < npart; i++) {
       NeutPart* part = fNUISANCEEvent->fNeutVect->PartInfo(i);
       if ((part->fIsAlive == false) && (part->fStatus == -1) &&
           std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
                      part->fPID)) {
         fNUISANCEEvent->probe_E = part->fP.T();
         fNUISANCEEvent->probe_pdg = part->fPID;
         break;
       } else {
         continue;
       }
     }
   }
 #endif
 
   // Setup Input scaling for joint inputs
   fNUISANCEEvent->InputWeight = GetInputWeight(entry);
 
   // Return event pointer
   return fNUISANCEEvent;
 }
 
 // From NEUT neutclass/neutpart.h
 //          Bool_t         fIsAlive; // Particle should be tracked or not
 //                          ( in the detector simulator )
 //
 //         Int_t          fStatus;  // Status flag of this particle
 //                            -2: Non existing particle
 //                            -1: Initial state particle
 //                             0: Normal
 //                             1: Decayed to the other particle
 //                             2: Escaped from the detector
 //                             3: Absorped
 //                             4: Charge exchanged
 //                             5: Pauli blocked
 //                             6: N/A
 //                             7: Produced child particles
 //                             8: Inelastically scattered
 //
 int NEUTInputHandler::GetNeutParticleStatus(NeutPart* part) {
   // State
   int state = kUndefinedState;
 
   // Remove Pauli blocked events, probably just single pion events
   if (part->fStatus == 5) {
     state = kFSIState;
 
   // fStatus == -1 means initial  state
   } else if (part->fIsAlive == false && part->fStatus == -1) {
     state = kInitialState;
 
     // NEUT has a bit of a strange convention for fIsAlive and fStatus
     // combinations
     // for NC and neutrino particle isAlive true/false and status 2 means
     // final state particle
     // for other particles in NC status 2 means it's an FSI particle
     // for CC it means it was an FSI particle
   } else if (part->fStatus == 2) {
     // NC case is a little strange... The outgoing neutrino might be alive or
     // not alive. Remaining particles with status 2 are FSI particles that
     // reinteracted
     if (abs(fNeutVect->Mode) > 30 &&
         (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) {
       state = kFinalState;
       // The usual CC case
     } else if (part->fIsAlive == true) {
       state = kFSIState;
     }
 
   } else if (part->fIsAlive == true && part->fStatus == 2 &&
       (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) {
         state = kFinalState;
 
   } else if (part->fIsAlive == true && part->fStatus == 0) {
     state = kFinalState;
 
   } else if (!part->fIsAlive && (part->fStatus == 1 || part->fStatus == 3 || part->fStatus == 4 || part->fStatus == 7 || part->fStatus == 8)) {
     state = kFSIState;
 
     // There's one hyper weird case where fStatus = -3. This apparently corresponds to a nucleon being ejected via pion FSI when there is "data available"
   } else if (!part->fIsAlive && (part->fStatus == -3)) {
     state = kUndefinedState;
     // NC neutrino outgoing
   } else if (!part->fIsAlive && part->fStatus == 0 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) {
     state = kFinalState;
 
   // Warn if we still find alive particles without classifying them
   } else if (part->fIsAlive == true) {
     ERR(WRN) << "Undefined NEUT state "
       << " Alive: " << part->fIsAlive << " Status: " << part->fStatus
       << " PDG: " << part->fPID << std::endl;
     throw;
     // Warn if we find dead particles that we haven't classified
   } else {
     ERR(WRN) << "Undefined NEUT state "
       << " Alive: " << part->fIsAlive << " Status: " << part->fStatus
       << " PDG: " << part->fPID << std::endl;
     throw;
   }
 
 
   return state;
 }
 
 void NEUTInputHandler::CalcNUISANCEKinematics() {
   // Reset all variables
   fNUISANCEEvent->ResetEvent();
 
   // Fill Globals
   fNUISANCEEvent->Mode = fNeutVect->Mode;
   fNUISANCEEvent->fEventNo = fNeutVect->EventNo;
   fNUISANCEEvent->fTargetA = fNeutVect->TargetA;
   fNUISANCEEvent->fTargetZ = fNeutVect->TargetZ;
   fNUISANCEEvent->fTargetH = fNeutVect->TargetH;
   fNUISANCEEvent->fBound = bool(fNeutVect->Ibound);
 
   if (fNUISANCEEvent->fBound) {
     fNUISANCEEvent->fTargetPDG = TargetUtils::GetTargetPDGFromZA(
         fNUISANCEEvent->fTargetZ, fNUISANCEEvent->fTargetA);
   } else {
     fNUISANCEEvent->fTargetPDG = 1000010010;
   }
 
   // Check Particle Stack
   UInt_t npart = fNeutVect->Npart();
   UInt_t kmax = fNUISANCEEvent->kMaxParticles;
   if (npart > kmax) {
     ERR(FTL) << "NEUT has too many particles. Expanding stack." << std::endl;
     fNUISANCEEvent->ExpandParticleStack(npart);
     throw;
   }
 
   int nprimary = fNeutVect->Nprimary();
   // Fill Particle Stack
   for (size_t i = 0; i < npart; i++) {
     // Get Current Count
     int curpart = fNUISANCEEvent->fNParticles;
 
     // Get NEUT Particle
     NeutPart* part = fNeutVect->PartInfo(i);
 
     // State
     int state = GetNeutParticleStatus(part);
 
     // Remove Undefined
     if (kRemoveUndefParticles && state == kUndefinedState) continue;
 
     // Remove FSI
     if (kRemoveFSIParticles && state == kFSIState) continue;
 
     // Remove Nuclear
     if (kRemoveNuclearParticles &&
         (state == kNuclearInitial || state == kNuclearRemnant))
       continue;
 
     // State
     fNUISANCEEvent->fParticleState[curpart] = state;
 
     // Is the paricle associated with the primary vertex?
     bool primary = false;
     // NEUT events are just popped onto the stack as primary, then continues to be non-primary
     if (i < nprimary) primary = true;
     fNUISANCEEvent->fPrimaryVertex[curpart] = primary;
 
     // Mom
     fNUISANCEEvent->fParticleMom[curpart][0] = part->fP.X();
     fNUISANCEEvent->fParticleMom[curpart][1] = part->fP.Y();
     fNUISANCEEvent->fParticleMom[curpart][2] = part->fP.Z();
     fNUISANCEEvent->fParticleMom[curpart][3] = part->fP.T();
 
     // PDG
     fNUISANCEEvent->fParticlePDG[curpart] = part->fPID;
 
     // Add up particle count
     fNUISANCEEvent->fNParticles++;
   }
 
   // Save Extra Generator Info
   if (fSaveExtra) {
     fNeutInfo->FillGeneratorInfo(fNeutVect);
   }
 
   // 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();
+  FitParticle* ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
+  if (ISAnyLepton) {
+    fNUISANCEEvent->probe_E = ISAnyLepton->E();
+    fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
   }
 
   return;
 }
 
 void NEUTUtils::FillNeutCommons(NeutVect* nvect) {
   // WARNING: This has only been implemented for a neuttree and not GENIE
   // This should be kept in sync with T2KNIWGUtils::GetNIWGEvent(TTree)
 
   // NEUT version info.  Can't get it to compile properly with this yet
   // neutversion_.corev  =   nvect->COREVer;
   // neutversion_.nucev  =   nvect->NUCEVer;
   // neutversion_.nuccv  =   nvect->NUCCVer;
 
   // Documentation: See nework.h
   nework_.modene = nvect->Mode;
   nework_.numne = nvect->Npart();
 
 #ifdef NEUT_COMMON_QEAV
   nemdls_.mdlqeaf = nvect->QEAVForm;
 #else
   nemdls_.mdlqeaf = nvect->QEVForm;
 #endif
   nemdls_.mdlqe = nvect->QEModel;
   nemdls_.mdlspi = nvect->SPIModel;
   nemdls_.mdldis = nvect->DISModel;
   nemdls_.mdlcoh = nvect->COHModel;
   neutcoh_.necohepi = nvect->COHModel;
 
   nemdls_.xmaqe = nvect->QEMA;
   nemdls_.xmvqe = nvect->QEMV;
   nemdls_.kapp = nvect->KAPPA;
 
   // nemdls_.sccfv = SCCFVdef;
   // nemdls_.sccfa = SCCFAdef;
   // nemdls_.fpqe = FPQEdef;
 
   nemdls_.xmaspi = nvect->SPIMA;
   nemdls_.xmvspi = nvect->SPIMV;
   nemdls_.xmares = nvect->RESMA;
   nemdls_.xmvres = nvect->RESMV;
 
   neut1pi_.xmanffres = nvect->SPIMA;
   neut1pi_.xmvnffres = nvect->SPIMV;
   neut1pi_.xmarsres = nvect->RESMA;
   neut1pi_.xmvrsres = nvect->RESMV;
   neut1pi_.neiff = nvect->SPIForm;
   neut1pi_.nenrtype = nvect->SPINRType;
   neut1pi_.rneca5i = nvect->SPICA5I;
   neut1pi_.rnebgscl = nvect->SPIBGScale;
 
   nemdls_.xmacoh = nvect->COHMA;
   nemdls_.rad0nu = nvect->COHR0;
   // nemdls_.fa1coh = nvect->COHA1err;
   // nemdls_.fb1coh = nvect->COHb1err;
 
   // neutdis_.nepdf = NEPDFdef;
   // neutdis_.nebodek = NEBODEKdef;
 
   neutcard_.nefrmflg = nvect->FrmFlg;
   neutcard_.nepauflg = nvect->PauFlg;
   neutcard_.nenefo16 = nvect->NefO16;
   neutcard_.nemodflg = nvect->ModFlg;
   // neutcard_.nenefmodl = 1;
   // neutcard_.nenefmodh = 1;
   // neutcard_.nenefkinh = 1;
   // neutpiabs_.neabspiemit = 1;
 
   nenupr_.iformlen = nvect->FormLen;
 
   neutpiless_.ipilessdcy = nvect->IPilessDcy;
   neutpiless_.rpilessdcy = nvect->RPilessDcy;
 
   neutpiless_.ipilessdcy = nvect->IPilessDcy;
   neutpiless_.rpilessdcy = nvect->RPilessDcy;
 
   neffpr_.fefqe = nvect->NuceffFactorPIQE;
   neffpr_.fefqeh = nvect->NuceffFactorPIQEH;
   neffpr_.fefinel = nvect->NuceffFactorPIInel;
   neffpr_.fefabs = nvect->NuceffFactorPIAbs;
   neffpr_.fefcx = nvect->NuceffFactorPICX;
   neffpr_.fefcxh = nvect->NuceffFactorPICXH;
 
   neffpr_.fefcoh = nvect->NuceffFactorPICoh;
   neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin;
   neffpr_.fefcxhf = nvect->NuceffFactorPICXKin;
   neffpr_.fefcohf = nvect->NuceffFactorPIQELKin;
 
   for (int i = 0; i < nework_.numne; i++) {
     nework_.ipne[i] = nvect->PartInfo(i)->fPID;
     nework_.pne[i][0] =
       (float)nvect->PartInfo(i)->fP.X() / 1000;  // VC(NE)WORK in M(G)eV
     nework_.pne[i][1] =
       (float)nvect->PartInfo(i)->fP.Y() / 1000;  // VC(NE)WORK in M(G)eV
     nework_.pne[i][2] =
       (float)nvect->PartInfo(i)->fP.Z() / 1000;  // VC(NE)WORK in M(G)eV
   }
   // fsihist.h
 
   // neutroot fills a dummy object for events with no FSI to prevent memory leak
   // when
   // reading the TTree, so check for it here
 
   if ((int)nvect->NfsiVert() ==
       1) {  // An event with FSI must have at least two vertices
     //    if (nvect->NfsiPart()!=1 || nvect->Fsiprob!=-1)
     //      ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or
     //      Fsiprob!=-1 when NfsiVert==1" << std::endl;
 
     fsihist_.nvert = 0;
     fsihist_.nvcvert = 0;
     fsihist_.fsiprob = 1;
   } else {  // Real FSI event
     fsihist_.nvert = (int)nvect->NfsiVert();
     for (int ivert = 0; ivert < fsihist_.nvert; ivert++) {
       fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID;
       fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X();
       fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y();
       fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z();
     }
 
     fsihist_.nvcvert = nvect->NfsiPart();
     for (int ip = 0; ip < fsihist_.nvcvert; ip++) {
       fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab;
       fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc;
       fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID;
       fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart;
       fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd;
       fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X();
       fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y();
       fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z();
     }
     fsihist_.fsiprob = nvect->Fsiprob;
   }
 
   neutcrscom_.crsx = nvect->Crsx;
   neutcrscom_.crsy = nvect->Crsy;
   neutcrscom_.crsz = nvect->Crsz;
   neutcrscom_.crsphi = nvect->Crsphi;
   neutcrscom_.crsq2 = nvect->Crsq2;
 
   neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ;
   neuttarget_.numbndp = nvect->TargetZ;
   neuttarget_.numfrep = nvect->TargetH;
   neuttarget_.numatom = nvect->TargetA;
   posinnuc_.ibound = nvect->Ibound;
 
   // put empty nucleon FSI history (since it is not saved in the NeutVect
   // format)
   // Comment out as NEUT does not have the necessary proton FSI information yet
   //  nucleonfsihist_.nfnvert = 0;
   //  nucleonfsihist_.nfnstep = 0;
 }
 
 #endif
diff --git a/src/InputHandler/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx
index 1853d95..d0d702c 100644
--- a/src/InputHandler/NuWroInputHandler.cxx
+++ b/src/InputHandler/NuWroInputHandler.cxx
@@ -1,507 +1,506 @@
 #ifdef __NUWRO_ENABLED__
 #include "NuWroInputHandler.h"
 #include "InputUtils.h"
 
 NuWroGeneratorInfo::~NuWroGeneratorInfo() { delete fNuWroParticlePDGs; }
 
 void NuWroGeneratorInfo::AddBranchesToTree(TTree* tn) {
   tn->Branch("NuWroParticlePDGs", &fNuWroParticlePDGs, "NuWroParticlePDGs/I");
 }
 
 void NuWroGeneratorInfo::SetBranchesFromTree(TTree* tn) {
   tn->SetBranchAddress("NuWroParticlePDGs", &fNuWroParticlePDGs);
 }
 
 void NuWroGeneratorInfo::AllocateParticleStack(int stacksize) {
   fNuWroParticlePDGs = new int[stacksize];
 }
 
 void NuWroGeneratorInfo::DeallocateParticleStack() {
   delete fNuWroParticlePDGs;
 }
 
 void NuWroGeneratorInfo::FillGeneratorInfo(event* e) { Reset(); }
 
 void NuWroGeneratorInfo::Reset() {
   for (int i = 0; i < kMaxParticles; i++) {
     fNuWroParticlePDGs[i] = 0;
   }
 }
 
 int event1_nof(event* e, int pdg) {
   int c = 0;
   for (size_t i = 0; i < e->out.size(); i++)
     if (e->out[i].pdg == pdg) c++;
   return c;
 }
 
 NuWroInputHandler::NuWroInputHandler(std::string const& handle,
                                      std::string const& rawinputs) {
   LOG(SAM) << "Creating NuWroInputHandler : " << handle << std::endl;
 
   // Run a joint input handling
   fName = handle;
   fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
   fSaveExtra = false;  // FitPar::Config().GetParB("NuWroSaveExtra");
   // Setup the TChain
   fNuWroTree = new TChain("treeout");
 
   // 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(inputs[inp_it].c_str(), "READ");
     if (!inp_file or inp_file->IsZombie()) {
       ERR(FTL) << "nuwro File IsZombie() at " << inputs[inp_it] << std::endl;
       throw;
     }
 
     // Get Flux/Event hist
     TH1D* fluxhist = (TH1D*)inp_file->Get(
         (PlotUtils::GetObjectWithName(inp_file, "FluxHist")).c_str());
     TH1D* eventhist = (TH1D*)inp_file->Get(
         (PlotUtils::GetObjectWithName(inp_file, "EvtHist")).c_str());
     if (!fluxhist or !eventhist) {
       ERR(FTL) << "nuwro FILE doesn't contain flux/xsec info" << std::endl;
       if (FitPar::Config().GetParB("regennuwro")) {
         ERR(FTL) << "Regen NuWro has not been added yet. Email the developers!"
                  << std::endl;
         // ProcessNuWroInputFlux(inputs[inp_it]);
         throw;
       } else {
         ERR(FTL) << "If you would like NUISANCE to generate these for you "
                  << "please set parameter regennuwro=1 and re-run."
                  << std::endl;
         throw;
       }
     }
 
     // Get N Events
     TTree* nuwrotree = (TTree*)inp_file->Get("treeout");
     if (!nuwrotree) {
       ERR(FTL) << "treeout not located in nuwro file! " << inputs[inp_it]
                << std::endl;
       throw;
     }
     int nevents = nuwrotree->GetEntries();
 
     // Register input to form flux/event rate hists
     RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
 
     // Add to TChain
     fNuWroTree->Add(inputs[inp_it].c_str());
   }
 
   // Registor all our file inputs
   SetupJointInputs();
 
   // Setup Events
   fNuWroEvent = NULL;
   fNuWroTree->SetBranchAddress("e", &fNuWroEvent);
   fNuWroTree->GetEntry(0);
 
   fNUISANCEEvent = new FitEvent();
   fNUISANCEEvent->fType = kNUWRO;
   fNUISANCEEvent->fNuwroEvent = fNuWroEvent;
 
   fNUISANCEEvent->HardReset();
 
   if (fSaveExtra) {
     fNuWroInfo = new NuWroGeneratorInfo();
     fNUISANCEEvent->AddGeneratorInfo(fNuWroInfo);
   }
 };
 
 NuWroInputHandler::~NuWroInputHandler() {
   if (fNuWroTree) delete fNuWroTree;
 }
 
 void NuWroInputHandler::CreateCache() {
   // fNuWroTree->SetCacheEntryRange(0, fNEvents);
   //    fNuWroTree->AddBranchToCache("*", 1);
   //    fNuWroTree->SetCacheSize(fCacheSize);
 }
 
 void NuWroInputHandler::RemoveCache() {
   // fNuWroTree->SetCacheEntryRange(0, fNEvents);
   //    fNuWroTree->AddBranchToCache("*", 0);
   //    fNuWroTree->SetCacheSize(0);
 }
 
 void NuWroInputHandler::ProcessNuWroInputFlux(const std::string file) {}
 
 FitEvent* NuWroInputHandler::GetNuisanceEvent(const UInt_t entry,
                                               const bool lightweight) {
   // Catch too large entries
   if (entry >= (UInt_t)fNEvents) return NULL;
 
   // Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
   fNuWroTree->GetEntry(entry);
 
   // Run NUISANCE Vector Filler
   if (!lightweight) {
     CalcNUISANCEKinematics();
   }
 #ifdef __PROB3PP_ENABLED__
   for (size_t i = 0; i < fNUISANCEEvent->fNuwroEvent->in.size(); i++) {
     if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
                    fNUISANCEEvent->fNuwroEvent->in[i].pdg)) {
       fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroEvent->in[i].t;
       fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroEvent->in[i].pdg;
       break;
     }
   }
 #endif
   // Setup Input scaling for joint inputs
   fNUISANCEEvent->InputWeight = GetInputWeight(entry);
 
 #ifdef __USE_NUWRO_SRW_EVENTS__
   if (!rwEvs.size()) {
     fNuwroParams = fNuWroEvent->par;
   }
 
   if (entry >= rwEvs.size()) {
     rwEvs.push_back(BaseFitEvt());
     rwEvs.back().fType = kNUWRO;
     rwEvs.back().Mode = fNUISANCEEvent->Mode;
     rwEvs.back().fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
     rwEvs.back().fNuwroEvent = NULL;
     rwEvs.back().fNuwroParams = &fNuwroParams;
     rwEvs.back().probe_E = rwEvs.back().fNuwroSRWEvent.NeutrinoEnergy;
     rwEvs.back().probe_pdg = rwEvs.back().fNuwroSRWEvent.NeutrinoPDG;
   }
 
   fNUISANCEEvent->fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
   fNUISANCEEvent->fNuwroParams = &fNuwroParams;
   fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoEnergy;
   fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoPDG;
 #endif
 
   return fNUISANCEEvent;
 }
 
 int NuWroInputHandler::ConvertNuwroMode(event* e) {
   Int_t proton_pdg, neutron_pdg, pion_pdg, pion_plus_pdg, pion_minus_pdg,
       lambda_pdg, eta_pdg, kaon_pdg, kaon_plus_pdg;
   proton_pdg = 2212;
   eta_pdg = 221;
   neutron_pdg = 2112;
   pion_pdg = 111;
   pion_plus_pdg = 211;
   pion_minus_pdg = -211;
   // O_16_pdg = 100069;   // oznacznie z Neuta
   lambda_pdg = 3122;
   kaon_pdg = 311;
   kaon_plus_pdg = 321;
 
   if (e->flag.qel)  // kwiazielastyczne oddziaływanie
   {
     if (e->flag.anty)  // jeśli jest to oddziaływanie z antyneutrinem
     {
       if (e->flag.cc)
         return -1;
       else {
         if (event1_nof(e, proton_pdg))
           return -51;
         else if (event1_nof(e, neutron_pdg))
           return -52;  // sprawdzam dodatkowo ?
       }
     } else  // oddziaływanie z neutrinem
     {
       if (e->flag.cc)
         return 1;
       else {
         if (event1_nof(e, proton_pdg))
           return 51;
         else if (event1_nof(e, neutron_pdg))
           return 52;
       }
     }
   }
 
   if (e->flag.mec) {
     if (e->flag.anty)
       return -2;
     else
       return 2;
   }
 
   if (e->flag.res)  // rezonansowa produkcja: pojedynczy pion, pojed.eta, kaon,
                     // multipiony
   {
     Int_t liczba_pionow, liczba_kaonow;
 
     liczba_pionow = event1_nof(e, pion_pdg) + event1_nof(e, pion_plus_pdg) +
                     event1_nof(e, pion_minus_pdg);
     liczba_kaonow = event1_nof(e, kaon_pdg) + event1_nof(e, kaon_pdg);
 
     if (liczba_pionow > 1 || liczba_pionow == 0)  // multipiony
     {
       if (e->flag.anty) {
         if (e->flag.cc)
           return -21;
         else
           return -41;
       } else {
         if (e->flag.cc)
           return 21;
         else
           return 41;
       }
     }
 
     if (liczba_pionow == 1) {
       if (e->flag.anty)  // jeśli jest to oddziaływanie z antyneutrinem
       {
         if (e->flag.cc) {
           if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_minus_pdg))
             return -11;
           if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_pdg)) return -12;
           if (event1_nof(e, proton_pdg) && event1_nof(e, pion_minus_pdg))
             return -13;
         } else {
           if (event1_nof(e, proton_pdg)) {
             if (event1_nof(e, pion_minus_pdg))
               return -33;
             else if (event1_nof(e, pion_pdg))
               return -32;
           } else if (event1_nof(e, neutron_pdg)) {
             if (event1_nof(e, pion_plus_pdg))
               return -34;
             else if (event1_nof(e, pion_pdg))
               return -31;
           }
         }
       } else  // oddziaływanie z neutrinem
       {
         if (e->flag.cc) {
           if (event1_nof(e, proton_pdg) && event1_nof(e, pion_plus_pdg))
             return 11;
           if (event1_nof(e, proton_pdg) && event1_nof(e, pion_pdg)) return 12;
           if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_plus_pdg))
             return 13;
         } else {
           if (event1_nof(e, proton_pdg)) {
             if (event1_nof(e, pion_minus_pdg))
               return 33;
             else if (event1_nof(e, pion_pdg))
               return 32;
           } else if (event1_nof(e, neutron_pdg)) {
             if (event1_nof(e, pion_plus_pdg))
               return 34;
             else if (event1_nof(e, pion_pdg))
               return 31;
           }
         }
       }
     }
 
     if (event1_nof(e, eta_pdg))  // produkcja rezonansowa ety
     {
       if (e->flag.anty)  // jeśli jest to oddziaływanie z antyneutrinem
       {
         if (e->flag.cc)
           return -22;
         else {
           if (event1_nof(e, neutron_pdg))
             return -42;
           else if (event1_nof(e, proton_pdg))
             return -43;  // sprawdzam dodatkowo ?
         }
       } else  // oddziaływanie z neutrinem
       {
         if (e->flag.cc)
           return 22;
         else {
           if (event1_nof(e, neutron_pdg))
             return 42;
           else if (event1_nof(e, proton_pdg))
             return 43;
         }
       }
     }
 
     if (event1_nof(e, lambda_pdg) == 1 &&
         liczba_kaonow == 1)  // produkcja rezonansowa kaonu
     {
       if (e->flag.anty)  // jeśli jest to oddziaływanie z antyneutrinem
       {
         if (e->flag.cc && event1_nof(e, kaon_pdg))
           return -23;
         else {
           if (event1_nof(e, kaon_pdg))
             return -44;
           else if (event1_nof(e, kaon_plus_pdg))
             return -45;
         }
       } else  // oddziaływanie z neutrinem
       {
         if (e->flag.cc && event1_nof(e, kaon_plus_pdg))
           return 23;
         else {
           if (event1_nof(e, kaon_pdg))
             return 44;
           else if (event1_nof(e, kaon_plus_pdg))
             return 45;
         }
       }
     }
   }
 
   if (e->flag.coh)  // koherentne  oddziaływanie tylko na O(16)
   {
     Int_t _target;
     _target = e->par.nucleus_p + e->par.nucleus_n;  // liczba masowa  O(16)
 
     if (_target == 16) {
       if (e->flag.anty)  // jeśli jest to oddziaływanie z antyneutrinem
       {
         if (e->flag.cc && event1_nof(e, pion_minus_pdg))
           return -16;
         else if (event1_nof(e, pion_pdg))
           return -36;
       } else  // oddziaływanie z neutrinem
       {
         if (e->flag.cc && event1_nof(e, pion_plus_pdg))
           return 16;
         else if (event1_nof(e, pion_pdg))
           return 36;
       }
     }
   }
 
   // gleboko nieelastyczne rozpraszanie
   if (e->flag.dis) {
     if (e->flag.anty) {
       if (e->flag.cc)
         return -26;
       else
         return -46;
     } else {
       if (e->flag.cc)
         return 26;
       else
         return 46;
     }
   }
 
   return 9999;
 }
 
 void NuWroInputHandler::CalcNUISANCEKinematics() {
   // std::cout << "NuWro Event Address " << fNuWroEvent << std::endl;
   // Reset all variables
   fNUISANCEEvent->ResetEvent();
   FitEvent* evt = fNUISANCEEvent;
 
   // Sort Event Info
   evt->Mode = ConvertNuwroMode(fNuWroEvent);
 
   if (abs(evt->Mode) > 60) {
     evt->Mode = 0;
   }
 
   evt->fEventNo = 0.0;
   evt->fTotCrs = 0.0;
   evt->fTargetA = fNuWroEvent->par.nucleus_p + fNuWroEvent->par.nucleus_n;
   evt->fTargetZ = fNuWroEvent->par.nucleus_p;
   evt->fTargetPDG = TargetUtils::GetTargetPDGFromZA(evt->fTargetZ, evt->fTargetA);
   evt->fTargetH = 0;
   evt->fBound = (evt->fTargetA != 1);
 
   // Check Particle Stack
   UInt_t npart_in = fNuWroEvent->in.size();
   UInt_t npart_out = fNuWroEvent->out.size();
   UInt_t npart_post = fNuWroEvent->post.size();
   UInt_t npart = npart_in + npart_out + npart_post;
   UInt_t kmax = evt->kMaxParticles;
 
   if (npart > kmax) {
     ERR(WRN) << "NUWRO has too many particles. Expanding stack." << std::endl;
     fNUISANCEEvent->ExpandParticleStack(npart);
   }
 
   evt->fNParticles = 0;
   std::vector<particle>::iterator p_iter;
 
   // Get the Initial State
   for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end(); p_iter++) {
     AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState, true);
   }
 
   // Try to find the FSI state particles
   // Loop over the primary vertex particles
   // If they match the post-FSI they haven't undergone FSI.
   // If they don't match post-FSI they have undergone FSI.
   for (p_iter = fNuWroEvent->out.begin(); p_iter != fNuWroEvent->out.end(); p_iter++) {
     // Get the particle
     particle p = (*p_iter);
     // Check against all the post particles, match them
     std::vector<particle>::iterator p2_iter;
     bool match = false;
     for (p2_iter = fNuWroEvent->post.begin(); p2_iter != fNuWroEvent->post.end(); p2_iter++) {
       particle p2 = (*p2_iter);
       // Check energy and pdg
       // A very small cascade which changes the energy by 1E-5 MeV should be matched
       match = (fabs(p2.E()-p.E()) < 1E-5 && p2.pdg == p.pdg);
       // If we match p to p2 break the loop
       if (match) break;
     }
     // If we've looped through the whole particle stack of post-FSI and haven't found a match it's a primary particle that has been FSIed
     if (!match) AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState, true);
   }
 
   // Loop over the final state particles
   for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end(); p_iter++) {
     particle p = (*p_iter);
     // To find if it's primary or not we have to loop through the primary ones and match, just like above
     bool match = false;
     std::vector<particle>::iterator p2_iter;
     for (p2_iter = fNuWroEvent->out.begin(); p2_iter != fNuWroEvent->out.end(); p2_iter++) {
       particle p2 = (*p2_iter);
       match = (fabs(p2.E()-p.E()) < 1E-5 && p2.pdg == p.pdg);
       if (match) break;
     }
     AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState, match);
   }
 
   // Fill Generator Info
   if (fSaveExtra) fNuWroInfo->FillGeneratorInfo(fNuWroEvent);
 
   // 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();
+  FitParticle* ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
+  if (ISAnyLepton) {
+    fNUISANCEEvent->probe_E = ISAnyLepton->E();
+    fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
   }
 
   return;
 }
 
 void NuWroInputHandler::AddNuWroParticle(FitEvent* evt, particle& p,
                                          int state, bool primary = false) {
   // Add Mom
   evt->fParticleMom[evt->fNParticles][0] = static_cast<vect&>(p).x;
   evt->fParticleMom[evt->fNParticles][1] = static_cast<vect&>(p).y;
   evt->fParticleMom[evt->fNParticles][2] = static_cast<vect&>(p).z;
   evt->fParticleMom[evt->fNParticles][3] = static_cast<vect&>(p).t;
 
   // For NuWro a particle that we've given a FSI state is a pre-FSI particle
   // An initial state particle is also a primary vertex praticle
   evt->fPrimaryVertex[evt->fNParticles] = primary;
 
   // Status/PDG
   evt->fParticleState[evt->fNParticles] = state;
   evt->fParticlePDG[evt->fNParticles] = p.pdg;
 
   // Add to particle count
   evt->fNParticles++;
 }
 
 void NuWroInputHandler::Print() {}
 
 #endif