No OneTemporary
Actions

Size

43 KB

Subscribers

None

View Options

	diff --git a/src/Reweight/CMakeLists.txt b/src/Reweight/CMakeLists.txt
	index 55b9fef..01296a8 100644
	--- a/src/Reweight/CMakeLists.txt
	+++ b/src/Reweight/CMakeLists.txt
	@@ -1,83 +1,84 @@
	# 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/>.
	################################################################################
	set(IMPLFILES
	GlobalDialList.cxx
	FitWeight.cxx
	WeightEngineBase.cxx
	NEUTWeightEngine.cxx
	NuWroWeightEngine.cxx
	GENIEWeightEngine.cxx
	WeightUtils.cxx
	SampleNormEngine.cxx
	LikelihoodWeightEngine.cxx
	SplineWeightEngine.cxx
	NUISANCESyst.cxx
	T2KWeightEngine.cxx
	NUISANCEWeightEngine.cxx
	NUISANCEWeightCalcs.cxx
	NIWGWeightEngine.cxx
	OscWeightEngine.cxx
	MINERvAWeightCalcs.cxx
	weightRPA.h
	)

	set(HEADERFILES
	GlobalDialList.h
	FitWeight.h
	WeightEngineBase.h
	NEUTWeightEngine.h
	NuWroWeightEngine.h
	GENIEWeightEngine.h
	WeightUtils.h
	SampleNormEngine.h
	LikelihoodWeightEngine.h
	SplineWeightEngine.h
	NUISANCESyst.h
	T2KWeightEngine.h
	NUISANCEWeightEngine.h
	NUISANCEWeightCalcs.h
	NIWGWeightEngine.h
	OscWeightEngine.h
	MINERvAWeightCalcs.h
	weightRPA.h
	+BeRPA.h
	)

	set(LIBNAME Reweight)

	if(CMAKE_BUILD_TYPE MATCHES DEBUG)
	add_library(${LIBNAME} STATIC ${IMPLFILES})
	else(CMAKE_BUILD_TYPE MATCHES RELEASE)
	add_library(${LIBNAME} SHARED ${IMPLFILES})
	endif()

	include_directories(${MINIMUM_INCLUDE_DIRECTORIES})

	set_target_properties(${LIBNAME} PROPERTIES VERSION
	"${ExtFit_VERSION_MAJOR}.${ExtFit_VERSION_MINOR}.${ExtFit_VERSION_REVISION}")
	#set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS})

	if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES)
	add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES})
	endif()

	install(TARGETS ${LIBNAME} DESTINATION lib)
	#Can uncomment this to install the headers... but is it really neccessary?
	install(FILES ${HEADERFILES} DESTINATION include)

	set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE)
	diff --git a/src/Reweight/NUISANCESyst.cxx b/src/Reweight/NUISANCESyst.cxx
	index 04c425c..d11114b 100644
	--- a/src/Reweight/NUISANCESyst.cxx
	+++ b/src/Reweight/NUISANCESyst.cxx
	@@ -1,74 +1,80 @@
	#include "NUISANCESyst.h"

	int Reweight::ConvertNUISANCEDial(std::string type) {

	for (int i = kUnkownNUISANCEDial + 1; i < kNUISANCEDial_LAST; i++) {
	if (!type.compare(ConvNUISANCEDial(i).c_str())) {
	return i;
	}
	}
	return kUnkownNUISANCEDial;
	};

	std::string Reweight::ConvNUISANCEDial(int type) {

	switch (type) {
	case kGaussianCorr_CCQE_norm: { return "GaussianCorr_CCQE_norm"; }
	case kGaussianCorr_CCQE_tilt: { return "GaussianCorr_CCQE_tilt"; }
	case kGaussianCorr_CCQE_Pq0: { return "GaussianCorr_CCQE_Pq0"; }
	case kGaussianCorr_CCQE_Wq0: { return "GaussianCorr_CCQE_Wq0"; }
	case kGaussianCorr_CCQE_Pq3: { return "GaussianCorr_CCQE_Pq3"; }
	case kGaussianCorr_CCQE_Wq3: { return "GaussianCorr_CCQE_Wq3"; }

	case kGaussianCorr_2p2h_norm: { return "GaussianCorr_2p2h_norm"; }
	case kGaussianCorr_2p2h_tilt: { return "GaussianCorr_2p2h_tilt"; }
	case kGaussianCorr_2p2h_Pq0: { return "GaussianCorr_2p2h_Pq0"; }
	case kGaussianCorr_2p2h_Wq0: { return "GaussianCorr_2p2h_Wq0"; }
	case kGaussianCorr_2p2h_Pq3: { return "GaussianCorr_2p2h_Pq3"; }
	case kGaussianCorr_2p2h_Wq3: { return "GaussianCorr_2p2h_Wq3"; }

	case kGaussianCorr_2p2h_PPandNN_norm: { return "GaussianCorr_2p2h_PPandNN_norm"; }
	case kGaussianCorr_2p2h_PPandNN_tilt: { return "GaussianCorr_2p2h_PPandNN_tilt"; }
	case kGaussianCorr_2p2h_PPandNN_Pq0: { return "GaussianCorr_2p2h_PPandNN_Pq0"; }
	case kGaussianCorr_2p2h_PPandNN_Wq0: { return "GaussianCorr_2p2h_PPandNN_Wq0"; }
	case kGaussianCorr_2p2h_PPandNN_Pq3: { return "GaussianCorr_2p2h_PPandNN_Pq3"; }
	case kGaussianCorr_2p2h_PPandNN_Wq3: { return "GaussianCorr_2p2h_PPandNN_Wq3"; }

	case kGaussianCorr_2p2h_NP_norm: { return "GaussianCorr_2p2h_NP_norm"; }
	case kGaussianCorr_2p2h_NP_tilt: { return "GaussianCorr_2p2h_NP_tilt"; }
	case kGaussianCorr_2p2h_NP_Pq0: { return "GaussianCorr_2p2h_NP_Pq0"; }
	case kGaussianCorr_2p2h_NP_Wq0: { return "GaussianCorr_2p2h_NP_Wq0"; }
	case kGaussianCorr_2p2h_NP_Pq3: { return "GaussianCorr_2p2h_NP_Pq3"; }
	case kGaussianCorr_2p2h_NP_Wq3: { return "GaussianCorr_2p2h_NP_Wq3"; }

	case kGaussianCorr_CC1pi_norm: { return "GaussianCorr_CC1pi_norm"; }
	case kGaussianCorr_CC1pi_tilt: { return "GaussianCorr_CC1pi_tilt"; }
	case kGaussianCorr_CC1pi_Pq0: { return "GaussianCorr_CC1pi_Pq0"; }
	case kGaussianCorr_CC1pi_Wq0: { return "GaussianCorr_CC1pi_Wq0"; }
	case kGaussianCorr_CC1pi_Pq3: { return "GaussianCorr_CC1pi_Pq3"; }
	case kGaussianCorr_CC1pi_Wq3: { return "GaussianCorr_CC1pi_Wq3"; }
	case kGaussianCorr_AllowSuppression: { return "GaussianCorr_AllowSuppression"; }

	+ case kBeRPA_A: { return "BeRPA_A"; }
	+ case kBeRPA_B: { return "BeRPA_B"; }
	+ case kBeRPA_D: { return "BeRPA_D"; }
	+ case kBeRPA_E: { return "BeRPA_E"; }
	+ case kBeRPA_U: { return "BeRPA_U"; }
	+
	case kModeNorm_NormRES: { return "NormRES"; }

	case kMINERvARW_NormCCQE: { return "MINERvARW_NormCCQE"; }
	case kMINERvARW_NormCCMEC: { return "MINERvARW_NormCCMEC"; }
	case kMINERvARW_NormCCRES: { return "MINERvARW_NormCCRES"; }

	case kMINERvARW_RikRPA_ApplyRPA: { return "MINERvARW_RikRPA_ApplyRPA"; }
	case kMINERvARW_RikRPA_LowQ2: { return "MINERvARW_RikRPA_LowQ2"; }
	case kMINERvARW_RikRPA_HighQ2: { return "MINERvARW_RikRPA_HighQ2"; }

	case kMINERvARW_RikRESRPA_ApplyRPA: { return "MINERvARW_RikRESRPA_ApplyRPA"; }
	case kMINERvARW_RikRESRPA_LowQ2: { return "MINERvARW_RikRESRPA_LowQ2"; }
	case kMINERvARW_RikRESRPA_HighQ2: { return "MINERvARW_RikRESRPA_HighQ2"; }


	case kSBLOsc_Distance: { return "SBLOsc_Distance"; }
	case kSBLOsc_MassSplitting: { return "SBLOsc_MassSplitting"; }
	case kSBLOsc_Sin2Theta: { return "SBLOsc_Sin2Theta"; }

	default: return "unknown_nuisance_dial";
	}

	};
	diff --git a/src/Reweight/NUISANCESyst.h b/src/Reweight/NUISANCESyst.h
	index 652e2ff..14825ef 100644
	--- a/src/Reweight/NUISANCESyst.h
	+++ b/src/Reweight/NUISANCESyst.h
	@@ -1,73 +1,78 @@
	#ifndef NUISANCESyst_H
	#define NUISANCESyst_H
	#include "GeneralUtils.h"

	namespace Reweight {

	enum NUISANCESyst {

	kUnkownNUISANCEDial = 0,

	kGaussianCorr_CCQE_norm,
	kGaussianCorr_CCQE_tilt,
	kGaussianCorr_CCQE_Pq0,
	kGaussianCorr_CCQE_Wq0,
	kGaussianCorr_CCQE_Pq3,
	kGaussianCorr_CCQE_Wq3,

	kGaussianCorr_2p2h_norm,
	kGaussianCorr_2p2h_tilt,
	kGaussianCorr_2p2h_Pq0,
	kGaussianCorr_2p2h_Wq0,
	kGaussianCorr_2p2h_Pq3,
	kGaussianCorr_2p2h_Wq3,

	kGaussianCorr_2p2h_PPandNN_norm,
	kGaussianCorr_2p2h_PPandNN_tilt,
	kGaussianCorr_2p2h_PPandNN_Pq0,
	kGaussianCorr_2p2h_PPandNN_Wq0,
	kGaussianCorr_2p2h_PPandNN_Pq3,
	kGaussianCorr_2p2h_PPandNN_Wq3,

	kGaussianCorr_2p2h_NP_norm,
	kGaussianCorr_2p2h_NP_tilt,
	kGaussianCorr_2p2h_NP_Pq0,
	kGaussianCorr_2p2h_NP_Wq0,
	kGaussianCorr_2p2h_NP_Pq3,
	kGaussianCorr_2p2h_NP_Wq3,

	kGaussianCorr_CC1pi_norm,
	kGaussianCorr_CC1pi_tilt,
	kGaussianCorr_CC1pi_Pq0,
	kGaussianCorr_CC1pi_Wq0,
	kGaussianCorr_CC1pi_Pq3,
	kGaussianCorr_CC1pi_Wq3,
	kGaussianCorr_AllowSuppression,

	+ kBeRPA_A,
	+ kBeRPA_B,
	+ kBeRPA_D,
	+ kBeRPA_E,
	+ kBeRPA_U,
	+
	kModeNorm_NormRES,

	kMINERvARW_NormCCQE,
	kMINERvARW_NormCCMEC,
	kMINERvARW_NormCCRES,
	//kMINERvARW_NormCCNonRES1pi,

	kMINERvARW_RikRPA_ApplyRPA,
	kMINERvARW_RikRPA_LowQ2,
	kMINERvARW_RikRPA_HighQ2,

	kMINERvARW_RikRESRPA_ApplyRPA,
	kMINERvARW_RikRESRPA_LowQ2,
	kMINERvARW_RikRESRPA_HighQ2,

	kSBLOsc_Distance,
	kSBLOsc_MassSplitting,
	kSBLOsc_Sin2Theta,

	kNUISANCEDial_LAST
	};

	int ConvertNUISANCEDial(std::string type);
	std::string ConvNUISANCEDial(int type);
	-
	};
	#endif
	diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx
	index be3c9ae..8c61eef 100644
	--- a/src/Reweight/NUISANCEWeightCalcs.cxx
	+++ b/src/Reweight/NUISANCEWeightCalcs.cxx
	@@ -1,444 +1,503 @@
	#include "NUISANCEWeightCalcs.h"

	#include "FitEvent.h"
	#include "GeneralUtils.h"
	#include "NUISANCESyst.h"
	#include "WeightUtils.h"

	using namespace Reweight;

	ModeNormCalc::ModeNormCalc() { fNormRES = 1.0; }

	double ModeNormCalc::CalcWeight(BaseFitEvt* evt) {
	int mode = abs(evt->Mode);
	double w = 1.0;

	if (mode == 11 or mode == 12 or mode == 13) {
	w *= fNormRES;
	}

	return w;
	}

	void ModeNormCalc::SetDialValue(std::string name, double val) {
	SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
	}

	void ModeNormCalc::SetDialValue(int rwenum, double val) {
	int curenum = rwenum % 1000;

	// Check Handled
	if (!IsHandled(curenum)) return;
	if (curenum == kModeNorm_NormRES) fNormRES = val;
	}

	bool ModeNormCalc::IsHandled(int rwenum) {
	int curenum = rwenum % 1000;
	switch (curenum) {
	case kModeNorm_NormRES:
	return true;
	default:
	return false;
	}
	}

	+BeRPACalc::BeRPACalc() : fBeRPA_A(0.59), fBeRPA_B(1.05), fBeRPA_D(1.13), fBeRPA_E(0.88), fBeRPA_U(1.2), nParams(0) {
	+ // A = 0.59 +/- 20%
	+ // B = 1.05 +/- 20%
	+ // D = 1.13 +/- 15%
	+ // E = 0.88 +/- 40%
	+ // U = 1.2
	+}
	+
	+double BeRPACalc::CalcWeight(BaseFitEvt* evt) {
	+ FitEvent* fevt = static_cast<FitEvent*>(evt);
	+ int mode = abs(evt->Mode);
	+ double w = 1.0;
	+ if (nParams == 0) {
	+ return w;
	+ }
	+
	+ // Get Q2
	+ // Get final state lepton
	+ if (mode == 1) {
	+ double Q2 = -1.0(fevt->GetHMFSAnyLeptons()->P4()-fevt->GetNeutrinoIn()->P4())(fevt->GetHMFSAnyLeptons()->P4()-fevt->GetNeutrinoIn()->P4())/1.E6;
	+ // Only CCQE events
	+ w *= calcRPA(Q2, fBeRPA_A, fBeRPA_B, fBeRPA_D, fBeRPA_E, fBeRPA_U);
	+ }
	+
	+ return w;
	+}
	+
	+void BeRPACalc::SetDialValue(std::string name, double val) {
	+ SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
	+}
	+
	+void BeRPACalc::SetDialValue(int rwenum, double val) {
	+ int curenum = rwenum % 1000;
	+
	+ // Check Handled
	+ if (!IsHandled(curenum)) return;
	+ // Need 4 or 5 reconfigures
	+ if (curenum == kBeRPA_A) fBeRPA_A = val;
	+ else if (curenum == kBeRPA_B) fBeRPA_B = val;
	+ else if (curenum == kBeRPA_D) fBeRPA_D = val;
	+ else if (curenum == kBeRPA_E) fBeRPA_E = val;
	+ else if (curenum == kBeRPA_U) fBeRPA_U = val;
	+ nParams++;
	+}
	+
	+bool BeRPACalc::IsHandled(int rwenum) {
	+ int curenum = rwenum % 1000;
	+ switch (curenum) {
	+ case kBeRPA_A:
	+ case kBeRPA_B:
	+ case kBeRPA_D:
	+ case kBeRPA_E:
	+ case kBeRPA_U:
	+ return true;
	+ default:
	+ return false;
	+ }
	+}
	+
	SBLOscWeightCalc::SBLOscWeightCalc() {
	fDistance = 0.0;
	fMassSplitting = 0.0;
	fSin2Theta = 0.0;
	}

	double SBLOscWeightCalc::CalcWeight(BaseFitEvt* evt) {
	FitEvent* fevt = static_cast<FitEvent*>(evt);

	FitParticle* pnu = fevt->PartInfo(0);
	double E = pnu->fP.E() / 1.E3;

	// Extract energy
	return GetSBLOscWeight(E);
	}

	void SBLOscWeightCalc::SetDialValue(std::string name, double val) {
	SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
	}

	void SBLOscWeightCalc::SetDialValue(int rwenum, double val) {
	int curenum = rwenum % 1000;
	if (!IsHandled(curenum)) return;
	if (curenum == kSBLOsc_Distance) fDistance = val;
	if (curenum == kSBLOsc_MassSplitting) fMassSplitting = val;
	if (curenum == kSBLOsc_Sin2Theta) fSin2Theta = val;
	}

	bool SBLOscWeightCalc::IsHandled(int rwenum) {
	int curenum = rwenum % 1000;
	switch (curenum) {
	case kSBLOsc_Distance:
	return true;
	case kSBLOsc_MassSplitting:
	return true;
	case kSBLOsc_Sin2Theta:
	return true;
	default:
	return false;
	}
	}

	double SBLOscWeightCalc::GetSBLOscWeight(double E) {
	if (E <= 0.0) return 1.0 - 0.5 * fSin2Theta;
	return 1.0 - fSin2Theta * pow(sin(1.267 * fMassSplitting * fDistance / E), 2);
	}

	GaussianModeCorr::GaussianModeCorr() {

	// Apply the tilt-shift Gauss by Patrick
	// Alternatively set in config
	fMethod = true;

	// Init
	fApply_CCQE = false;
	fGausVal_CCQE[kPosNorm] = 0.0;
	fGausVal_CCQE[kPosTilt] = 0.0;
	fGausVal_CCQE[kPosPq0] = 1.0;
	fGausVal_CCQE[kPosWq0] = 1.0;
	fGausVal_CCQE[kPosPq3] = 1.0;
	fGausVal_CCQE[kPosWq3] = 1.0;

	fApply_2p2h = false;
	fGausVal_2p2h[kPosNorm] = 0.0;
	fGausVal_2p2h[kPosTilt] = 0.0;
	fGausVal_2p2h[kPosPq0] = 1.0;
	fGausVal_2p2h[kPosWq0] = 1.0;
	fGausVal_2p2h[kPosPq3] = 1.0;
	fGausVal_2p2h[kPosWq3] = 1.0;

	fApply_2p2h_PPandNN = false;
	fGausVal_2p2h_PPandNN[kPosNorm] = 0.0;
	fGausVal_2p2h_PPandNN[kPosTilt] = 0.0;
	fGausVal_2p2h_PPandNN[kPosPq0] = 1.0;
	fGausVal_2p2h_PPandNN[kPosWq0] = 1.0;
	fGausVal_2p2h_PPandNN[kPosPq3] = 1.0;
	fGausVal_2p2h_PPandNN[kPosWq3] = 1.0;

	fApply_2p2h_NP = false;
	fGausVal_2p2h_NP[kPosNorm] = 0.0;
	fGausVal_2p2h_NP[kPosTilt] = 0.0;
	fGausVal_2p2h_NP[kPosPq0] = 1.0;
	fGausVal_2p2h_NP[kPosWq0] = 1.0;
	fGausVal_2p2h_NP[kPosPq3] = 1.0;
	fGausVal_2p2h_NP[kPosWq3] = 1.0;

	fApply_CC1pi = false;
	fGausVal_CC1pi[kPosNorm] = 0.0;
	fGausVal_CC1pi[kPosTilt] = 0.0;
	fGausVal_CC1pi[kPosPq0] = 1.0;
	fGausVal_CC1pi[kPosWq0] = 1.0;
	fGausVal_CC1pi[kPosPq3] = 1.0;
	fGausVal_CC1pi[kPosWq3] = 1.0;

	fAllowSuppression = false;

	fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG");
	//fDebugStatements = true;
	}

	double GaussianModeCorr::CalcWeight(BaseFitEvt* evt) {
	FitEvent* fevt = static_cast<FitEvent*>(evt);
	double rw_weight = 1.0;

	// Get Neutrino
	if (!fevt->Npart()) {
	THROW("NO particles found in stack!");
	}
	FitParticle* pnu = fevt->GetHMISAnyLeptons();
	if (!pnu) {
	THROW("NO Starting particle found in stack!");
	}
	int pdgnu = pnu->fPID;

	int expect_fsleppdg = 0;

	if (pdgnu & 1) {
	expect_fsleppdg = pdgnu;
	} else {
	expect_fsleppdg = abs(pdgnu) - 1;
	}

	FitParticle* plep = fevt->GetHMFSParticle(expect_fsleppdg);
	if (!plep) return 1.0;

	TLorentzVector q = pnu->fP - plep->fP;

	// Extra q0,q3
	double q0 = fabs(q.E()) / 1.E3;
	double q3 = fabs(q.Vect().Mag()) / 1.E3;

	int initialstate = -1; // Undef
	if (abs(fevt->Mode) == 2) {
	int npr = 0;
	int nne = 0;

	for (UInt_t j = 0; j < fevt->Npart(); j++) {
	if ((fevt->PartInfo(j))->fIsAlive) continue;

	if (fevt->PartInfo(j)->fPID == 2212)
	npr++;
	else if (fevt->PartInfo(j)->fPID == 2112)
	nne++;
	}

	if (fevt->Mode == 2 && npr == 1 && nne == 1) {
	initialstate = 2;

	} else if (fevt->Mode == 2 && ((npr == 0 && nne == 2) \|\| (npr == 2 && nne == 0))) {
	initialstate = 1;
	}
	}


	// Apply weighting
	if (fApply_CCQE && abs(fevt->Mode) == 1) {
	if (fDebugStatements) std::cout << "Getting CCQE Weight" << std::endl;
	double g = GetGausWeight(q0, q3, fGausVal_CCQE);
	if (g < 1.0) g = 1.0;
	rw_weight *= g;
	}

	if (fApply_2p2h && abs(fevt->Mode) == 2) {
	if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl;
	if (fDebugStatements) std::cout << "Got q0 q3 = " << q0 << " " << q3 << " mode = " << fevt->Mode << std::endl;
	rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h);
	if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl;
	}

	if (fApply_2p2h_PPandNN && abs(fevt->Mode) == 2 && initialstate == 1) {
	if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl;
	rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN);
	}

	if (fApply_2p2h_NP && abs(fevt->Mode) == 2 && initialstate == 2) {
	if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl;
	rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP);
	}

	if (fApply_CC1pi && abs(fevt->Mode) >= 11 && abs(fevt->Mode) <= 13) {
	if (fDebugStatements) std::cout << "Getting CC1pi Weight" << std::endl;
	rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi);
	}

	return rw_weight;
	}

	void GaussianModeCorr::SetMethod(bool method) {
	fMethod = method;
	if (fMethod == true) {
	LOG(FIT) << " Using tilt-shift Gaussian parameters for Gaussian enhancement..." << std::endl;
	} else {
	LOG(FIT) << " Using Normal Gaussian parameters for Gaussian enhancement..." << std::endl;
	}
	};

	double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) {
	// The weight
	double w = 1.0;

	// Use tilt-shift method by Patrick
	if (fMethod) {
	if (fDebugStatements) {
	std::cout << "Using Patrick gaussian" << std::endl;
	}
	// // CCQE Without Suppression
	// double Norm = 4.82788679036;
	// double Tilt = 2.3501416116;
	// double Pq0 = 0.363964889702;
	// double Wq0 = 0.133976806938;
	// double Pq3 = 0.431769740224;
	// double Wq3 = 0.207666663434;

	// // Also add for CCQE at the end
	// return (w > 1.0) ? w : 1.0;

	// // 2p2h with suppression
	// double Norm = 15.967;
	// double Tilt = -0.455655;
	// double Pq0 = 0.214598;
	// double Wq0 = 0.0291061;
	// double Pq3 = 0.480194;
	// double Wq3 = 0.134588;

	double Norm = vals[kPosNorm];
	double Tilt = vals[kPosTilt];
	double Pq0 = vals[kPosPq0];
	double Wq0 = vals[kPosWq0];
	double Pq3 = vals[kPosPq3];
	double Wq3 = vals[kPosWq3];

	double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0);
	a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3);

	double b = -sin(2 * Tilt) / (4 * Wq0 * Wq0);
	b += sin(2 * Tilt) / (4 * Wq3 * Wq3);

	double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0);
	c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3);

	w = Norm;
	w = exp(-a (q0 - Pq0) * (q0 - Pq0));
	w = exp(+2.0 b * (q0 - Pq0) * (q3 - Pq3));
	w = exp(-c (q3 - Pq3) * (q3 - Pq3));

	if (fDebugStatements) {
	std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " " << sin(Tilt) << std::endl;
	std::cout << "abc = " << a << " " << b << " " << c << std::endl;
	std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " " << Pq3 << " " << Wq3 << " " << w << std::endl;
	}

	if (w != w \|\| std::isnan(w) \|\| w < 0.0) {
	w = 0.0;
	}

	if (w < 1.0 and !fAllowSuppression) {
	w = 1.0;
	}

	// Use the MINERvA Gaussian method
	} else {
	/*
	* From MINERvA and Daniel Ruterbories:
	* Old notes here: * http://cdcvs.fnal.gov/cgi-bin/public-cvs/cvsweb-public.cgi/AnalysisFramework/Ana/CCQENu/ana_common/data/?cvsroot=mnvsoft
	* These parameters are slightly altered
	*
	* FRESH:
	* 10.5798
	* 0.254032
	* 0.50834
	* 0.0571035
	* 0.129051
	* 0.875287
	*/
	if (fDebugStatements) {
	std::cout << "Using MINERvA Gaussian" << std::endl;
	}

	double norm = vals[kPosNorm];
	double meanq0 = vals[kPosTilt];
	double meanq3 = vals[kPosPq0];
	double sigmaq0 = vals[kPosWq0];
	double sigmaq3 = vals[kPosPq3];
	double corr = vals[kPosWq3];

	double z = (q0 - meanq0)*(q0 - meanq0) /sigmaq0/sigmaq0
	+ (q3 - meanq3)*(q3 - meanq3) / sigmaq3/sigmaq3
	- 2corr(q0-meanq0)(q3-meanq3)/ (sigmaq0 sigmaq3);

	double ret = normexp( -0.5 z / (1 - corr*corr) );
	//Need to add 1 to the results
	w = 1.0 + ret;
	}

	return w;
	}

	void GaussianModeCorr::SetDialValue(std::string name, double val) {
	SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
	}

	void GaussianModeCorr::SetDialValue(int rwenum, double val) {
	int curenum = rwenum % 1000;

	// Check Handled
	if (!IsHandled(curenum)) return;

	// CCQE Setting
	for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) {
	if (i == curenum) {
	int index = i - kGaussianCorr_CCQE_norm;
	fGausVal_CCQE[index] = val;
	fApply_CCQE = true;
	}
	}

	// 2p2h Setting
	for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) {
	if (i == curenum) {
	int index = i - kGaussianCorr_2p2h_norm;
	fGausVal_2p2h[index] = val;
	fApply_2p2h = true;
	}
	}

	// 2p2h_PPandNN Setting
	for (int i = kGaussianCorr_2p2h_PPandNN_norm; i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) {
	if (i == curenum) {
	int index = i - kGaussianCorr_2p2h_PPandNN_norm;
	fGausVal_2p2h_PPandNN[index] = val;
	fApply_2p2h_PPandNN = true;
	}
	}

	// 2p2h_NP Setting
	for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; i++) {
	if (i == curenum) {
	int index = i - kGaussianCorr_2p2h_NP_norm;
	fGausVal_2p2h_NP[index] = val;
	fApply_2p2h_NP = true;
	}
	}

	// CC1pi Setting
	for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) {
	if (i == curenum) {
	int index = i - kGaussianCorr_CC1pi_norm;
	fGausVal_CC1pi[index] = val;
	fApply_CC1pi = true;
	}
	}

	if (curenum == kGaussianCorr_AllowSuppression) {
	fAllowSuppression = (val > 0.5);
	}
	}

	bool GaussianModeCorr::IsHandled(int rwenum) {
	int curenum = rwenum % 1000;
	switch (curenum) {
	case kGaussianCorr_CCQE_norm:
	case kGaussianCorr_CCQE_tilt:
	case kGaussianCorr_CCQE_Pq0:
	case kGaussianCorr_CCQE_Wq0:
	case kGaussianCorr_CCQE_Pq3:
	case kGaussianCorr_CCQE_Wq3:

	case kGaussianCorr_2p2h_norm:
	case kGaussianCorr_2p2h_tilt:
	case kGaussianCorr_2p2h_Pq0:
	case kGaussianCorr_2p2h_Wq0:
	case kGaussianCorr_2p2h_Pq3:
	case kGaussianCorr_2p2h_Wq3:

	case kGaussianCorr_2p2h_PPandNN_norm:
	case kGaussianCorr_2p2h_PPandNN_tilt:
	case kGaussianCorr_2p2h_PPandNN_Pq0:
	case kGaussianCorr_2p2h_PPandNN_Wq0:
	case kGaussianCorr_2p2h_PPandNN_Pq3:
	case kGaussianCorr_2p2h_PPandNN_Wq3:

	case kGaussianCorr_2p2h_NP_norm:
	case kGaussianCorr_2p2h_NP_tilt:
	case kGaussianCorr_2p2h_NP_Pq0:
	case kGaussianCorr_2p2h_NP_Wq0:
	case kGaussianCorr_2p2h_NP_Pq3:
	case kGaussianCorr_2p2h_NP_Wq3:

	case kGaussianCorr_CC1pi_norm:
	case kGaussianCorr_CC1pi_tilt:
	case kGaussianCorr_CC1pi_Pq0:
	case kGaussianCorr_CC1pi_Wq0:
	case kGaussianCorr_CC1pi_Pq3:
	case kGaussianCorr_CC1pi_Wq3:
	case kGaussianCorr_AllowSuppression:
	return true;
	default:
	return false;
	}
	}
	diff --git a/src/Reweight/NUISANCEWeightCalcs.h b/src/Reweight/NUISANCEWeightCalcs.h
	index c09b424..a1d2b9e 100644
	--- a/src/Reweight/NUISANCEWeightCalcs.h
	+++ b/src/Reweight/NUISANCEWeightCalcs.h
	@@ -1,106 +1,129 @@
	#ifndef NUISANCE_WEIGHT_CALCS
	#define NUISANCE_WEIGHT_CALCS

	#include "BaseFitEvt.h"
	+#include "BeRPA.h"

	class NUISANCEWeightCalc {
	public:
	NUISANCEWeightCalc() {};
	virtual ~NUISANCEWeightCalc() {};

	virtual double CalcWeight(BaseFitEvt* evt){return 1.0;};
	virtual void SetDialValue(std::string name, double val){};
	virtual void SetDialValue(int rwenum, double val){};
	virtual bool IsHandled(int rwenum){return false;};

	virtual void Print(){};

	std::map<std::string, int> fDialNameIndex;
	std::map<int, int> fDialEnumIndex;
	std::vector<double> fDialValues;

	std::string fName;
	};

	class ModeNormCalc : public NUISANCEWeightCalc {
	public:
	ModeNormCalc();
	~ModeNormCalc(){};

	double CalcWeight(BaseFitEvt* evt);
	void SetDialValue(std::string name, double val);
	void SetDialValue(int rwenum, double val);
	bool IsHandled(int rwenum);

	double fNormRES;
	};

	+class BeRPACalc : public NUISANCEWeightCalc {
	+ public:
	+ BeRPACalc();
	+ ~BeRPACalc(){};
	+
	+ double CalcWeight(BaseFitEvt* evt);
	+ void SetDialValue(std::string name, double val);
	+ void SetDialValue(int rwenum, double val);
	+ bool IsHandled(int rwenum);
	+
	+ private:
	+ // Parameter values
	+ double fBeRPA_A;
	+ double fBeRPA_B;
	+ double fBeRPA_D;
	+ double fBeRPA_E;
	+ double fBeRPA_U;
	+ // Counts of enabled parameters
	+ int nParams;
	+
	+};
	+
	class SBLOscWeightCalc : public NUISANCEWeightCalc {
	public:
	SBLOscWeightCalc();
	~SBLOscWeightCalc(){};

	double CalcWeight(BaseFitEvt* evt);
	void SetDialValue(std::string name, double val);
	void SetDialValue(int rwenum, double val);
	bool IsHandled(int rwenum);

	double GetSBLOscWeight(double E);

	double fDistance;
	double fMassSplitting;
	double fSin2Theta;

	};

	class GaussianModeCorr : public NUISANCEWeightCalc {
	public:

	GaussianModeCorr();
	~GaussianModeCorr(){};

	double CalcWeight(BaseFitEvt* evt);
	void SetDialValue(std::string name, double val);
	void SetDialValue(int rwenum, double val);
	bool IsHandled(int rwenum);
	double GetGausWeight(double q0, double q3, double vals[]);
	// Set the Gaussian method (tilt-shift or normal Gaussian parameters)
	void SetMethod(bool method);

	// 5 pars describe the Gaussain
	// 0 norm.
	// 1 q0 pos
	// 2 q0 width
	// 3 q3 pos
	// 4 q3 width
	static const int kPosNorm = 0;
	static const int kPosTilt = 1;
	static const int kPosPq0 = 2;
	static const int kPosWq0 = 3;
	static const int kPosPq3 = 4;
	static const int kPosWq3 = 5;

	bool fApply_CCQE;
	double fGausVal_CCQE[6];

	bool fApply_2p2h;
	double fGausVal_2p2h[6];

	bool fApply_2p2h_PPandNN;
	double fGausVal_2p2h_PPandNN[6];

	bool fApply_2p2h_NP;
	double fGausVal_2p2h_NP[6];

	bool fApply_CC1pi;
	double fGausVal_CC1pi[6];

	bool fAllowSuppression;

	bool fDebugStatements;

	bool fMethod;
	};


	#endif
	diff --git a/src/Reweight/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx
	index ebb62fc..39be11c 100644
	--- a/src/Reweight/NUISANCEWeightEngine.cxx
	+++ b/src/Reweight/NUISANCEWeightEngine.cxx
	@@ -1,118 +1,119 @@
	#include "NUISANCEWeightEngine.h"
	#include "NUISANCEWeightCalcs.h"

	#ifdef __MINERVA_RW_ENABLED__
	#ifdef __GENIE_ENABLED__
	#include "MINERvAWeightCalcs.h"
	#endif
	#endif


	NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) {
	// Setup the NUISANCE Reweight engine
	fCalcName = name;
	LOG(FIT) << "Setting up NUISANCE Custom RW : " << fCalcName << std::endl;

	// Load in all Weight Calculations
	GaussianModeCorr* GaussianMode = new GaussianModeCorr();
	std::string Gaussian_Method = FitPar::Config().GetParS("Gaussian_Enhancement");
	if (Gaussian_Method == "Tilt-Shift") {
	GaussianMode->SetMethod(true);
	} else if (Gaussian_Method == "Normal") {
	GaussianMode->SetMethod(false);
	} else {
	ERR(FTL) << "I do not recognise method " << Gaussian_Method << " for the Gaussian enhancement, so will die now..." << std::endl;
	throw;
	}
	fWeightCalculators.push_back(GaussianMode);
	fWeightCalculators.push_back(new ModeNormCalc());
	fWeightCalculators.push_back(new SBLOscWeightCalc());
	+ fWeightCalculators.push_back(new BeRPACalc());

	#ifdef __MINERVA_RW_ENABLED__
	#ifdef __GENIE_ENABLED__
	fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_MEC());
	fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_RES());
	fWeightCalculators.push_back(new nuisance::reweight::RikRPA());
	#endif
	#endif

	// Set Abs Twk Config
	fIsAbsTwk = true;
	};

	void NUISANCEWeightEngine::IncludeDial(std::string name, double startval) {
	// Get First enum
	int nuisenum = Reweight::ConvDial(name, kCUSTOM);

	// Setup Maps
	fEnumIndex[nuisenum]; // = std::vector<size_t>(0);
	fNameIndex[name]; // = std::vector<size_t>(0);

	// Split by commas
	std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
	for (uint i = 0; i < allnames.size(); i++) {
	std::string singlename = allnames[i];

	// Get RW
	int singleenum = Reweight::ConvDial(singlename, kCUSTOM);

	// Fill Maps
	int index = fValues.size();
	fValues.push_back(0.0);
	fNUISANCEEnums.push_back(singleenum);

	// Initialize dial
	LOG(FIT) << "Registering " << singlename << " from " << name << std::endl;

	// Setup index
	fEnumIndex[nuisenum].push_back(index);
	fNameIndex[name].push_back(index);
	}

	// Set Value if given
	if (startval != -999.9) {
	SetDialValue(nuisenum, startval);
	}
	};

	void NUISANCEWeightEngine::SetDialValue(int nuisenum, double val) {
	std::vector<size_t> indices = fEnumIndex[nuisenum];
	for (uint i = 0; i < indices.size(); i++) {
	fValues[indices[i]] = val;
	}
	}

	void NUISANCEWeightEngine::SetDialValue(std::string name, double val) {
	std::vector<size_t> indices = fNameIndex[name];
	for (uint i = 0; i < indices.size(); i++) {
	fValues[indices[i]] = val;
	}
	}

	void NUISANCEWeightEngine::Reconfigure(bool silent) {
	for (size_t i = 0; i < fNUISANCEEnums.size(); i++) {
	for (std::vector<NUISANCEWeightCalc*>::iterator calciter = fWeightCalculators.begin(); calciter != fWeightCalculators.end(); calciter++) {
	NUISANCEWeightCalc* nuiscalc = static_cast<NUISANCEWeightCalc>(calciter);

	if (nuiscalc->IsHandled(fNUISANCEEnums[i])) {
	nuiscalc->SetDialValue(fNUISANCEEnums[i], fValues[i]);
	}
	}
	}
	}

	double NUISANCEWeightEngine::CalcWeight(BaseFitEvt* evt) {
	double rw_weight = 1.0;

	// Cast as usable class
	for (std::vector<NUISANCEWeightCalc*>::iterator iter =
	fWeightCalculators.begin();
	iter != fWeightCalculators.end(); iter++) {
	NUISANCEWeightCalc* nuiscalc = static_cast<NUISANCEWeightCalc>(iter);

	rw_weight *= nuiscalc->CalcWeight(evt);
	}

	// Return rw_weight
	return rw_weight;
	}
	diff --git a/src/Reweight/T2KWeightEngine.cxx b/src/Reweight/T2KWeightEngine.cxx
	index ac0123b..6fadd76 100644
	--- a/src/Reweight/T2KWeightEngine.cxx
	+++ b/src/Reweight/T2KWeightEngine.cxx
	@@ -1,149 +1,140 @@
	#include "T2KWeightEngine.h"

	T2KWeightEngine::T2KWeightEngine(std::string name) {
	#ifdef __T2KREW_ENABLED__

	- // Setup the NEUT Reweight engien
	- fCalcName = name;
	- LOG(FIT) << "Setting up T2K RW : " << fCalcName << std::endl;
	+ // Setup the NEUT Reweight engien
	+ fCalcName = name;
	+ LOG(FIT) << "Setting up T2K RW : " << fCalcName << std::endl;

	- // Create RW Engine suppressing cout
	- StopTalking();
	+ // Create RW Engine suppressing cout
	+ StopTalking();

	- // Create Main RW Engine
	- fT2KRW = new t2krew::T2KReWeight();
	+ // Create Main RW Engine
	+ fT2KRW = new t2krew::T2KReWeight();

	- // Setup Sub RW Engines (Only activated for neut and niwg)
	- fT2KNeutRW = new t2krew::T2KNeutReWeight();
	- fT2KNIWGRW = new t2krew::T2KNIWGReWeight();
	+ // Setup Sub RW Engines (Only activated for neut and niwg)
	+ fT2KNeutRW = new t2krew::T2KNeutReWeight();
	+ fT2KNIWGRW = new t2krew::T2KNIWGReWeight();

	- fT2KRW->AdoptWghtEngine("fNeutRW", fT2KNeutRW);
	- fT2KRW->AdoptWghtEngine("fNIWGRW", fT2KNIWGRW);
	+ fT2KRW->AdoptWghtEngine("fNeutRW", fT2KNeutRW);
	+ fT2KRW->AdoptWghtEngine("fNIWGRW", fT2KNIWGRW);

	- fT2KRW->Reconfigure();
	+ fT2KRW->Reconfigure();

	- // allow cout again
	- StartTalking();
	+ // allow cout again
	+ StartTalking();

	- // Set Abs Twk Config
	- fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
	+ // Set Abs Twk Config
	+ fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));

	#else
	- ERR(FTL) << "T2K RW NOT ENABLED" << std::endl;
	- throw;
	+ ERR(FTL) << "T2K RW NOT ENABLED" << std::endl;
	+ throw;
	#endif
	};

	void T2KWeightEngine::IncludeDial(std::string name, double startval) {
	#ifdef __T2KREW_ENABLED__


	- // Get First enum
	- int nuisenum = Reweight::ConvDial(name, kT2K);
	+ // Get First enum
	+ int nuisenum = Reweight::ConvDial(name, kT2K);

	- // Setup Maps
	- fEnumIndex[nuisenum];// = std::vector<size_t>(0);
	- fNameIndex[name]; // = std::vector<size_t>(0);
	+ // Setup Maps
	+ fEnumIndex[nuisenum];// = std::vector<size_t>(0);
	+ fNameIndex[name]; // = std::vector<size_t>(0);

	- // Split by commas
	- std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
	- for (uint i = 0; i < allnames.size(); i++) {
	- std::string singlename = allnames[i];
	+ // Split by commas
	+ std::vector<std::string> allnames = GeneralUtils::ParseToStr(name, ",");
	+ for (uint i = 0; i < allnames.size(); i++) {
	+ std::string singlename = allnames[i];

	- // Get RW
	- t2krew::T2KSyst_t gensyst = t2krew::T2KSyst::FromString(name);
	+ // Get RW
	+ t2krew::T2KSyst_t gensyst = t2krew::T2KSyst::FromString(name);

	- // Fill Maps
	- int index = fValues.size();
	- fValues.push_back(0.0);
	- fT2KSysts.push_back(gensyst);
	+ // Fill Maps
	+ int index = fValues.size();
	+ fValues.push_back(0.0);
	+ fT2KSysts.push_back(gensyst);

	- // Initialize dial
	- std::cout << "Registering " << singlename << " from " << name << std::endl;
	- fT2KRW->Systematics().Include(gensyst);
	+ // Initialize dial
	+ std::cout << "Registering " << singlename << " from " << name << std::endl;
	+ fT2KRW->Systematics().Include(gensyst);

	- // If Absolute
	- if (fIsAbsTwk) {
	- fT2KRW->Systematics().SetAbsTwk(gensyst);
	- }
	+ // If Absolute
	+ if (fIsAbsTwk) {
	+ fT2KRW->Systematics().SetAbsTwk(gensyst);
	+ }

	- // Setup index
	- fEnumIndex[nuisenum].push_back(index);
	- fNameIndex[name].push_back(index);
	+ // Setup index
	+ fEnumIndex[nuisenum].push_back(index);
	+ fNameIndex[name].push_back(index);

	- }
	+ }

	- // Set Value if given
	- if (startval != -999.9) {
	- SetDialValue(nuisenum, startval);
	- }
	+ // Set Value if given
	+ if (startval != -999.9) {
	+ SetDialValue(nuisenum, startval);
	+ }
	#endif
	}

	void T2KWeightEngine::SetDialValue(int nuisenum, double val) {
	#ifdef __T2KREW_ENABLED__
	- std::vector<size_t> indices = fEnumIndex[nuisenum];
	- for (uint i = 0; i < indices.size(); i++) {
	- fValues[indices[i]] = val;
	- fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
	- }
	+ std::vector<size_t> indices = fEnumIndex[nuisenum];
	+ for (uint i = 0; i < indices.size(); i++) {
	+ fValues[indices[i]] = val;
	+ fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
	+ }
	#endif
	}

	void T2KWeightEngine::SetDialValue(std::string name, double val) {
	#ifdef __T2KREW_ENABLED__
	- std::vector<size_t> indices = fNameIndex[name];
	- for (uint i = 0; i < indices.size(); i++) {
	- fValues[indices[i]] = val;
	- fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
	- }
	+ std::vector<size_t> indices = fNameIndex[name];
	+ for (uint i = 0; i < indices.size(); i++) {
	+ fValues[indices[i]] = val;
	+ fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
	+ }
	#endif
	}


	void T2KWeightEngine::Reconfigure(bool silent) {
	#ifdef __T2KREW_ENABLED__
	- // Hush now...
	- if (silent) StopTalking();
	+ // Hush now...
	+ if (silent) StopTalking();

	- // Reconf
	- StopTalking();
	- fT2KRW->Reconfigure();
	- StartTalking();
	+ // Reconf
	+ StopTalking();
	+ fT2KRW->Reconfigure();
	+ StartTalking();

	- // Shout again
	- if (silent) StartTalking();
	+ // Shout again
	+ if (silent) StartTalking();
	#endif
	}


	double T2KWeightEngine::CalcWeight(BaseFitEvt* evt) {
	- double rw_weight = 1.0;
	+ double rw_weight = 1.0;

	#ifdef __T2KREW_ENABLED__
	- // Skip Non GENIE
	- if (evt->fType != kNEUT) return 1.0;
	+ // Skip Non GENIE
	+ if (evt->fType != kNEUT) return 1.0;

	- // Hush now
	- StopTalking();
	+ // Hush now
	+ StopTalking();

	- // Get Weight For NEUT
	- rw_weight = fT2KRW->CalcWeight(evt->fNeutVect);
	+ // Get Weight For NEUT
	+ rw_weight = fT2KRW->CalcWeight(evt->fNeutVect);

	- // Speak Now
	- StartTalking();
	+ // Speak Now
	+ StartTalking();
	#endif

	- // Return rw_weight
	- return rw_weight;
	+ // Return rw_weight
	+ return rw_weight;
	}

	-
	-
	-
	-
	-
	-
	-
	-
	-
	diff --git a/src/Utils/FitUtils.h b/src/Utils/FitUtils.h
	index 9510e20..ccfe8bd 100644
	--- a/src/Utils/FitUtils.h
	+++ b/src/Utils/FitUtils.h
	@@ -1,186 +1,187 @@
	// 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 FITUTILS_H_SEEN
	#define FITUTILS_H_SEEN

	#include <math.h>
	#include <stdlib.h>
	#include <unistd.h>
	#include <ctime>
	#include <iostream>
	#include <numeric>

	#include <TChain.h>
	#include <TFile.h>
	#include <TH1D.h>
	#include <TH2D.h>
	#include <THStack.h>
	#include <TKey.h>
	#include <TLegend.h>
	#include <TList.h>
	#include <TLorentzVector.h>
	#include <TObjArray.h>
	#include <TROOT.h>
	#include <TRandom3.h>
	#include <TTree.h>
	#include "FitEvent.h"
	#include "TGraph.h"
	#include "TH2Poly.h"
	#include "FitEvent.h"

	#include "FitLogger.h"

	/*!
	* \addtogroup Utils
	* @{
	*/

	/// Functions needed by individual samples for calculating kinematic quantities.
	namespace FitUtils {

	/// Return a vector of all values saved in map
	double *GetArrayFromMap(std::vector<std::string> invals,
	std::map<std::string, double> inmap);

	/// Returns kinetic energy of particle
	double T(TLorentzVector part);

	/// Returns momentum of particle
	double p(TLorentzVector part);
	double p(FitParticle* part);

	/// Returns angle between particles (_NOT_ cosine!)
	double th(TLorentzVector part, TLorentzVector part2);
	double th(FitParticle* part1, FitParticle* part2);

	/// Hadronic mass reconstruction
	double Wrec(TLorentzVector pnu, TLorentzVector pmu);

	/// Hadronic mass true from initial state particles and muon; useful if the full
	/// FSI vectors aren't not saved and we for some reasons need W_true
	double Wtrue(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector pnuc);

	double SumKE_PartVect(std::vector<FitParticle *> const fps);
	double SumTE_PartVect(std::vector<FitParticle *> const fps);

	/// Return E Hadronic for all FS Particles in Hadronic System
	double GetErecoil_TRUE(FitEvent *event);

	/// Return E Hadronic for all Charged FS Particles in Hadronic System
	double GetErecoil_CHARGED(FitEvent *event);
	double Eavailable(FitEvent *event);

	/*
	CCQE MiniBooNE/MINERvA
	*/
	/// Function to calculate the reconstructed Q^{2}_{QE}
	double Q2QErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);

	/// Function returns the reconstructed E_{nu} values
	double EnuQErec(TLorentzVector pmu, double costh, double binding, bool neutrino = true);

	//! Function to calculate the reconstructed Q^{2}_{QE}
	double Q2QErec(double pl, double costh, double binding, bool neutrino = true);

	//! Function to calculate the reconstructed Q^{2}_{QE}
	double Q2QErec(TLorentzVector Pmu, TLorentzVector Pnu, double binding, bool neutrino = true);

	//! Function returns the reconstructed E_{nu} values
	double EnuQErec(double pl, double costh, double binding,
	bool neutrino = true);

	/*
	CCQE1p MINERvA
	*/
	/// Reconstruct Q2QE given just the maximum energy proton.
	double ProtonQ2QErec(double pE, double binding);

	/*
	E Recoil MINERvA
	*/
	double GetErecoil_MINERvA_LowRecoil(FitEvent *event);

	/*
	CC1pi0 MiniBooNE
	*/
	/// Reconstruct Enu from CCpi0 vectors and binding energy
	double EnuCC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
	TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));

	/// Reconstruct Q2 from CCpi0 vectors and binding energy
	double Q2CC1pi0rec(TLorentzVector pnu, TLorentzVector pmu,
	TLorentzVector ppi0 = TLorentzVector(0, 0, 0, 0));

	/*
	CC1pi+ MiniBooNE
	*/

	/// returns reconstructed Enu a la MiniBooNE CCpi+
	/// returns reconstructed Enu a la MiniBooNE CCpi+
	// Also for when not having pion info (so when we have a Michel tag in T2K)
	double EnuCC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
	bool pionInfo = true);

	/// returns reconstructed Enu assumming resonance interaction where intermediate
	/// resonance was a Delta
	double EnuCC1piprecDelta(TLorentzVector pnu, TLorentzVector pmu);

	/// returns reconstructed in a variety of flavours
	double Q2CC1piprec(TLorentzVector pnu, TLorentzVector pmu, TLorentzVector ppip,
	int enuType = 0, bool pionInfo = true);

	/*
	T2K CC1pi+ on CH
	*/
	double thq3pi_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
	TLorentzVector ppi);
	double q3_CC1pip_T2K(TLorentzVector pnu, TLorentzVector pmu,
	TLorentzVector ppi);
	double WrecCC1pip_T2K_MB(TLorentzVector pnu, TLorentzVector pmu,
	TLorentzVector ppip);
	double EnuCC1piprec_T2K_eMB(TLorentzVector pnu, TLorentzVector pmu,
	TLorentzVector ppi);

	/*
	nucleon single pion
	*/
	double MpPi(TLorentzVector pp, TLorentzVector ppi);

	/// Gets delta p T as defined in Phys.Rev. C94 (2016) no.1, 015503
	double Get_STV_dpt(FitEvent *event, int ISPDG, bool Is0pi);
	/// Gets delta phi T as defined in Phys.Rev. C94 (2016) no.1, 015503
	double Get_STV_dphit(FitEvent *event, int ISPDG, bool Is0pi);
	/// Gets delta alpha T as defined in Phys.Rev. C94 (2016) no.1, 015503
	double Get_STV_dalphat(FitEvent *event, int ISPDG, bool Is0pi);

	// As defined in PhysRevC.95.065501
	// Using prescription from arXiv 1805.05486
	double Get_pn_reco_C(FitEvent *event, int ISPDG, bool Is0pi);
	+double Get_pn_reco_Ar(FitEvent *event, int ISPDG, bool Is0pi);

	//For T2K inferred kinematics analyis - variables defined as on page 7 of T2K TN287v11 (and now arXiv 1802.05078)
	double ppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
	TVector3 tppInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);
	double cthpInfK(TLorentzVector pmu, double costh, double binding, bool neutrino);

	double CosThAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
	double PhiAdler(TLorentzVector Pnu, TLorentzVector Pmu, TLorentzVector Ppi, TLorentzVector Pprot);
	}

	/! @} /
	#endif

File Metadata

Mime Type: text/x-diff
Expires: Sat, May 3, 6:27 AM (1 d, 17 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4902781
Default Alt Text: (43 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions