No OneTemporary
Actions

Size

62 KB

Subscribers

None

View Options

	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

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 4:02 PM (1 d, 21 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4022569
Default Alt Text: (62 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions