No OneTemporary
Actions

Size

82 KB

Subscribers

None

View Options

	diff --git a/src/FitBase/Measurement1D.cxx b/src/FitBase/Measurement1D.cxx
	index 1579746..831dc72 100644
	--- a/src/FitBase/Measurement1D.cxx
	+++ b/src/FitBase/Measurement1D.cxx
	@@ -1,1781 +1,1787 @@
	// Copyright 2016 L. Pickering, P caltowell, 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;

	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";
	+ "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 (fCovar) delete fCovar;
	if (fInvert) delete fInvert;
	if (fDecomp) delete fDecomp;

	}

	//********************************************************************
	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
	LOG(SAM) << "Finalising Sample Settings: " << fName << std::endl;

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

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

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

	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) {
	ERR(WRN) << "Norm error for class " << fName << " is 0.0!" << std::endl;
	ERR(WRN) << "If you want to use it please add fNormError=VAL" << std::endl;
	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) {
	//********************************************************************

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

	}

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

	LOG(SAM) << "Reading data from root file: " << datafile << ";" << histname << std::endl;
	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) {
	ERR(FTL) << "Need a data hist to setup possion errors! " << std::endl;
	ERR(FTL) << "Setup Data First!" << std::endl;
	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) {
	LOG(SAM) << "Setting diagonal covariance for: " << data->GetName() << std::endl;
	fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
	covar = StatUtils::GetInvert(fFullCovar);
	fDecomp = StatUtils::GetDecomp(fFullCovar);
	} else {
	ERR(FTL) << "No data input provided to set diagonal covar from!" << std::endl;

	}

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

	}

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

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

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

	}

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

	LOG(SAM) << "Reading covariance from text file: " << covfile << ";" << histname << std::endl;
	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();
	}

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

	}

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

	LOG(SAM) << "Reading inverted covariance from text file: " << covfile << ";" << histname << std::endl;
	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();
	LOG(SAM) << "Reading data correlations from text file: " << covfile << ";" << dim << std::endl;
	TMatrixDSym* correlation = StatUtils::GetCovarFromTextFile(covfile, dim);

	if (!fDataHist) {
	ERR(FTL) << "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" << std::endl;
	throw;
	}

	// 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::SetCorrelationFromRootFile(std::string covfile, std::string histname) {
	//********************************************************************

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

	if (!fDataHist) {
	ERR(FTL) << "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" << std::endl;
	throw;
	}

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

	LOG(SAM) << "Reading cholesky from text file: " << covfile << std::endl;
	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) {
	//********************************************************************

	LOG(SAM) << "Reading cholesky decomp from root file: " << covfile << ";" << histname << std::endl;
	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::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;
	LOG(SAM) << "Reading bin mask from file: " << maskfile << std::endl;

	// 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(), ifstream::in);

	if (!mask.is_open()) {
	LOG(FTL) << " Cannot find mask file." << std::endl;
	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) {
	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.
	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() {
	//********************************************************************

	LOG(SAM) << "Finalising Measurement: " << fName << std::endl;

	// Make sure data is setup
	if (!fDataHist) {
	ERR(FTL) << "No data has been setup inside " << fName << " constructor!" << std::endl;
	throw;
	}

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

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

	if (!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(),
	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) {
	ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" << __LINE__ << std::endl;
	ERR(FTL) << "fScaleFactor = " << fScaleFactor << std::endl;
	ERR(FTL) << "EXITING" << std::endl;
	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) {
	ERR(FTL) << "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." << std::endl;
	throw;
	}
	}
	}

	// 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) {
	ERR(FTL) << "ERROR: Fit Option '" << fit_options_input.at(i)
	<< "' Provided is not allowed for this measurement."
	<< std::endl;
	ERR(FTL) << "Fit Options should be provided as a '/' seperated list "
	"(e.g. FREE/DIAG/NORM)"
	<< std::endl;
	ERR(FTL) << "Available options for " << fName << " are '" << fAllowedTypes
	<< "'" << std::endl;

	throw;
	}
	}

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

	ERR(FTL) << "No other LIKELIHOODS properly supported!" << std::endl;
	ERR(FTL) << "Try to use a chi2!" << std::endl;
	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(), ifstream::in);

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

	if (smear.is_open())
	LOG(SAM) << "Reading smearing matrix from file: " << smearfile << std::endl;
	else
	ERR(FTL) << "Smearing matrix provided is incorrect: " << smearfile
	<< std::endl;

	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) {
	ERR(WRN) << fName
	<< ": attempted to apply smearing matrix, but none was set"
	<< std::endl;
	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) {

	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) {
	// 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) 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) {
	std::cout << "Getting likelihood from covariance " << std::endl;
	stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist);
	}

	}

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

	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;
	LOG(SAM) << "Setting fake data from : " << fFakeDataInput << std::endl;

	// 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 (fDataHist) delete fDataHist;
	fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);

	return;
	};

	/*
	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 Data/MC
	GetDataList().at(0)->Write();
	GetMCList().at(0)->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())
	GetEventHistogram()->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());
	}

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

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

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

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

	// Returning
	LOG(SAM) << "Written Histograms: " << fName << std::endl;
	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());

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

	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.AddS("name", fName);
	// samplekey.AddS("type",type);
	// samplekey.AddS("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;
	LOG(SAM) << "Found event rate measurement but fIsRawEvents == false!"
	<< std::endl;
	LOG(SAM) << "Overriding this and setting fIsRawEvents == true!"
	<< std::endl;
	}

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

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

	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
	LOG(SAM) << "Reading data from: " << dataFile.c_str() << std::endl;
	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) {
	//********************************************************************

	LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
	<< std::endl;
	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) {
	//********************************************************************

	LOG(SAM) << "Filling histogram from " << inhistfile << "->" << histname
	<< std::endl;
	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* decmpPlot = new TH2D();
	TH2D* covarInvPlot = 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());
	covarInvPlot = (TH2D*)tempFile->Get((covName + "covinv").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
	ERR(WRN) << "Incorrect thrown_covariance option in parameters."
	<< std::endl;

	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(), ifstream::in);

	this->covar = new TMatrixDSym(dim);
	fFullCovar = new TMatrixDSym(dim);
	if (covarread.is_open())
	LOG(SAM) << "Reading covariance matrix from file: " << covarFile
	<< std::endl;
	else
	ERR(FTL) << "Covariance matrix provided is incorrect: " << covarFile
	<< std::endl;

	// 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) {
	ERR(WRN) << "SetCovarMatrixFromText -> Covariance matrix only has <= 1 "
	"entries on this line: "
	<< row << std::endl;
	}

	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(), ifstream::in);

	this->covar = new TMatrixDSym(dim);
	this->fFullCovar = new TMatrixDSym(dim);
	if (corr.is_open())
	LOG(SAM) << "Reading and converting correlation matrix from file: "
	<< corrFile << std::endl;
	else {
	ERR(FTL) << "Correlation matrix provided is incorrect: " << corrFile
	<< std::endl;
	exit(-1);
	}

	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) {
	ERR(FTL) << "Found a zero value in the covariance matrix, assuming "
	"this is an error!"
	<< std::endl;
	exit(-1);
	}

	(*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) {
	//********************************************************************

	LOG(SAM) << "Getting covariance from " << covarFile << "->" << covName
	<< std::endl;

	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(), 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/MeasurementBase.cxx b/src/FitBase/MeasurementBase.cxx
	index e0d4af3..d5675e4 100644
	--- a/src/FitBase/MeasurementBase.cxx
	+++ b/src/FitBase/MeasurementBase.cxx
	@@ -1,538 +1,542 @@
	// 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 "MeasurementBase.h"

	/*
	Constructor/Destructors
	*/

	//********************************************************************
	// 2nd Level Constructor (Inherits From MeasurementBase.h)
	MeasurementBase::MeasurementBase(void) {
	//********************************************************************

	fScaleFactor = 1.0;
	fMCFilled = false;
	fNoData = false;
	fInput = NULL;

	// Set the default values
	// After-wards this gets set in SetupMeasurement
	EnuMin = 0.;
	EnuMax = 1.E5;

	fMeasurementSpeciesType = kSingleSpeciesMeasurement;
	fEventVariables = NULL;
	fIsJoint = false;
	};

	void MeasurementBase::FinaliseMeasurement() {

	// Used to setup default data hists, covars, etc.



	}

	//********************************************************************
	// 2nd Level Destructor (Inherits From MeasurementBase.h)
	MeasurementBase::~MeasurementBase() {
	//********************************************************************

	};

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

	// Set Energy Limits
	if (low == -9999.9) low = this->EnuMin;
	if (high == -9999.9) high = this->EnuMax;
	return GetInput()->TotalIntegratedFlux(low, high, intOpt);
	};

	//********************************************************************
	double MeasurementBase::PredictedEventRate(std::string intOpt, double low,
	double high) {
	//********************************************************************

	// Set Energy Limits
	if (low == -9999.9) low = this->EnuMin;
	if (high == -9999.9) high = this->EnuMax;

	return GetInput()->PredictedEventRate(low, high, intOpt) * 1E-38;
	};

	//********************************************************************
	void MeasurementBase::SetupInputs(std::string inputfile) {
	//********************************************************************


	// Add this infile to the global manager
	if (FitPar::Config().GetParB("EventManager")) {
	fInput = FitBase::AddInput(fName, inputfile);
	} else {
	std::vector<std::string> file_descriptor =
	GeneralUtils::ParseToStr(inputfile, ":");
	if (file_descriptor.size() != 2) {
	ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfile
	<< "\". expected \"FILETYPE:file.root\"" << std::endl;
	throw;
	}
	InputUtils::InputType inpType =
	InputUtils::ParseInputType(file_descriptor[0]);

	fInput = InputUtils::CreateInputHandler(fName, inpType, file_descriptor[1]);
	}


	fNEvents = fInput->GetNEvents();

	// Expect INPUTTYPE:FileLocation(s)
	std::vector<std::string> file_descriptor =
	GeneralUtils::ParseToStr(inputfile, ":");
	if (file_descriptor.size() != 2) {
	ERR(FTL) << "File descriptor had no filetype declaration: \"" << inputfile
	<< "\". expected \"FILETYPE:file.root\"" << std::endl;
	throw;
	}
	fInputType = InputUtils::ParseInputType(file_descriptor[0]);

	fInputFileName = file_descriptor[1];
	if (EnuMin == 0 && EnuMax == 1.E5) {
	EnuMin = fInput->GetFluxHistogram()->GetBinLowEdge(1);
	EnuMax = fInput->GetFluxHistogram()->GetBinLowEdge(
	fInput->GetFluxHistogram()->GetNbinsX() + 1);
	}

	fFluxHist = fInput->GetFluxHistogram();
	fEventHist = fInput->GetEventHistogram();
	}

	//***********************************************
	int MeasurementBase::GetInputID() {
	//***********************************************
	return FitBase::GetInputID(fInputFileName);
	}

	//***********************************************
	SampleSettings MeasurementBase::LoadSampleSettings(nuiskey samplekey) {
	//***********************************************
	SampleSettings setting = SampleSettings(samplekey);
	fName = setting.GetS("name");

	// Used as an initial setup function incase we need to do anything here.
	LOG(SAM) << "Loading Sample : " << setting.GetName() << std::endl;
	if (!fIsJoint) SetupInputs( setting.GetS("input") );

	return setting;
	}

	//***********************************************
	SampleSettings MeasurementBase::LoadSampleSettings(std::string name, std::string input, std::string type) {
	//***********************************************

	nuiskey samplekey = Config::CreateKey("sample");
	samplekey.SetS("name",name);
	samplekey.SetS("input",input);
	samplekey.SetS("type",type);

	return LoadSampleSettings(samplekey);
	}

	void MeasurementBase::FinaliseSampleSettings() {

	EnuMin = fSettings.GetD("enu_min");
	EnuMax = fSettings.GetD("enu_max");

	}


	//***********************************************
	void MeasurementBase::Reconfigure() {
	//***********************************************

	LOG(REC) << " Reconfiguring sample " << fName << std::endl;

	// Reset Histograms
	ResetExtraHistograms();
	AutoResetExtraTH1();
	this->ResetAll();

	// FitEvent* cust_event = fInput->GetEventPointer();
	int fNEvents = fInput->GetNEvents();
	int countwidth = (fNEvents / 5);


	// MAIN EVENT LOOP
	FitEvent* cust_event = fInput->FirstNuisanceEvent();
	int i = 0;
	int npassed = 0;
	while(cust_event){

	cust_event->RWWeight = fRW->CalcWeight(cust_event);
	cust_event->Weight = cust_event->RWWeight * cust_event->InputWeight;

	Weight = cust_event->Weight;

	// Initialize
	fXVar = -999.9;
	fYVar = -999.9;
	fZVar = -999.9;
	Signal = false;
	Mode = cust_event->Mode;

	// Extract Measurement Variables
	this->FillEventVariables(cust_event);
	Signal = this->isSignal(cust_event);
	if (Signal) npassed++;

	GetBox()->SetX(fXVar);
	GetBox()->SetY(fYVar);
	GetBox()->SetZ(fZVar);
	GetBox()->SetMode(Mode);
	// GetBox()->fSignal = Signal;

	// Fill Histogram Values
	GetBox()->FillBoxFromEvent(cust_event);
	// this->FillExtraHistograms(GetBox(), Weight);
	this->FillHistogramsFromBox(GetBox(), Weight);

	// Print Out
	if (LOG_LEVEL(REC) && countwidth > 0 && !(i % countwidth)) {
	std::stringstream ss("");
	ss.unsetf(std::ios_base::fixed);
	ss << std::setw(7) << std::right << i << "/" << fNEvents << " events ("
	<< std::setw(2) << double(i) / double(fNEvents) * 100. << std::left
	<< std::setw(5) << "%) "
	<< "[S,X,Y,Z,M,W] = [" << std::fixed << std::setprecision(2)
	<< std::right << Signal << ", " << std::setw(5) << fXVar << ", "
	<< std::setw(5) << fYVar << ", " << std::setw(5) << fYVar << ", "
	<< std::setw(3) << (int)Mode << ", " << std::setw(5) << Weight << "] "
	<< std::endl;
	LOG(SAM) << ss.str();
	}

	// iterate
	cust_event = fInput->NextNuisanceEvent();
	i++;
	}

	LOG(SAM) << npassed << "/" << fNEvents << " passed selection " << std::endl;
	if (npassed == 0) {
	LOG(SAM) << "WARNING: NO EVENTS PASSED SELECTION!" << std::endl;
	}
	LOG(REC) << std::setw(10) << std::right << NSignal << "/"
	<< fNEvents << " events passed selection + binning after reweight"
	<< std::endl;

	// Finalise Histograms
	fMCFilled = true;
	this->ConvertEventRates();
	}

	void MeasurementBase::FillHistogramsFromBox(MeasurementVariableBox* var, double weight) {

	fXVar = var->GetX();
	fYVar = var->GetY();
	fZVar = var->GetZ();
	// Signal = var->fSignal;
	// Mode = var->fMode;
	Weight = weight;

	FillHistograms();
	FillExtraHistograms(var, weight);

	}

	void MeasurementBase::FillHistograms(double weight){
	Weight = weight * GetBox()->GetSampleWeight();
	FillHistograms();
	FillExtraHistograms(GetBox(), Weight);
	}


	MeasurementVariableBox* MeasurementBase::FillVariableBox(FitEvent* event) {

	GetBox()->Reset();
	Mode = event->Mode;
	Weight = 1.0; //event->Weight;

	this->FillEventVariables(event);
	Signal = this->isSignal(event);

	GetBox()->FillBoxFromEvent(event);

	GetBox()->SetX(fXVar);
	GetBox()->SetY(fYVar);
	GetBox()->SetZ(fZVar);
	GetBox()->SetMode(event->Mode);
	GetBox()->SetSampleWeight(Weight);
	// GetBox()->fSignal = Signal;

	return GetBox();
	}

	MeasurementVariableBox* MeasurementBase::GetBox() {
	if (!fEventVariables) fEventVariables = CreateBox();
	return fEventVariables;
	}

	//***********************************************
	void MeasurementBase::ReconfigureFast() {
	//***********************************************
	this->Reconfigure();
	}

	//***********************************************
	void MeasurementBase::ConvertEventRates() {
	//***********************************************

	AutoScaleExtraTH1();
	ScaleExtraHistograms(GetBox());
	this->ScaleEvents();

	double normval = fRW->GetSampleNorm(this->fName);
	if (normval < 0.01 or normval > 10.0){
	ERR(WRN) << "Norm Value inside MeasurementBase::ConvertEventRates() looks off!" << std::endl;
	ERR(WRN) << "It could have become out of sync with the minimizer norm list." << std::endl;
	ERR(WRN) << "Setting it to 1.0" << std::endl;
	normval = 1.0;
	}
	AutoNormExtraTH1(normval);
	NormExtraHistograms(GetBox(), normval);
	this->ApplyNormScale(normval);

	}

	//***********************************************
	InputHandlerBase* MeasurementBase::GetInput() {
	//***********************************************

	if (!fInput) {
	ERR(FTL) << "MeasurementBase::fInput not set. Please submit your command "
	"line options and input cardfile with a bug report to: "
	"nuisance@projects.hepforge.org"
	<< std::endl;
	throw;
	}
	return fInput;
	};

	//***********************************************
	void MeasurementBase::Renormalise() {
	//***********************************************

	// Called when the fitter has changed a measurements normalisation but not any
	// reweight dials
	// Means we don't have to call the time consuming reconfigure when this
	// happens.
	double norm = fRW->GetDialValue(this->fName + "_norm");

	if ((this->fCurrentNorm == 0.0 and norm != 0.0) or not fMCFilled) {
	this->ReconfigureFast();
	return;
	}

	if (this->fCurrentNorm == norm) return;

	this->ApplyNormScale(1.0 / this->fCurrentNorm);
	this->ApplyNormScale(norm);

	return;
	};

	//***********************************************
	void MeasurementBase::SetSignal(bool sig) {
	//***********************************************
	Signal = sig;
	}

	//***********************************************
	void MeasurementBase::SetSignal(FitEvent* evt) {
	//***********************************************
	Signal = this->isSignal(evt);
	}

	//***********************************************
	void MeasurementBase::SetWeight(double wght) {
	//***********************************************
	Weight = wght;
	}

	//***********************************************
	void MeasurementBase::SetMode(int md) {
	//***********************************************
	Mode = md;
	}

	//***********************************************
	std::vector<TH1*> MeasurementBase::GetFluxList() {
	//***********************************************
	return GetInput()->GetFluxList();
	}

	//***********************************************
	std::vector<TH1*> MeasurementBase::GetEventRateList() {
	//***********************************************
	return GetInput()->GetEventList();
	}

	//***********************************************
	std::vector<TH1*> MeasurementBase::GetXSecList() {
	//***********************************************
	return GetInput()->GetXSecList();
	}


	void MeasurementBase::ProcessExtraHistograms(int cmd,
	MeasurementVariableBox* vars,
	double weight) {
	// This should be overriden if we have extra histograms!!!
	// Add a flag to tell user this...
	return;
	}

	void MeasurementBase::FillExtraHistograms(MeasurementVariableBox* vars,
	double weight) {
	ProcessExtraHistograms(kCMD_Fill, vars, weight);
	}

	void MeasurementBase::ScaleExtraHistograms(MeasurementVariableBox* vars) {
	ProcessExtraHistograms(kCMD_Scale, vars, 1.0);
	}
	void MeasurementBase::ResetExtraHistograms() {
	ProcessExtraHistograms(kCMD_Reset, NULL, 1.0);
	}
	void MeasurementBase::NormExtraHistograms(MeasurementVariableBox* vars,
	double norm) {
	ProcessExtraHistograms(kCMD_Norm, vars, norm);
	}
	void MeasurementBase::WriteExtraHistograms() {
	ProcessExtraHistograms(kCMD_Write, NULL, 1.00);
	}

	void MeasurementBase::SetAutoProcessTH1(TH1* hist, int c1, int c2, int c3, int c4, int c5) {
	FakeStack* fake = new FakeStack(hist);
	SetAutoProcessTH1(fake, c1, c2, c3, c4, c5); // Need to add a destroy command!
	}

	void MeasurementBase::SetAutoProcessTH1(StackBase* hist, int c1, int c2, int c3, int c4, int c5) {

	// Set Defaults
	// int ncommands = kCMD_extraplotflags;
	bool autoflags[5];
	autoflags[0] = false;
	autoflags[1] = false;
	autoflags[2] = false;
	autoflags[3] = false;
	autoflags[4] = false;

	int givenflags[5];
	givenflags[0] = c1;
	givenflags[1] = c2;
	givenflags[2] = c3;
	givenflags[3] = c4;
	givenflags[4] = c5;
	fExtraTH1s[hist] = std::vector<int>(5,0);

	// Setup a default one.
	if (c1 == -1 && c2 == -1 && c3 == -1 && c4 == -1 && c5 == -1){
	fExtraTH1s[hist][kCMD_Reset] = 1;
	fExtraTH1s[hist][kCMD_Scale] = 1;
	fExtraTH1s[hist][kCMD_Norm] = 1;
	fExtraTH1s[hist][kCMD_Write] = 1;
	}

	for (int i = 0; i < 5; i++) {
	switch (givenflags[i]) {

	// Skip over...
	case -1:
	break;

	case kCMD_Reset:
	case kCMD_Scale:
	case kCMD_Norm:
	case kCMD_Write:
	fExtraTH1s[hist][givenflags[i]] = 1;
	break;

	case kCMD_Fill:
	ERR(FTL) << "Can't auto fill yet!" << std::endl;
	autoflags[givenflags[i]] = 1;
	break;

	default:
	break;
	}
	}

	// LOG(SAM) << "AutoProcessing " << hist->GetName() << std::endl;
	};

	void MeasurementBase::AutoFillExtraTH1() {
	ERR(FTL) << "Can't auto fill yet! it's too inefficent!" << std::endl;
	return;
	}

	void MeasurementBase::AutoResetExtraTH1() {

	for (std::map<StackBase*, std::vector<int> >::iterator iter = fExtraTH1s.begin();
	iter != fExtraTH1s.end(); iter++) {

	if (!((*iter).second)[kCMD_Reset]) continue;
	(*iter).first->Reset();
	}
	};

	void MeasurementBase::AutoScaleExtraTH1() {
	for (std::map<StackBase*, std::vector<int> >::iterator iter = fExtraTH1s.begin();
	iter != fExtraTH1s.end(); iter++) {

	if (!((*iter).second)[kCMD_Scale]) continue;
	- (*iter).first->Scale(fScaleFactor, "width");
	+ if (fIsNoWidth){
	+ (*iter).first->Scale(fScaleFactor);
	+ } else {
	+ (*iter).first->Scale(fScaleFactor, "width");
	+ }
	}
	};

	void MeasurementBase::AutoNormExtraTH1(double norm) {
	double sfactor = 0.0;
	if (norm != 0.0) sfactor = 1.0 / norm;

	for (std::map<StackBase*, std::vector<int> >::iterator iter = fExtraTH1s.begin();
	iter != fExtraTH1s.end(); iter++) {

	if (!((*iter).second)[kCMD_Norm]) continue;
	(*iter).first->Scale(sfactor);
	}
	};

	void MeasurementBase::AutoWriteExtraTH1() {
	for (std::map<StackBase*, std::vector<int> >::iterator iter = fExtraTH1s.begin();
	iter != fExtraTH1s.end(); iter++) {

	if (!(((*iter).second)[kCMD_Write])) continue;
	(*iter).first->Write();
	}
	};



	diff --git a/src/FitBase/MeasurementBase.h b/src/FitBase/MeasurementBase.h
	index ff86ab8..fb456f9 100644
	--- a/src/FitBase/MeasurementBase.h
	+++ b/src/FitBase/MeasurementBase.h
	@@ -1,343 +1,344 @@
	// 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 INPUTHANDLER_H_SEEN
	#define INPUTHANDLER_H_SEEN

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

	// C Includes
	#include <math.h>
	#include <stdlib.h>
	#include <time.h>
	#include <deque>
	#include <iomanip>
	#include <iostream>
	#include <list>
	#include <numeric>
	#include <sstream>
	#include <string>

	// ROOT includes
	#include <TArrayF.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 "FitParameters.h"
	#include "FitUtils.h"
	#include "GeneralUtils.h"
	#include "PlotUtils.h"
	#include "StatUtils.h"
	#include "InputFactory.h"
	#include "FitWeight.h"

	#include "TMultiDimFit.h"

	#ifdef __GENIE_ENABLED__
	#include "Conventions/Units.h"
	#endif

	#include "EventManager.h"
	#include "TObject.h"
	#include "InputHandler.h"
	#include "NuisConfig.h"
	#include "NuisKey.h"
	#include "SampleSettings.h"
	#include "StackBase.h"
	#include "StandardStacks.h"

	/// Enumerations to help with extra plot functions
	enum extraplotflags {
	kCMD_Reset = 0,
	kCMD_Fill,
	kCMD_Scale,
	kCMD_Norm,
	kCMD_Write,
	kCMD_Error,
	kCMD_extraplotflags
	};

	enum MeasurementSpeciesClass {
	kSingleSpeciesMeasurement = 0,
	kNumuWithWrongSignMeasurement,
	kNueWithWrongSignMeasurement,
	kFourSpeciesMeasurement,
	};

	/// InputHandler Class
	///
	/// Inherits from Measurement base to handle whatever input is throwna t the
	/// fitter automatically.
	/// All functions here handle how files are read in, converted to custom formats
	/// and reconfigures are called.
	/// Used generally for the MC inputs.

	//! 2nd level experiment class that handles converting MC into a common format
	//! and calling reconfigure


	class MeasurementBase {
	public:
	/*
	Constructor/Destructors
	*/
	//! Default Constructor. Set everything to NULL
	MeasurementBase();

	//! Default virtual destructor
	virtual ~MeasurementBase(void);
	virtual void InitialSetup(void) {};

	/*
	Reconfigure Functions
	*/
	//! Function called if MC tuning dials haven't been changed and all we want to
	//! do is update the normalisation.
	virtual void Renormalise(void);

	//! Call reconfigure only looping over signal events to save time.
	virtual void ReconfigureFast(void);

	virtual void FillHistograms(double weight);


	//! Call reconfigure looping over all MC events including background
	virtual void Reconfigure(void);

	// virtual TH2D GetCovarMatrix(void) = 0;
	virtual double GetLikelihood(void) { return 0.0; };
	virtual int GetNDOF(void) { return 0; };
	virtual void ThrowCovariance(void) = 0;
	virtual void SetFakeDataValues(std::string fkdt) = 0;

	//! Get the total integrated flux between this samples energy range
	virtual double TotalIntegratedFlux(std::string intOpt = "width",
	double low = -9999.9,
	double high = -9999.9);

	//! Get the predicted event rate for this sample
	virtual double PredictedEventRate(std::string intOpt = "width",
	double low = -9999.9,
	double high = -9999.9);

	virtual SampleSettings LoadSampleSettings(nuiskey samplekey);
	virtual SampleSettings LoadSampleSettings(std::string name, std::string input, std::string type);

	virtual void FinaliseSampleSettings();
	virtual void FinaliseMeasurement();
	virtual void ProcessExtraHistograms(int cmd, MeasurementVariableBox* vars,
	double weight = 1.0);

	virtual void FillExtraHistograms(MeasurementVariableBox* vars, double weight = 1.0);
	virtual void ScaleExtraHistograms(MeasurementVariableBox* vars);
	virtual void ResetExtraHistograms();
	virtual void NormExtraHistograms(MeasurementVariableBox* vars, double norm = 1.0);
	virtual void WriteExtraHistograms();
	virtual MeasurementVariableBox* CreateBox() {return new MeasurementVariableBox();};

	int GetPassed() {
	int signalSize = 0;
	return signalSize;
	}

	int GetTotal() { return fNEvents; }

	/*
	Reconfigure LOOP
	*/
	// All these should be virtual
	///! Reset Histograms (Handled at Measurement Stage)
	virtual void ResetAll(void) = 0;

	///! Fill the event variables for this sample (Handled in each inherited
	/// sample)
	virtual void FillEventVariables(FitEvent* event) { (void)event; };

	///! Check whether this event is signle (Handled in each inherited sample)
	virtual bool isSignal(FitEvent* event) {
	(void)event;
	return false;
	};

	///! Fill the histogram for this event using fXVar and fYVar (Handled in each
	/// inherited sample)
	virtual void FillHistograms(void) {};

	///! Convert event rates to whatever distributions you need.
	virtual void ConvertEventRates(void);

	///! Call scale events after the plots have been filled at the end of
	/// reconfigure.
	virtual void ScaleEvents(void) {};

	///! Apply the scale factor at the end of reconfigure.
	virtual void ApplyNormScale(double norm) { (void)norm; };

	///! Save Histograms
	virtual void Write(std::string drawOpt = "") = 0;

	virtual MeasurementVariableBox* FillVariableBox(FitEvent* event);

	virtual MeasurementVariableBox* GetBox();

	void FillHistogramsFromBox(MeasurementVariableBox* var, double weight);
	/*
	Histogram Access Functions
	*/

	///! Virtual function to get data histogram
	virtual std::vector<TH1*> GetDataList(void) = 0;

	///! Virtual function to get MC histogram
	virtual std::vector<TH1*> GetMCList(void) = 0;
	virtual std::vector<TH1*> GetFineList(void) = 0;
	virtual std::vector<TH1*> GetMaskList(void) = 0;

	///! Return flux histograms in a vector
	virtual std::vector<TH1*> GetFluxList(void);
	virtual std::vector<TH1*> GetEventRateList(void);
	virtual std::vector<TH1*> GetXSecList(void);

	virtual TH1D* GetEventHistogram() { return fInput->GetEventHistogram(); };
	virtual TH1D* GetXSecHistogram() { return fInput->GetXSecHistogram(); };
	virtual TH1D* GetFluxHistogram() { return fInput->GetFluxHistogram(); };

	///! Return input for this sample
	InputHandlerBase* GetInput(void);

	std::string GetName(void) { return fName; };
	double GetScaleFactor(void) { return fScaleFactor; };

	double GetXVar(void) { return fXVar; };
	double GetYVar(void) { return fYVar; };
	double GetZVar(void) { return fZVar; };
	double GetMode(void) { return this->Mode; };
	double GetEnu(void) { return this->Enu; };

	void SetupInputs(std::string inputfile);
	int GetInputID(void);
	std::string GetInputFileName() { return fInputFileName; };
	void SetSignal(bool sig);
	void SetSignal(FitEvent* evt);
	void SetWeight(double wght);
	void SetMode(int md);
	void SetNoData(bool isTrue = true) { fNoData = isTrue; };

	inline void SetXVar(double xvar) { fXVar = xvar; };
	inline void SetYVar(double yvar) { fYVar = yvar; };
	inline void SetZVar(double zvar) { fZVar = zvar; };

	virtual std::vector<MeasurementBase*> GetSubSamples() {
	return std::vector<MeasurementBase*>(1, this);
	}


	void SetAutoProcessTH1(TH1* hist, int c1 = -1,
	int c2 = -1, int c3 = -1,
	int c4 = -1, int c5 = -1);
	void SetAutoProcessTH1(StackBase* hist, int c1 = -1,
	int c2 = -1, int c3 = -1,
	int c4 = -1, int c5 = -1);
	void AutoFillExtraTH1();
	void AutoResetExtraTH1();
	void AutoScaleExtraTH1();
	void AutoNormExtraTH1(double norm);
	void AutoWriteExtraTH1();


	// functions that need to be added.
	// - Initial Check
	// - Check Target/Beam loop.
	// - Check flux shape if suggested one given.
	// - Return MeasurementList (returns )


	protected:
	// Minimum and maximum energies
	double Enu; //!< Neutrino Energy
	double EnuMin; //!< Minimum incoming particle energy of events to include
	double EnuMax; //!< Maximum incoming particle energy of events to include

	BaseFitEvt* signal_event;
	FitEvent* cust_event;

	FitWeight* fRW; //!< Pointer to the rw engine
	InputHandlerBase* fInput; //!< Instance of the input handler

	std::string fName; //!< Name of the sample
	int fEventType;

	double fBeamDistance; //!< Incoming Particle flight distance (for oscillation
	//! analysis)
	double fScaleFactor; //!< fScaleFactor applied to events to convert from
	//! eventrate to final distribution
	double
	fCurrentNorm; //!< current normalisation factor applied if fit is "FREE"
	bool fMCFilled; //!< flag whether MC plots have been filled (For
	//! ApplyNormalisation)
	bool fNoData; //!< flag whether data plots do not exist (for ratios)
	+ bool fIsNoWidth; ///< Flag : Don't scale by bin width

	// TEMP OBJECTS TO HANDLE MERGE
	double fXVar, fYVar, fZVar, Mode, Weight;
	bool Signal;
	int ievt;
	int fNEvents;
	double Enu_rec, ThetaMu, CosThetaMu;

	InputUtils::InputType fInputType;
	std::string fInputFileName;
	TH1D* fFluxHist;
	TH1D* fEventHist;

	MeasurementSpeciesClass fMeasurementSpeciesType;
	SampleSettings fSettings;

	MeasurementVariableBox* fEventVariables;

	std::map<StackBase*, std::vector<int> > fExtraTH1s;
	int NSignal;
	// std::map<TH1*, bool[6] > fExtaStacks;

	bool fIsJoint;

	};




	// Class TypeDefs
	typedef std::list<MeasurementBase*>::const_iterator MeasListConstIter;
	typedef std::list<MeasurementBase*>::iterator MeasListIter;
	typedef std::vector<MeasurementBase*>::const_iterator MeasVectConstIter;
	typedef std::vector<MeasurementBase*>::iterator MeasVectIter;

	/! @} /
	#endif

File Metadata

Mime Type: text/x-diff
Expires: Sat, Dec 21, 5:35 PM (12 h, 5 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023660
Default Alt Text: (82 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions