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