No OneTemporary
Actions

Size

73 KB

Referenced Files

None

Subscribers

None

View Options

	diff --git a/src/Routines/SplineRoutines.cxx b/src/Routines/SplineRoutines.cxx
	index 996122e1..ea67d2df 100755
	--- a/src/Routines/SplineRoutines.cxx
	+++ b/src/Routines/SplineRoutines.cxx
	@@ -1,1616 +1,1637 @@
	// 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 "StatusMessage.h"

	#include "SplineRoutines.h"

	/*
	Constructor/Destructor
	*/
	//************************
	void SplineRoutines::Init() {
	//************************

	fStrategy = "SaveEvents";
	fRoutines.clear();

	fCardFile = "";

	fSampleFCN = NULL;
	fRW = NULL;

	fAllowedRoutines = ("SaveEvents,TestEvents,SaveSplineEvents");
	};

	//*************************************
	SplineRoutines::~SplineRoutines() {
	//*************************************
	};

	/*
	Input Functions
	*/
	//*************************************
	SplineRoutines::SplineRoutines(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;

	// 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, "-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
	nuiskey fCompKey;
	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);
	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")) {
	std::cout << "[ NUISANCE ] : Overriding " << "MAXEVENTS=" + maxevents << std::endl;
	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;
	FitPar::log_verb = verbocount;
	LOG_VERB(verbocount);
	ERR_VERB(errorcount);

	// Starting Setup
	// ---------------------------
	SetupRWEngine();

	return;
	};


	/*
	Setup Functions
	*/
	//*************************************
	void SplineRoutines::SetupRWEngine() {
	//*************************************

	fRW = new FitWeight("splineweight");
	// std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
	std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");

	// Add Parameters
	for (size_t i = 0; i < parameterkeys.size(); i++) {
	nuiskey key = parameterkeys[i];

	std::string parname = key.GetS("name");
	std::string partype = key.GetS("type");
	double nom = key.GetD("nominal");

	fRW->IncludeDial( key.GetS("name"),
	FitBase::ConvDialType(key.GetS("type")), nom);
	fRW->SetDialValue( key.GetS("name"), key.GetD("nominal") );

	}
	fRW->Reconfigure();

	return;
	}


	/*
	Fitting Functions
	*/
	//*************************************
	void SplineRoutines::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;
	fRW->SetDialValue(name, updateVals.at(name));
	}

	fRW->Reconfigure();
	return;
	}

	//*************************************
	void SplineRoutines::Run() {
	//*************************************
	std::cout << "Running " << std::endl;

	// Parse given routines
	fRoutines = GeneralUtils::ParseToStr(fStrategy, ",");
	if (fRoutines.empty()) {
	ERR(FTL) << "Trying to run ComparisonRoutines with no routines given!" << std::endl;
	throw;
	}

	for (size_t i = 0; i < fRoutines.size(); i++) {

	LOG(FIT) << "Running Routine: " << fRoutines[i] << std::endl;
	std::string rout = fRoutines[i];
	if (!rout.compare("SaveEvents")) SaveEvents();
	else if (!rout.compare("TestEvents")) TestEvents();
	else if (!rout.compare("GenerateEventSplines")){ GenerateEventSplines(); }
	// else if (!rout.compare("GenerateEventWeights")) { GenerateEventWeights(); }
	// else if (!rout.compare("BuildEventSplines")) { BuildEventSplines(); }
	else if (!rout.compare("TestSplines_1DEventScan")) TestSplines_1DEventScan();
	else if (!rout.compare("TestSplines_NDEventThrow")) TestSplines_NDEventThrow();
	else if (!rout.compare("SaveSplinePlots")) SaveSplinePlots();
	else if (!rout.compare("TestSplines_1DLikelihoodScan")) TestSplines_1DLikelihoodScan();
	else if (!rout.compare("TestSplines_NDLikelihoodThrow")) TestSplines_NDLikelihoodThrow();

	}


	}

	//*************************************
	void SplineRoutines::SaveEvents() {
	//*************************************

	if (fRW) delete fRW;
	SetupRWEngine();
	fRW->Reconfigure();
	fRW->Print();

	// Generate a set of nominal events
	// Method, Loop over inputs, create input handler, then create a ttree
	std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
	for (size_t i = 0; i < eventkeys.size(); i++) {
	nuiskey key = eventkeys.at(i);

	// Get I/O
	std::string inputfilename = key.GetS("input");
	if (inputfilename.empty()) {
	ERR(FTL) << "No input given for set of input events!" << std::endl;
	throw;
	}

	std::string outputfilename = key.GetS("output");
	if (outputfilename.empty()) {
	outputfilename = inputfilename + ".nuisance.root";
	ERR(FTL) << "No output give for set of output events! Saving to "
	<< outputfilename << std::endl;
	}

	// Make new outputfile
	TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
	outputfile->cd();

	// Make a new input handler
	std::vector<std::string> file_descriptor =
	GeneralUtils::ParseToStr(inputfilename, ":");
	if (file_descriptor.size() != 2) {
	ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
	<< "\". expected \"FILETYPE:file.root\"" << std::endl;
	throw;
	}
	InputUtils::InputType inptype =
	InputUtils::ParseInputType(file_descriptor[0]);

	InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);

	// Get info from inputhandler
	int nevents = input->GetNEvents();
	int countwidth = (nevents / 10);
	FitEvent* nuisevent = input->FirstNuisanceEvent();

	// Setup a TTree to save the event
	outputfile->cd();
	TTree* eventtree = new TTree("nuisance_events", "nuisance_events");
	nuisevent->AddBranchesToTree(eventtree);

	// Loop over all events and fill the TTree
	int i = 0;
	// int countwidth = nevents / 5;

	while (nuisevent) {

	// Get Event Weight
	nuisevent->RWWeight = fRW->CalcWeight(nuisevent);
	// if (nuisevent->RWWeight != 1.0){
	// std::cout << "Weight = " << nuisevent->RWWeight << std::endl;
	// }
	// Save everything
	eventtree->Fill();

	// Logging
	if (i % countwidth == 0) {
	LOG(REC) << "Saved " << i << "/" << nevents
	<< " nuisance events. [M, W] = ["
	<< nuisevent->Mode << ", " << nuisevent->RWWeight << "]" << std::endl;
	}

	// iterate
	nuisevent = input->NextNuisanceEvent();
	i++;
	}

	// Save flux and close file
	outputfile->cd();
	eventtree->Write();
	input->GetFluxHistogram()->Write("nuisance_fluxhist");
	input->GetEventHistogram()->Write("nuisance_eventhist");

	// Close Output
	outputfile->Close();

	// Delete Inputs
	delete input;
	}

	// remove Keys
	eventkeys.clear();

	// Finished
	LOG(FIT) << "Finished processing all nuisance events." << std::endl;
	}

	//*************************************
	void SplineRoutines::TestEvents() {
	//*************************************

	LOG(FIT) << "Testing events." << std::endl;

	// Create a new file for the test samples
	if (!fOutputRootFile) {
	fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
	}

	// Loop over all tests
	int count = 0;
	std::vector<nuiskey> testkeys = Config::QueryKeys("sampletest");
	for (std::vector<nuiskey>::iterator iter = testkeys.begin();
	iter != testkeys.end(); iter++) {
	nuiskey key = (*iter);

	// 0. Create new measurement list
	std::list<MeasurementBase*> samplelist;

	// 1. Build Sample From Events
	std::string samplename = key.GetS("name");
	std::string eventsid = key.GetS("inputid");
	nuiskey eventskey = Config::QueryLastKey("events", "id=" + eventsid);
	std::string rawfile = eventskey.GetS("input");
	LOG(FIT) << "Creating sample " << samplename << std::endl;
	MeasurementBase* rawsample = SampleUtils::CreateSample(samplename, rawfile, "", "", FitBase::GetRW());

	// 2. Build Sample From Nuisance Events
	std::string eventsfile = eventskey.GetS("output");
	LOG(FIT) << "Creating Fit Eevnt Sample " << samplename << " " << eventsfile << std::endl;
	MeasurementBase* nuissample = SampleUtils::CreateSample(samplename, "FEVENT:" + eventsfile, "", "", FitBase::GetRW());

	// 3. Make some folders to save stuff
	TDirectory* sampledir = (TDirectory*) fOutputRootFile->mkdir(Form((samplename + "_test_%d").c_str(), count));
	TDirectory* rawdir = (TDirectory*) sampledir->mkdir("raw");
	TDirectory* nuisancedir = (TDirectory*) sampledir->mkdir("nuisance");
	TDirectory* difdir = (TDirectory*) sampledir->mkdir("difference");

	// 4. Reconfigure both
	rawdir->cd();
	rawsample->Reconfigure();
	rawsample->Write();

	nuisancedir->cd();
	nuissample->Reconfigure();
	nuissample->Write();

	// 4. Compare Raw to Nuisance Events

	// Loop over all keyse
	TIter next(rawdir->GetListOfKeys());
	TKey *dirkey;
	while ((dirkey = (TKey*)next())) {

	// If not a 1D/2D histogram skip
	TClass *cl = gROOT->GetClass(dirkey->GetClassName());
	if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;

	// Get TH1* from both dir
	TH1 rawplot = (TH1)rawdir->Get(dirkey->GetName());
	TH1 nuisanceplot = (TH1)nuisancedir->Get(dirkey->GetName());

	// Take Difference
	nuisanceplot->Add(rawplot, -1.0);

	// Save to dif folder
	difdir->cd();
	nuisanceplot->Write();
	}

	// 5. Tidy Up
	samplelist.clear();

	// Iterator
	count++;
	}
	}



	//*************************************
	void SplineRoutines::GenerateEventSplines() {
	//*************************************
	if (fRW) delete fRW;
	SetupRWEngine();

	// Setup the spline reader
	SplineWriter* splwrite = new SplineWriter(fRW);
	std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");

	// Add splines to splinewriter
	for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
	iter != splinekeys.end(); iter++) {
	nuiskey splkey = (*iter);

	// Add Spline Info To Reader
	splwrite->AddSpline(splkey);
	}
	splwrite->SetupSplineSet();


	// Make an ugly list for N cores
	int ncores = FitPar::Config().GetParI("NCORES");//omp_get_max_threads();
	std::vector<SplineWriter*> splwriterlist;

	for (int i = 0; i < ncores; i++){
	SplineWriter* tmpwriter = new SplineWriter(fRW);

	for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
	iter != splinekeys.end(); iter++) {
	nuiskey splkey = (*iter);

	// Add Spline Info To Reader
	tmpwriter->AddSpline(splkey);
	}
	tmpwriter->SetupSplineSet();

	splwriterlist.push_back(tmpwriter);
	}


	// Event Loop
	// Loop over all events and calculate weights for each parameter set.

	// Generate a set of nominal events
	// Method, Loop over inputs, create input handler, then create a ttree
	std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
	for (size_t i = 0; i < eventkeys.size(); i++) {
	nuiskey key = eventkeys.at(i);

	// Get I/O
	std::string inputfilename = key.GetS("input");
	if (inputfilename.empty()) {
	ERR(FTL) << "No input given for set of input events!" << std::endl;
	throw;
	}

	std::string outputfilename = key.GetS("output");
	if (outputfilename.empty()) {
	outputfilename = inputfilename + ".nuisance.root";
	ERR(FTL) << "No output give for set of output events! Saving to "
	<< outputfilename << std::endl;
	}

	// Make new outputfile
	TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
	outputfile->cd();

	// Make a new input handler
	std::vector<std::string> file_descriptor =
	GeneralUtils::ParseToStr(inputfilename, ":");
	if (file_descriptor.size() != 2) {
	ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
	<< "\". expected \"FILETYPE:file.root\"" << std::endl;
	throw;
	}
	InputUtils::InputType inptype =
	InputUtils::ParseInputType(file_descriptor[0]);

	InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);

	// Get info from inputhandler
	int nevents = input->GetNEvents();
	int countwidth = (nevents / 1000);
	FitEvent* nuisevent = input->FirstNuisanceEvent();

	// Setup a TTree to save the event
	outputfile->cd();
	TTree* eventtree = new TTree("nuisance_events", "nuisance_events");

	// Add a flag that allows just splines to be saved.
	nuisevent->AddBranchesToTree(eventtree);

	// Save the spline reader
	splwrite->Write("spline_reader");

	// Setup the spline TTree
	TTree* weighttree = new TTree("weight_tree", "weight_tree");
	splwrite->AddWeightsToTree(weighttree);

	// Make container for all weights
	- int nweights = splwrite->GetNWeights()+1;
	+ int nweights = splwrite->GetNWeights();
	double** weightcont = new double*[nevents];
	for (int k = 0; k < nevents; k++) {
	weightcont[k] = new double[nweights];
	}

	int npar = splwrite->GetNPars();


	int lasttime = time(NULL);

	// Could reorder this to save the weightconts in order instead of reconfiguring per event.
	// Loop over all events and fill the TTree
	while (nuisevent) {

	// std::cout << "Fitting event " << i << std::endl;
	// Calculate the weights for each parameter set
	// splwrite->FitSplinesForEvent(nuisevent);
	splwrite->GetWeightsForEvent(nuisevent, weightcont[i]);
	+ bool hasresponse = false;
	+ for (int j = 0; j < nweights; j++){
	+
	+ if (weightcont[i][j] != 1.0){
	+ // std::cout << "Non Zero Weight at " << i << " " << j << std::endl;
	+ hasresponse = true;
	+ } else {
	+ // std::cout << "Empty Weight at " << i << " " << j << std::endl;
	+ }
	+ }
	+ if (!hasresponse){
	+ // std::cout << "Deleting flat response " << nuisevent->Mode << std::endl;
	+ delete weightcont[i];
	+ weightcont[i] = NULL;
	+ }

	// Save everything
	eventtree->Fill();
	weighttree->Fill();
	// splinetree->Fill();

	// nuisevent->Print();
	// std::cout << "Done with event " << i << std::endl;

	// Push weight sets into a the array

	// sleep(4);
	// Logging
	if (i % countwidth == 0) {

	std::ostringstream timestring;
	int timeelapsed = time(NULL) - lasttime;
	if (i != 0 and timeelapsed) {
	lasttime = time(NULL);

	int eventsleft = nevents - i;
	float speed = float(countwidth) / float(timeelapsed);
	float proj = (float(eventsleft) / float(speed)) / 60 / 60;
	timestring << proj << " hours remaining.";

	}
	LOG(REC) << "Saved " << i << "/" << nevents << " nuisance spline weights. " << timestring.str() << std::endl;
	}

	// Iterate
	i++;
	nuisevent = input->NextNuisanceEvent();
	}

	outputfile->cd();
	eventtree->Write();
	weighttree->Write();
	input->GetFluxHistogram()->Write("nuisance_fluxhist");
	input->GetEventHistogram()->Write("nuisance_eventhist");
	outputfile->Close();

	outputfile = new TFile(outputfilename.c_str(), "UPDATE");
	outputfile->cd();

	weighttree = (TTree*) outputfile->Get("weight_tree");
	splwrite->ReadWeightsFromTree(weighttree);


	// Divide weights container into Ncores.
	// Parrallelise this loop checking for what core we are on.
	// for (int i = 0; i < nevents; i++){
	// splwriterlist[int(i / (nevents/4))]->FitSplinesForEvent(coeff);
	// }

	// // Now loop over weights tree
	// for (int i = 0; i < weighttree->GetEntries(); i++) {
	// weighttree->GetEntry(i);
	// splwrite->FitSplinesForEvent();
	// splinetree->Fill();

	// if (i % countwidth == 0) {

	// std::ostringstream timestring;
	// int timeelapsed = time(NULL) - lasttime;
	// if (i != 0 and timeelapsed) {
	// lasttime = time(NULL);

	// int eventsleft = nevents - i;
	// float speed = float(countwidth) / float(timeelapsed);
	// float proj = (float(eventsleft) / float(speed)) / 60 / 60;
	// timestring << proj << " hours remaining.";

	// }
	// LOG(REC) << "Built " << i << "/" << nevents << " nuisance spline events. " << timestring.str() << std::endl;
	// }
	// }

	// Get Splines
	float** allcoeff = new float*[nevents];
	for (int k = 0; k < nevents; k++){
	allcoeff[k] = new float[npar];
	}



	#pragma omp parallel for num_threads(ncores)
	for (int i = 0; i < nevents; i++) {

	//#pragma omp atomic
	// printf("Using Thread %d to build event %d \n", int(omp_get_thread_num()), (int)i );
	// std::cout<< " -> Writer = " << splwriterlist[ i / (nevents/ncores) ] << std::endl;

	// #pragma omp atomic
	- splwriterlist[ int(omp_get_thread_num()) ]->FitSplinesForEvent(weightcont[i], allcoeff[i]);
	+ if (weightcont[i]){
	+ splwriterlist[ int(omp_get_thread_num()) ]->FitSplinesForEvent(weightcont[i], allcoeff[i]);
	+ } else {
	+ for (int j = 0; j < npar; j++){
	+ allcoeff[i][j] = float(0.0);
	+ }
	+ }

	// splwrite->FitSplinesForEvent(weightcont[i], allcoeff[i]);


	if (i % 500 == 0) {

	if (LOG_LEVEL(REC)){
	printf("Using Thread %d to build event %d \n", int(omp_get_thread_num()), (int)i );
	}
	}
	/*

	std::ostringstream timestring;
	int timeelapsed = time(NULL) - lasttime;
	if (i != 0 and timeelapsed) {
	lasttime = time(NULL);

	int eventsleft = nevents - i;
	float speed = float(countwidth) / float(timeelapsed);
	float proj = (float(eventsleft) / float(speed)) / 60 / 60;
	timestring << proj << " hours remaining.";
	timestring << " Using Writer at " << i / (nevents/ncores) << " = " << splwriterlist[ i / (nevents/ncores) ] << std::endl;

	}
	LOG(REC) << "Built " << i << "/" << nevents << " nuisance spline events. " << timestring.str() << std::endl;
	}
	*/
	}

	// Save Splines into TTree
	float* coeff = new float[npar];
	outputfile->cd();
	TTree* splinetree = new TTree("spline_tree", "spline_tree");

	splinetree->Branch("SplineCoeff", coeff, Form("SplineCoeff[%d]/F",npar));

	std::cout << "Saving to the allcoeff" << std::endl;
	for (int k = 0; k < nevents; k++) {
	for (int l = 0; l < npar; l++) {
	coeff[l] = allcoeff[k][l];
	}
	std::cout << "Coeff 0, 1, 2 = " << coeff[0] << " " << coeff[1] << " " << coeff[2] << std::endl;
	splinetree->Fill();
	}

	// Save flux and close file
	outputfile->cd();
	splinetree->Write();

	// Delete the container.
	for (int k = 0; k < nevents; k++) {
	delete weightcont[k];
	}
	delete weightcont;
	delete coeff;

	// Close Output
	outputfile->Close();

	// Delete Inputs
	delete input;
	}

	// remove Keys
	eventkeys.clear();

	}


	//*************************************
	void SplineRoutines::MergeSplines() {
	//*************************************
	// Loop over all 'splinemerge' keys.
	// Add them to the Merger.
	// Call setup splines.

	// Get the key with eventinput
	// - remaining keys should have splineinput
	// - Loop over number of entries.
	// - FillEntry in merger.
	// - Fill NUISANCEEvent into a new TTree.

	SplineMerger* splmerge = new SplineMerger();
	std::vector<nuiskey> splinekeys = Config::QueryKeys("splinemerge");
	for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
	iter != splinekeys.end(); iter++) {
	nuiskey splkey = (*iter);

	TFile* infile = new TFile(splkey.GetS("input").c_str(), "READ");
	splmerge->AddSplineSetFromFile(infile);

	}
	splmerge->SetupSplineSet();

	// Now get Event File
	std::vector<nuiskey> eventkeys = Config::QueryKeys("eventmerge");
	nuiskey key = eventkeys[0];

	std::string inputfilename = key.GetS("input");

	// Make a new input handler
	std::vector<std::string> file_descriptor =
	GeneralUtils::ParseToStr(inputfilename, ":");
	if (file_descriptor.size() != 2) {
	ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
	<< "\". expected \"FILETYPE:file.root\"" << std::endl;
	throw;
	}
	InputUtils::InputType inptype =
	InputUtils::ParseInputType(file_descriptor[0]);

	InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);

	std::string outputfilename = key.GetS("output");
	if (outputfilename.empty()) {
	outputfilename = inputfilename + ".nuisance.root";
	ERR(FTL) << "No output give for set of output events! Saving to "
	<< outputfilename << std::endl;
	}

	// Make new outputfile
	TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
	outputfile->cd();


	// Get info from inputhandler
	int nevents = input->GetNEvents();
	int countwidth = (nevents / 1000);
	FitEvent* nuisevent = input->FirstNuisanceEvent();

	// Setup a TTree to save the event
	outputfile->cd();
	TTree* eventtree = new TTree("nuisance_events", "nuisance_events");

	// Add a flag that allows just splines to be saved.
	nuisevent->AddBranchesToTree(eventtree);

	// Save the spline reader
	splmerge->Write("spline_reader");

	// Setup the spline TTree
	TTree* splinetree = new TTree("spline_tree", "spline_tree");
	splmerge->AddCoefficientsToTree(splinetree);

	int lasttime = time(NULL);
	int i = 0;
	// Loop over all events and fill the TTree
	while (nuisevent) {

	// Calculate the weights for each parameter set
	splmerge->FillMergedSplines(i);

	// Save everything
	eventtree->Fill();
	splinetree->Fill();

	// Logging
	if (i % countwidth == 0) {

	std::ostringstream timestring;
	int timeelapsed = time(NULL) - lasttime;
	if (i != 0 and timeelapsed) {
	lasttime = time(NULL);

	int eventsleft = nevents - i;
	float speed = float(countwidth) / float(timeelapsed);
	float proj = (float(eventsleft) / float(speed)) / 60 / 60;
	timestring << proj << " hours remaining.";

	}
	LOG(REC) << "Saved " << i << "/" << nevents << " nuisance spline events. " << timestring.str() << std::endl;
	}

	// Iterate
	i++;
	nuisevent = input->NextNuisanceEvent();
	}

	// Save flux and close file
	outputfile->cd();
	eventtree->Write();
	splinetree->Write();

	input->GetFluxHistogram()->Write("nuisance_fluxhist");
	input->GetEventHistogram()->Write("nuisance_eventhist");

	// Close Output
	outputfile->Close();

	// Delete Inputs
	delete input;
	}



	//*************************************
	void SplineRoutines::TestSplines_1DEventScan() {
	//*************************************

	// Setup RW Engine
	if (fRW) delete fRW;
	SetupRWEngine();

	// Make a spline RW Engine too.
	FitWeight* splweight = new FitWeight("splinerwaweight");
	// std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
	std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
	TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));

	// Add Parameters
	for (size_t i = 0; i < parameterkeys.size(); i++) {
	nuiskey key = parameterkeys[i];

	std::string parname = key.GetS("name");
	std::string partype = key.GetS("type");
	double nom = key.GetD("nominal");

	parhisttemplate->SetBinContent(i + 1, nom);
	parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());

	splweight->IncludeDial( key.GetS("name"),
	kSPLINEPARAMETER, nom);
	splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );

	}
	splweight->Reconfigure();

	// Make a high resolution spline set.
	std::vector<double> nomvals = fRW->GetDialValues();
	int testres = FitPar::Config().GetParI("spline_test_resolution");

	std::vector< std::vector<double> > scanparset_vals;
	std::vector< TH1D* > scanparset_hists;

	// Loop over all params
	// Add Parameters
	for (size_t i = 0; i < parameterkeys.size(); i++) {
	nuiskey key = parameterkeys[i];

	// Get Par Name
	std::string name = key.GetS("name");

	if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
	continue;
	}

	// Push Back Scan
	double low = key.GetD("low");
	double high = key.GetD("high");
	double cur = low;
	double step = key.GetD("step");

	while (cur <= high) {

	// Make new set
	std::vector<double> newvals = nomvals;
	newvals[i] = cur;

	// Add to vects
	scanparset_vals.push_back(newvals);

	TH1D* parhist = (TH1D*)parhisttemplate->Clone();
	for (size_t j = 0; j < newvals.size(); j++) {
	parhist->SetBinContent(j + 1, newvals[j]);
	}
	scanparset_hists.push_back(parhist);


	// Move to next one
	cur += step;
	}
	}

	// Print out the parameter set to test
	for (int i = 0; i < scanparset_vals.size(); i++) {
	std::cout << "Parset " << i;
	for (int j = 0 ; j < scanparset_vals[i].size(); j++) {
	std::cout << " " << scanparset_vals[i][j];
	}
	std::cout << std::endl;
	}


	// Weight holders
	double* rawweights = new double[scanparset_vals.size()];
	double* splweights = new double[scanparset_vals.size()];
	double* difweights = new double[scanparset_vals.size()];

	int NParSets = scanparset_vals.size();

	// Loop over all event I/O
	std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
	for (size_t i = 0; i < eventkeys.size(); i++) {
	nuiskey key = eventkeys.at(i);

	// Get I/O
	std::string inputfilename = key.GetS("input");
	if (inputfilename.empty()) {
	ERR(FTL) << "No input given for set of input events!" << std::endl;
	throw;
	}

	std::string outputfilename = key.GetS("output");
	if (outputfilename.empty()) {
	outputfilename = inputfilename + ".nuisance.root";
	ERR(FTL) << "No output give for set of output events! Saving to "
	<< outputfilename << std::endl;
	}

	// Make a new input handler
	std::vector<std::string> file_descriptor =
	GeneralUtils::ParseToStr(inputfilename, ":");
	if (file_descriptor.size() != 2) {
	ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
	<< "\". expected \"FILETYPE:file.root\"" << std::endl;
	throw;
	}
	InputUtils::InputType inptype =
	InputUtils::ParseInputType(file_descriptor[0]);


	// Make handlers for input and output
	InputHandlerBase* input = InputUtils::CreateInputHandler("rawevents", inptype, file_descriptor[1]);
	InputHandlerBase* output = InputUtils::CreateInputHandler("splineevents", InputUtils::kEVSPLN_Input, outputfilename);

	// Get Base Events for each case.
	FitEvent* rawevent = input->FirstNuisanceEvent();
	FitEvent* splevent = output->FirstNuisanceEvent();


	// Setup outputfile
	std::string outputtest = outputfilename + ".splinetest.1DEventScan.root";
	TFile* outputtestfile = new TFile(outputtest.c_str(), "RECREATE");
	outputtestfile->cd();

	// Save Parameter Sets
	for (size_t i = 0; i < scanparset_hists.size(); i++) {
	scanparset_hists[i]->Write(Form("Paramater_Set_%i", (int)i));
	}

	// Save a TTree of weights and differences.
	TTree* weighttree = new TTree("weightscan", "weightscan");

	// Make a branch for each weight set
	for (size_t i = 0; i < scanparset_hists.size(); i++) {
	weighttree->Branch(Form("RawWeights_Set_%i", (int)i), &rawweights[i], Form("RawWeights_Set_%i/D", (int)i) );
	weighttree->Branch(Form("SplineWeights_Set_%i", (int)i), &splweights[i], Form("SplineWeights_Set_%i/D", (int)i) );
	weighttree->Branch(Form("DifWeights_Set_%i", (int)i), &difweights[i], Form("DifWeights_Set_%i/D", (int)i) );

	}

	// Count
	int i = 0;
	int nevents = input->GetNEvents();
	while (rawevent and splevent) {

	// Loop over 1D parameter sets.
	for (size_t j = 0; j < scanparset_vals.size(); j++) {

	// Reconfigure
	fRW->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
	fRW->Reconfigure();

	// Reconfigure spline RW
	splweight->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
	splweight->Reconfigure();

	splevent->fSplineRead->SetNeedsReconfigure(true);

	// Calc weight for both events
	rawweights[j] = fRW->CalcWeight(rawevent);
	splweights[j] = splweight->CalcWeight(splevent);
	difweights[j] = splweights[j] - rawweights[j];
	}


	if (i % 1000 == 0) {
	LOG(FIT) << "Processed " << i << "/" << nevents << std::endl;
	}

	// Fill Array
	weighttree->Fill();

	// Iterate to next event.
	i++;
	rawevent = input->NextNuisanceEvent();
	splevent = output->NextNuisanceEvent();
	}

	outputtestfile->cd();
	weighttree->Write();
	outputtestfile->Close();
	}
	}



	//*************************************
	void SplineRoutines::TestSplines_NDEventThrow() {
	//*************************************

	// Setup RW Engine
	if (fRW) delete fRW;
	SetupRWEngine();

	// Make a spline RW Engine too.
	FitWeight* splweight = new FitWeight("splinerwaweight");
	std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
	TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));

	// Add Parameters
	for (size_t i = 0; i < parameterkeys.size(); i++) {
	nuiskey key = parameterkeys[i];

	std::string parname = key.GetS("name");
	std::string partype = key.GetS("type");
	double nom = key.GetD("nominal");

	parhisttemplate->SetBinContent(i + 1, nom);
	parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());

	splweight->IncludeDial( key.GetS("name"),
	kSPLINEPARAMETER, nom);
	splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );

	}
	splweight->Reconfigure();

	// Make a high resolution spline set.
	std::vector<double> nomvals = fRW->GetDialValues();
	int testres = FitPar::Config().GetParI("spline_test_resolution");

	std::vector< std::string > scanparset_names;
	std::vector< std::vector<double> > scanparset_vals;
	std::vector< TH1D* > scanparset_hists;

	// Loop over all params
	// Add Parameters
	int nthrows = FitPar::Config().GetParI("spline_test_throws");
	for (int i = 0; i < nthrows; i++) {

	std::vector<double> newvals = nomvals;

	for (size_t j = 0; j < parameterkeys.size(); j++) {
	nuiskey key = parameterkeys[j];

	if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
	continue;
	}

	// Push Back Scan
	double low = key.GetD("low");
	double high = key.GetD("high");
	newvals[j] = gRandom->Uniform(low, high);

	}
	// Add to vects
	scanparset_vals.push_back(newvals);

	TH1D* parhist = (TH1D*)parhisttemplate->Clone();
	for (size_t j = 0; j < newvals.size(); j++) {
	parhist->SetBinContent(j + 1, newvals[j]);
	}
	scanparset_hists.push_back(parhist);
	}

	// Print out the parameter set to test
	for (int i = 0; i < scanparset_vals.size(); i++) {
	std::cout << "Parset " << i;
	for (int j = 0 ; j < scanparset_vals[i].size(); j++) {
	std::cout << " " << scanparset_vals[i][j];
	}
	std::cout << std::endl;
	}


	// Weight holders
	double* rawweights = new double[scanparset_vals.size()];
	double* splweights = new double[scanparset_vals.size()];
	double* difweights = new double[scanparset_vals.size()];

	int NParSets = scanparset_vals.size();

	// Loop over all event I/O
	std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
	for (size_t i = 0; i < eventkeys.size(); i++) {
	nuiskey key = eventkeys.at(i);

	// Get I/O
	std::string inputfilename = key.GetS("input");
	if (inputfilename.empty()) {
	ERR(FTL) << "No input given for set of input events!" << std::endl;
	throw;
	}

	std::string outputfilename = key.GetS("output");
	if (outputfilename.empty()) {
	outputfilename = inputfilename + ".nuisance.root";
	ERR(FTL) << "No output give for set of output events! Saving to "
	<< outputfilename << std::endl;
	}

	// Make a new input handler
	std::vector<std::string> file_descriptor =
	GeneralUtils::ParseToStr(inputfilename, ":");
	if (file_descriptor.size() != 2) {
	ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
	<< "\". expected \"FILETYPE:file.root\"" << std::endl;
	throw;
	}
	InputUtils::InputType inptype =
	InputUtils::ParseInputType(file_descriptor[0]);


	// Make handlers for input and output
	InputHandlerBase* input = InputUtils::CreateInputHandler("rawevents", inptype, file_descriptor[1]);
	InputHandlerBase* output = InputUtils::CreateInputHandler("splineevents", InputUtils::kEVSPLN_Input, outputfilename);

	// Get Base Events for each case.
	FitEvent* rawevent = input->FirstNuisanceEvent();
	FitEvent* splevent = output->FirstNuisanceEvent();


	// Setup outputfile
	std::string outputtest = outputfilename + ".splinetest.NDEventThrow.root";
	TFile* outputtestfile = new TFile(outputtest.c_str(), "RECREATE");
	outputtestfile->cd();

	// Save Parameter Sets
	for (size_t i = 0; i < scanparset_hists.size(); i++) {
	scanparset_hists[i]->Write(Form("Paramater_Set_%i", (int)i));
	}

	// Save a TTree of weights and differences.
	TTree* weighttree = new TTree("weightscan", "weightscan");

	// Make a branch for each weight set
	for (size_t i = 0; i < scanparset_hists.size(); i++) {
	weighttree->Branch(Form("RawWeights_Set_%i", (int)i), &rawweights[i], Form("RawWeights_Set_%i/D", (int)i) );
	weighttree->Branch(Form("SplineWeights_Set_%i", (int)i), &splweights[i], Form("SplineWeights_Set_%i/D", (int)i) );
	weighttree->Branch(Form("DifWeights_Set_%i", (int)i), &difweights[i], Form("DifWeights_Set_%i/D", (int)i) );

	}

	// Count
	int i = 0;
	int nevents = input->GetNEvents();
	while (rawevent and splevent) {

	// Loop over 1D parameter sets.
	for (size_t j = 0; j < scanparset_vals.size(); j++) {

	// Reconfigure
	fRW->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
	fRW->Reconfigure();

	// Reconfigure spline RW
	splweight->SetAllDials(&scanparset_vals[j][0], scanparset_vals[j].size());
	splweight->Reconfigure();

	splevent->fSplineRead->SetNeedsReconfigure(true);

	// Calc weight for both events
	rawweights[j] = fRW->CalcWeight(rawevent);
	splweights[j] = splweight->CalcWeight(splevent);
	difweights[j] = splweights[j] - rawweights[j];
	}


	if (i % 1000 == 0) {
	LOG(FIT) << "Processed " << i << "/" << nevents << std::endl;
	}

	// Fill Array
	weighttree->Fill();

	// Iterate to next event.
	i++;
	rawevent = input->NextNuisanceEvent();
	splevent = output->NextNuisanceEvent();
	}

	outputtestfile->cd();
	weighttree->Write();
	outputtestfile->Close();
	}
	}


	void SplineRoutines::SaveSplinePlots() {

	if (fRW) delete fRW;
	SetupRWEngine();

	// Setup the spline reader
	SplineWriter* splwrite = new SplineWriter(fRW);
	std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");

	// Add splines to splinewriter
	for (std::vector<nuiskey>::iterator iter = splinekeys.begin();
	iter != splinekeys.end(); iter++) {
	nuiskey splkey = (*iter);

	// Add Spline Info To Reader
	splwrite->AddSpline(splkey);
	}
	splwrite->SetupSplineSet();



	// Event Loop
	// Loop over all events and calculate weights for each parameter set.

	// Generate a set of nominal events
	// Method, Loop over inputs, create input handler, then create a ttree
	std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
	for (size_t i = 0; i < eventkeys.size(); i++) {
	nuiskey key = eventkeys.at(i);

	// Get I/O
	std::string inputfilename = key.GetS("input");
	if (inputfilename.empty()) {
	ERR(FTL) << "No input given for set of input events!" << std::endl;
	throw;
	}

	std::string outputfilename = key.GetS("output");
	if (outputfilename.empty()) {
	outputfilename = inputfilename + ".nuisance.root";
	ERR(FTL) << "No output give for set of output events! Saving to "
	<< outputfilename << std::endl;
	}

	// Make new outputfile
	outputfilename += ".SplinePlots.root";
	TFile* outputfile = new TFile(outputfilename.c_str(), "RECREATE");
	outputfile->cd();

	// Make a new input handler
	std::vector<std::string> file_descriptor =
	GeneralUtils::ParseToStr(inputfilename, ":");
	if (file_descriptor.size() != 2) {
	ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfilename
	<< "\". expected \"FILETYPE:file.root\"" << std::endl;
	throw;
	}
	InputUtils::InputType inptype =
	InputUtils::ParseInputType(file_descriptor[0]);

	InputHandlerBase* input = InputUtils::CreateInputHandler("eventsaver", inptype, file_descriptor[1]);

	// Get info from inputhandler
	int nevents = input->GetNEvents();
	int countwidth = (nevents / 1000);
	FitEvent* nuisevent = input->FirstNuisanceEvent();

	outputfile->cd();

	int lasttime = time(NULL);
	TCanvas* fitcanvas = NULL;

	// Loop over all events and fill the TTree
	while (nuisevent) {

	// std::cout << "Fitting event " << i << std::endl;
	// Calculate the weights for each parameter set
	splwrite->GetWeightsForEvent(nuisevent);
	splwrite->FitSplinesForEvent(fitcanvas, true);

	if (fitcanvas) {
	outputfile->cd();
	fitcanvas->Write(Form("Event_SplineCanvas_%i", (int)i));
	}

	// Logging
	if (i % countwidth == 0) {
	LOG(REC) << "Saved " << i << "/" << nevents << " nuisance spline plots. " << std::endl;
	}

	// Iterate
	i++;
	nuisevent = input->NextNuisanceEvent();
	}
	// Save flux and close file
	outputfile->cd();

	// Close Output
	outputfile->Close();

	// Delete Inputs
	delete input;
	}

	// remove Keys
	eventkeys.clear();

	}

	void SplineRoutines::TestSplines_NDLikelihoodThrow() {

	// Setup RW Engine
	if (fRW) delete fRW;
	SetupRWEngine();

	// Make a spline RW Engine too.
	FitWeight* splweight = new FitWeight("splinerwaweight");
	std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
	TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));

	// Add Parameters
	for (size_t i = 0; i < parameterkeys.size(); i++) {
	nuiskey key = parameterkeys[i];

	std::string parname = key.GetS("name");
	std::string partype = key.GetS("type");
	double nom = key.GetD("nominal");

	parhisttemplate->SetBinContent(i + 1, nom);
	parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());

	splweight->IncludeDial( key.GetS("name"),
	kSPLINEPARAMETER, nom);
	splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );

	}
	splweight->Reconfigure();

	// Make a high resolution spline set.
	std::vector<double> nomvals = fRW->GetDialValues();
	int testres = FitPar::Config().GetParI("spline_test_resolution");

	std::vector< std::string > scanparset_names;
	std::vector< std::vector<double> > scanparset_vals;
	std::vector< TH1D* > scanparset_hists;

	// Loop over all params
	// Add Parameters
	int nthrows = FitPar::Config().GetParI("spline_test_throws");
	for (int i = 0; i < nthrows; i++) {

	std::vector<double> newvals = nomvals;

	for (size_t j = 0; j < parameterkeys.size(); j++) {
	nuiskey key = parameterkeys[j];

	if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
	continue;
	}

	// Push Back Scan
	double low = key.GetD("low");
	double high = key.GetD("high");
	newvals[j] = gRandom->Uniform(low, high);

	}
	// Add to vects
	scanparset_vals.push_back(newvals);

	TH1D* parhist = (TH1D*)parhisttemplate->Clone();
	for (size_t j = 0; j < newvals.size(); j++) {
	parhist->SetBinContent(j + 1, newvals[j]);
	}
	scanparset_hists.push_back(parhist);
	}

	// Print out the parameter set to test
	for (int i = 0; i < scanparset_vals.size(); i++) {
	std::cout << "Parset " << i;
	for (int j = 0 ; j < scanparset_vals[i].size(); j++) {
	std::cout << " " << scanparset_vals[i][j];
	}
	std::cout << std::endl;
	}

	// Make a new set of Raw/Spline Sample Keys
	std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
	std::vector<nuiskey> testkeys = Config::QueryKeys("sampletest");

	std::vector<nuiskey> rawkeys;
	std::vector<nuiskey> splkeys;


	for (std::vector<nuiskey>::iterator iter = testkeys.begin();
	iter != testkeys.end(); iter++) {
	nuiskey key = (*iter);

	std::string samplename = key.GetS("name");
	std::string eventsid = key.GetS("inputid");
	nuiskey eventskey = Config::QueryLastKey("events", "id=" + eventsid);
	std::string rawfile = eventskey.GetS("input");
	std::string splfile = eventskey.GetS("output");

	nuiskey rawkeytemp = Config::CreateKey("sample");
	rawkeytemp.SetS("name", samplename);
	rawkeytemp.SetS("input", rawfile);

	nuiskey splkeytemp = Config::CreateKey("sample");
	splkeytemp.SetS("name", samplename);
	splkeytemp.SetS("input", "EVSPLN:" + splfile);

	rawkeys.push_back(rawkeytemp);
	splkeys.push_back(splkeytemp);
	}

	if (fOutputRootFile) delete fOutputRootFile;
	fOutputRootFile = new TFile(fOutputFile.c_str(), "RECREATE");

	fOutputRootFile->ls();
	// Make two new JointFCN
	JointFCN* rawfcn = new JointFCN(rawkeys, fOutputRootFile);
	JointFCN* splfcn = new JointFCN(splkeys, fOutputRootFile);

	// Create iteration tree in output file
	fOutputRootFile->cd();
	rawfcn->CreateIterationTree("raw_iterations", fRW);
	splfcn->CreateIterationTree("spl_iterations", splweight);

	// Loop over parameter sets.
	for (size_t j = 0; j < scanparset_vals.size(); j++) {

	FitBase::SetRW(fRW);
	double rawtotal = rawfcn->DoEval(&scanparset_vals[j][0]);

	FitBase::SetRW(splweight);
	double spltotal = splfcn->DoEval(&scanparset_vals[j][0]);

	LOG(FIT) << "RAW SPLINE DIF = " << rawtotal << " " << spltotal << " " << spltotal - rawtotal << std::endl;
	}

	fOutputRootFile->cd();

	rawfcn->WriteIterationTree();
	splfcn->WriteIterationTree();

	}


	void SplineRoutines::TestSplines_1DLikelihoodScan() {

	// Setup RW Engine.
	if (fRW) delete fRW;
	SetupRWEngine();

	// Setup Parameter Set.
	// Make a spline RW Engine too.
	FitWeight* splweight = new FitWeight("splinerwaweight");
	// std::vector<nuiskey> splinekeys = Config::QueryKeys("spline");
	std::vector<nuiskey> parameterkeys = Config::QueryKeys("parameter");
	TH1D* parhisttemplate = new TH1D("parhist", "parhist", parameterkeys.size(), 0.0, float(parameterkeys.size()));

	// Add Parameters
	for (size_t i = 0; i < parameterkeys.size(); i++) {
	nuiskey key = parameterkeys[i];

	std::string parname = key.GetS("name");
	std::string partype = key.GetS("type");
	double nom = key.GetD("nominal");

	parhisttemplate->SetBinContent(i + 1, nom);
	parhisttemplate->GetXaxis()->SetBinLabel(i + 1, parname.c_str());

	splweight->IncludeDial( key.GetS("name"),
	kSPLINEPARAMETER, nom);
	splweight->SetDialValue( key.GetS("name"), key.GetD("nominal") );

	}
	splweight->Reconfigure();

	// Make a high resolution spline set.
	std::vector<double> nomvals = fRW->GetDialValues();
	int testres = FitPar::Config().GetParI("spline_test_resolution");

	std::vector< std::vector<double> > scanparset_vals;
	std::vector< TH1D* > scanparset_hists;

	// Loop over all params
	// Add Parameters
	for (size_t i = 0; i < parameterkeys.size(); i++) {
	nuiskey key = parameterkeys[i];

	// Get Par Name
	std::string name = key.GetS("name");

	if (!key.Has("low") or !key.Has("high") or !key.Has("step")) {
	continue;
	}

	// Push Back Scan
	double low = key.GetD("low");
	double high = key.GetD("high");
	double cur = low;
	double step = key.GetD("step");

	while (cur <= high) {

	// Make new set
	std::vector<double> newvals = nomvals;
	newvals[i] = cur;

	// Add to vects
	scanparset_vals.push_back(newvals);

	TH1D* parhist = (TH1D*)parhisttemplate->Clone();
	for (size_t j = 0; j < newvals.size(); j++) {
	parhist->SetBinContent(j + 1, newvals[j]);
	}
	scanparset_hists.push_back(parhist);


	// Move to next one
	cur += step;
	}
	}

	// Print out the parameter set to test
	for (int i = 0; i < scanparset_vals.size(); i++) {
	std::cout << "Parset " << i;
	for (int j = 0 ; j < scanparset_vals[i].size(); j++) {
	std::cout << " " << scanparset_vals[i][j];
	}
	std::cout << std::endl;
	}

	// Make a new set of Raw/Spline Sample Keys
	std::vector<nuiskey> eventkeys = Config::QueryKeys("events");
	std::vector<nuiskey> testkeys = Config::QueryKeys("sampletest");

	std::vector<nuiskey> rawkeys;
	std::vector<nuiskey> splkeys;


	for (std::vector<nuiskey>::iterator iter = testkeys.begin();
	iter != testkeys.end(); iter++) {
	nuiskey key = (*iter);

	std::string samplename = key.GetS("name");
	std::string eventsid = key.GetS("inputid");
	nuiskey eventskey = Config::QueryLastKey("events", "id=" + eventsid);
	std::string rawfile = eventskey.GetS("input");
	std::string splfile = eventskey.GetS("output");

	nuiskey rawkeytemp = Config::CreateKey("sample");
	rawkeytemp.SetS("name", samplename);
	rawkeytemp.SetS("input", rawfile);

	nuiskey splkeytemp = Config::CreateKey("sample");
	splkeytemp.SetS("name", samplename);
	splkeytemp.SetS("input", "EVSPLN:" + splfile);

	rawkeys.push_back(rawkeytemp);
	splkeys.push_back(splkeytemp);
	}

	if (fOutputRootFile) delete fOutputRootFile;
	fOutputRootFile = new TFile(fOutputFile.c_str(), "RECREATE");

	fOutputRootFile->ls();
	// Make two new JointFCN
	JointFCN* rawfcn = new JointFCN(rawkeys, fOutputRootFile);
	JointFCN* splfcn = new JointFCN(splkeys, fOutputRootFile);

	// Create iteration tree in output file
	fOutputRootFile->cd();
	rawfcn->CreateIterationTree("raw_iterations", fRW);
	splfcn->CreateIterationTree("spl_iterations", splweight);

	// Loop over parameter sets.
	for (size_t j = 0; j < scanparset_vals.size(); j++) {

	FitBase::SetRW(fRW);
	double rawtotal = rawfcn->DoEval(&scanparset_vals[j][0]);

	FitBase::SetRW(splweight);
	double spltotal = splfcn->DoEval(&scanparset_vals[j][0]);

	LOG(FIT) << "RAW SPLINE DIF = " << rawtotal << " " << spltotal << " " << spltotal - rawtotal << std::endl;
	}

	fOutputRootFile->cd();

	rawfcn->WriteIterationTree();
	splfcn->WriteIterationTree();

	}


	/*
	MISC Functions
	*/
	//*************************************
	int SplineRoutines::GetStatus() {
	//*************************************

	return 0;
	}
	diff --git a/src/Splines/SplineWriter.cxx b/src/Splines/SplineWriter.cxx
	index 5055b477..d4f8be87 100644
	--- a/src/Splines/SplineWriter.cxx
	+++ b/src/Splines/SplineWriter.cxx
	@@ -1,765 +1,766 @@
	#include "SplineWriter.h"
	using namespace SplineUtils;

	// Spline reader should have access to every spline.
	// Should know when reconfigure is called what its limits are and adjust accordingly.
	// Then should pass it the index required in the stack as &xvals[index], and &pars[parindex]

	void SplineWriter::Write(std::string name) {
	// Create a TTree with each form and scan points in it.
	TTree* tr = new TTree(name.c_str(), name.c_str());

	// Loop over all splines and add a TString in the TTree for its inputs
	tr->Branch("Spline", &fSpline);
	tr->Branch("Type", &fType);
	tr->Branch("Form", &fForm);
	tr->Branch("Points", &fPoints);
	tr->Fill();
	tr->Write();

	delete tr;
	}

	void SplineWriter::AddWeightsToTree(TTree* tr) {
	// Add only the fitted spline coefficients to the ttree
	std::cout << "Saving Spline Weights to TTree = " << Form("SplineWeights[%d]/F", (int)fValList.size()) << std::endl;
	// sleep(1);
	tr->Branch("SplineWeights", fWeightList, Form("SplineWeights[%d]/D", (int)fValList.size()));
	}

	void SplineWriter::ReadWeightsFromTree(TTree* tr) {
	tr->SetBranchAddress("SplineWeights", fWeightList);
	}


	void SplineWriter::AddCoefficientsToTree(TTree* tr) {
	// Add only the fitted spline coefficients to the ttree
	std::cout << "Saving Spline Coeff to TTree = " << Form("SplineCoeff[%d]/F", fNCoEff) << std::endl;
	// sleep(1);
	tr->Branch("SplineCoeff", fCoEffStorer, Form("SplineCoeff[%d]/F", fNCoEff));
	}


	void SplineWriter::SetupSplineSet() {
	std::cout << "Setting up spline set" << std::endl;

	// Create the coefficients double*
	fNCoEff = 0;
	for (size_t i = 0; i < fAllSplines.size(); i++) {
	fNCoEff += fAllSplines[i].GetNPar();
	}
	fCoEffStorer = new float[fNCoEff];

	std::cout << "NCoeff = " << fNCoEff << std::endl;

	// Calculate the grid of parsets
	// Setup the list of parameter coefficients.
	std::vector<double> nomvals = fRW->GetDialValues();
	fParVect.clear();
	fSetIndex.clear();

	fParVect.push_back(nomvals);
	fSetIndex.push_back(0);
	// fWeightList.push_back(1.0);
	fValList.push_back(std::vector<double>(1, 0.0));

	// Loop over all splines.
	for (size_t i = 0; i < fAllSplines.size(); i++) {

	// For each dial loop over all points within a given position
	std::vector<double> newvals = nomvals;

	std::vector<int> pos; // = SplineUtils::GetSplitDialPositions(fRW, fSpline[i]);
	std::vector<std::string> splitnames = GeneralUtils::ParseToStr(fSpline[i], ";");

	for (size_t j = 0; j < splitnames.size(); j++) {
	int temppos = fRW->GetDialPos(splitnames[j]);
	pos.push_back(temppos);
	}

	std::vector< std::vector<double> > vals = SplineUtils::GetSplitDialPoints(fPoints[i]);

	for (size_t j = 0; j < vals.size(); j++) {

	for (size_t k = 0; k < pos.size(); k++) {
	newvals[ pos[k] ] = vals[j][k];
	}
	fParVect.push_back(newvals);
	fValList.push_back(vals[j]);
	// fWeightList.push_back(1.0);
	fSetIndex.push_back(i + 1);
	}
	}

	fWeightList = new double[fValList.size()];
	for (int i = 0; i < fValList.size(); i++) {
	fWeightList[i] = 1.0;
	}

	// Print out the parameter set
	LOG(FIT) << "Parset \| Index \| Pars --- " << std::endl;
	for (size_t i = 0; i < fSetIndex.size(); i++) {
	LOG(FIT) << " Set " << i << ". \| " << fSetIndex[i] << " \| ";
	if (LOG_LEVEL(FIT)) {
	for (size_t j = 0; j < fParVect[i].size(); j++) {
	std::cout << " " << fParVect[i][j];
	}
	std::cout << std::endl;
	}
	}
	}

	void SplineWriter::GetWeightsForEvent(FitEvent* event) {
	// Get Starting Weight
	fRW->SetAllDials(&fParVect[0][0], fParVect[0].size());
	double nomweight = fRW->CalcWeight(event);
	event->RWWeight = nomweight;

	if (fDrawSplines) {
	std::cout << "Nominal Spline Weight = " << nomweight << std::endl;
	}
	fWeightList[0] = nomweight;

	// Loop over parameter sets
	for (size_t i = 1; i < fParVect.size(); i++) {
	// Update FRW
	fRW->SetAllDials(&fParVect[i][0], fParVect[i].size());

	// Calculate a weight for event
	double weight = fRW->CalcWeight(event);


	if (weight >= 0.0 and weight < 200) {
	// Fill Weight Set
	fWeightList[i] = weight / nomweight;
	if (fDrawSplines) std::cout << "Calculating values from weight set " << i << " " << fParVect[i][0] << " = " << weight << " " << weight / nomweight << std::endl;
	} else {
	fWeightList[i] = 1.0;
	}
	}
	}

	void SplineWriter::GetWeightsForEvent(FitEvent* event, double* weights) {
	// Get Starting Weight
	fRW->SetAllDials(&fParVect[0][0], fParVect[0].size());
	double nomweight = fRW->CalcWeight(event);
	event->RWWeight = nomweight;

	if (fDrawSplines) {
	std::cout << "Nominal Spline Weight = " << nomweight << std::endl;
	}
	weights[0] = nomweight;

	// Loop over parameter sets
	for (size_t i = 1; i < fParVect.size(); i++) {
	// Update FRW
	fRW->SetAllDials(&fParVect[i][0], fParVect[i].size());

	// Calculate a weight for event
	double weight = fRW->CalcWeight(event);


	if (weight >= 0.0 and weight < 200) {
	// Fill Weight Set
	weights[i] = weight / nomweight;
	if (fDrawSplines) std::cout << "Calculating values from weight set " << i << " " << fParVect[i][0] << " = " << weight << " " << weight / nomweight << std::endl;
	} else {
	weights[i] = 1.0;
	}
	+ fWeightList[i] = weights[i];
	}
	}


	void SplineWriter::FitSplinesForEvent(double* inputweights, float* coeff) {

	int n = fAllSplines.size();
	int coeffcount = 0;

	for (int i = 0; i < n; i++) {

	// DialVals
	std::vector<std::vector<double> > dialvals;
	std::vector<double> weightvals;
	bool hasresponse = false;

	for (size_t j = 0; j < fSetIndex.size(); j++) {
	if (fSetIndex[j] != i + 1) continue;

	dialvals.push_back(fValList[j]);
	double tempw = inputweights[j];
	weightvals.push_back(tempw);

	if (tempw != 1.0) hasresponse = true;
	}

	// Perform Fit
	if (hasresponse) {
	// std::cout << "Fitting Coeff" << std::endl;
	FitCoeff(&fAllSplines[i], dialvals, weightvals, &coeff[coeffcount], fDrawSplines);
	} else {
	for (int j = 0; coeffcount + j < fNCoEff; j++) {
	// std::cout << "Setting 0.0 response " << coeffcount + i << " " << fNCoEff << std::endl;
	coeff[coeffcount + j] = 0.0;
	}
	}

	// std::cout << "Adding coeffcount" << std::endl;
	// Offset coeffcount
	coeffcount += (fAllSplines[i]).GetNPar();
	}
	// std::cout << "FitEvent" << std::endl;
	return;
	}


	void SplineWriter::FitSplinesForEvent(TCanvas * fitcanvas, bool saveplot) {

	// Loop over splines
	// int count = 0;
	int coeffcount = 0;
	int n = fAllSplines.size();

	std::vector<int> splinecounter;
	for (int i = 0; i < n; i++) {
	splinecounter.push_back(coeffcount);
	coeffcount += fAllSplines[i].GetNPar();
	}

	// #pragma omp parallel for
	for (int i = 0; i < n; i++) {

	// Store X/Y Vals
	std::vector<std::vector<double> > dialvals;
	std::vector<double> weightvals;
	bool hasresponse = false;
	int npar = (fAllSplines[i]).GetNPar();
	coeffcount = splinecounter[i];

	for (size_t j = 0; j < fSetIndex.size(); j++) {
	if (fSetIndex[j] != i + 1) continue;

	dialvals.push_back(fValList[j]);
	double tempw = fWeightList[j];
	weightvals.push_back(tempw);

	if (tempw != 1.0) hasresponse = true;
	}

	// Make a new graph and fit coeff if response
	if (hasresponse) {
	FitCoeff(&fAllSplines[i], dialvals, weightvals, &fCoEffStorer[coeffcount], fDrawSplines);
	} else {
	for (int i = 0; i < npar; i++) {
	fCoEffStorer[coeffcount + i] = 0.0;
	}
	}
	}


	// Check overrides
	// if (saveplot) {
	// coeffcount = 0;

	// // Make new canvas to save stuff into
	// if (fitcanvas) delete fitcanvas;
	// fitcanvas = new TCanvas("c1", "c1", fAllSplines.size() * 400 , 600);
	// fitcanvas->Divide(fAllSplines.size(), 1);

	// // Clear out points
	// for (size_t i = 0; i < fAllDrawnGraphs.size(); i++) {
	// if (fAllDrawnGraphs[i]) delete fAllDrawnGraphs[i];
	// if (fAllDrawnHists[i]) delete fAllDrawnHists[i];
	// }
	// fAllDrawnGraphs.clear();
	// fAllDrawnHists.clear();


	// for (size_t i = 0; i < fAllSplines.size(); i++) {

	// fitcanvas->cd(i + 1);

	// // Store X/Y Vals
	// std::vector<std::vector<double> > dialvals;
	// std::vector<double> weightvals;
	// bool hasresponse = false;
	// int npar = (fAllSplines[i]).GetNPar();

	// for (size_t j = 0; j < fSetIndex.size(); j++) {
	// if ((UInt_t)fSetIndex[j] != (UInt_t)i + 1) continue;
	// dialvals.push_back(fValList[j]);
	// weightvals.push_back(fWeightList[j] - 0.0);
	// if (fWeightList[j] != 1.0) hasresponse = true;
	// }

	// if (hasresponse) {

	// TGraph* gr = new TGraph(dialvals.size(), dialvals, weightvals);
	// fAllDrawnGraphs.push_back(gr);

	// // Get XMax Min
	// int n = xvals.size();
	// double xmin = 99999.9;
	// double xmax = -99999.9;
	// for (int j = 0; j < n; j++) {
	// if (xvals[j] > xmax) xmax = xvals[j];
	// if (xvals[j] < xmin) xmin = xvals[j];
	// }

	// TH1D* hist = new TH1D("temp", "temp", 100, xmin, xmax);
	// fAllDrawnHists.push_back(hist);

	// for (int k = 0; k < 100; k++) {
	// double xtemp = hist->GetXaxis()->GetBinCenter(k + 1);
	// fAllSplines[i].Reconfigure(xtemp);
	// double ytemp = fAllSplines[i].DoEval(&fCoEffStorer[coeffcount]);
	// hist->SetBinContent(k + 1, ytemp);
	// }

	// // gr->Draw("APL");
	// hist->SetLineColor(kRed);
	// hist->Draw("HIST C");
	// hist->SetTitle("Spline Response");
	// // hist->GetYaxis()->SetRangeUser(0.0, 3.0);

	// // gStyle->SetOptStat(0);
	// hist->SetStats(0);
	// gr->SetMarkerStyle(20);
	// gr->SetTitle("True Weight Points");
	// gr->Draw("P SAME");
	// gPad->BuildLegend(0.6, 0.6, 0.85, 0.85);
	// gPad->Update();

	// hist->SetTitle(fSpline[i].c_str());
	// hist->GetXaxis()->SetTitle("Dial Variation");
	// hist->GetYaxis()->SetTitle("Event Weight");
	// fitcanvas->Update();
	// }
	// coeffcount += npar;
	// }
	// }
	}




	// Fitting
	// ----------------------------------------------
	void SplineWriter::FitCoeff(Spline * spl, std::vector< std::vector<double> >& v, std::vector<double>& w, float * coeff, bool draw) {

	std::vector<double> x;
	std::vector<double> y;
	std::vector<double> z;
	for (size_t i = 0; i < v.size(); i++) {
	x.push_back(v[i][0]);
	if (v[i].size() > 1) y.push_back(v[i][1]);
	if (v[i].size() > 2) z.push_back(v[i][2]);
	}

	switch (spl->GetType()) {

	// Polynominal Graph Fits
	case k1DPol1:
	case k1DPol2:
	case k1DPol3:
	case k1DPol4:
	case k1DPol5:
	case k1DPol6:



	FitCoeff1DGraph(spl, v.size(), &x[0], &w[0], coeff, draw);
	break;

	case k1DTSpline3:

	GetCoeff1DTSpline3(spl, x.size(), &x[0], &w[0], coeff, draw);
	break;

	case k2DPol6:
	case k2DGaus:
	case k2DTSpline3:
	FitCoeff2DGraph(spl, v.size(), &x[0], &y[0], &w[0], coeff, draw);
	break;

	default:
	break;
	}

	// fSplineFCNs[spl] = new SplineFCN(spl, v, w);
	// fSplineFCNs[spl]->SaveAs("mysplinetest_" + spl->GetName() + ".pdf", coeff);
	// sleep(1);
	// delete fSplineFCNs[spl];

	}

	void SplineWriter::FitCoeff1DGraph(Spline * spl, int n, double * x, double * y, float * coeff, bool draw) {

	TGraph* gr = new TGraph(n, x, y);
	double xmin = 99999.9;
	double xmax = -99999.9;
	for (int i = 0; i < n; i++) {
	if (x[i] > xmax) xmax = x[i];
	if (x[i] < xmin) xmin = x[i];
	}

	double xwidth = xmax - xmin;
	xmin = xmin - xwidth * 0.01;
	xmax = xmax + xwidth * 0.01;

	// Create a new function for fitting.
	TF1* func = new TF1("f1", spl, -30.0, 30.0, spl->GetNPar());
	func->SetNpx(400);
	func->FixParameter(0, 1.0); // Fix so 1.0 at nominal

	// Run the actual spline fit
	StopTalking();


	// If linear fit with two points
	if (n == 2 and spl->GetType() == k1DPol1) {

	float m = (y[1] - y[0]) / (x[1] - x[0]);
	float c = y[0] - (0.0 - x[0]) * m;

	func->SetParameter(0, c);
	func->SetParameter(1, m);

	} else if (spl->GetType() == k1DPol1) {
	gr->Fit(func, "WQ");
	} else {
	gr->Fit(func, "FMWQ");
	}


	StartTalking();

	for (int i = 0; i < spl->GetNPar(); i++) {
	coeff[i] = func->GetParameter(i);
	}

	if (draw) {
	gr->Draw("APL");
	gPad->Update();
	gPad->SaveAs(("plot_test_" + spl->GetName() + ".pdf").c_str());
	std::cout << "Saving Graph" << std::endl;
	sleep(3);
	}

	delete func;
	delete gr;
	}

	double SplineFCN::operator()(const double * x) const {
	return DoEval(x);
	}

	double SplineFCN::DoEval(const double * x) const {

	float* fx = new float[fSpl->GetNPar()];
	for (int i = 0; i < fSpl->GetNPar(); i++) {
	fx[i] = x[i];
	}//

	double tot = 0;
	for (size_t i = 0; i < fVal.size(); i++) {

	int nonzero = 0;
	for (size_t j = 0; j < fVal[i].size(); j++) {
	fSpl->Reconfigure(fVal[i][j], j);
	// std::cout << "Reconfiguring " << fVal[i][j] << " " << j << std::endl;
	if (fVal[i][j] != 0.0) nonzero++;
	}
	if (uncorrelated and nonzero > 1) continue;

	double w = fSpl->DoEval(fx);
	double wdif = (w - fWeight[i]);// / (fWeight[i] * 0.05);
	tot += sqrt(wdif * wdif);
	}

	// delete fx;

	return tot;
	}

	void SplineFCN::UpdateWeights(std::vector<double>& w) {
	for (int i = 0; i < w.size(); i++) {
	fWeight[i] = w[i];
	}
	}

	void SplineFCN::SetCorrelated(bool state) {
	uncorrelated = !state;
	}



	void SplineFCN::SaveAs(std::string name, const float * fx) {


	if (fVal[0].size() != 2) {
	TCanvas* c1 = new TCanvas("c1", "c1", 800, 600);
	c1->Divide(2, 1);
	TH1D* histmc = new TH1D("hist", "hist", fVal.size(), 0.0, double(fVal.size()));
	TH1D* histdt = new TH1D("histdt", "histdt", fVal.size(), 0.0, double(fVal.size()));

	for (size_t i = 0; i < fVal.size(); i++) {
	for (size_t j = 0; j < fVal[i].size(); j++) {
	fSpl->Reconfigure(fVal[i][j], j);
	}

	histmc->SetBinContent( i + 1, fSpl->DoEval(fx) );
	histdt->SetBinContent( i + 1, fWeight[i] );
	}

	// histmc->Add(histdt,-1.0);

	c1->cd(1);
	histmc->SetTitle("Spline;Parameter Set;Weight");
	histmc->Draw("HIST");
	histdt->SetLineColor(kRed);
	histdt->Draw("SAME HIST");

	c1->cd(2);
	TH1D* histdif = (TH1D*) histmc->Clone();
	histdif->Add(histdt, -1.0);
	histdif->Draw("HIST");

	c1->Update();
	c1->SaveAs(name.c_str());
	delete c1;
	} else {

	TGraph2D* histmc = new TGraph2D();
	TGraph2D* histdt = new TGraph2D();
	TGraph2D* histdif = new TGraph2D();

	for (size_t i = 0; i < fVal.size(); i++) {
	for (size_t j = 0; j < fVal[i].size(); j++) {
	fSpl->Reconfigure(fVal[i][j], j);
	}

	histmc->SetPoint(histmc->GetN(), fVal[i][0], fVal[i][1], fSpl->DoEval(fx));
	histdt->SetPoint(histdt->GetN(), fVal[i][0], fVal[i][1], fWeight[i]);
	histdif->SetPoint(histdif->GetN(), fVal[i][0], fVal[i][1], fabs(fSpl->DoEval(fx) - fWeight[i]));


	}
	TCanvas* c1 = new TCanvas("c1", "c1", 800, 600);
	c1->Divide(3, 1);
	c1->cd(1);
	histmc->SetTitle(("Spline;" + fSpl->GetName()).c_str());
	histmc->Draw("COLZ");
	gPad->Update();
	c1->cd(2);
	histdt->SetTitle(("Raw;" + fSpl->GetName()).c_str());
	histdt->Draw("COLZ");
	c1->cd(3);
	histdif->SetTitle(("Dif;" + fSpl->GetName()).c_str());
	histdif->Draw("COLZ");
	gPad->SetRightMargin(0.15);
	// gPad->SetLogz(1);

	gPad->Update();
	c1->SaveAs(name.c_str());
	delete c1;
	}



	// gPad->SaveAs(name.c_str());
	}

	namespace SplineUtils {
	Spline* gSpline = NULL;
	}

	double SplineUtils::Func2DWrapper(double * x, double * p) {
	// std::cout << "2D Func Wrapper " << x[0] << " " << x[1] << std::endl;
	return (*gSpline)(x, p);
	}

	void SplineWriter::FitCoeff2DGraph(Spline * spl, int n, double * x, double * y, double * w, float * coeff, bool draw) {

	#pragma omp critical
	{
	// if (gSpline != spl) gSpline = spl;
	if (SplineUtils::gSpline != spl)
	SplineUtils::gSpline = spl;

	TF2* f2 = new TF2("f2", Func2DWrapper, -30.0, 30.0, -30.0, 30.0, spl->GetNPar());
	TGraph2D* histmc = new TGraph2D(n, x, y, w);
	f2->SetNpx(400);

	StopTalking();
	histmc->Fit(f2, "FMWQ");
	StartTalking();

	for (int i = 0; i < spl->GetNPar(); i++) {
	coeff[i] = f2->GetParameter(i);
	// std::cout << "Fit 2D Func " << i << " = " << coeff[i] << std::endl;
	}

	delete f2;
	delete histmc;
	}
	}


	void SplineWriter::FitCoeffNDGraph(Spline * spl, std::vector< std::vector<double> >& v, std::vector<double>& w, float * coeff, bool draw) {

	// TGraph2D* gr = new TGraph2D(n, x, y);
	// double xmin = 99999.9;
	// double xmax = -99999.9;
	// double ymin = 99999.9;
	// double ymax = -99999.9;
	// for (int i = 0; i < n; i++) {
	// if (x[i] > xmax) xmax = x[i];
	// if (x[i] < xmin) xmin = x[i];
	// if (y[i] > ymax) ymax = y[i];
	// if (y[i] < ymin) ymin = y[i];
	// }
	// double xwidth = xmax - xmin;
	// xmin = xmin - xwidth * 0.01;
	// xmax = xmax + xwidth * 0.01;
	// double ywidth = ymax - ymin;
	// ymin = ymin - ywidth * 0.01;
	// ymax = ymax + ywidth * 0.01;


	// StopTalking();

	// Build if not already in writer
	// if (fSplineFCNs.find(spl) != fSplineFCNs.end()) {
	// fSplineFCNs.erase(spl);
	// }
	// fSplineFCNs[spl] = new SplineFCN(spl, v, w);


	if (fSplineFunctors.find(spl) != fSplineFunctors.end()) {
	delete fSplineFunctors[spl];
	fSplineFunctors.erase(spl);
	}

	if (fSplineFCNs.find(spl) != fSplineFCNs.end()) {
	delete fSplineFCNs[spl];
	fSplineFCNs.erase(spl);
	}

	if (fSplineMinimizers.find(spl) != fSplineMinimizers.end()) {
	delete fSplineMinimizers[spl];
	fSplineMinimizers.erase(spl);
	}



	if (fSplineMinimizers.find(spl) == fSplineMinimizers.end()) {
	std::cout << "Building new ND minimizer for " << spl << std::endl;
	fSplineFCNs[spl] = new SplineFCN(spl, v, w);

	// fSplineFCNs[spl] = new SplineFCN(spl, v, w);
	fSplineFunctors[spl] = new ROOT::Math::Functor( *fSplineFCNs[spl], spl->GetNPar() );
	// fitclass = "GSLSimAn"; fittype = "";
	fSplineMinimizers[spl] = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Combined");
	fSplineMinimizers[spl]->SetMaxFunctionCalls(1E8);
	fSplineMinimizers[spl]->SetMaxIterations(1E8);
	fSplineMinimizers[spl]->SetTolerance(1.E-6);
	fSplineMinimizers[spl]->SetStrategy(2);

	for (int j = 0; j < spl->GetNPar(); j++) {
	if (j != 0) {
	fSplineMinimizers[spl]->SetVariable(j, Form("Par_%i", j), 0.1, 0.1);
	// fSplineMinimizers[spl]->SetParameter(j, Form("Par_%i", j), 0.1, 0.1, -100.0, 100.0 );
	} else {
	fSplineMinimizers[spl]->SetVariable(j, Form("Par_%i", j), 1.0, 0.1);
	// fSplineMinimizers[spl]->SetParameter(j, Form("Par_%i", j), 1.0, 0.1, -100.0, 100.0 );
	// fSplineMinimizers[spl]->FixVariable(j);
	}
	}
	// fSplineMinimizers[spl]->SetFunction(*fSplineFunctors[spl]);
	sleep(1);
	}
	// Create a new function for fitting.
	// StopTalking();

	// Update FCN
	fSplineFCNs[spl]->UpdateWeights(w);
	// fSplineMinimizers[spl]->SetDefaultOptions();
	// fSplineMinimizers[spl]->Clear();
	for (int j = 0; j < spl->GetNPar(); j++) {
	if (j != 0) {
	fSplineMinimizers[spl]->SetVariableValue(j, 0.1);
	// fSplineMinimizers[spl]->SetParameter(j, Form("Par_%i", j), 0.1, 0.1, -100.0, 100.0 );

	// fSplineMinimizers[spl]->SetVariableValue(j, 0.1);
	}
	}
	// fSplineFCNs[spl]->SetCorrelated(false);
	// fSplineFunctors[spl] = new ROOT::Math::Functor( *fSplineFCNs[spl], spl->GetNPar() );
	// fSplineMinimizers[spl]->SetFunction(*fSplineFunctors[spl]);
	// fSplineMinimizers[spl]->Minimize();

	fSplineFCNs[spl]->SetCorrelated(true);
	delete fSplineFunctors[spl];
	fSplineFunctors[spl] = new ROOT::Math::Functor( *fSplineFCNs[spl], spl->GetNPar() );
	fSplineMinimizers[spl]->SetFunction(*fSplineFunctors[spl]);
	// ((TFitterMinuit*)fSplineMinimizers[spl])->CreateMinimizer(TFitterMinuit::kMigrad);
	fSplineMinimizers[spl]->Minimize();
	fSplineMinimizers[spl]->Minimize();

	// ((TFitterMinuit*)fSplineMinimizers[spl])->CreateMinimizer(TFitterMinuit::kMigrad);
	// fSplineMinimizers[spl]->Minimize();

	// delete minimizer;
	// StartTalking();

	// Now Get Parameters
	const double *values = fSplineMinimizers[spl]->X();
	for (int i = 0; i < spl->GetNPar(); i++) {

	std::cout << "Updated Coeff " << i << " " << values[i] << std::endl;
	coeff[i] = values[i];
	}

	// Save a sample
	fSplineFCNs[spl]->SaveAs("mysplinetest_" + spl->GetName() + ".pdf", coeff);
	sleep(1);

	// delete values;
	// delete minimizer;
	}



	// Spline extraction Functions
	void SplineWriter::GetCoeff1DTSpline3(Spline * spl, int n, double * x, double * y, float * coeff, bool draw) {

	StopTalking();
	TSpline3 temp_spline = TSpline3("temp_spline", x, y, n);
	StartTalking();

	for (size_t i = 0; i < (UInt_t)n; i++) {

	double a, b, c, d, e;
	temp_spline.GetCoeff(i, a, b, c, d, e);

	coeff[i * 4] = y[i];
	coeff[i * 4 + 1] = c;
	coeff[i * 4 + 2] = d;
	coeff[i * 4 + 3] = e;

	// std::cout << "Setting Spline Coeff Set " << i << " to " << y[i] << " " << c << " " << d << " " << e << std::endl;
	}

	if (draw) {
	TGraph* gr = new TGraph(n, x, y);
	temp_spline.Draw("CA");
	gr->Draw("PL SAME");
	gPad->Update();
	gPad->SaveAs(("plot_test_" + spl->GetName() + ".pdf").c_str());

	delete gr;
	}

	return;
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Sep 30, 5:44 AM (7 h, 16 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6566327
Default Alt Text: (73 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions