No OneTemporary
Actions

Size

47 KB

Subscribers

None

View Options

	diff --git a/parameters/config.xml b/parameters/config.xml
	index 8ad5bc2..ef6b88e 100644
	--- a/parameters/config.xml
	+++ b/parameters/config.xml
	@@ -1,200 +1,210 @@
	<nuisance>
	<!-- # ###################################################### -->
	<!-- # NUISANCE CONFIGURATION OPTIONS -->
	<!-- # This file is read in by default at runtime -->
	<!-- # If you want to override on a case by case bases use -q at runtime -->
	<!-- # ###################################################### -->

	<!-- # MAIN Configs -->
	<!-- # ###################################################### -->

	<!-- # Logger goes from -->
	<!-- # 1 Quiet -->
	<!-- # 2 Fitter -->
	<!-- # 3 Samples -->
	<!-- # 4 Reconfigure Loops -->
	<!-- # 5 Every Event print out (SHOUT) -->
	<!-- # -1 DEBUGGING -->
	<config verbosity='5'/>
	<config VERBOSITY='5'/>

	<!-- # ERROR goes from -->
	<!-- # 0 NONE -->
	<!-- # 1 FATAL -->
	<!-- # 2 WARN -->
	<config ERROR='2'/>
	<config TRACE='0'/>

	<config cores='1' />
	<config spline_test_throws='50' />
	<config spline_cores='1' />
	<config spline_chunks='20' />
	<config spline_procchunk='-1' />

	<config Electron_NThetaBins='4' />
	<config Electron_NEnergyBins='4' />
	<config Electron_ThetaWidth='1.0' />
	<config Electron_EnergyWidth='0.10' />

	<config RemoveFSIParticles='0' />
	<config RemoveUndefParticles='0' />
	<config RemoveNuclearParticles='0'/>
	<config logging.JointFCN.cxx='4'/>
	<config MINERvASaveExtraCCQE='0' />

	<!-- # Input Configs -->
	<!-- # ###################################################### -->

	<!-- # Default Requirements file for the externalDataFitter Package -->
	<!-- # MAX Events : -1 is no cut. Events will be scaled automatically to give good xsec predictions. -->
	<config input.maxevents='-1'/>
	<config MAXEVENTS='-1'/>
	<config input.MAXEVENTS='-1'/>
	<config includeemptystackhists='0'/>
	<!-- # Turn on/off event manager -->
	<!-- # EventManager enables us to only loop number of events once for multiple projections of the same measurements -->
	<!-- # e.g. MiniBooNE CC1pi+ Q2 and MiniBooNE CC1pi+ Tmu would ordinarily require 2 reconfigures, but with this enabled it requires only one -->
	<config input.eventmanager='1'/>
	<config EventManager='1'/>

	<!-- # Event Directories -->
	<!-- # Can setup default directories and use @EVENT_DIR/path to link to it -->
	<config EVENT_DIR='/data2/stowell/NIWG/'/>
	<config NEUT_EVENT_DIR='/data2/stowell/NIWG/neut/fit_samples_neut5.3.3/'/>
	<config GENIE_EVENT_DIR='/data2/stowell/NIWG/genie/fit_samples_R.2.10.0/'/>
	<config NUWRO_EVENT_DIR='/data2/stowell/NIWG/nuwro/fit_samples/'/>
	<config GIBUU_EVENT_DIR='/data/GIBUU/DIR/'/>
	<config SaveNuWroExtra='0' />

	<!-- # In PrepareGENIE the reconstructed splines can be saved into the file -->
	<config save_genie_splines='1'/>

	<!-- # In InputHandler the option to regenerate NuWro flux/xsec plots is available -->
	<!-- # Going to move this to its own app soon -->
	<config input.regen_nuwro_plots='0'/>

	<!-- # DEVEL CONFIG OPTION, don't touch! -->
	<config CacheSize='0'/>

	<!-- # ReWeighting Configuration Options -->
	<!-- # ###################################################### -->

	<!-- # Set absolute twkdial for parameters -->
	<config params.setabstwk='0'/>

	<!-- # Convert Dials in output statements using dial conversion card -->
	<config convert_dials='0'/>

	<!-- # Make RW Calculations be quiet -->
	<config params.silentweighting='0'/>

	<!-- # Vetos can be used to specify RW dials NOT to be loaded in -->
	<!-- # Useful if one specific one has an issue -->
	<config FitWeight.fNIWGRW_veto=''/>
	<config FitWeight.fNuwroRW_veto=''/>
	<config FitWeight.fNeutRW_veto=''/>
	<config FitWeight.fGenieRW_veto=''/>


	<!-- # Output Options -->
	<!-- # ###################################################### -->

	<!-- # Save Nominal prediction with all rw engines at default -->
	<config savenominal='0'/>

	<!-- # Save prefit with values at starting values -->
	<config saveprefit='0'/>

	<!-- # Here's the full list of drawing options -->
	<!-- # See src/FitBase/Measurement1D::Write for more info -->
	<!-- #config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/RESIDUAL/MATRIX/FLUX/MASK/MAP -->
	<!-- #config drawopts DATA/MC -->
	<config drawopts='DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC/SETTINGS'/>

	<config InterpolateSigmaQ0Histogram='1' />
	<config InterpolateSigmaQ0HistogramRes='100' />
	<config InterpolateSigmaQ0HistogramThrow='1' />
	<config InterpolateSigmaQ0HistogramNTHROWS='100000' />

	<!-- # Save the shape scaling applied with option SHAPE into the main MC hist -->
	<config saveshapescaling='0'/>

	<config CorrectGENIEMECNorm='1'/>

	<!-- # Set style of 1D output histograms -->
	<config linecolour='1'/>
	<config linestyle='1'/>
	<config linewidth='1'/>

	<!-- # For GenericFlux -->
	<config isLiteMode='0'/>

	<!-- # Statistical Options -->
	<!-- # ###################################################### -->

	<!-- # Add MC Statistical error to likelihoods -->
	<config statutils.addmcerror='0'/>

	<!-- # NUISMIN Configurations -->
	<!-- # ###################################################### -->

	<config minimizer.maxcalls='1000000'/>
	<config minimizer.maxiterations='1000000'/>
	<config minimizer.tolerance='0.001'/>

	<!-- # Number of events required in low stats routines -->
	<config minimizer.lowstatevents='25000'/>


	<!-- # Error band generator configs -->
	<!-- # ###################################################### -->

	<!-- # For -f ErrorBands creates error bands for given measurements -->
	<!-- # How many throws do we want (The higher the better precision) -->
	<config error_throws='250'/>

	<!-- # Are we throwing uniform or according to Gaussian? -->
	<!-- # Only use uniform if wanting to study the limits of a dial. -->
	<config error_uniform='0'/>
	<config WriteSeperateStacks='1'/>

	<!-- # Other Individual Case Configs -->
	<!-- # ###################################################### -->

	<!-- # Covariance throw options for fake data studies with MiniBooNE data. -->
	<config thrown_covariance='FULL'/>
	<config throw_mc_stat='0.0'/>
	<config throw_diag_syst='0'/>
	<config throw_corr_syst='0'/>
	<config throw_mc_stat='0'/>

	<!-- # Apply a shift to the muon momentum before calculation of Q2 -->
	<config muon_momentum_shift='0.0'/>
	<config muon_momentum_throw='0'/>

	<!-- # MINERvA Specific Configs -->
	<config MINERvA_XSec_CCinc_2DEavq3_nu.hadron_cut='0'/>
	<config MINERvA_CCinc_XSec_2DEavq3_nu.useq3true='0'/>
	<config Modes.split_PN_NN='0'/>
	<config SignalReconfigures='0'/>

	<!-- # SciBooNE specific -->
	<config SciBarDensity='1.04'/>
	<config SciBarRecoDist='10.0'/>
	<config PenetratingMuonEnergy='1.4'/>
	<config NumRangeSteps='50'/>

	<config MINERvADensity='1.04'/>
	<config MINERvARecoDist='10.0'/>
	<config PenetratingMuonEnergy='1.4'/>
	<config NumRangeSteps='50'/>



	<config GENIEWeightEngine_CCRESMode="kModeMaMv"/>
	<config GENIEInputHandler.SavePrimary="0" />

	<config NToyThrows='100000' />

	+<config GENIEXSecModelCCRES="genie::ReinSehgalRESPXSec" />
	+<config GENIEXSecModelCOH="genie::ReinSehgalCOHPiPXSec" />
	+<config GENIEXSecModelCCQE="genie::LwlynSmithQELCCPXSec" />
	+
	+<!--
	+<config GENIEXSecModelCCQE="genie::NievesQELCCPXSec" />
	+<config GENIEXSecModelCCRES="genie::BergerSehgalRESPXSec2014" />
	+<config GENIEXSecModelCOH="genie::BergerSehgalCOHPiPXSec2015" />
	+-->
	+
	</nuisance>

	diff --git a/src/MCStudies/MCStudy_CCQEHistograms.cxx b/src/MCStudies/MCStudy_CCQEHistograms.cxx
	index bde660b..13298ee 100644
	--- a/src/MCStudies/MCStudy_CCQEHistograms.cxx
	+++ b/src/MCStudies/MCStudy_CCQEHistograms.cxx
	@@ -1,225 +1,230 @@
	// 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 "MCStudy_CCQEHistograms.h"

	#include "T2K_SignalDef.h"
	#include "MINERvA_SignalDef.h"

	//********************************************************************
	/// @brief Class to perform MC Studies on a custom measurement
	MCStudy_CCQEHistograms::MCStudy_CCQEHistograms(std::string name, std::string inputfile,
	FitWeight *rw, std::string type,
	std::string fakeDataFile) {
	//********************************************************************

	// Measurement Details
	fName = name;
	fEventTree = NULL;

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

	// 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);
	fEventTree = NULL;

	// Setup fDataHist as a placeholder
	this->fDataHist = new TH1D(("approximate_data"), ("kaon_data"), 5, 1.0, 6.0);

	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.
	this->fScaleFactor = (this->fEventHist->Integral("width") * 1E-38 / (fNEvents + 0.)) /
	this->TotalIntegratedFlux();


	hist_Enu = new TH1D("MCStudy_CCQE_Enu","MCStudy_CCQE_Enu",30,0.0,2.0);
	hist_TLep = new TH1D("MCStudy_CCQE_TLep","MCStudy_CCQE_TLep",30,0.0,4.0);
	hist_CosLep = new TH1D("MCStudy_CCQE_CosLep","MCStudy_CCQE_CosLep",30,-1.0,1.0);
	- hist_Q2 = new TH1D("MCStudy_CCQE_Q2","MCStudy_CCQE_Q2",30,0.0,2.0);
	- hist_Q2QE = new TH1D("MCStudy_CCQE_Q2QE","MCStudy_CCQE_Q2QE",30,0.0,2.0);
	+ hist_Q2 = new TH1D("MCStudy_CCQE_Q2;Q^{2} (GeV^{2});d#sigma/dQ^{2} (cm^{2}/nucleon/GeV^{2})","MCStudy_CCQE_Q2",30,0.0,3.0);
	+ hist_Q2QE = new TH1D("MCStudy_CCQE_Q2QE","MCStudy_CCQE_Q2QE",30,0.0,3.0);
	hist_EQE = new TH1D("MCStudy_CCQE_EQE","MCStudy_CCQE_EQE",30,0.0,5.0);
	hist_q0 = new TH1D("MCStudy_CCQE_q0","MCStudy_CCQE_q0",30,0.0,2.0);
	hist_q3 = new TH1D("MCStudy_CCQE_q3","MCStudy_CCQE_q3",30,0.0,2.0);
	hist_TLepCosLep = new TH2D("MCStudy_CCQE_TLepCosLep","MCStudy_CCQE_TLepCosLep",15,0.0,5.0,15,-1.0,1.0);
	hist_Total = new TH1D("MCStudy_CCQE_TotalXSec","MXStudy_CCQE_TotalXSec",1,0.0,1.0);

	+ hist_q0q3 = new TH2D("MCStudy_CCQE_q0q3","MCStudy_CCQE_q0q3;q_{3} (GeV); q_{0} (GeV); d#sigma/dq_{0}dq_{3} (cm^{2}/nucleon/GeV^{2})",40,0.0,2.0,40,0.0,2.0);
	+
	return;
	}

	//********************************************************************
	void MCStudy_CCQEHistograms::FillEventVariables(FitEvent *event) {
	//********************************************************************

	// std::cout << "Event fBound = " << event->fBound << " " << event->Mode << std::endl;
	- if (event->fBound > 0) return;
	+// if (event->fBound > 0) return;
	if (abs(event->Mode) != 1) return;
	- std::cout << "Event fBound = " << event->fBound << " " << event->Mode << "-> Signal " << std::endl;
	+ // std::cout << "Event fBound = " << event->fBound << " " << event->Mode << "-> Signal " << std::endl;


	FitParticle* muon = NULL;
	FitParticle* nu = event->GetNeutrinoIn();
	bool IsNuMu = event->PDGnu() > 0;

	if (IsNuMu) muon = event->GetHMFSParticle(13);
	else muon = event->GetHMFSParticle(-13);

	// Reset Variables
	Enu = -999.9;
	TLep = -999.9;
	CosLep = -999.9;
	Q2 = -999.9;
	Q2QE = -999.9;
	EQE = -999.9;
	q0 = -999.9;
	q3 = -999.9;

	// Fill Variables
	if (muon){

	Enu = event->Enu() / 1.E3;

	TLep = (muon->fP.E() - muon->fP.Mag()) / 1.E3;
	CosLep = cos(muon->fP.Vect().Angle( nu->fP.Vect() ));

	Q2 = fabs((muon->fP - nu->fP).Mag2() / 1.E6);

	Q2QE = FitUtils::Q2QErec(muon->fP, CosLep, 34., IsNuMu);
	EQE = FitUtils::EnuQErec(muon->fP, CosLep, 34., IsNuMu);

	q0 = fabs((muon->fP - nu->fP).E()) / 1.E3;
	q3 = fabs((muon->fP - nu->fP).Vect().Mag()) / 1.E3;

	LocalRWWeight = event->RWWeight;
	LocalInputWeight = event->InputWeight;

	}

	// Fill Tree
	if (abs(Mode) == 1 and Signal){
	hist_Enu->Fill(Enu,event->Weight);
	hist_TLep->Fill(TLep,event->Weight);
	hist_CosLep->Fill(CosLep,event->Weight);
	hist_Q2->Fill(Q2,event->Weight);
	hist_Q2QE->Fill(Q2QE,event->Weight);
	hist_EQE->Fill(EQE,event->Weight);
	hist_q0->Fill(q0,event->Weight);
	hist_q3->Fill(q3,event->Weight);
	hist_TLepCosLep->Fill(TLep,CosLep,event->Weight);
	-
	+ hist_q0q3->Fill(q3,q0,event->Weight);
	hist_Total->Fill(0.5,event->Weight);

	fXVar = Q2;
	}

	return;
	};

	//********************************************************************
	void MCStudy_CCQEHistograms::Write(std::string drawOpt) {
	//********************************************************************
	// Measurement1D::Write(drawOpt);

	LOG(FIT) << "Writing MCStudy_CCQEHistograms " << std::endl;
	// FitPar::Config().out->cd();
	hist_Enu->Write();
	hist_TLep->Write();
	hist_CosLep->Write();
	hist_Q2->Write();
	hist_Q2QE->Write();
	hist_EQE->Write();
	hist_q0->Write();
	hist_q3->Write();
	hist_TLepCosLep->Write();
	+ hist_q0q3->Write();
	hist_Total->Write();

	return;
	}

	//********************************************************************
	void MCStudy_CCQEHistograms::ResetAll(){
	//********************************************************************
	hist_Enu->Reset();
	hist_TLep->Reset();
	hist_CosLep->Reset();
	hist_Q2->Reset();
	hist_Q2QE->Reset();
	hist_EQE->Reset();
	hist_q0->Reset();
	hist_q3->Reset();
	+ hist_q0q3->Reset();
	hist_TLepCosLep->Reset();
	hist_Total->Reset();

	return;
	}


	//********************************************************************
	void MCStudy_CCQEHistograms::ScaleEvents(){
	//********************************************************************

	hist_Enu->Scale(fScaleFactor,"width");
	hist_TLep->Scale(fScaleFactor,"width");
	hist_CosLep->Scale(fScaleFactor,"width");
	hist_Q2->Scale(fScaleFactor,"width");
	hist_Q2QE->Scale(fScaleFactor,"width");
	hist_EQE->Scale(fScaleFactor,"width");
	hist_q0->Scale(fScaleFactor,"width");
	hist_q3->Scale(fScaleFactor,"width");
	+ hist_q0q3->Scale(fScaleFactor,"width");
	hist_TLepCosLep->Scale(fScaleFactor,"width");
	hist_Total->Scale(fScaleFactor,"width");

	return;
	}



	//********************************************************************
	/// Select only events with final state Muons
	bool MCStudy_CCQEHistograms::isSignal(FitEvent *event) {
	//********************************************************************


	if (abs(event->Mode) != 1) return false;
	- if (event->fBound > 0) return false;
	+ //if (event->fBound > 0) return false;
	// if (!event->HasFSMuon()) return false;

	// Do we want any other signal?
	return true;
	};



	diff --git a/src/MCStudies/MCStudy_CCQEHistograms.h b/src/MCStudies/MCStudy_CCQEHistograms.h
	index 455a9f2..4799ef3 100644
	--- a/src/MCStudies/MCStudy_CCQEHistograms.h
	+++ b/src/MCStudies/MCStudy_CCQEHistograms.h
	@@ -1,75 +1,76 @@
	// 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 MCStudy_CCQEHistograms_H_SEEN
	#define MCStudy_CCQEHistograms_H_SEEN
	#include "Measurement1D.h"

	//********************************************************************
	class MCStudy_CCQEHistograms : public Measurement1D {
	//********************************************************************

	public:

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

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

	void ScaleEvents();
	void ResetAll();

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

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

	private:

	double fEventScaleFactor;
	TTree* fEventTree;
	TH1D* hist_Enu;
	float Enu;
	TH1D* hist_TLep;
	float TLep ;
	TH1D* hist_CosLep;
	float CosLep;
	TH1D* hist_Q2;
	float Q2 ;
	TH1D* hist_Q2QE;
	float Q2QE ;
	TH1D* hist_EQE;
	float EQE ;
	TH1D* hist_q0;
	float q0 ;
	TH1D* hist_q3;
	float q3 ;
	+ TH2D* hist_q0q3;

	TH1D* hist_Total;

	TH2D* hist_TLepCosLep;

	double LocalRWWeight;
	double LocalInputWeight;

	};

	#endif
	diff --git a/src/Reweight/CMakeLists.txt b/src/Reweight/CMakeLists.txt
	index 685133c..55b9fef 100644
	--- a/src/Reweight/CMakeLists.txt
	+++ b/src/Reweight/CMakeLists.txt
	@@ -1,81 +1,83 @@
	# 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
	)

	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/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx
	index 4b4ceb1..f5d769b 100644
	--- a/src/Reweight/GENIEWeightEngine.cxx
	+++ b/src/Reweight/GENIEWeightEngine.cxx
	@@ -1,231 +1,242 @@
	#include "GENIEWeightEngine.h"

	GENIEWeightEngine::GENIEWeightEngine(std::string name) {
	#ifdef __GENIE_ENABLED__

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

	// Create RW Engine suppressing cout
	StopTalking();
	fGenieRW = new genie::rew::GReWeight();

	// Get List of Vetos (Just for debugging)
	std::string rw_engine_list = FitPar::Config().GetParS("FitWeight.fGenieRW_veto");
	bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
	bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
	bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
	bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos;
	bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos;
	bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos;
	bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos;
	bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos;
	bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos;
	bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos;
	bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
	bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
	bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
	bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
	bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos;
	bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos;

	// Now actually add the RW Calcs
	if (xsec_ncel)
	fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL);
	- if (xsec_ccqe)
	- fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE);
	- if (xsec_coh)
	- fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH);
	+ if (xsec_ccqe){
	+ fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE);
	+ // (dynamic_cast<GReWeightNuXSecCCQE*> (fGenieRW->WghtCalc("xsec_ccqe")))
	+ // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") );
	+ }
	+ if (xsec_coh){
	+ fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH());
	+ // (dynamic_cast<GReWeightNuXSecCOH*> (fGenieRW->WghtCalc("xsec_coh")))
	+ // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") );
	+ }
	+
	if (xsec_nnres)
	fGenieRW->AdoptWghtCalc("xsec_nonresbkg",
	new genie::rew::GReWeightNonResonanceBkg);
	if (xsec_nudis)
	fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod);
	if (xsec_resdec)
	fGenieRW->AdoptWghtCalc("hadro_res_decay",
	new genie::rew::GReWeightResonanceDecay);
	if (xsec_fzone)
	fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone);
	if (xsec_intra)
	fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke);
	if (xsec_agky)
	fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY);
	if (xsec_qevec)
	fGenieRW->AdoptWghtCalc("xsec_ccqe_vec",
	new genie::rew::GReWeightNuXSecCCQEvec);
	#if __GENIE_VERSION__ >= 212
	if (xsec_qeaxial)
	fGenieRW->AdoptWghtCalc("xsec_ccqe_axial",
	new genie::rew::GReWeightNuXSecCCQEaxial);
	#endif

	if (xsec_dis)
	fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS);
	if (xsec_nc)
	fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC);
	- if (xsec_ccres)
	+ if (xsec_ccres){
	fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES);
	+ // (dynamic_cast<GReWeightNuXSecCCRES*> (fGenieRW->WghtCalc("xsec_ccres")))
	+ // ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCRES") );
	+ }
	+
	if (xsec_ncres)
	fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES);
	if (xsec_nucqe)
	fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM);

	GReWeightNuXSecCCQE * rwccqe =
	dynamic_cast<GReWeightNuXSecCCQE *> (fGenieRW->WghtCalc("xsec_ccqe"));
	rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);

	// Default to include shape and normalization changes for CCRES (can be changed downstream if desired)
	GReWeightNuXSecCCRES * rwccres =
	dynamic_cast<GReWeightNuXSecCCRES *> (fGenieRW->WghtCalc("xsec_ccres"));
	std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode");
	if (!marestype.compare("kModeNormAndMaMvShape")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); }
	else if (!marestype.compare("kModeMaMv")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); }
	else {
	THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype );
	}

	// Default to include shape and normalization changes for NCRES (can be changed downstream if desired)
	GReWeightNuXSecNCRES * rwncres =
	dynamic_cast<GReWeightNuXSecNCRES *> (fGenieRW->WghtCalc("xsec_ncres"));
	rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);

	// Default to include shape and normalization changes for DIS (can be changed downstream if desired)
	GReWeightNuXSecDIS * rwdis =
	dynamic_cast<GReWeightNuXSecDIS *> (fGenieRW->WghtCalc("xsec_dis"));
	rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u);

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

	// allow cout again
	StartTalking();

	#else
	ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl;
	#endif
	};


	void GENIEWeightEngine::IncludeDial(std::string name, double startval) {
	#ifdef __GENIE_ENABLED__

	// Get First enum
	int nuisenum = Reweight::ConvDial(name, kGENIE);

	// 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
	genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename);

	// Fill Maps
	int index = fValues.size();
	fValues.push_back(0.0);
	fGENIESysts.push_back(rwsyst);

	// Initialize dial
	std::cout << "Registering " << singlename << " from " << name << std::endl;
	fGenieRW->Systematics().Init( fGENIESysts[index] );

	// If Absolute
	if (fIsAbsTwk) {
	GSystUncertainty::Instance()->SetUncertainty( rwsyst, 1.0, 1.0 );
	}

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

	}

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


	void GENIEWeightEngine::SetDialValue(int nuisenum, double val) {
	#ifdef __GENIE_ENABLED__
	std::vector<size_t> indices = fEnumIndex[nuisenum];
	for (uint i = 0; i < indices.size(); i++) {
	fValues[indices[i]] = val;
	fGenieRW->Systematics().Set( fGENIESysts[indices[i]], val);
	}
	#endif
	}

	void GENIEWeightEngine::SetDialValue(std::string name, double val) {
	#ifdef __GENIE_ENABLED__
	std::vector<size_t> indices = fNameIndex[name];
	for (uint i = 0; i < indices.size(); i++) {
	fValues[indices[i]] = val;
	fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
	}
	#endif
	}

	void GENIEWeightEngine::Reconfigure(bool silent) {
	#ifdef __GENIE_ENABLED__
	// Hush now...
	if (silent) StopTalking();

	// Reconf
	fGenieRW->Reconfigure();
	fGenieRW->Print();

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


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

	#ifdef __GENIE_ENABLED__
	// Skip Non GENIE
	if (evt->fType != kGENIE) return 1.0;

	// Make nom weight
	if (!evt) {
	THROW("evt not found : " << evt);
	}

	if (!(evt->genie_event)) {
	THROW("evt->genie_event not found!" << evt->genie_event);
	}

	if (!(evt->genie_event->event)) {
	THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event));
	}

	if (!fGenieRW) {
	THROW("GENIE RW Not Found!" << fGenieRW);
	}

	rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event));
	// std::cout << "Returning GENIE Weight for electron scattering = " << rw_weight << std::endl;
	#endif

	// Return rw_weight
	return rw_weight;
	}










	diff --git a/src/Reweight/MINERvAWeightCalcs.cxx b/src/Reweight/MINERvAWeightCalcs.cxx
	index ebd453a..6cbc7ef 100644
	--- a/src/Reweight/MINERvAWeightCalcs.cxx
	+++ b/src/Reweight/MINERvAWeightCalcs.cxx
	@@ -1,469 +1,478 @@
	#ifdef __MINERVA_RW_ENABLED__
	#ifdef __GENIE_ENABLED__

	#include "MINERvAWeightCalcs.h"
	#include "BaseFitEvt.h"

	namespace nuisance {
	namespace reweight {


	//*******************************************************
	MINERvAReWeight_QE::MINERvAReWeight_QE() {
	//*******************************************************
	fTwk_NormCCQE = 0.0;
	fDef_NormCCQE = 1.0;
	fCur_NormCCQE = fDef_NormCCQE;
	}

	MINERvAReWeight_QE::~MINERvAReWeight_QE(){};

	double MINERvAReWeight_QE::CalcWeight(BaseFitEvt* evt) {
	// Check GENIE
	if (evt->fType != kGENIE) return 1.0;

	// Extract the GENIE Record
	GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
	const Interaction* interaction = ghep->Summary();
	const InitialState& init_state = interaction->InitState();
	const ProcessInfo& proc_info = interaction->ProcInfo();
	const Target& tgt = init_state.Tgt();

	// If the event is not QE this Calc doesn't handle it
	if (!proc_info.IsQuasiElastic()) return 1.0;

	// WEIGHT CALCULATIONS -------------
	double w = 1.0;

	// CCQE Dial
	if (!proc_info.IsWeakCC()) w *= fCur_NormCCQE;

	// Return Combined Weight
	return w;
	}

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

	void MINERvAReWeight_QE::SetDialValue(int rwenum, double val) {
	// Check Handled
	int curenum = rwenum % 1000;
	if (!IsHandled(curenum)) return;

	// Set Values
	if (curenum == kMINERvARW_NormCCQE) {
	fTwk_NormCCQE = val;
	fCur_NormCCQE = fDef_NormCCQE + fTwk_NormCCQE;
	}

	// Define Tweaked
	fTweaked = ((fTwk_NormCCQE != 0.0));
	}

	bool MINERvAReWeight_QE::IsHandled(int rwenum) {
	int curenum = rwenum % 1000;

	switch (curenum) {
	case kMINERvARW_NormCCQE:
	return true;
	default:
	return false;
	}
	}


	//*******************************************************
	MINERvAReWeight_MEC::MINERvAReWeight_MEC() {
	//*******************************************************
	fTwk_NormCCMEC = 0.0;
	fDef_NormCCMEC = 1.0;
	fCur_NormCCMEC = fDef_NormCCMEC;
	}

	MINERvAReWeight_MEC::~MINERvAReWeight_MEC(){};

	double MINERvAReWeight_MEC::CalcWeight(BaseFitEvt* evt) {
	// Check GENIE
	if (evt->fType != kGENIE) return 1.0;

	// Extract the GENIE Record
	GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
	const Interaction* interaction = ghep->Summary();
	const InitialState& init_state = interaction->InitState();
	const ProcessInfo& proc_info = interaction->ProcInfo();
	const Target& tgt = init_state.Tgt();

	// If the event is not MEC this Calc doesn't handle it
	if (!proc_info.IsMEC()) return 1.0;

	// WEIGHT CALCULATIONS -------------
	double w = 1.0;

	// CCMEC Dial
	if (!proc_info.IsWeakCC()) w *= fCur_NormCCMEC;

	// Return Combined Weight
	return w;
	}

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

	void MINERvAReWeight_MEC::SetDialValue(int rwenum, double val) {
	// Check Handled
	int curenum = rwenum % 1000;
	if (!IsHandled(curenum)) return;

	// Set Values
	if (curenum == kMINERvARW_NormCCMEC) {
	fTwk_NormCCMEC = val;
	fCur_NormCCMEC = fDef_NormCCMEC + fTwk_NormCCMEC;
	}

	// Define Tweaked
	fTweaked = ((fTwk_NormCCMEC != 0.0));
	}

	bool MINERvAReWeight_MEC::IsHandled(int rwenum) {
	int curenum = rwenum % 1000;

	switch (curenum) {
	case kMINERvARW_NormCCMEC:
	return true;
	default:
	return false;
	}
	}

	//*******************************************************
	MINERvAReWeight_RES::MINERvAReWeight_RES() {
	//*******************************************************
	fTwk_NormCCRES = 0.0;
	fDef_NormCCRES = 1.0;
	fCur_NormCCRES = fDef_NormCCRES;
	}

	MINERvAReWeight_RES::~MINERvAReWeight_RES(){};

	double MINERvAReWeight_RES::CalcWeight(BaseFitEvt* evt) {

	// std::cout << "Caculating RES" << std::endl;
	// Check GENIE
	if (evt->fType != kGENIE) return 1.0;

	// Extract the GENIE Record
	GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
	const Interaction* interaction = ghep->Summary();
	const InitialState& init_state = interaction->InitState();
	const ProcessInfo& proc_info = interaction->ProcInfo();
	const Target& tgt = init_state.Tgt();

	// If the event is not RES this Calc doesn't handle it
	if (!proc_info.IsResonant()) return 1.0;

	// WEIGHT CALCULATIONS -------------
	double w = 1.0;

	// CCRES Dial
	if (proc_info.IsWeakCC()) w *= fCur_NormCCRES;

	// Return Combined Weight
	return w;
	}

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

	void MINERvAReWeight_RES::SetDialValue(int rwenum, double val) {
	// Check Handled
	int curenum = rwenum % 1000;
	if (!IsHandled(curenum)) return;

	// Set Values
	if (curenum == kMINERvARW_NormCCRES) {
	fTwk_NormCCRES = val;
	fCur_NormCCRES = fDef_NormCCRES + fTwk_NormCCRES;
	}

	// Define Tweaked
	fTweaked = ((fTwk_NormCCRES != 0.0));
	}

	bool MINERvAReWeight_RES::IsHandled(int rwenum) {
	int curenum = rwenum % 1000;

	switch (curenum) {
	case kMINERvARW_NormCCRES:
	return true;
	default:
	return false;
	}
	}

	//*******************************************************
	RikRPA::RikRPA() {
	//*******************************************************

	// - Syst : kMINERvA_RikRPA_ApplyRPA
	// - Type : Binary
	// - Limits : 0.0 (false) -> 1.0 (true)
	// - Default : 0.0
	fApplyDial_RPACorrection = false;

	// - Syst : kMINERvA_RikRPA_LowQ2
	// - Type : Absolute
	// - Limits : 1.0 -> 1.0
	// - Default : 0.0
	// - Frac Error : 100%
	fDefDial_RPALowQ2 = 0.0;
	fCurDial_RPALowQ2 = fDefDial_RPALowQ2;
	fErrDial_RPALowQ2 = 0.0;

	// - Syst : kMINERvA_RikRPA_HighQ2
	// - Type : Absolute
	// - Limits : 1.0 -> 1.0
	// - Default : 0.0
	// - Frac Error : 100%
	fDefDial_RPAHighQ2 = 0.0;
	fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2;
	fErrDial_RPAHighQ2 = 1.0;

	fEventWeights = new double[5];

	for (size_t i = 0; i < kMaxCalculators; i++) {
	- // fRPACalculators[i] = NULL;
	+ fRPACalculators[i] = NULL;
	}

	fTweaked = false;
	}

	RikRPA::~RikRPA() {
	- delete fEventWeights;
	+ // delete fEventWeights;

	- for (size_t i = 0; i < kMaxCalculators; i++) {
	- // if (fRPACalculators[i]) delete fRPACalculators[i];
	- // fRPACalculators[i] = NULL;
	- }
	+ // for (size_t i = 0; i < kMaxCalculators; i++) {
	+ // if (fRPACalculators[i]) delete fRPACalculators[i];
	+ // fRPACalculators[i] = NULL;
	+ // }
	}

	double RikRPA::CalcWeight(BaseFitEvt* evt) {
	// LOG(FIT) << "Calculating RikRPA" << std::endl;
	// Return 1.0 if not tweaked
	if (!fTweaked) return 1.0;

	double w = 1.0;

	// Extract the GENIE Record
	GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
	const Interaction* interaction = ghep->Summary();
	const InitialState& init_state = interaction->InitState();
	const ProcessInfo& proc_info = interaction->ProcInfo();
	// const Kinematics & kine = interaction->Kine();
	// const XclsTag & xcls = interaction->ExclTag();
	const Target& tgt = init_state.Tgt();

	// If not QE return 1.0
	// LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() <<
	// std::endl;
	if (!tgt.IsNucleus()) return 1.0;
	if (!proc_info.IsQuasiElastic()) return 1.0;

	// Extract Beam and Target PDG
	GHepParticle* neutrino = ghep->Probe();
	int bpdg = neutrino->Pdg();

	GHepParticle* target = ghep->Particle(1);
	assert(target);
	int tpdg = target->Pdg();

	// Find the enum we need
	int calcenum = GetRPACalcEnum(bpdg, tpdg);
	if (calcenum == -1) return 1.0;

	// Check we have the RPA Calc setup for this enum
	// if not, set it up at that point
	- // if (!fRPACalculators[calcenum]) SetupRPACalculator(calcenum);
	- // weightRPA* rpacalc = fRPACalculators[calcenum];
	- // if (!rpacalc) {
	- // THROW("Failed to grab the RPA Calculator : " << calcenum);
	- // }
	+ if (!fRPACalculators[calcenum]) SetupRPACalculator(calcenum);
	+ weightRPA* rpacalc = fRPACalculators[calcenum];
	+ if (!rpacalc) {
	+ THROW("Failed to grab the RPA Calculator : " << calcenum);
	+ }

	// Extract Q0-Q3
	GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
	const TLorentzVector& k1 = *(neutrino->P4());
	const TLorentzVector& k2 = *(fsl->P4());
	double q0 = fabs((k1 - k2).E());
	double q3 = fabs((k1 - k2).Vect().Mag());
	+ double Q2 = fabs((k1 - k2).Mag2());

	// Now use q0-q3 and RPA Calculator to fill fWeights
	- // LOG(FIT) << "Getting Weights = " << q0 << " " << q3 << std::endl;
	- // rpacalc->getWeight(q0, q3, fEventWeights);
	+ //LOG(FIT) << "Getting Weights = " << q0 << " " << q3 << std::endl;
	+ rpacalc->getWeight(q0, q3, fEventWeights);

	// Apply Interpolation (for the time being simple linear)

	// Syst Application : kMINERvA_RikRPA_ApplyRPA
	if (fApplyDial_RPACorrection) {
	w *= fEventWeights[0]; // CV
	}

	+ /*
	LOG(FIT) << " fCurDial_RPALowQ2 = " << fCurDial_RPALowQ2
	<< " fCurDial_RPAHighQ2 = " << fCurDial_RPAHighQ2 << " Weights "
	<< fEventWeights[0] << " " << fEventWeights[1] << " "
	<< fEventWeights[2] << " " << fEventWeights[3] << " "
	<< fEventWeights[4] << std::endl;
	-
	+ */
	// Syst Application : kMINERvA_RikRPA_LowQ2
	if (fabs(fCurDial_RPALowQ2) > 0.0) {
	double interpw = fEventWeights[0];

	- if (fCurDial_RPALowQ2 > 0.0) {
	+ if (fCurDial_RPALowQ2 > 0.0 && Q2 < 2.0) {
	interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[1]) *
	fCurDial_RPALowQ2; // WLow+ } else if
	- } else if (fCurDial_RPALowQ2 < 0.0) {
	+ } else if (fCurDial_RPALowQ2 < 0.0 && Q2 < 2.0) {
	interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) *
	fCurDial_RPALowQ2; // WLow-
	}
	w *= interpw / fEventWeights[0]; // Div by CV again
	}

	// Syst Application : kMINERvA_RikRPA_HighQ2
	if (fabs(fCurDial_RPAHighQ2) > 0.0) {
	double interpw = fEventWeights[0];
	+
	if (fCurDial_RPAHighQ2 > 0.0) {
	interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
	- fCurDial_RPAHighQ2; // WHigh+ } else
	- if (fCurDial_RPAHighQ2 < 0.0) {
	+ fCurDial_RPAHighQ2; // WHigh+
	+
	+ } else if (fCurDial_RPAHighQ2 < 0.0) {
	interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
	fCurDial_RPAHighQ2; // WHigh-
	- }
	- w *= interpw / fEventWeights[0]; // Div by CV again
	}
	+ w *= interpw / fEventWeights[0]; // Div by CV again
	}

	// LOG(FIT) << "RPA Weight = " << w << std::endl;
	return w;
	} // namespace reweight

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

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

	// Check Handled
	if (!IsHandled(curenum)) return;
	if (curenum == kMINERvARW_RikRPA_ApplyRPA)
	fApplyDial_RPACorrection = (val > 0.5);
	if (curenum == kMINERvARW_RikRPA_LowQ2) fCurDial_RPALowQ2 = val;
	if (curenum == kMINERvARW_RikRPA_HighQ2) fCurDial_RPAHighQ2 = val;

	// Assign flag to say stuff has changed
	fTweaked = (fApplyDial_RPACorrection \|\|
	fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 \|\|
	fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0);
	}

	bool RikRPA::IsHandled(int rwenum) {
	int curenum = rwenum % 1000;
	switch (curenum) {
	case kMINERvARW_RikRPA_ApplyRPA:
	return true;
	case kMINERvARW_RikRPA_LowQ2:
	return true;
	case kMINERvARW_RikRPA_HighQ2:
	return true;
	default:
	return false;
	}
	}

	void RikRPA::SetupRPACalculator(int calcenum) {
	std::string rwdir = FitPar::GetDataBase() + "/reweight/MINERvA/RikRPA/";
	std::string fidir = "";
	switch (calcenum) {
	case kNuMuC12:
	fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root";
	break;

	case kNuMuO16:
	fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root";
	break;

	case kNuMuAr40:
	fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root";
	break;

	case kNuMuCa40:
	fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root";
	break;

	case kNuMuFe56:
	fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root";
	break;

	case kNuMuBarC12:
	fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root";
	break;

	case kNuMuBarO16:
	fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root";
	break;

	case kNuMuBarAr40:
	fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root";
	break;

	case kNuMuBarCa40:
	fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root";
	break;

	case kNuMuBarFe56:
	fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root";
	break;
	}

	LOG(FIT) << "Loading RPA CALC : " << fidir << std::endl;
	TDirectory* olddir = gDirectory;

	- // fRPACalculators[calcenum] = new weightRPA(rwdir + "/" + fidir);
	+ std::cout << "***********************************************" << std::endl;
	+ std::cout << "Loading a new weightRPA calculator" << std::endl;
	+ std::cout << "Authors: Rik Gran, Heidi Schellman" << std::endl;
	+ std::cout << "Citation: arXiv:1705.02932 [hep-ex]" << std::endl;
	+ std::cout << "***********************************************" << std::endl;
	+
	+ fRPACalculators[calcenum] = new weightRPA(rwdir + "/" + fidir);
	olddir->cd();
	return;
	}

	int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) {
	if (bpdg == 14 && tpdg == 1000060120)
	return kNuMuC12;
	else if (bpdg == 14 && tpdg == 1000080160)
	return kNuMuO16;
	else if (bpdg == 14 && tpdg == 1000180400)
	return kNuMuAr40;
	else if (bpdg == 14 && tpdg == 1000200400)
	return kNuMuCa40;
	else if (bpdg == 14 && tpdg == 1000280560)
	return kNuMuFe56;
	else if (bpdg == -14 && tpdg == 1000060120)
	return kNuMuBarC12;
	else if (bpdg == -14 && tpdg == 1000080160)
	return kNuMuBarO16;
	else if (bpdg == -14 && tpdg == 1000180400)
	return kNuMuBarAr40;
	else if (bpdg == -14 && tpdg == 1000200400)
	return kNuMuBarCa40;
	else if (bpdg == -14 && tpdg == 1000280560)
	return kNuMuBarFe56;
	else {
	ERROR(WRN, "Unknown beam and target combination for RPA Calcs! "
	<< bpdg << " " << tpdg);
	}

	return -1;
	}

	} // namespace reweight
	} // namespace nuisance
	#endif
	#endif
	diff --git a/src/Reweight/MINERvAWeightCalcs.h b/src/Reweight/MINERvAWeightCalcs.h
	index d3ec9a1..5bedfd3 100644
	--- a/src/Reweight/MINERvAWeightCalcs.h
	+++ b/src/Reweight/MINERvAWeightCalcs.h
	@@ -1,136 +1,136 @@
	#ifndef MINERVA_WEIGHT_CALCS
	#define MINERVA_WEIGHT_CALCS

	#include <string>

	#ifdef __MINERVA_RW_ENABLED__
	#ifdef __GENIE_ENABLED__
	#include "Conventions/Units.h"
	#include "EVGCore/EventRecord.h"
	#include "FitEvent.h"
	#include "FitParameters.h"
	#include "GHEP/GHepParticle.h"
	#include "GHEP/GHepRecord.h"
	#include "GHEP/GHepUtils.h"
	#include "GeneralUtils.h"
	#include "NUISANCESyst.h"
	#include "NUISANCEWeightCalcs.h"
	#include "Ntuple/NtpMCEventRecord.h"
	#include "PDG/PDGUtils.h"
	#include "WeightUtils.h"
	-// #include "weightRPA.h"
	+#include "weightRPA.h"
	using namespace genie;
	class BaseFitEvt;


	namespace nuisance {
	namespace reweight {

	// MEC Dials
	class MINERvAReWeight_QE : public NUISANCEWeightCalc {
	public:
	MINERvAReWeight_QE();
	virtual ~MINERvAReWeight_QE();

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

	double fTwk_NormCCQE;
	double fCur_NormCCQE;
	double fDef_NormCCQE;
	bool fTweaked;

	};


	// MEC Dials
	class MINERvAReWeight_MEC : public NUISANCEWeightCalc {
	public:
	MINERvAReWeight_MEC();
	virtual ~MINERvAReWeight_MEC();

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

	double fTwk_NormCCMEC;
	double fCur_NormCCMEC;
	double fDef_NormCCMEC;
	bool fTweaked;

	};

	// RES Dials
	class MINERvAReWeight_RES : public NUISANCEWeightCalc {
	public:
	MINERvAReWeight_RES();
	virtual ~MINERvAReWeight_RES();

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

	double fTwk_NormCCRES;
	double fCur_NormCCRES;
	double fDef_NormCCRES;
	bool fTweaked;

	};

	/// RPA Weight Calculator that applies RPA systematics
	/// to GENIE events. GENIE EVENTS ONLY!
	class RikRPA : public NUISANCEWeightCalc {
	public:
	RikRPA();
	~RikRPA();

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

	void SetupRPACalculator(int calcenum);
	int GetRPACalcEnum(int bpdg, int tpdg);

	bool fApplyDial_RPACorrection;

	double fTwkDial_RPALowQ2;
	double fDefDial_RPALowQ2;
	double fCurDial_RPALowQ2;
	double fErrDial_RPALowQ2;

	double fTwkDial_RPAHighQ2;
	double fDefDial_RPAHighQ2;
	double fCurDial_RPAHighQ2;
	double fErrDial_RPAHighQ2;

	double* fEventWeights;
	bool fTweaked;

	const static int kMaxCalculators = 10;
	enum rpacalcenums {
	kNuMuC12,
	kNuMuO16,
	kNuMuAr40,
	kNuMuCa40,
	kNuMuFe56,
	kNuMuBarC12,
	kNuMuBarO16,
	kNuMuBarAr40,
	kNuMuBarCa40,
	kNuMuBarFe56
	};
	- // weightRPA* fRPACalculators[kMaxCalculators];
	+ weightRPA* fRPACalculators[kMaxCalculators];
	};

	}; // namespace reweight
	}; // namespace nuisance

	#endif // __GENIE_ENABLED__
	#endif //__MINERVA_RW_ENABLED__

	#endif

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 5:37 PM (1 d, 13 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805424
Default Alt Text: (47 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions