No OneTemporary
Actions

Size

209 KB

Subscribers

None

View Options

	diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx
	index a1ebaa7..d9e09cb 100644
	--- a/src/FitBase/Measurement1D.cxx
	+++ b/src/FitBase/Measurement1D.cxx
	@@ -1,1898 +1,1940 @@
	// Copyright 2016 L. Pickering, P. Stowell, R. Terri, C. Wilkinson, C. Wret

	/*******************************************************************************
	* This ile 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 "Measurement1D.h"

	//********************************************************************
	Measurement1D::Measurement1D(void) {
	//********************************************************************

	// XSec Scalings
	fScaleFactor = -1.0;
	fCurrentNorm = 1.0;

	// Histograms
	fDataHist = NULL;
	fDataTrue = NULL;

	fMCHist = NULL;
	fMCFine = NULL;
	fMCWeighted = NULL;

	fMaskHist = NULL;

	// Covar
	covar = NULL;
	fFullCovar = NULL;
	fShapeCovar = NULL;

	fCovar = NULL;
	fInvert = NULL;
	fDecomp = NULL;

	+ fResidualHist = NULL;
	+ fChi2LessBinHist = NULL;
	+
	// Fake Data
	fFakeDataInput = "";
	fFakeDataFile = NULL;

	// Options
	fDefaultTypes = "FIX/FULL/CHI2";
	fAllowedTypes =
	"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK/NOWIDTH";

	fIsFix = false;
	fIsShape = false;
	fIsFree = false;
	fIsDiag = false;
	fIsFull = false;
	fAddNormPen = false;
	fIsMask = false;
	fIsChi2SVD = false;
	fIsRawEvents = false;
	fIsNoWidth = false;
	fIsDifXSec = false;
	fIsEnu1D = false;

	// Inputs
	fInput = NULL;
	fRW = NULL;

	// Extra Histograms
	fMCHist_Modes = NULL;
	}

	//********************************************************************
	Measurement1D::~Measurement1D(void) {
	//********************************************************************

	if (fDataHist)
	delete fDataHist;
	if (fDataTrue)
	delete fDataTrue;
	if (fMCHist)
	delete fMCHist;
	if (fMCFine)
	delete fMCFine;
	if (fMCWeighted)
	delete fMCWeighted;
	if (fMaskHist)
	delete fMaskHist;
	if (covar)
	delete covar;
	if (fFullCovar)
	delete fFullCovar;
	if (fShapeCovar)
	delete fShapeCovar;
	if (fCovar)
	delete fCovar;
	if (fInvert)
	delete fInvert;
	if (fDecomp)
	delete fDecomp;
	+
	+ delete fResidualHist;
	+ delete fChi2LessBinHist;
	}

	//********************************************************************
	void Measurement1D::FinaliseSampleSettings() {
	//********************************************************************

	MeasurementBase::FinaliseSampleSettings();

	// Setup naming + renaming
	fName = fSettings.GetName();
	fSettings.SetS("originalname", fName);
	if (fSettings.Has("rename")) {
	fName = fSettings.GetS("rename");
	fSettings.SetS("name", fName);
	}

	// Setup all other options
	NUIS_LOG(SAM, "Finalising Sample Settings: " << fName);

	if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) {
	fIsRawEvents = true;
	- NUIS_LOG(SAM, "Found event rate measurement but using poisson likelihoods.");
	+ NUIS_LOG(SAM,
	+ "Found event rate measurement but using poisson likelihoods.");
	}

	if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) {
	fIsEnu1D = true;
	NUIS_LOG(SAM, "::" << fName << "::");
	- NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, "
	- << "not flux averaged!");
	+ NUIS_LOG(SAM,
	+ "Found XSec Enu measurement, applying flux integrated scaling, "
	+ << "not flux averaged!");
	}

	if (fIsEnu1D && fIsRawEvents) {
	NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
	- "really correct?!");
	+ "really correct?!");
	NUIS_ERR(FTL, "Check experiment constructor for " << fName
	- << " and correct this!");
	+ << " and correct this!");
	NUIS_ERR(FTL, "I live in " << __FILE__ << ":" << __LINE__);
	throw;
	}

	if (!fRW)
	fRW = FitBase::GetRW();
	if (!fInput and !fIsJoint)
	SetupInputs(fSettings.GetS("input"));

	// Setup options
	SetFitOptions(fDefaultTypes); // defaults
	SetFitOptions(fSettings.GetS("type")); // user specified

	EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min"));
	EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max"));

	if (fAddNormPen) {
	if (fNormError <= 0.0) {
	NUIS_ERR(FTL, "Norm error for class " << fName << " is 0.0!");
	NUIS_ERR(FTL, "If you want to use it please add fNormError=VAL");
	throw;
	}
	}
	}

	//********************************************************************
	void Measurement1D::CreateDataHistogram(int dimx, double *binx) {
	//********************************************************************

	if (fDataHist)
	delete fDataHist;

	fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(),
	(fSettings.GetFullTitles()).c_str(), dimx, binx);
	}

	//********************************************************************
	void Measurement1D::SetDataFromTextFile(std::string datafile) {
	//********************************************************************

	NUIS_LOG(SAM, "Reading data from text file: " << datafile);
	fDataHist = PlotUtils::GetTH1DFromFile(
	datafile, fSettings.GetName() + "_data", fSettings.GetFullTitles());
	}

	//********************************************************************
	void Measurement1D::SetDataFromRootFile(std::string datafile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM, "Reading data from root file: " << datafile << ";" << histname);
	fDataHist = PlotUtils::GetTH1DFromRootFile(datafile, histname);
	fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
	(fSettings.GetFullTitles()).c_str());

	return;
	};

	//********************************************************************
	void Measurement1D::SetEmptyData() {
	//********************************************************************

	fDataHist = new TH1D("EMPTY_DATA", "EMPTY_DATA", 1, 0.0, 1.0);
	}

	//********************************************************************
	void Measurement1D::SetPoissonErrors() {
	//********************************************************************

	if (!fDataHist) {
	NUIS_ERR(FTL, "Need a data hist to setup possion errors! ");
	NUIS_ERR(FTL, "Setup Data First!");
	throw;
	}

	for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
	fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1)));
	}
	}

	//********************************************************************
	void Measurement1D::SetCovarFromDiagonal(TH1D *data) {
	//********************************************************************

	if (!data and fDataHist) {
	data = fDataHist;
	}

	if (data) {
	NUIS_LOG(SAM, "Setting diagonal covariance for: " << data->GetName());
	fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	} else {
	NUIS_ABORT("No data input provided to set diagonal covar from!");
	}

	// if (!fIsDiag) {
	// ERR(FTL) << "SetCovarMatrixFromDiag called for measurement "
	// << "that is not set as diagonal." );
	// throw;
	// }
	}

	//********************************************************************
	void Measurement1D::SetCovarFromTextFile(std::string covfile, int dim) {
	//********************************************************************

	if (dim == -1) {
	dim = fDataHist->GetNbinsX();
	}

	NUIS_LOG(SAM, "Reading covariance from text file: " << covfile);
	fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement1D::SetCovarFromMultipleTextFiles(std::string covfiles,
	int dim) {
	//********************************************************************

	if (dim == -1) {
	dim = fDataHist->GetNbinsX();
	}

	std::vector<std::string> covList = GeneralUtils::ParseToStr(covfiles, ";");

	fFullCovar = new TMatrixDSym(dim);
	for (uint i = 0; i < covList.size(); ++i) {
	NUIS_LOG(SAM, "Reading covariance from text file: " << covList[i]);
	TMatrixDSym *temp_cov = StatUtils::GetCovarFromTextFile(covList[i], dim);
	(fFullCovar) += (temp_cov);
	delete temp_cov;
	}
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement1D::SetCovarFromRootFile(std::string covfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM,
	- "Reading covariance from text file: " << covfile << ";" << histname);
	+ "Reading covariance from text file: " << covfile << ";" << histname);
	fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname);
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement1D::SetCovarInvertFromTextFile(std::string covfile, int dim) {
	//********************************************************************

	if (dim == -1) {
	dim = fDataHist->GetNbinsX();
	}

	NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile);
	covar = StatUtils::GetCovarFromTextFile(covfile, dim);
	fFullCovar = StatUtils::GetInvert(covar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement1D::SetCovarInvertFromRootFile(std::string covfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile << ";"
	- << histname);
	+ << histname);
	covar = StatUtils::GetCovarFromRootFile(covfile, histname);
	fFullCovar = StatUtils::GetInvert(covar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement1D::SetCorrelationFromTextFile(std::string covfile, int dim) {
	//********************************************************************

	if (dim == -1)
	dim = fDataHist->GetNbinsX();
	- NUIS_LOG(SAM,
	- "Reading data correlations from text file: " << covfile << ";" << dim);
	+ NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";"
	+ << dim);
	TMatrixDSym *correlation = StatUtils::GetCovarFromTextFile(covfile, dim);

	if (!fDataHist) {
	NUIS_ABORT("Trying to set correlations from text file but there is no "
	- "data to build it from. \n"
	- << "In constructor make sure data is set before "
	- "SetCorrelationFromTextFile is called. \n");
	+ "data to build it from. \n"
	+ << "In constructor make sure data is set before "
	+ "SetCorrelationFromTextFile is called. \n");
	}

	// Fill covar from data errors and correlations
	fFullCovar = new TMatrixDSym(dim);
	for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
	for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
	(fFullCovar)(i, j) = (correlation)(i, j) *
	fDataHist->GetBinError(i + 1) *
	fDataHist->GetBinError(j + 1) * 1.E76;
	}
	}

	// Fill other covars.
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);

	delete correlation;
	}

	//********************************************************************
	void Measurement1D::SetCorrelationFromMultipleTextFiles(std::string corrfiles,
	int dim) {
	//********************************************************************

	if (dim == -1) {
	dim = fDataHist->GetNbinsX();
	}

	std::vector<std::string> corrList = GeneralUtils::ParseToStr(corrfiles, ";");

	fFullCovar = new TMatrixDSym(dim);
	for (uint i = 0; i < corrList.size(); ++i) {
	NUIS_LOG(SAM, "Reading covariance from text file: " << corrList[i]);
	TMatrixDSym *temp_cov = StatUtils::GetCovarFromTextFile(corrList[i], dim);

	for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
	for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
	(temp_cov)(i, j) = (temp_cov)(i, j) * fDataHist->GetBinError(i + 1) *
	fDataHist->GetBinError(j + 1) * 1.E76;
	}
	}

	(fFullCovar) += (temp_cov);
	delete temp_cov;
	}
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement1D::SetCorrelationFromRootFile(std::string covfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";"
	- << histname);
	+ << histname);
	TMatrixDSym *correlation = StatUtils::GetCovarFromRootFile(covfile, histname);

	if (!fDataHist) {
	NUIS_ABORT("Trying to set correlations from text file but there is no "
	- "data to build it from. \n"
	- << "In constructor make sure data is set before "
	- "SetCorrelationFromTextFile is called. \n");
	+ "data to build it from. \n"
	+ << "In constructor make sure data is set before "
	+ "SetCorrelationFromTextFile is called. \n");
	}

	// Fill covar from data errors and correlations
	fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
	for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
	for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
	(fFullCovar)(i, j) = (correlation)(i, j) *
	fDataHist->GetBinError(i + 1) *
	fDataHist->GetBinError(j + 1) * 1.E76;
	}
	}

	// Fill other covars.
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);

	delete correlation;
	}

	//********************************************************************
	void Measurement1D::SetCholDecompFromTextFile(std::string covfile, int dim) {
	//********************************************************************

	if (dim == -1) {
	dim = fDataHist->GetNbinsX();
	}

	NUIS_LOG(SAM, "Reading cholesky from text file: " << covfile);
	TMatrixD *temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim);

	TMatrixD trans = (TMatrixD )temp->Clone();
	trans->T();
	(trans) = (*temp);

	fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);

	delete temp;
	delete trans;
	}

	//********************************************************************
	void Measurement1D::SetCholDecompFromRootFile(std::string covfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM, "Reading cholesky decomp from root file: " << covfile << ";"
	- << histname);
	+ << histname);
	TMatrixD *temp = StatUtils::GetMatrixFromRootFile(covfile, histname);

	TMatrixD trans = (TMatrixD )temp->Clone();
	trans->T();
	(trans) = (*temp);

	fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), "");
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);

	delete temp;
	delete trans;
	}

	void Measurement1D::SetShapeCovar() {

	// Return if this is missing any pre-requisites
	if (!fFullCovar)
	return;
	if (!fDataHist)
	return;

	// Also return if it's bloody stupid under the circumstances
	if (fIsDiag)
	return;

	fShapeCovar = StatUtils::ExtractShapeOnlyCovar(fFullCovar, fDataHist);
	return;
	}

	//********************************************************************
	void Measurement1D::ScaleData(double scale) {
	//********************************************************************
	fDataHist->Scale(scale);
	}

	//********************************************************************
	void Measurement1D::ScaleDataErrors(double scale) {
	//********************************************************************
	for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
	fDataHist->SetBinError(i + 1, fDataHist->GetBinError(i + 1) * scale);
	}
	}

	//********************************************************************
	void Measurement1D::ScaleCovar(double scale) {
	//********************************************************************
	(fFullCovar) = scale;
	(covar) = 1.0 / scale;
	(fDecomp) = sqrt(scale);
	}

	//********************************************************************
	void Measurement1D::SetBinMask(std::string maskfile) {
	//********************************************************************

	if (!fIsMask)
	return;
	NUIS_LOG(SAM, "Reading bin mask from file: " << maskfile);

	// Create a mask histogram with dim of data
	int nbins = fDataHist->GetNbinsX();
	fMaskHist = new TH1I((fSettings.GetName() + "_BINMASK").c_str(),
	(fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(),
	nbins, 0, nbins);
	std::string line;
	std::ifstream mask(maskfile.c_str(), std::ifstream::in);

	if (!mask.is_open()) {
	NUIS_ABORT("Cannot find mask file.");
	}

	while (std::getline(mask >> std::ws, line, '\n')) {
	std::vector<int> entries = GeneralUtils::ParseToInt(line, " ");

	// Skip lines with poorly formatted lines
	if (entries.size() < 2) {
	- NUIS_LOG(WRN, "Measurement1D::SetBinMask(), couldn't parse line: " << line);
	+ NUIS_LOG(WRN,
	+ "Measurement1D::SetBinMask(), couldn't parse line: " << line);
	continue;
	}

	// The first index should be the bin number, the second should be the mask
	// value.
	int val = 0;
	if (entries[1] > 0)
	val = 1;
	fMaskHist->SetBinContent(entries[0], val);
	}

	// Apply masking by setting masked data bins to zero
	PlotUtils::MaskBins(fDataHist, fMaskHist);

	return;
	}

	//********************************************************************
	void Measurement1D::FinaliseMeasurement() {
	//********************************************************************

	NUIS_LOG(SAM, "Finalising Measurement: " << fName);

	if (fSettings.GetB("onlymc")) {
	if (fDataHist)
	delete fDataHist;
	fDataHist = new TH1D("empty_data", "empty_data", 1, 0.0, 1.0);
	}

	// Make sure data is setup
	if (!fDataHist) {
	NUIS_ABORT("No data has been setup inside " << fName << " constructor!");
	}

	// Make sure covariances are setup
	if (!fFullCovar) {
	fIsDiag = true;
	SetCovarFromDiagonal(fDataHist);
	}

	if (!covar) {
	covar = StatUtils::GetInvert(fFullCovar);
	}

	if (!fDecomp) {
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	// Push the diagonals of fFullCovar onto the data histogram
	// Comment this out until the covariance/data scaling is consistent!
	StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1E-38);

	// If shape only, set covar and fDecomp using the shape-only matrix (if set)
	if (fIsShape && fShapeCovar && FitPar::Config().GetParB("UseShapeCovar")) {
	if (covar)
	delete covar;
	covar = StatUtils::GetInvert(fShapeCovar);
	if (fDecomp)
	delete fDecomp;
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	// Setup fMCHist from data
	fMCHist = (TH1D *)fDataHist->Clone();
	fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
	(fSettings.GetFullTitles()).c_str());
	fMCHist->Reset();

	// Setup fMCFine
	fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX() * 8,
	fMCHist->GetBinLowEdge(1),
	fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1));
	fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(),
	(fSettings.GetFullTitles()).c_str());
	fMCFine->Reset();

	// Setup MC Stat
	fMCStat = (TH1D *)fMCHist->Clone();
	fMCStat->Reset();

	// Search drawopts for possible types to include by default
	std::string drawopts = FitPar::Config().GetParS("drawopts");
	if (drawopts.find("MODES") != std::string::npos) {
	fMCHist_Modes = new TrueModeStack((fSettings.GetName() + "_MODES").c_str(),
	("True Channels"), fMCHist);
	SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write);
	}

	// Setup bin masks using sample name
	if (fIsMask) {

	std::string curname = fName;
	std::string origname = fSettings.GetS("originalname");

	// Check rename.mask
	std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask");

	// Check origname.mask
	if (maskloc.empty())
	maskloc = FitPar::Config().GetParDIR(origname + ".mask");

	// Check database
	if (maskloc.empty()) {
	maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask";
	}

	// Setup Bin Mask
	SetBinMask(maskloc);
	}

	if (fScaleFactor < 0) {
	NUIS_ERR(FTL, "I found a negative fScaleFactor in " << __FILE__ << ":"
	- << __LINE__);
	+ << __LINE__);
	NUIS_ERR(FTL, "fScaleFactor = " << fScaleFactor);
	NUIS_ERR(FTL, "EXITING");
	throw;
	}

	// Create and fill Weighted Histogram
	if (!fMCWeighted) {

	fMCWeighted = (TH1D *)fMCHist->Clone();
	fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(),
	(fName + "_MCWGHTS" + fPlotTitles).c_str());
	fMCWeighted->GetYaxis()->SetTitle("Weighted Events");
	}
	}

	//********************************************************************
	void Measurement1D::SetFitOptions(std::string opt) {
	//********************************************************************

	// Do nothing if default given
	if (opt == "DEFAULT")
	return;

	// CHECK Conflicting Fit Options
	std::vector<std::string> fit_option_allow =
	GeneralUtils::ParseToStr(fAllowedTypes, "/");

	for (UInt_t i = 0; i < fit_option_allow.size(); i++) {
	std::vector<std::string> fit_option_section =
	GeneralUtils::ParseToStr(fit_option_allow.at(i), ",");

	bool found_option = false;

	for (UInt_t j = 0; j < fit_option_section.size(); j++) {
	std::string av_opt = fit_option_section.at(j);

	if (!found_option and opt.find(av_opt) != std::string::npos) {
	found_option = true;

	} else if (found_option and opt.find(av_opt) != std::string::npos) {
	- NUIS_ABORT("ERROR: Conflicting fit options provided: "
	- << opt << std::endl
	- << "Conflicting group = " << fit_option_section.at(i)
	- << std::endl
	- << "You should only supply one of these options in card file.");
	+ NUIS_ABORT(
	+ "ERROR: Conflicting fit options provided: "
	+ << opt << std::endl
	+ << "Conflicting group = " << fit_option_section.at(i) << std::endl
	+ << "You should only supply one of these options in card file.");
	}
	}
	}

	// Check all options are allowed
	std::vector<std::string> fit_options_input =
	GeneralUtils::ParseToStr(opt, "/");
	for (UInt_t i = 0; i < fit_options_input.size(); i++) {
	if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) {
	NUIS_ERR(WRN, "ERROR: Fit Option '"
	- << fit_options_input.at(i)
	- << "' Provided is not allowed for this measurement.");
	+ << fit_options_input.at(i)
	+ << "' Provided is not allowed for this measurement.");
	NUIS_ERR(WRN, "Fit Options should be provided as a '/' seperated list "
	- "(e.g. FREE/DIAG/NORM)");
	+ "(e.g. FREE/DIAG/NORM)");
	NUIS_ABORT("Available options for " << fName << " are '" << fAllowedTypes
	- << "'");
	+ << "'");
	}
	}

	// Set TYPE
	fFitType = opt;

	// FIX,SHAPE,FREE
	if (opt.find("FIX") != std::string::npos) {
	fIsFree = fIsShape = false;
	fIsFix = true;
	} else if (opt.find("SHAPE") != std::string::npos) {
	fIsFree = fIsFix = false;
	fIsShape = true;
	} else if (opt.find("FREE") != std::string::npos) {
	fIsFix = fIsShape = false;
	fIsFree = true;
	}

	// DIAG,FULL (or default to full)
	if (opt.find("DIAG") != std::string::npos) {
	fIsDiag = true;
	fIsFull = false;
	} else if (opt.find("FULL") != std::string::npos) {
	fIsDiag = false;
	fIsFull = true;
	}

	// CHI2/LL (OTHERS?)
	if (opt.find("LOG") != std::string::npos) {
	fIsChi2 = false;

	NUIS_ERR(FTL, "No other LIKELIHOODS properly supported!");
	NUIS_ERR(FTL, "Try to use a chi2!");
	throw;

	} else {
	fIsChi2 = true;
	}

	// EXTRAS
	if (opt.find("RAW") != std::string::npos)
	fIsRawEvents = true;
	if (opt.find("NOWIDTH") != std::string::npos)
	fIsNoWidth = true;
	if (opt.find("DIF") != std::string::npos)
	fIsDifXSec = true;
	if (opt.find("ENU1D") != std::string::npos)
	fIsEnu1D = true;
	if (opt.find("NORM") != std::string::npos)
	fAddNormPen = true;
	if (opt.find("MASK") != std::string::npos)
	fIsMask = true;

	return;
	};

	//********************************************************************
	void Measurement1D::SetSmearingMatrix(std::string smearfile, int truedim,
	int recodim) {
	//********************************************************************

	// The smearing matrix describes the migration from true bins (rows) to reco
	// bins (columns)
	// Counter over the true bins!
	int row = 0;

	std::string line;
	std::ifstream smear(smearfile.c_str(), std::ifstream::in);

	// Note that the smearing matrix may be rectangular.
	fSmearMatrix = new TMatrixD(truedim, recodim);

	if (smear.is_open()) {
	NUIS_LOG(SAM, "Reading smearing matrix from file: " << smearfile);
	} else {
	NUIS_ABORT("Smearing matrix provided is incorrect: " << smearfile);
	}

	while (std::getline(smear >> std::ws, line, '\n')) {
	int column = 0;

	std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
	for (std::vector<double>::iterator iter = entries.begin();
	iter != entries.end(); iter++) {
	(fSmearMatrix)(row, column) = (iter) / 100.; // Convert to fraction from
	// percentage (this may not be
	// general enough)
	column++;
	}
	row++;
	}
	return;
	}

	//********************************************************************
	void Measurement1D::ApplySmearingMatrix() {
	//********************************************************************

	if (!fSmearMatrix) {
	NUIS_ERR(WRN,
	- fName << ": attempted to apply smearing matrix, but none was set");
	+ fName << ": attempted to apply smearing matrix, but none was set");
	return;
	}

	TH1D unsmeared = (TH1D )fMCHist->Clone();
	TH1D smeared = (TH1D )fMCHist->Clone();
	smeared->Reset();

	// Loop over reconstructed bins
	// true = row; reco = column
	for (int rbin = 0; rbin < fSmearMatrix->GetNcols(); ++rbin) {
	// Sum up the constributions from all true bins
	double rBinVal = 0;

	// Loop over true bins
	for (int tbin = 0; tbin < fSmearMatrix->GetNrows(); ++tbin) {
	rBinVal +=
	(fSmearMatrix)(tbin, rbin) unsmeared->GetBinContent(tbin + 1);
	}
	smeared->SetBinContent(rbin + 1, rBinVal);
	}
	fMCHist = (TH1D *)smeared->Clone();

	return;
	}

	/*
	Reconfigure LOOP
	*/
	//********************************************************************
	void Measurement1D::ResetAll() {
	//********************************************************************

	fMCHist->Reset();
	fMCFine->Reset();
	fMCStat->Reset();

	return;
	};

	//********************************************************************
	void Measurement1D::FillHistograms() {
	//********************************************************************

	if (Signal) {

	NUIS_LOG(DEB, "Fill MCHist: " << fXVar << ", " << Weight);

	fMCHist->Fill(fXVar, Weight);
	fMCFine->Fill(fXVar, Weight);
	fMCStat->Fill(fXVar, 1.0);

	if (fMCHist_Modes)
	fMCHist_Modes->Fill(Mode, fXVar, Weight);
	}

	return;
	};

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

	// Fill MCWeighted;
	// for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
	// fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1));
	// fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1));
	// }

	// Setup Stat ratios for MC and MC Fine
	double *statratio = new double[fMCHist->GetNbinsX()];
	for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
	if (fMCHist->GetBinContent(i + 1) != 0) {
	statratio[i] =
	fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1);
	} else {
	statratio[i] = 0.0;
	}
	}

	double *statratiofine = new double[fMCFine->GetNbinsX()];
	for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
	if (fMCFine->GetBinContent(i + 1) != 0) {
	statratiofine[i] =
	fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1);
	} else {
	statratiofine[i] = 0.0;
	}
	}

	// Scaling for raw event rates
	if (fIsRawEvents) {
	double datamcratio = fDataHist->Integral() / fMCHist->Integral();

	fMCHist->Scale(datamcratio);
	fMCFine->Scale(datamcratio);

	if (fMCHist_Modes)
	fMCHist_Modes->Scale(datamcratio);

	// Scaling for XSec as function of Enu
	} else if (fIsEnu1D) {

	PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(),
	GetEventHistogram(), fScaleFactor, fNEvents);
	PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(),
	GetEventHistogram(), fScaleFactor, fNEvents);

	if (fMCHist_Modes) {
	// Loop over the modes
	fMCHist_Modes->FluxUnfold(GetFluxHistogram(), GetEventHistogram(),
	fScaleFactor, fNEvents);
	// PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(),
	// GetEventHistogram(), fScaleFactor,
	// fNEvents);
	}

	} else if (fIsNoWidth) {
	fMCHist->Scale(fScaleFactor);
	fMCFine->Scale(fScaleFactor);
	if (fMCHist_Modes)
	fMCHist_Modes->Scale(fScaleFactor);
	// Any other differential scaling
	} else {
	fMCHist->Scale(fScaleFactor, "width");
	fMCFine->Scale(fScaleFactor, "width");

	if (fMCHist_Modes)
	fMCHist_Modes->Scale(fScaleFactor, "width");
	}

	// Proper error scaling - ROOT Freaks out with xsec weights sometimes
	for (int i = 0; i < fMCStat->GetNbinsX(); i++) {
	fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]);
	}

	for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
	fMCFine->SetBinError(i + 1,
	fMCFine->GetBinContent(i + 1) * statratiofine[i]);
	}

	// Clean up
	delete[] statratio;
	delete[] statratiofine;

	return;
	};

	//********************************************************************
	void Measurement1D::ApplyNormScale(double norm) {
	//********************************************************************

	fCurrentNorm = norm;

	fMCHist->Scale(1.0 / norm);
	fMCFine->Scale(1.0 / norm);

	return;
	};

	/*
	Statistic Functions - Outsources to StatUtils
	*/

	//********************************************************************
	int Measurement1D::GetNDOF() {
	//********************************************************************
	int ndof = fDataHist->GetNbinsX();
	if (fMaskHist and fIsMask)
	ndof -= fMaskHist->Integral();
	return ndof;
	}

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

	// If this is for a ratio, there is no data histogram to compare to!
	if (fNoData \|\| !fDataHist)
	return 0.;

	// Apply Masking to MC if Required.
	if (fIsMask and fMaskHist) {
	PlotUtils::MaskBins(fMCHist, fMaskHist);
	}

	// Sort Shape Scaling
	double scaleF = 0.0;
	// TODO Include !fIsRawEvents
	if (fIsShape) {
	if (fMCHist->Integral(1, fMCHist->GetNbinsX(), "width")) {
	scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX(), "width") /
	fMCHist->Integral(1, fMCHist->GetNbinsX(), "width");
	fMCHist->Scale(scaleF);
	fMCFine->Scale(scaleF);
	}
	}

	// Likelihood Calculation
	double stat = 0.;
	if (fIsChi2) {

	if (fIsRawEvents) {
	stat = StatUtils::GetChi2FromEventRate(fDataHist, fMCHist, fMaskHist);
	} else if (fIsDiag) {
	stat = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMaskHist);
	} else if (!fIsDiag and !fIsRawEvents) {
	stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist);
	+
	+ stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist, 1,
	+ 1E76, fResidualHist);
	+ if (fChi2LessBinHist) {
	+ TH1I *binmask = new TH1I("mask", "", fDataHist->GetNbinsX(), 0,
	+ fDataHist->GetNbinsX());
	+ binmask->SetDirectory(NULL);
	+ for (int xi = 0; xi < fDataHist->GetNbinsX(); ++xi) {
	+ binmask->Reset();
	+ binmask->SetBinContent(xi + 1, 1);
	+ fChi2LessBinHist->SetBinContent(
	+ xi + 1,
	+ StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, binmask));
	+ }
	+ delete binmask;
	+ }
	}
	}

	// Sort Penalty Terms
	if (fAddNormPen) {
	double penalty =
	(1. - fCurrentNorm) * (1. - fCurrentNorm) / (fNormError * fNormError);

	stat += penalty;
	}

	// Return to normal scaling
	if (fIsShape) { // and !FitPar::Config().GetParB("saveshapescaling")) {
	fMCHist->Scale(1. / scaleF);
	fMCFine->Scale(1. / scaleF);
	}

	fLikelihood = stat;

	return stat;
	}

	/*
	Fake Data Functions
	*/
	//********************************************************************
	void Measurement1D::SetFakeDataValues(std::string fakeOption) {
	//********************************************************************

	// Setup original/datatrue
	TH1D tempdata = (TH1D )fDataHist->Clone();

	if (!fIsFakeData) {
	fIsFakeData = true;

	// Make a copy of the original data histogram.
	if (!fDataOrig)
	fDataOrig = (TH1D *)fDataHist->Clone((fName + "_data_original").c_str());

	} else {
	ResetFakeData();
	}

	// Setup Inputs
	fFakeDataInput = fakeOption;
	NUIS_LOG(SAM, "Setting fake data from : " << fFakeDataInput);

	// From MC
	if (fFakeDataInput.compare("MC") == 0) {
	fDataHist = (TH1D *)fMCHist->Clone((fName + "_MC").c_str());

	// Fake File
	} else {
	if (!fFakeDataFile)
	fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ");
	fDataHist = (TH1D *)fFakeDataFile->Get((fName + "_MC").c_str());
	}

	// Setup Data Hist
	fDataHist->SetNameTitle((fName + "_FAKE").c_str(),
	(fName + fPlotTitles).c_str());

	// Replace Data True
	if (fDataTrue)
	delete fDataTrue;
	fDataTrue = (TH1D *)fDataHist->Clone();
	fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(),
	(fName + fPlotTitles).c_str());

	// Make a new covariance for fake data hist.
	int nbins = fDataHist->GetNbinsX();
	double alpha_i = 0.0;
	double alpha_j = 0.0;

	for (int i = 0; i < nbins; i++) {
	for (int j = 0; j < nbins; j++) {
	alpha_i =
	fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1);
	alpha_j =
	fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1);

	(fFullCovar)(i, j) = alpha_i alpha_j * (*fFullCovar)(i, j);
	}
	}

	// Setup Covariances
	if (covar)
	delete covar;
	covar = StatUtils::GetInvert(fFullCovar);

	if (fDecomp)
	delete fDecomp;
	fDecomp = StatUtils::GetInvert(fFullCovar);

	delete tempdata;

	return;
	};

	//********************************************************************
	void Measurement1D::ResetFakeData() {
	//********************************************************************

	if (fIsFakeData) {
	if (fDataHist)
	delete fDataHist;
	fDataHist =
	(TH1D *)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str());
	}
	}

	//********************************************************************
	void Measurement1D::ResetData() {
	//********************************************************************

	if (fIsFakeData) {
	if (fDataHist)
	delete fDataHist;
	fDataHist =
	(TH1D *)fDataOrig->Clone((fSettings.GetName() + "_data").c_str());
	}

	fIsFakeData = false;
	}

	//********************************************************************
	void Measurement1D::ThrowCovariance() {
	//********************************************************************

	// Take a fDecomposition and use it to throw the current dataset.
	// Requires fDataTrue also be set incase used repeatedly.

	if (!fDataTrue)
	fDataTrue = (TH1D *)fDataHist->Clone();
	if (fDataHist)
	delete fDataHist;
	fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);

	return;
	};

	//********************************************************************
	void Measurement1D::ThrowDataToy() {
	//********************************************************************
	if (!fDataTrue)
	fDataTrue = (TH1D *)fDataHist->Clone();
	if (fMCHist)
	delete fMCHist;
	fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
	}

	/*
	Access Functions
	*/

	//********************************************************************
	TH1D *Measurement1D::GetMCHistogram() {
	//********************************************************************

	if (!fMCHist)
	return fMCHist;

	std::ostringstream chi2;
	chi2 << std::setprecision(5) << this->GetLikelihood();

	int linecolor = kRed;
	int linestyle = 1;
	int linewidth = 1;

	int fillcolor = 0;
	int fillstyle = 1001;

	// if (fSettings.Has("linecolor")) linecolor = fSettings.GetI("linecolor");
	// if (fSettings.Has("linestyle")) linestyle = fSettings.GetI("linestyle");
	// if (fSettings.Has("linewidth")) linewidth = fSettings.GetI("linewidth");

	// if (fSettings.Has("fillcolor")) fillcolor = fSettings.GetI("fillcolor");
	// if (fSettings.Has("fillstyle")) fillstyle = fSettings.GetI("fillstyle");

	fMCHist->SetTitle(chi2.str().c_str());

	fMCHist->SetLineColor(linecolor);
	fMCHist->SetLineStyle(linestyle);
	fMCHist->SetLineWidth(linewidth);

	fMCHist->SetFillColor(fillcolor);
	fMCHist->SetFillStyle(fillstyle);

	return fMCHist;
	};

	//********************************************************************
	TH1D *Measurement1D::GetDataHistogram() {
	//********************************************************************

	if (!fDataHist)
	return fDataHist;

	int datacolor = kBlack;
	int datastyle = 1;
	int datawidth = 1;

	// if (fSettings.Has("datacolor")) datacolor = fSettings.GetI("datacolor");
	// if (fSettings.Has("datastyle")) datastyle = fSettings.GetI("datastyle");
	// if (fSettings.Has("datawidth")) datawidth = fSettings.GetI("datawidth");

	fDataHist->SetLineColor(datacolor);
	fDataHist->SetLineWidth(datawidth);
	fDataHist->SetMarkerStyle(datastyle);

	return fDataHist;
	};

	/*
	Write Functions
	*/

	// Save all the histograms at once
	//********************************************************************
	void Measurement1D::Write(std::string drawOpt) {
	//********************************************************************

	// Get Draw Options
	drawOpt = FitPar::Config().GetParS("drawopts");

	// Write Settigns
	if (drawOpt.find("SETTINGS") != std::string::npos) {
	fSettings.Set("#chi^{2}", fLikelihood);
	fSettings.Set("NDOF", this->GetNDOF());
	fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF());
	fSettings.Write();
	}

	// Write Data/MC
	if (drawOpt.find("DATA") != std::string::npos)
	GetDataList().at(0)->Write();
	if (drawOpt.find("MC") != std::string::npos) {
	GetMCList().at(0)->Write();
	if ((fEvtRateScaleFactor != 0xdeadbeef) && GetMCList().at(0)) {
	TH1D PredictedEvtRate = static_cast<TH1D >(GetMCList().at(0)->Clone());
	PredictedEvtRate->Scale(fEvtRateScaleFactor);
	PredictedEvtRate->GetYaxis()->SetTitle("Predicted event rate");
	PredictedEvtRate->Write();
	}
	}

	// Write Fine Histogram
	if (drawOpt.find("FINE") != std::string::npos)
	GetFineList().at(0)->Write();

	// Write Weighted Histogram
	if (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted)
	fMCWeighted->Write();

	// Save Flux/Evt if no event manager
	if (!FitPar::Config().GetParB("EventManager")) {
	if (drawOpt.find("FLUX") != std::string::npos && GetFluxHistogram())
	GetFluxHistogram()->Write();

	if (drawOpt.find("EVT") != std::string::npos && GetEventHistogram())
	GetEventHistogram()->Write();

	if (drawOpt.find("XSEC") != std::string::npos && GetEventHistogram())
	GetXSecHistogram()->Write();
	}

	// Write Mask
	if (fIsMask && (drawOpt.find("MASK") != std::string::npos)) {
	fMaskHist->Write();
	}

	// Write Covariances
	if (drawOpt.find("COV") != std::string::npos && fFullCovar) {
	PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName())->Write();
	}

	if (drawOpt.find("INVCOV") != std::string::npos && covar) {
	PlotUtils::GetInvCovarPlot(covar, fSettings.GetName())->Write();
	}

	if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) {
	PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName())->Write();
	}

	// // Likelihood residual plots
	// if (drawOpt.find("RESIDUAL") != std::string::npos) {
	// WriteResidualPlots();
	// }

	// Ratio and Shape Plots
	if (drawOpt.find("RATIO") != std::string::npos) {
	WriteRatioPlot();
	}

	if (drawOpt.find("SHAPE") != std::string::npos) {
	WriteShapePlot();
	if (drawOpt.find("RATIO") != std::string::npos)
	WriteShapeRatioPlot();
	}

	// // RATIO
	// if (drawOpt.find("CANVMC") != std::string::npos) {
	// TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist);
	// c1->Write();
	// delete c1;
	// }

	// // PDG
	// if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) {
	// TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes);
	// c2->Write();
	// delete c2;
	// }

	+ if (fIsChi2 && !fIsDiag) {
	+ fResidualHist = (TH1D *)fMCHist->Clone((fName + "_RESIDUAL").c_str());
	+ fResidualHist->GetYaxis()->SetTitle("#Delta#chi^{2}");
	+ fResidualHist->Reset();
	+
	+ fChi2LessBinHist =
	+ (TH1D *)fMCHist->Clone((fName + "_Chi2NMinusOne").c_str());
	+ fChi2LessBinHist->GetYaxis()->SetTitle("Total #chi^{2} without bin_{i}");
	+ fChi2LessBinHist->Reset();
	+
	+ (void)GetLikelihood();
	+
	+ fResidualHist->Write();
	+ fChi2LessBinHist->Write();
	+ }
	+
	// Write Extra Histograms
	AutoWriteExtraTH1();
	WriteExtraHistograms();

	// Returning
	NUIS_LOG(SAM, "Written Histograms: " << fName);
	return;
	}

	//********************************************************************
	void Measurement1D::WriteRatioPlot() {
	//********************************************************************

	// Setup mc data ratios
	TH1D dataRatio = (TH1D )fDataHist->Clone((fName + "_data_RATIO").c_str());
	TH1D mcRatio = (TH1D )fMCHist->Clone((fName + "_MC_RATIO").c_str());

	// Extra MC Data Ratios
	for (int i = 0; i < mcRatio->GetNbinsX(); i++) {

	dataRatio->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) /
	fMCHist->GetBinContent(i + 1));
	dataRatio->SetBinError(i + 1, fDataHist->GetBinError(i + 1) /
	fMCHist->GetBinContent(i + 1));

	mcRatio->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1) /
	fMCHist->GetBinContent(i + 1));
	mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) /
	fMCHist->GetBinContent(i + 1));
	}

	// Write ratios
	mcRatio->Write();
	dataRatio->Write();

	delete mcRatio;
	delete dataRatio;
	}

	//********************************************************************
	void Measurement1D::WriteShapePlot() {
	//********************************************************************

	TH1D mcShape = (TH1D )fMCHist->Clone((fName + "_MC_SHAPE").c_str());

	TH1D dataShape = (TH1D )fDataHist->Clone((fName + "_data_SHAPE").c_str());
	// Don't check error
	if (fShapeCovar)
	StatUtils::SetDataErrorFromCov(dataShape, fShapeCovar, 1E-38, false);

	double shapeScale = 1.0;
	if (fIsRawEvents) {
	shapeScale = fDataHist->Integral() / fMCHist->Integral();
	} else {
	shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
	}

	mcShape->Scale(shapeScale);

	std::stringstream ss;
	ss << shapeScale;
	mcShape->SetTitle(ss.str().c_str());

	mcShape->SetLineWidth(3);
	mcShape->SetLineStyle(7);

	mcShape->Write();
	dataShape->Write();

	delete mcShape;
	}

	//********************************************************************
	void Measurement1D::WriteShapeRatioPlot() {
	//********************************************************************

	// Get a mcshape histogram
	TH1D mcShape = (TH1D )fMCHist->Clone((fName + "_MC_SHAPE").c_str());

	double shapeScale = 1.0;
	if (fIsRawEvents) {
	shapeScale = fDataHist->Integral() / fMCHist->Integral();
	} else {
	shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
	}

	mcShape->Scale(shapeScale);

	// Create shape ratio histograms
	TH1D *mcShapeRatio =
	(TH1D *)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str());
	TH1D *dataShapeRatio =
	(TH1D *)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str());

	// Divide the histograms
	mcShapeRatio->Divide(mcShape);
	dataShapeRatio->Divide(mcShape);

	// Colour the shape ratio plots
	mcShapeRatio->SetLineWidth(3);
	mcShapeRatio->SetLineStyle(7);

	mcShapeRatio->Write();
	dataShapeRatio->Write();

	delete mcShapeRatio;
	delete dataShapeRatio;
	}

	//// CRAP TO BE REMOVED

	//********************************************************************
	void Measurement1D::SetupMeasurement(std::string inputfile, std::string type,
	FitWeight *rw, std::string fkdt) {
	//********************************************************************

	nuiskey samplekey = Config::CreateKey("sample");
	samplekey.Set("name", fName);
	samplekey.Set("type", type);
	samplekey.Set("input", inputfile);
	fSettings = LoadSampleSettings(samplekey);

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

	// Check if name contains Evt, indicating that it is a raw number of events
	// measurements and should thus be treated as once
	fIsRawEvents = false;
	if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) {
	fIsRawEvents = true;
	NUIS_LOG(SAM, "Found event rate measurement but fIsRawEvents == false!");
	NUIS_LOG(SAM, "Overriding this and setting fIsRawEvents == true!");
	}

	fIsEnu1D = false;
	if (fName.find("XSec_1DEnu") != std::string::npos) {
	fIsEnu1D = true;
	NUIS_LOG(SAM, "::" << fName << "::");
	- NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, "
	- "not flux averaged!");
	+ NUIS_LOG(SAM,
	+ "Found XSec Enu measurement, applying flux integrated scaling, "
	+ "not flux averaged!");
	}

	if (fIsEnu1D && fIsRawEvents) {
	NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
	- "really correct?!");
	+ "really correct?!");
	NUIS_ERR(FTL, "Check experiment constructor for " << fName
	- << " and correct this!");
	+ << " and correct this!");
	NUIS_ERR(FTL, "I live in " << __FILE__ << ":" << __LINE__);
	throw;
	}

	fRW = rw;

	if (!fInput and !fIsJoint)
	SetupInputs(inputfile);

	// Set Default Options
	SetFitOptions(fDefaultTypes);

	// Set Passed Options
	SetFitOptions(type);

	// Still adding support for flat flux inputs
	// // Set Enu Flux Scaling
	// if (isFlatFluxFolding) this->Input()->ApplyFluxFolding(
	// this->defaultFluxHist );

	// FinaliseMeasurement();
	}

	//********************************************************************
	void Measurement1D::SetupDefaultHist() {
	//********************************************************************

	// Setup fMCHist
	fMCHist = (TH1D *)fDataHist->Clone();
	fMCHist->SetNameTitle((fName + "_MC").c_str(),
	(fName + "_MC" + fPlotTitles).c_str());

	// Setup fMCFine
	Int_t nBins = fMCHist->GetNbinsX();
	fMCFine = new TH1D(
	(fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(),
	nBins * 6, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(nBins + 1));

	fMCStat = (TH1D *)fMCHist->Clone();
	fMCStat->Reset();

	fMCHist->Reset();
	fMCFine->Reset();

	// Setup the NEUT Mode Array
	PlotUtils::CreateNeutModeArray((TH1D )fMCHist, (TH1 *)fMCHist_PDG);
	PlotUtils::ResetNeutModeArray((TH1 **)fMCHist_PDG);

	// Setup bin masks using sample name
	if (fIsMask) {
	std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask");
	if (maskloc.empty()) {
	maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask";
	}

	SetBinMask(maskloc);
	}

	fMCHist_Modes =
	new TrueModeStack((fName + "_MODES").c_str(), ("True Channels"), fMCHist);
	SetAutoProcessTH1(fMCHist_Modes, kCMD_Reset, kCMD_Norm, kCMD_Write);

	return;
	}

	//********************************************************************
	void Measurement1D::SetDataValues(std::string dataFile) {
	//********************************************************************

	// Override this function if the input file isn't in a suitable format
	NUIS_LOG(SAM, "Reading data from: " << dataFile.c_str());
	fDataHist =
	PlotUtils::GetTH1DFromFile(dataFile, (fName + "_data"), fPlotTitles);
	fDataTrue = (TH1D *)fDataHist->Clone();

	// Number of data points is number of bins
	fNDataPointsX = fDataHist->GetXaxis()->GetNbins();

	return;
	};

	//********************************************************************
	void Measurement1D::SetDataFromDatabase(std::string inhistfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM, "Filling histogram from " << inhistfile << "->" << histname);
	fDataHist = PlotUtils::GetTH1DFromRootFile(
	(GeneralUtils::GetTopLevelDir() + "/data/" + inhistfile), histname);
	fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());

	return;
	};

	//********************************************************************
	void Measurement1D::SetDataFromFile(std::string inhistfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM, "Filling histogram from " << inhistfile << "->" << histname);
	fDataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname);
	fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());

	return;
	};

	//********************************************************************
	void Measurement1D::SetCovarMatrix(std::string covarFile) {
	//********************************************************************

	// Covariance function, only really used when reading in the MB Covariances.

	TFile *tempFile = new TFile(covarFile.c_str(), "READ");

	TH2D *covarPlot = new TH2D();
	TH2D *fFullCovarPlot = new TH2D();
	std::string covName = "";
	std::string covOption = FitPar::Config().GetParS("thrown_covariance");

	if (fIsShape \|\| fIsFree)
	covName = "shp_";
	if (fIsDiag)
	covName += "diag";
	else
	covName += "full";

	covarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str());

	if (!covOption.compare("SUB"))
	fFullCovarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str());
	else if (!covOption.compare("FULL"))
	fFullCovarPlot = (TH2D *)tempFile->Get("fullcov");
	else {
	NUIS_ERR(WRN, "Incorrect thrown_covariance option in parameters.");
	}

	int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral());
	int covdim = int(fDataHist->GetNbinsX());

	this->covar = new TMatrixDSym(dim);
	fFullCovar = new TMatrixDSym(dim);
	fDecomp = new TMatrixDSym(dim);

	int row, column = 0;
	row = 0;
	column = 0;
	for (Int_t i = 0; i < covdim; i++) {
	// if (this->masked->GetBinContent(i+1) > 0) continue;

	for (Int_t j = 0; j < covdim; j++) {
	// if (this->masked->GetBinContent(j+1) > 0) continue;

	(*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1);
	(*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1);

	column++;
	}
	column = 0;
	row++;
	}

	// Set bin errors on data
	if (!fIsDiag) {
	StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar);
	}

	// Get Deteriminant and inverse matrix
	// fCovDet = this->covar->Determinant();

	TDecompSVD LU = TDecompSVD(*this->covar);
	this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");

	return;
	};

	//********************************************************************
	// Sets the covariance matrix from a provided file in a text format
	// scale is a multiplicative pre-factor to apply in the case where the
	// covariance is given in some unit (e.g. 1E-38)
	void Measurement1D::SetCovarMatrixFromText(std::string covarFile, int dim,
	double scale) {
	//********************************************************************

	// Make a counter to track the line number
	int row = 0;

	std::string line;
	std::ifstream covarread(covarFile.c_str(), std::ifstream::in);

	this->covar = new TMatrixDSym(dim);
	fFullCovar = new TMatrixDSym(dim);
	if (covarread.is_open()) {
	NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile);
	} else {
	NUIS_ABORT("Covariance matrix provided is incorrect: " << covarFile);
	}

	// Loop over the lines in the file
	while (std::getline(covarread >> std::ws, line, '\n')) {
	int column = 0;

	// Loop over entries and insert them into matrix
	std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");

	if (entries.size() <= 1) {
	NUIS_ERR(WRN, "SetCovarMatrixFromText -> Covariance matrix only has <= 1 "
	- "entries on this line: "
	- << row);
	+ "entries on this line: "
	+ << row);
	}

	for (std::vector<double>::iterator iter = entries.begin();
	iter != entries.end(); iter++) {
	(covar)(row, column) = iter;
	(fFullCovar)(row, column) = iter;

	column++;
	}

	row++;
	}
	covarread.close();

	// Scale the actualy covariance matrix by some multiplicative factor
	(fFullCovar) = scale;

	// Robust matrix inversion method
	TDecompSVD LU = TDecompSVD(*this->covar);
	// THIS IS ACTUALLY THE INVERSE COVARIANCE MATRIXA AAAAARGH
	delete this->covar;
	this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");

	// Now need to multiply by the scaling factor
	// If the covariance
	(this->covar) = 1. / (scale);

	return;
	};

	//********************************************************************
	void Measurement1D::SetCovarMatrixFromCorrText(std::string corrFile, int dim) {
	//********************************************************************

	// Make a counter to track the line number
	int row = 0;

	std::string line;
	std::ifstream corr(corrFile.c_str(), std::ifstream::in);

	this->covar = new TMatrixDSym(dim);
	this->fFullCovar = new TMatrixDSym(dim);
	if (corr.is_open()) {
	- NUIS_LOG(SAM,
	- "Reading and converting correlation matrix from file: " << corrFile);
	+ NUIS_LOG(SAM, "Reading and converting correlation matrix from file: "
	+ << corrFile);
	} else {
	NUIS_ABORT("Correlation matrix provided is incorrect: " << corrFile);
	}

	while (std::getline(corr >> std::ws, line, '\n')) {
	int column = 0;

	// Loop over entries and insert them into matrix
	// Multiply by the errors to get the covariance, rather than the correlation
	// matrix
	std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
	for (std::vector<double>::iterator iter = entries.begin();
	iter != entries.end(); iter++) {
	double val = (iter) this->fDataHist->GetBinError(row + 1) * 1E38 *
	this->fDataHist->GetBinError(column + 1) * 1E38;
	if (val == 0) {
	NUIS_ABORT("Found a zero value in the covariance matrix, assuming "
	- "this is an error!");
	+ "this is an error!");
	}

	(*this->covar)(row, column) = val;
	(*this->fFullCovar)(row, column) = val;

	column++;
	}

	row++;
	}

	// Robust matrix inversion method
	TDecompSVD LU = TDecompSVD(*this->covar);
	delete this->covar;
	this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");

	return;
	};

	//********************************************************************
	// FullUnits refers to if we have "real" unscaled units in the covariance
	// matrix, e.g. 1E-76. If this is the case we need to scale it so that the chi2
	// contribution is correct NUISANCE internally assumes the covariance matrix has
	// units of 1E76
	void Measurement1D::SetCovarFromDataFile(std::string covarFile,
	std::string covName, bool FullUnits) {
	//********************************************************************

	NUIS_LOG(SAM, "Getting covariance from " << covarFile << "->" << covName);

	TFile *tempFile = new TFile(covarFile.c_str(), "READ");
	TH2D covPlot = (TH2D )tempFile->Get(covName.c_str());
	covPlot->SetDirectory(0);
	// Scale the covariance matrix if it comes in normal units
	if (FullUnits) {
	covPlot->Scale(1.E76);
	}

	int dim = covPlot->GetNbinsX();
	fFullCovar = new TMatrixDSym(dim);

	for (int i = 0; i < dim; i++) {
	for (int j = 0; j < dim; j++) {
	(*fFullCovar)(i, j) = covPlot->GetBinContent(i + 1, j + 1);
	}
	}

	this->covar = (TMatrixDSym *)fFullCovar->Clone();
	fDecomp = (TMatrixDSym *)fFullCovar->Clone();

	TDecompSVD LU = TDecompSVD(*this->covar);
	this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");

	TDecompChol LUChol = TDecompChol(*fDecomp);
	LUChol.Decompose();
	fDecomp = new TMatrixDSym(dim, LU.GetU().GetMatrixArray(), "");

	return;
	};

	// //********************************************************************
	// void Measurement1D::SetBinMask(std::string maskFile) {
	// //********************************************************************

	// // Create a mask histogram.
	// int nbins = fDataHist->GetNbinsX();
	// fMaskHist =
	// new TH1I((fName + "_fMaskHist").c_str(),
	// (fName + "_fMaskHist; Bin; Mask?").c_str(), nbins, 0, nbins);
	// std::string line;
	// std::ifstream mask(maskFile.c_str(), std::ifstream::in);

	// if (mask.is_open())
	// LOG(SAM) << "Reading bin mask from file: " << maskFile << std::endl;
	// else
	// LOG(FTL) << " Cannot find mask file." << std::endl;

	// while (std::getline(mask >> std::ws, line, '\n')) {
	// std::vector<int> entries = GeneralUtils::ParseToInt(line, " ");

	// // Skip lines with poorly formatted lines
	// if (entries.size() < 2) {
	// LOG(WRN) << "Measurement1D::SetBinMask(), couldn't parse line: " <<
	// line
	// << std::endl;
	// continue;
	// }

	// // The first index should be the bin number, the second should be the
	// mask
	// // value.
	// fMaskHist->SetBinContent(entries[0], entries[1]);
	// }

	// // Set masked data bins to zero
	// PlotUtils::MaskBins(fDataHist, fMaskHist);

	// return;
	// }

	// //********************************************************************
	// void Measurement1D::GetBinContents(std::vector<double>& cont,
	// std::vector<double>& err) {
	// //********************************************************************

	// // Return a vector of the main bin contents
	// for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
	// cont.push_back(fMCHist->GetBinContent(i + 1));
	// err.push_back(fMCHist->GetBinError(i + 1));
	// }

	// return;
	// };

	/*
	XSec Functions
	*/

	// //********************************************************************
	// void Measurement1D::SetFluxHistogram(std::string fluxFile, int minE, int
	// maxE,
	// double fluxNorm) {
	// //********************************************************************

	// // Note this expects the flux bins to be given in terms of MeV
	// LOG(SAM) << "Reading flux from file: " << fluxFile << std::endl;

	// TGraph f(fluxFile.c_str(), "%lg %lg");

	// fFluxHist =
	// new TH1D((fName + "_flux").c_str(), (fName + "; E_{#nu} (GeV)").c_str(),
	// f.GetN() - 1, minE, maxE);

	// Double_t* yVal = f.GetY();

	// for (int i = 0; i < fFluxHist->GetNbinsX(); ++i)
	// fFluxHist->SetBinContent(i + 1, yVal[i] * fluxNorm);
	// };

	// //********************************************************************
	// double Measurement1D::TotalIntegratedFlux(std::string intOpt, double low,
	// double high) {
	// //********************************************************************

	// if (fInput->GetType() == kGiBUU) {
	// return 1.0;
	// }

	// // The default case of low = -9999.9 and high = -9999.9
	// if (low == -9999.9) low = this->EnuMin;
	// if (high == -9999.9) high = this->EnuMax;

	// int minBin = fFluxHist->GetXaxis()->FindBin(low);
	// int maxBin = fFluxHist->GetXaxis()->FindBin(high);

	// // Get integral over custom range
	// double integral = fFluxHist->Integral(minBin, maxBin + 1, intOpt.c_str());

	// return integral;
	// };
	diff --git a/src/FitBase/Measurement1D.h b/src/FitBase/Measurement1D.h
	index dbd6471..6c6b0c1 100644
	--- a/src/FitBase/Measurement1D.h
	+++ b/src/FitBase/Measurement1D.h
	@@ -1,656 +1,659 @@
	// Copyright 2016 L. Pickering, P towell, 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 MEASUREMENT_1D_H_SEEN
	#define MEASUREMENT_1D_H_SEEN

	/*!
	* \addtogroup FitBase
	* @{
	*/

	#include <math.h>
	#include <stdlib.h>
	#include <deque>
	#include <iomanip>
	#include <iostream>
	#include <numeric>
	#include <sstream>
	#include <string>

	// ROOT includes
	#include <TArrayF.h>
	#include <TCanvas.h>
	#include <TCut.h>
	#include <TDecompChol.h>
	#include <TDecompSVD.h>
	#include <TGraph.h>
	#include <TGraphErrors.h>
	#include <TH1D.h>
	#include <TH2D.h>
	#include <TMatrixDSym.h>
	#include <TROOT.h>
	#include <TSystem.h>

	// External data fit includes
	#include "FitEvent.h"

	#include "FitUtils.h"
	#include "MeasurementBase.h"
	#include "PlotUtils.h"
	#include "StatUtils.h"

	#include "SignalDef.h"
	#include "MeasurementVariableBox.h"
	#include "MeasurementVariableBox1D.h"

	namespace NUISANCE {
	namespace FitBase {

	}
	}

	//********************************************************************
	/// 1D Measurement base class. Histogram handling is done in this base layer.
	class Measurement1D : public MeasurementBase {
	//********************************************************************

	public:
	/*
	Constructor/Deconstuctor
	*/
	Measurement1D(void);
	virtual ~Measurement1D(void);

	/*
	Setup Functions
	*/

	/// \brief Setup all configs once initialised
	///
	/// Should be called after all configs have been setup inside fSettings container.
	/// Handles the processing of inputs and setting up of types.
	/// Replaces the old 'SetupMeasurement' function.
	void FinaliseSampleSettings();

	/// \brief Creates the 1D data distribution given the binning provided.
	virtual void CreateDataHistogram(int dimx, double* binx);

	/// \brief Read 1D data inputs from a text file.
	///
	/// Inputfile should have the format: \n
	/// low_binedge_1 bin_content_1 bin_error_1 \n
	/// low_binedge_2 bin_content_2 bin_error_2 \n
	/// .... .... .... \n
	/// high_bin_edge_N 0.0 0.0
	virtual void SetDataFromTextFile(std::string datafile);

	/// \brief Read 1D data inputs from a root file.
	///
	/// inhistfile specifies the path to the root file
	/// histname specifies the name of the histogram.
	///
	/// If no histogram name is given the inhistfile value
	/// is automatically parsed with ';' so that: \n
	/// 'myhistfile.root;myhistname' \n
	/// will also work.
	virtual void SetDataFromRootFile(std::string inhistfile, std::string histname = "");


	/// \brief Setup a default empty data histogram
	///
	/// Only used for flattree creators.
	virtual void SetEmptyData();

	/// \brief Set data bin errors to sqrt(entries)
	///
	/// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST
	///
	/// Sets the data errors as the sqrt of the bin contents
	/// Should be use for counting experiments
	virtual void SetPoissonErrors();

	/// \brief Make diagonal covariance from data
	///
	/// \warning If no histogram passed, data must be setup first!
	/// Setup the covariance inputs by taking the data histogram
	/// errors and setting up a diagonal covariance matrix.
	///
	/// If no data is supplied, fDataHist is used if already set.
	virtual void SetCovarFromDiagonal(TH1D* data = NULL);

	/// \brief Read the data covariance from a text file.
	///
	/// Inputfile should have the format: \n
	/// covariance_11 covariance_12 covariance_13 ... \n
	/// covariance_21 covariance_22 covariance_23 ... \n
	/// ... ... ... ... \n
	///
	/// If no dimensions are given, it is assumed from the number
	/// entries in the first line of covfile.
	virtual void SetCovarFromTextFile(std::string covfile, int dim = -1);

	virtual void SetCovarFromMultipleTextFiles(std::string covfiles, int dim = -1);

	/// \brief Read the data covariance from a ROOT file.
	///
	/// - covfile specifies the full path to the file
	/// - histname specifies the name of the covariance object. Both TMatrixDSym and TH2D are supported.
	///
	/// If no histogram name is given the inhistfile value
	/// is automatically parsed with ; so that: \n
	/// mycovfile.root;myhistname \n
	/// will also work.
	virtual void SetCovarFromRootFile(std::string covfile, std::string histname="");

	/// \brief Read the inverted data covariance from a text file.
	///
	/// Inputfile should have similar format to that shown
	/// in SetCovarFromTextFile.
	///
	/// If no dimensions are given, it is assumed from the number
	/// entries in the first line of covfile.
	virtual void SetCovarInvertFromTextFile(std::string covfile, int dim = -1);


	/// \brief Read the inverted data covariance from a ROOT file.
	///
	/// Inputfile should have similar format to that shown
	/// in SetCovarFromRootFile.
	///
	/// If no histogram name is given the inhistfile value
	/// is automatically parsed with ; so that: \n
	/// mycovfile.root;myhistname \n
	/// will also work.
	virtual void SetCovarInvertFromRootFile(std::string covfile, std::string histname="");

	/// \brief Read the data correlations from a text file.
	///
	/// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST
	///
	/// Inputfile should have similar format to that shown
	/// in SetCovarFromTextFile.
	///
	/// If no dimensions are given, it is assumed from the number
	/// entries in the first line of covfile.
	virtual void SetCorrelationFromTextFile(std::string covfile, int dim = -1);

	/// \brief Read the data correlations from multiple text files.
	///
	/// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST
	///
	/// Inputfile should have similar format to that shown
	/// in SetCovarFromTextFile.
	///
	/// If no dimensions are given, it is assumed from the number
	/// entries in the first line of the first corrfile.
	virtual void SetCorrelationFromMultipleTextFiles(std::string corrfiles, int dim = -1);

	/// \brief Read the data correlations from a ROOT file.
	///
	/// \warning REQUIRES DATA TO BE SET FIRST
	///
	/// Inputfile should have similar format to that shown
	/// in SetCovarFromRootFile.
	///
	/// If no histogram name is given the inhistfile value
	/// is automatically parsed with ; so that: \n
	/// mycovfile.root;myhistname \n
	/// will also work.
	virtual void SetCorrelationFromRootFile(std::string covfile, std::string histname="");


	/// \brief Read the cholesky decomposed covariance from a text file and turn it into a covariance
	///
	/// Inputfile should have similar format to that shown
	/// in SetCovarFromTextFile.
	///
	/// If no dimensions are given, it is assumed from the number
	/// entries in the first line of covfile.
	virtual void SetCholDecompFromTextFile(std::string covfile, int dim = -1);


	/// \brief Read the cholesky decomposed covariance from a ROOT file and turn it into a covariance
	///
	/// Inputfile should have similar format to that shown
	/// in SetCovarFromRootFile.
	///
	/// If no histogram name is given the inhistfile value
	/// is automatically parsed with ; so that: \n
	/// mycovfile.root;myhistname \n
	/// will also work.
	virtual void SetCholDecompFromRootFile(std::string covfile, std::string histname="");

	/// \brief Try to extract a shape-only matrix from the existing covariance
	virtual void SetShapeCovar();

	/// \brief Scale the data by some scale factor
	virtual void ScaleData(double scale);


	/// \brief Scale the data error bars by some scale factor
	virtual void ScaleDataErrors(double scale);


	/// \brief Scale the covariaince and its invert/decomp by some scale factor.
	virtual void ScaleCovar(double scale);



	/// \brief Setup a bin masking histogram and apply masking to data
	///
	/// \warning REQUIRES DATA HISTOGRAM TO BE SET FIRST
	///
	/// Reads in a list of bins in a text file to be masked. Format is: \n
	/// bin_index_1 1 \n
	/// bin_index_2 1 \n
	/// bin_index_3 1 \n
	///
	/// If 0 is given then a bin entry will NOT be masked. So for example: \n\n
	/// 1 1 \n
	/// 2 1 \n
	/// 3 0 \n
	/// 4 1 \n\n
	/// Will mask only the 1st, 2nd, and 4th bins.
	///
	/// Masking can be turned on by specifiying the MASK option when creating a sample.
	/// When this is passed NUISANCE will look in the following locations for the mask file:
	/// - FitPar::Config().GetParS(fName + ".mask")
	/// - "data/masks/" + fName + ".mask";
	virtual void SetBinMask(std::string maskfile);


	/// \brief Set the current fit options from a string.
	///
	/// This is called twice for each sample, once to set the default
	/// and once to set the current setting (if anything other than default given)
	///
	/// For this to work properly it requires the default and allowed types to be
	/// set correctly. These should be specified as a string listing options.
	///
	/// To split up options so that NUISANCE can automatically detect ones that
	/// are conflicting. Any options seperated with the '/' symbol are non conflicting
	/// and can be given together, whereas any seperated with the ',' symbol cannot
	/// be specified by the end user at the same time.
	///
	/// Default Type Examples:
	/// - DIAG/FIX = Default option will be a diagonal covariance, with FIXED norm.
	/// - MASK/SHAPE = Default option will be a masked hist, with SHAPE always on.
	///
	/// Allowed Type examples:
	/// - 'FULL/DIAG/NORM/MASK' = Any of these options can be specified.
	/// - 'FULL,FREE,SHAPE/MASK/NORM' = User can give either FULL, FREE, or SHAPE as on option.
	/// MASK and NORM can also be included as options.
	virtual void SetFitOptions(std::string opt);


	/// \brief Final constructor setup
	/// \warning Should be called right at the end of the constructor.
	///
	/// Contains a series of checks to ensure the data and inputs have been setup.
	/// Also creates the MC histograms needed for fitting.
	void FinaliseMeasurement();




	/*
	Smearing
	*/
	/// \brief Read in smearing matrix from file
	///
	/// Set the smearing matrix from a text file given the size of the matrix
	virtual void SetSmearingMatrix(std::string smearfile, int truedim,
	int recodim);

	/// \brief Apply smearing to MC true to get MC reco
	///
	/// Apply smearing matrix to fMCHist using fSmearingMatrix
	virtual void ApplySmearingMatrix(void);



	/*
	Reconfigure Functions
	*/

	/// \brief Create a Measurement1D box
	///
	/// Creates a new 1D variable box containing just fXVar.
	///
	/// This box is the bare minimum required by the JointFCN when
	/// running fast reconfigures during a routine.
	/// If for some reason a sample needs extra variables to be saved then
	/// it should override this function creating its own MeasurementVariableBox
	/// that contains the extra variables.
	virtual MeasurementVariableBox* CreateBox() {return new MeasurementVariableBox1D();};

	/// \brief Reset all MC histograms
	///
	/// Resets all standard histograms and those registered to auto
	/// process to zero.
	///
	/// If extra histograms are not included in auto processing, then they must be reset
	/// by overriding this function and doing it manually if required.
	virtual void ResetAll(void);

	/// \brief Fill MC Histograms from XVar
	///
	/// Fill standard histograms using fXVar, Weight read from the variable box.
	///
	/// WARNING : Any extra MC histograms need to be filled by overriding this function,
	/// even if they have been set to auto process.
	virtual void FillHistograms(void);

	// \brief Convert event rates to final histogram
	///
	/// Apply standard scaling procedure to standard mc histograms to convert from
	/// raw events to xsec prediction.
	///
	/// If any distributions have been set to auto process
	/// that is done during this function call, and a differential xsec is assumed.
	/// If that is not the case this function must be overriden.
	virtual void ScaleEvents(void);

	/// \brief Scale MC by a factor=1/norm
	///
	/// Apply a simple normalisation scaling if the option FREE or a norm_parameter
	/// has been specified in the NUISANCE routine.
	virtual void ApplyNormScale(double norm);


	/*
	Statistical Functions
	*/

	/// \brief Get Number of degrees of freedom
	///
	/// Returns the number bins inside the data histogram accounting for
	/// any bin masking applied.
	virtual int GetNDOF(void);

	/// \brief Return Data/MC Likelihood at current state
	///
	/// Returns the likelihood of the data given the current MC prediction.
	/// Diferent likelihoods definitions are used depending on the FitOptions.
	virtual double GetLikelihood(void);


	/*
	Fake Data
	*/

	/// \brief Set the fake data values from either a file, or MC
	///
	/// - Setting from a file "path": \n
	/// When reading from a file the full path must be given to a standard
	/// nuisance output. The standard MC histogram should have a name that matches
	/// this sample for it to be read in.
	/// \n\n
	/// - Setting from "MC": \n
	/// If the MC option is given the current MC prediction is used as fake data.
	virtual void SetFakeDataValues(std::string fakeOption);

	/// \brief Reset fake data back to starting fake data
	///
	/// Reset the fake data back to original fake data (Reset back to before
	/// ThrowCovariance was first called)
	virtual void ResetFakeData(void);

	/// \brief Reset fake data back to original data
	///
	/// Reset the data histogram back to the true original dataset for this sample
	/// before any fake data was defined.
	virtual void ResetData(void);

	/// \brief Generate fake data by throwing the covariance.
	///
	/// Can be used on fake MC data or just the original dataset.
	/// Call ResetFakeData or ResetData to return to values before the throw.
	virtual void ThrowCovariance(void);

	/// \brief Throw the data by its assigned errors and assign this to MC
	///
	/// Used when creating data toys by assign the MC to this thrown data
	/// so that the likelihood is calculated between data and thrown data
	virtual void ThrowDataToy(void);


	/*
	Access Functions
	*/

	/// \brief Returns nicely formatted MC Histogram
	///
	/// Format options can also be given in the samplesettings:
	/// - linecolor
	/// - linestyle
	/// - linewidth
	/// - fillcolor
	/// - fillstyle
	///
	/// So to have a sample line colored differently in the xml cardfile put: \n
	/// <sample name="MiniBooNE_CCQE_XSec_1DQ2_nu" input="NEUT:input.root"
	/// linecolor="2" linestyle="7" linewidth="2" />
	virtual TH1D* GetMCHistogram(void);

	/// \brief Returns nicely formatted data Histogram
	///
	/// Format options can also be given in the samplesettings:
	/// - datacolor
	/// - datastyle
	/// - datawidth
	///
	/// So to have a sample data colored differently in the xml cardfile put: \n
	/// <sample name="MiniBooNE_CCQE_XSec_1DQ2_nu" input="NEUT:input.root"
	/// datacolor="2" datastyle="7" datawidth="2" />
	virtual TH1D* GetDataHistogram(void);

	/// \brief Returns a list of all MC histograms.
	///
	/// Override this if you have extra histograms that need to be
	/// accessed outside of the Measurement1D class.
	inline virtual std::vector<TH1*> GetMCList(void) {
	return std::vector<TH1*>(1, GetMCHistogram());
	}

	/// \brief Returns a list of all Data histograms.
	///
	/// Override this if you have extra histograms that need to be
	/// accessed outside of the Measurement1D class.
	inline virtual std::vector<TH1*> GetDataList(void) {
	return std::vector<TH1*>(1, GetDataHistogram());
	}

	/// \brief Returns a list of all Mask histograms.
	///
	/// Override this if you have extra histograms that need to be
	/// accessed outside of the Measurement1D class.
	inline virtual std::vector<TH1*> GetMaskList(void) {
	return std::vector<TH1*>(1, fMaskHist);
	};

	/// \brief Returns a list of all Fine histograms.
	///
	/// Override this if you have extra histograms that need to be
	/// accessed outside of the Measurement1D class.
	inline virtual std::vector<TH1*> GetFineList(void) {
	return std::vector<TH1*>(1, fMCFine);
	};


	/*
	Write Functions
	*/
	/// \brief Save the current state to the current TFile directory \n
	///
	/// Data/MC are both saved by default.
	/// A range of other histograms can be saved by setting the
	/// config option 'drawopts'.
	///
	/// Possible options: \n
	/// - FINE = Write Fine Histogram \n
	/// - WEIGHTS = Write Weighted MC Histogram (before scaling) \n
	/// - FLUX = Write Flux Histogram from MC Input \n
	/// - EVT = Write Event Histogram from MC Input \n
	/// - XSEC = Write XSec Histogram from MC Input \n
	/// - MASK = Write Mask Histogram \n
	/// - COV = Write Covariance Histogram \n
	/// - INVCOV = Write Inverted Covariance Histogram \n
	/// - DECMOP = Write Decomp. Covariance Histogram \n
	/// - RESIDUAL= Write Resudial Histograms \n
	/// - RATIO = Write Data/MC Ratio Histograms \n
	/// - SHAPE = Write MC Shape Histograms norm. to Data \n
	/// - CANVMC = Build MC Canvas Showing Data, MC, Shape \n
	/// - MODES = Write PDG Stack \n
	/// - CANVPDG = Build MC Canvas Showing Data, PDGStack \n
	///
	/// So to save a range of these in parameters/config.xml set: \n
	/// <config drawopts='FINE/COV/SHAPE/RATIO' />
	virtual void Write(std::string drawOpt);



	virtual void WriteRatioPlot();




	virtual void WriteShapePlot();
	virtual void WriteShapeRatioPlot();




	//*
	// OLD DEFUNCTIONS
	//

	/// OLD FUNCTION
	virtual void SetupMeasurement(std::string input, std::string type,
	FitWeight* rw, std::string fkdt);

	/// OLD FUNCTION
	virtual void SetupDefaultHist(void);

	/// OLD FUNCTION
	virtual void SetDataValues(std::string dataFile);

	/// OLD FUNCTION
	virtual void SetDataFromFile(std::string inhistfile, std::string histname);
	/// OLD FUNCTION
	virtual void SetDataFromDatabase(std::string inhistfile,
	std::string histname);

	/// OLD FUNCTION
	virtual void SetCovarMatrix(std::string covarFile);
	/// OLD FUNCTION
	virtual void SetCovarMatrixFromText(std::string covarFile, int dim,
	double scale = 1.0);
	/// OLD FUNCTION
	virtual void SetCovarMatrixFromCorrText(std::string covarFile, int dim);


	/// OLD FUNCTION
	virtual void SetCovarFromDataFile(std::string covarFile, std::string covName,
	bool FullUnits = false);


	/// OLD FUNCTION
	// virtual THStack GetModeStack(void);

	protected:

	// Data
	TH1D* fDataHist; ///< default data histogram
	TH1D* fDataOrig; ///< histogram to store original data before throws.
	TH1D* fDataTrue; ///< histogram to store true dataset
	std::string fPlotTitles; ///< Plot title x and y for the histograms


	// MC
	TH1D* fMCHist; ///< default MC Histogram used in the chi2 fits
	TH1D* fMCFine; ///< finely binned MC histogram
	TH1D* fMCStat; ///< histogram with unweighted events to properly calculate
	TH1D* fMCWeighted; ///< Weighted histogram before xsec scaling

	TH1I* fMaskHist; ///< Mask histogram for neglecting specific bins
	TMatrixD* fSmearMatrix; ///< Smearing matrix (note, this is not symmetric)

	+ TH1D *fResidualHist;
	+ TH1D *fChi2LessBinHist;
	+
	TrueModeStack* fMCHist_Modes; ///< Optional True Mode Stack


	// Statistical
	TMatrixDSym* covar; ///< Inverted Covariance
	TMatrixDSym* fFullCovar; ///< Full Covariance
	TMatrixDSym* fDecomp; ///< Decomposed Covariance
	TMatrixDSym* fCorrel; ///< Correlation Matrix

	TMatrixDSym* fShapeCovar; ///< Shape-only covariance
	TMatrixDSym* fShapeDecomp; ///< Decomposed shape-only covariance
	TMatrixDSym* fShapeInvert; ///< Inverted shape-only covariance

	TMatrixDSym* fCovar; ///< New FullCovar
	TMatrixDSym* fInvert; ///< New covar

	double fNormError; ///< Sample norm error

	double fLikelihood; ///< Likelihood value

	// Fake Data
	bool fIsFakeData; ///< Flag: is the current data fake from MC
	std::string fFakeDataInput; ///< Input fake data file path
	TFile* fFakeDataFile; ///< Input fake data file


	// Fit specific flags
	std::string fFitType; ///< Current fit type
	std::string fAllowedTypes; ///< Fit Types Possible
	std::string fDefaultTypes; ///< Starting Default Fit Types

	bool fIsShape; ///< Flag : Perform Shape-only fit
	bool fIsFree; ///< Flag : Perform normalisation free fit
	bool fIsDiag; ///< Flag : only include uncorrelated diagonal errors
	bool fIsMask; ///< Flag : Apply bin masking
	bool fIsRawEvents; ///< Flag : Are bin contents just event rates
	bool fIsEnu1D; ///< Flag : Perform Flux Unfolded Scaling
	bool fIsChi2SVD; ///< Flag : Use alternative Chi2 SVD Method (Do not use)
	bool fAddNormPen; ///< Flag : Add a normalisation penalty term to the chi2.
	bool fIsFix; ///< Flag : keeping norm fixed
	bool fIsFull; ///< Flag : using full covariaince
	bool fIsDifXSec; ///< Flag : creating a dif xsec
	bool fIsChi2; ///< Flag : using Chi2 over LL methods
	bool fIsSmeared; ///< Flag : Apply smearing?




	/// OLD STUFF TO REMOVE
	TH1D* fMCHist_PDG[61]; ///< REMOVE OLD MC PDG Plot

	// Arrays for data entries
	Double_t* fXBins; ///< REMOVE xBin edges
	Double_t* fDataValues; ///< REMOVE data bin contents
	Double_t* fDataErrors; ///< REMOVE data bin errors
	Int_t fNDataPointsX; ///< REMOVE number of data points

	};

	/! @} /
	#endif
	diff --git a/src/FitBase/Measurement2D.cxx b/src/FitBase/Measurement2D.cxx
	index 0cc89f3..fd2cb59 100644
	--- a/src/FitBase/Measurement2D.cxx
	+++ b/src/FitBase/Measurement2D.cxx
	@@ -1,2018 +1,2020 @@
	// 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 "Measurement2D.h"
	#include "TDecompChol.h"

	//********************************************************************
	Measurement2D::Measurement2D(void) {
	//********************************************************************

	covar = NULL;
	fDecomp = NULL;
	fFullCovar = NULL;

	fMCHist = NULL;
	fMCFine = NULL;
	fDataHist = NULL;

	fMCHist_X = NULL;
	fMCHist_Y = NULL;
	fDataHist_X = NULL;
	fDataHist_Y = NULL;

	fMaskHist = NULL;
	fMapHist = NULL;
	fDataOrig = NULL;
	fDataTrue = NULL;
	fMCWeighted = NULL;

	fResidualHist = NULL;
	+ fChi2LessBinHist = NULL;

	fDefaultTypes = "FIX/FULL/CHI2";
	fAllowedTypes =
	"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/FITPROJX/"
	"FITPROJY";

	fIsFix = false;
	fIsShape = false;
	fIsFree = false;

	fIsDiag = false;
	fIsFull = false;

	fAddNormPen = false;
	fIsMask = false;
	fIsChi2SVD = false;

	fIsRawEvents = false;
	fIsDifXSec = false;
	fIsEnu = false;

	// XSec Scalings
	fScaleFactor = -1.0;
	fCurrentNorm = 1.0;

	// Histograms
	fDataHist = NULL;
	fDataTrue = NULL;

	fMCHist = NULL;
	fMCFine = NULL;
	fMCWeighted = NULL;

	fMaskHist = NULL;

	// Covar
	covar = NULL;
	fFullCovar = NULL;

	fCovar = NULL;
	fInvert = NULL;
	fDecomp = NULL;

	// Fake Data
	fFakeDataInput = "";
	fFakeDataFile = NULL;

	// Options
	fDefaultTypes = "FIX/FULL/CHI2";
	fAllowedTypes =
	"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK";

	fIsFix = false;
	fIsShape = false;
	fIsFree = false;
	fIsDiag = false;
	fIsFull = false;
	fAddNormPen = false;
	fIsMask = false;
	fIsChi2SVD = false;
	fIsRawEvents = false;
	fIsDifXSec = false;
	fIsEnu1D = false;

	// Inputs
	fInput = NULL;
	fRW = NULL;

	// Extra Histograms
	fMCHist_Modes = NULL;
	}

	//********************************************************************
	Measurement2D::~Measurement2D(void) {
	//********************************************************************

	if (fDataHist)
	delete fDataHist;
	if (fDataTrue)
	delete fDataTrue;
	if (fMCHist)
	delete fMCHist;
	if (fMCFine)
	delete fMCFine;
	if (fMCWeighted)
	delete fMCWeighted;
	if (fMaskHist)
	delete fMaskHist;
	if (covar)
	delete covar;
	if (fFullCovar)
	delete fFullCovar;
	if (fCovar)
	delete fCovar;
	if (fInvert)
	delete fInvert;
	if (fDecomp)
	delete fDecomp;

	delete fResidualHist;
	+ delete fChi2LessBinHist;
	}

	//********************************************************************
	void Measurement2D::FinaliseSampleSettings() {
	//********************************************************************

	MeasurementBase::FinaliseSampleSettings();

	// Setup naming + renaming
	fName = fSettings.GetName();
	fSettings.SetS("originalname", fName);
	if (fSettings.Has("rename")) {
	fName = fSettings.GetS("rename");
	fSettings.SetS("name", fName);
	}

	// Setup all other options
	NUIS_LOG(SAM, "Finalising Sample Settings: " << fName);

	if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) {
	fIsRawEvents = true;
	NUIS_LOG(SAM, "Found event rate measurement but using poisson likelihoods.");
	}

	if (fSettings.GetS("originalname").find("Enu") != std::string::npos) {
	fIsEnu1D = true;
	NUIS_LOG(SAM, "::" << fName << "::");
	NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, "
	<< "not flux averaged!");
	}

	if (fIsEnu1D && fIsRawEvents) {
	NUIS_ERR(FTL, "Found 2D Enu XSec distribution AND fIsRawEvents, is this "
	"really correct?!");
	NUIS_ERR(FTL, "Check experiment constructor for " << fName
	<< " and correct this!");
	NUIS_ABORT("I live in " << __FILE__ << ":" << __LINE__);
	}

	if (!fRW)
	fRW = FitBase::GetRW();
	if (!fInput)
	SetupInputs(fSettings.GetS("input"));

	// Setup options
	SetFitOptions(fDefaultTypes); // defaults
	SetFitOptions(fSettings.GetS("type")); // user specified

	EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min"));
	EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max"));

	if (fAddNormPen) {
	fNormError = fSettings.GetNormError();
	if (fNormError <= 0.0) {
	NUIS_ERR(FTL, "Norm error for class " << fName << " is 0.0!");
	NUIS_ABORT("If you want to use it please add fNormError=VAL");
	}
	}
	}

	void Measurement2D::CreateDataHistogram(int dimx, double *binx, int dimy,
	double *biny) {
	if (fDataHist)
	delete fDataHist;

	NUIS_LOG(SAM, "Creating Data Histogram dim : " << dimx << " " << dimy);

	fDataHist = new TH2D((fSettings.GetName() + "_data").c_str(),
	(fSettings.GetFullTitles()).c_str(), dimx - 1, binx,
	dimy - 1, biny);
	}

	void Measurement2D::SetDataFromTextFile(std::string data, std::string binx,
	std::string biny) {
	// Get the data hist
	fDataHist = PlotUtils::GetTH2DFromTextFile(data, binx, biny);
	// Set the name properly
	fDataHist->SetName((fSettings.GetName() + "_data").c_str());
	fDataHist->SetTitle(fSettings.GetFullTitles().c_str());
	}

	void Measurement2D::SetDataFromRootFile(std::string datfile,
	std::string histname) {
	NUIS_LOG(SAM, "Reading data from root file: " << datfile << ";" << histname);
	fDataHist = PlotUtils::GetTH2DFromRootFile(datfile, histname);
	fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
	(fSettings.GetFullTitles()).c_str());
	}

	void Measurement2D::SetDataValuesFromTextFile(std::string datfile, TH2D *hist) {

	NUIS_LOG(SAM, "Setting data values from text file");
	if (!hist)
	hist = fDataHist;

	// Read TH2D From textfile
	TH2D valhist = (TH2D )hist->Clone();
	valhist->Reset();
	PlotUtils::Set2DHistFromText(datfile, valhist, 1.0, true);

	NUIS_LOG(SAM, " -> Filling values from read hist.");
	for (int i = 0; i < valhist->GetNbinsX(); i++) {
	for (int j = 0; j < valhist->GetNbinsY(); j++) {
	hist->SetBinContent(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1));
	}
	}
	NUIS_LOG(SAM, " --> Done");
	}

	void Measurement2D::SetDataErrorsFromTextFile(std::string datfile, TH2D *hist) {
	NUIS_LOG(SAM, "Setting data errors from text file");

	if (!hist)
	hist = fDataHist;

	// Read TH2D From textfile
	TH2D valhist = (TH2D )hist->Clone();
	valhist->Reset();
	PlotUtils::Set2DHistFromText(datfile, valhist, 1.0);

	// Fill Errors
	NUIS_LOG(SAM, " -> Filling errors from read hist.");

	for (int i = 0; i < valhist->GetNbinsX(); i++) {
	for (int j = 0; j < valhist->GetNbinsY(); j++) {
	hist->SetBinError(i + 1, j + 1, valhist->GetBinContent(i + 1, j + 1));
	}
	}
	NUIS_LOG(SAM, " --> Done");
	}

	void Measurement2D::SetMapValuesFromText(std::string dataFile) {

	TH2D *hist = fDataHist;
	std::vector<double> edgex;
	std::vector<double> edgey;

	for (int i = 0; i <= hist->GetNbinsX(); i++)
	edgex.push_back(hist->GetXaxis()->GetBinLowEdge(i + 1));
	for (int i = 0; i <= hist->GetNbinsY(); i++)
	edgey.push_back(hist->GetYaxis()->GetBinLowEdge(i + 1));

	fMapHist = new TH2I((fName + "_map").c_str(), (fName + fPlotTitles).c_str(),
	edgex.size() - 1, &edgex[0], edgey.size() - 1, &edgey[0]);

	NUIS_LOG(SAM, "Reading map from: " << dataFile);
	PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0);
	}

	//********************************************************************
	void Measurement2D::SetPoissonErrors() {
	//********************************************************************

	if (!fDataHist) {
	NUIS_ERR(FTL, "Need a data hist to setup possion errors! ");
	NUIS_ABORT("Setup Data First!");
	}

	for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
	fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1)));
	}
	}

	//********************************************************************
	void Measurement2D::SetCovarFromDiagonal(TH2D *data) {
	//********************************************************************

	if (!data and fDataHist) {
	data = fDataHist;
	}

	if (data) {
	NUIS_LOG(SAM, "Setting diagonal covariance for: " << data->GetName());
	fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	} else {
	NUIS_ABORT("No data input provided to set diagonal covar from!");
	}

	// if (!fIsDiag) {
	// ERR(FTL) << "SetCovarMatrixFromDiag called for measurement "
	// << "that is not set as diagonal." << std::endl;
	// throw;
	// }
	}

	//********************************************************************
	void Measurement2D::SetCovarFromTextFile(std::string covfile, int dim) {
	//********************************************************************

	if (dim == -1) {
	dim = this->GetNDOF();
	}

	NUIS_LOG(SAM, "Reading covariance from text file: " << covfile << " " << dim);
	fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement2D::SetCovarFromRootFile(std::string covfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM,
	"Reading covariance from text file: " << covfile << ";" << histname);
	fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname);
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement2D::SetCovarInvertFromTextFile(std::string covfile, int dim) {
	//********************************************************************

	if (dim == -1) {
	dim = this->GetNDOF();
	}

	NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile);
	covar = StatUtils::GetCovarFromTextFile(covfile, dim);
	fFullCovar = StatUtils::GetInvert(covar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement2D::SetCovarInvertFromRootFile(std::string covfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile << ";"
	<< histname);
	covar = StatUtils::GetCovarFromRootFile(covfile, histname);
	fFullCovar = StatUtils::GetInvert(covar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	//********************************************************************
	void Measurement2D::SetCorrelationFromTextFile(std::string covfile, int dim) {
	//********************************************************************

	if (dim == -1)
	dim = this->GetNDOF();
	NUIS_LOG(SAM,
	"Reading data correlations from text file: " << covfile << ";" << dim);
	TMatrixDSym *correlation = StatUtils::GetCovarFromTextFile(covfile, dim);

	if (!fDataHist) {
	NUIS_ABORT("Trying to set correlations from text file but there is no "
	"data to build it from. \n"
	<< "In constructor make sure data is set before "
	"SetCorrelationFromTextFile is called. \n");
	}

	// Fill covar from data errors and correlations
	fFullCovar = new TMatrixDSym(dim);
	for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
	for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
	(fFullCovar)(i, j) = (correlation)(i, j) *
	fDataHist->GetBinError(i + 1) *
	fDataHist->GetBinError(j + 1) * 1.E76;
	}
	}

	// Fill other covars.
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);

	delete correlation;
	}

	//********************************************************************
	void Measurement2D::SetCorrelationFromRootFile(std::string covfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";"
	<< histname);
	TMatrixDSym *correlation = StatUtils::GetCovarFromRootFile(covfile, histname);

	if (!fDataHist) {
	NUIS_ABORT("Trying to set correlations from text file but there is no "
	"data to build it from. \n"
	<< "In constructor make sure data is set before "
	"SetCorrelationFromTextFile is called. \n");
	}

	// Fill covar from data errors and correlations
	fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
	for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
	for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
	(fFullCovar)(i, j) = (correlation)(i, j) *
	fDataHist->GetBinError(i + 1) *
	fDataHist->GetBinError(j + 1) * 1.E76;
	}
	}

	// Fill other covars.
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);

	delete correlation;
	}

	//********************************************************************
	void Measurement2D::SetCholDecompFromTextFile(std::string covfile, int dim) {
	//********************************************************************

	if (dim == -1) {
	dim = this->GetNDOF();
	}

	NUIS_LOG(SAM, "Reading cholesky from text file: " << covfile << " " << dim);
	TMatrixD *temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim);

	TMatrixD trans = (TMatrixD )temp->Clone();
	trans->T();
	(trans) = (*temp);

	fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);

	delete temp;
	delete trans;
	}

	//********************************************************************
	void Measurement2D::SetCholDecompFromRootFile(std::string covfile,
	std::string histname) {
	//********************************************************************

	NUIS_LOG(SAM, "Reading cholesky decomp from root file: " << covfile << ";"
	<< histname);
	TMatrixD *temp = StatUtils::GetMatrixFromRootFile(covfile, histname);

	TMatrixD trans = (TMatrixD )temp->Clone();
	trans->T();
	(trans) = (*temp);

	fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), "");
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);

	delete temp;
	delete trans;
	}

	//********************************************************************
	void Measurement2D::ScaleData(double scale) {
	//********************************************************************
	fDataHist->Scale(scale);
	}

	//********************************************************************
	void Measurement2D::ScaleDataErrors(double scale) {
	//********************************************************************
	for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
	for (int j = 0; j < fDataHist->GetNbinsY(); j++) {
	fDataHist->SetBinError(i + 1, j + 1,
	fDataHist->GetBinError(i + 1, j + 1) * scale);
	}
	}
	}

	//********************************************************************
	void Measurement2D::ScaleCovar(double scale) {
	//********************************************************************
	(fFullCovar) = scale;
	(covar) = 1.0 / scale;
	(fDecomp) = sqrt(scale);
	}

	//********************************************************************
	void Measurement2D::SetBinMask(std::string maskfile) {
	//********************************************************************

	if (!fIsMask)
	return;
	NUIS_LOG(SAM, "Reading bin mask from file: " << maskfile);

	// Create a mask histogram with dim of data
	int nbinsx = fDataHist->GetNbinsX();
	int nbinxy = fDataHist->GetNbinsY();
	fMaskHist = new TH2I((fSettings.GetName() + "_BINMASK").c_str(),
	(fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(),
	nbinsx, 0, nbinsx, nbinxy, 0, nbinxy);
	std::string line;
	std::ifstream mask(maskfile.c_str(), std::ifstream::in);

	if (!mask.is_open()) {
	NUIS_LOG(FTL, " Cannot find mask file.");
	throw;
	}

	while (std::getline(mask >> std::ws, line, '\n')) {
	std::vector<int> entries = GeneralUtils::ParseToInt(line, " ");

	// Skip lines with poorly formatted lines
	if (entries.size() < 2) {
	NUIS_LOG(WRN, "Measurement2D::SetBinMask(), couldn't parse line: " << line);
	continue;
	}

	// The first index should be the bin number, the second should be the mask
	// value.
	int val = 0;
	if (entries[2] > 0)
	val = 1;
	fMaskHist->SetBinContent(entries[0], entries[1], val);
	}

	// Apply masking by setting masked data bins to zero
	PlotUtils::MaskBins(fDataHist, fMaskHist);

	return;
	}

	//********************************************************************
	void Measurement2D::FinaliseMeasurement() {
	//********************************************************************

	NUIS_LOG(SAM, "Finalising Measurement: " << fName);
	if (fSettings.GetB("onlymc")) {
	if (fDataHist)
	delete fDataHist;
	fDataHist = new TH2D("empty_data", "empty_data", 1, 0.0, 1.0, 1, 0.0, 1.0);
	}
	// Make sure data is setup
	if (!fDataHist) {
	NUIS_ABORT("No data has been setup inside " << fName << " constructor!");
	}

	// Make sure covariances are setup
	if (!fFullCovar) {
	fIsDiag = true;
	SetCovarFromDiagonal(fDataHist);
	}

	if (!covar) {
	covar = StatUtils::GetInvert(fFullCovar);
	}

	if (!fDecomp) {
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	}

	// Setup fMCHist from data
	fMCHist = (TH2D *)fDataHist->Clone();
	fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
	(fSettings.GetFullTitles()).c_str());
	fMCHist->Reset();

	// Setup fMCFine
	fMCFine = new TH2D(
	"mcfine", "mcfine", fDataHist->GetNbinsX() * 6,
	fMCHist->GetXaxis()->GetBinLowEdge(1),
	fMCHist->GetXaxis()->GetBinLowEdge(fDataHist->GetNbinsX() + 1),
	fDataHist->GetNbinsY() * 6, fMCHist->GetYaxis()->GetBinLowEdge(1),
	fMCHist->GetYaxis()->GetBinLowEdge(fDataHist->GetNbinsY() + 1));

	fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(),
	(fSettings.GetFullTitles()).c_str());
	fMCFine->Reset();

	// Setup MC Stat
	fMCStat = (TH2D *)fMCHist->Clone();
	fMCStat->Reset();

	// Search drawopts for possible types to include by default
	std::string drawopts = FitPar::Config().GetParS("drawopts");
	if (drawopts.find("MODES") != std::string::npos) {
	fMCHist_Modes = new TrueModeStack((fSettings.GetName() + "_MODES").c_str(),
	("True Channels"), fMCHist);
	SetAutoProcessTH1(fMCHist_Modes);
	}

	// Setup bin masks using sample name
	if (fIsMask) {

	std::string curname = fName;
	std::string origname = fSettings.GetS("originalname");

	// Check rename.mask
	std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask");

	// Check origname.mask
	if (maskloc.empty())
	maskloc = FitPar::Config().GetParDIR(origname + ".mask");

	// Check database
	if (maskloc.empty()) {
	maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask";
	}

	// Setup Bin Mask
	SetBinMask(maskloc);
	}

	if (fScaleFactor < 0) {
	NUIS_ERR(FTL, "I found a negative fScaleFactor in " << __FILE__ << ":"
	<< __LINE__);
	NUIS_ERR(FTL, "fScaleFactor = " << fScaleFactor);
	NUIS_ABORT("EXITING");
	}

	// Create and fill Weighted Histogram
	if (!fMCWeighted) {

	fMCWeighted = (TH2D *)fMCHist->Clone();
	fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(),
	(fName + "_MCWGHTS" + fPlotTitles).c_str());
	fMCWeighted->GetYaxis()->SetTitle("Weighted Events");
	}

	if (!fMapHist)
	fMapHist = StatUtils::GenerateMap(fDataHist);
	}

	//********************************************************************
	void Measurement2D::SetFitOptions(std::string opt) {
	//********************************************************************

	// Do nothing if default given
	if (opt == "DEFAULT")
	return;

	// CHECK Conflicting Fit Options
	std::vector<std::string> fit_option_allow =
	GeneralUtils::ParseToStr(fAllowedTypes, "/");

	for (UInt_t i = 0; i < fit_option_allow.size(); i++) {
	std::vector<std::string> fit_option_section =
	GeneralUtils::ParseToStr(fit_option_allow.at(i), ",");

	bool found_option = false;

	for (UInt_t j = 0; j < fit_option_section.size(); j++) {
	std::string av_opt = fit_option_section.at(j);

	if (!found_option and opt.find(av_opt) != std::string::npos) {
	found_option = true;

	} else if (found_option and opt.find(av_opt) != std::string::npos) {
	NUIS_ABORT("ERROR: Conflicting fit options provided: "
	<< opt << std::endl
	<< "Conflicting group = " << fit_option_section.at(i)
	<< std::endl
	<< "You should only supply one of these options in card file.");
	}
	}
	}

	// Check all options are allowed
	std::vector<std::string> fit_options_input =
	GeneralUtils::ParseToStr(opt, "/");
	for (UInt_t i = 0; i < fit_options_input.size(); i++) {
	if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) {
	NUIS_ERR(FTL, "ERROR: Fit Option '"
	<< fit_options_input.at(i)
	<< "' Provided is not allowed for this measurement.");
	NUIS_ERR(FTL, "Fit Options should be provided as a '/' seperated list "
	"(e.g. FREE/DIAG/NORM)");
	NUIS_ABORT("Available options for " << fName << " are '" << fAllowedTypes
	<< "'");
	}
	}

	// Set TYPE
	fFitType = opt;

	// FIX,SHAPE,FREE
	if (opt.find("FIX") != std::string::npos) {
	fIsFree = fIsShape = false;
	fIsFix = true;
	} else if (opt.find("SHAPE") != std::string::npos) {
	fIsFree = fIsFix = false;
	fIsShape = true;
	} else if (opt.find("FREE") != std::string::npos) {
	fIsFix = fIsShape = false;
	fIsFree = true;
	}

	// DIAG,FULL (or default to full)
	if (opt.find("DIAG") != std::string::npos) {
	fIsDiag = true;
	fIsFull = false;
	} else if (opt.find("FULL") != std::string::npos) {
	fIsDiag = false;
	fIsFull = true;
	}

	// CHI2/LL (OTHERS?)
	if (opt.find("LOG") != std::string::npos) {
	fIsChi2 = false;

	NUIS_ERR(FTL, "No other LIKELIHOODS properly supported!");
	NUIS_ABORT("Try to use a chi2!");

	} else {
	fIsChi2 = true;
	}

	// EXTRAS
	if (opt.find("RAW") != std::string::npos)
	fIsRawEvents = true;
	if (opt.find("DIF") != std::string::npos)
	fIsDifXSec = true;
	if (opt.find("ENU1D") != std::string::npos)
	fIsEnu1D = true;
	if (opt.find("NORM") != std::string::npos)
	fAddNormPen = true;
	if (opt.find("MASK") != std::string::npos)
	fIsMask = true;

	// Set TYPE
	fFitType = opt;

	// FIX,SHAPE,FREE
	if (opt.find("FIX") != std::string::npos) {
	fIsFree = fIsShape = false;
	fIsFix = true;
	} else if (opt.find("SHAPE") != std::string::npos) {
	fIsFree = fIsFix = false;
	fIsShape = true;
	} else if (opt.find("FREE") != std::string::npos) {
	fIsFix = fIsShape = false;
	fIsFree = true;
	}

	// DIAG,FULL (or default to full)
	if (opt.find("DIAG") != std::string::npos) {
	fIsDiag = true;
	fIsFull = false;
	} else if (opt.find("FULL") != std::string::npos) {
	fIsDiag = false;
	fIsFull = true;
	}

	// CHI2/LL (OTHERS?)
	if (opt.find("LOG") != std::string::npos)
	fIsChi2 = false;
	else
	fIsChi2 = true;

	// EXTRAS
	if (opt.find("RAW") != std::string::npos)
	fIsRawEvents = true;
	if (opt.find("DIF") != std::string::npos)
	fIsDifXSec = true;
	if (opt.find("ENU1D") != std::string::npos)
	fIsEnu = true;
	if (opt.find("NORM") != std::string::npos)
	fAddNormPen = true;
	if (opt.find("MASK") != std::string::npos)
	fIsMask = true;

	fIsProjFitX = (opt.find("FITPROJX") != std::string::npos);
	fIsProjFitY = (opt.find("FITPROJY") != std::string::npos);

	return;
	};

	/*
	Reconfigure LOOP
	*/
	//********************************************************************
	void Measurement2D::ResetAll() {
	//********************************************************************

	fMCHist->Reset();
	fMCFine->Reset();
	fMCStat->Reset();

	return;
	};

	//********************************************************************
	void Measurement2D::FillHistograms() {
	//********************************************************************

	if (Signal) {
	fMCHist->Fill(fXVar, fYVar, Weight);
	fMCFine->Fill(fXVar, fYVar, Weight);
	fMCStat->Fill(fXVar, fYVar, 1.0);

	if (fMCHist_Modes)
	fMCHist_Modes->Fill(Mode, fXVar, fYVar, Weight);
	}

	return;
	};

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

	// Fill MCWeighted;
	// for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
	// fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1));
	// fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1));
	// }

	// Setup Stat ratios for MC and MC Fine
	double *statratio = new double[fMCHist->GetNbinsX()];
	for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
	if (fMCHist->GetBinContent(i + 1) != 0) {
	statratio[i] =
	fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1);
	} else {
	statratio[i] = 0.0;
	}
	}

	double *statratiofine = new double[fMCFine->GetNbinsX()];
	for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
	if (fMCFine->GetBinContent(i + 1) != 0) {
	statratiofine[i] =
	fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1);
	} else {
	statratiofine[i] = 0.0;
	}
	}

	// Scaling for raw event rates
	if (fIsRawEvents) {
	double datamcratio = fDataHist->Integral() / fMCHist->Integral();

	fMCHist->Scale(datamcratio);
	fMCFine->Scale(datamcratio);

	if (fMCHist_Modes)
	fMCHist_Modes->Scale(datamcratio);

	// Scaling for XSec as function of Enu
	} else if (fIsEnu1D) {

	PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(),
	GetEventHistogram(), fScaleFactor);

	PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(),
	GetEventHistogram(), fScaleFactor);

	// if (fMCHist_Modes) {
	// PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(),
	// GetEventHistogram(), fScaleFactor,
	// fNEvents);
	// }

	// Any other differential scaling
	} else {
	fMCHist->Scale(fScaleFactor, "width");
	fMCFine->Scale(fScaleFactor, "width");

	// if (fMCHist_Modes) fMCHist_Modes->Scale(fScaleFactor, "width");
	}

	// Proper error scaling - ROOT Freaks out with xsec weights sometimes
	for (int i = 0; i < fMCStat->GetNbinsX(); i++) {
	fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]);
	}

	for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
	fMCFine->SetBinError(i + 1,
	fMCFine->GetBinContent(i + 1) * statratiofine[i]);
	}

	// Clean up
	delete statratio;
	delete statratiofine;

	return;
	};

	//********************************************************************
	void Measurement2D::ApplyNormScale(double norm) {
	//********************************************************************

	fCurrentNorm = norm;

	fMCHist->Scale(1.0 / norm);
	fMCFine->Scale(1.0 / norm);

	return;
	};

	/*
	Statistic Functions - Outsources to StatUtils
	*/

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

	// Just incase it has gone...
	if (!fDataHist)
	return -1;

	int nDOF = 0;

	// If datahist has no errors make sure we don't include those bins as they are
	// not data points
	for (int xBin = 0; xBin < fDataHist->GetNbinsX(); ++xBin) {
	for (int yBin = 0; yBin < fDataHist->GetNbinsY(); ++yBin) {
	if (fDataHist->GetBinError(xBin + 1, yBin + 1) != 0)
	++nDOF;
	}
	}

	// Account for possible bin masking
	int nMasked = 0;
	if (fMaskHist and fIsMask)
	if (fMaskHist->Integral() > 0)
	for (int xBin = 0; xBin < fMaskHist->GetNbinsX() + 1; ++xBin)
	for (int yBin = 0; yBin < fMaskHist->GetNbinsY() + 1; ++yBin)
	if (fMaskHist->GetBinContent(xBin, yBin) > 0.5)
	++nMasked;

	// Take away those masked DOF
	if (fIsMask) {
	nDOF -= nMasked;
	}

	return nDOF;
	}

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

	// If this is for a ratio, there is no data histogram to compare to!
	if (fNoData \|\| !fDataHist)
	return 0.;

	// Fix weird masking bug
	if (!fIsMask) {
	if (fMaskHist) {
	fMaskHist = NULL;
	}
	} else {
	if (fMaskHist) {
	PlotUtils::MaskBins(fMCHist, fMaskHist);
	}
	}

	// if (fIsProjFitX or fIsProjFitY) return GetProjectedChi2();

	// Scale up the results to match each other (Not using width might be
	// inconsistent with Meas1D)
	double scaleF = fDataHist->Integral() / fMCHist->Integral();
	if (fIsShape) {
	fMCHist->Scale(scaleF);
	fMCFine->Scale(scaleF);
	// PlotUtils::ScaleNeutModeArray((TH1**)fMCHist_PDG, scaleF);
	}

	if (!fMapHist)
	fMapHist = StatUtils::GenerateMap(fDataHist);

	// Get the chi2 from either covar or diagonals
	double chi2 = 0.0;

	if (fIsChi2) {
	if (fIsDiag) {
	chi2 =
	StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMapHist, fMaskHist);
	} else {
	chi2 = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMapHist,
	fMaskHist, fResidualHist);
	if (fChi2LessBinHist) {
	TH2I *binmask = new TH2I("mask", "", fDataHist->GetNbinsX(), 0,
	fDataHist->GetNbinsX(), fDataHist->GetNbinsY(),
	0, fDataHist->GetNbinsY());
	binmask->SetDirectory(NULL);
	for (int xi = 0; xi < fDataHist->GetNbinsX(); ++xi) {
	for (int yi = 0; yi < fDataHist->GetNbinsY(); ++yi) {
	binmask->Reset();
	binmask->SetBinContent(xi + 1, yi + 1, 1);
	fChi2LessBinHist->SetBinContent(
	xi + 1, yi + 1,
	StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMapHist,
	binmask));
	}
	}
	delete binmask;
	}
	}
	}

	// Add a normal penalty term
	if (fAddNormPen) {
	chi2 +=
	(1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) / (fNormError * fNormError);
	NUIS_LOG(REC, "Norm penalty = " << (1 - (fCurrentNorm)) * (1 - (fCurrentNorm)) /
	(fNormError * fNormError));
	}

	// Adjust the shape back to where it was.
	if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) {
	fMCHist->Scale(1. / scaleF);
	fMCFine->Scale(1. / scaleF);
	}

	fLikelihood = chi2;

	return chi2;
	}

	/*
	Fake Data Functions
	*/
	//********************************************************************
	void Measurement2D::SetFakeDataValues(std::string fakeOption) {
	//********************************************************************

	// Setup original/datatrue
	TH2D tempdata = (TH2D )fDataHist->Clone();

	if (!fIsFakeData) {
	fIsFakeData = true;

	// Make a copy of the original data histogram.
	if (!fDataOrig)
	fDataOrig = (TH2D *)fDataHist->Clone((fName + "_data_original").c_str());

	} else {
	ResetFakeData();
	}

	// Setup Inputs
	fFakeDataInput = fakeOption;
	NUIS_LOG(SAM, "Setting fake data from : " << fFakeDataInput);

	// From MC
	if (fFakeDataInput.compare("MC") == 0) {
	fDataHist = (TH2D *)fMCHist->Clone((fName + "_MC").c_str());

	// Fake File
	} else {
	if (!fFakeDataFile)
	fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ");
	fDataHist = (TH2D *)fFakeDataFile->Get((fName + "_MC").c_str());
	}

	// Setup Data Hist
	fDataHist->SetNameTitle((fName + "_FAKE").c_str(),
	(fName + fPlotTitles).c_str());

	// Replace Data True
	if (fDataTrue)
	delete fDataTrue;
	fDataTrue = (TH2D *)fDataHist->Clone();
	fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(),
	(fName + fPlotTitles).c_str());

	// Make a new covariance for fake data hist.
	int nbins = fDataHist->GetNbinsX() * fDataHist->GetNbinsY();
	double alpha_i = 0.0;
	double alpha_j = 0.0;

	for (int i = 0; i < nbins; i++) {
	for (int j = 0; j < nbins; j++) {
	if (tempdata->GetBinContent(i + 1) && tempdata->GetBinContent(j + 1)) {
	alpha_i =
	fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1);
	alpha_j =
	fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1);
	} else {
	alpha_i = 0.0;
	alpha_j = 0.0;
	}

	(fFullCovar)(i, j) = alpha_i alpha_j * (*fFullCovar)(i, j);
	}
	}

	// Setup Covariances
	if (covar)
	delete covar;
	covar = StatUtils::GetInvert(fFullCovar);

	if (fDecomp)
	delete fDecomp;
	fDecomp = StatUtils::GetInvert(fFullCovar);

	delete tempdata;

	return;
	};

	//********************************************************************
	void Measurement2D::ResetFakeData() {
	//********************************************************************

	if (fIsFakeData) {
	if (fDataHist)
	delete fDataHist;
	fDataHist =
	(TH2D *)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str());
	}
	}

	//********************************************************************
	void Measurement2D::ResetData() {
	//********************************************************************

	if (fIsFakeData) {
	if (fDataHist)
	delete fDataHist;
	fDataHist =
	(TH2D *)fDataOrig->Clone((fSettings.GetName() + "_data").c_str());
	}

	fIsFakeData = false;
	}

	//********************************************************************
	void Measurement2D::ThrowCovariance() {
	//********************************************************************

	// Take a fDecomposition and use it to throw the current dataset.
	// Requires fDataTrue also be set incase used repeatedly.

	if (fDataHist)
	delete fDataHist;
	fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);

	return;
	};

	//********************************************************************
	void Measurement2D::ThrowDataToy() {
	//********************************************************************
	if (!fDataTrue)
	fDataTrue = (TH2D *)fDataHist->Clone();
	if (fMCHist)
	delete fMCHist;
	fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
	}

	/*
	Access Functions
	*/

	//********************************************************************
	TH2D *Measurement2D::GetMCHistogram() {
	//********************************************************************

	if (!fMCHist)
	return fMCHist;

	std::ostringstream chi2;
	chi2 << std::setprecision(5) << this->GetLikelihood();

	int linecolor = kRed;
	int linestyle = 1;
	int linewidth = 1;

	int fillcolor = 0;
	int fillstyle = 1001;

	if (fSettings.Has("linecolor"))
	linecolor = fSettings.GetI("linecolor");
	if (fSettings.Has("linestyle"))
	linestyle = fSettings.GetI("linestyle");
	if (fSettings.Has("linewidth"))
	linewidth = fSettings.GetI("linewidth");

	if (fSettings.Has("fillcolor"))
	fillcolor = fSettings.GetI("fillcolor");
	if (fSettings.Has("fillstyle"))
	fillstyle = fSettings.GetI("fillstyle");

	fMCHist->SetTitle(chi2.str().c_str());

	fMCHist->SetLineColor(linecolor);
	fMCHist->SetLineStyle(linestyle);
	fMCHist->SetLineWidth(linewidth);

	fMCHist->SetFillColor(fillcolor);
	fMCHist->SetFillStyle(fillstyle);

	return fMCHist;
	};

	//********************************************************************
	TH2D *Measurement2D::GetDataHistogram() {
	//********************************************************************

	if (!fDataHist)
	return fDataHist;

	int datacolor = kBlack;
	int datastyle = 1;
	int datawidth = 1;

	if (fSettings.Has("datacolor"))
	datacolor = fSettings.GetI("datacolor");
	if (fSettings.Has("datastyle"))
	datastyle = fSettings.GetI("datastyle");
	if (fSettings.Has("datawidth"))
	datawidth = fSettings.GetI("datawidth");

	fDataHist->SetLineColor(datacolor);
	fDataHist->SetLineWidth(datawidth);
	fDataHist->SetMarkerStyle(datastyle);

	return fDataHist;
	};

	/*
	Write Functions
	*/

	// Save all the histograms at once
	//********************************************************************
	void Measurement2D::Write(std::string drawOpt) {
	//********************************************************************

	// Get Draw Options
	drawOpt = FitPar::Config().GetParS("drawopts");

	// Write Settigns
	if (drawOpt.find("SETTINGS") != std::string::npos) {
	fSettings.Set("#chi^{2}", fLikelihood);
	fSettings.Set("NDOF", this->GetNDOF());
	fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF());
	fSettings.Write();
	}

	if (fIsChi2 && !fIsDiag) {
	fResidualHist = (TH2D *)fMCHist->Clone((fName + "_RESIDUAL").c_str());
	fResidualHist->GetYaxis()->SetTitle("#Delta#chi^{2}");
	fResidualHist->Reset();

	fChi2LessBinHist = (TH2D *)fMCHist->Clone((fName + "_Chi2NMinusOne").c_str());
	fChi2LessBinHist->GetYaxis()->SetTitle("Total #chi^{2} without bin_{i}");
	fChi2LessBinHist->Reset();

	(void)GetLikelihood();

	fResidualHist->Write();
	fChi2LessBinHist->Write();
	}

	// // Likelihood residual plots
	// if (drawOpt.find("RESIDUAL") != std::string::npos) {
	// WriteResidualPlots();
	//}

	// // RATIO
	// if (drawOpt.find("CANVMC") != std::string::npos) {
	// TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist);
	// c1->Write();
	// delete c1;
	// }

	// // PDG
	// if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) {
	// TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes);
	// c2->Write();
	// delete c2;
	// }

	/// 2D VERSION
	// If null pointer return
	if (!fMCHist and !fDataHist) {
	NUIS_LOG(SAM, fName << "Incomplete histogram set!");
	return;
	}

	// Config::Get().out->cd();

	// Get Draw Options
	drawOpt = FitPar::Config().GetParS("drawopts");
	bool drawData = (drawOpt.find("DATA") != std::string::npos);
	bool drawNormal = (drawOpt.find("MC") != std::string::npos);
	bool drawEvents = (drawOpt.find("EVT") != std::string::npos);
	bool drawFine = (drawOpt.find("FINE") != std::string::npos);
	bool drawRatio = (drawOpt.find("RATIO") != std::string::npos);
	// bool drawModes = (drawOpt.find("MODES") != std::string::npos);
	bool drawShape = (drawOpt.find("SHAPE") != std::string::npos);
	bool residual = (drawOpt.find("RESIDUAL") != std::string::npos);
	bool drawMatrix = (drawOpt.find("MATRIX") != std::string::npos);
	bool drawFlux = (drawOpt.find("FLUX") != std::string::npos);
	bool drawMask = (drawOpt.find("MASK") != std::string::npos);
	bool drawMap = (drawOpt.find("MAP") != std::string::npos);
	bool drawProj = (drawOpt.find("PROJ") != std::string::npos);
	// bool drawCanvPDG = (drawOpt.find("CANVPDG") != std::string::npos);
	bool drawCov = (drawOpt.find("COV") != std::string::npos);
	bool drawSliceMC = (drawOpt.find("CANVSLICEMC") != std::string::npos);
	bool drawWeighted =
	(drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted);

	if (FitPar::Config().GetParB("EventManager")) {
	drawFlux = false;
	drawEvents = false;
	}

	// Save standard plots
	if (drawData) {
	GetDataList().at(0)->Write();
	// Generate a simple map
	if (!fMapHist)
	fMapHist = StatUtils::GenerateMap(fDataHist);
	// Convert to 1D Lists
	TH1D *data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist);
	data_1D->Write();
	delete data_1D;
	}
	if (drawNormal) {
	GetMCList().at(0)->Write();
	if (!fMapHist)
	fMapHist = StatUtils::GenerateMap(fDataHist);
	TH1D *mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist);
	mc_1D->SetLineColor(kRed);
	mc_1D->Write();
	delete mc_1D;
	}

	// Write Weighted Histogram
	if (drawWeighted)
	fMCWeighted->Write();

	if (drawCov) {
	TH2D(*fFullCovar).Write((fName + "_COV").c_str());
	}

	if (drawOpt.find("INVCOV") != std::string::npos) {
	TH2D(*covar).Write((fName + "_INVCOV").c_str());
	}

	// Save only mc and data if splines
	if (fEventType == 4 or fEventType == 3) {
	return;
	}

	// Draw Extra plots
	if (drawFine)
	this->GetFineList().at(0)->Write();
	if (drawFlux and GetFluxHistogram()) {
	GetFluxHistogram()->Write();
	}
	if (drawEvents and GetEventHistogram()) {
	GetEventHistogram()->Write();
	}
	if (fIsMask and drawMask) {
	fMaskHist->Write((fName + "_MSK").c_str()); //< save mask
	TH1I *mask_1D = StatUtils::MapToMask(fMaskHist, fMapHist);
	if (mask_1D) {
	mask_1D->Write();

	TMatrixDSym *calc_cov =
	StatUtils::ApplyInvertedMatrixMasking(covar, mask_1D);
	TH1D *data_1D = StatUtils::MapToTH1D(fDataHist, fMapHist);
	TH1D *mc_1D = StatUtils::MapToTH1D(fMCHist, fMapHist);
	TH1D *calc_data = StatUtils::ApplyHistogramMasking(data_1D, mask_1D);
	TH1D *calc_mc = StatUtils::ApplyHistogramMasking(mc_1D, mask_1D);

	TH2D bin_cov = new TH2D(calc_cov);

	bin_cov->Write();
	calc_data->Write();
	calc_mc->Write();

	delete mask_1D;
	delete calc_cov;
	delete calc_data;
	delete calc_mc;
	delete bin_cov;
	delete data_1D;
	delete mc_1D;
	}
	}
	if (drawMap)
	fMapHist->Write((fName + "_MAP").c_str()); //< save map

	// // Save neut stack
	// if (drawModes) {
	// THStack combo_fMCHist_PDG = PlotUtils::GetNeutModeStack(
	// (fName + "_MC_PDG").c_str(),
	// (TH1**)fMCHist_PDG, 0);
	// combo_fMCHist_PDG.Write();
	// }

	// Save Matrix plots
	if (drawMatrix and fFullCovar and covar and fDecomp) {
	TH2D cov = TH2D((*fFullCovar));
	cov.SetNameTitle((fName + "_cov").c_str(),
	(fName + "_cov;Bins; Bins;").c_str());
	cov.Write();

	TH2D covinv = TH2D((*this->covar));
	covinv.SetNameTitle((fName + "_covinv").c_str(),
	(fName + "_cov;Bins; Bins;").c_str());
	covinv.Write();

	TH2D covdec = TH2D((*fDecomp));
	covdec.SetNameTitle((fName + "_covdec").c_str(),
	(fName + "_cov;Bins; Bins;").c_str());
	covdec.Write();
	}

	// Save ratio plots if required
	if (drawRatio) {
	// Needed for error bars
	for (int i = 0; i < fMCHist->GetNbinsX() * fMCHist->GetNbinsY(); i++)
	fMCHist->SetBinError(i + 1, 0.0);

	fDataHist->GetSumw2();
	fMCHist->GetSumw2();

	// Create Ratio Histograms
	TH2D dataRatio = (TH2D )fDataHist->Clone((fName + "_data_RATIO").c_str());
	TH2D mcRatio = (TH2D )fMCHist->Clone((fName + "_MC_RATIO").c_str());

	mcRatio->Divide(fMCHist);
	dataRatio->Divide(fMCHist);

	// Cancel bin errors on MC
	for (int i = 0; i < mcRatio->GetNbinsX() * mcRatio->GetNbinsY(); i++) {
	mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) /
	fMCHist->GetBinContent(i + 1));
	}

	mcRatio->SetMinimum(0);
	mcRatio->SetMaximum(2);
	dataRatio->SetMinimum(0);
	dataRatio->SetMaximum(2);

	mcRatio->Write();
	dataRatio->Write();

	delete mcRatio;
	delete dataRatio;
	}

	// Save Shape Plots if required
	if (drawShape) {
	// Create Shape Histogram
	TH2D mcShape = (TH2D )fMCHist->Clone((fName + "_MC_SHAPE").c_str());

	double shapeScale = 1.0;
	if (fIsRawEvents) {
	shapeScale = fDataHist->Integral() / fMCHist->Integral();
	} else {
	shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
	}

	mcShape->Scale(shapeScale);

	mcShape->SetLineWidth(3);
	mcShape->SetLineStyle(7); // dashes

	mcShape->Write();

	// Save shape ratios
	if (drawRatio) {
	// Needed for error bars
	mcShape->GetSumw2();

	// Create shape ratio histograms
	TH2D *mcShapeRatio =
	(TH2D *)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str());
	TH2D *dataShapeRatio =
	(TH2D *)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str());

	// Divide the histograms
	mcShapeRatio->Divide(mcShape);
	dataShapeRatio->Divide(mcShape);

	// Colour the shape ratio plots
	mcShapeRatio->SetLineWidth(3);
	mcShapeRatio->SetLineStyle(7); // dashes

	mcShapeRatio->Write();
	dataShapeRatio->Write();

	delete mcShapeRatio;
	delete dataShapeRatio;
	}

	delete mcShape;
	}

	// Save residual calculations of what contributed to the chi2 values.
	if (residual) {
	}

	if (fIsProjFitX \|\| fIsProjFitY \|\| drawProj) {
	// If not already made, make the projections
	if (!fMCHist_X) {
	PlotUtils::MatchEmptyBins(fDataHist, fMCHist);

	fMCHist_X = PlotUtils::GetProjectionX(fMCHist, fMaskHist);
	fMCHist_Y = PlotUtils::GetProjectionY(fMCHist, fMaskHist);

	fDataHist_X = PlotUtils::GetProjectionX(fDataHist, fMaskHist);
	fDataHist_Y = PlotUtils::GetProjectionY(fDataHist, fMaskHist);

	// This is not the correct way of doing it
	// double chi2X = StatUtils::GetChi2FromDiag(fDataHist_X, fMCHist_X);
	// double chi2Y = StatUtils::GetChi2FromDiag(fDataHist_Y, fMCHist_Y);

	// fMCHist_X->SetTitle(Form("%f", chi2X));
	// fMCHist_Y->SetTitle(Form("%f", chi2Y));
	}

	// Save the histograms
	fDataHist_X->Write();
	fMCHist_X->Write();

	fDataHist_Y->Write();
	fMCHist_Y->Write();
	}

	if (drawSliceMC) {
	TCanvas *c1 = new TCanvas((fName + "_MC_CANV_Y").c_str(),
	(fName + "_MC_CANV_Y").c_str(), 1024, 1024);

	c1->Divide(2, int(fDataHist->GetNbinsY() / 3. + 1));
	TH2D mcShape = (TH2D )fMCHist->Clone((fName + "_MC_SHAPE").c_str());
	double shapeScale =
	fDataHist->Integral("width") / fMCHist->Integral("width");
	mcShape->Scale(shapeScale);
	mcShape->SetLineStyle(7);

	c1->cd(1);
	TLegend *leg = new TLegend(0.0, 0.0, 1.0, 1.0);
	leg->AddEntry(fDataHist, (fName + " Data").c_str(), "lep");
	leg->AddEntry(fMCHist, (fName + " MC").c_str(), "l");
	leg->AddEntry(mcShape, (fName + " Shape").c_str(), "l");
	leg->SetLineColor(0);
	leg->SetLineStyle(0);
	leg->SetFillColor(0);
	leg->SetLineStyle(0);
	leg->Draw("SAME");

	// Make Y slices
	for (int i = 1; i < fDataHist->GetNbinsY() + 1; i++) {
	c1->cd(i + 1);
	TH1D *fDataHist_SliceY = PlotUtils::GetSliceY(fDataHist, i);
	fDataHist_SliceY->Draw("E1");

	TH1D *fMCHist_SliceY = PlotUtils::GetSliceY(fMCHist, i);
	fMCHist_SliceY->Draw("SAME");

	TH1D *mcShape_SliceY = PlotUtils::GetSliceY(mcShape, i);
	mcShape_SliceY->Draw("SAME");
	mcShape_SliceY->SetLineStyle(mcShape->GetLineStyle());
	}
	c1->Write();
	delete c1;
	delete leg;

	TCanvas *c2 = new TCanvas((fName + "_MC_CANV_X").c_str(),
	(fName + "_MC_CANV_X").c_str(), 1024, 1024);

	c2->Divide(2, int(fDataHist->GetNbinsX() / 3. + 1));
	mcShape = (TH2D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
	shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
	mcShape->Scale(shapeScale);
	mcShape->SetLineStyle(7);

	c2->cd(1);
	TLegend *leg2 = new TLegend(0.0, 0.0, 1.0, 1.0);
	leg2->AddEntry(fDataHist, (fName + " Data").c_str(), "lep");
	leg2->AddEntry(fMCHist, (fName + " MC").c_str(), "l");
	leg2->AddEntry(mcShape, (fName + " Shape").c_str(), "l");
	leg2->SetLineColor(0);
	leg2->SetLineStyle(0);
	leg2->SetFillColor(0);
	leg2->SetLineStyle(0);
	leg2->Draw("SAME");

	// Make Y slices
	for (int i = 1; i < fDataHist->GetNbinsX() + 1; i++) {
	c2->cd(i + 1);
	TH1D *fDataHist_SliceX = PlotUtils::GetSliceX(fDataHist, i);
	fDataHist_SliceX->Draw("E1");

	TH1D *fMCHist_SliceX = PlotUtils::GetSliceX(fMCHist, i);
	fMCHist_SliceX->Draw("SAME");

	TH1D *mcShape_SliceX = PlotUtils::GetSliceX(mcShape, i);
	mcShape_SliceX->Draw("SAME");
	mcShape_SliceX->SetLineStyle(mcShape->GetLineStyle());
	}

	c2->Write();
	delete c2;
	delete leg2;
	}

	// Write Extra Histograms
	AutoWriteExtraTH1();
	WriteExtraHistograms();

	// Returning
	NUIS_LOG(SAM, "Written Histograms: " << fName);
	return;
	}

	/*
	Setup Functions
	*/
	//********************************************************************
	void Measurement2D::SetupMeasurement(std::string inputfile, std::string type,
	FitWeight *rw, std::string fkdt) {
	//********************************************************************

	// Check if name contains Evt, indicating that it is a raw number of events
	// measurements and should thus be treated as once
	fIsRawEvents = false;
	if ((fName.find("Evt") != std::string::npos) && fIsRawEvents == false) {
	fIsRawEvents = true;
	NUIS_LOG(SAM, "Found event rate measurement but fIsRawEvents == false!");
	NUIS_LOG(SAM, "Overriding this and setting fIsRawEvents == true!");
	}

	fIsEnu = false;
	if ((fName.find("XSec") != std::string::npos) &&
	(fName.find("Enu") != std::string::npos)) {
	fIsEnu = true;
	NUIS_LOG(SAM, "::" << fName << "::");
	NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, "
	"not flux averaged!");
	if (FitPar::Config().GetParB("EventManager")) {
	NUIS_ERR(FTL, "Enu Measurements do not yet work with the Event Manager!");
	NUIS_ERR(FTL, "If you want decent flux unfolded results please run in "
	"series mode (-q EventManager=0)");
	sleep(2);
	throw;
	}
	}

	if (fIsEnu && fIsRawEvents) {
	NUIS_ERR(FTL, "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
	"really correct?!");
	NUIS_ERR(FTL, "Check experiment constructor for " << fName
	<< " and correct this!");
	NUIS_ABORT("I live in " << __FILE__ << ":" << __LINE__);
	}

	// Reset everything to NULL
	fRW = rw;

	// Setting up 2D Inputs
	this->SetupInputs(inputfile);

	// Set Default Options
	SetFitOptions(fDefaultTypes);

	// Set Passed Options
	SetFitOptions(type);
	}

	//********************************************************************
	void Measurement2D::SetupDefaultHist() {
	//********************************************************************

	// Setup fMCHist
	fMCHist = (TH2D *)fDataHist->Clone();
	fMCHist->SetNameTitle((fName + "_MC").c_str(),
	(fName + "_MC" + fPlotTitles).c_str());

	// Setup fMCFine
	Int_t nBinsX = fMCHist->GetNbinsX();
	Int_t nBinsY = fMCHist->GetNbinsY();
	fMCFine = new TH2D((fName + "_MC_FINE").c_str(),
	(fName + "_MC_FINE" + fPlotTitles).c_str(), nBinsX * 3,
	fMCHist->GetXaxis()->GetBinLowEdge(1),
	fMCHist->GetXaxis()->GetBinLowEdge(nBinsX + 1), nBinsY * 3,
	fMCHist->GetYaxis()->GetBinLowEdge(1),
	fMCHist->GetYaxis()->GetBinLowEdge(nBinsY + 1));

	// Setup MC Stat
	fMCStat = (TH2D *)fMCHist->Clone();
	fMCStat->Reset();

	// Setup the NEUT Mode Array
	// PlotUtils::CreateNeutModeArray(fMCHist, (TH1**)fMCHist_PDG);

	// Setup bin masks using sample name
	if (fIsMask) {
	std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask");
	if (maskloc.empty()) {
	maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask";
	}

	SetBinMask(maskloc);
	}

	return;
	}

	//********************************************************************
	void Measurement2D::SetDataValues(std::string dataFile, std::string TH2Dname) {
	//********************************************************************

	if (dataFile.find(".root") == std::string::npos) {
	NUIS_ERR(FTL, "Error! " << dataFile << " is not a .root file");
	NUIS_ERR(FTL, "Currently only .root file reading is supported (MiniBooNE "
	"CC1pi+ 2D), but implementing .txt should be dirt easy");
	NUIS_ABORT("See me at " << __FILE__ << ":" << __LINE__);

	} else {
	TFile *inFile = new TFile(dataFile.c_str(), "READ");
	fDataHist = (TH2D *)(inFile->Get(TH2Dname.c_str())->Clone());
	fDataHist->SetDirectory(0);

	fDataHist->SetNameTitle((fName + "_data").c_str(),
	(fName + "_MC" + fPlotTitles).c_str());

	delete inFile;
	}

	return;
	}

	//********************************************************************
	void Measurement2D::SetDataValues(std::string dataFile, double dataNorm,
	std::string errorFile, double errorNorm) {
	//********************************************************************

	// Make a counter to track the line number
	int yBin = 0;

	std::string line;
	std::ifstream data(dataFile.c_str(), std::ifstream::in);

	fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(),
	fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);

	if (data.is_open()) {
	NUIS_LOG(SAM, "Reading data from: " << dataFile.c_str());
	}

	while (std::getline(data >> std::ws, line, '\n')) {
	int xBin = 0;

	// Loop over entries and insert them into the histogram
	std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
	for (std::vector<double>::iterator iter = entries.begin();
	iter != entries.end(); iter++) {
	fDataHist->SetBinContent(xBin + 1, yBin + 1, (iter) dataNorm);
	xBin++;
	}
	yBin++;
	}

	yBin = 0;
	std::ifstream error(errorFile.c_str(), std::ifstream::in);

	if (error.is_open()) {
	NUIS_LOG(SAM, "Reading errors from: " << errorFile.c_str());
	}

	while (std::getline(error >> std::ws, line, '\n')) {
	int xBin = 0;

	// Loop over entries and insert them into the histogram
	std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
	for (std::vector<double>::iterator iter = entries.begin();
	iter != entries.end(); iter++) {
	fDataHist->SetBinError(xBin + 1, yBin + 1, (iter) errorNorm);
	xBin++;
	}
	yBin++;
	}

	return;
	};

	//********************************************************************
	void Measurement2D::SetDataValuesFromText(std::string dataFile,
	double dataNorm) {
	//********************************************************************

	fDataHist = new TH2D((fName + "_data").c_str(), (fName + fPlotTitles).c_str(),
	fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);

	NUIS_LOG(SAM, "Reading data from: " << dataFile);
	PlotUtils::Set2DHistFromText(dataFile, fDataHist, dataNorm, true);

	return;
	};

	//********************************************************************
	void Measurement2D::SetCovarMatrix(std::string covarFile) {
	//********************************************************************

	// Used to read a covariance matrix from a root file
	TFile *tempFile = new TFile(covarFile.c_str(), "READ");

	// Make plots that we want
	TH2D *covarPlot = new TH2D();
	TH2D *fFullCovarPlot = new TH2D();

	// Get covariance options for fake data studies
	std::string covName = "";
	std::string covOption = FitPar::Config().GetParS("throw_covariance");

	// Which matrix to get?
	if (fIsShape \|\| fIsFree)
	covName = "shp_";
	if (fIsDiag)
	covName += "diag";
	else
	covName += "full";

	covarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str());

	// Throw either the sub matrix or the full matrix
	if (!covOption.compare("SUB"))
	fFullCovarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str());
	else if (!covOption.compare("FULL"))
	fFullCovarPlot = (TH2D *)tempFile->Get("fullcov");
	else {
	NUIS_ERR(WRN, " Incorrect thrown_covariance option in parameters.");
	}

	// Bin masking?
	int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral());
	int covdim = int(fDataHist->GetNbinsX());

	// Make new covars
	this->covar = new TMatrixDSym(dim);
	fFullCovar = new TMatrixDSym(dim);
	fDecomp = new TMatrixDSym(dim);

	// Full covariance values
	int row, column = 0;
	row = 0;
	column = 0;
	for (Int_t i = 0; i < covdim; i++) {
	// masking can be dodgy
	// if (this->masked->GetBinContent(i+1) > 0) continue;
	for (Int_t j = 0; j < covdim; j++) {
	// if (this->masked->GetBinContent(j+1) > 0) continue;
	(*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1);
	(*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1);

	column++;
	}
	column = 0;
	row++;
	}

	// Set bin errors on data
	if (!fIsDiag) {
	for (Int_t i = 0; i < fDataHist->GetNbinsX(); i++) {
	fDataHist->SetBinError(
	i + 1, sqrt((covarPlot->GetBinContent(i + 1, i + 1))) * 1E-38);
	}
	}

	TDecompSVD LU = TDecompSVD(*this->covar);
	this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");

	tempFile->Close();
	delete tempFile;

	return;
	};

	//********************************************************************
	void Measurement2D::SetCovarMatrixFromText(std::string covarFile, int dim) {
	//********************************************************************

	// Make a counter to track the line number
	int row = 0;

	std::string line;
	std::ifstream covar(covarFile.c_str(), std::ifstream::in);

	this->covar = new TMatrixDSym(dim);
	fFullCovar = new TMatrixDSym(dim);
	if (covar.is_open()) {
	NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile);
	}
	while (std::getline(covar >> std::ws, line, '\n')) {
	int column = 0;

	// Loop over entries and insert them into matrix
	// Multiply by the errors to get the covariance, rather than the correlation
	// matrix
	std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
	for (std::vector<double>::iterator iter = entries.begin();
	iter != entries.end(); iter++) {
	double val = (iter) fDataHist->GetBinError(row + 1) * 1E38 *
	fDataHist->GetBinError(column + 1) * 1E38;
	(*this->covar)(row, column) = val;
	(*fFullCovar)(row, column) = val;

	column++;
	}

	row++;
	}

	// Robust matrix inversion method
	TDecompSVD LU = TDecompSVD(*this->covar);
	this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");

	return;
	};

	//********************************************************************
	void Measurement2D::SetCovarMatrixFromChol(std::string covarFile, int dim) {
	//********************************************************************

	// Make a counter to track the line number
	int row = 0;

	std::string line;
	std::ifstream covarread(covarFile.c_str(), std::ifstream::in);

	TMatrixD *newcov = new TMatrixD(dim, dim);

	if (covarread.is_open()) {
	NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile);
	}
	while (std::getline(covarread >> std::ws, line, '\n')) {
	int column = 0;

	// Loop over entries and insert them into matrix
	// Multiply by the errors to get the covariance, rather than the correlation
	// matrix
	std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
	for (std::vector<double>::iterator iter = entries.begin();
	iter != entries.end(); iter++) {
	(newcov)(row, column) = iter;
	column++;
	}

	row++;
	}
	covarread.close();

	// Form full covariance
	TMatrixD trans = (TMatrixD )(newcov)->Clone();
	trans->T();
	(trans) = (*newcov);

	fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");

	delete newcov;
	delete trans;

	// Robust matrix inversion method
	TDecompChol LU = TDecompChol(*this->fFullCovar);
	this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");

	return;
	};

	// //********************************************************************
	// void Measurement2D::SetMapValuesFromText(std::string dataFile) {
	// //********************************************************************

	// fMapHist = new TH2I((fName + "_map").c_str(), (fName +
	// fPlotTitles).c_str(),
	// fNDataPointsX - 1, fXBins, fNDataPointsY - 1, fYBins);

	// LOG(SAM) << "Reading map from: " << dataFile << std::endl;
	// PlotUtils::Set2DHistFromText(dataFile, fMapHist, 1.0);

	// return;
	// };
	diff --git a/src/Statistical/StatUtils.cxx b/src/Statistical/StatUtils.cxx
	index 586cbed..a12a28a 100644
	--- a/src/Statistical/StatUtils.cxx
	+++ b/src/Statistical/StatUtils.cxx
	@@ -1,1448 +1,1451 @@
	// 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 "StatUtils.h"
	#include "GeneralUtils.h"
	#include "NuisConfig.h"
	#include "TH1D.h"

	//*******************************************************************
	Double_t StatUtils::GetChi2FromDiag(TH1D data, TH1D mc, TH1I *mask) {
	//*******************************************************************

	Double_t Chi2 = 0.0;
	TH1D calc_data = (TH1D )data->Clone("calc_data");
	calc_data->SetDirectory(NULL);
	TH1D calc_mc = (TH1D )mc->Clone("calc_mc");
	calc_mc->SetDirectory(NULL);

	// Add MC Error to data if required
	if (FitPar::Config().GetParB("addmcerror")) {
	for (int i = 0; i < calc_data->GetNbinsX(); i++) {
	double dterr = calc_data->GetBinError(i + 1);
	double mcerr = calc_mc->GetBinError(i + 1);

	if (dterr > 0.0) {
	calc_data->SetBinError(i + 1, sqrt(dterr * dterr + mcerr * mcerr));
	}
	}
	}

	// Apply masking if required
	if (mask) {
	calc_data = ApplyHistogramMasking(data, mask);
	calc_data->SetDirectory(NULL);
	calc_mc = ApplyHistogramMasking(mc, mask);
	calc_mc->SetDirectory(NULL);
	}

	// Iterate over bins in X
	for (int i = 0; i < calc_data->GetNbinsX(); i++) {
	// Ignore bins with zero data or zero bin error
	if (calc_data->GetBinError(i + 1) <= 0.0 \|\|
	calc_data->GetBinContent(i + 1) == 0.0)
	continue;

	// Take mc data difference
	double diff =
	calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1);
	double err = calc_data->GetBinError(i + 1);
	Chi2 += (diff * diff) / (err * err);
	}

	// cleanup
	delete calc_data;
	delete calc_mc;

	return Chi2;
	};

	//*******************************************************************
	Double_t StatUtils::GetChi2FromDiag(TH2D data, TH2D mc, TH2I *map,
	TH2I *mask) {
	//*******************************************************************

	// Generate a simple map
	if (!map)
	map = GenerateMap(data);

	// Convert to 1D Histograms
	TH1D *data_1D = MapToTH1D(data, map);
	TH1D *mc_1D = MapToTH1D(mc, map);
	TH1I *mask_1D = MapToMask(mask, map);

	// Calculate 1D chi2 from 1D Plots
	Double_t Chi2 = StatUtils::GetChi2FromDiag(data_1D, mc_1D, mask_1D);

	// CleanUp
	delete data_1D;
	delete mc_1D;
	delete mask_1D;

	return Chi2;
	};

	//*******************************************************************
	Double_t StatUtils::GetChi2FromCov(TH1D data, TH1D mc, TMatrixDSym *invcov,
	TH1I *mask, double data_scale,
	double covar_scale, TH1D *outchi2perbin) {
	//*******************************************************************

	Double_t Chi2 = 0.0;
	- TMatrixDSym calc_cov = (TMatrixDSym )invcov->Clone();
	- TH1D calc_data = (TH1D )data->Clone();
	- TH1D calc_mc = (TH1D )mc->Clone();
	+ TMatrixDSym calc_cov = (TMatrixDSym )invcov->Clone("local_invcov");
	+ TH1D calc_data = (TH1D )data->Clone("local_data");
	+ TH1D calc_mc = (TH1D )mc->Clone("local_mc");
	+
	+ calc_data->SetDirectory(NULL);
	+ calc_mc->SetDirectory(NULL);

	// If a mask if applied we need to apply it before the matrix is inverted
	if (mask) {
	calc_cov = ApplyInvertedMatrixMasking(invcov, mask);
	calc_data = ApplyHistogramMasking(data, mask);
	calc_mc = ApplyHistogramMasking(mc, mask);
	}

	// Add MC Error to data if required
	if (FitPar::Config().GetParB("statutils.addmcerror")) {
	// Make temp cov
	TMatrixDSym *newcov = StatUtils::GetInvert(calc_cov);

	// Add MC err to diag
	for (int i = 0; i < calc_data->GetNbinsX(); i++) {
	double mcerr = calc_mc->GetBinError(i + 1) * sqrt(covar_scale);
	double oldval = (*newcov)(i, i);

	- NUIS_LOG(FIT,
	- "Adding cov stat " << mcerr * mcerr << " to " << (*newcov)(i, i));
	+ NUIS_LOG(FIT, "Adding cov stat " << mcerr * mcerr << " to "
	+ << (*newcov)(i, i));
	(newcov)(i, i) = oldval + mcerr mcerr;
	}

	// Reset the calc_cov to new invert
	delete calc_cov;
	calc_cov = GetInvert(newcov);

	// Delete the tempcov
	delete newcov;
	}

	calc_data->Scale(data_scale);
	calc_mc->Scale(data_scale);
	(calc_cov) = covar_scale;

	// iterate over bins in X (i,j)
	NUIS_LOG(DEB, "START Chi2 Calculation=================");
	for (int i = 0; i < calc_data->GetNbinsX(); i++) {
	double ibin_contrib = 0;
	- NUIS_LOG(DEB,
	- "[CHI2] i = " << i << " ["
	- << calc_data->GetXaxis()->GetBinLowEdge(i + 1) << " -- "
	- << calc_data->GetXaxis()->GetBinUpEdge(i + 1) << "].");
	+ NUIS_LOG(DEB, "[CHI2] i = "
	+ << i << " ["
	+ << calc_data->GetXaxis()->GetBinLowEdge(i + 1) << " -- "
	+ << calc_data->GetXaxis()->GetBinUpEdge(i + 1) << "].");
	for (int j = 0; j < calc_data->GetNbinsX(); j++) {
	NUIS_LOG(DEB, "[CHI2]\t j = "
	- << i << " [" << calc_data->GetXaxis()->GetBinLowEdge(j + 1)
	- << " -- " << calc_data->GetXaxis()->GetBinUpEdge(j + 1)
	- << "].");
	+ << i << " ["
	+ << calc_data->GetXaxis()->GetBinLowEdge(j + 1) << " -- "
	+ << calc_data->GetXaxis()->GetBinUpEdge(j + 1) << "].");
	if ((calc_data->GetBinContent(i + 1) != 0 \|\|
	calc_mc->GetBinContent(i + 1) != 0) &&
	((*calc_cov)(i, j) != 0)) {
	- NUIS_LOG(DEB,
	- "[CHI2]\t\t Chi2 contribution (i,j) = (" << i << "," << j << ")");
	+ NUIS_LOG(DEB, "[CHI2]\t\t Chi2 contribution (i,j) = (" << i << "," << j
	+ << ")");
	NUIS_LOG(DEB, "[CHI2]\t\t Data - MC(i) = "
	- << calc_data->GetBinContent(i + 1) << " - "
	- << calc_mc->GetBinContent(i + 1) << " = "
	- << (calc_data->GetBinContent(i + 1) -
	- calc_mc->GetBinContent(i + 1)));
	+ << calc_data->GetBinContent(i + 1) << " - "
	+ << calc_mc->GetBinContent(i + 1) << " = "
	+ << (calc_data->GetBinContent(i + 1) -
	+ calc_mc->GetBinContent(i + 1)));

	NUIS_LOG(DEB, "[CHI2]\t\t Data - MC(j) = "
	- << calc_data->GetBinContent(j + 1) << " - "
	- << calc_mc->GetBinContent(j + 1) << " = "
	- << (calc_data->GetBinContent(j + 1) -
	- calc_mc->GetBinContent(j + 1)));
	+ << calc_data->GetBinContent(j + 1) << " - "
	+ << calc_mc->GetBinContent(j + 1) << " = "
	+ << (calc_data->GetBinContent(j + 1) -
	+ calc_mc->GetBinContent(j + 1)));

	NUIS_LOG(DEB, "[CHI2]\t\t Covar = " << (*calc_cov)(i, j));

	- NUIS_LOG(DEB,
	- "[CHI2]\t\t Cont chi2 = " << ((calc_data->GetBinContent(i + 1) -
	- calc_mc->GetBinContent(i + 1)) *
	- (calc_cov)(i, j)
	- (calc_data->GetBinContent(j + 1) -
	- calc_mc->GetBinContent(j + 1)))
	- << " " << Chi2);
	+ NUIS_LOG(DEB, "[CHI2]\t\t Cont chi2 = "
	+ << ((calc_data->GetBinContent(i + 1) -
	+ calc_mc->GetBinContent(i + 1)) *
	+ (calc_cov)(i, j)
	+ (calc_data->GetBinContent(j + 1) -
	+ calc_mc->GetBinContent(j + 1)))
	+ << " " << Chi2);
	double bin_cont =
	((calc_data->GetBinContent(i + 1) - calc_mc->GetBinContent(i + 1)) *
	(calc_cov)(i, j)
	(calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1)));

	Chi2 += bin_cont;
	ibin_contrib += bin_cont;
	} else {
	NUIS_LOG(DEB, "Skipping chi2 contribution (i,j) = ("
	- << i << "," << j
	- << "), Data = " << calc_data->GetBinContent(i + 1)
	- << ", MC = " << calc_mc->GetBinContent(i + 1)
	- << ", Cov = " << (*calc_cov)(i, j));
	+ << i << "," << j
	+ << "), Data = " << calc_data->GetBinContent(i + 1)
	+ << ", MC = " << calc_mc->GetBinContent(i + 1)
	+ << ", Cov = " << (*calc_cov)(i, j));
	Chi2 += 0.;
	}
	}
	if (outchi2perbin) {
	outchi2perbin->SetBinContent(i + 1, ibin_contrib);
	}
	}

	// Cleanup
	delete calc_cov;
	delete calc_data;
	delete calc_mc;

	return Chi2;
	}

	//*******************************************************************
	Double_t StatUtils::GetChi2FromCov(TH2D data, TH2D mc, TMatrixDSym *invcov,
	TH2I map, TH2I mask, TH2D *outchi2perbin) {
	//*******************************************************************

	// Generate a simple map
	if (!map) {
	map = StatUtils::GenerateMap(data);
	}

	// Convert to 1D Histograms
	TH1D *data_1D = MapToTH1D(data, map);
	TH1D *mc_1D = MapToTH1D(mc, map);
	TH1I *mask_1D = MapToMask(mask, map);
	TH1D *outchi2perbin_1D = NULL;
	TH1D *outchi2perbin_map_1D = NULL;

	if (outchi2perbin) {
	outchi2perbin_1D = MapToTH1D(outchi2perbin, map);
	for (Int_t xbi_it = 0; xbi_it < outchi2perbin->GetXaxis()->GetNbins();
	++xbi_it) {
	for (Int_t ybi_it = 0; ybi_it < outchi2perbin->GetYaxis()->GetNbins();
	++ybi_it) {
	int gbin = outchi2perbin->GetBin(xbi_it + 1, ybi_it + 1);
	// std::cout << " gbin " << gbin << " corresponds to "
	// << " x: " << (xbi_it + 1) << ", y: " << (ybi_it + 1)
	// << std::endl;
	outchi2perbin->SetBinContent(xbi_it + 1, ybi_it + 1, gbin);
	}
	}
	outchi2perbin_map_1D = MapToTH1D(outchi2perbin, map);
	}

	// Calculate 1D chi2 from 1D Plots
	Double_t Chi2 = StatUtils::GetChi2FromCov(data_1D, mc_1D, invcov, mask_1D, 1,
	1E76, outchi2perbin_1D);
	if (outchi2perbin) {
	for (int xbi_it = 0; xbi_it < outchi2perbin_1D->GetXaxis()->GetNbins();
	++xbi_it) {
	// std::cout << " adding chi2 "
	// << outchi2perbin_1D->GetBinContent(xbi_it + 1)
	// << " from 1d bin " << (xbi_it + 1) << " to gbin "
	- // << outchi2perbin_map_1D->GetBinContent(xbi_it + 1) << std::endl;
	+ // << outchi2perbin_map_1D->GetBinContent(xbi_it + 1) <<
	+ // std::endl;
	outchi2perbin->SetBinContent(
	outchi2perbin_map_1D->GetBinContent(xbi_it + 1),
	outchi2perbin_1D->GetBinContent(xbi_it + 1));
	}
	}

	// CleanUp
	delete data_1D;
	delete mc_1D;
	delete mask_1D;
	delete outchi2perbin_1D;
	delete outchi2perbin_map_1D;

	return Chi2;
	}

	//*******************************************************************
	Double_t StatUtils::GetChi2FromSVD(TH1D data, TH1D mc, TMatrixDSym *cov,
	TH1I *mask) {
	//*******************************************************************

	Double_t Chi2 = 0.0;
	TMatrixDSym calc_cov = (TMatrixDSym )cov->Clone();
	TH1D calc_data = (TH1D )data->Clone();
	TH1D calc_mc = (TH1D )mc->Clone();

	// If a mask if applied we need to apply it before the matrix is inverted
	if (mask) {
	calc_cov = StatUtils::ApplyMatrixMasking(cov, mask);
	calc_data = StatUtils::ApplyHistogramMasking(data, mask);
	calc_mc = StatUtils::ApplyHistogramMasking(mc, mask);
	}

	// Decompose matrix
	TDecompSVD LU = TDecompSVD((*calc_cov));
	LU.Decompose();
	TMatrixDSym *cov_U =
	new TMatrixDSym(calc_data->GetNbinsX(), LU.GetU().GetMatrixArray(), "");
	TVectorD *cov_S = new TVectorD(LU.GetSig());

	// Apply basis rotation before adding up chi2
	Double_t rotated_difference = 0.0;

	for (int i = 0; i < calc_data->GetNbinsX(); i++) {
	rotated_difference = 0.0;

	// Rotate basis of Data - MC
	for (int j = 0; j < calc_data->GetNbinsY(); j++)
	rotated_difference +=
	(calc_data->GetBinContent(j + 1) - calc_mc->GetBinContent(j + 1)) *
	(*cov_U)(j, i);

	// Divide by rotated error cov_S
	Chi2 += rotated_difference * rotated_difference * 1E76 / (*cov_S)(i);
	}

	// Cleanup
	delete calc_cov;
	delete calc_data;
	delete calc_mc;
	delete cov_U;
	delete cov_S;

	return Chi2;
	}

	//*******************************************************************
	Double_t StatUtils::GetChi2FromSVD(TH2D data, TH2D mc, TMatrixDSym *cov,
	TH2I map, TH2I mask) {
	//*******************************************************************

	// Generate a simple map
	if (!map)
	map = StatUtils::GenerateMap(data);

	// Convert to 1D Histograms
	TH1D *data_1D = MapToTH1D(data, map);
	TH1D *mc_1D = MapToTH1D(mc, map);
	TH1I *mask_1D = MapToMask(mask, map);

	// Calculate from 1D
	Double_t Chi2 = StatUtils::GetChi2FromSVD(data_1D, mc_1D, cov, mask_1D);

	// CleanUp
	delete data_1D;
	delete mc_1D;
	delete mask_1D;

	return Chi2;
	}

	//*******************************************************************
	double StatUtils::GetChi2FromEventRate(TH1D data, TH1D mc, TH1I *mask) {
	//*******************************************************************

	// If just an event rate, for chi2 just use Poission Likelihood to calculate
	// the chi2 component
	double chi2 = 0.0;
	TH1D calc_data = (TH1D )data->Clone();
	TH1D calc_mc = (TH1D )mc->Clone();

	// Apply masking if required
	if (mask) {
	calc_data = ApplyHistogramMasking(data, mask);
	calc_mc = ApplyHistogramMasking(mc, mask);
	}

	// Iterate over bins in X
	for (int i = 0; i < calc_data->GetNbinsX(); i++) {
	double dt = calc_data->GetBinContent(i + 1);
	double mc = calc_mc->GetBinContent(i + 1);

	if (mc <= 0)
	continue;

	if (dt <= 0) {
	// Only add difference
	chi2 += 2 * (mc - dt);
	} else {
	// Do the chi2 for Poisson distributions
	chi2 += 2 * (mc - dt + (dt * log(dt / mc)));
	}

	/*
	LOG(REC)<<"Evt Chi2 cont = "<<i<<" "
	<<mc<<" "<<dt<<" "
	<<2 * (mc - dt + (dt+0.) * log((dt+0.) / (mc+0.)))
	<<" "<<Chi2<<std::endl;
	*/
	}

	// cleanup
	delete calc_data;
	delete calc_mc;

	return chi2;
	}

	//*******************************************************************
	Double_t StatUtils::GetChi2FromEventRate(TH2D data, TH2D mc, TH2I *map,
	TH2I *mask) {
	//*******************************************************************

	// Generate a simple map
	if (!map)
	map = StatUtils::GenerateMap(data);

	// Convert to 1D Histograms
	TH1D *data_1D = MapToTH1D(data, map);
	TH1D *mc_1D = MapToTH1D(mc, map);
	TH1I *mask_1D = MapToMask(mask, map);

	// Calculate from 1D
	Double_t Chi2 = StatUtils::GetChi2FromEventRate(data_1D, mc_1D, mask_1D);

	// CleanUp
	delete data_1D;
	delete mc_1D;
	delete mask_1D;

	return Chi2;
	}

	//*******************************************************************
	Double_t StatUtils::GetLikelihoodFromDiag(TH1D data, TH1D mc, TH1I *mask) {
	//*******************************************************************
	// Currently just a placeholder!
	(void)data;
	(void)mc;
	(void)mask;
	return 0.0;
	};

	//*******************************************************************
	Double_t StatUtils::GetLikelihoodFromDiag(TH2D data, TH2D mc, TH2I *map,
	TH2I *mask) {
	//*******************************************************************

	// Generate a simple map
	if (!map)
	map = StatUtils::GenerateMap(data);

	// Convert to 1D Histograms
	TH1D *data_1D = MapToTH1D(data, map);
	TH1D *mc_1D = MapToTH1D(mc, map);
	TH1I *mask_1D = MapToMask(mask, map);

	// Calculate from 1D
	Double_t MLE = StatUtils::GetLikelihoodFromDiag(data_1D, mc_1D, mask_1D);

	// CleanUp
	delete data_1D;
	delete mc_1D;
	delete mask_1D;

	return MLE;
	};

	//*******************************************************************
	Double_t StatUtils::GetLikelihoodFromCov(TH1D data, TH1D mc,
	TMatrixDSym invcov, TH1I mask) {
	//*******************************************************************
	// Currently just a placeholder !
	(void)data;
	(void)mc;
	(void)invcov;
	(void)mask;

	return 0.0;
	};

	//*******************************************************************
	Double_t StatUtils::GetLikelihoodFromCov(TH2D data, TH2D mc,
	TMatrixDSym invcov, TH2I map,
	TH2I *mask) {
	//*******************************************************************

	// Generate a simple map
	if (!map)
	map = StatUtils::GenerateMap(data);

	// Convert to 1D Histograms
	TH1D *data_1D = MapToTH1D(data, map);
	TH1D *mc_1D = MapToTH1D(mc, map);
	TH1I *mask_1D = MapToMask(mask, map);

	// Calculate from 1D
	Double_t MLE =
	StatUtils::GetLikelihoodFromCov(data_1D, mc_1D, invcov, mask_1D);

	// CleanUp
	delete data_1D;
	delete mc_1D;
	delete mask_1D;

	return MLE;
	};

	//*******************************************************************
	Double_t StatUtils::GetLikelihoodFromSVD(TH1D data, TH1D mc, TMatrixDSym *cov,
	TH1I *mask) {
	//*******************************************************************
	// Currently just a placeholder!
	(void)data;
	(void)mc;
	(void)cov;
	(void)mask;

	return 0.0;
	};

	//*******************************************************************
	Double_t StatUtils::GetLikelihoodFromSVD(TH2D data, TH2D mc, TMatrixDSym *cov,
	TH2I map, TH2I mask) {
	//*******************************************************************

	// Generate a simple map
	if (!map)
	map = StatUtils::GenerateMap(data);

	// Convert to 1D Histograms
	TH1D *data_1D = MapToTH1D(data, map);
	TH1D *mc_1D = MapToTH1D(mc, map);
	TH1I *mask_1D = MapToMask(mask, map);

	// Calculate from 1D
	Double_t MLE = StatUtils::GetLikelihoodFromSVD(data_1D, mc_1D, cov, mask_1D);

	// CleanUp
	delete data_1D;
	delete mc_1D;
	delete mask_1D;

	return MLE;
	};

	//*******************************************************************
	Double_t StatUtils::GetLikelihoodFromEventRate(TH1D data, TH1D mc,
	TH1I *mask) {
	//*******************************************************************
	// Currently just a placeholder!
	(void)data;
	(void)mc;
	(void)mask;

	return 0.0;
	};

	//*******************************************************************
	Double_t StatUtils::GetLikelihoodFromEventRate(TH2D data, TH2D mc, TH2I *map,
	TH2I *mask) {
	//*******************************************************************

	// Generate a simple map
	if (!map)
	map = StatUtils::GenerateMap(data);

	// Convert to 1D Histograms
	TH1D *data_1D = MapToTH1D(data, map);
	TH1D *mc_1D = MapToTH1D(mc, map);
	TH1I *mask_1D = MapToMask(mask, map);

	// Calculate from 1D
	Double_t MLE = StatUtils::GetChi2FromEventRate(data_1D, mc_1D, mask_1D);

	// CleanUp
	delete data_1D;
	delete mc_1D;
	delete mask_1D;

	return MLE;
	};

	//*******************************************************************
	Int_t StatUtils::GetNDOF(TH1D hist, TH1I mask) {
	//*******************************************************************

	TH1D calc_hist = (TH1D )hist->Clone();

	// If a mask is provided we need to apply it before getting NDOF
	if (mask) {
	calc_hist = StatUtils::ApplyHistogramMasking(hist, mask);
	}

	// NDOF is defined as total number of bins with non-zero errors
	Int_t NDOF = 0;
	for (int i = 0; i < calc_hist->GetNbinsX(); i++) {
	if (calc_hist->GetBinError(i + 1) > 0.0)
	NDOF++;
	}

	delete calc_hist;

	return NDOF;
	};

	//*******************************************************************
	Int_t StatUtils::GetNDOF(TH2D hist, TH2I map, TH2I *mask) {
	//*******************************************************************

	Int_t NDOF = 0;
	if (!map)
	map = StatUtils::GenerateMap(hist);

	for (int i = 0; i < hist->GetNbinsX(); i++) {
	for (int j = 0; j < hist->GetNbinsY(); j++) {
	if (mask->GetBinContent(i + 1, j + 1))
	continue;
	if (map->GetBinContent(i + 1, j + 1) <= 0)
	continue;

	NDOF++;
	}
	}

	return NDOF;
	};

	//*******************************************************************
	TH1D StatUtils::ThrowHistogram(TH1D hist, TMatrixDSym *cov, bool throwdiag,
	TH1I *mask) {
	//*******************************************************************

	TH1D *calc_hist =
	(TH1D *)hist->Clone((std::string(hist->GetName()) + "_THROW").c_str());
	TMatrixDSym calc_cov = (TMatrixDSym )cov->Clone();
	Double_t correl_val = 0.0;

	// If a mask if applied we need to apply it before the matrix is decomposed
	if (mask) {
	calc_cov = ApplyMatrixMasking(cov, mask);
	calc_hist = ApplyHistogramMasking(calc_hist, mask);
	}

	// If a covariance is provided we need a preset random vector and a decomp
	std::vector<Double_t> rand_val;
	TMatrixDSym *decomp_cov = NULL;

	if (cov) {
	for (int i = 0; i < hist->GetNbinsX(); i++) {
	rand_val.push_back(gRandom->Gaus(0.0, 1.0));
	}

	// Decomp the matrix
	decomp_cov = StatUtils::GetDecomp(calc_cov);
	}

	// iterate over bins
	for (int i = 0; i < hist->GetNbinsX(); i++) {
	// By Default the errors on the histogram are thrown uncorrelated to the
	// other errors
	/*
	if (throwdiag) {
	calc_hist->SetBinContent(i + 1, (calc_hist->GetBinContent(i + 1) + \
	gRandom->Gaus(0.0, 1.0) * calc_hist->GetBinError(i + 1)) );
	}
	*/

	// If a covariance is provided that is also thrown
	if (cov) {
	correl_val = 0.0;

	for (int j = 0; j < hist->GetNbinsX(); j++) {
	correl_val += rand_val[j] * (*decomp_cov)(j, i);
	}
	calc_hist->SetBinContent(
	i + 1, (calc_hist->GetBinContent(i + 1) + correl_val * 1E-38));
	}
	}

	delete calc_cov;
	delete decomp_cov;

	// return this new thrown data
	return calc_hist;
	};

	//*******************************************************************
	TH2D StatUtils::ThrowHistogram(TH2D hist, TMatrixDSym cov, TH2I map,
	bool throwdiag, TH2I *mask) {
	//*******************************************************************

	// PLACEHOLDER!!!!!!!!!
	// Currently no support for throwing 2D Histograms from a covariance
	(void)hist;
	(void)cov;
	(void)map;
	(void)throwdiag;
	(void)mask;

	// /todo
	// Sort maps if required
	// Throw the covariance for a 1D plot
	// Unmap back to 2D Histogram

	return hist;
	}

	//*******************************************************************
	TH1D StatUtils::ApplyHistogramMasking(TH1D hist, TH1I *mask) {
	//*******************************************************************

	if (!mask)
	return ((TH1D *)hist->Clone());

	// This masking is only sufficient for chi2 calculations, and will have dodgy
	// bin edges.

	// Get New Bin Count
	Int_t NBins = 0;
	for (int i = 0; i < hist->GetNbinsX(); i++) {
	if (mask->GetBinContent(i + 1))
	continue;
	NBins++;
	}

	// Make new hist
	std::string newmaskname = std::string(hist->GetName()) + "_MSKD";
	TH1D *calc_hist =
	new TH1D(newmaskname.c_str(), newmaskname.c_str(), NBins, 0, NBins);

	// fill new hist
	int binindex = 0;
	for (int i = 0; i < hist->GetNbinsX(); i++) {
	if (mask->GetBinContent(i + 1)) {
	- NUIS_LOG(REC, "Applying mask to bin " << i + 1 << " " << hist->GetName());
	+ NUIS_LOG(DEB, "Applying mask to bin " << i + 1 << " " << hist->GetName());
	continue;
	}
	calc_hist->SetBinContent(binindex + 1, hist->GetBinContent(i + 1));
	calc_hist->SetBinError(binindex + 1, hist->GetBinError(i + 1));
	binindex++;
	}

	return calc_hist;
	};

	//*******************************************************************
	TH2D StatUtils::ApplyHistogramMasking(TH2D hist, TH2I *mask) {
	//*******************************************************************

	TH2D newhist = (TH2D )hist->Clone();

	if (!mask)
	return newhist;

	for (int i = 0; i < hist->GetNbinsX(); i++) {
	for (int j = 0; j < hist->GetNbinsY(); j++) {
	if (mask->GetBinContent(i + 1, j + 1) > 0) {
	newhist->SetBinContent(i + 1, j + 1, 0.0);
	newhist->SetBinContent(i + 1, j + 1, 0.0);
	}
	}
	}

	return newhist;
	}

	//*******************************************************************
	TMatrixDSym StatUtils::ApplyMatrixMasking(TMatrixDSym mat, TH1I *mask) {
	//*******************************************************************

	if (!mask)
	return (TMatrixDSym *)(mat->Clone());

	// Get New Bin Count
	Int_t NBins = 0;
	for (int i = 0; i < mask->GetNbinsX(); i++) {
	if (mask->GetBinContent(i + 1))
	continue;
	NBins++;
	}

	// make new matrix
	TMatrixDSym *calc_mat = new TMatrixDSym(NBins);
	int col, row;

	// Need to mask out bins in the current matrix
	row = 0;
	for (int i = 0; i < mask->GetNbinsX(); i++) {
	col = 0;

	// skip if masked
	if (mask->GetBinContent(i + 1) > 0.5)
	continue;

	for (int j = 0; j < mask->GetNbinsX(); j++) {
	// skip if masked
	if (mask->GetBinContent(j + 1) > 0.5)
	continue;

	(calc_mat)(row, col) = (mat)(i, j);
	col++;
	}
	row++;
	}

	return calc_mat;
	};

	//*******************************************************************
	TMatrixDSym StatUtils::ApplyMatrixMasking(TMatrixDSym mat, TH2D *data,
	TH2I mask, TH2I map) {
	//*******************************************************************

	if (!map)
	map = StatUtils::GenerateMap(data);
	TH1I *mask_1D = StatUtils::MapToMask(mask, map);

	TMatrixDSym *newmat = StatUtils::ApplyMatrixMasking(mat, mask_1D);

	delete mask_1D;
	return newmat;
	}

	//*******************************************************************
	TMatrixDSym StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym mat,
	TH1I *mask) {
	//*******************************************************************

	TMatrixDSym *new_mat = GetInvert(mat);
	TMatrixDSym *masked_mat = ApplyMatrixMasking(new_mat, mask);

	TMatrixDSym *inverted_mat = GetInvert(masked_mat);

	delete masked_mat;
	delete new_mat;

	return inverted_mat;
	};

	//*******************************************************************
	TMatrixDSym StatUtils::ApplyInvertedMatrixMasking(TMatrixDSym mat, TH2D *data,
	TH2I mask, TH2I map) {
	//*******************************************************************

	if (!map)
	map = StatUtils::GenerateMap(data);
	TH1I *mask_1D = StatUtils::MapToMask(mask, map);

	TMatrixDSym *newmat = ApplyInvertedMatrixMasking(mat, mask_1D);

	delete mask_1D;
	return newmat;
	}

	//*******************************************************************
	TMatrixDSym StatUtils::GetInvert(TMatrixDSym mat) {
	//*******************************************************************

	TMatrixDSym new_mat = (TMatrixDSym )mat->Clone();

	// Check for diagonal
	bool non_diagonal = false;
	for (int i = 0; i < new_mat->GetNrows(); i++) {
	for (int j = 0; j < new_mat->GetNrows(); j++) {
	if (i == j)
	continue;

	if ((*new_mat)(i, j) != 0.0) {
	non_diagonal = true;
	break;
	}
	}
	}

	// If diag, just flip the diag
	if (!non_diagonal or new_mat->GetNrows() == 1) {
	for (int i = 0; i < new_mat->GetNrows(); i++) {
	if ((*new_mat)(i, i) != 0.0)
	(new_mat)(i, i) = 1.0 / (new_mat)(i, i);
	else
	(*new_mat)(i, i) = 0.0;
	}
	return new_mat;
	}

	// Invert full matrix
	TDecompSVD LU = TDecompSVD((*new_mat));
	new_mat =
	new TMatrixDSym(new_mat->GetNrows(), LU.Invert().GetMatrixArray(), "");

	return new_mat;
	}

	//*******************************************************************
	TMatrixDSym StatUtils::GetDecomp(TMatrixDSym mat) {
	//*******************************************************************

	TMatrixDSym new_mat = (TMatrixDSym )mat->Clone();
	int nrows = new_mat->GetNrows();

	// Check for diagonal
	bool diagonal = true;
	for (int i = 0; i < nrows; i++) {
	for (int j = 0; j < nrows; j++) {
	if (i == j)
	continue;

	if ((*new_mat)(i, j) != 0.0) {
	diagonal = false;
	break;
	}
	}
	}

	// If diag, just flip the diag
	if (diagonal or nrows == 1) {
	for (int i = 0; i < nrows; i++) {
	if ((*new_mat)(i, i) > 0.0)
	(new_mat)(i, i) = sqrt((new_mat)(i, i));
	else
	(*new_mat)(i, i) = 0.0;
	}
	return new_mat;
	}

	TDecompChol LU = TDecompChol(*new_mat);
	LU.Decompose();
	delete new_mat;

	TMatrixDSym *dec_mat = new TMatrixDSym(nrows, LU.GetU().GetMatrixArray(), "");

	return dec_mat;
	}

	//*******************************************************************
	void StatUtils::ForceNormIntoCovar(TMatrixDSym &mat, TH1D hist, double norm) {
	//*******************************************************************

	if (!mat)
	mat = MakeDiagonalCovarMatrix(hist);

	int nbins = mat->GetNrows();
	TMatrixDSym *new_mat = new TMatrixDSym(nbins);

	for (int i = 0; i < nbins; i++) {
	for (int j = 0; j < nbins; j++) {
	double valx = hist->GetBinContent(i + 1) * 1E38;
	double valy = hist->GetBinContent(j + 1) * 1E38;

	(new_mat)(i, j) = (mat)(i, j) + norm * norm * valx * valy;
	}
	}

	// Swap the two
	delete mat;
	mat = new_mat;

	return;
	};

	//*******************************************************************
	void StatUtils::ForceNormIntoCovar(TMatrixDSym mat, TH2D data, double norm,
	TH2I *map) {
	//*******************************************************************

	if (!map)
	map = StatUtils::GenerateMap(data);
	TH1D *data_1D = MapToTH1D(data, map);

	StatUtils::ForceNormIntoCovar(mat, data_1D, norm);
	delete data_1D;

	return;
	}

	//*******************************************************************
	TMatrixDSym StatUtils::MakeDiagonalCovarMatrix(TH1D data, double scaleF) {
	//*******************************************************************

	TMatrixDSym *newmat = new TMatrixDSym(data->GetNbinsX());

	for (int i = 0; i < data->GetNbinsX(); i++) {
	(*newmat)(i, i) =
	data->GetBinError(i + 1) * data->GetBinError(i + 1) * scaleF * scaleF;
	}

	return newmat;
	}

	//*******************************************************************
	TMatrixDSym StatUtils::MakeDiagonalCovarMatrix(TH2D data, TH2I *map,
	double scaleF) {
	//*******************************************************************

	if (!map)
	map = StatUtils::GenerateMap(data);
	TH1D *data_1D = MapToTH1D(data, map);

	return StatUtils::MakeDiagonalCovarMatrix(data_1D, scaleF);
	};

	//*******************************************************************
	void StatUtils::SetDataErrorFromCov(TH1D DataHist, TMatrixDSym cov,
	double scale, bool ErrorCheck) {
	//*******************************************************************

	// Check
	if (ErrorCheck) {
	if (cov->GetNrows() != DataHist->GetNbinsX()) {
	NUIS_ERR(
	FTL,
	"Nrows in cov don't match nbins in DataHist for SetDataErrorFromCov");
	NUIS_ERR(FTL, "Nrows = " << cov->GetNrows());
	NUIS_ABORT("Nbins = " << DataHist->GetNbinsX());
	}
	}

	// Set bin errors form cov diag
	// Check if the errors are set
	bool ErrorsSet = false;
	for (int i = 0; i < DataHist->GetNbinsX(); i++) {
	if (ErrorsSet == true)
	break;
	if (DataHist->GetBinError(i + 1) != 0 && DataHist->GetBinContent(i + 1) > 0)
	ErrorsSet = true;
	}

	// Now loop over
	if (ErrorsSet && ErrorCheck) {
	for (int i = 0; i < DataHist->GetNbinsX(); i++) {
	double DataHisterr = DataHist->GetBinError(i + 1);
	double coverr = sqrt((cov)(i, i)) scale;
	// Check that the errors are within 1% of eachother
	if (fabs(DataHisterr - coverr) / DataHisterr > 0.01) {
	NUIS_ERR(WRN, "Data error does not match covariance error for bin "
	- << i + 1 << " ("
	- << DataHist->GetXaxis()->GetBinLowEdge(i + 1) << "-"
	- << DataHist->GetXaxis()->GetBinLowEdge(i + 2) << ")");
	+ << i + 1 << " ("
	+ << DataHist->GetXaxis()->GetBinLowEdge(i + 1) << "-"
	+ << DataHist->GetXaxis()->GetBinLowEdge(i + 2) << ")");
	NUIS_ERR(WRN, "Data error: " << DataHisterr);
	NUIS_ERR(WRN, "Cov error: " << coverr);
	}
	}
	// Else blindly trust the covariance
	} else {
	for (int i = 0; i < DataHist->GetNbinsX(); i++) {
	DataHist->SetBinError(i + 1, sqrt((cov)(i, i)) scale);
	}
	}

	return;
	}

	//*******************************************************************
	void StatUtils::SetDataErrorFromCov(TH2D data, TMatrixDSym cov, TH2I *map,
	double scale, bool ErrorCheck) {
	//*******************************************************************

	// Check
	if (ErrorCheck) {
	if (cov->GetNrows() != data->GetNbinsX() * data->GetNbinsY()) {
	- NUIS_ERR(
	- FTL,
	- "Nrows in cov don't match nbins in data for SetDataNUIS_ERRorFromCov");
	+ NUIS_ERR(FTL, "Nrows in cov don't match nbins in data for "
	+ "SetDataNUIS_ERRorFromCov");
	NUIS_ERR(FTL, "Nrows = " << cov->GetNrows());
	NUIS_ABORT("Nbins = " << data->GetNbinsX());
	}
	}

	// Set bin errors form cov diag
	// Check if the errors are set
	bool ErrorsSet = false;
	for (int i = 0; i < data->GetNbinsX(); i++) {
	for (int j = 0; j < data->GetNbinsX(); j++) {
	if (ErrorsSet == true)
	break;
	if (data->GetBinError(i + 1, j + 1) != 0)
	ErrorsSet = true;
	}
	}

	// Create map if required
	if (!map)
	map = StatUtils::GenerateMap(data);

	// Set Bin Errors from cov diag
	int count = 0;
	for (int i = 0; i < data->GetNbinsX(); i++) {
	for (int j = 0; j < data->GetNbinsY(); j++) {
	if (data->GetBinContent(i + 1, j + 1) == 0.0)
	continue;
	// If we have errors on our histogram the map is good
	count = map->GetBinContent(i + 1, j + 1) - 1;
	double dataerr = data->GetBinError(i + 1, j + 1);
	double coverr = sqrt((cov)(count, count)) scale;
	// Check that the errors are within 1% of eachother
	if (ErrorsSet && ErrorCheck) {
	if (fabs(dataerr - coverr) / dataerr > 0.01) {
	NUIS_ERR(WRN, "Data error does not match covariance error for bin "
	- << i + 1 << " ("
	- << data->GetXaxis()->GetBinLowEdge(i + 1) << "-"
	- << data->GetXaxis()->GetBinLowEdge(i + 2) << ")");
	+ << i + 1 << " ("
	+ << data->GetXaxis()->GetBinLowEdge(i + 1) << "-"
	+ << data->GetXaxis()->GetBinLowEdge(i + 2) << ")");
	NUIS_ERR(WRN, "Data error: " << dataerr);
	NUIS_ERR(WRN, "Cov error: " << coverr);
	}
	} else {
	data->SetBinError(i + 1, j + 1, sqrt((cov)(count, count)) scale);
	}
	}
	}
	// Delete the map now that we don't need it
	// Woops, it's needed elsewhere there! (Grrrrr)
	// map->Delete();
	}

	TMatrixDSym StatUtils::ExtractShapeOnlyCovar(TMatrixDSym full_covar,
	TH1 *data_hist,
	double data_scale) {
	int nbins = full_covar->GetNrows();
	TMatrixDSym *shape_covar = new TMatrixDSym(nbins);

	// Check nobody is being silly
	if (data_hist->GetNbinsX() != nbins) {
	NUIS_ERR(WRN, "Inconsistent matrix and data histogram passed to "
	- "StatUtils::ExtractShapeOnlyCovar!");
	+ "StatUtils::ExtractShapeOnlyCovar!");
	NUIS_ERR(WRN, "data_hist has " << data_hist->GetNbinsX() << " matrix has "
	- << nbins);
	+ << nbins);
	int err_bins = data_hist->GetNbinsX();
	if (nbins > err_bins)
	err_bins = nbins;
	for (int i = 0; i < err_bins; ++i) {
	NUIS_ERR(WRN, "Matrix diag. = " << (*full_covar)(i, i) << " data = "
	- << data_hist->GetBinContent(i + 1));
	+ << data_hist->GetBinContent(i + 1));
	}
	return NULL;
	}

	double total_data = 0;
	double total_covar = 0;

	// Initial loop to calculate some constants
	for (int i = 0; i < nbins; ++i) {
	total_data += data_hist->GetBinContent(i + 1) * data_scale;
	for (int j = 0; j < nbins; ++j) {
	total_covar += (*full_covar)(i, j);
	}
	}

	if (total_data == 0 \|\| total_covar == 0) {
	NUIS_ERR(WRN, "Stupid matrix or data histogram passed to "
	- "StatUtils::ExtractShapeOnlyCovar! Ignoring...");
	+ "StatUtils::ExtractShapeOnlyCovar! Ignoring...");
	return NULL;
	}

	NUIS_LOG(SAM, "Norm error = " << sqrt(total_covar) / total_data);

	// Now loop over and calculate the shape-only matrix
	for (int i = 0; i < nbins; ++i) {
	double data_i = data_hist->GetBinContent(i + 1) * data_scale;

	for (int j = 0; j < nbins; ++j) {
	double data_j = data_hist->GetBinContent(j + 1) * data_scale;

	double norm_term =
	data_i * data_j * total_covar / total_data / total_data;
	double mix_sum1 = 0;
	double mix_sum2 = 0;

	for (int k = 0; k < nbins; ++k) {
	mix_sum1 += (*full_covar)(k, j);
	mix_sum2 += (*full_covar)(i, k);
	}

	double mix_term1 =
	data_i * (mix_sum1 / total_data -
	total_covar * data_j / total_data / total_data);
	double mix_term2 =
	data_j * (mix_sum2 / total_data -
	total_covar * data_i / total_data / total_data);

	(*shape_covar)(i, j) =
	(*full_covar)(i, j) - mix_term1 - mix_term2 - norm_term;
	}
	}
	return shape_covar;
	}

	//*******************************************************************
	TH2I StatUtils::GenerateMap(TH2D hist) {
	//*******************************************************************

	std::string maptitle = std::string(hist->GetName()) + "_MAP";

	TH2I *map =
	new TH2I(maptitle.c_str(), maptitle.c_str(), hist->GetNbinsX(), 0,
	hist->GetNbinsX(), hist->GetNbinsY(), 0, hist->GetNbinsY());

	Int_t index = 1;

	for (int i = 0; i < hist->GetNbinsX(); i++) {
	for (int j = 0; j < hist->GetNbinsY(); j++) {
	if (hist->GetBinContent(i + 1, j + 1) > 0) {
	map->SetBinContent(i + 1, j + 1, index);
	index++;
	} else {
	map->SetBinContent(i + 1, j + 1, 0);
	}
	}
	}

	return map;
	}

	//*******************************************************************
	TH1D StatUtils::MapToTH1D(TH2D hist, TH2I *map) {
	//*******************************************************************

	if (!hist)
	return NULL;

	// Get N bins for 1D plot
	Int_t Nbins = map->GetMaximum();

	std::string name1D = std::string(hist->GetName()) + "_1D";

	// Make new 1D Hist
	TH1D *newhist = new TH1D(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins);

	// map bin contents
	for (int i = 0; i < map->GetNbinsX(); i++) {
	for (int j = 0; j < map->GetNbinsY(); j++) {
	if (map->GetBinContent(i + 1, j + 1) == 0)
	continue;
	newhist->SetBinContent(map->GetBinContent(i + 1, j + 1),
	hist->GetBinContent(i + 1, j + 1));
	newhist->SetBinError(map->GetBinContent(i + 1, j + 1),
	hist->GetBinError(i + 1, j + 1));
	}
	}

	// return
	return newhist;
	}

	//*******************************************************************
	TH1I StatUtils::MapToMask(TH2I hist, TH2I *map) {
	//*******************************************************************

	TH1I *newhist = NULL;
	if (!hist)
	return newhist;

	// Get N bins for 1D plot
	Int_t Nbins = map->GetMaximum();
	std::string name1D = std::string(hist->GetName()) + "_1D";

	// Make new 1D Hist
	newhist = new TH1I(name1D.c_str(), name1D.c_str(), Nbins, 0, Nbins);

	// map bin contents
	for (int i = 0; i < map->GetNbinsX(); i++) {
	for (int j = 0; j < map->GetNbinsY(); j++) {
	if (map->GetBinContent(i + 1, j + 1) == 0)
	continue;

	newhist->SetBinContent(map->GetBinContent(i + 1, j + 1),
	hist->GetBinContent(i + 1, j + 1));
	}
	}

	// return
	return newhist;
	}

	TMatrixDSym StatUtils::GetCovarFromCorrel(TMatrixDSym correl, TH1D *data) {
	int nbins = correl->GetNrows();
	TMatrixDSym *covar = new TMatrixDSym(nbins);

	for (int i = 0; i < nbins; i++) {
	for (int j = 0; j < nbins; j++) {
	(*covar)(i, j) =
	(correl)(i, j) data->GetBinError(i + 1) * data->GetBinError(j + 1);
	}
	}

	return covar;
	}

	//*******************************************************************
	TMatrixD *StatUtils::GetMatrixFromTextFile(std::string covfile, int dimx,
	int dimy) {
	//*******************************************************************

	// Determine dim
	if (dimx == -1 and dimy == -1) {
	std::string line;
	std::ifstream covar(covfile.c_str(), std::ifstream::in);

	int row = 0;
	while (std::getline(covar >> std::ws, line, '\n')) {
	int column = 0;

	std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");

	if (entries.size() <= 1) {
	NUIS_ERR(WRN, "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 "
	- "entries on this line: "
	- << row);
	+ "entries on this line: "
	+ << row);
	}
	for (std::vector<double>::iterator iter = entries.begin();
	iter != entries.end(); iter++) {
	column++;

	if (column > dimx)
	dimx = column;
	}
	row++;
	if (row > dimy)
	dimy = row;
	}
	}

	// Or assume symmetric
	if (dimx != -1 and dimy == -1) {
	dimy = dimx;
	}
	assert(dimy != -1 && " matrix dimy not set.");

	// Make new matrix
	TMatrixD *mat = new TMatrixD(dimx, dimy);
	std::string line;
	std::ifstream covar(covfile.c_str(), std::ifstream::in);

	int row = 0;
	while (std::getline(covar >> std::ws, line, '\n')) {
	int column = 0;

	std::vector<double> entries = GeneralUtils::ParseToDbl(line, " ");
	if (entries.size() <= 1) {
	NUIS_ERR(WRN, "StatUtils::GetMatrixFromTextFile, matrix only has <= 1 "
	- "entries on this line: "
	- << row);
	+ "entries on this line: "
	+ << row);
	}
	for (std::vector<double>::iterator iter = entries.begin();
	iter != entries.end(); iter++) {
	// Check Rows
	// assert(row > mat->GetNrows() && " covar rows doesn't match matrix
	// rows.");
	// assert(column > mat->GetNcols() && " covar cols doesn't match matrix
	// cols.");

	// Fill Matrix
	(mat)(row, column) = (iter);
	column++;
	}
	row++;
	}

	return mat;
	}

	//*******************************************************************
	TMatrixD *StatUtils::GetMatrixFromRootFile(std::string covfile,
	std::string histname) {
	//*******************************************************************

	std::string inputfile = covfile + ";" + histname;
	std::vector<std::string> splitfile = GeneralUtils::ParseToStr(inputfile, ";");

	if (splitfile.size() < 2) {
	NUIS_ABORT("No object name given!");
	}

	// Get file
	TFile *tempfile = new TFile(splitfile[0].c_str(), "READ");

	// Get Object
	TObject *obj = tempfile->Get(splitfile[1].c_str());
	if (!obj) {
	NUIS_ABORT("Object " << splitfile[1] << " doesn't exist!");
	}

	// Try casting
	TMatrixD mat = dynamic_cast<TMatrixD >(obj);
	if (mat) {
	TMatrixD newmat = (TMatrixD )mat->Clone();

	delete mat;
	tempfile->Close();

	return newmat;
	}

	TMatrixDSym matsym = dynamic_cast<TMatrixDSym >(obj);
	if (matsym) {
	TMatrixD *newmat = new TMatrixD(matsym->GetNrows(), matsym->GetNrows());
	for (int i = 0; i < matsym->GetNrows(); i++) {
	for (int j = 0; j < matsym->GetNrows(); j++) {
	(newmat)(i, j) = (matsym)(i, j);
	}
	}

	delete matsym;
	tempfile->Close();

	return newmat;
	}

	TH2D mathist = dynamic_cast<TH2D >(obj);
	if (mathist) {
	TMatrixD *newmat = new TMatrixD(mathist->GetNbinsX(), mathist->GetNbinsX());
	for (int i = 0; i < mathist->GetNbinsX(); i++) {
	for (int j = 0; j < mathist->GetNbinsX(); j++) {
	(*newmat)(i, j) = mathist->GetBinContent(i + 1, j + 1);
	}
	}

	delete mathist;
	tempfile->Close();

	return newmat;
	}

	return NULL;
	}

	//*******************************************************************
	TMatrixDSym *StatUtils::GetCovarFromTextFile(std::string covfile, int dim) {
	//*******************************************************************

	// Delete TempMat
	TMatrixD *tempmat = GetMatrixFromTextFile(covfile, dim, dim);

	// Make a symmetric covariance
	TMatrixDSym *newmat = new TMatrixDSym(tempmat->GetNrows());
	for (int i = 0; i < tempmat->GetNrows(); i++) {
	for (int j = 0; j < tempmat->GetNrows(); j++) {
	(newmat)(i, j) = (tempmat)(i, j);
	}
	}

	delete tempmat;
	return newmat;
	}

	//*******************************************************************
	TMatrixDSym *StatUtils::GetCovarFromRootFile(std::string covfile,
	std::string histname) {
	//*******************************************************************

	TMatrixD *tempmat = GetMatrixFromRootFile(covfile, histname);
	TMatrixDSym *newmat = new TMatrixDSym(tempmat->GetNrows());

	for (int i = 0; i < tempmat->GetNrows(); i++) {
	for (int j = 0; j < tempmat->GetNrows(); j++) {
	(newmat)(i, j) = (tempmat)(i, j);
	}
	}

	delete tempmat;
	return newmat;
	}
	diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
	index 9f3e0c8..ea0706d 100644
	--- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
	+++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
	@@ -1,229 +1,270 @@
	// 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 "T2K_SignalDef.h"

	#include "T2K_CC0pi_XSec_2DPcos_nu_I.h"

	//********************************************************************
	T2K_CC0pi_XSec_2DPcos_nu_I::T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey) {
	//********************************************************************

	// Sample overview ---------------------------------------------------
	std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu_I sample. \n"
	"Target: CH \n"
	"Flux: T2K 2.5 degree off-axis (ND280) \n"
	"Signal: CC0pi\n"
	"https://journals.aps.org/prd/abstract/10.1103/"
	"PhysRevD.93.112012 Analysis I";

	// Setup common settings
	fSettings = LoadSampleSettings(samplekey);
	fSettings.SetDescription(descrip);
	fSettings.SetXTitle("P_{#mu} (GeV)");
	fSettings.SetYTitle("cos#theta_{#mu}");
	fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
	fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV",
	"FIX/FULL");
	fSettings.SetEnuRange(0.0, 10.0);
	fSettings.DefineAllowedTargets("C,H");

	// CCQELike plot information
	fSettings.SetTitle("T2K_CC0pi_XSec_2DPcos_nu_I");
	fSettings.DefineAllowedSpecies("numu");

	FinaliseSampleSettings();

	// Scaling Setup ---------------------------------------------------
	// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
	fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) *
	1E-38 / (TotalIntegratedFlux()));

	// Setup Histograms
	SetHistograms();

	// Final setup ---------------------------------------------------
	FinaliseMeasurement();
	};

	bool T2K_CC0pi_XSec_2DPcos_nu_I::isSignal(FitEvent *event) {
	return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_I);
	};

	void T2K_CC0pi_XSec_2DPcos_nu_I::FillEventVariables(FitEvent *event) {

	if (event->NumFSParticle(13) == 0)
	return;

	TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
	TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;

	double pmu = Pmu.Vect().Mag() / 1000.;
	double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect()));

	fXVar = pmu;
	fYVar = CosThetaMu;

	return;
	};

	void T2K_CC0pi_XSec_2DPcos_nu_I::FillHistograms() {

	Measurement1D::FillHistograms();
	if (Signal) {
	fMCHist_Fine2D->Fill(fXVar, fYVar, Weight);
	FillMCSlice(fXVar, fYVar, Weight);
	}
	}

	// Modification is needed after the full reconfigure to move bins around
	// Otherwise this would need to be replaced by a TH2Poly which is too awkward.
	void T2K_CC0pi_XSec_2DPcos_nu_I::ConvertEventRates() {

	for (int i = 0; i < 9; i++) {
	fMCHist_Slices[i]->GetSumw2();
	}

	// Do standard conversion.
	Measurement1D::ConvertEventRates();

	// First scale MC slices also by their width in Y
	fMCHist_Slices[0]->Scale(1.0 / 1.00);
	fMCHist_Slices[1]->Scale(1.0 / 0.60);
	fMCHist_Slices[2]->Scale(1.0 / 0.10);
	fMCHist_Slices[3]->Scale(1.0 / 0.10);
	fMCHist_Slices[4]->Scale(1.0 / 0.05);
	fMCHist_Slices[5]->Scale(1.0 / 0.05);
	fMCHist_Slices[6]->Scale(1.0 / 0.04);
	fMCHist_Slices[7]->Scale(1.0 / 0.04);
	fMCHist_Slices[8]->Scale(1.0 / 0.02);

	// Now Convert into 1D list
	fMCHist->Reset();
	int bincount = 0;
	for (int i = 0; i < 9; i++) {
	for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
	fMCHist->SetBinContent(bincount + 1,
	fMCHist_Slices[i]->GetBinContent(j + 1));
	fMCHist->SetBinError(bincount + 1, fMCHist_Slices[i]->GetBinError(j + 1));
	bincount++;
	}
	}

	return;
	}

	void T2K_CC0pi_XSec_2DPcos_nu_I::FillMCSlice(double x, double y, double w) {

	if (y >= -1.0 and y < 0.0)
	fMCHist_Slices[0]->Fill(x, w);
	else if (y >= 0.0 and y < 0.6)
	fMCHist_Slices[1]->Fill(x, w);
	else if (y >= 0.6 and y < 0.7)
	fMCHist_Slices[2]->Fill(x, w);
	else if (y >= 0.7 and y < 0.8)
	fMCHist_Slices[3]->Fill(x, w);
	else if (y >= 0.8 and y < 0.85)
	fMCHist_Slices[4]->Fill(x, w);
	else if (y >= 0.85 and y < 0.90)
	fMCHist_Slices[5]->Fill(x, w);
	else if (y >= 0.90 and y < 0.94)
	fMCHist_Slices[6]->Fill(x, w);
	else if (y >= 0.94 and y < 0.98)
	fMCHist_Slices[7]->Fill(x, w);
	else if (y >= 0.98 and y <= 1.00)
	fMCHist_Slices[8]->Fill(x, w);
	}

	void T2K_CC0pi_XSec_2DPcos_nu_I::SetHistograms() {

	// Read in 1D Data Histograms
	fInputFile = new TFile(
	(FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root")
	.c_str(),
	"READ");
	// fInputFile->ls();

	// Read in 1D Data
	fDataHist = (TH1D *)fInputFile->Get("datahist");

	fMCHist_Fine2D = new TH2D("T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D",
	"T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", 400, 0.0, 30.0,
	100, -1.0, 1.0);
	SetAutoProcessTH1(fMCHist_Fine2D);

	TH2D tempcov = (TH2D )fInputFile->Get("analysis1_totcov");

	fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
	for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
	for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
	(*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
	}
	}
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);

	// Read in 2D Data
	fDataPoly = (TH2Poly *)fInputFile->Get("datapoly");
	fDataPoly->SetNameTitle("T2K_CC0pi_XSec_2DPcos_nu_I_datapoly",
	"T2K_CC0pi_XSec_2DPcos_nu_I_datapoly");
	SetAutoProcessTH1(fDataPoly, kCMD_Write);
	fDataHist->Reset();

	// Read in 2D Data Slices and Make MC Slices
	int bincount = 0;
	for (int i = 0; i < 9; i++) {

	// Get Data Histogram
	// fInputFile->ls();
	fDataHist_Slices.push_back(
	(TH1D *)fInputFile->Get(Form("dataslice_%i", i))->Clone());
	fDataHist_Slices[i]->SetNameTitle(
	Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%i", i),
	(Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%i", i)));

	// Loop over nbins and set errors from covar
	for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
	fDataHist_Slices[i]->SetBinError(
	j + 1, sqrt((fFullCovar)(bincount, bincount)) 1E-38);

	// std::cout << "Setting data hist " <<
	// fDataHist_Slices[i]->GetBinContent(j+1) << " " <<
	// fDataHist_Slices[i]->GetBinError(j+1) << std::endl;
	fDataHist->SetBinContent(bincount + 1,
	fDataHist_Slices[i]->GetBinContent(j + 1));
	fDataHist->SetBinError(bincount + 1,
	fDataHist_Slices[i]->GetBinError(j + 1));

	bincount++;
	}

	// Make MC Clones
	fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone());
	fMCHist_Slices[i]->SetNameTitle(
	Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%i", i),
	(Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%i", i)));

	SetAutoProcessTH1(fDataHist_Slices[i], kCMD_Write);
	SetAutoProcessTH1(fMCHist_Slices[i]);
	// fMCHist_Slices[i]->Reset();
	}

	return;
	};
	+
	+void T2K_CC0pi_XSec_2DPcos_nu_I::Write(std::string drawopt) {
	+ this->Measurement1D::Write(drawopt);
	+
	+ if (fResidualHist) {
	+ std::vector<TH1D *> MCResidual_Slices;
	+ size_t tb_it = 0;
	+ for (size_t i = 0; i < fMCHist_Slices.size(); ++i) {
	+ std::string name = Form("T2K_CC0pi_XSec_2DPcos_nu_I_RESIDUAL_Slice%i", i);
	+ MCResidual_Slices.push_back(
	+ static_cast<TH1D *>(fMCHist_Slices[i]->Clone(name.c_str())));
	+ MCResidual_Slices.back()->Reset();
	+
	+ for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) {
	+ double bc = fResidualHist->GetBinContent(tb_it + 1);
	+ MCResidual_Slices.back()->SetBinContent(j + 1, bc);
	+ tb_it++;
	+ }
	+ MCResidual_Slices.back()->Write();
	+ }
	+ }
	+
	+ if (fChi2LessBinHist) {
	+ std::vector<TH1D *> MCChi2LessBin_Slices;
	+ size_t tb_it = 0;
	+ for (size_t i = 0; i < fMCHist_Slices.size(); ++i) {
	+ std::string name =
	+ Form("T2K_CC0pi_XSec_2DPcos_nu_I_Chi2NMinusOne_Slice%i", i);
	+ MCChi2LessBin_Slices.push_back(
	+ static_cast<TH1D *>(fMCHist_Slices[i]->Clone(name.c_str())));
	+ MCChi2LessBin_Slices.back()->Reset();
	+
	+ for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) {
	+ double bc = fChi2LessBinHist->GetBinContent(tb_it + 1);
	+ MCChi2LessBin_Slices.back()->SetBinContent(j + 1, bc);
	+ tb_it++;
	+ }
	+ MCChi2LessBin_Slices.back()->Write();
	+ }
	+ }
	+}
	\ No newline at end of file
	diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h
	index 269a6de..a1f63b6 100644
	--- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h
	+++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h
	@@ -1,72 +1,74 @@
	// 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 T2K_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN
	#define T2K_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN

	#include "Measurement1D.h"
	#include "TH2Poly.h"
	#include "MeasurementVariableBox2D.h"

	class T2K_CC0pi_XSec_2DPcos_nu_I : public Measurement1D {
	public:

	T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey);

	/// Virtual Destructor
	~T2K_CC0pi_XSec_2DPcos_nu_I() {};

	/// Numu CC0PI Signal Definition
	///
	/// /item
	bool isSignal(FitEvent *nvect);

	/// Read histograms in a special way because format is different.
	/// Read from FitPar::GetDataBase()+"/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root"
	void SetHistograms();

	/// Bin Tmu CosThetaMu
	void FillEventVariables(FitEvent* customEvent);

	// Fill Histograms
	void FillHistograms();

	/// Have to do a weird event scaling for analysis 1
	void ConvertEventRates();

	+ void Write(std::string drawopt);
	+
	/// \brief Create Q2 Box to save correction info
	inline MeasurementVariableBox* CreateBox(){ return new MeasurementVariableBox2D(); };

	private:

	bool fIsSystCov, fIsStatCov, fIsNormCov;

	TH2Poly* fDataPoly;
	TH2Poly* fMCPoly;

	TFile* fInputFile;
	TH2D* fMCHist_Fine2D;

	std::vector<TH1D*> fMCHist_Slices;
	std::vector<TH1D*> fDataHist_Slices;

	void FillMCSlice(double x, double y, double w);

	};

	#endif

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 9:08 PM (1 d, 27 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3806206
Default Alt Text: (209 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions