No OneTemporary
Actions

Size

125 KB

Subscribers

None

View Options

	diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx
	index 6ce228b..85ed5a1 100755
	--- a/src/FCN/JointFCN.cxx
	+++ b/src/FCN/JointFCN.cxx
	@@ -1,1123 +1,1154 @@
	#include "JointFCN.h"
	#include "FitUtils.h"
	#include <stdio.h>

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

	fOutputDir = gDirectory;
	if (outfile)
	Config::Get().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)
	Config::Get().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) {
	//***************************************************

	NUIS_LOG(FIT, " Creating new iteration container! ");
	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
	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++) {
	fNameValues.push_back(dials[i]);
	fCurrentValues.push_back(0.0);
	}
	fNDials = dials.size();
	fDialVals = new double[fNDials];

	// Set IterationTree Flag
	fIterationTree = true;
	}

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

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

	//***************************************************
	void JointFCN::WriteIterationTree() {
	//***************************************************
	NUIS_LOG(FIT, "Writing iteration tree");

	// 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");
	for (size_t 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++) {
	std::vector<double> itervals = fIterationValues[i];

	// Fill iteration state
	count = fIterationCount[i];
	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++) {
	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++) {
	fCurrentValues[count++] = double(fDialVals[i]);
	}

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

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

	+ double *par_vals = new double[fNPars];
	+
	+ for (int i = 0; i < fNPars; ++i) {
	+ if (fMirroredParams.count(i)) {
	+ if (!fMirroredParams[i].mirror_above &&
	+ (x[i] < fMirroredParams[i].mirror_value)) {
	+ double xabove = fMirroredParams[i].mirror_value - x[i];
	+ par_vals[i] = fMirroredParams[i].mirror_value + xabove;
	+ std::cout << "\t--Parameter " << i << " mirrored from " << x[i]
	+ << " -> " << par_vals[i] << std::endl;
	+ } else if (fMirroredParams[i].mirror_above &&
	+ (x[i] >= fMirroredParams[i].mirror_value)) {
	+ double xabove = x[i] - fMirroredParams[i].mirror_value;
	+ par_vals[i] = fMirroredParams[i].mirror_value - xabove;
	+ std::cout << "\t--Parameter " << i << " mirrored from " << x[i]
	+ << " -> " << par_vals[i] << std::endl;
	+ } else {
	+ par_vals[i] = x[i];
	+ }
	+ } else {
	+ par_vals[i] = x[i];
	+ }
	+ }
	+
	// WEIGHT ENGINE
	- fDialChanged = FitBase::GetRW()->HasRWDialChanged(x);
	- FitBase::GetRW()->UpdateWeightEngine(x);
	+ fDialChanged = FitBase::GetRW()->HasRWDialChanged(par_vals);
	+ FitBase::GetRW()->UpdateWeightEngine(par_vals);
	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
	- NUIS_LOG(FIT, "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood);
	+ NUIS_LOG(FIT,
	+ "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood);

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

	+ delete[] par_vals;
	+
	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;
	}
	return totaldof;
	}

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

	NUIS_LOG(MIN, std::left << std::setw(43) << "Getting likelihoods..."
	- << " : "
	- << "-2logL");
	+ << " : "
	+ << "-2logL");

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

	NUIS_LOG(MIN, "-> " << std::left << std::setw(40) << exp->GetName() << " : "
	- << newlike << "/" << ndof);
	+ << newlike << "/" << ndof);

	// 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) {
	fSampleLikes[count] = fLikelihood;
	}

	return like;
	};

	void JointFCN::LoadSamples(std::vector<nuiskey> samplekeys) {
	NUIS_LOG(MIN, "Loading Samples : " << samplekeys.size());
	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 = "";

	NUIS_LOG(MIN, "Loading Sample : " << samplename);

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

	if (!NewLoadedSample) {
	NUIS_ERR(FTL, "Could not load sample provided: " << samplename);
	NUIS_ERR(FTL, "Check spelling with that in src/FCN/SampleList.cxx");
	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);
	NUIS_LOG(REC, "Starting Reconfigure iter. " << this->fCurIter);
	- // std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled <<
	- // std::endl;
	- // Event Manager Reconf
	+ // std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled
	+ // << std::endl; Event Manager Reconf
	if (fUsingEventManager) {
	if (!fullconfig && 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;
	NUIS_LOG(MIN, "Finished Reconfigure iter. " << fCurIter << " in "
	- << time(NULL) - starttime << "s");
	+ << time(NULL) - starttime << "s");

	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
	NUIS_LOG(REC, "Event Manager Reconfigure");
	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.
	+ // 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 / 10;

	// Start event loop iterating until we get a NULL pointer.
	while (curevent) {
	// Get Event Weight
	// The reweighting weight
	curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
	// The Custom weight and reweight
	curevent->Weight =
	curevent->RWWeight * curevent->InputWeight * curevent->CustomWeight;

	if (LOGGING(REC)) {
	if (countwidth && (i % countwidth == 0)) {
	NUIS_LOG(REC, curinput->GetName()
	- << " : Processed " << i << " events. [M, W] = ["
	- << curevent->Mode << ", " << curevent->Weight << "]");
	+ << " : Processed " << i << " events. [M, W] = ["
	+ << curevent->Mode << ", " << curevent->Weight
	+ << "]");
	}
	}

	// 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 && 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 && savesignal && 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++)
	+ // 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.
	NUIS_LOG(REC, "Filled " << fillcount << " signal events.");
	if (savesignal) {
	int mem = ( // sizeof(fSignalEventBoxes) +
	// fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
	sizeof(MeasurementVariableBox1D) * fillcount) *
	1E-6;
	- NUIS_LOG(REC, " -> Saved " << fillcount << " signal boxes for faster access. (~"
	- << mem << " MB)");
	+ NUIS_LOG(REC, " -> Saved " << fillcount
	+ << " signal boxes for faster access. (~" << mem
	+ << " MB)");
	if (fIsAllSplines and !fSignalEventSplines.empty()) {
	int splmem = sizeof(float) * fSignalEventSplines.size() *
	fSignalEventSplines[0].size() * 1E-6;
	- NUIS_LOG(REC, " -> Saved " << fillcount << " " << fSignalEventSplines.size()
	- << " spline sets into memory. (~" << splmem
	- << " MB)");
	+ NUIS_LOG(REC, " -> Saved "
	+ << fillcount << " " << fSignalEventSplines.size()
	+ << " spline sets into memory. (~" << splmem << " MB)");
	}
	}

	NUIS_LOG(REC,
	- "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart);
	+ "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart);

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

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

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

	NUIS_LOG(FIT, " -> Doing FAST using manager");
	// 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()) {
	NUIS_ERR(WRN, "Signal Flags Empty! Using normal manager.");
	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;
	// This is just the total number of events
	// int nevents = fSignalEventFlags.size();
	// This is the number of events that are signal
	int nevents = fSignalEventBoxes.size();
	int countwidth = nevents / 10;

	// If All Splines tell splines they need a reconfigure.
	std::vector<InputHandlerBase *>::iterator inp_iter = fInputList.begin();
	if (fIsAllSplines) {
	NUIS_LOG(REC, "All Spline Inputs so using fast spline loop.");
	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;

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

	// Loop over the events in each input
	for (int i = 0; i < curinput->GetNEvents(); i++) {
	double rwweight = 0.0;

	// If the event is a signal event
	if (fSignalEventFlags[sigcount]) {
	// Get Event Info
	if (!fIsAllSplines) {
	if (fFillNuisanceEvent) {
	curevent = 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 * curevent->CustomWeight;
	rwweight = curevent->Weight;

	coreeventweights[splinecount] = rwweight;
	if (countwidth && ((splinecount % countwidth) == 0)) {
	NUIS_LOG(REC, curinput->GetName()
	- << " : Processed " << i
	- << " events. W = " << curevent->Weight << std::endl);
	+ << " : Processed " << i << " events. W = "
	+ << curevent->Weight << std::endl);
	}

	// #pragma omp atomic
	splinecount++;
	}

	// #pragma omp atomic
	sigcount++;
	}
	}

	NUIS_LOG(SAM, "Processed event weights.");

	// #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) {
	NUIS_LOG(REC, "Filled " << ispline << " sample weights.");
	}

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

	NUIS_LOG(SAM, "Filled sample distributions.");

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

	// Cleanup coreeventweights
	delete coreeventweights;

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

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

	// Save a likelihood/ndof plot
	NUIS_LOG(MIN, "Writing likelihood plot...");
	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);
	}
	if (likes.size()) {
	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 (int 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
	NUIS_LOG(MIN, "Writing each of the data classes...");
	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) {
	//***************************************************

	NUIS_LOG(MIN, "Setting fake data from " << fakeinput);
	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;
	NUIS_LOG(MIN, "Likelihoods : ");

	// 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
	NUIS_LOG(MIN, "-> " << std::left << std::setw(40) << exp->GetName() << " : "
	- << singlelike);
	+ << singlelike);
	}

	// 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
	- NUIS_LOG(MIN, "-> " << std::left << std::setw(40) << pull->GetName() << " : "
	- << singlelike);
	+ NUIS_LOG(MIN, "-> " << std::left << std::setw(40) << pull->GetName()
	+ << " : " << singlelike);
	}

	// 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 1f47f28..895bb9c 100755
	--- a/src/FCN/JointFCN.h
	+++ b/src/FCN/JointFCN.h
	@@ -1,180 +1,196 @@
	#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:

	//! 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
	//! 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; }

	//! Return list of pointers to all the pulls
	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<std::string> GetAllNames();
	std::vector<double> GetAllLikelihoods();
	std::vector<int> GetAllNDOF();

	+
	+ void SetVariableMirrored(int ipar, double mirror_value, bool mirror_above){
	+ mirror_param mir_par;
	+ mir_par.mirror_value = mirror_value;
	+ mir_par.mirror_above = mirror_above;
	+ fMirroredParams[ipar] = mir_par;
	+ }
	+ void SetNParams(int npar){
	+ fNPars = npar;
	+ }
	+
	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;


	-
	-
	+ struct mirror_param {
	+ double mirror_value;
	+ bool mirror_above;
	+ };
	+ std::map<int, mirror_param> fMirroredParams;
	+ //the number of pars added to the minimizer, should be the same as fNDials
	+ int fNPars;
	};

	/! @} /
	#endif // _JOINT_FCN_H_
	diff --git a/src/Routines/ComparisonRoutines.cxx b/src/Routines/ComparisonRoutines.cxx
	index 78807a2..a651d24 100755
	--- a/src/Routines/ComparisonRoutines.cxx
	+++ b/src/Routines/ComparisonRoutines.cxx
	@@ -1,526 +1,637 @@
	// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret

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

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

	fOutputFile = "";
	fOutputRootFile = NULL;

	fStrategy = "Compare";

	fRoutines.clear();

	fCardFile = "";

	fFakeDataInput = "";

	fSampleFCN = NULL;

	fAllowedRoutines = ("Compare");
	};

	//*************************************
	ComparisonRoutines::~ComparisonRoutines() {
	//*************************************
	delete fOutputRootFile;
	};

	/*
	Input Functions
	*/
	//*************************************
	ComparisonRoutines::ComparisonRoutines(int argc, char *argv[]) {
	//*************************************

	// Initialise Defaults
	Init();
	nuisconfig configuration = Config::Get();

	// Default containers
	std::string cardfile = "";
	std::string maxevents = "-1";
	std::string skipevents = "0";
	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, "-d", fFakeDataInput, false, false);
	ParserUtils::ParseArgument(args, "-i", xmlcmds);
	ParserUtils::ParseArgument(args, "-q", configargs);
	ParserUtils::ParseArgument(args, "-s", skipevents);
	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()) {
	NUIS_ABORT("No input supplied!");
	}

	if (fOutputFile.empty() and !fCardFile.empty()) {
	fOutputFile = fCardFile + ".root";
	NUIS_ERR(WRN, "No output supplied so saving it to: " << fOutputFile);

	} else if (fOutputFile.empty()) {
	NUIS_ABORT("No output file or cardfile supplied!");
	}

	// 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.Set("cardfile", fCardFile);
	if (!fOutputFile.empty())
	fCompKey.Set("outputfile", fOutputFile);
	if (!fStrategy.empty())
	fCompKey.Set("strategy", fStrategy);

	// Load XML Cardfile
	configuration.LoadSettings(fCompKey.GetS("cardfile"), "");

	// Add Config Args
	for (size_t i = 0; i < configargs.size(); i++) {
	configuration.OverrideConfig(configargs[i]);
	}
	if (maxevents.compare("-1")) {
	configuration.OverrideConfig("MAXEVENTS=" + maxevents);
	}
	configuration.OverrideConfig("NSKIPEVENTS=" + skipevents);

	// Finish configuration XML
	configuration.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml");

	// Add Error Verbo Lines
	verbocount += Config::GetParI("VERBOSITY");
	errorcount += Config::GetParI("ERROR");
	bool trace = Config::GetParB("TRACE");
	std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
	std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
	SETVERBOSITY(verbocount);
	SETTRACE(trace);

	// Comparison Setup ========================================

	// Proper Setup
	fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
	SetupComparisonsFromXML();

	SetupRWEngine();
	SetupFCN();

	return;
	};

	//*************************************
	void ComparisonRoutines::SetupComparisonsFromXML() {
	//*************************************

	NUIS_LOG(FIT, "Setting up nuiscomp");

	// Setup Parameters ------------------------------------------
	std::vector<nuiskey> parkeys = Config::QueryKeys("parameter");
	if (!parkeys.empty()) {
	NUIS_LOG(FIT, "Number of parameters : " << parkeys.size());
	}

	for (size_t i = 0; i < parkeys.size(); i++) {
	nuiskey key = parkeys.at(i);

	// Check for type,name,nom
	if (!key.Has("type")) {
	NUIS_ERR(FTL, "No type given for parameter " << i);
	NUIS_ABORT("type='PARAMETER_TYPE'");
	} else if (!key.Has("name")) {
	NUIS_ERR(FTL, "No name given for parameter " << i);
	NUIS_ABORT("name='SAMPLE_NAME'");
	} else if (!key.Has("nominal")) {
	NUIS_ERR(FTL, "No nominal given for parameter " << i);
	NUIS_ABORT("nominal='NOMINAL_VALUE'");
	}

	// 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");

	// Check for incomplete limtis
	int limdef =
	((int)key.Has("low") + (int)key.Has("high") + (int)key.Has("step"));

	if (limdef > 0 and limdef < 3) {
	NUIS_ERR(FTL, "Incomplete limit set given for parameter : " << parname);
	NUIS_ABORT(
	"Requires: low='LOWER_LIMIT' high='UPPER_LIMIT' step='STEP_SIZE' ");
	}

	// Extra limits
	if (key.Has("low")) {

	parlow = key.GetD("low");
	parhigh = key.GetD("high");
	parstep = key.GetD("step");

	NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom
	<< " : " << parlow << " < p < " << parhigh << " : "
	<< parstate);
	} else {
	NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom
	<< " : " << parstate);
	}

	+ bool ismirr = false;
	+ if (key.Has("mirror_point")) {
	+ ismirr = true;
	+ mirror_param mir;
	+ mir.mirror_value = key.GetD("mirror_point");
	+ mir.mirror_above = key.GetB("mirror_above");
	+ fMirroredParams[parname] = mir;
	+ NUIS_LOG(FIT,
	+ "\t\t" << parname << " is mirrored at " << mir.mirror_value
	+ << " "
	+ << (mir.mirror_above ? "from above" : "from below"));
	+ }
	+
	// Convert if required
	if (parstate.find("ABS") != std::string::npos) {
	+ if (ismirr) {
	+ NUIS_ABORT("Cannot mirror parameters with ABS state!");
	+ }
	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) {
	+ if (ismirr) {
	+ NUIS_ABORT("Cannot mirror parameters with FRAC state!");
	+ }
	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);
	fCurVals[parname] = parnom;
	fStateVals[parname] = parstate;
	}

	// Setup Samples ----------------------------------------------
	std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
	if (!samplekeys.empty()) {
	NUIS_LOG(FIT, "Number of samples : " << samplekeys.size());
	}

	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
	NUIS_LOG(FIT, "Read Sample " << i << ". : " << samplename << " ("
	<< sampletype << ") [Norm=" << samplenorm
	<< "]" << std::endl
	<< " -> input='"
	<< samplefile << "'");

	// If FREE add to parameters otherwise continue
	if (sampletype.find("FREE") == std::string::npos) {
	if (samplenorm != 1.0) {
	NUIS_ERR(FTL, "You provided a sample normalisation but did not specify "
	"that the sample is free");
	NUIS_ABORT("Change so sample contains type=\"FREE\" and re-run");
	}
	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;
	}

	// Setup Fake Parameters -----------------------------
	std::vector<nuiskey> fakekeys = Config::QueryKeys("fakeparameter");
	if (!fakekeys.empty()) {
	NUIS_LOG(FIT, "Number of fake parameters : " << fakekeys.size());
	}

	for (size_t i = 0; i < fakekeys.size(); i++) {
	nuiskey key = fakekeys.at(i);

	// Check for type,name,nom
	if (!key.Has("name")) {
	NUIS_ABORT("No name given for fakeparameter " << i);

	} else if (!key.Has("nominal")) {
	NUIS_ABORT("No nominal given for fakeparameter " << i);
	}

	// Get Inputs
	std::string parname = key.GetS("name");
	double parnom = key.GetD("nominal");

	// Push into vectors
	fFakeVals[parname] = parnom;
	}
	}

	//*************************************
	void ComparisonRoutines::SetupRWEngine() {
	//*************************************

	NUIS_LOG(FIT, "Setting up FitWeight Engine");
	for (UInt_t i = 0; i < fParams.size(); i++) {
	std::string name = fParams[i];
	FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name));
	}

	return;
	}

	//*************************************
	void ComparisonRoutines::SetupFCN() {
	//*************************************

	NUIS_LOG(FIT, "Building the SampleFCN");
	if (fSampleFCN)
	delete fSampleFCN;
	Config::Get().out = fOutputRootFile;
	fOutputRootFile->cd();
	fSampleFCN = new JointFCN(fOutputRootFile);
	SetFakeData();

	return;
	}

	//*************************************
	void ComparisonRoutines::SetFakeData() {
	//*************************************

	if (fFakeDataInput.empty())
	return;

	if (fFakeDataInput.compare("MC") == 0) {
	NUIS_LOG(FIT, "Setting fake data from MC starting prediction.");
	UpdateRWEngine(fFakeVals);

	FitBase::GetRW()->Reconfigure();
	fSampleFCN->ReconfigureAllEvents();
	fSampleFCN->SetFakeData("MC");

	- UpdateRWEngine(fCurVals);
	+ std::map<std::string, double> CurVals_wmirr;
	+ for (size_t i = 0; i < fParams.size(); ++i) {
	+ std::string const &pname = fParams[i];
	+ if (fMirroredParams.count(pname)) {
	+ if (!fMirroredParams[pname].mirror_above &&
	+ (fCurVals[pname] < fMirroredParams[pname].mirror_value)) {
	+ double xabove = fMirroredParams[pname].mirror_value - fCurVals[pname];
	+ CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value + xabove;
	+ std::cout << "\t--Parameter " << pname << " mirrored from "
	+ << fCurVals[pname] << " -> " << CurVals_wmirr[pname]
	+ << std::endl;
	+ } else if (fMirroredParams[pname].mirror_above &&
	+ (fCurVals[pname] >= fMirroredParams[pname].mirror_value)) {
	+ double xabove = fCurVals[pname] - fMirroredParams[pname].mirror_value;
	+ CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value - xabove;
	+ std::cout << "\t--Parameter " << pname << " mirrored from "
	+ << fCurVals[pname] << " -> " << CurVals_wmirr[pname]
	+ << std::endl;
	+ } else {
	+ CurVals_wmirr[pname] = fCurVals[pname];
	+ }
	+ } else {
	+ CurVals_wmirr[pname] = fCurVals[pname];
	+ }
	+ }
	+
	+ UpdateRWEngine(CurVals_wmirr);

	NUIS_LOG(FIT, "Set all data to fake MC predictions.");
	} else {
	NUIS_LOG(FIT, "Setting fake data from: " << fFakeDataInput);
	fSampleFCN->SetFakeData(fFakeDataInput);
	}

	return;
	}

	/*
	Fitting Functions
	*/
	//*************************************
	void ComparisonRoutines::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 ComparisonRoutines::Run() {
	//*************************************

	NUIS_LOG(FIT, "Running ComparisonRoutines : " << fStrategy);

	if (FitPar::Config().GetParB("save_nominal")) {
	SaveNominal();
	}

	// Parse given routines
	fRoutines = GeneralUtils::ParseToStr(fStrategy, ",");
	if (fRoutines.empty()) {
	NUIS_ABORT("Trying to run ComparisonRoutines with no routines given!");
	}

	for (UInt_t i = 0; i < fRoutines.size(); i++) {
	std::string routine = fRoutines.at(i);

	NUIS_LOG(FIT, "Routine: " << routine);
	if (!routine.compare("Compare")) {
	- UpdateRWEngine(fCurVals);
	+
	+ std::map<std::string, double> CurVals_wmirr;
	+ for (size_t i = 0; i < fParams.size(); ++i) {
	+ std::string const &pname = fParams[i];
	+ if (fMirroredParams.count(pname)) {
	+ std::cout << "MA? " << fMirroredParams[pname].mirror_above
	+ << ", OVal: " << fCurVals[pname]
	+ << ", MPoint: " << fMirroredParams[pname].mirror_value
	+ << std::endl;
	+ if (!fMirroredParams[pname].mirror_above &&
	+ (fCurVals[pname] < fMirroredParams[pname].mirror_value)) {
	+ double xabove =
	+ fMirroredParams[pname].mirror_value - fCurVals[pname];
	+ CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value + xabove;
	+ std::cout << "\t--Parameter " << pname << " mirrored from "
	+ << fCurVals[pname] << " -> " << CurVals_wmirr[pname]
	+ << std::endl;
	+ } else if (fMirroredParams[pname].mirror_above &&
	+ (fCurVals[pname] >= fMirroredParams[pname].mirror_value)) {
	+ double xabove =
	+ fCurVals[pname] - fMirroredParams[pname].mirror_value;
	+ CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value - xabove;
	+ std::cout << "\t--Parameter " << pname << " mirrored from "
	+ << fCurVals[pname] << " -> " << CurVals_wmirr[pname]
	+ << std::endl;
	+ } else {
	+ CurVals_wmirr[pname] = fCurVals[pname];
	+ }
	+ } else {
	+ CurVals_wmirr[pname] = fCurVals[pname];
	+ }
	+ std::cout << "~~~~~~~" << pname << " : " << fCurVals[pname] << " -> "
	+ << CurVals_wmirr[pname] << std::endl;
	+ }
	+
	+ UpdateRWEngine(CurVals_wmirr);
	GenerateComparison();
	PrintState();
	SaveCurrentState();
	}
	}

	return;
	}

	//*************************************
	void ComparisonRoutines::GenerateComparison() {
	//*************************************
	NUIS_LOG(FIT, "Generating Comparison.");
	// Main Event Loop from event Manager
	fSampleFCN->ReconfigureAllEvents();
	return;
	}

	//*************************************
	void ComparisonRoutines::PrintState() {
	//*************************************
	NUIS_LOG(FIT, "------------");

	// Count max size
	int maxcount = 0;
	for (UInt_t i = 0; i < fParams.size(); i++) {
	maxcount = max(int(fParams[i].size()), maxcount);
	}

	// Header
	NUIS_LOG(FIT, " # " << left << setw(maxcount) << "Parameter "
	<< " = " << setw(10) << "Value"
	<< " +- " << setw(10) << "Error"
	<< " " << setw(8) << "(Units)"
	<< " " << setw(10) << "Conv. Val"
	<< " +- " << setw(10) << "Conv. Err"
	<< " " << setw(8) << "(Units)");

	+ std::map<std::string, double> CurVals_wmirr;
	+ for (size_t i = 0; i < fParams.size(); ++i) {
	+ std::string const &pname = fParams[i];
	+ if (fMirroredParams.count(pname)) {
	+ if (!fMirroredParams[pname].mirror_above &&
	+ (fCurVals[pname] < fMirroredParams[pname].mirror_value)) {
	+ double xabove = fMirroredParams[pname].mirror_value - fCurVals[pname];
	+ CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value + xabove;
	+ std::cout << "\t--Parameter " << pname << " mirrored from "
	+ << fCurVals[pname] << " -> " << CurVals_wmirr[pname]
	+ << std::endl;
	+ } else if (fMirroredParams[pname].mirror_above &&
	+ (fCurVals[pname] >= fMirroredParams[pname].mirror_value)) {
	+ double xabove = fCurVals[pname] - fMirroredParams[pname].mirror_value;
	+ CurVals_wmirr[pname] = fMirroredParams[pname].mirror_value - xabove;
	+ std::cout << "\t--Parameter " << pname << " mirrored from "
	+ << fCurVals[pname] << " -> " << CurVals_wmirr[pname]
	+ << std::endl;
	+ } else {
	+ CurVals_wmirr[pname] = fCurVals[pname];
	+ }
	+ } else {
	+ CurVals_wmirr[pname] = fCurVals[pname];
	+ }
	+ }
	+
	// Parameters
	for (UInt_t i = 0; i < fParams.size(); i++) {
	std::string syst = fParams.at(i);

	std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
	std::string curunits = "(sig.)";
	double curval = fCurVals[syst];
	double curerr = 0.0;

	if (fStateVals[syst].find("ABS") != std::string::npos) {
	curval = FitBase::RWSigmaToAbs(typestr, syst, curval);
	curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
	FitBase::RWSigmaToAbs(typestr, syst, 0.0));
	curunits = "(Abs.)";
	} else if (fStateVals[syst].find("FRAC") != std::string::npos) {
	curval = FitBase::RWSigmaToFrac(typestr, syst, curval);
	curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) -
	FitBase::RWSigmaToFrac(typestr, syst, 0.0));
	curunits = "(Frac)";
	}

	std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")";
	double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
	double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
	FitBase::RWSigmaToAbs(typestr, syst, 0.0));

	std::ostringstream curparstring;

	curparstring << " " << setw(3) << left << i << ". " << setw(maxcount)
	<< syst << " = " << setw(10) << curval << " +- " << setw(10)
	<< curerr << " " << setw(8) << curunits << " " << setw(10)
	<< convval << " +- " << setw(10) << converr << " " << setw(8)
	<< convunits;
	-
	NUIS_LOG(FIT, curparstring.str());
	+
	+ if (fMirroredParams.count(syst)) {
	+ NUIS_LOG(FIT, "\t\t--> Mirrored at " << fMirroredParams[syst].mirror_value
	+ << " to effective value "
	+ << CurVals_wmirr[syst]);
	+ }
	}

	NUIS_LOG(FIT, "------------");
	double like = fSampleFCN->GetLikelihood();
	NUIS_LOG(FIT,
	std::left << std::setw(46) << "Likelihood for JointFCN: " << like);
	NUIS_LOG(FIT, "------------");
	}

	/*
	Write Functions
	*/
	//*************************************
	void ComparisonRoutines::SaveCurrentState(std::string subdir) {
	//*************************************

	NUIS_LOG(FIT, "Saving current full FCN predictions");

	// Setup DIRS
	TDirectory *curdir = gDirectory;
	if (!subdir.empty()) {
	TDirectory newdir = (TDirectory )gDirectory->mkdir(subdir.c_str());
	newdir->cd();
	}

	fSampleFCN->Write();

	// Change back to current DIR
	curdir->cd();

	return;
	}

	//*************************************
	void ComparisonRoutines::SaveNominal() {
	//*************************************

	fOutputRootFile->cd();

	NUIS_LOG(FIT, "Saving Nominal Predictions (be cautious with this)");
	FitBase::GetRW()->Reconfigure();
	GenerateComparison();
	SaveCurrentState("nominal");
	};
	diff --git a/src/Routines/ComparisonRoutines.h b/src/Routines/ComparisonRoutines.h
	index 1d967ac..9111ef8 100755
	--- a/src/Routines/ComparisonRoutines.h
	+++ b/src/Routines/ComparisonRoutines.h
	@@ -1,168 +1,175 @@
	// 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 COMPARISON_ROUTINES_H
	#define COMPARISON_ROUTINES_H
	/! \addtogroup Routines @{ /

	#include "TH1.h"
	#include "TF1.h"
	#include "TMatrixD.h"
	#include "TVectorD.h"
	#include "TSystem.h"
	#include "TFile.h"
	#include "TProfile.h"

	#include <vector>
	#include <string>
	#include <iostream>
	#include <sstream>
	#include <cstring>

	#include "FitEvent.h"
	#include "JointFCN.h"

	#include "GeneralUtils.h"
	#include "NuisConfig.h"
	#include "NuisKey.h"
	#include "FitLogger.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 ComparisonRoutines {
	//*************************************

	public:

	/*
	Constructor/Destructor
	*/

	/// Constructor reads in arguments given at the command line for the fit here.
	ComparisonRoutines(int argc, char* argv[]);

	/// Default destructor
	~ComparisonRoutines();

	/// Reset everything to default/NULL
	void Init();

	/*
	Input Functions
	*/

	/// Queries configuration keys to setup Parameters/Samples/FakeParameters
	void SetupComparisonsFromXML();

	/*
	Setup Functions
	*/

	/// Setups up our custom RW engine with all the parameters passed in the card file
	void SetupRWEngine();

	/// Setups up the jointFCN.
	void SetupFCN();

	/// Set the current data histograms in each sample to the fake data.
	void SetFakeData();

	/*
	Fitting Functions
	*/

	/// Main function to actually start iterating over the different required fit routines
	void Run();

	/// Creates a comparison from FCN
	void GenerateComparison();

	/// Given a new map change the values that the RW engine is currently set to
	void UpdateRWEngine(std::map<std::string,double>& updateVals);

	/// Print current value
	void PrintState();

	/*
	Write Functions
	*/

	/// Save the sample plots for current MC
	/// dir if not empty forces plots to be saved in a subdirectory of outputfile
	void SaveCurrentState(std::string subdir="");

	/// Save starting predictions into a seperate folder
	void SaveNominal();

	/*
	MISC Functions
	*/

	/// Get previous fit status from a file
	Int_t GetStatus();

	protected:

	//! Our Custom ReWeight Object
	FitWeight* rw;

	std::string fOutputFile; ///< Output file name
	// std::string fInputFile; ///< Input file name

	// TFile* fInputRootFile; ///<
	TFile* fOutputRootFile; ///< Output ROOT TFile

	JointFCN* fSampleFCN; ///< Joint Samples Container that handles reconfigures.

	std::string fCardFile; ///< Input card/XML file.

	std::string fStrategy; ///< Comparison routine selection.
	std::vector<std::string> fRoutines; ///< Split vector of comparison routine selection.
	std::string fAllowedRoutines; ///< Hard coded list of allowed routines.

	/// Fake data flag. Can be 'MC' to use 'fake_parameter'
	/// or 'path_to_file.root' to use previous NUISANCE MC predictions.
	std::string fFakeDataInput;

	// Input Dial Vals
	std::vector<std::string> fParams; ///< Vector of dial names.
	std::map<std::string, std::string> fStateVals; ///< Map of dial states
	std::map<std::string, double> fCurVals; ///< Map of dial values
	std::map<std::string, int> fTypeVals; ///< Map of dial type enums.

	+
	+ struct mirror_param {
	+ double mirror_value;
	+ bool mirror_above;
	+ };
	+ std::map<std::string, mirror_param> fMirroredParams;
	+
	// Fake Dial Vals
	std::map<std::string,double> fFakeVals; ///< Map of fake data settings.

	// Configuration
	nuiskey fCompKey; ///< Configuration Key for this Comparison Instance

	};

	/! @} /
	#endif
	diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx
	index 4e20a4a..9293e51 100755
	--- a/src/Routines/MinimizerRoutines.cxx
	+++ b/src/Routines/MinimizerRoutines.cxx
	@@ -1,1535 +1,1561 @@
	// 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 "MinimizerRoutines.h"

	#include "Simple_MH_Sampler.h"

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

	fInputFile = "";
	fInputRootFile = NULL;

	fOutputFile = "";
	fOutputRootFile = NULL;

	fCovar = NULL;
	fCovFree = NULL;
	fCorrel = NULL;
	fCorFree = NULL;
	fDecomp = NULL;
	fDecFree = NULL;

	fStrategy = "Migrad,FixAtLimBreak,Migrad";
	fRoutines.clear();

	fCardFile = "";

	fFakeDataInput = "";

	fSampleFCN = NULL;

	fMinimizer = NULL;
	fMinimizerFCN = NULL;
	fCallFunctor = NULL;

	fAllowedRoutines = ("Migrad,Simplex,Combined,"
	"Brute,Fumili,ConjugateFR,"
	"ConjugatePR,BFGS,BFGS2,"
	"SteepDesc,GSLSimAn,FixAtLim,FixAtLimBreak,"
	"Chi2Scan1D,Chi2Scan2D,Contours,ErrorBands,"
	"DataToys,MCMC");
	};

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

	/*
	Input Functions
	*/
	//*************************************
	MinimizerRoutines::MinimizerRoutines() {
	//*************************************
	Init();
	}

	//*************************************
	MinimizerRoutines::MinimizerRoutines(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, "-d", fFakeDataInput, 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()) {
	NUIS_ABORT("No input supplied!");
	}

	if (fOutputFile.empty() and !fCardFile.empty()) {
	fOutputFile = fCardFile + ".root";
	NUIS_ERR(WRN, "No output supplied so saving it to: " << fOutputFile);

	} else if (fOutputFile.empty()) {
	NUIS_ABORT("No output file or cardfile supplied!");
	}

	// 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.Set("cardfile", fCardFile);
	if (!fOutputFile.empty())
	fCompKey.Set("outputfile", fOutputFile);
	if (!fStrategy.empty())
	fCompKey.Set("strategy", fStrategy);

	// Load XML Cardfile
	configuration.LoadSettings(fCompKey.GetS("cardfile"), "");

	// 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.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml");

	// Add Error Verbo Lines
	verbocount += Config::GetParI("VERBOSITY");
	errorcount += Config::GetParI("ERROR");
	NUIS_LOG(FIT, "[ NUISANCE ]: Setting VERBOSITY=" << verbocount);
	NUIS_LOG(FIT, "[ NUISANCE ]: Setting ERROR=" << errorcount);
	// FitPar::log_verb = verbocount;
	SETVERBOSITY(verbocount);
	// ERR_VERB(errorcount);

	// Minimizer Setup ========================================
	fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
	SetupMinimizerFromXML();

	SetupCovariance();
	SetupRWEngine();
	SetupFCN();

	return;
	};

	//*************************************
	void MinimizerRoutines::SetupMinimizerFromXML() {
	//*************************************

	NUIS_LOG(FIT, "Setting up nuismin");

	// Setup Parameters ------------------------------------------
	std::vector<nuiskey> parkeys = Config::QueryKeys("parameter");
	if (!parkeys.empty()) {
	NUIS_LOG(FIT, "Number of parameters : " << parkeys.size());
	}

	for (size_t i = 0; i < parkeys.size(); i++) {
	nuiskey key = parkeys.at(i);

	// Check for type,name,nom
	if (!key.Has("type")) {
	NUIS_ABORT("No type given for parameter " << i);
	} else if (!key.Has("name")) {
	NUIS_ABORT("No name given for parameter " << i);
	} else if (!key.Has("nominal")) {
	NUIS_ABORT("No nominal given for parameter " << i);
	}

	// 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 state if none 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");

	NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom
	- << " : " << parlow << " < p < " << parhigh << " : "
	- << parstate);
	+ << " : " << parlow << " < p < " << parhigh << " : "
	+ << parstate);
	} else {
	NUIS_LOG(FIT, "Read " << partype << " : " << parname << " = " << parnom
	- << " : " << parstate);
	+ << " : " << parstate);
	+ }
	+
	+ bool ismirr = false;
	+ if (key.Has("mirror_point")) {
	+ ismirr=true;
	+ mirror_param mir;
	+ mir.mirror_value = key.GetD("mirror_point");
	+ mir.mirror_above = key.GetB("mirror_above");
	+ fMirroredParams[parname] = mir;
	+ NUIS_LOG(FIT,
	+ "\t\t" << parname << " is mirrored at " << mir.mirror_value
	+ << " "
	+ << (mir.mirror_above ? "from above" : "from below"));
	}

	// Run Parameter Conversion if needed
	if (parstate.find("ABS") != std::string::npos) {
	+ if(ismirr){
	+ NUIS_ABORT("Cannot mirror parameters with ABS state!");
	+ }
	double opnom = parnom;
	double oparstep = parstep;
	parnom = FitBase::RWAbsToSigma(partype, parname, parnom);
	parlow = FitBase::RWAbsToSigma(partype, parname, parlow);
	parhigh = FitBase::RWAbsToSigma(partype, parname, parhigh);
	parstep =
	FitBase::RWAbsToSigma(partype, parname, opnom + parstep) - parnom;
	NUIS_LOG(FIT, "ParStep: " << parstep << " (" << oparstep << ").");
	} else if (parstate.find("FRAC") != std::string::npos) {
	+ if(ismirr){
	+ NUIS_ABORT("Cannot mirror parameters with FRAC state!");
	+ }
	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()) {
	NUIS_LOG(FIT, "Number of samples : " << samplekeys.size());
	}

	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
	- NUIS_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);
	+ NUIS_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);

	// If FREE add to parameters otherwise continue
	if (sampletype.find("FREE") == std::string::npos) {
	if (samplenorm != 1.0) {
	NUIS_ERR(FTL, "You provided a sample normalisation but did not specify "
	- "that the sample is free");
	+ "that the sample is free");
	NUIS_ABORT("Change so sample contains type=\"FREE\" and re-run");
	}
	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;
	fStartVals[normname] = samplenorm;
	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()) {
	NUIS_LOG(FIT, "Number of fake parameters : " << fakekeys.size());
	}

	for (size_t i = 0; i < fakekeys.size(); i++) {
	nuiskey key = fakekeys.at(i);

	// Check for type,name,nom
	if (!key.Has("name")) {
	NUIS_ABORT("No name given for fakeparameter " << i);
	} else if (!key.Has("nom")) {
	NUIS_ABORT("No nominal given for fakeparameter " << i);
	}

	// Get Inputs
	std::string parname = key.GetS("name");
	double parnom = key.GetD("nom");

	// Push into vectors
	fFakeVals[parname] = parnom;
	}
	}

	/*
	Setup Functions
	*/
	//*************************************
	void MinimizerRoutines::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 MinimizerRoutines::SetupFCN() {
	//*************************************

	NUIS_LOG(FIT, "Making the jointFCN");
	if (fSampleFCN)
	delete fSampleFCN;
	// fSampleFCN = new JointFCN(fCardFile, fOutputRootFile);
	fSampleFCN = new JointFCN(fOutputRootFile);

	SetFakeData();

	fMinimizerFCN = new MinimizerFCN(fSampleFCN);
	fCallFunctor = new ROOT::Math::Functor(*fMinimizerFCN, fParams.size());

	fSampleFCN->CreateIterationTree("fit_iterations", FitBase::GetRW());

	return;
	}

	//******************************************
	void MinimizerRoutines::SetupFitter(std::string routine) {
	//******************************************

	// Make the fitter
	std::string fitclass = "";
	std::string fittype = "";

	bool UseMCMC = false;
	// Get correct types
	if (!routine.compare("Migrad")) {
	fitclass = "Minuit2";
	fittype = "Migrad";
	} else if (!routine.compare("Simplex")) {
	fitclass = "Minuit2";
	fittype = "Simplex";
	} else if (!routine.compare("Combined")) {
	fitclass = "Minuit2";
	fittype = "Combined";
	} else if (!routine.compare("Brute")) {
	fitclass = "Minuit2";
	fittype = "Scan";
	} else if (!routine.compare("Fumili")) {
	fitclass = "Minuit2";
	fittype = "Fumili";
	} else if (!routine.compare("ConjugateFR")) {
	fitclass = "GSLMultiMin";
	fittype = "ConjugateFR";
	} else if (!routine.compare("ConjugatePR")) {
	fitclass = "GSLMultiMin";
	fittype = "ConjugatePR";
	} else if (!routine.compare("BFGS")) {
	fitclass = "GSLMultiMin";
	fittype = "BFGS";
	} else if (!routine.compare("BFGS2")) {
	fitclass = "GSLMultiMin";
	fittype = "BFGS2";
	} else if (!routine.compare("SteepDesc")) {
	fitclass = "GSLMultiMin";
	fittype = "SteepestDescent";
	// } else if (!routine.compare("GSLMulti")) { fitclass = "GSLMultiFit";
	// fittype = ""; // Doesn't work out of the box
	} else if (!routine.compare("GSLSimAn")) {
	fitclass = "GSLSimAn";
	fittype = "";
	} else if (!routine.compare("MCMC")) {
	UseMCMC = true;
	}

	// make minimizer
	if (fMinimizer)
	delete fMinimizer;

	if (UseMCMC) {
	fMinimizer = new Simple_MH_Sampler();
	} else {
	fMinimizer = ROOT::Math::Factory::CreateMinimizer(fitclass, fittype);
	}

	fMinimizer->SetMaxFunctionCalls(FitPar::Config().GetParI("MAXCALLS"));

	if (!routine.compare("Brute")) {
	fMinimizer->SetMaxFunctionCalls(fParams.size() * fParams.size() * 4);
	fMinimizer->SetMaxIterations(fParams.size() * fParams.size() * 4);
	}

	fMinimizer->SetMaxIterations(FitPar::Config().GetParI("MAXITERATIONS"));
	fMinimizer->SetTolerance(FitPar::Config().GetParD("TOLERANCE"));
	fMinimizer->SetStrategy(FitPar::Config().GetParI("STRATEGY"));
	fMinimizer->SetFunction(*fCallFunctor);

	int ipar = 0;
	// Add Fit Parameters
	for (UInt_t i = 0; i < fParams.size(); i++) {
	std::string syst = fParams.at(i);

	bool fixed = true;
	double vstart, vstep, vlow, vhigh;
	vstart = vstep = vlow = vhigh = 0.0;

	if (fCurVals.find(syst) != fCurVals.end())
	vstart = fCurVals.at(syst);
	if (fMinVals.find(syst) != fMinVals.end())
	vlow = fMinVals.at(syst);
	if (fMaxVals.find(syst) != fMaxVals.end())
	vhigh = fMaxVals.at(syst);
	if (fStepVals.find(syst) != fStepVals.end())
	vstep = fStepVals.at(syst);
	if (fFixVals.find(syst) != fFixVals.end())
	fixed = fFixVals.at(syst);

	// fix for errors
	if (vhigh == vlow)
	vhigh += 1.0;

	fMinimizer->SetVariable(ipar, syst, vstart, vstep);
	fMinimizer->SetVariableLimits(ipar, vlow, vhigh);
	+ if (fMirroredParams.count(syst)) {
	+ fSampleFCN->SetVariableMirrored(ipar, fMirroredParams[syst].mirror_value,
	+ fMirroredParams[syst].mirror_above);
	+ }

	if (fixed) {
	fMinimizer->FixVariable(ipar);
	NUIS_LOG(FIT, "Fixed Param: " << syst);

	} else {
	NUIS_LOG(FIT, "Free Param: " << syst << " Start:" << vstart
	- << " Range:" << vlow << " to " << vhigh
	- << " Step:" << vstep);
	+ << " Range:" << vlow << " to " << vhigh
	+ << " Step:" << vstep);
	}

	ipar++;
	}
	+ fSampleFCN->SetNParams(ipar);

	NUIS_LOG(FIT, "Setup Minimizer: " << fMinimizer->NDim() << "(NDim) "
	- << fMinimizer->NFree() << "(NFree)");
	+ << fMinimizer->NFree() << "(NFree)");

	return;
	}

	//*************************************
	// Set fake data from user input
	void MinimizerRoutines::SetFakeData() {
	//*************************************

	// If the fake data input field (-d) isn't provided, return to caller
	if (fFakeDataInput.empty())
	return;

	// If user specifies -d MC we set the data to the MC
	// User can also specify fake data parameters to reweight by doing
	// "fake_parameter" in input card file
	// "fake_parameter" gets read in ReadCard function (reads to fFakeVals)
	if (fFakeDataInput.compare("MC") == 0) {
	NUIS_LOG(FIT, "Setting fake data from MC starting prediction.");
	// fFakeVals get read in in ReadCard
	UpdateRWEngine(fFakeVals);

	// Reconfigure the reweight engine
	FitBase::GetRW()->Reconfigure();
	// Reconfigure all the samples to the new reweight
	fSampleFCN->ReconfigureAllEvents();
	// Feed on and set the fake-data in each measurement class
	fSampleFCN->SetFakeData("MC");

	// Changed the reweight engine values back to the current values
	// So we start the fit at a different value than what we set the fake-data
	// to
	UpdateRWEngine(fCurVals);

	NUIS_LOG(FIT, "Set all data to fake MC predictions.");
	} else {
	fSampleFCN->SetFakeData(fFakeDataInput);
	}

	return;
	}

	/*
	Fitting Functions
	*/
	//*************************************
	void MinimizerRoutines::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 MinimizerRoutines::Run() {
	//*************************************

	NUIS_LOG(FIT, "Running MinimizerRoutines : " << fStrategy);
	if (FitPar::Config().GetParB("save_nominal")) {
	SaveNominal();
	}

	// Parse given routines
	fRoutines = GeneralUtils::ParseToStr(fStrategy, ",");
	if (fRoutines.empty()) {
	NUIS_ABORT("Trying to run MinimizerRoutines with no routines given!");
	}

	for (UInt_t i = 0; i < fRoutines.size(); i++) {
	std::string routine = fRoutines.at(i);
	int fitstate = kFitUnfinished;
	NUIS_LOG(FIT, "Running Routine: " << routine);

	// Try Routines
	if (routine.find("LowStat") != std::string::npos)
	LowStatRoutine(routine);
	else if (routine == "FixAtLim")
	FixAtLimit();
	else if (routine == "FixAtLimBreak")
	fitstate = FixAtLimit();
	else if (routine.find("ErrorBands") != std::string::npos)
	GenerateErrorBands();
	else if (routine.find("DataToys") != std::string::npos)
	ThrowDataToys();
	else if (!routine.compare("Chi2Scan1D"))
	Create1DScans();
	else if (!routine.compare("Chi2Scan2D"))
	Chi2Scan2D();
	else
	fitstate = RunFitRoutine(routine);

	// If ending early break here
	if (fitstate == kFitFinished \|\| fitstate == kNoChange) {
	NUIS_LOG(FIT, "Ending fit routines loop.");
	break;
	}
	}

	return;
	}

	//*************************************
	int MinimizerRoutines::RunFitRoutine(std::string routine) {
	//*************************************
	int endfits = kFitUnfinished;

	// set fitter at the current start values
	fOutputRootFile->cd();
	SetupFitter(routine);

	// choose what to do with the minimizer depending on routine.
	if (!routine.compare("Migrad") or !routine.compare("Simplex") or
	!routine.compare("Combined") or !routine.compare("Brute") or
	!routine.compare("Fumili") or !routine.compare("ConjugateFR") or
	!routine.compare("ConjugatePR") or !routine.compare("BFGS") or
	!routine.compare("BFGS2") or !routine.compare("SteepDesc") or
	// !routine.compare("GSLMulti") or
	!routine.compare("GSLSimAn") or !routine.compare("MCMC")) {
	if (fMinimizer->NFree() > 0) {
	NUIS_LOG(FIT, fMinimizer->Minimize());
	GetMinimizerState();
	}
	}

	// other otptions
	else if (!routine.compare("Contour")) {
	CreateContours();
	}

	return endfits;
	}

	//*************************************
	void MinimizerRoutines::PrintState() {
	//*************************************
	NUIS_LOG(FIT, "------------");

	// Count max size
	int maxcount = 0;
	for (UInt_t i = 0; i < fParams.size(); i++) {
	maxcount = max(int(fParams[i].size()), maxcount);
	}

	// Header
	NUIS_LOG(FIT, " # " << left << setw(maxcount) << "Parameter "
	- << " = " << setw(10) << "Value"
	- << " +- " << setw(10) << "Error"
	- << " " << setw(8) << "(Units)"
	- << " " << setw(10) << "Conv. Val"
	- << " +- " << setw(10) << "Conv. Err"
	- << " " << setw(8) << "(Units)");
	+ << " = " << setw(10) << "Value"
	+ << " +- " << setw(10) << "Error"
	+ << " " << setw(8) << "(Units)"
	+ << " " << setw(10) << "Conv. Val"
	+ << " +- " << setw(10) << "Conv. Err"
	+ << " " << setw(8) << "(Units)");

	// Parameters
	for (UInt_t i = 0; i < fParams.size(); i++) {
	std::string syst = fParams.at(i);

	std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
	std::string curunits = "(sig.)";
	double curval = fCurVals[syst];
	double curerr = fErrorVals[syst];

	if (fStateVals[syst].find("ABS") != std::string::npos) {
	curval = FitBase::RWSigmaToAbs(typestr, syst, curval);
	curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
	FitBase::RWSigmaToAbs(typestr, syst, 0.0));
	curunits = "(Abs.)";
	} else if (fStateVals[syst].find("FRAC") != std::string::npos) {
	curval = FitBase::RWSigmaToFrac(typestr, syst, curval);
	curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) -
	FitBase::RWSigmaToFrac(typestr, syst, 0.0));
	curunits = "(Frac)";
	}

	std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")";
	double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
	double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
	FitBase::RWSigmaToAbs(typestr, syst, 0.0));

	std::ostringstream curparstring;

	curparstring << " " << setw(3) << left << i << ". " << setw(maxcount)
	<< syst << " = " << setw(10) << curval << " +- " << setw(10)
	<< curerr << " " << setw(8) << curunits << " " << setw(10)
	<< convval << " +- " << setw(10) << converr << " " << setw(8)
	<< convunits;

	NUIS_LOG(FIT, curparstring.str());
	}

	NUIS_LOG(FIT, "------------");
	double like = fSampleFCN->GetLikelihood();
	- NUIS_LOG(FIT, std::left << std::setw(46) << "Likelihood for JointFCN: " << like);
	+ NUIS_LOG(FIT,
	+ std::left << std::setw(46) << "Likelihood for JointFCN: " << like);
	NUIS_LOG(FIT, "------------");
	}

	//*************************************
	void MinimizerRoutines::GetMinimizerState() {
	//*************************************

	NUIS_LOG(FIT, "Minimizer State: ");
	// Get X and Err
	const double *values = fMinimizer->X();
	const double *errors = fMinimizer->Errors();
	// int ipar = 0;

	for (UInt_t i = 0; i < fParams.size(); i++) {
	std::string syst = fParams.at(i);

	fCurVals[syst] = values[i];
	fErrorVals[syst] = errors[i];
	}

	PrintState();

	// Covar
	SetupCovariance();
	if (fMinimizer->CovMatrixStatus() > 0) {
	// Fill Full Covar
	NUIS_LOG(FIT, "Filling covariance");
	for (int i = 0; i < fCovar->GetNbinsX(); i++) {
	for (int j = 0; j < fCovar->GetNbinsY(); j++) {
	fCovar->SetBinContent(i + 1, j + 1, fMinimizer->CovMatrix(i, j));
	}
	}

	int freex = 0;
	int freey = 0;
	for (int i = 0; i < fCovar->GetNbinsX(); i++) {
	if (fMinimizer->IsFixedVariable(i))
	continue;
	freey = 0;

	for (int j = 0; j < fCovar->GetNbinsY(); j++) {
	if (fMinimizer->IsFixedVariable(j))
	continue;

	fCovFree->SetBinContent(freex + 1, freey + 1,
	fMinimizer->CovMatrix(i, j));
	freey++;
	}
	freex++;
	}

	fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation");
	fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition");
	if (fMinimizer->NFree() > 0) {
	fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free");
	fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free");
	}
	}

	return;
	};

	//*************************************
	void MinimizerRoutines::LowStatRoutine(std::string routine) {
	//*************************************

	NUIS_LOG(FIT, "Running Low Statistics Routine: " << routine);
	int lowstatsevents = FitPar::Config().GetParI("LOWSTATEVENTS");
	int maxevents = FitPar::Config().GetParI("MAXEVENTS");
	int verbosity = FitPar::Config().GetParI("VERBOSITY");

	std::string trueroutine = routine;
	std::string substring = "LowStat";
	trueroutine.erase(trueroutine.find(substring), substring.length());

	// Set MAX EVENTS=1000
	Config::SetPar("MAXEVENTS", lowstatsevents);
	Config::SetPar("VERBOSITY", 3);
	SetupFCN();

	RunFitRoutine(trueroutine);

	Config::SetPar("MAXEVENTS", maxevents);
	SetupFCN();

	Config::SetPar("VERBOSITY", verbosity);
	return;
	}

	//*************************************
	void MinimizerRoutines::Create1DScans() {
	//*************************************

	// 1D Scan Routine
	// Steps through all free parameters about nominal using the step size
	// Creates a graph for each free parameter

	// At the current point create a 1D Scan for all parametes (Uncorrelated)
	for (UInt_t i = 0; i < fParams.size(); i++) {
	if (fFixVals[fParams[i]])
	continue;

	NUIS_LOG(FIT, "Running 1D Scan for " << fParams[i]);
	fSampleFCN->CreateIterationTree(fParams[i] + "_scan1D_iterations",
	FitBase::GetRW());

	double scanmiddlepoint = fCurVals[fParams[i]];

	// Determine N points needed
	double limlow = fMinVals[fParams[i]];
	double limhigh = fMaxVals[fParams[i]];
	double step = fStepVals[fParams[i]];

	int npoints = int(fabs(limhigh - limlow) / (step + 0.));

	TH1D *contour =
	new TH1D(("Chi2Scan1D_" + fParams[i]).c_str(),
	("Chi2Scan1D_" + fParams[i] + ";" + fParams[i]).c_str(),
	npoints, limlow, limhigh);

	// Fill bins
	for (int x = 0; x < contour->GetNbinsX(); x++) {
	// Set X Val
	fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1);

	// Run Eval
	double *vals = FitUtils::GetArrayFromMap(fParams, fCurVals);
	double chi2 = fSampleFCN->DoEval(vals);
	delete vals;

	// Fill Contour
	contour->SetBinContent(x + 1, chi2);
	}

	// Save contour
	contour->Write();

	// Reset Parameter
	fCurVals[fParams[i]] = scanmiddlepoint;

	// Save TTree
	fSampleFCN->WriteIterationTree();
	}

	return;
	}

	//*************************************
	void MinimizerRoutines::Chi2Scan2D() {
	//*************************************

	// Chi2 Scan 2D
	// Creates a 2D chi2 scan by stepping through all free parameters
	// Works for all pairwise combos of free parameters

	// Scan I
	for (UInt_t i = 0; i < fParams.size(); i++) {
	if (fFixVals[fParams[i]])
	continue;

	// Scan J
	for (UInt_t j = 0; j < i; j++) {
	if (fFixVals[fParams[j]])
	continue;

	fSampleFCN->CreateIterationTree(fParams[i] + "_" + fParams[j] + "_" +
	"scan2D_iterations",
	FitBase::GetRW());

	double scanmid_i = fCurVals[fParams[i]];
	double scanmid_j = fCurVals[fParams[j]];

	double limlow_i = fMinVals[fParams[i]];
	double limhigh_i = fMaxVals[fParams[i]];
	double step_i = fStepVals[fParams[i]];

	double limlow_j = fMinVals[fParams[j]];
	double limhigh_j = fMaxVals[fParams[j]];
	double step_j = fStepVals[fParams[j]];

	int npoints_i = int(fabs(limhigh_i - limlow_i) / (step_i + 0.)) + 1;
	int npoints_j = int(fabs(limhigh_j - limlow_j) / (step_j + 0.)) + 1;

	TH2D *contour = new TH2D(
	("Chi2Scan2D_" + fParams[i] + "_" + fParams[j]).c_str(),
	("Chi2Scan2D_" + fParams[i] + "_" + fParams[j] + ";" + fParams[i] +
	";" + fParams[j])
	.c_str(),
	npoints_i, limlow_i, limhigh_i, npoints_j, limlow_j, limhigh_j);

	// Begin Scan
	NUIS_LOG(FIT, "Running scan for " << fParams[i] << " " << fParams[j]);

	// Fill bins
	for (int x = 0; x < contour->GetNbinsX(); x++) {
	// Set X Val
	fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1);

	// Loop Y
	for (int y = 0; y < contour->GetNbinsY(); y++) {
	// Set Y Val
	fCurVals[fParams[j]] = contour->GetYaxis()->GetBinCenter(y + 1);

	// Run Eval
	double *vals = FitUtils::GetArrayFromMap(fParams, fCurVals);
	double chi2 = fSampleFCN->DoEval(vals);
	delete vals;

	// Fill Contour
	contour->SetBinContent(x + 1, y + 1, chi2);

	fCurVals[fParams[j]] = scanmid_j;
	}

	fCurVals[fParams[i]] = scanmid_i;
	fCurVals[fParams[j]] = scanmid_j;
	}

	// Save contour
	contour->Write();

	// Save Iterations
	fSampleFCN->WriteIterationTree();
	}
	}

	return;
	}

	//*************************************
	void MinimizerRoutines::CreateContours() {
	//*************************************

	// Use MINUIT for this if possible
	NUIS_ABORT("Contours not yet implemented as it is really slow!");

	return;
	}

	//*************************************
	int MinimizerRoutines::FixAtLimit() {
	//*************************************

	bool fixedparam = false;
	for (UInt_t i = 0; i < fParams.size(); i++) {
	std::string syst = fParams.at(i);
	if (fFixVals[syst])
	continue;

	double curVal = fCurVals.at(syst);
	double minVal = fMinVals.at(syst);
	double maxVal = fMinVals.at(syst);

	if (fabs(curVal - minVal) < 0.0001) {
	fCurVals[syst] = minVal;
	fFixVals[syst] = true;
	fixedparam = true;
	}

	if (fabs(maxVal - curVal) < 0.0001) {
	fCurVals[syst] = maxVal;
	fFixVals[syst] = true;
	fixedparam = true;
	}
	}

	if (!fixedparam) {
	NUIS_LOG(FIT, "No dials needed fixing!");
	return kNoChange;
	} else
	return kStateChange;
	}

	/*
	Write Functions
	*/
	//*************************************
	void MinimizerRoutines::SaveResults() {
	//*************************************

	fOutputRootFile->cd();

	if (fMinimizer) {
	SetupCovariance();
	SaveMinimizerState();
	}

	SaveCurrentState();
	}

	//*************************************
	void MinimizerRoutines::SaveMinimizerState() {
	//*************************************

	NUIS_LOG(FIT, "Saving Minimizer State");
	if (!fMinimizer) {
	NUIS_ABORT("Can't save minimizer state without min object");
	}

	// Save main fit tree
	fSampleFCN->WriteIterationTree();

	// Get Vals and Errors
	GetMinimizerState();

	// Save tree with fit status
	std::vector<std::string> nameVect;
	std::vector<double> valVect;
	std::vector<double> errVect;
	std::vector<double> minVect;
	std::vector<double> maxVect;
	std::vector<double> startVect;
	std::vector<int> endfixVect;
	std::vector<int> startfixVect;

	// int NFREEPARS = fMinimizer->NFree();
	int NPARS = fMinimizer->NDim();

	int ipar = 0;
	// Dial Vals
	for (UInt_t i = 0; i < fParams.size(); i++) {
	std::string name = fParams.at(i);
	nameVect.push_back(name);

	valVect.push_back(fCurVals.at(name));

	errVect.push_back(fErrorVals.at(name));

	minVect.push_back(fMinVals.at(name));

	maxVect.push_back(fMaxVals.at(name));

	startVect.push_back(fStartVals.at(name));

	endfixVect.push_back(fFixVals.at(name));

	startfixVect.push_back(fStartFixVals.at(name));

	ipar++;
	}

	int NFREE = fMinimizer->NFree();
	int NDIM = fMinimizer->NDim();

	double CHI2 = fSampleFCN->GetLikelihood();
	int NBINS = fSampleFCN->GetNDOF();
	int NDOF = NBINS - NFREE;

	// Write fit results
	TTree *fit_tree = new TTree("fit_result", "fit_result");
	fit_tree->Branch("parameter_names", &nameVect);
	fit_tree->Branch("parameter_values", &valVect);
	fit_tree->Branch("parameter_errors", &errVect);
	fit_tree->Branch("parameter_min", &minVect);
	fit_tree->Branch("parameter_max", &maxVect);
	fit_tree->Branch("parameter_start", &startVect);
	fit_tree->Branch("parameter_fix", &endfixVect);
	fit_tree->Branch("parameter_startfix", &startfixVect);
	fit_tree->Branch("CHI2", &CHI2, "CHI2/D");
	fit_tree->Branch("NDOF", &NDOF, "NDOF/I");
	fit_tree->Branch("NBINS", &NBINS, "NBINS/I");
	fit_tree->Branch("NDIM", &NDIM, "NDIM/I");
	fit_tree->Branch("NFREE", &NFREE, "NFREE/I");
	fit_tree->Fill();
	fit_tree->Write();

	// Make dial variables
	TH1D dialvar = TH1D("fit_dials", "fit_dials", NPARS, 0, NPARS);
	TH1D startvar = TH1D("start_dials", "start_dials", NPARS, 0, NPARS);
	TH1D minvar = TH1D("min_dials", "min_dials", NPARS, 0, NPARS);
	TH1D maxvar = TH1D("max_dials", "max_dials", NPARS, 0, NPARS);

	TH1D dialvarfree = TH1D("fit_dials_free", "fit_dials_free", NFREE, 0, NFREE);
	TH1D startvarfree =
	TH1D("start_dials_free", "start_dials_free", NFREE, 0, NFREE);
	TH1D minvarfree = TH1D("min_dials_free", "min_dials_free", NFREE, 0, NFREE);
	TH1D maxvarfree = TH1D("max_dials_free", "max_dials_free", NFREE, 0, NFREE);

	int freecount = 0;

	for (UInt_t i = 0; i < nameVect.size(); i++) {
	std::string name = nameVect.at(i);

	dialvar.SetBinContent(i + 1, valVect.at(i));
	dialvar.SetBinError(i + 1, errVect.at(i));
	dialvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());

	startvar.SetBinContent(i + 1, startVect.at(i));
	startvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());

	minvar.SetBinContent(i + 1, minVect.at(i));
	minvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());

	maxvar.SetBinContent(i + 1, maxVect.at(i));
	maxvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());

	if (NFREE > 0) {
	if (!startfixVect.at(i)) {
	freecount++;

	dialvarfree.SetBinContent(freecount, valVect.at(i));
	dialvarfree.SetBinError(freecount, errVect.at(i));
	dialvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());

	startvarfree.SetBinContent(freecount, startVect.at(i));
	startvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());

	minvarfree.SetBinContent(freecount, minVect.at(i));
	minvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());

	maxvarfree.SetBinContent(freecount, maxVect.at(i));
	maxvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());
	}
	}
	}

	// Save Dial Plots
	dialvar.Write();
	startvar.Write();
	minvar.Write();
	maxvar.Write();

	if (NFREE > 0) {
	dialvarfree.Write();
	startvarfree.Write();
	minvarfree.Write();
	maxvarfree.Write();
	}

	// Save fit_status plot
	TH1D statusplot = TH1D("fit_status", "fit_status", 8, 0, 8);
	std::string fit_labels[8] = {"status", "cov_status", "maxiter",
	"maxfunc", "iter", "func",
	"precision", "tolerance"};
	double fit_vals[8];
	fit_vals[0] = fMinimizer->Status() + 0.;
	fit_vals[1] = fMinimizer->CovMatrixStatus() + 0.;
	fit_vals[2] = fMinimizer->MaxIterations() + 0.;
	fit_vals[3] = fMinimizer->MaxFunctionCalls() + 0.;
	fit_vals[4] = fMinimizer->NIterations() + 0.;
	fit_vals[5] = fMinimizer->NCalls() + 0.;
	fit_vals[6] = fMinimizer->Precision() + 0.;
	fit_vals[7] = fMinimizer->Tolerance() + 0.;

	for (int i = 0; i < 8; i++) {
	statusplot.SetBinContent(i + 1, fit_vals[i]);
	statusplot.GetXaxis()->SetBinLabel(i + 1, fit_labels[i].c_str());
	}

	statusplot.Write();

	// Save Covars
	if (fCovar)
	fCovar->Write();
	if (fCovFree)
	fCovFree->Write();
	if (fCorrel)
	fCorrel->Write();
	if (fCorFree)
	fCorFree->Write();
	if (fDecomp)
	fDecomp->Write();
	if (fDecFree)
	fDecFree->Write();

	return;
	}

	//*************************************
	void MinimizerRoutines::SaveCurrentState(std::string subdir) {
	//*************************************

	NUIS_LOG(FIT, "Saving current full FCN predictions");

	// Setup DIRS
	TDirectory *curdir = gDirectory;
	if (!subdir.empty()) {
	TDirectory newdir = (TDirectory )gDirectory->mkdir(subdir.c_str());
	newdir->cd();
	}

	FitBase::GetRW()->Reconfigure();
	fSampleFCN->ReconfigureAllEvents();
	fSampleFCN->Write();

	// Change back to current DIR
	curdir->cd();

	return;
	}

	//*************************************
	void MinimizerRoutines::SaveNominal() {
	//*************************************

	fOutputRootFile->cd();

	NUIS_LOG(FIT, "Saving Nominal Predictions (be cautious with this)");
	FitBase::GetRW()->Reconfigure();
	SaveCurrentState("nominal");
	};

	//*************************************
	void MinimizerRoutines::SavePrefit() {
	//*************************************

	fOutputRootFile->cd();

	NUIS_LOG(FIT, "Saving Prefit Predictions");
	UpdateRWEngine(fStartVals);
	SaveCurrentState("prefit");
	UpdateRWEngine(fCurVals);
	};

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

	return 0;
	}

	//*************************************
	void MinimizerRoutines::SetupCovariance() {
	//*************************************

	// Remove covares if they exist
	if (fCovar)
	delete fCovar;
	if (fCovFree)
	delete fCovFree;
	if (fCorrel)
	delete fCorrel;
	if (fCorFree)
	delete fCorFree;
	if (fDecomp)
	delete fDecomp;
	if (fDecFree)
	delete fDecFree;

	NUIS_LOG(FIT, "Building covariance matrix..");

	int NFREE = 0;
	int NDIM = 0;

	// Get NFREE from min or from vals (for cases when doing throws)
	if (fMinimizer) {
	NFREE = fMinimizer->NFree();
	NDIM = fMinimizer->NDim();
	} else {
	NDIM = fParams.size();
	for (UInt_t i = 0; i < fParams.size(); i++) {
	NUIS_LOG(FIT, "Getting Param " << fParams[i]);

	if (!fFixVals[fParams[i]])
	NFREE++;
	}
	}

	if (NDIM == 0)
	return;
	NUIS_LOG(FIT, "NFREE == " << NFREE);
	fCovar = new TH2D("covariance", "covariance", NDIM, 0, NDIM, NDIM, 0, NDIM);
	if (NFREE > 0) {
	fCovFree = new TH2D("covariance_free", "covariance_free", NFREE, 0, NFREE,
	NFREE, 0, NFREE);
	} else {
	fCovFree = NULL;
	}

	// Set Bin Labels
	int countall = 0;
	int countfree = 0;
	for (UInt_t i = 0; i < fParams.size(); i++) {
	fCovar->GetXaxis()->SetBinLabel(countall + 1, fParams[i].c_str());
	fCovar->GetYaxis()->SetBinLabel(countall + 1, fParams[i].c_str());
	countall++;

	if (!fFixVals[fParams[i]] and NFREE > 0) {
	fCovFree->GetXaxis()->SetBinLabel(countfree + 1, fParams[i].c_str());
	fCovFree->GetYaxis()->SetBinLabel(countfree + 1, fParams[i].c_str());
	countfree++;
	}
	}

	fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation");
	fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition");

	if (NFREE > 0) {
	fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free");
	fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free");
	} else {
	fCorFree = NULL;
	fDecFree = NULL;
	}

	return;
	};

	//*************************************
	void MinimizerRoutines::ThrowCovariance(bool uniformly) {
	//*************************************
	std::vector<double> rands;

	if (!fDecFree) {
	NUIS_ERR(WRN, "Trying to throw 0 free parameters");
	return;
	}

	// Generate Random Gaussians
	for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) {
	rands.push_back(gRandom->Gaus(0.0, 1.0));
	}

	// Reset Thrown Values
	for (UInt_t i = 0; i < fParams.size(); i++) {
	fThrownVals[fParams[i]] = fCurVals[fParams[i]];
	}

	// Loop and get decomp
	for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) {
	std::string parname = std::string(fDecFree->GetXaxis()->GetBinLabel(i + 1));
	double mod = 0.0;

	if (!uniformly) {
	for (Int_t j = 0; j < fDecFree->GetNbinsY(); j++) {
	mod += rands[j] * fDecFree->GetBinContent(j + 1, i + 1);
	}
	}

	if (fCurVals.find(parname) != fCurVals.end()) {
	if (uniformly)
	fThrownVals[parname] =
	gRandom->Uniform(fMinVals[parname], fMaxVals[parname]);
	else {
	fThrownVals[parname] = fCurVals[parname] + mod;
	}
	}
	}

	// Check Limits
	for (UInt_t i = 0; i < fParams.size(); i++) {
	std::string syst = fParams[i];
	if (fFixVals[syst])
	continue;
	if (fThrownVals[syst] < fMinVals[syst])
	fThrownVals[syst] = fMinVals[syst];
	if (fThrownVals[syst] > fMaxVals[syst])
	fThrownVals[syst] = fMaxVals[syst];
	}

	return;
	};

	//*************************************
	void MinimizerRoutines::GenerateErrorBands() {
	//*************************************

	TDirectory errorDIR = (TDirectory )fOutputRootFile->mkdir("error_bands");
	errorDIR->cd();

	// Make a second file to store throws
	std::string tempFileName = fOutputFile;
	if (tempFileName.find(".root") != std::string::npos)
	tempFileName.erase(tempFileName.find(".root"), 5);
	tempFileName += ".throws.root";
	TFile *tempfile = new TFile(tempFileName.c_str(), "RECREATE");

	tempfile->cd();
	int nthrows = FitPar::Config().GetParI("error_throws");

	UpdateRWEngine(fCurVals);
	fSampleFCN->ReconfigureAllEvents();

	TDirectory nominal = (TDirectory )tempfile->mkdir("nominal");
	nominal->cd();
	fSampleFCN->Write();

	TDirectory *outnominal =
	(TDirectory *)fOutputRootFile->mkdir("nominal_throw");
	outnominal->cd();
	fSampleFCN->Write();

	errorDIR->cd();
	TTree *parameterTree = new TTree("throws", "throws");
	double chi2;
	for (UInt_t i = 0; i < fParams.size(); i++)
	parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]],
	(fParams[i] + "/D").c_str());
	parameterTree->Branch("chi2", &chi2, "chi2/D");

	bool uniformly = FitPar::Config().GetParB("error_uniform");

	// Run Throws and save
	for (Int_t i = 0; i < nthrows; i++) {
	TDirectory *throwfolder =
	(TDirectory *)tempfile->mkdir(Form("throw_%i", i));
	throwfolder->cd();

	// Generate Random Parameter Throw
	ThrowCovariance(uniformly);

	// Run Eval
	double *vals = FitUtils::GetArrayFromMap(fParams, fThrownVals);
	chi2 = fSampleFCN->DoEval(vals);
	delete vals;

	// Save the FCN
	fSampleFCN->Write();

	parameterTree->Fill();
	}

	errorDIR->cd();
	fDecFree->Write();
	fCovFree->Write();
	parameterTree->Write();

	delete parameterTree;

	// Now go through the keys in the temporary file and look for TH1D, and TH2D
	// plots
	TIter next(nominal->GetListOfKeys());
	TKey *key;
	while ((key = (TKey *)next())) {
	TClass *cl = gROOT->GetClass(key->GetClassName());
	if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D"))
	continue;
	TH1D baseplot = (TH1D )key->ReadObj();
	std::string plotname = std::string(baseplot->GetName());

	int nbins = baseplot->GetNbinsX() * baseplot->GetNbinsY();

	// Setup TProfile with RMS option
	TProfile *tprof =
	new TProfile((plotname + "_prof").c_str(), (plotname + "_prof").c_str(),
	nbins, 0, nbins, "S");

	// Setup The TTREE
	double *bincontents;
	bincontents = new double[nbins];

	double *binlowest;
	binlowest = new double[nbins];

	double *binhighest;
	binhighest = new double[nbins];

	errorDIR->cd();
	TTree *bintree =
	new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str());
	for (Int_t i = 0; i < nbins; i++) {
	bincontents[i] = 0.0;
	binhighest[i] = 0.0;
	binlowest[i] = 0.0;
	bintree->Branch(Form("content_%i", i), &bincontents[i],
	Form("content_%i/D", i));
	}

	for (Int_t i = 0; i < nthrows; i++) {
	TH1 *newplot =
	(TH1 *)tempfile->Get(Form(("throw_%i/" + plotname).c_str(), i));

	for (Int_t j = 0; j < nbins; j++) {
	tprof->Fill(j + 0.5, newplot->GetBinContent(j + 1));
	bincontents[j] = newplot->GetBinContent(j + 1);

	if (bincontents[j] < binlowest[j] or i == 0)
	binlowest[j] = bincontents[j];
	if (bincontents[j] > binhighest[j] or i == 0)
	binhighest[j] = bincontents[j];
	}

	errorDIR->cd();
	bintree->Fill();

	delete newplot;
	}

	errorDIR->cd();

	for (Int_t j = 0; j < nbins; j++) {
	if (!uniformly) {
	baseplot->SetBinError(j + 1, tprof->GetBinError(j + 1));

	} else {
	baseplot->SetBinContent(j + 1, (binlowest[j] + binhighest[j]) / 2.0);
	baseplot->SetBinError(j + 1, (binhighest[j] - binlowest[j]) / 2.0);
	}
	}

	errorDIR->cd();
	baseplot->Write();
	tprof->Write();
	bintree->Write();

	delete baseplot;
	delete tprof;
	delete bintree;
	delete[] bincontents;
	}

	return;
	};

	void MinimizerRoutines::ThrowDataToys() {
	NUIS_LOG(FIT, "Generating Toy Data Throws");
	int verb = Config::GetParI("VERBOSITY");
	SETVERBOSITY(FIT);

	int nthrows = FitPar::Config().GetParI("NToyThrows");
	double maxlike = -1.0;
	double minlike = -1.0;
	std::vector<double> values;
	for (int i = 0; i < 1.E4; i++) {
	fSampleFCN->ThrowDataToy();
	double like = fSampleFCN->GetLikelihood();
	values.push_back(like);
	if (maxlike == -1.0 or like > maxlike)
	maxlike = like;
	if (minlike == -1.0 or like < minlike)
	minlike = like;
	}
	SETVERBOSITY(verb);

	// Fill Histogram
	TH1D *likes = new TH1D("toydatalikelihood", "toydatalikelihood",
	int(sqrt(nthrows)), minlike, maxlike);
	for (size_t i = 0; i < values.size(); i++) {
	likes->Fill(values[i]);
	}

	// Save to file
	NUIS_LOG(FIT, "Writing toy data throws");
	fOutputRootFile->cd();
	likes->Write();
	}
	diff --git a/src/Routines/MinimizerRoutines.h b/src/Routines/MinimizerRoutines.h
	index 3918e7a..dcff46b 100755
	--- a/src/Routines/MinimizerRoutines.h
	+++ b/src/Routines/MinimizerRoutines.h
	@@ -1,278 +1,284 @@
	// 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 MINIMIZER_ROUTINES_H
	#define MINIMIZER_ROUTINES_H

	/*!
	* \addtogroup Minimizer
	* @{
	*/

	#include "TH1.h"
	#include "TF1.h"
	#include "TMatrixD.h"
	#include "TVectorD.h"

	#ifdef ROOT6_USE_FIT_FITTER_INTERFACE
	#include "Fit/Fitter.h"
	#else
	#include "Minuit2/FCNBase.h"
	#include "TFitterMinuit.h"
	#endif

	#include "TSystem.h"
	#include "TFile.h"
	#include "TProfile.h"


	#include <vector>
	#include <string>
	#include <iostream>
	#include <sstream>
	#include <cstring>

	#include "FitEvent.h"
	#include "JointFCN.h"
	#include "MinimizerFCN.h"


	#include "Math/Minimizer.h"
	#include "Math/Factory.h"
	#include "Math/Functor.h"
	#include "FitLogger.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 MinimizerRoutines{
	//*************************************

	public:

	/*
	Constructor/Destructor
	*/

	MinimizerRoutines();

	//! Constructor reads in arguments given at the command line for the fit here.
	MinimizerRoutines(int argc, char* argv[]);

	//! Default destructor
	~MinimizerRoutines();

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

	//! Loops through each line of the card file and passes it to other read functions
	void ReadCard(std::string cardfile);

	//! Check for parameter string in the line and assign the correct type.
	//! Fills maps for each of the parameters
	int ReadParameters(std::string parstring);

	//! Reads in fake parameters and assigns them (Requires the parameter to be included as a normal parameter as well)
	int ReadFakeDataPars(std::string parstring);

	//! Read in the samples so we can set up the free normalisation dials if required
	int ReadSamples(std::string sampleString);
	void SetupMinimizerFromXML();

	/*
	Setup Functions
	*/

	//! Setup the configuration given the arguments passed at the commandline and card file
	void SetupConfig();

	//! Setups up our custom RW engine with all the parameters passed in the card file
	void SetupRWEngine();

	//! Setups up the jointFCN.
	void SetupFCN();

	//! Sets up the minimizerObj for ROOT. there are cases where this is called repeatedly, e.g. If you are using a brute force scan before using Migrad.
	void SetupFitter(std::string routine);

	//! Set the current data histograms in each sample to the fake data.
	void SetFakeData();

	//! Setup the covariances with the correct dimensions. At the start this is either uncorrelated or merged given all the input covariances.
	//! At the end of the fit this produces the blank covariances which can then be filled by the minimizerObj with best fit covariances.
	void SetupCovariance();

	/*
	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);

	//! Get the current state of minimizerObj and fill it into currentVals and currentNorms
	void GetMinimizerState();

	//! Print current value
	void PrintState();

	//! Performs a fit routine where the input.maxevents is set to a much lower value to try and move closer to the best fit minimum.
	void LowStatRoutine(std::string routine);

	//! Perform a chi2 scan in 1D around the current point
	void Create1DScans();

	//! Perform a chi2 scan in 2D around the current point
	void Chi2Scan2D();

	//! Currently a placeholder NEEDS UPDATING
	void CreateContours();

	//! If any currentVals are close to the limits set them to the limit and fix them
	int FixAtLimit();

	//! 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 ThrowCovariance(bool uniformly);

	//! Given the covariance we currently have generate error bands by throwing the covariance.
	//! The FitPar config "error_uniform" defines whether to throw using the covariance or uniformly.
	//! The FitPar config "error_throws" defines how many throws are needed.
	//! Currently only supports TH1D plots.
	void GenerateErrorBands();

	/*
	Write Functions
	*/

	//! Write plots and TTrees listing the minimizerObj result of the fit to file
	void SaveMinimizerState();

	//! Save the sample plots for current MC
	//! dir if not empty forces plots to be saved in a subdirectory of outputfile
	void SaveCurrentState(std::string subdir="");

	//! Save starting predictions into a seperate folder
	void SaveNominal();

	//! Save predictions before the fit is ran into a seperate folder
	void SavePrefit();

	void SaveResults();
	/*
	MISC Functions
	*/

	//! Get previous fit status from a file
	Int_t GetStatus();

	/// Makes a histogram of likelihoods when throwing the data according to its statistics
	void ThrowDataToys();

	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
	ROOT::Math::Minimizer* fMinimizer;

	JointFCN* fSampleFCN;
	MinimizerFCN* fMinimizerFCN;
	ROOT::Math::Functor* fCallFunctor;

	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;

	+ struct mirror_param {
	+ double mirror_value;
	+ bool mirror_above;
	+ };
	+ std::map<std::string, mirror_param> fMirroredParams;
	+
	//! 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;

	TH2D* fCorrel;
	TH2D* fDecomp;
	TH2D* fCovar;

	TH2D* fCorFree;
	TH2D* fDecFree;
	TH2D* fCovFree;

	nuiskey fCompKey;

	};

	/! @} /
	#endif

File Metadata

Mime Type: text/x-diff
Expires: Thu, Apr 24, 6:37 AM (1 d, 21 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4887791
Default Alt Text: (125 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions