No OneTemporary
Actions

Size

81 KB

Subscribers

None

View Options

	diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
	index c260ba6..878957e 100644
	--- a/app/CMakeLists.txt
	+++ b/app/CMakeLists.txt
	@@ -1,192 +1,201 @@
	# 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(TARGETS_TO_BUILD)

	if(USE_MINIMIZER)
	add_executable(nuismin nuismin.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuismin)
	target_link_libraries(nuismin ${MODULETargets})
	target_link_libraries(nuismin ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(nuismin ${ROOT_LIBS})
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(nuismin PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()

	add_executable(nuissplines nuissplines.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissplines)
	target_link_libraries(nuissplines ${MODULETargets})
	target_link_libraries(nuissplines ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(nuissplines ${ROOT_LIBS})
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(nuissplines PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()

	endif()

	include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
	include_directories(${CMAKE_SOURCE_DIR}/src/Routines)
	include_directories(${CMAKE_SOURCE_DIR}/src/InputHandler)
	include_directories(${CMAKE_SOURCE_DIR}/src/Genie)
	include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
	include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
	include_directories(${CMAKE_SOURCE_DIR}/src/Splines)
	include_directories(${CMAKE_SOURCE_DIR}/src/Reweight)
	include_directories(${CMAKE_SOURCE_DIR}/src/FCN)
	include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies)
	include_directories(${CMAKE_SOURCE_DIR}/src/Smearceptance)
	include_directories(${EXP_INCLUDE_DIRECTORIES})

	if (USE_NuWro AND NOT NUWRO_BUILT_FROM_FILE)
	add_executable(nuwro_nuisance nuwro_NUISANCE.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuwro_nuisance)
	target_link_libraries(nuwro_nuisance ${MODULETargets})
	target_link_libraries(nuwro_nuisance ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(nuwro_nuisance ${ROOT_LIBS})
	include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(nuwro_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()
	endif()

	if (USE_NEUT)
	add_executable(neut_nuisance neut_NUISANCE.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};neut_nuisance)
	target_link_libraries(neut_nuisance ${MODULETargets})
	target_link_libraries(neut_nuisance ${CMAKE_DEPENDLIB_FLAGS})
	target_link_libraries(neut_nuisance ${ROOT_LIBS})
	include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(neut_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()
	endif()

	if (BUILD_GEVGEN)
	add_executable(gevgen_nuisance gEvGen_NUISANCE.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance)
	target_link_libraries(gevgen_nuisance ${MODULETargets})
	target_link_libraries(gevgen_nuisance ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(gevgen_nuisance ${ROOT_LIBS})
	include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
	include_directories(${GENIE_INCLUDES}/Apps)
	include_directories(${GENIE_INCLUDES}/FluxDrivers)
	include_directories(${GENIE_INCLUDES}/EVGDrivers)
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(gevgen_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()


	add_executable(gevgen_nuisance_mixed gEvGen_NUISANCE_MIXED.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance_mixed)
	target_link_libraries(gevgen_nuisance_mixed ${MODULETargets})
	target_link_libraries(gevgen_nuisance_mixed ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(gevgen_nuisance_mixed ${ROOT_LIBS})
	include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
	include_directories(${GENIE_INCLUDES}/Apps)
	include_directories(${GENIE_INCLUDES}/FluxDrivers)
	include_directories(${GENIE_INCLUDES}/EVGDrivers)
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(gevgen_nuisance_mixed PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()
	endif()

	if (USE_GiBUU)
	add_executable(DumpGiBUUEvents DumpGiBUUEvents.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};DumpGiBUUEvents)
	target_link_libraries(DumpGiBUUEvents ${MODULETargets})
	target_link_libraries(DumpGiBUUEvents ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(DumpGiBUUEvents ${ROOT_LIBS})
	include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(DumpGiBUUEvents PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()
	endif()

	add_executable(nuiscomp nuiscomp.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuiscomp)
	target_link_libraries(nuiscomp ${MODULETargets})
	target_link_libraries(nuiscomp ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(nuiscomp ${ROOT_LIBS})
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(nuiscomp PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()

	add_executable(nuisflat nuisflat.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisflat)
	target_link_libraries(nuisflat ${MODULETargets})
	target_link_libraries(nuisflat ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(nuisflat ${ROOT_LIBS})
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(nuisflat PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()

	add_executable(nuissmear nuissmear.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissmear)
	target_link_libraries(nuissmear ${MODULETargets})
	target_link_libraries(nuissmear ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(nuissmear ${ROOT_LIBS})
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(nuissmear PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()

	add_executable(nuissyst nuissyst.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissyst)
	target_link_libraries(nuissyst ${MODULETargets})
	target_link_libraries(nuissyst ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(nuissyst ${ROOT_LIBS})
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(nuissyst PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()

	+add_executable(nuisbayes nuisbayes.cxx)
	+set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisbayes)
	+target_link_libraries(nuisbayes ${MODULETargets})
	+target_link_libraries(nuisbayes ${CMAKE_DEPENDLIB_FLAGS})
	+# target_link_libraries(nuisbayes ${ROOT_LIBS})
	+if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	+ set_target_properties(nuisbayes PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	+endif()
	+
	if(USE_GENIE)
	add_executable(PrepareGENIE PrepareGENIE.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareGENIE)
	target_link_libraries(PrepareGENIE ${MODULETargets})
	target_link_libraries(PrepareGENIE ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(PrepareGENIE ${ROOT_LIBS})
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(PrepareGENIE PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()
	endif()

	if(USE_NEUT)
	add_executable(PrepareNEUT PrepareNEUT.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNEUT)
	target_link_libraries(PrepareNEUT ${MODULETargets})
	target_link_libraries(PrepareNEUT ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(PrepareNEUT ${ROOT_LIBS})
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(PrepareNEUT PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()
	endif()

	# PREPARE NUWRO
	# Commented out for the time being until it is finished..
	if(USE_NuWro)
	add_executable(PrepareNuwro PrepareNuwroEvents.cxx)
	set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNuwro)
	target_link_libraries(PrepareNuwro ${MODULETargets})
	target_link_libraries(PrepareNuwro ${CMAKE_DEPENDLIB_FLAGS})
	# target_link_libraries(PrepareNuwro ${ROOT_LIBS})
	if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
	set_target_properties(PrepareNuwro PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
	endif()
	endif()

	install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin)
	diff --git a/app/nuisbayes.cxx b/app/nuisbayes.cxx
	new file mode 100644
	index 0000000..44a3df0
	--- /dev/null
	+++ b/app/nuisbayes.cxx
	@@ -0,0 +1,68 @@
	+// 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/>.
	+*******************************************************************************/
	+
	+// Author: Patrick Stowell 09/2016
	+/*
	+ Usage: ./nuissyst -c card file -o output file, where the results of the throws are stored
	+ where:
	+*/
	+
	+#include "BayesianRoutines.h"
	+
	+//*******************************
	+void printInputCommands(){
	+//*******************************
	+
	+ exit(-1);
	+
	+};
	+
	+//*******************************
	+int main(int argc, char* argv[]){
	+//*******************************
	+
	+ // Program status;
	+ int status = 0;
	+
	+ // If No Arguments print commands
	+ if (argc == 1) printInputCommands();
	+
	+ for (int i = 1; i< argc; ++i){
	+ // Cardfile
	+ if (!std::strcmp(argv[i], "-h")) printInputCommands();
	+ else break;
	+ }
	+
	+ // Read input arguments such as card file, parameter arguments, and fit routines
	+ LOG(FIT)<<"Starting nuissyst"<<std::endl;
	+
	+ // Make systematic class and run fit
	+ BayesianRoutines* min = new BayesianRoutines(argc, argv);
	+ min->Run();
	+
	+
	+ // Show Final Status
	+ LOG(FIT)<<"-------------------------------------"<<std::endl;
	+ LOG(FIT)<<"FINISHED" << std::endl;
	+ LOG(FIT)<<"-------------------------------------"<<std::endl;
	+ return status;
	+}
	+
	+
	+
	diff --git a/parameters/config.xml b/parameters/config.xml
	index f10a590..4e103d0 100644
	--- a/parameters/config.xml
	+++ b/parameters/config.xml
	@@ -1,213 +1,214 @@
	<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" />
	-->

	<config UseShapeCovar="0" />
	+<config CC0piNBINS="156" />

	</nuisance>

	diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx
	index 753ac96..5102cbd 100755
	--- a/src/FCN/JointFCN.cxx
	+++ b/src/FCN/JointFCN.cxx
	@@ -1,1014 +1,1136 @@
	#include "JointFCN.h"
	#include <stdio.h>
	#include "FitUtils.h"


	//***************************************************
	JointFCN::JointFCN(TFile* outfile) {
	//***************************************************

	fOutputDir = gDirectory;
	if (outfile) FitPar::Config().out = outfile;

	std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
	LoadSamples(samplekeys);

	std::vector<nuiskey> covarkeys = Config::QueryKeys("covar");
	LoadPulls(covarkeys);

	fCurIter = 0;
	fMCFilled = false;

	fIterationTree = false;
	fDialVals = NULL;
	fNDials = 0;

	fUsingEventManager = FitPar::Config().GetParB("EventManager");
	fOutputDir->cd();
	}

	//***************************************************
	JointFCN::JointFCN(std::vector<nuiskey> samplekeys, TFile* outfile) {
	//***************************************************

	fOutputDir = gDirectory;
	if (outfile) FitPar::Config().out = outfile;

	LoadSamples(samplekeys);

	fCurIter = 0;
	fMCFilled = false;

	fOutputDir->cd();

	fIterationTree = false;
	fDialVals = NULL;
	fNDials = 0;

	fUsingEventManager = FitPar::Config().GetParB("EventManager");
	fOutputDir->cd();
	}

	//***************************************************
	JointFCN::~JointFCN() {
	//***************************************************

	// Delete Samples
	for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
	iter++) {
	MeasurementBase* exp = *iter;
	delete exp;
	}

	for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
	ParamPull* pull = *iter;
	delete pull;
	}

	// Sort Tree
	if (fIterationTree) DestroyIterationTree();
	if (fDialVals) delete fDialVals;
	if (fSampleLikes) delete fSampleLikes;
	};

	//***************************************************
	void JointFCN::CreateIterationTree(std::string name, FitWeight* rw) {
	//***************************************************

	LOG(FIT) << " Creating new iteration container! " << std::endl;
	DestroyIterationTree();
	fIterationTreeName = name;

	// Add sample likelihoods and ndof
	for (MeasListConstIter iter = fSamples.begin();
	iter != fSamples.end();
	iter++) {

	MeasurementBase* exp = *iter;
	std::string name = exp->GetName();

	std::string liketag = name + "_likelihood";
	fNameValues.push_back(liketag);
	fCurrentValues.push_back(0.0);

	std::string ndoftag = name + "_ndof";
	fNameValues.push_back(ndoftag);
	fCurrentValues.push_back(0.0);
	}

	// Add Pull terms
	for (PullListConstIter iter = fPulls.begin();
	iter != fPulls.end(); iter++) {

	ParamPull* pull = *iter;
	std::string name = pull->GetName();

	std::string liketag = name + "_likelihood";
	fNameValues.push_back(liketag);
	fCurrentValues.push_back(0.0);

	std::string ndoftag = name + "_ndof";
	fNameValues.push_back(ndoftag);
	fCurrentValues.push_back(0.0);
	}

	// Add Likelihoods
	fNameValues.push_back("total_likelihood");
	fCurrentValues.push_back(0.0);

	fNameValues.push_back("total_ndof");
	fCurrentValues.push_back(0.0);

	- // Setup Containers
	+ // Setup Containers
	fSampleN = fSamples.size() + fPulls.size();
	fSampleLikes = new double[fSampleN];
	fSampleNDOF = new int[fSampleN];

	// Add Dials
	std::vector<std::string> dials = rw->GetDialNames();
	- for (size_t i = 0; i < dials.size(); i++){
	+ for (size_t i = 0; i < dials.size(); i++) {
	fNameValues.push_back( dials[i] );
	fCurrentValues.push_back( 0.0 );
	}
	fNDials = dials.size();
	- fDialVals = new double[fNDials];
	+ fDialVals = new double[fNDials];

	// Set IterationTree Flag
	fIterationTree = true;
	-
	+
	}

	//***************************************************
	void JointFCN::DestroyIterationTree() {
	//***************************************************

	fIterationCount.clear();
	fCurrentValues.clear();
	fNameValues.clear();
	fIterationValues.clear();

	}

	//***************************************************
	void JointFCN::WriteIterationTree() {
	//***************************************************
	LOG(FIT) << "Writing iteration tree" << std::endl;

	// Make a new TTree
	TTree* itree = new TTree(fIterationTreeName.c_str(),
	fIterationTreeName.c_str());

	double* vals = new double[fNameValues.size()];
	int count = 0;

	- itree->Branch("iteration",&count,"Iteration/I");
	+ itree->Branch("iteration", &count, "Iteration/I");
	for (int i = 0; i < fNameValues.size(); i++) {
	itree->Branch( fNameValues[i].c_str(),
	&vals[i],
	(fNameValues[i] + "/D").c_str() );
	}

	// Fill Iterations
	- for (size_t i = 0; i < fIterationValues.size(); i++){
	+ for (size_t i = 0; i < fIterationValues.size(); i++) {
	std::vector<double> itervals = fIterationValues[i];

	// Fill iteration state
	count = fIterationCount[i];
	- for (size_t j = 0; j < itervals.size(); j++){
	+ for (size_t j = 0; j < itervals.size(); j++) {
	vals[j] = itervals[j];
	}

	// Save to TTree
	itree->Fill();
	}

	// Write to file
	itree->Write();
	}

	//***************************************************
	void JointFCN::FillIterationTree(FitWeight* rw) {
	//***************************************************

	// Loop over samples count
	int count = 0;
	- for (int i = 0; i < fSampleN; i++){
	+ for (int i = 0; i < fSampleN; i++) {
	fCurrentValues[count++] = fSampleLikes[i];
	fCurrentValues[count++] = double(fSampleNDOF[i]);
	}

	// Fill Totals
	fCurrentValues[count++] = fLikelihood;
	fCurrentValues[count++] = double(fNDOF);

	// Loop Over Parameter Counts
	rw->GetAllDials(fDialVals, fNDials);
	- for (int i = 0; i < fNDials; i++){
	+ for (int i = 0; i < fNDials; i++) {
	fCurrentValues[count++] = double(fDialVals[i]);
	}

	// Push Back Into Container
	fIterationCount.push_back( fCurIter );
	fIterationValues.push_back(fCurrentValues);
	}

	//***************************************************
	double JointFCN::DoEval(const double* x) {
	//***************************************************

	// WEIGHT ENGINE
	fDialChanged = FitBase::GetRW()->HasRWDialChanged(x);
	FitBase::GetRW()->UpdateWeightEngine(x);
	if (fDialChanged) {
	FitBase::GetRW()->Reconfigure();
	FitBase::EvtManager().ResetWeightFlags();
	}
	if (LOG_LEVEL(REC)) {
	FitBase::GetRW()->Print();
	}

	// SORT SAMPLES
	ReconfigureSamples();

	// GET TEST STAT
	fLikelihood = GetLikelihood();
	fNDOF = GetNDOF();
	-
	+
	// PRINT PROGRESS
	LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood
	<< std::endl;

	// UPDATE TREE
	if (fIterationTree) FillIterationTree(FitBase::GetRW());

	return fLikelihood;
	}

	//***************************************************
	int JointFCN::GetNDOF() {
	//***************************************************

	int totaldof = 0;
	int count = 0;

	// Total number of Free bins in each MC prediction
	for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
	iter++) {
	MeasurementBase* exp = *iter;
	int dof = exp->GetNDOF();

	// Save Seperate DOF
	if (fIterationTree) {
	fSampleNDOF[count] = dof;
	}

	// Add to total
	totaldof += dof;
	count++;
	}

	// Loop over pulls
	for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
	ParamPull* pull = *iter;
	double dof = pull->GetLikelihood();

	// Save seperate DOF
	if (fIterationTree) {
	fSampleNDOF[count] = dof;
	}

	// Add to total
	totaldof += dof;
	count++;
	}

	// Set Data Variable
	- if (fIterationTree){
	- fSampleNDOF[count] = totaldof;
	+ if (fIterationTree) {
	+ fSampleNDOF[count] = totaldof;
	}
	return totaldof;
	}

	//***************************************************
	double JointFCN::GetLikelihood() {
	//***************************************************

	LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..."
	<< " : "
	<< "-2logL" << std::endl;

	// Loop and add up likelihoods in an uncorrelated way
	double like = 0.0;
	int count = 0;
	for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
	iter++) {
	MeasurementBase* exp = *iter;
	double newlike = exp->GetLikelihood();
	int ndof = exp->GetNDOF();
	// Save seperate likelihoods
	if (fIterationTree) {
	fSampleLikes[count] = newlike;
	}

	LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : "
	<< newlike << "/" << ndof << std::endl;

	// Add Weight Scaling
	// like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName());

	// Add to total
	like += newlike;
	count++;
	}

	// Loop over pulls
	for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
	ParamPull* pull = *iter;
	double newlike = pull->GetLikelihood();

	// Save seperate likelihoods
	if (fIterationTree) {
	fSampleLikes[count] = newlike;
	}

	// Add to total
	like += newlike;
	count++;
	}

	// Set Data Variable
	fLikelihood = like;
	- if (fIterationTree){
	+ if (fIterationTree) {
	fSampleLikes[count] = fLikelihood;
	}

	return like;
	};

	void JointFCN::LoadSamples(std::vector<nuiskey> samplekeys) {
	LOG(MIN) << "Loading Samples : " << samplekeys.size() << std::endl;
	for (size_t i = 0; i < samplekeys.size(); i++) {
	nuiskey key = samplekeys[i];

	// Get Sample Options
	std::string samplename = key.GetS("name");
	std::string samplefile = key.GetS("input");
	std::string sampletype = key.GetS("type");
	std::string fakeData = "";

	LOG(MIN) << "Loading Sample : " << samplename << std::endl;

	fOutputDir->cd();
	MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(key);

	if (!NewLoadedSample) {
	ERR(FTL) << "Could not load sample provided: " << samplename << std::endl;
	ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
	<< std::endl;
	throw;
	} else {
	fSamples.push_back(NewLoadedSample);
	}
	}
	}

	//***************************************************
	void JointFCN::LoadPulls(std::vector<nuiskey> pullkeys) {
	//***************************************************
	for (size_t i = 0; i < pullkeys.size(); i++) {
	nuiskey key = pullkeys[i];

	std::string pullname = key.GetS("name");
	std::string pullfile = key.GetS("input");
	std::string pulltype = key.GetS("type");

	fOutputDir->cd();
	fPulls.push_back(new ParamPull(pullname, pullfile, pulltype));
	}
	}

	//***************************************************
	void JointFCN::ReconfigureSamples(bool fullconfig) {
	//***************************************************

	int starttime = time(NULL);
	LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << std::endl;
	// std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled <<
	// std::endl;
	// Event Manager Reconf
	if (fUsingEventManager) {
	if (!fullconfig and fMCFilled)
	ReconfigureFastUsingManager();
	else
	ReconfigureUsingManager();

	} else {
	// Loop over all Measurement Classes
	for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
	iter++) {
	MeasurementBase* exp = *iter;

	// If RW Either do signal or full reconfigure.
	if (fDialChanged or !fMCFilled or fullconfig) {
	if (!fullconfig and fMCFilled)
	exp->ReconfigureFast();
	else
	exp->Reconfigure();

	// If RW Not needed just do normalisation
	} else {
	exp->Renormalise();
	}
	}
	}

	// Loop over pulls and update
	for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
	ParamPull* pull = *iter;
	pull->Reconfigure();
	}

	fMCFilled = true;
	LOG(MIN) << "Finished Reconfigure iter. " << fCurIter << " in "
	<< time(NULL) - starttime << "s" << std::endl;

	fCurIter++;
	}

	//***************************************************
	void JointFCN::ReconfigureSignal() {
	//***************************************************
	ReconfigureSamples(false);
	}

	//***************************************************
	void JointFCN::ReconfigureAllEvents() {
	//***************************************************
	FitBase::GetRW()->Reconfigure();
	FitBase::EvtManager().ResetWeightFlags();
	ReconfigureSamples(true);
	}

	std::vector<InputHandlerBase*> JointFCN::GetInputList() {
	std::vector<InputHandlerBase*> InputList;
	fIsAllSplines = true;

	MeasListConstIter iterSam = fSamples.begin();
	for (; iterSam != fSamples.end(); iterSam++) {
	MeasurementBase* exp = (*iterSam);

	std::vector<MeasurementBase*> subsamples = exp->GetSubSamples();
	for (size_t i = 0; i < subsamples.size(); i++) {
	InputHandlerBase* inp = subsamples[i]->GetInput();
	if (std::find(InputList.begin(), InputList.end(), inp) ==
	InputList.end()) {
	if (subsamples[i]->GetInput()->GetType() != kSPLINEPARAMETER)
	fIsAllSplines = false;

	InputList.push_back(subsamples[i]->GetInput());
	}
	}
	}

	return InputList;
	}

	std::vector<MeasurementBase*> JointFCN::GetSubSampleList() {
	std::vector<MeasurementBase*> SampleList;

	MeasListConstIter iterSam = fSamples.begin();
	for (; iterSam != fSamples.end(); iterSam++) {
	MeasurementBase* exp = (*iterSam);

	std::vector<MeasurementBase*> subsamples = exp->GetSubSamples();
	for (size_t i = 0; i < subsamples.size(); i++) {
	SampleList.push_back(subsamples[i]);
	}
	}

	return SampleList;
	}

	//***************************************************
	void JointFCN::ReconfigureUsingManager() {
	//***************************************************

	// 'Slow' Event Manager Reconfigure
	LOG(REC) << "Event Manager Reconfigure" << std::endl;
	int timestart = time(NULL);

	// Reset all samples
	MeasListConstIter iterSam = fSamples.begin();
	for (; iterSam != fSamples.end(); iterSam++) {
	MeasurementBase* exp = (*iterSam);
	exp->ResetAll();
	}

	// If we are siving signal, reset all containers.
	bool savesignal = (FitPar::Config().GetParB("SignalReconfigures"));

	if (savesignal) {
	// Reset all of our event signal vectors
	fSignalEventBoxes.clear();
	fSignalEventFlags.clear();
	fSampleSignalFlags.clear();
	fSignalEventSplines.clear();
	}

	// Make sure we have a list of inputs
	if (fInputList.empty()) {
	fInputList = GetInputList();
	fSubSampleList = GetSubSampleList();
	}

	// If all inputs are splines make sure the readers are told
	// they need to be reconfigured.
	std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();

	if (fIsAllSplines) {
	for (; inp_iter != fInputList.end(); inp_iter++) {
	InputHandlerBase* curinput = (*inp_iter);

	// Tell reader in each BaseEvent it needs a Reconfigure next weight calc.
	BaseFitEvt* curevent = curinput->FirstBaseEvent();
	if (curevent->fSplineRead) {
	curevent->fSplineRead->SetNeedsReconfigure(true);
	}
	}
	}

	// MAIN INPUT LOOP ====================

	int fillcount = 0;
	int inputcount = 0;
	inp_iter = fInputList.begin();

	// Loop over each input in manager
	for (; inp_iter != fInputList.end(); inp_iter++) {
	InputHandlerBase* curinput = (*inp_iter);

	// Get event information
	FitEvent* curevent = curinput->FirstNuisanceEvent();
	curinput->CreateCache();

	int i = 0;
	int nevents = curinput->GetNEvents();
	int countwidth = nevents / 5;

	// Start event loop iterating until we get a NULL pointer.
	while (curevent) {
	// Get Event Weight
	curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
	curevent->Weight = curevent->RWWeight * curevent->InputWeight;
	double rwweight = curevent->Weight;
	// std::cout << "RWWeight = " << curevent->RWWeight << " " <<
	// curevent->InputWeight << std::endl;

	// Logging
	// std::cout << CHECKLOG(1) << std::endl;
	if (LOGGING(REC)) {
	if (i % countwidth == 0) {
	QLOG(REC, curinput->GetName()
	<< " : Processed " << i << " events. [M, W] = ["
	<< curevent->Mode << ", " << rwweight << "]");
	}
	}

	// Setup flag for if signal found in at least one sample
	bool foundsignal = false;

	// Create a new signal bitset for this event
	std::vector<bool> signalbitset(fSubSampleList.size());

	// Create a new signal box vector for this event
	std::vector<MeasurementVariableBox*> signalboxes;

	// Start measurement iterator
	size_t measitercount = 0;
	std::vector<MeasurementBase*>::iterator meas_iter =
	fSubSampleList.begin();

	// Loop over all subsamples (sub in JointMeas)
	for (; meas_iter != fSubSampleList.end(); meas_iter++) {
	MeasurementBase* curmeas = (*meas_iter);

	// Compare input pointers, to current input, skip if not.
	// Pointer tells us if it matches without doing ID checks.
	if (curinput != curmeas->GetInput()) {
	if (savesignal) {
	// Set bit to 0 as definitely not signal
	signalbitset[measitercount] = 0;
	}

	// Count up what measurement we are on.
	measitercount++;

	// Skip sample as input not signal.
	continue;
	}

	// Fill events for matching inputs.
	MeasurementVariableBox* box = curmeas->FillVariableBox(curevent);

	bool signal = curmeas->isSignal(curevent);
	curmeas->SetSignal(signal);
	curmeas->FillHistograms(curevent->Weight);

	// If its Signal tally up fills
	if (signal) {
	fillcount++;
	}

	// If we are saving signal/splines fill the bitset
	if (savesignal) {
	signalbitset[measitercount] = signal;
	}

	// If signal save a clone of the event box for use later.
	if (savesignal and signal) {
	foundsignal = true;
	signalboxes.push_back(box->CloneSignalBox());
	}

	// Keep track of Measurement we are on.
	measitercount++;
	}

	// Once we've filled the measurements, if saving signal
	// push back if any sample flagged this event as signal
	if (savesignal) {
	fSignalEventFlags.push_back(foundsignal);
	}

	// Save the vector of signal boxes for this event
	if (savesignal and foundsignal) {
	fSignalEventBoxes.push_back(signalboxes);
	fSampleSignalFlags.push_back(signalbitset);
	}

	// If all inputs are splines we can save the spline coefficients
	// for fast in memory reconfigures later.
	if (fIsAllSplines and savesignal and foundsignal) {
	// Make temp vector to push back with
	std::vector<float> coeff;
	for (size_t l = 0; l < (UInt_t)curevent->fSplineRead->GetNPar(); l++) {
	coeff.push_back(curevent->fSplineCoeff[l]);
	}

	// Push back to signal event splines. Kept in sync with
	// fSignalEventBoxes size.
	// int splinecount = fSignalEventSplines.size();
	fSignalEventSplines.push_back(coeff);

	// if (splinecount % 1000 == 0) {
	// std::cout << "Pushed Back Coeff " << splinecount << " : ";
	// for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++)
	// {
	// std::cout << " " << fSignalEventSplines[splinecount][l];
	// }
	// std::cout << std::endl;
	// }
	}

	// Clean up vectors once done with this event
	signalboxes.clear();
	signalbitset.clear();

	// Iterate to the next event.
	curevent = curinput->NextNuisanceEvent();
	i++;
	}

	// curinput->RemoveCache();

	// Keep track of what input we are on.
	inputcount++;
	}

	// End of Event Loop ===============================

	// Now event loop is finished loop over all Measurements
	// Converting Binned events to XSec Distributions
	iterSam = fSamples.begin();
	for (; iterSam != fSamples.end(); iterSam++) {
	MeasurementBase* exp = (*iterSam);
	exp->ConvertEventRates();
	}

	// Print out statements on approximate memory usage for profiling.
	LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
	if (savesignal) {
	int mem =
	( // sizeof(fSignalEventBoxes) +
	// fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
	sizeof(MeasurementVariableBox1D) * fillcount) *
	1E-6;
	LOG(REC) << " -> Saved " << fillcount
	<< " signal boxes for faster access. (~" << mem << " MB)"
	<< std::endl;
	if (fIsAllSplines and !fSignalEventSplines.empty()) {
	int splmem = sizeof(float) * fSignalEventSplines.size() *
	fSignalEventSplines[0].size() * 1E-6;
	LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size()
	<< " spline sets into memory. (~" << splmem << " MB)"
	<< std::endl;
	}
	}

	LOG(REC) << "Time taken ReconfigureUsingManager() : "
	<< time(NULL) - timestart << std::endl;

	// Check SignalReconfigures works for all samples
	if (savesignal) {
	double likefull = GetLikelihood();
	ReconfigureFastUsingManager();
	double likefast = GetLikelihood();

	if (fabs(likefull - likefast) > 0.0001)
	{
	ERROR(FTL, "Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast);
	ERROR(FTL, "This means some samples you are using are not setup to use SignalReconfigures=1");
	ERROR(FTL, "Please turn OFF signal reconfigures.");
	throw;
	} else {
	LOG(FIT) << "Likelihoods for FULL and FAST match. Will use FAST next time." << std::endl;
	}
	}

	// End of reconfigure
	return;
	};

	//***************************************************
	void JointFCN::ReconfigureFastUsingManager() {
	//***************************************************

	LOG(FIT) << " -> Doing FAST using manager" << std::endl;
	// Get Start time for profilling
	int timestart = time(NULL);

	// Reset all samples
	MeasListConstIter iterSam = fSamples.begin();
	for (; iterSam != fSamples.end(); iterSam++) {
	MeasurementBase* exp = (*iterSam);
	exp->ResetAll();
	}

	// Check for saved variables if not do a full reconfigure.
	if (fSignalEventFlags.empty()) {
	ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl;
	ReconfigureUsingManager();
	return;
	}

	bool fFillNuisanceEvent =
	FitPar::Config().GetParB("FullEventOnSignalReconfigure");

	// Setup fast vector iterators.
	std::vector<bool>::iterator inpsig_iter = fSignalEventFlags.begin();
	std::vector<std::vector<MeasurementVariableBox*> >::iterator box_iter =
	fSignalEventBoxes.begin();
	std::vector<std::vector<float> >::iterator spline_iter =
	fSignalEventSplines.begin();
	std::vector<std::vector<bool> >::iterator samsig_iter =
	fSampleSignalFlags.begin();
	int splinecount = 0;

	// Setup stuff for logging
	int fillcount = 0;
	int nevents = fSignalEventFlags.size();
	int countwidth = nevents / 20;

	// If All Splines tell splines they need a reconfigure.
	std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
	if (fIsAllSplines) {
	LOG(REC) << "All Spline Inputs so using fast spline loop." << std::endl;
	for (; inp_iter != fInputList.end(); inp_iter++) {
	InputHandlerBase* curinput = (*inp_iter);

	// Tell each fSplineRead in BaseFitEvent to reconf next weight calc
	BaseFitEvt* curevent = curinput->FirstBaseEvent();
	if (curevent->fSplineRead)
	curevent->fSplineRead->SetNeedsReconfigure(true);
	}
	}

	// Loop over all possible spline inputs
	double* coreeventweights = new double[fSignalEventBoxes.size()];
	splinecount = 0;

	inp_iter = fInputList.begin();
	inpsig_iter = fSignalEventFlags.begin();
	spline_iter = fSignalEventSplines.begin();

	// Loop over all signal flags
	// For each valid signal flag add one to splinecount
	// Get Splines from that count and add to weight
	// Add splinecount
	int sigcount = 0;
	splinecount = 0;

	// #pragma omp parallel for shared(splinecount,sigcount)
	for (uint iinput = 0; iinput < fInputList.size(); iinput++) {
	InputHandlerBase* curinput = fInputList[iinput];
	BaseFitEvt* curevent = curinput->FirstBaseEvent();

	for (int i = 0; i < curinput->GetNEvents(); i++) {
	double rwweight = 0.0;
	if (fSignalEventFlags[sigcount]) {
	// Get Event Info
	if (!fIsAllSplines) {
	if (fFillNuisanceEvent)
	curinput->GetNuisanceEvent(i);
	else
	curevent = curinput->GetBaseEvent(i);
	} else {
	curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0];
	}

	curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
	curevent->Weight = curevent->RWWeight * curevent->InputWeight;
	rwweight = curevent->Weight;

	coreeventweights[splinecount] = rwweight;
	if (splinecount % countwidth == 0) {
	LOG(REC) << "Processed " << splinecount
	<< " event weights. W = " << rwweight << std::endl;
	}

	// #pragma omp atomic
	splinecount++;
	}

	// #pragma omp atomic
	sigcount++;
	}
	}
	LOG(SAM) << "Processed event weights." << std::endl;

	// #pragma omp barrier

	// Reset Iterators
	inpsig_iter = fSignalEventFlags.begin();
	spline_iter = fSignalEventSplines.begin();
	box_iter = fSignalEventBoxes.begin();
	samsig_iter = fSampleSignalFlags.begin();
	int nsplineweights = splinecount;
	splinecount = 0;

	// Start of Fast Event Loop ============================

	// Start input iterators
	// Loop over number of inputs
	for (int ispline = 0; ispline < nsplineweights; ispline++) {
	double rwweight = coreeventweights[ispline];

	// Get iterators for this event
	std::vector<bool>::iterator subsamsig_iter = (*samsig_iter).begin();
	std::vector<MeasurementVariableBox*>::iterator subbox_iter =
	(*box_iter).begin();

	// Loop over all sub measurements.
	std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
	for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) {
	MeasurementBase* curmeas = (*meas_iter);

	// If event flagged as signal for this sample fill from the box.
	if (*subsamsig_iter) {
	curmeas->SetSignal(true);
	curmeas->FillHistogramsFromBox((*subbox_iter), rwweight);

	// Move onto next box if there is one.
	subbox_iter++;
	fillcount++;
	}
	}

	if (ispline % countwidth == 0) {
	LOG(REC) << "Filled " << ispline << " sample weights." << std::endl;
	}

	// Iterate over the main signal event containers.
	samsig_iter++;
	box_iter++;
	spline_iter++;
	splinecount++;
	}
	// End of Fast Event Loop ===================

	LOG(SAM) << "Filled sample distributions." << std::endl;

	// Now loop over all Measurements
	// Convert Binned events
	iterSam = fSamples.begin();
	for (; iterSam != fSamples.end(); iterSam++) {
	MeasurementBase* exp = (*iterSam);
	exp->ConvertEventRates();
	}

	// Cleanup coreeventweights
	if (fIsAllSplines) {
	delete coreeventweights;
	}

	// Print some reconfigure profiling.
	LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
	LOG(REC) << "Time taken ReconfigureFastUsingManager() : "
	<< time(NULL) - timestart << std::endl;
	}

	//***************************************************
	void JointFCN::Write() {
	//***************************************************

	// Save a likelihood/ndof plot
	LOG(MIN) << "Writing likelihood plot.." << std::endl;
	std::vector<double> likes;
	std::vector<double> ndofs;
	std::vector<std::string> names;
	for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
	iter++) {
	MeasurementBase* exp = *iter;
	double like = exp->GetLikelihood();
	double ndof = exp->GetNDOF();
	std::string name = exp->GetName();
	likes.push_back(like);
	ndofs.push_back(ndof);
	names.push_back(name);
	}
	TH1D likehist = TH1D("likelihood_hist", "likelihood_hist",
	likes.size(), 0.0, double(likes.size()));
	TH1D ndofhist = TH1D("ndof_hist", "ndof_hist",
	ndofs.size(), 0.0, double(ndofs.size()));
	TH1D divhist = TH1D("likedivndof_hist", "likedivndof_hist",
	likes.size(), 0.0, double(likes.size()));
	for (size_t i = 0; i < likehist.GetNbinsX(); i++) {
	likehist.SetBinContent(i + 1, likes[i]);
	ndofhist.SetBinContent(i + 1, ndofs[i]);
	if (ndofs[i] != 0.0) {
	divhist.SetBinContent(i + 1, likes[i] / ndofs[i]);
	}
	likehist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
	ndofhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
	divhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str());
	}
	likehist.Write();
	ndofhist.Write();
	divhist.Write();

	// Loop over individual experiments and call Write
	LOG(MIN) << "Writing each of the data classes..." << std::endl;
	for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
	iter++) {
	MeasurementBase* exp = *iter;
	exp->Write();
	}

	// Save Pull Terms
	for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
	ParamPull* pull = *iter;
	pull->Write();
	}

	if (FitPar::Config().GetParB("EventManager")) {
	// Get list of inputs
	std::map<int, InputHandlerBase*> fInputs =
	FitBase::EvtManager().GetInputs();
	std::map<int, InputHandlerBase*>::const_iterator iterInp;

	for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) {
	InputHandlerBase* input = (iterInp->second);

	input->GetFluxHistogram()->Write();
	input->GetXSecHistogram()->Write();
	input->GetEventHistogram()->Write();
	}
	}
	};

	//***************************************************
	void JointFCN::SetFakeData(std::string fakeinput) {
	//***************************************************

	LOG(MIN) << "Setting fake data from " << fakeinput << std::endl;
	for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
	iter++) {
	MeasurementBase* exp = *iter;
	exp->SetFakeDataValues(fakeinput);
	}

	return;
	}

	//***************************************************
	void JointFCN::ThrowDataToy() {
	//***************************************************

	for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
	iter++) {
	MeasurementBase* exp = *iter;
	exp->ThrowDataToy();
	}

	return;
	}
	+
	+//***************************************************
	+std::vector<std::string> JointFCN::GetAllNames() {
	+//***************************************************
	+
	+ // Vect of all likelihoods and total
	+ std::vector<std::string> namevect;
	+
	+ // Loop over samples first
	+ for (MeasListConstIter iter = fSamples.begin();
	+ iter != fSamples.end();
	+ iter++) {
	+ MeasurementBase* exp = *iter;
	+
	+ // Get Likelihoods and push to vector
	+ namevect.push_back(exp->GetName());
	+ }
	+
	+
	+ // Loop over pulls second
	+ for (PullListConstIter iter = fPulls.begin();
	+ iter != fPulls.end();
	+ iter++) {
	+ ParamPull* pull = *iter;
	+
	+ // Push back to vector
	+ namevect.push_back(pull->GetName());
	+ }
	+
	+ // Finally add the total
	+ namevect.push_back("total");
	+
	+ return namevect;
	+}
	+
	+//***************************************************
	+std::vector<double> JointFCN::GetAllLikelihoods() {
	+//***************************************************
	+
	+ // Vect of all likelihoods and total
	+ std::vector<double> likevect;
	+ double total_likelihood = 0.0;
	+ LOG(MIN) << "Likelihoods : " << std::endl;
	+
	+ // Loop over samples first
	+ for (MeasListConstIter iter = fSamples.begin();
	+ iter != fSamples.end();
	+ iter++) {
	+ MeasurementBase* exp = *iter;
	+
	+ // Get Likelihoods and push to vector
	+ double singlelike = exp->GetLikelihood();
	+ likevect.push_back(singlelike);
	+ total_likelihood += singlelike;
	+
	+ // Print Out
	+ LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : "
	+ << singlelike << std::endl;
	+ }
	+
	+
	+ // Loop over pulls second
	+ for (PullListConstIter iter = fPulls.begin();
	+ iter != fPulls.end();
	+ iter++) {
	+ ParamPull* pull = *iter;
	+
	+ // Push back to vector
	+ double singlelike = pull->GetLikelihood();
	+ likevect.push_back(singlelike);
	+ total_likelihood += singlelike;
	+
	+ // Print Out
	+ LOG(MIN) << "-> " << std::left << std::setw(40) << pull->GetName() << " : "
	+ << singlelike << std::endl;
	+
	+ }
	+
	+ // Finally add the total likelihood
	+ likevect.push_back(total_likelihood);
	+
	+ return likevect;
	+}
	+
	+//***************************************************
	+std::vector<int> JointFCN::GetAllNDOF() {
	+//***************************************************
	+
	+ // Vect of all ndof and total
	+ std::vector<int> ndofvect;
	+ int total_ndof = 0;
	+
	+ // Loop over samples first
	+ for (MeasListConstIter iter = fSamples.begin();
	+ iter != fSamples.end();
	+ iter++) {
	+ MeasurementBase* exp = *iter;
	+
	+ // Get Likelihoods and push to vector
	+ int singlendof = exp->GetNDOF();
	+ ndofvect.push_back(singlendof);
	+ total_ndof += singlendof;
	+ }
	+
	+
	+ // Loop over pulls second
	+ for (PullListConstIter iter = fPulls.begin();
	+ iter != fPulls.end();
	+ iter++) {
	+ ParamPull* pull = *iter;
	+
	+ // Push back to vector
	+ int singlendof = pull->GetNDOF();
	+ ndofvect.push_back(singlendof);
	+ total_ndof += singlendof;
	+ }
	+
	+ // Finally add the total ndof
	+ ndofvect.push_back(total_ndof);
	+
	+ return ndofvect;
	+}
	diff --git a/src/FCN/JointFCN.h b/src/FCN/JointFCN.h
	index bf87424..1f47f28 100755
	--- a/src/FCN/JointFCN.h
	+++ b/src/FCN/JointFCN.h
	@@ -1,176 +1,180 @@
	#ifndef _JOINT_FCN_H_
	#define _JOINT_FCN_H_
	/*!
	* \addtogroup FCN
	* @{
	*/

	#include <iostream>
	#include <vector>
	#include <fstream>
	#include <list>

	// ROOT headers
	#include "TTree.h"
	#include "TH1D.h"
	#include "TH2D.h"
	#include <TMatrixDSym.h>
	#include "TGraphErrors.h"
	#include <TVectorD.h>
	#include <TMath.h>

	#include "FitUtils.h"

	// External fitter headers
	#include "MeasurementBase.h"
	#include "SampleList.h"

	#define GetCurrentDir getcwd
	#include "EventManager.h"
	#include "Math/Functor.h"
	#include "ParamPull.h"

	#include "NuisConfig.h"
	#include "NuisKey.h"
	#include "MeasurementVariableBox.h"
	#include "MeasurementVariableBox1D.h"

	using namespace FitUtils;
	using namespace FitBase;
	//! Main FCN Class which ROOT's joint function needs to evaulate the chi2 at each stage of the fit.
	class JointFCN
	{
	- public:
	+public:

	//! Constructor
	//! cardfile = Path to input card file listing samples
	- JointFCN(std::vector<nuiskey> samplekeys, TFile* outfile=NULL);
	- JointFCN(TFile* outfile=NULL); // Loads from global config
	+ JointFCN(std::vector<nuiskey> samplekeys, TFile* outfile = NULL);
	+ JointFCN(TFile* outfile = NULL); // Loads from global config
	//! Destructor
	~JointFCN();

	//! Create sample list from cardfile
	void LoadSamples(std::vector<nuiskey> samplekeys);
	void LoadPulls(std::vector<nuiskey> pullkeys);

	//! Main Likelihood evaluation FCN
	double DoEval(const double *x);

	//! Func Wrapper for ROOT
	inline double operator() (const std::vector<double> & x) {
	double* x_array = new double[x.size()];
	return this->DoEval(x_array);
	};

	//! Func Wrapper for ROOT
	inline double operator() (const double *x) {
	return this->DoEval(x);
	};

	//! Create a TTree to save all dial value iterations for this FCN
	void CreateIterationTree(std::string name, FitWeight* rw);

	//! Fills dial values and sample likelihoods into tree
	void FillIterationTree(FitWeight* rw);

	//! Writes TTree to fOutput directory
	void WriteIterationTree();

	//! Deletes TTree
	void DestroyIterationTree();

	//! Get Degrees of Freedom for samples (NBins)
	int GetNDOF();

	//! Return NDOF wrapper
	inline unsigned int NDim() {return this->GetNDOF();};

	//! Reconfigure samples where we force all events to be looped over.
	void ReconfigureAllEvents() ;

	//! Call Reconfigure on samples.
	//! Choice of reconfigure type depends on whether dials have changed
	//! and the MC has been filled.
	void ReconfigureSamples(bool fullconfig = false);

	//! Call reconfigure on only signal events (defaults to all events if CurIter=0)
	void ReconfigureSignal();

	//! Gets likelihood for all samples in FCN (assuming uncorrelated)
	double GetLikelihood();

	//! Returns list of pointers to the samples
	- inline std::list<MeasurementBase*> GetSampleList(){ return fSamples; }
	+ inline std::list<MeasurementBase*> GetSampleList() { return fSamples; }

	//! Return list of pointers to all the pulls
	- inline std::list<ParamPull*> GetPullList(){ return fPulls; };
	+ inline std::list<ParamPull*> GetPullList() { return fPulls; };

	//! Write all samples to output DIR
	void Write();

	//! Set Fake data from file/MC
	void SetFakeData(std::string fakeinput);

	//! Reconfigure looping over duplicate inputs
	void ReconfigureUsingManager();

	//! Reconfigure Fast looping over duplicate inputs
	void ReconfigureFastUsingManager();


	/// Throws data according to current stats
	void ThrowDataToy();

	-std::vector<MeasurementBase*> GetSubSampleList();
	-std::vector<InputHandlerBase*> GetInputList();
	+ std::vector<MeasurementBase*> GetSubSampleList();
	+ std::vector<InputHandlerBase*> GetInputList();

	- private:
	+ std::vector<std::string> GetAllNames();
	+ std::vector<double> GetAllLikelihoods();
	+ std::vector<int> GetAllNDOF();
	+
	+private:

	//! Append the experiments to include in the fit to this list
	std::list<MeasurementBase*> fSamples;

	//! Append parameter pull terms to include penalties in the fit to this list
	std::list<ParamPull*> fPulls;

	TDirectory *fOutputDir; //!< Directory to save contents

	std::string fCardFile; //!< Input Card text file

	bool fDialChanged; //!< Flag for if RW dials changed
	UInt_t fCurIter; //!< Counter for how many times reconfigure called
	bool fMCFilled; //!< Check MC has at least been filled once

	bool fIterationTree; //!< Tree to save RW values on each function call
	int fNDials; //!< Number of RW Dials in FitWeight
	double* fDialVals; //!< Current Values of RW Dials
	double fLikelihood; //!< Current likelihood for joint sample likelihood
	double fNDOF; //!< Total NDOF
	double* fSampleLikes; //!< Likelihoods for each individual measurement in list
	int * fSampleNDOF; //!< NDOF for each individual measurement in list

	bool fUsingEventManager; //!< Flag for doing joint comparisons

	std::vector< std::vector<float> > fSignalEventSplines;
	std::vector< std::vector<MeasurementVariableBox*> > fSignalEventBoxes;
	std::vector< bool > fSignalEventFlags;
	std::vector< std::vector<bool> > fSampleSignalFlags;

	std::vector<InputHandlerBase*> fInputList;
	std::vector<MeasurementBase*> fSubSampleList;
	bool fIsAllSplines;


	std::vector< int > fIterationCount;
	std::vector< double > fCurrentValues;
	std::vector< std::string > fNameValues;
	std::vector< std::vector<double> > fIterationValues;
	int fSampleN;
	std::string fIterationTreeName;




	};

	/! @} /
	#endif // _JOINT_FCN_H_
	diff --git a/src/Routines/BayesianRoutines.cxx b/src/Routines/BayesianRoutines.cxx
	new file mode 100755
	index 0000000..05bfa3f
	--- /dev/null
	+++ b/src/Routines/BayesianRoutines.cxx
	@@ -0,0 +1,524 @@
	+// 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 "BayesianRoutines.h"
	+
	+void BayesianRoutines::Init() {
	+
	+ fInputFile = "";
	+ fInputRootFile = NULL;
	+
	+ fOutputFile = "";
	+ fOutputRootFile = NULL;
	+
	+ fStrategy = "BayesianThrows";
	+ fRoutines.clear();
	+ fRoutines.push_back("BayesianThrows");
	+
	+ fCardFile = "";
	+
	+ fFakeDataInput = "";
	+
	+ fSampleFCN = NULL;
	+
	+ fAllowedRoutines = ("f");
	+
	+};
	+
	+BayesianRoutines::~BayesianRoutines() {
	+};
	+
	+BayesianRoutines::BayesianRoutines(int argc, char* argv[]) {
	+
	+ // Initialise Defaults
	+ Init();
	+ nuisconfig configuration = Config::Get();
	+
	+ // Default containers
	+ std::string cardfile = "";
	+ std::string maxevents = "-1";
	+ int errorcount = 0;
	+ int verbocount = 0;
	+ std::vector<std::string> xmlcmds;
	+ std::vector<std::string> configargs;
	+ fNThrows = 250;
	+ fStartThrows = 0;
	+ fThrowString = "";
	+ // Make easier to handle arguments.
	+ std::vector<std::string> args = GeneralUtils::LoadCharToVectStr(argc, argv);
	+ ParserUtils::ParseArgument(args, "-c", fCardFile, true);
	+ ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
	+ ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
	+ ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
	+ ParserUtils::ParseArgument(args, "-t", fNThrows, false, false);
	+ ParserUtils::ParseArgument(args, "-i", xmlcmds);
	+ ParserUtils::ParseArgument(args, "-q", configargs);
	+ ParserUtils::ParseCounter(args, "e", errorcount);
	+ ParserUtils::ParseCounter(args, "v", verbocount);
	+ ParserUtils::CheckBadArguments(args);
	+
	+ // Add extra defaults if none given
	+ if (fCardFile.empty() and xmlcmds.empty()) {
	+ ERR(FTL) << "No input supplied!" << std::endl;
	+ throw;
	+ }
	+
	+ if (fOutputFile.empty() and !fCardFile.empty()) {
	+ fOutputFile = fCardFile + ".root";
	+ ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
	+
	+ } else if (fOutputFile.empty()) {
	+ ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
	+ throw;
	+ }
	+
	+ // Configuration Setup =============================
	+
	+ // Check no comp key is available
	+ if (Config::Get().GetNodes("nuiscomp").empty()) {
	+ fCompKey = Config::Get().CreateNode("nuiscomp");
	+ } else {
	+ fCompKey = Config::Get().GetNodes("nuiscomp")[0];
	+ }
	+
	+ if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
	+ if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile);
	+ if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
	+
	+ // Load XML Cardfile
	+ configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
	+
	+ // Add CMD XML Structs
	+ for (size_t i = 0; i < xmlcmds.size(); i++) {
	+ configuration.AddXMLLine(xmlcmds[i]);
	+ }
	+
	+ // Add Config Args
	+ for (size_t i = 0; i < configargs.size(); i++) {
	+ configuration.OverrideConfig(configargs[i]);
	+ }
	+ if (maxevents.compare("-1")) {
	+ configuration.OverrideConfig("MAXEVENTS=" + maxevents);
	+ }
	+
	+ // Finish configuration XML
	+ configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
	+
	+ // Add Error Verbo Lines
	+ verbocount += Config::Get().GetParI("VERBOSITY");
	+ errorcount += Config::Get().GetParI("ERROR");
	+ std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
	+ std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
	+ SETVERBOSITY(verbocount);
	+
	+ // Proper Setup
	+ if (fStrategy.find("ErrorBands") != std::string::npos \|\|
	+ fStrategy.find("MergeErrors") != std::string::npos) {
	+ fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
	+ }
	+
	+ // fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
	+ SetupSystematicsFromXML();
	+
	+ SetupRWEngine();
	+ SetupFCN();
	+
	+ return;
	+};
	+
	+void BayesianRoutines::SetupSystematicsFromXML() {
	+
	+ LOG(FIT) << "Setting up nuismin" << std::endl;
	+
	+ // Setup Parameters ------------------------------------------
	+ std::vector<nuiskey> parkeys = Config::QueryKeys("parameter");
	+ if (!parkeys.empty()) {
	+ LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
	+ }
	+
	+ for (size_t i = 0; i < parkeys.size(); i++) {
	+ nuiskey key = parkeys.at(i);
	+
	+ // Check for type,name,nom
	+ if (!key.Has("type")) {
	+ ERR(FTL) << "No type given for parameter " << i << std::endl;
	+ throw;
	+ } else if (!key.Has("name")) {
	+ ERR(FTL) << "No name given for parameter " << i << std::endl;
	+ throw;
	+ } else if (!key.Has("nominal")) {
	+ ERR(FTL) << "No nominal given for parameter " << i << std::endl;
	+ throw;
	+ }
	+
	+ // Get Inputs
	+ std::string partype = key.GetS("type");
	+ std::string parname = key.GetS("name");
	+ double parnom = key.GetD("nominal");
	+ double parlow = parnom - 1;
	+ double parhigh = parnom + 1;
	+ double parstep = 1;
	+
	+
	+ // Override if state not given
	+ if (!key.Has("state")) {
	+ key.SetS("state", "FIX");
	+ }
	+
	+ std::string parstate = key.GetS("state");
	+
	+ // Extra limits
	+ if (key.Has("low")) {
	+ parlow = key.GetD("low");
	+ parhigh = key.GetD("high");
	+ parstep = key.GetD("step");
	+
	+ LOG(FIT) << "Read " << partype << " : "
	+ << parname << " = "
	+ << parnom << " : "
	+ << parlow << " < p < " << parhigh
	+ << " : " << parstate << std::endl;
	+ } else {
	+ LOG(FIT) << "Read " << partype << " : "
	+ << parname << " = "
	+ << parnom << " : "
	+ << parstate << std::endl;
	+ }
	+
	+ // Run Parameter Conversion if needed
	+ if (parstate.find("ABS") != std::string::npos) {
	+ parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
	+ parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
	+ parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
	+ parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
	+ } else if (parstate.find("FRAC") != std::string::npos) {
	+ parnom = FitBase::RWFracToSigma( partype, parname, parnom );
	+ parlow = FitBase::RWFracToSigma( partype, parname, parlow );
	+ parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
	+ parstep = FitBase::RWFracToSigma( partype, parname, parstep );
	+ }
	+
	+ // Push into vectors
	+ fParams.push_back(parname);
	+
	+ fTypeVals[parname] = FitBase::ConvDialType(partype);;
	+ fStartVals[parname] = parnom;
	+ fCurVals[parname] = parnom;
	+
	+ fErrorVals[parname] = 0.0;
	+
	+ fStateVals[parname] = parstate;
	+ bool fixstate = parstate.find("FIX") != std::string::npos;
	+ fFixVals[parname] = fixstate;
	+ fStartFixVals[parname] = fFixVals[parname];
	+
	+ fMinVals[parname] = parlow;
	+ fMaxVals[parname] = parhigh;
	+ fStepVals[parname] = parstep;
	+
	+ }
	+
	+ // Setup Samples ----------------------------------------------
	+ std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
	+ if (!samplekeys.empty()) {
	+ LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
	+ }
	+
	+ for (size_t i = 0; i < samplekeys.size(); i++) {
	+ nuiskey key = samplekeys.at(i);
	+
	+ // Get Sample Options
	+ std::string samplename = key.GetS("name");
	+ std::string samplefile = key.GetS("input");
	+
	+ std::string sampletype =
	+ key.Has("type") ? key.GetS("type") : "DEFAULT";
	+
	+ double samplenorm =
	+ key.Has("norm") ? key.GetD("norm") : 1.0;
	+
	+ // Print out
	+ LOG(FIT) << "Read sample info " << i << " : "
	+ << samplename << std::endl
	+ << "\t\t input -> " << samplefile << std::endl
	+ << "\t\t state -> " << sampletype << std::endl
	+ << "\t\t norm -> " << samplenorm << std::endl;
	+
	+ // If FREE add to parameters otherwise continue
	+ if (sampletype.find("FREE") == std::string::npos) {
	+ continue;
	+ }
	+
	+ // Form norm dial from samplename + sampletype + "_norm";
	+ std::string normname = samplename + "_norm";
	+
	+ // Check normname not already present
	+ if (fTypeVals.find(normname) != fTypeVals.end()) {
	+ continue;
	+ }
	+
	+ // Add new norm dial to list if its passed above checks
	+ fParams.push_back(normname);
	+
	+ fTypeVals[normname] = kNORM;
	+ fStateVals[normname] = sampletype;
	+ fCurVals[normname] = samplenorm;
	+
	+ fErrorVals[normname] = 0.0;
	+
	+ fMinVals[normname] = 0.1;
	+ fMaxVals[normname] = 10.0;
	+ fStepVals[normname] = 0.5;
	+
	+ bool state = sampletype.find("FREE") == std::string::npos;
	+ fFixVals[normname] = state;
	+ fStartFixVals[normname] = state;
	+ }
	+
	+ // Setup Fake Parameters -----------------------------
	+ std::vector<nuiskey> fakekeys = Config::QueryKeys("fakeparameter");
	+ if (!fakekeys.empty()) {
	+ LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
	+ }
	+
	+ for (size_t i = 0; i < fakekeys.size(); i++) {
	+ nuiskey key = fakekeys.at(i);
	+
	+ // Check for type,name,nom
	+ if (!key.Has("name")) {
	+ ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
	+ throw;
	+ } else if (!key.Has("nom")) {
	+ ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
	+ throw;
	+ }
	+
	+ // Get Inputs
	+ std::string parname = key.GetS("name");
	+ double parnom = key.GetD("nom");
	+
	+ // Push into vectors
	+ fFakeVals[parname] = parnom;
	+ }
	+}
	+
	+/*
	+ Setup Functions
	+*/
	+//*************************************
	+void BayesianRoutines::SetupRWEngine() {
	+//*************************************
	+
	+ for (UInt_t i = 0; i < fParams.size(); i++) {
	+ std::string name = fParams[i];
	+ FitBase::GetRW() -> IncludeDial(name, fTypeVals.at(name) );
	+ }
	+ UpdateRWEngine(fStartVals);
	+
	+ return;
	+}
	+
	+//*************************************
	+void BayesianRoutines::SetupFCN() {
	+//*************************************
	+
	+ LOG(FIT) << "Making the jointFCN" << std::endl;
	+ if (fSampleFCN) delete fSampleFCN;
	+ fSampleFCN = new JointFCN(fOutputRootFile);
	+
	+ fInputThrows = fSampleFCN->GetPullList();
	+
	+ return;
	+}
	+
	+/*
	+ Fitting Functions
	+*/
	+//*************************************
	+void BayesianRoutines::UpdateRWEngine(std::map<std::string, double>& updateVals) {
	+//*************************************
	+
	+ for (UInt_t i = 0; i < fParams.size(); i++) {
	+ std::string name = fParams[i];
	+
	+ if (updateVals.find(name) == updateVals.end()) continue;
	+ FitBase::GetRW()->SetDialValue(name, updateVals.at(name));
	+ }
	+
	+ FitBase::GetRW()->Reconfigure();
	+ return;
	+}
	+
	+//*************************************
	+void BayesianRoutines::ThrowParameters() {
	+//*************************************
	+
	+ // Set fThrownVals to all values in currentVals
	+ for (UInt_t i = 0; i < fParams.size(); i++) {
	+ std::string name = fParams.at(i);
	+ fThrownVals[name] = fCurVals[name];
	+ }
	+
	+ for (PullListConstIter iter = fInputThrows.begin();
	+ iter != fInputThrows.end(); iter++) {
	+ ParamPull* pull = *iter;
	+
	+ pull->ThrowCovariance();
	+ TH1D dialhist = pull->GetDataHist();
	+
	+ for (int i = 0; i < dialhist.GetNbinsX(); i++) {
	+ std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i + 1));
	+ if (fCurVals.find(name) != fCurVals.end()) {
	+ fThrownVals[name] = dialhist.GetBinContent(i + 1);
	+ }
	+ }
	+
	+ // Reset throw incase pulls are calculated.
	+ pull->ResetToy();
	+
	+ }
	+
	+ // Now update Parameters
	+ UpdateRWEngine(fThrownVals);
	+
	+ // Update Pulls
	+ for (PullListConstIter iter = fInputThrows.begin();
	+ iter != fInputThrows.end(); iter++) {
	+ ParamPull* pull = *iter;
	+ pull->Reconfigure();
	+ }
	+
	+ return;
	+};
	+
	+//*************************************
	+void BayesianRoutines::Run() {
	+//*************************************
	+
	+ std::cout << "Running routines " << std::endl;
	+ fRoutines = GeneralUtils::ParseToStr(fStrategy, ",");
	+
	+ for (UInt_t i = 0; i < fRoutines.size(); i++) {
	+
	+ std::string routine = fRoutines.at(i);
	+ LOG(FIT) << "Running Routine: " << routine << std::endl;
	+
	+ if (!routine.compare("BayesianThrows")) GenerateThrows();
	+ else THROW("UNKNOWN ROUTINE " << routine);
	+ }
	+
	+ return;
	+}
	+
	+//*************************************
	+void BayesianRoutines::GenerateThrows() {
	+//*************************************
	+
	+ // Create a new output file
	+ TFile* outfile = new TFile((fOutputFile + ".throws.root").c_str(), "RECREATE");
	+ outfile->cd();
	+
	+ int nthrows = fNThrows;
	+
	+ // Setting Seed
	+ // Matteo Mazzanti's Fix
	+ struct timeval mytime;
	+ gettimeofday(&mytime, NULL);
	+ Double_t seed = time(NULL) + int(getpid()) + (mytime.tv_sec * 1000.) + (mytime.tv_usec / 1000.);
	+ gRandom->SetSeed(seed);
	+ LOG(FIT) << "Using Seed : " << seed << std::endl;
	+ LOG(FIT) << "nthrows = " << nthrows << std::endl;
	+
	+ // Run the Initial Reconfigure
	+ LOG(FIT) << "Making nominal prediction " << std::endl;
	+ TDirectory* nominal = (TDirectory*) outfile->mkdir("nominal");
	+ nominal->cd();
	+ UpdateRWEngine(fStartVals);
	+ fSampleFCN->ReconfigureUsingManager();
	+ fSampleFCN->Write();
	+
	+ // Create an iteration tree inside SampleFCN
	+ fSampleFCN->CreateIterationTree("error_iterations", FitBase::GetRW());
	+
	+ // Create a new iteration TTree
	+ TTree* LIKETREE = new TTree("likelihood", "likelihood");
	+ std::vector<std::string> likenames = fSampleFCN->GetAllNames();
	+ std::vector<double> likevals = fSampleFCN->GetAllLikelihoods();
	+ std::vector<int> likendof = fSampleFCN->GetAllNDOF();
	+ double* LIKEVALS = new double[likevals.size()];
	+ int* LIKENDOF = new int[likendof.size()];
	+
	+ for (size_t i = 0; i < likendof.size(); i++) {
	+ LIKETREE->Branch( (likenames[i] + "_likelihood" ).c_str(), &LIKEVALS[i],
	+ (likenames[i] + "_likelihood/D").c_str() );
	+ LIKETREE->Branch( (likenames[i] + "_ndof" ).c_str(), &LIKENDOF[i],
	+ (likenames[i] + "_ndof/I").c_str() );
	+ LIKENDOF[i] = likendof[i];
	+ }
	+
	+ likenames .clear();
	+ likevals .clear();
	+ likendof .clear();
	+
	+ double* PARAMVALS = new double[fParams.size()];
	+ for (size_t i = 0; i < fParams.size(); i++){
	+ LIKETREE->Branch( fParams[i].c_str(), &PARAMVALS[i], (fParams[i] + "/D").c_str() );
	+ }
	+
	+ // Run Throws and save
	+ for (Int_t i = 0; i < nthrows; i++) {
	+
	+ // Skip the start throw
	+ if (i == 0) continue;
	+ LOG(FIT) << "Throw " << i << " ================================" << std::endl;
	+
	+ // Throw Parameters
	+ ThrowParameters();
	+ FitBase::GetRW()->Print();
	+
	+ // Get Parameter Values
	+ for (size_t i = 0; i < fParams.size(); i++){
	+ PARAMVALS[i] = fThrownVals[fParams[i]];
	+ }
	+
	+ // Run Sample Prediction
	+ fSampleFCN->ReconfigureFastUsingManager();
	+
	+ // Get vector of likelihoods/ndof
	+ std::vector<double> likevals = fSampleFCN->GetAllLikelihoods();
	+ for (size_t i = 0; i < likevals.size(); i++) {
	+ LIKEVALS[i] = likevals[i];
	+ }
	+
	+ // Save to TTree
	+ LIKETREE->Fill();
	+
	+ // Save the FCN
	+ // if (fSavePredictions){ SaveSamplePredictions(); }
	+ LOG(FIT) << "END OF THROW ================================" << std::endl;
	+
	+ }
	+
	+ // Finish up
	+ outfile->cd();
	+ LIKETREE->Write();
	+ outfile->Close();
	+ delete LIKEVALS;
	+ delete LIKENDOF;
	+ delete PARAMVALS;
	+}
	diff --git a/src/Routines/BayesianRoutines.h b/src/Routines/BayesianRoutines.h
	new file mode 100755
	index 0000000..3e56888
	--- /dev/null
	+++ b/src/Routines/BayesianRoutines.h
	@@ -0,0 +1,183 @@
	+// 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 BAYESIAN_ROUTINES_H
	+#define BAYESIAN_ROUTINES_H
	+
	+/*!
	+ * \addtogroup Minimizer
	+ * @{
	+ */
	+
	+#include "TH1.h"
	+#include "TF1.h"
	+#include "TMatrixD.h"
	+#include "TVectorD.h"
	+#include "TSystem.h"
	+#include "TFile.h"
	+#include "TProfile.h"
	+
	+#include <sys/time.h>
	+#include <vector>
	+#include <string>
	+#include <iostream>
	+#include <sstream>
	+#include <cstring>
	+
	+#include "FitEvent.h"
	+#include "JointFCN.h"
	+#include "FitParameters.h"
	+#include "ParserUtils.h"
	+
	+enum minstate {
	+ kErrorStatus = -1,
	+ kGoodStatus,
	+ kFitError,
	+ kNoChange,
	+ kFitFinished,
	+ kFitUnfinished,
	+ kStateChange,
	+};
	+
	+//*************************************
	+//! Collects all possible fit routines into a single class to avoid repeated code
	+class BayesianRoutines{
	+//*************************************
	+
	+public:
	+
	+ /*
	+ Constructor/Destructor
	+ */
	+
	+ //! Constructor reads in arguments given at the command line for the fit here.
	+ BayesianRoutines(int argc, char* argv[]);
	+
	+ //! Default destructor
	+ ~BayesianRoutines();
	+
	+ //! Reset everything to default/NULL
	+ void Init();
	+
	+ /*
	+ Input Functions
	+ */
	+
	+ //! Splits the arguments ready for initial setup
	+ void ParseArgs(int argc, char* argv[]);
	+
	+ //! Sorts out configuration and verbosity right at the very start.
	+ //! Calls readCard to set everything else up.
	+ void InitialSetup();
	+
	+ /*
	+ Setup Functions
	+ */
	+ void SetupSystematicsFromXML();
	+
	+ //! Setups up our custom RW engine with all the parameters passed in the card file
	+ void SetupRWEngine();
	+
	+ //! Setups up the jointFCN.
	+ void SetupFCN();
	+
	+ /*
	+ Fitting Functions
	+ */
	+
	+ //! Main function to actually start iterating over the different required fit routines
	+ void Run();
	+
	+ //! Given a new map change the values that the RW engine is currently set to
	+ void UpdateRWEngine(std::map<std::string,double>& updateVals);
	+
	+ //! Given a single routine (see tutorial for options) run that fit routine now.
	+ int RunFitRoutine(std::string routine);
	+
	+ //! Throw the current covariance of dial values we have, and fill the thrownVals and thrownNorms maps.
	+ //! If uniformly is true parameters will be thrown uniformly between their upper and lower limits.
	+ void ThrowParameters();
	+
	+ //! Run Throws
	+ void GenerateThrows();
	+
	+protected:
	+
	+ //! Our Custom ReWeight Object
	+ FitWeight* rw;
	+
	+ std::string fOutputFile;
	+ std::string fInputFile;
	+
	+ TFile* fInputRootFile;
	+ TFile* fOutputRootFile;
	+
	+ //! Flag for whether the fit should be continued if an output file is already found.
	+ bool fitContinue;
	+
	+ //! Minimizer Object for handling roots different minimizer methods
	+ JointFCN* fSampleFCN;
	+
	+ int nfreepars;
	+
	+ std::string fCardFile;
	+
	+ std::string fStrategy;
	+ std::vector<std::string> fRoutines;
	+ std::string fAllowedRoutines;
	+
	+ std::string fFakeDataInput;
	+
	+ // Input Dial Vals
	+ //! Vector of dial names
	+ std::vector<std::string> fParams;
	+ std::map<std::string, std::string> fStateVals;
	+ std::map<std::string, double> fStartVals;
	+ std::map<std::string, double> fCurVals;
	+ std::map<std::string, double> fErrorVals;
	+ std::map<std::string, double> fMinVals;
	+ std::map<std::string, double> fMaxVals;
	+ std::map<std::string, double> fStepVals;
	+ std::map<std::string, int> fTypeVals;
	+ std::map<std::string, bool> fFixVals;
	+ std::map<std::string, bool> fStartFixVals;
	+
	+ //! Vector of fake parameter names
	+ std::map<std::string,double> fFakeVals;
	+
	+ //! Map of thrown parameter names and values (After ThrowCovariance)
	+ std::map<std::string,double> fThrownVals;
	+
	+
	+ std::list <ParamPull*> fInputThrows; //!< Pointers to pull terms
	+ std::vector <TH1D> fInputDials; //!< Vector of Input Histograms
	+ std::vector <TMatrixDSym> fInputCovar; //!< Vector of Input Covariances
	+
	+ nuiskey fCompKey;
	+ std::vector<std::string> fThrowList;
	+ std::string fThrowString;
	+
	+ int fNThrows;
	+ int fStartThrows;
	+
	+
	+};
	+
	+/! @} /
	+#endif
	diff --git a/src/Routines/CMakeLists.txt b/src/Routines/CMakeLists.txt
	index 2c5cf1c..a87e4a0 100644
	--- a/src/Routines/CMakeLists.txt
	+++ b/src/Routines/CMakeLists.txt
	@@ -1,66 +1,68 @@
	# 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
	ComparisonRoutines.cxx
	SystematicRoutines.cxx
	SplineRoutines.cxx
	+BayesianRoutines.cxx
	)

	if(USE_MINIMIZER)
	set(IMPLFILES ${IMPLFILES};MinimizerRoutines.cxx)
	endif()

	set(HEADERFILES
	ComparisonRoutines.h
	SystematicRoutines.h
	+BayesianRoutines.h
	SplineRoutines.h
	)

	if(USE_MINIMIZER)
	set(HEADERFILES ${HEADERFILES};MinimizerRoutines.h)
	endif()

	set(LIBNAME Routines)

	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(${EXP_INCLUDE_DIRECTORIES})

	include_directories(${MINIMUM_INCLUDE_DIRECTORIES})

	include_directories(${CMAKE_SOURCE_DIR}/src/FCN)
	include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies)

	set_target_properties(${LIBNAME} PROPERTIES VERSION
	"${NUISANCE_VERSION_MAJOR}.${NUISANCE_VERSION_MINOR}.${NUISANCE_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)

File Metadata

Mime Type: text/x-diff
Expires: Sun, Feb 23, 3:04 PM (31 m, 23 s)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4486116
Default Alt Text: (81 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions