No OneTemporary
Actions

Size

82 KB

Subscribers

None

View Options

	diff --git a/src/InputHandler/FitEvent.h b/src/InputHandler/FitEvent.h
	index ddc5f99..f0b3db6 100644
	--- a/src/InputHandler/FitEvent.h
	+++ b/src/InputHandler/FitEvent.h
	@@ -1,636 +1,637 @@
	// 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/>.
	*******************************************************************************/
	#ifndef FITEVENT2_H_SEEN
	#define FITEVENT2_H_SEEN
	/*!
	* \addtogroup InputHandler
	* @{
	*/

	#include <algorithm>
	#include <iterator>
	#include <vector>
	#include "FitParticle.h"
	#include "TLorentzVector.h"
	#include "TSpline.h"

	#include "BaseFitEvt.h"
	#include "FitLogger.h"
	#include "TArrayD.h"
	#include "TTree.h"
	#include "TChain.h"
	+#include "TargetUtils.h"

	#include "PhysConst.h"

	/// Common container for event particles
	class FitEvent : public BaseFitEvt {
	public:

	FitEvent();
	~FitEvent() {};

	void FreeFitParticles();
	void ClearFitParticles();
	void ResetEvent(void);
	void ResetParticleList(void);
	void HardReset();
	void OrderStack();
	void SetBranchAddress(TChain* tn);
	void AddBranchesToTree(TTree* tn);
	void Print();
	void PrintChris();
	void DeallocateParticleStack();
	void AllocateParticleStack(int stacksize);
	void ExpandParticleStack(int stacksize);
	void AddGeneratorInfo(GeneratorInfoBase* gen);


	// ---- HELPER/ACCESS FUNCTIONS ---- //
	/// Return True Interaction ID
	inline int GetMode (void) const { return Mode; };
	/// Return Target Atomic Number
	inline int GetTargetA (void) const { return fTargetA; };
	/// Return Target Nuclear Charge
	inline int GetTargetZ (void) const { return fTargetZ; };
	/// Get Event Total Cross-section
	inline int GetTotCrs (void) const { return fTotCrs; };

	/// Is Event Charged Current?
	inline bool IsCC() const { return (abs(Mode) <= 30); };
	/// Is Event Neutral Current?
	inline bool IsNC() const { return (abs(Mode) > 30); };

	/// Return Particle 4-momentum for given index in particle stack
	TLorentzVector GetParticleP4 (int index) const;
	/// Return Particle 3-momentum for given index in particle stack
	TVector3 GetParticleP3 (int index) const;
	/// Return Particle absolute momentum for given index in particle stack
	double GetParticleMom (int index) const;
	/// Return Particle absolute momentum-squared for given index in particle stack
	double GetParticleMom2 (int index) const;
	/// Return Particle energy for given index in particle stack
	double GetParticleE (int index) const;
	/// Return Particle State for given index in particle stack
	int GetParticleState (int index) const;
	/// Return Particle PDG for given index in particle stack
	int GetParticlePDG (int index) const;


	/// Allows the removal of KE up to total KE.
	inline void RemoveKE(int index, double KE){

	FitParticle *fp = GetParticle(index);

	double mass = fp->M();
	double oKE = fp->KE();
	double nE = mass + (oKE - KE);
	if(nE < mass){ // Can't take more KE than it has
	nE = mass;
	}
	double n3Mom = sqrt(nEnE - massmass);
	TVector3 np3 = fp->P3().Unit()*n3Mom;

	fParticleMom[index][0] = np3[0];
	fParticleMom[index][1] = np3[1];
	fParticleMom[index][2] = np3[2];
	fParticleMom[index][3] = nE;

	}

	/// Allows the removal of KE up to total KE.
	inline void GiveKE(int index, double KE){
	RemoveKE(index,-KE);
	}
	/// Return Particle for given index in particle stack
	FitParticle* GetParticle(int const index);
	/// Get Total Number of Particles in stack
	inline uint NParticles (void) const { return fNParticles; };

	/// Check if event contains a particle given a pdg and state.
	/// If no state is passed all states are considered.
	bool HasParticle (int const pdg = 0, int const state = -1) const ;

	template <size_t N>
	inline bool HasParticle(int const (&pdgs)[N], int const state = -1) const {
	for (size_t i = 0; i < N; i++) {
	if (HasParticle( pdgs[i], state )) {
	return false;
	}
	}
	return false;
	}

	/// Get total number of particles given a pdg and state.
	/// If no state is passed all states are considered.
	int NumParticle (int const pdg = 0, int const state = -1) const;

	template <size_t N>
	inline int NumParticle(int const (&pdgs)[N], int const state = -1) const {
	int ncount = 0;
	for (size_t i = 0; i < N; i++) {
	ncount += NumParticle( pdgs[i], state );
	}
	return ncount;
	}

	/// Return a vector of particle indices that can be used with 'GetParticle'
	/// functions given a particle pdg and state.
	/// If no state is passed all states are considered.
	std::vector<int> GetAllParticleIndices (int const pdg = -1, int const state = -1) const;

	template <size_t N>
	inline std::vector<int> GetAllParticleIndices(int const (&pdgs)[N], const int state = -1) const {
	std::vector<int> plist;
	for (size_t i = 0; i < N; i++) {
	std::vector<int> plisttemp = GetAllParticleIndices(pdgs[i], state);
	plist.insert( plist.end(), plisttemp.begin(), plisttemp.end() );
	}
	return plist;
	}

	/// Return a vector of FitParticles given a particle pdg and state.
	/// This is memory intensive and slow than GetAllParticleIndices,
	/// but is slightly easier to use.
	std::vector<FitParticle*> GetAllParticle (int const pdg = -1, int const state = -1);

	template <size_t N>
	inline std::vector<FitParticle*> GetAllParticle(int const (&pdgs)[N], int const state = -1) {
	std::vector<FitParticle*> plist;
	for (size_t i = 0; i < N; i++) {
	std::vector<FitParticle*> plisttemp = GetAllParticle(pdgs[i], state);
	plist.insert( plist.end(), plisttemp.begin(), plisttemp.end() );
	}
	return plist;
	}

	inline std::vector<int> GetAllNuElectronIndices (void) { return GetAllParticleIndices(12); };
	inline std::vector<int> GetAllNuMuonIndices (void) { return GetAllParticleIndices(14); };
	inline std::vector<int> GetAllNuTauIndices (void) { return GetAllParticleIndices(16); };
	inline std::vector<int> GetAllElectronIndices (void) { return GetAllParticleIndices(11); };
	inline std::vector<int> GetAllMuonIndices (void) { return GetAllParticleIndices(13); };
	inline std::vector<int> GetAllTauIndices (void) { return GetAllParticleIndices(15); };
	inline std::vector<int> GetAllProtonIndices (void) { return GetAllParticleIndices(2212); };
	inline std::vector<int> GetAllNeutronIndices (void) { return GetAllParticleIndices(2112); };
	inline std::vector<int> GetAllPiZeroIndices (void) { return GetAllParticleIndices(111); };
	inline std::vector<int> GetAllPiPlusIndices (void) { return GetAllParticleIndices(211); };
	inline std::vector<int> GetAllPiMinusIndices (void) { return GetAllParticleIndices(-211); };
	inline std::vector<int> GetAllPhotonIndices (void) { return GetAllParticleIndices(22); };

	inline std::vector<FitParticle*> GetAllNuElectron (void) { return GetAllParticle(12); };
	inline std::vector<FitParticle*> GetAllNuMuon (void) { return GetAllParticle(14); };
	inline std::vector<FitParticle*> GetAllNuTau (void) { return GetAllParticle(16); };
	inline std::vector<FitParticle*> GetAllElectron (void) { return GetAllParticle(11); };
	inline std::vector<FitParticle*> GetAllMuon (void) { return GetAllParticle(13); };
	inline std::vector<FitParticle*> GetAllTau (void) { return GetAllParticle(15); };
	inline std::vector<FitParticle*> GetAllProton (void) { return GetAllParticle(2212); };
	inline std::vector<FitParticle*> GetAllNeutron (void) { return GetAllParticle(2112); };
	inline std::vector<FitParticle*> GetAllPiZero (void) { return GetAllParticle(111); };
	inline std::vector<FitParticle*> GetAllPiPlus (void) { return GetAllParticle(211); };
	inline std::vector<FitParticle*> GetAllPiMinus (void) { return GetAllParticle(-211); };
	inline std::vector<FitParticle*> GetAllPhoton (void) { return GetAllParticle(22); };

	// --- Highest Momentum Search Functions --- //
	/// Returns the Index of the highest momentum particle given a pdg and state.
	/// If no state is given all states are considered, but that will just return the
	/// momentum of the beam in most cases so is not advised.
	int GetHMParticleIndex (int const pdg = 0, int const state = -1) const;

	template <size_t N>
	inline int GetHMParticleIndex (int const (&pdgs)[N], int const state = -1) const {

	double mom = -999.9;
	int rtnindex = -1;

	for (size_t i = 0; i < N; ++i) {
	// Use ParticleMom as doesn't require Particle Mem alloc
	int pindex = GetHMParticleIndex(pdgs[i], state);
	if (pindex != -1){
	double momnew = GetParticleMom2(pindex);
	if (momnew > mom) {
	rtnindex = pindex;
	mom = momnew;
	}
	}
	}
	return rtnindex;
	};

	/// Returns the highest momentum particle given a pdg and state.
	/// If no state is given all states are considered, but that will just return the
	/// momentum of the beam in most cases so is not advised.
	inline FitParticle* GetHMParticle(int const pdg = 0, int const state = -1) {
	return GetParticle( GetHMParticleIndex(pdg, state) );
	}

	template <size_t N>
	inline FitParticle* GetHMParticle(int const (&pdgs)[N], int const state) {
	return GetParticle(GetHMParticleIndex(pdgs, state));
	};


	// ---- Initial State Helpers --- //
	/// Checks the event has a particle of a given pdg in the initial state.
	inline bool HasISParticle(int const pdg) const {
	return HasParticle(pdg, kInitialState);
	};
	template <size_t N>
	inline bool HasISParticle(int const (&pdgs)[N]) const {
	return HasParticle(pdgs, kInitialState);
	};

	/// Returns the number of particles with a given pdg in the initial state.
	inline int NumISParticle(int const pdg = 0) const {
	return NumParticle(pdg, kInitialState);
	};
	template <size_t N>
	inline int NumISParticle(int const (&pdgs)[N]) const {
	return NumParticle(pdgs, kInitialState);
	};

	/// Returns a list of indices for all particles with a given pdg
	/// in the initial state. These can be used with the 'GetParticle' functions.
	inline std::vector<int> GetAllISParticleIndices(int const pdg = -1) const {
	return GetAllParticleIndices(pdg, kInitialState);
	};
	template <size_t N>
	inline std::vector<int> GetAllISParticleIndices(int const (&pdgs)[N]) const {
	return GetAllParticleIndices(pdgs, kInitialState);
	};

	/// Returns a list of particles with a given pdg in the initial state.
	/// This function is more memory intensive and slower than the Indices function.
	inline std::vector<FitParticle*> GetAllISParticle(int const pdg = -1) {
	return GetAllParticle(pdg, kInitialState);
	};
	template <size_t N>
	inline std::vector<FitParticle*> GetAllISParticle(int const (&pdgs)[N]) {
	return GetAllParticle(pdgs, kInitialState);
	};

	/// Returns the highest momentum particle with a given pdg in the initial state.
	inline FitParticle* GetHMISParticle(int const pdg) {
	return GetHMParticle(pdg, kInitialState);
	};
	template <size_t N>
	inline FitParticle* GetHMISParticle(int const (&pdgs)[N]) {
	return GetHMParticle(pdgs, kInitialState);
	};

	/// Returns the highest momentum particle index with a given pdg in the initial state.
	inline int GetHMISParticleIndex(int const pdg) const {
	return GetHMParticleIndex(pdg, kInitialState);
	};
	template <size_t N>
	inline int GetHMISParticleIndex(int const (&pdgs)[N]) const {
	return GetHMParticleIndex(pdgs, kInitialState);
	};

	inline bool HasISNuElectron (void) const { return HasISParticle(12); };
	inline bool HasISNuMuon (void) const { return HasISParticle(14); };
	inline bool HasISNuTau (void) const { return HasISParticle(16); };
	inline bool HasISElectron (void) const { return HasISParticle(11); };
	inline bool HasISMuon (void) const { return HasISParticle(13); };
	inline bool HasISTau (void) const { return HasISParticle(15); };
	inline bool HasISProton (void) const { return HasISParticle(2212); };
	inline bool HasISNeutron (void) const { return HasISParticle(2112); };
	inline bool HasISPiZero (void) const { return HasISParticle(111); };
	inline bool HasISPiPlus (void) const { return HasISParticle(211); };
	inline bool HasISPiMinus (void) const { return HasISParticle(-211); };
	inline bool HasISPhoton (void) const { return HasISParticle(22); };
	inline bool HasISLeptons (void) const { return HasISParticle(PhysConst::pdg_leptons); };
	inline bool HasISPions (void) const { return HasISParticle(PhysConst::pdg_pions); };
	inline bool HasISChargePions (void) const { return HasISParticle(PhysConst::pdg_charged_pions); };

	inline int NumISNuElectron (void) const { return NumISParticle(12); };
	inline int NumISNuMuon (void) const { return NumISParticle(14); };
	inline int NumISNuTau (void) const { return NumISParticle(16); };
	inline int NumISElectron (void) const { return NumISParticle(11); };
	inline int NumISMuon (void) const { return NumISParticle(13); };
	inline int NumISTau (void) const { return NumISParticle(15); };
	inline int NumISProton (void) const { return NumISParticle(2212); };
	inline int NumISNeutron (void) const { return NumISParticle(2112); };
	inline int NumISPiZero (void) const { return NumISParticle(111); };
	inline int NumISPiPlus (void) const { return NumISParticle(211); };
	inline int NumISPiMinus (void) const { return NumISParticle(-211); };
	inline int NumISPhoton (void) const { return NumISParticle(22); };
	inline int NumISLeptons (void) const { return NumISParticle(PhysConst::pdg_leptons); };
	inline int NumISPions (void) const { return NumISParticle(PhysConst::pdg_pions); };
	inline int NumISChargePions (void) const { return NumISParticle(PhysConst::pdg_charged_pions); };

	inline std::vector<int> GetAllISNuElectronIndices (void) const { return GetAllISParticleIndices(12); };
	inline std::vector<int> GetAllISNuMuonIndices (void) const { return GetAllISParticleIndices(14); };
	inline std::vector<int> GetAllISNuTauIndices (void) const { return GetAllISParticleIndices(16); };
	inline std::vector<int> GetAllISElectronIndices (void) const { return GetAllISParticleIndices(11); };
	inline std::vector<int> GetAllISMuonIndices (void) const { return GetAllISParticleIndices(13); };
	inline std::vector<int> GetAllISTauIndices (void) const { return GetAllISParticleIndices(15); };
	inline std::vector<int> GetAllISProtonIndices (void) const { return GetAllISParticleIndices(2212); };
	inline std::vector<int> GetAllISNeutronIndices (void) const { return GetAllISParticleIndices(2112); };
	inline std::vector<int> GetAllISPiZeroIndices (void) const { return GetAllISParticleIndices(111); };
	inline std::vector<int> GetAllISPiPlusIndices (void) const { return GetAllISParticleIndices(211); };
	inline std::vector<int> GetAllISPiMinusIndices (void) const { return GetAllISParticleIndices(-211); };
	inline std::vector<int> GetAllISPhotonIndices (void) const { return GetAllISParticleIndices(22); };
	inline std::vector<int> GetAllISLeptonsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_leptons); };
	inline std::vector<int> GetAllISPionsIndices (void) const { return GetAllISParticleIndices(PhysConst::pdg_pions); };
	inline std::vector<int> GetAllISChargePionsIndices(void) const { return GetAllISParticleIndices(PhysConst::pdg_charged_pions); };

	inline std::vector<FitParticle*> GetAllISNuElectron (void) { return GetAllISParticle(12); };
	inline std::vector<FitParticle*> GetAllISNuMuon (void) { return GetAllISParticle(14); };
	inline std::vector<FitParticle*> GetAllISNuTau (void) { return GetAllISParticle(16); };
	inline std::vector<FitParticle*> GetAllISElectron (void) { return GetAllISParticle(11); };
	inline std::vector<FitParticle*> GetAllISMuon (void) { return GetAllISParticle(13); };
	inline std::vector<FitParticle*> GetAllISTau (void) { return GetAllISParticle(15); };
	inline std::vector<FitParticle*> GetAllISProton (void) { return GetAllISParticle(2212); };
	inline std::vector<FitParticle*> GetAllISNeutron (void) { return GetAllISParticle(2112); };
	inline std::vector<FitParticle*> GetAllISPiZero (void) { return GetAllISParticle(111); };
	inline std::vector<FitParticle*> GetAllISPiPlus (void) { return GetAllISParticle(211); };
	inline std::vector<FitParticle*> GetAllISPiMinus (void) { return GetAllISParticle(-211); };
	inline std::vector<FitParticle*> GetAllISPhoton (void) { return GetAllISParticle(22); };
	inline std::vector<FitParticle*> GetAllISLeptons (void) { return GetAllISParticle(PhysConst::pdg_leptons); };
	inline std::vector<FitParticle*> GetAllISPions (void) { return GetAllISParticle(PhysConst::pdg_pions); };
	inline std::vector<FitParticle*> GetAllISChargePions(void) { return GetAllISParticle(PhysConst::pdg_charged_pions); };

	inline FitParticle* GetHMISNuElectron (void) { return GetHMISParticle(12); };
	inline FitParticle* GetHMISNuMuon (void) { return GetHMISParticle(14); };
	inline FitParticle* GetHMISNuTau (void) { return GetHMISParticle(16); };
	inline FitParticle* GetHMISElectron (void) { return GetHMISParticle(11); };
	inline FitParticle* GetHMISMuon (void) { return GetHMISParticle(13); };
	inline FitParticle* GetHMISTau (void) { return GetHMISParticle(15); };
	inline FitParticle* GetHMISAnyLeptons (void) { return GetHMISParticle(PhysConst::pdg_all_leptons); };
	inline FitParticle* GetHMISProton (void) { return GetHMISParticle(2212); };
	inline FitParticle* GetHMISNeutron (void) { return GetHMISParticle(2112); };
	inline FitParticle* GetHMISPiZero (void) { return GetHMISParticle(111); };
	inline FitParticle* GetHMISPiPlus (void) { return GetHMISParticle(211); };
	inline FitParticle* GetHMISPiMinus (void) { return GetHMISParticle(-211); };
	inline FitParticle* GetHMISPhoton (void) { return GetHMISParticle(22); };
	inline FitParticle* GetHMISLepton (void) { return GetHMISParticle(PhysConst::pdg_leptons); };
	inline FitParticle* GetHMISPions (void) { return GetHMISParticle(PhysConst::pdg_pions); };
	inline FitParticle* GetHMISChargePions(void) { return GetHMISParticle(PhysConst::pdg_charged_pions); };

	inline int GetHMISNuElectronIndex (void) { return GetHMISParticleIndex(12); };
	inline int GetHMISNuMuonIndex (void) { return GetHMISParticleIndex(14); };
	inline int GetHMISNuTauIndex (void) { return GetHMISParticleIndex(16); };
	inline int GetHMISElectronIndex (void) { return GetHMISParticleIndex(11); };
	inline int GetHMISMuonIndex (void) { return GetHMISParticleIndex(13); };
	inline int GetHMISTauIndex (void) { return GetHMISParticleIndex(15); };
	inline int GetHMISProtonIndex (void) { return GetHMISParticleIndex(2212); };
	inline int GetHMISNeutronIndex (void) { return GetHMISParticleIndex(2112); };
	inline int GetHMISPiZeroIndex (void) { return GetHMISParticleIndex(111); };
	inline int GetHMISPiPlusIndex (void) { return GetHMISParticleIndex(211); };
	inline int GetHMISPiMinusIndex (void) { return GetHMISParticleIndex(-211); };
	inline int GetHMISPhotonIndex (void) { return GetHMISParticleIndex(22); };
	inline int GetHMISLeptonsIndex (void) { return GetHMISParticleIndex(PhysConst::pdg_leptons); };
	inline int GetHMISPionsIndex (void) { return GetHMISParticleIndex(PhysConst::pdg_pions); };
	inline int GetHMISChargePionsIndex(void) { return GetHMISParticleIndex(PhysConst::pdg_charged_pions); };

	// ---- Final State Helpers --- //
	inline bool HasFSParticle(int const pdg) const {
	return HasParticle(pdg, kFinalState);
	};
	template <size_t N>
	inline bool HasFSParticle(int const (&pdgs)[N]) const {
	return HasParticle(pdgs, kFinalState);
	};

	inline int NumFSParticle(int const pdg = 0) const {
	return NumParticle(pdg, kFinalState);
	};
	template <size_t N>
	inline int NumFSParticle(int const (&pdgs)[N]) const {
	return NumParticle(pdgs, kFinalState);
	};

	inline std::vector<int> GetAllFSParticleIndices(int const pdg = -1) const {
	return GetAllParticleIndices(pdg, kFinalState);
	};
	template <size_t N>
	inline std::vector<int> GetAllFSParticleIndices(int const (&pdgs)[N]) const {
	return GetAllParticleIndices(pdgs, kFinalState);
	};

	inline std::vector<FitParticle*> GetAllFSParticle(int const pdg = -1) {
	return GetAllParticle(pdg, kFinalState);
	};
	template <size_t N>
	inline std::vector<FitParticle*> GetAllFSParticle(int const (&pdgs)[N]) {
	return GetAllParticle(pdgs, kFinalState);
	};

	inline FitParticle* GetHMFSParticle(int const pdg) {
	return GetHMParticle(pdg, kFinalState);
	};
	template <size_t N>
	inline FitParticle* GetHMFSParticle(int const (&pdgs)[N]) {
	return GetHMParticle(pdgs, kFinalState);
	};

	inline int GetHMFSParticleIndex(int const pdg) const {
	return GetHMParticleIndex(pdg, kFinalState);
	};
	template <size_t N>
	inline int GetHMFSParticleIndex(int const (&pdgs)[N]) const {
	return GetHMParticleIndex(pdgs, kFinalState);
	};

	inline bool HasFSNuElectron (void) const { return HasFSParticle(12); };
	inline bool HasFSNuMuon (void) const { return HasFSParticle(14); };
	inline bool HasFSNuTau (void) const { return HasFSParticle(16); };
	inline bool HasFSElectron (void) const { return HasFSParticle(11); };
	inline bool HasFSMuon (void) const { return HasFSParticle(13); };
	inline bool HasFSTau (void) const { return HasFSParticle(15); };
	inline bool HasFSProton (void) const { return HasFSParticle(2212); };
	inline bool HasFSNeutron (void) const { return HasFSParticle(2112); };
	inline bool HasFSPiZero (void) const { return HasFSParticle(111); };
	inline bool HasFSPiPlus (void) const { return HasFSParticle(211); };
	inline bool HasFSPiMinus (void) const { return HasFSParticle(-211); };
	inline bool HasFSPhoton (void) const { return HasFSParticle(22); };
	inline bool HasFSLeptons (void) const { return HasFSParticle(PhysConst::pdg_leptons); };
	inline bool HasFSPions (void) const { return HasFSParticle(PhysConst::pdg_pions); };
	inline bool HasFSChargePions (void) const { return HasFSParticle(PhysConst::pdg_charged_pions); };

	inline int NumFSNuElectron (void) const { return NumFSParticle(12); };
	inline int NumFSNuMuon (void) const { return NumFSParticle(14); };
	inline int NumFSNuTau (void) const { return NumFSParticle(16); };
	inline int NumFSElectron (void) const { return NumFSParticle(11); };
	inline int NumFSMuon (void) const { return NumFSParticle(13); };
	inline int NumFSTau (void) const { return NumFSParticle(15); };
	inline int NumFSProton (void) const { return NumFSParticle(2212); };
	inline int NumFSNeutron (void) const { return NumFSParticle(2112); };
	inline int NumFSPiZero (void) const { return NumFSParticle(111); };
	inline int NumFSPiPlus (void) const { return NumFSParticle(211); };
	inline int NumFSPiMinus (void) const { return NumFSParticle(-211); };
	inline int NumFSPhoton (void) const { return NumFSParticle(22); };
	int NumFSLeptons (void) const; // { return NumFSParticle(PhysConst::pdg_leptons); };
	inline int NumFSPions (void) const { return NumFSParticle(PhysConst::pdg_pions); };
	inline int NumFSChargePions (void) const { return NumFSParticle(PhysConst::pdg_charged_pions); };

	inline std::vector<int> GetAllFSNuElectronIndices (void) const { return GetAllFSParticleIndices(12); };
	inline std::vector<int> GetAllFSNuMuonIndices (void) const { return GetAllFSParticleIndices(14); };
	inline std::vector<int> GetAllFSNuTauIndices (void) const { return GetAllFSParticleIndices(16); };
	inline std::vector<int> GetAllFSElectronIndices (void) const { return GetAllFSParticleIndices(11); };
	inline std::vector<int> GetAllFSMuonIndices (void) const { return GetAllFSParticleIndices(13); };
	inline std::vector<int> GetAllFSTauIndices (void) const { return GetAllFSParticleIndices(15); };
	inline std::vector<int> GetAllFSProtonIndices (void) const { return GetAllFSParticleIndices(2212); };
	inline std::vector<int> GetAllFSNeutronIndices (void) const { return GetAllFSParticleIndices(2112); };
	inline std::vector<int> GetAllFSPiZeroIndices (void) const { return GetAllFSParticleIndices(111); };
	inline std::vector<int> GetAllFSPiPlusIndices (void) const { return GetAllFSParticleIndices(211); };
	inline std::vector<int> GetAllFSPiMinusIndices (void) const { return GetAllFSParticleIndices(-211); };
	inline std::vector<int> GetAllFSPhotonIndices (void) const { return GetAllFSParticleIndices(22); };
	inline std::vector<int> GetAllFSLeptonsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_leptons); };
	inline std::vector<int> GetAllFSPionsIndices (void) const { return GetAllFSParticleIndices(PhysConst::pdg_pions); };
	inline std::vector<int> GetAllFSChargePionsIndices(void) const { return GetAllFSParticleIndices(PhysConst::pdg_charged_pions); };

	inline std::vector<FitParticle*> GetAllFSNuElectron (void) { return GetAllFSParticle(12); };
	inline std::vector<FitParticle*> GetAllFSNuMuon (void) { return GetAllFSParticle(14); };
	inline std::vector<FitParticle*> GetAllFSNuTau (void) { return GetAllFSParticle(16); };
	inline std::vector<FitParticle*> GetAllFSElectron (void) { return GetAllFSParticle(11); };
	inline std::vector<FitParticle*> GetAllFSMuon (void) { return GetAllFSParticle(13); };
	inline std::vector<FitParticle*> GetAllFSTau (void) { return GetAllFSParticle(15); };
	inline std::vector<FitParticle*> GetAllFSProton (void) { return GetAllFSParticle(2212); };
	inline std::vector<FitParticle*> GetAllFSNeutron (void) { return GetAllFSParticle(2112); };
	inline std::vector<FitParticle*> GetAllFSPiZero (void) { return GetAllFSParticle(111); };
	inline std::vector<FitParticle*> GetAllFSPiPlus (void) { return GetAllFSParticle(211); };
	inline std::vector<FitParticle*> GetAllFSPiMinus (void) { return GetAllFSParticle(-211); };
	inline std::vector<FitParticle*> GetAllFSPhoton (void) { return GetAllFSParticle(22); };
	inline std::vector<FitParticle*> GetAllFSLeptons (void) { return GetAllFSParticle(PhysConst::pdg_leptons); };
	inline std::vector<FitParticle*> GetAllFSPions (void) { return GetAllFSParticle(PhysConst::pdg_pions); };
	inline std::vector<FitParticle*> GetAllFSChargePions (void) { return GetAllFSParticle(PhysConst::pdg_charged_pions); };

	inline FitParticle* GetHMFSNuElectron (void) { return GetHMFSParticle(12); };
	inline FitParticle* GetHMFSNuMuon (void) { return GetHMFSParticle(14); };
	inline FitParticle* GetHMFSNuTau (void) { return GetHMFSParticle(16); };
	inline FitParticle* GetHMFSElectron (void) { return GetHMFSParticle(11); };
	inline FitParticle* GetHMFSMuon (void) { return GetHMFSParticle(13); };
	inline FitParticle* GetHMFSTau (void) { return GetHMFSParticle(15); };
	inline FitParticle* GetHMFSAnyLeptons (void) { return GetHMFSParticle(PhysConst::pdg_all_leptons); };
	inline FitParticle* GetHMFSProton (void) { return GetHMFSParticle(2212); };
	inline FitParticle* GetHMFSNeutron (void) { return GetHMFSParticle(2112); };
	inline FitParticle* GetHMFSPiZero (void) { return GetHMFSParticle(111); };
	inline FitParticle* GetHMFSPiPlus (void) { return GetHMFSParticle(211); };
	inline FitParticle* GetHMFSPiMinus (void) { return GetHMFSParticle(-211); };
	inline FitParticle* GetHMFSPhoton (void) { return GetHMFSParticle(22); };

	inline FitParticle* GetHMFSLeptons (void) { return GetHMFSParticle(PhysConst::pdg_leptons); };
	inline FitParticle* GetHMFSAnyLepton (void) { return GetHMFSParticle(PhysConst::pdg_all_leptons); };

	inline FitParticle* GetHMFSPions (void) { return GetHMFSParticle(PhysConst::pdg_pions); };
	inline FitParticle* GetHMFSChargePions(void) { return GetHMFSParticle(PhysConst::pdg_charged_pions); };

	inline int GetHMFSNuElectronIndex (void) const { return GetHMFSParticleIndex(12); };
	inline int GetHMFSNuMuonIndex (void) const { return GetHMFSParticleIndex(14); };
	inline int GetHMFSNuTauIndex (void) const { return GetHMFSParticleIndex(16); };
	inline int GetHMFSElectronIndex (void) const { return GetHMFSParticleIndex(11); };
	inline int GetHMFSMuonIndex (void) const { return GetHMFSParticleIndex(13); };
	inline int GetHMFSTauIndex (void) const { return GetHMFSParticleIndex(15); };
	inline int GetHMFSProtonIndex (void) const { return GetHMFSParticleIndex(2212); };
	inline int GetHMFSNeutronIndex (void) const { return GetHMFSParticleIndex(2112); };
	inline int GetHMFSPiZeroIndex (void) const { return GetHMFSParticleIndex(111); };
	inline int GetHMFSPiPlusIndex (void) const { return GetHMFSParticleIndex(211); };
	inline int GetHMFSPiMinusIndex (void) const { return GetHMFSParticleIndex(-211); };
	inline int GetHMFSPhotonIndex (void) const { return GetHMFSParticleIndex(22); };

	inline int GetHMFSLeptonsIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_leptons); };
	inline int GetHMFSAnyLeptonIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_all_leptons); };

	inline int GetHMFSPionsIndex (void) const { return GetHMFSParticleIndex(PhysConst::pdg_pions); };
	inline int GetHMFSChargePionsIndex(void) const { return GetHMFSParticleIndex(PhysConst::pdg_charged_pions); };

	// ---- NEUTRINO INCOMING Related Functions
	int GetBeamNeutrinoIndex (void) const;
	inline TLorentzVector GetBeamNeutrinoP4 (void) const { return GetParticleP4(GetBeamNeutrinoIndex()); };
	inline TVector3 GetBeamNeutrinoP3 (void) const { return GetParticleP3(GetBeamNeutrinoIndex()); };
	inline double GetBeamNeutrinoMom (void) const { return GetParticleMom(GetBeamNeutrinoIndex()); };
	inline double GetBeamNeutrinoMom2 (void) const { return GetParticleMom2(GetBeamNeutrinoIndex()); };
	inline double GetBeamNeutrinoE (void) const { return GetParticleE(GetBeamNeutrinoIndex()); };
	inline double Enu (void) const { return GetBeamNeutrinoE(); };
	inline int GetBeamNeutrinoPDG (void) const { return GetParticlePDG(GetBeamNeutrinoIndex()); };
	inline int PDGnu (void) const { return GetBeamNeutrinoPDG(); };
	inline int GetNeutrinoInPos (void) const { return GetBeamNeutrinoIndex(); };
	inline FitParticle* GetBeamNeutrino (void) { return GetParticle(GetBeamNeutrinoIndex()); };
	inline FitParticle* GetNeutrinoIn (void) { return GetParticle(GetBeamNeutrinoIndex()); };


	// ---- Electron INCOMING Related Functions
	int GetBeamElectronIndex (void) const;
	inline TLorentzVector GetBeamElectronP4 (void) const { return GetParticleP4(GetBeamElectronIndex()); };
	inline TVector3 GetBeamElectronP3 (void) const { return GetParticleP3(GetBeamElectronIndex()); };
	inline double GetBeamElectronMom (void) const { return GetParticleMom(GetBeamElectronIndex()); };
	inline double GetBeamElectronMom2 (void) const { return GetParticleMom2(GetBeamElectronIndex()); };
	inline double GetBeamElectronE (void) const { return GetParticleE(GetBeamElectronIndex()); };
	inline FitParticle* GetBeamElectron (void) { return GetParticle(GetBeamElectronIndex()); };

	// ---- Pion INCOMING Functions
	int GetBeamPionIndex (void) const;
	inline TLorentzVector GetBeamPionP4 (void) const { return GetParticleP4(GetBeamPionIndex()); };
	inline TVector3 GetBeamPionP3 (void) const { return GetParticleP3(GetBeamPionIndex()); };
	inline double GetBeamPionMom (void) const { return GetParticleMom(GetBeamPionIndex()); };
	inline double GetBeamPionMom2 (void) const { return GetParticleMom2(GetBeamPionIndex()); };
	inline double GetBeamPionE (void) const { return GetParticleE(GetBeamPionIndex()); };
	inline FitParticle* GetBeamPion (void) { return GetParticle(GetBeamPionIndex()); };

	/// Legacy Functions
	inline FitParticle* PartInfo(uint i) { return GetParticle(i); };
	inline UInt_t Npart (void) const { return NPart(); };
	inline UInt_t NPart (void) const { return fNParticles; };



	// Other Functions
	int NumFSMesons();
	inline int GetBeamPartPos(void) const { return 0; };
	FitParticle* GetBeamPart(void);

	int GetLeptonOutPos(void) const;
	FitParticle* GetLeptonOut(void);


	// Event Information
	UInt_t fEventNo;
	double fTotCrs;
	int fTargetA;
	int fTargetZ;
	int fTargetH;
	bool fBound;
	int fDistance;
	int fTargetPDG;

	// Reduced Particle Stack
	UInt_t kMaxParticles;
	int fNParticles;
	double** fParticleMom;
	UInt_t* fParticleState;
	int* fParticlePDG;
	FitParticle** fParticleList;

	double** fOrigParticleMom;
	UInt_t* fOrigParticleState;
	int* fOrigParticlePDG;

	double* fNEUT_ParticleStatusCode;
	double* fNEUT_ParticleAliveCode;
	GeneratorInfoBase* fGenInfo;

	// Config Options for this class
	bool kRemoveFSIParticles;
	bool kRemoveUndefParticles;



	};
	/! @} /
	#endif
	diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
	index 600d076..0b1ac2a 100644
	--- a/src/InputHandler/GENIEInputHandler.cxx
	+++ b/src/InputHandler/GENIEInputHandler.cxx
	@@ -1,523 +1,526 @@
	// 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();
	- fNUISANCEEvent->fTargetA = 0.0;
	- fNUISANCEEvent->fTargetZ = 0.0;
	+ // Set the TargetPDG
	+ fNUISANCEEvent->fTargetPDG = fGenieGHep->TargetNucleus()->Pdg();
	+ // Set the A and Z and H from the target PDG
	+ fNUISANCEEvent->fTargetA = TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
	+ fNUISANCEEvent->fTargetZ = TargetUtils::GetTargetZFromPDG(fNUISANCEEvent->fTargetPDG);
	fNUISANCEEvent->fTargetH = 0;
	- fNUISANCEEvent->fBound = 0.0;
	+ fNUISANCEEvent->fBound = (fNUISANCEEvent->fTargetA != 1);
	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();

	// Add to N particle count
	fNUISANCEEvent->fNParticles++;

	// Extra Check incase GENIE fails.
	if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
	ERR(WRN) << "Number of GENIE Particles exceeds maximum!" << std::endl;
	ERR(WRN) << "Extend kMax, or run without including FSI particles!"
	<< std::endl;
	break;
	}
	}

	// Fill Extra Stack
	if (fSaveExtra) fGenieInfo->FillGeneratorInfo(fGenieNtpl);

	// Run Initial, FSI, Final, Other ordering.
	fNUISANCEEvent->OrderStack();

	FitParticle* ISNeutralLepton =
	fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
	if (ISNeutralLepton) {
	fNUISANCEEvent->probe_E = ISNeutralLepton->E();
	fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
	}

	return;
	}

	void GENIEInputHandler::Print() {}

	#endif
	diff --git a/src/InputHandler/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx
	index 384a326..9c7bf05 100644
	--- a/src/InputHandler/NuWroInputHandler.cxx
	+++ b/src/InputHandler/NuWroInputHandler.cxx
	@@ -1,479 +1,480 @@
	#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;
	+ 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);
	}

	// Sort Particles
	evt->fNParticles = 0;
	std::vector<particle>::iterator p_iter;

	// Initial State
	for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end();
	p_iter++) {
	AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState);
	}

	// FSI State
	// for (size_t i = 0; i < npart_in; i++ ) {
	// AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState);
	// }

	// Final State
	for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end();
	p_iter++) {
	AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState);
	}

	// 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();
	}

	return;
	}

	void NuWroInputHandler::AddNuWroParticle(FitEvent* evt, particle& p,
	int state) {
	// 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;

	// Status/PDG
	evt->fParticleState[evt->fNParticles] = state;
	evt->fParticlePDG[evt->fNParticles] = p.pdg;

	// Add to particle count
	evt->fNParticles++;
	}

	void NuWroInputHandler::Print() {}

	#endif

	diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
	index 0543dc5..6ab61b7 100644
	--- a/src/MCStudies/GenericFlux_Vectors.cxx
	+++ b/src/MCStudies/GenericFlux_Vectors.cxx
	@@ -1,363 +1,368 @@
	// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret

	/*******************************************************************************
	* This file is part of NUISANCE.
	*
	* NUISANCE is free software: you can redistribute it and/or modify
	* it under the terms of the GNU General Public License as published by
	* the Free Software Foundation, either version 3 of the License, or
	* (at your option) any later version.
	*
	* NUISANCE is distributed in the hope that it will be useful,
	* but WITHOUT ANY WARRANTY; without even the implied warranty of
	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	* GNU General Public License for more details.
	*
	* You should have received a copy of the GNU General Public License
	* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
	*******************************************************************************/

	#include "GenericFlux_Vectors.h"

	GenericFlux_Vectors::GenericFlux_Vectors(std::string name,
	std::string inputfile, FitWeight *rw,
	std::string type,
	std::string fakeDataFile) {
	// Measurement Details
	fName = name;
	eventVariables = NULL;

	// Define our energy range for flux calcs
	EnuMin = 0.;
	EnuMax = 1E10; // Arbritrarily high energy limit

	if (Config::HasPar("EnuMin")) {
	EnuMin = Config::GetParD("EnuMin");
	}

	if (Config::HasPar("EnuMax")) {
	EnuMax = Config::GetParD("EnuMax");
	}

	// Set default fitter flags
	fIsDiag = true;
	fIsShape = false;
	fIsRawEvents = false;

	// This function will sort out the input files automatically and parse all the
	// inputs,flags,etc.
	// There may be complex cases where you have to do this by hand, but usually
	// this will do.
	Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);

	eventVariables = NULL;

	// Setup fDataHist as a placeholder
	this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
	this->SetupDefaultHist();
	fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
	covar = StatUtils::GetInvert(fFullCovar);

	// 1. The generator is organised in SetupMeasurement so it gives the
	// cross-section in "per nucleon" units.
	// So some extra scaling for a specific measurement may be required. For
	// Example to get a "per neutron" measurement on carbon
	// which we do here, we have to multiple by the number of nucleons 12 and
	// divide by the number of neutrons 6.
	// N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is
	// often included here in other classes that directly integrate the event
	// histogram. This method is used here as it now respects EnuMin and EnuMax
	// correctly.
	this->fScaleFactor =
	(this->PredictedEventRate("width", 0, EnuMax) / double(fNEvents)) /
	this->TotalIntegratedFlux("width");

	LOG(SAM) << " Generic Flux Scaling Factor = " << fScaleFactor
	<< " [= " << (GetEventHistogram()->Integral("width") * 1E-38) << "/("
	<< (fNEvents + 0.) << "*" << TotalIntegratedFlux("width") << ")]"
	<< std::endl;

	if (fScaleFactor <= 0.0) {
	ERR(WRN) << "SCALE FACTOR TOO LOW " << std::endl;
	throw;
	}

	// Setup our TTrees
	this->AddEventVariablesToTree();
	this->AddSignalFlagsToTree();
	}

	void GenericFlux_Vectors::AddEventVariablesToTree() {
	// Setup the TTree to save everything
	if (!eventVariables) {
	Config::Get().out->cd();
	eventVariables = new TTree((this->fName + "_VARS").c_str(),
	(this->fName + "_VARS").c_str());
	}

	LOG(SAM) << "Adding Event Variables" << std::endl;

	eventVariables->Branch("Mode", &Mode, "Mode/I");
	eventVariables->Branch("cc", &cc, "cc/B");
	eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I");
	eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F");
	eventVariables->Branch("tgt", &tgt, "tgt/I");
	+ eventVariables->Branch("tgta", &tgta, "tgta/I");
	+ eventVariables->Branch("tgtz", &tgtz, "tgtz/I");
	eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I");
	eventVariables->Branch("ELep", &ELep, "ELep/F");
	eventVariables->Branch("CosLep", &CosLep, "CosLep/F");

	// Basic interaction kinematics
	eventVariables->Branch("Q2", &Q2, "Q2/F");
	eventVariables->Branch("q0", &q0, "q0/F");
	eventVariables->Branch("q3", &q3, "q3/F");
	eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F");
	eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F");
	eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F");
	eventVariables->Branch("W", &W, "W/F");
	eventVariables->Branch("W_genie", &W_genie, "W_genie/F");
	eventVariables->Branch("x", &x, "x/F");
	eventVariables->Branch("y", &y, "y/F");
	eventVariables->Branch("Eav", &Eav, "Eav/F");
	eventVariables->Branch("EavAlt", &EavAlt, "EavAlt/F");

	eventVariables->Branch("dalphat", &dalphat, "dalphat/F");
	eventVariables->Branch("dpt", &dpt, "dpt/F");
	eventVariables->Branch("dphit", &dphit, "dphit/F");
	eventVariables->Branch("pnreco_C", &pnreco_C, "pnreco_C/F");

	// Save outgoing particle vectors
	eventVariables->Branch("nfsp", &nfsp, "nfsp/I");
	eventVariables->Branch("px", px, "px[nfsp]/F");
	eventVariables->Branch("py", py, "py[nfsp]/F");
	eventVariables->Branch("pz", pz, "pz[nfsp]/F");
	eventVariables->Branch("E", E, "E[nfsp]/F");
	eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I");

	// Event Scaling Information
	eventVariables->Branch("Weight", &Weight, "Weight/F");
	eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
	eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
	// Should be a double because may be 1E-39 and less
	eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D");

	// The customs
	eventVariables->Branch("CustomWeight", &CustomWeight, "CustomWeight/F");
	eventVariables->Branch("CustomWeightArray", CustomWeightArray, "CustomWeightArray[6]/F");

	return;
	}

	void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {

	ResetVariables();

	// Fill Signal Variables
	FillSignalFlags(event);
	LOG(DEB) << "Filling signal" << std::endl;

	// Now fill the information
	Mode = event->Mode;
	cc = (abs(event->Mode) < 30);

	// Get the incoming neutrino and outgoing lepton
	FitParticle *nu = event->GetNeutrinoIn();
	FitParticle *lep = event->GetHMFSAnyLepton();

	PDGnu = nu->fPID;
	Enu_true = nu->fP.E() / 1E3;
	tgt = event->fTargetPDG;
	+ tgta = event->fTargetA;
	+ tgtz = event->fTargetZ;
	+
	if (lep != NULL) {
	PDGLep = lep->fPID;
	ELep = lep->fP.E() / 1E3;
	CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));

	// Basic interaction kinematics
	Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6;
	q0 = (nu->fP - lep->fP).E() / 1E3;
	q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3;

	// These assume C12 binding from MINERvA... not ideal
	Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
	Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);

	Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event)/1.E3;
	EavAlt = FitUtils::Eavailable(event)/1.E3;

	// Get W_true with assumption of initial state nucleon at rest
	float m_n = (float)PhysConst::mass_proton;
	// Q2 assuming nucleon at rest
	W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
	// True Q2
	W = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
	x = Q2 / (2 * m_n * q0);
	y = 1 - ELep / Enu_true;

	dalphat = FitUtils::Get_STV_dalphat(event, PDGnu, true);
	dpt = FitUtils::Get_STV_dpt(event, PDGnu, true);
	dphit = FitUtils::Get_STV_dphit(event, PDGnu, true);
	pnreco_C = FitUtils::Get_pn_reco_C(event, PDGnu, true);
	}

	// Loop over the particles and store all the final state particles in a vector
	for (UInt_t i = 0; i < event->Npart(); ++i) {

	bool part_alive = event->PartInfo(i)->fIsAlive &&
	event->PartInfo(i)->Status() == kFinalState;
	if (!part_alive) continue;

	partList.push_back(event->PartInfo(i));
	}

	// Save outgoing particle vectors
	nfsp = (int)partList.size();

	for (int i = 0; i < nfsp; ++i) {
	px[i] = partList[i]->fP.X() / 1E3;
	py[i] = partList[i]->fP.Y() / 1E3;
	pz[i] = partList[i]->fP.Z() / 1E3;
	E[i] = partList[i]->fP.E() / 1E3;
	pdg[i] = partList[i]->fPID;
	}

	#ifdef __GENIE_ENABLED__
	if (event->fType == kGENIE) {
	EventRecord * gevent = static_cast<EventRecord*>(event->genie_event->event);
	const Interaction * interaction = gevent->Summary();
	const Kinematics & kine = interaction->Kine();
	double W_genie = kine.W();
	}
	#endif

	// Fill event weights
	Weight = event->RWWeight * event->InputWeight;
	RWWeight = event->RWWeight;
	InputWeight = event->InputWeight;
	// And the Customs
	CustomWeight = event->CustomWeight;
	for (int i = 0; i < 6; ++i) {
	CustomWeightArray[i] = event->CustomWeightArray[i];
	}

	// Fill the eventVariables Tree
	eventVariables->Fill();
	return;
	};

	//********************************************************************
	void GenericFlux_Vectors::ResetVariables() {
	//********************************************************************

	cc = false;

	// Reset all Function used to extract any variables of interest to the event
	- Mode = PDGnu = tgt = PDGLep = 0;
	+ Mode = PDGnu = tgt = tgta = tgtz = PDGLep = 0;

	Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W = x = y = Eav = EavAlt = -999.9;

	W_genie = -999;
	// Other fun variables
	// MINERvA-like ones
	dalphat = dpt = dphit = pnreco_C = -999.99;

	nfsp = 0;
	for (int i = 0; i < kMAX; ++i){
	px[i] = py[i] = pz[i] = E[i] = -999;
	pdg[i] = 0;
	}

	Weight = InputWeight = RWWeight = 0.0;

	CustomWeight = 0.0;
	for (int i = 0; i < 6; ++i) CustomWeightArray[i] = 0.0;

	partList.clear();

	flagCCINC = flagNCINC = flagCCQE = flagCC0pi = flagCCQELike = flagNCEL = flagNC0pi = flagCCcoh = flagNCcoh = flagCC1pip = flagNC1pip = flagCC1pim = flagNC1pim = flagCC1pi0 = flagNC1pi0 = false;
	}

	//********************************************************************
	void GenericFlux_Vectors::FillSignalFlags(FitEvent *event) {
	//********************************************************************

	// Some example flags are given from SignalDef.
	// See src/Utils/SignalDef.cxx for more.
	int nuPDG = event->PartInfo(0)->fPID;

	// Generic signal flags
	flagCCINC = SignalDef::isCCINC(event, nuPDG);
	flagNCINC = SignalDef::isNCINC(event, nuPDG);
	flagCCQE = SignalDef::isCCQE(event, nuPDG);
	flagCCQELike = SignalDef::isCCQELike(event, nuPDG);
	flagCC0pi = SignalDef::isCC0pi(event, nuPDG);
	flagNCEL = SignalDef::isNCEL(event, nuPDG);
	flagNC0pi = SignalDef::isNC0pi(event, nuPDG);
	flagCCcoh = SignalDef::isCCCOH(event, nuPDG, 211);
	flagNCcoh = SignalDef::isNCCOH(event, nuPDG, 111);
	flagCC1pip = SignalDef::isCC1pi(event, nuPDG, 211);
	flagNC1pip = SignalDef::isNC1pi(event, nuPDG, 211);
	flagCC1pim = SignalDef::isCC1pi(event, nuPDG, -211);
	flagNC1pim = SignalDef::isNC1pi(event, nuPDG, -211);
	flagCC1pi0 = SignalDef::isCC1pi(event, nuPDG, 111);
	flagNC1pi0 = SignalDef::isNC1pi(event, nuPDG, 111);
	}

	void GenericFlux_Vectors::AddSignalFlagsToTree() {
	if (!eventVariables) {
	Config::Get().out->cd();
	eventVariables = new TTree((this->fName + "_VARS").c_str(),
	(this->fName + "_VARS").c_str());
	}

	LOG(SAM) << "Adding signal flags" << std::endl;

	// Signal Definitions from SignalDef.cxx
	eventVariables->Branch("flagCCINC", &flagCCINC, "flagCCINC/O");
	eventVariables->Branch("flagNCINC", &flagNCINC, "flagNCINC/O");
	eventVariables->Branch("flagCCQE", &flagCCQE, "flagCCQE/O");
	eventVariables->Branch("flagCC0pi", &flagCC0pi, "flagCC0pi/O");
	eventVariables->Branch("flagCCQELike", &flagCCQELike, "flagCCQELike/O");
	eventVariables->Branch("flagNCEL", &flagNCEL, "flagNCEL/O");
	eventVariables->Branch("flagNC0pi", &flagNC0pi, "flagNC0pi/O");
	eventVariables->Branch("flagCCcoh", &flagCCcoh, "flagCCcoh/O");
	eventVariables->Branch("flagNCcoh", &flagNCcoh, "flagNCcoh/O");
	eventVariables->Branch("flagCC1pip", &flagCC1pip, "flagCC1pip/O");
	eventVariables->Branch("flagNC1pip", &flagNC1pip, "flagNC1pip/O");
	eventVariables->Branch("flagCC1pim", &flagCC1pim, "flagCC1pim/O");
	eventVariables->Branch("flagNC1pim", &flagNC1pim, "flagNC1pim/O");
	eventVariables->Branch("flagCC1pi0", &flagCC1pi0, "flagCC1pi0/O");
	eventVariables->Branch("flagNC1pi0", &flagNC1pi0, "flagNC1pi0/O");
	};

	void GenericFlux_Vectors::Write(std::string drawOpt) {

	// First save the TTree
	eventVariables->Write();

	// Save Flux and Event Histograms too
	GetInput()->GetFluxHistogram()->Write();
	GetInput()->GetEventHistogram()->Write();

	return;
	}

	// Override functions which aren't really necessary
	bool GenericFlux_Vectors::isSignal(FitEvent *event) {
	(void)event;
	return true;
	};

	void GenericFlux_Vectors::ScaleEvents() { return; }

	void GenericFlux_Vectors::ApplyNormScale(float norm) {
	this->fCurrentNorm = norm;
	return;
	}

	void GenericFlux_Vectors::FillHistograms() { return; }

	void GenericFlux_Vectors::ResetAll() {
	eventVariables->Reset();
	return;
	}

	float GenericFlux_Vectors::GetChi2() { return 0.0; }
	diff --git a/src/MCStudies/GenericFlux_Vectors.h b/src/MCStudies/GenericFlux_Vectors.h
	index 101335b..181a283 100644
	--- a/src/MCStudies/GenericFlux_Vectors.h
	+++ b/src/MCStudies/GenericFlux_Vectors.h
	@@ -1,133 +1,135 @@
	// 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/>.
	*******************************************************************************/

	#ifndef GenericFlux_Vectors_H_SEEN
	#define GenericFlux_Vectors_H_SEEN
	#include "Measurement1D.h"
	#include "FitEvent.h"

	class GenericFlux_Vectors : public Measurement1D {

	public:

	GenericFlux_Vectors(std::string name, std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile);
	virtual ~GenericFlux_Vectors() {};

	//! Grab info from event
	void FillEventVariables(FitEvent *event);

	//! Fill signal flags
	void FillSignalFlags(FitEvent *event);

	void ResetVariables();

	//! Fill Custom Histograms
	void FillHistograms();

	//! ResetAll
	void ResetAll();

	//! Scale
	void ScaleEvents();

	//! Norm
	void ApplyNormScale(float norm);

	//! Define this samples signal
	bool isSignal(FitEvent *nvect);

	//! Write Files
	void Write(std::string drawOpt);

	//! Get Chi2
	float GetChi2();

	void AddEventVariablesToTree();
	void AddSignalFlagsToTree();

	private:

	TTree* eventVariables;
	std::vector<FitParticle*> partList;

	int Mode;
	bool cc;
	int PDGnu;
	int tgt;
	+ int tgta;
	+ int tgtz;
	int PDGLep;
	float ELep;
	float CosLep;

	// Basic interaction kinematics
	float Q2;
	float q0;
	float q3;
	float Enu_QE;
	float Enu_true;
	float Q2_QE;
	float W_nuc_rest;
	float W;
	float x;
	float y;
	float Eav;
	float EavAlt;
	float dalphat;
	float W_genie;
	float dpt;
	float dphit;
	float pnreco_C;

	// Save outgoing particle vectors
	int nfsp;
	static const int kMAX = 200;
	float px[kMAX];
	float py[kMAX];
	float pz[kMAX];
	float E[kMAX];
	int pdg[kMAX];

	// Basic event info
	float Weight;
	float InputWeight;
	float RWWeight;
	double fScaleFactor;

	// Custom weights
	float CustomWeight;
	float CustomWeightArray[6];

	// Generic signal flags
	bool flagCCINC;
	bool flagNCINC;
	bool flagCCQE;
	bool flagCC0pi;
	bool flagCCQELike;
	bool flagNCEL;
	bool flagNC0pi;
	bool flagCCcoh;
	bool flagNCcoh;
	bool flagCC1pip;
	bool flagNC1pip;
	bool flagCC1pim;
	bool flagNC1pim;
	bool flagCC1pi0;
	bool flagNC1pi0;
	};

	#endif

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 12:01 PM (56 m, 3 s)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5063563
Default Alt Text: (82 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions