diff --git a/src/FitBase/ParamPull.cxx b/src/FitBase/ParamPull.cxx
index cd101bc..449789a 100644
--- a/src/FitBase/ParamPull.cxx
+++ b/src/FitBase/ParamPull.cxx
@@ -1,849 +1,857 @@
// 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 .
*******************************************************************************/
#include "ParamPull.h"
//*******************************************************************************
ParamPull::ParamPull(std::string name, std::string inputfile, std::string type, std::string dials) {
//*******************************************************************************
fMinHist = NULL;
fMaxHist = NULL;
fTypeHist = NULL;
fDialSelection = dials;
fLimitHist = NULL;
fName = name;
fInput = inputfile;
fType = type;
// Set the pull type
SetType(fType);
- std::cout << fType << std::endl;
// Setup Histograms from input file
SetupHistograms(fInput);
};
//*******************************************************************************
void ParamPull::SetType(std::string type) {
//*******************************************************************************
fType = type;
// Assume Default if empty
if (type.empty() || type == "DEFAULT") {
ERR(WRN) << "No type specified for ParmPull class " << fName << std::endl;
ERR(WRN) << "Assuming GAUSTHROW/GAUSPULL" << std::endl;
type = "GAUSTHROW/GAUSPULL";
}
// Set Dial options
if (type.find("FRAC") != std::string::npos) {
fDialOptions = "FRAC";
fPlotTitles = ";; Fractional RW Value";
} else if (type.find("ABS") != std::string::npos) {
fDialOptions = "ABS";
fPlotTitles = ";; ABS RW Value";
} else {
fDialOptions = "";
fPlotTitles = ";; RW Value";
}
// Parse throw types
if (type.find("GAUSPULL") != std::string::npos) fCalcType = kGausPull;
else fCalcType = kNoPull;
if (type.find("GAUSTHROW") != std::string::npos) fThrowType = kGausThrow;
else if (type.find("FLATTHROW") != std::string::npos) fThrowType = kFlatThrow;
else fThrowType = kNoThrow;
// Extra check to see if throws or pulls are turned off
if (type.find("NOPULL") != std::string::npos) fCalcType = kNoPull;
if (type.find("NOTHROW") != std::string::npos) fThrowType = kNoThrow;
}
//*******************************************************************************
void ParamPull::SetupHistograms(std::string input) {
//*******************************************************************************
// Extract Types from Input
fFileType = "";
const int nfiletypes = 4;
const std::string filetypes[nfiletypes] = {"FIT", "ROOT", "TXT", "DIAL"};
for (int i = 0; i < nfiletypes; i++) {
std::string tempTypes = filetypes[i] + ":";
if (input.find(tempTypes) != std::string::npos) {
fFileType = filetypes[i];
input.replace(input.find(tempTypes), tempTypes.size(), "");
break;
}
}
// Read Files
if (!fFileType.compare("FIT")) ReadFitFile(input);
else if (!fFileType.compare("ROOT")) ReadRootFile(input);
else if (!fFileType.compare("VECT")) ReadVectFile(input);
else if (!fFileType.compare("DIAL")) ReadDialInput(input);
else {
ERR(FTL) << "Unknown ParamPull Type: " << input << std::endl;
+ ERR(FTL) << "Need FIT, ROOT, VECT or DIAL" << std::endl;
throw;
}
// Check Dials are all good
- CheckDialsValid();
+ if (!CheckDialsValid()) {
+ ERR(FTL) << "DIALS NOT VALID" << std::endl;
+ throw;
+ }
// Setup MC Histogram
fMCHist = (TH1D*) fDataHist->Clone();
fMCHist->Reset();
fMCHist->SetNameTitle( (fName + "_MC").c_str(),
(fName + " MC" + fPlotTitles).c_str() );
// If no Covar input make an uncorrelated one
if (!fCovar) {
fCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist, 1.0);
}
// If no types or limits are provided give them a default option
if (!fMinHist) {
+ LOG(FIT) << "No minimum histogram found for pull parameters, setting to be content - 1E6..." << std::endl;
fMinHist = (TH1D*) fDataHist->Clone();
fMinHist->SetNameTitle( (fName + "_min").c_str(),
(fName + " min" + fPlotTitles).c_str() );
for (int i = 0; i < fMinHist->GetNbinsX(); i++) {
// TODO (P.Stowell) Change this to a NULL system where limits are actually free!
fMinHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6);
}
}
if (!fMaxHist) {
+ LOG(FIT) << "No maximum histogram found for pull parameters, setting to be content - 1E6..." << std::endl;
fMaxHist = (TH1D*) fDataHist->Clone();
fMaxHist->SetNameTitle( (fName + "_min").c_str(),
(fName + " min" + fPlotTitles).c_str() );
for (int i = 0; i < fMaxHist->GetNbinsX(); i++) {
fMaxHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6);
}
}
// Set types from state, or to unknown
+ // Not really sure when or if this is ever used
if (!fTypeHist) {
-
int deftype = -1;
if (fType.find("T2K") != std::string::npos) { deftype = kT2K; }
else if (fType.find("NEUT") != std::string::npos) { deftype = kNEUT; }
else if (fType.find("NIWG") != std::string::npos) { deftype = kNIWG; }
else if (fType.find("GENIE") != std::string::npos) { deftype = kGENIE; }
else if (fType.find("NORM") != std::string::npos) { deftype = kNORM; }
else if (fType.find("NUWRO") != std::string::npos) { deftype = kNUWRO; }
fTypeHist = new TH1I( (fName + "_type").c_str(),
(fName + " type" + fPlotTitles).c_str(),
fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX() );
for (int i = 0; i < fTypeHist->GetNbinsX(); i++) {
fTypeHist->SetBinContent(i + 1, deftype );
}
}
// Sort Covariances
fInvCovar = StatUtils::GetInvert(fCovar);
fDecomp = StatUtils::GetDecomp(fCovar);
// Create DataTrue for Throws
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
(fName + " truedata" + fPlotTitles).c_str() );
fDataOrig = (TH1D*) fDataHist->Clone();
fDataOrig->SetNameTitle( (fName + "_origdata").c_str(),
(fName + " origdata" + fPlotTitles).c_str() );
// Select only dials we want
if (!fDialSelection.empty()) {
(*fDataHist) = RemoveBinsNotInString(*fDataHist, fDialSelection);
-
-
-
}
}
//*******************************************************************************
TH1D ParamPull::RemoveBinsNotInString(TH1D hist, std::string mystr) {
//*******************************************************************************
// Make list of allowed bins
std::vector allowedbins;
for (int i = 0; i < hist.GetNbinsX(); i++) {
std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1));
if (mystr.find(syst) != std::string::npos) {
allowedbins.push_back(syst);
}
}
// Make new histogram
UInt_t nbins = allowedbins.size();
TH1D newhist = TH1D( hist.GetName(), hist.GetTitle(),
(Int_t)nbins, 0.0, (Double_t)nbins);
// Setup bins
for (UInt_t i = 0; i < nbins; i++) {
// Set Labels
newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str());
// Copy Values
for (Int_t j = 0; j < hist.GetNbinsX(); j++) {
if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) {
newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) );
newhist.SetBinError(i + 1, hist.GetBinError(j + 1) );
}
}
}
return newhist;
}
//*******************************************************************************
TH1I ParamPull::RemoveBinsNotInString(TH1I hist, std::string mystr) {
//*******************************************************************************
// Make list of allowed bins
std::vector allowedbins;
for (int i = 0; i < hist.GetNbinsX(); i++) {
std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1));
if (mystr.find(syst) != std::string::npos) {
allowedbins.push_back(syst);
}
}
// Make new histogram
UInt_t nbins = allowedbins.size();
TH1I newhist = TH1I( hist.GetName(), hist.GetTitle(),
(Int_t)nbins, 0.0, (Int_t)nbins);
// Setup bins
for (UInt_t i = 0; i < nbins; i++) {
// Set Labels
newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str());
// Copy Values
for (Int_t j = 0; j < hist.GetNbinsX(); j++) {
if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) {
newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) );
newhist.SetBinError(i + 1, hist.GetBinError(j + 1) );
}
}
}
return newhist;
}
//*******************************************************************************
void ParamPull::ReadFitFile(std::string input) {
//*******************************************************************************
TFile* tempfile = new TFile(input.c_str(), "READ");
// Read Data
fDataHist = (TH1D*) tempfile->Get("fit_dials_free");
- if (!fDataHist) {
- ERR(FTL) << "Can't find TH1D hist fit_dials in " << fName << std::endl;
- ERR(FTL) << "File Entries:" << std::endl;
- tempfile->ls();
-
- throw;
- }
-
+ CheckHist(fDataHist);
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
+ fMinHist = (TH1D*)tempfile->Get("min_dials_free");
+ CheckHist(fMinHist);
+ fMinHist->SetDirectory(0);
+ fMinHist->SetNameTitle( (fName + "_min").c_str(),
+ (fName + " min" + fPlotTitles).c_str() );
+
+ fMaxHist = (TH1D*)tempfile->Get("max_dials_free");
+ CheckHist(fMaxHist);
+ fMaxHist->SetDirectory(0);
+ fMaxHist->SetNameTitle( (fName + "_max").c_str(),
+ (fName + " max" + fPlotTitles).c_str() );
+
// Read Covar
TH2D* tempcov = (TH2D*) tempfile->Get("covariance_free");
if (!tempcov) {
ERR(FTL) << "Can't find TH2D covariance_free in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
// Setup Covar
int nbins = fDataHist->GetNbinsX();
fCovar = new TMatrixDSym( nbins );
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
(*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
return;
}
//*******************************************************************************
void ParamPull::ReadRootFile(std::string input) {
//*******************************************************************************
std::vector inputlist = GeneralUtils::ParseToStr(input, ";");
// Check all given
if (inputlist.size() < 2) {
ERR(FTL) << "Covar supplied in 'ROOT' format should have 3 semi-colon seperated entries!" << std::endl
<< "ROOT:filename;histname[;covarname]" << std::endl;
ERR(FTL) << "histname = TH1D, covarname = TH2D" << std::endl;
throw;
}
// Get Entries
std::string filename = inputlist[0];
- LOG(DEB) << filename << std::endl;
std::string histname = inputlist[1];
- LOG(DEB) << histname << std::endl;
- LOG(DEB) << input << std::endl;
// Read File
TFile* tempfile = new TFile(filename.c_str(), "READ");
if (tempfile->IsZombie()){
LOG(FIT) << "Looking for ParamPull input inside database" << std::endl;
filename = FitPar::GetDataBase() + "/" + filename;
tempfile = new TFile(filename.c_str(), "READ");
}
if (tempfile->IsZombie()) {
ERR(FTL) << "Can't find file in " << fName << std::endl;
ERR(FTL) << "location = " << filename << std::endl;
throw;
}
// Read Hist
fDataHist = (TH1D*) tempfile->Get(histname.c_str());
if (!fDataHist) {
ERR(FTL) << "Can't find TH1D hist " << histname << " in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
fDataHist->SetDirectory(0);
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
LOG(DEB) << "READING COVAR" << std::endl;
// Read Covar
if (inputlist.size() > 2) {
std::string covarname = inputlist[2];
LOG(DEB) << "COVARNAME = " << covarname << std::endl;
TH2D* tempcov = (TH2D*) tempfile->Get(covarname.c_str());
if (!tempcov) {
ERR(FTL) << "Can't find TH2D covar " << covarname << " in " << fName << std::endl;
ERR(FTL) << "File Entries:" << std::endl;
tempfile->ls();
throw;
}
// Setup Covar
int nbins = fDataHist->GetNbinsX();
fCovar = new TMatrixDSym( nbins );
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
(*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
// Uncorrelated
} else {
LOG(SAM) << "No Covar provided so using diagonal errors for "
<< fName << std::endl;
fCovar = NULL;
}
}
//*******************************************************************************
void ParamPull::ReadVectFile(std::string input) {
//*******************************************************************************
std::vector inputlist = GeneralUtils::ParseToStr(input, ";");
if (inputlist.size() < 4) {
ERR(FTL) << "Need 3 inputs for vector input in " << fName << std::endl;
ERR(FTL) << "Inputs: " << input << std::endl;
throw;
}
// Open File
std::string rootname = inputlist[0];
TFile* tempfile = new TFile(rootname.c_str(), "READ");
if (tempfile->IsZombie()) {
ERR(FTL) << "Can't find file in " << fName << std::endl;
ERR(FTL) << "location = " << rootname << std::endl;
throw;
}
// Get Name
std::string tagname = inputlist[1];
// TVector dialtags = tempfile->Get(tagname.c_str());
// if (!dialtags){
// ERR(FTL) << "Can't find list of dial names!" << std::endl;
// }
// Get Values
std::string valuename = inputlist[2];
TVectorD* dialvals = (TVectorD*)tempfile->Get(valuename.c_str());
if (!dialvals) {
ERR(FTL) << "Can't find dial values" << std::endl;
}
// Get Matrix
std::string matrixname = inputlist[3];
TMatrixD* matrixvals = (TMatrixD*)tempfile->Get(matrixname.c_str());
if (!matrixvals) {
ERR(FTL) << "Can't find matirx values" << std::endl;
}
// Get Types
if (inputlist.size() > 4) {
std::string typesname = inputlist[3];
}
// Get Minimum
if (inputlist.size() > 5) {
std::string minname = inputlist[4];
}
// Get Maximum
if (inputlist.size() > 6) {
std::string maxname = inputlist[5];
}
}
//*******************************************************************************
void ParamPull::ReadDialInput(std::string input) {
//*******************************************************************************
std::vector inputlist = GeneralUtils::ParseToStr(input, ";");
if (inputlist.size() < 3) {
ERR(FTL) << "Need 3 inputs for dial input in " << fName << std::endl;
ERR(FTL) << "Inputs: " << input << std::endl;
throw;
}
std::vector inputvals = GeneralUtils::ParseToDbl(input, ";");
std::string dialname = inputlist[0];
double val = inputvals[1];
double err = inputvals[2];
fDataHist = new TH1D( (fName + "_data").c_str(),
(fName + "_data" + fPlotTitles).c_str(), 1, 0, 1);
fDataHist->SetBinContent(1, val);
fDataHist->SetBinError(1, err);
fDataHist->GetXaxis()->SetBinLabel(1, dialname.c_str());
fLimitHist = new TH1D( (fName + "_limits").c_str(),
(fName + "_limits" + fPlotTitles).c_str(), 1, 0, 1);
fLimitHist->Reset();
if (inputvals.size() > 4) {
fLimitHist->SetBinContent(1, (inputvals[3] + inputvals[4]) / 2.0);
fLimitHist->SetBinError(1, (inputvals[4] - inputvals[3]) / 2.0);
}
fCovar = NULL;
}
//*******************************************************************************
std::map ParamPull::GetAllDials() {
//*******************************************************************************
-
std::map dialtypemap;
-
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
-
std::string name = fDataHist->GetXaxis()->GetBinLabel(i + 1);
int type = fTypeHist->GetBinContent(i + 1);
-
dialtypemap[name] = type;
-
}
-
return dialtypemap;
}
//*******************************************************************************
bool ParamPull::CheckDialsValid() {
//*******************************************************************************
- return true;
std::string helpstring = "";
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
std::string name = std::string(fDataHist->GetXaxis()->GetBinLabel(i + 1));
// If dial exists its all good
- if (FitBase::GetRW()->DialIncluded(name)) continue;
+ if (FitBase::GetRW()->DialIncluded(name)) {
+ LOG(DEB) << "Found dial " << name << " in covariance " << fInput << " and matched to reweight engine " << std::endl;
+ continue;
+ }
// If it doesn't but its a sample norm also continue
if (name.find("_norm") != std::string::npos) {
ERR(WRN) << "Norm dial included in covar but not set in FitWeight." << std::endl;
ERR(WRN) << "Assuming its a sample norm and skipping..." << std::endl;
}
// Dial unknown so print a help statement
std::ostringstream tempstr;
tempstr << "unknown_parameter " << name << " "
<< fDataHist->GetBinContent(i + 1) << " "
<< fDataHist->GetBinContent(i + 1) - fDataHist->GetBinError(i + 1) << " "
<< fDataHist->GetBinContent(i + 1) + fDataHist->GetBinError(i + 1) << " "
<< fDataHist->GetBinError(i + 1) << " ";
if (!fType.empty()) tempstr << fType << std::endl;
else tempstr << "FREE" << std::endl;
helpstring += tempstr.str();
}
// Show statement before failing
if (!helpstring.empty()) {
- ERR(WRN) << "Dial(s) included in covar but not set in FitWeight." << std::endl
- << "ParamPulls needs to know how you want it to be treated." << std::endl
- << "Include the following lines into your card:" << std::endl;
+ ERR(FTL) << "Dial(s) included in covar but not set in FitWeight." << std::endl;
+ ERR(FTL) << "ParamPulls needs to know how you want it to be treated." << std::endl;
+ ERR(FTL) << "Include the following lines into your card to throw UNCORRELATED:" << std::endl
- ERR(WRN) << helpstring << std::endl;
+ << helpstring << std::endl;
throw;
return false;
} else {
return true;
}
+ return false;
}
//*******************************************************************************
void ParamPull::Reconfigure() {
//*******************************************************************************
FitWeight* rw = FitBase::GetRW();
// Get Dial Names that are valid
std::vector namevec = rw->GetDialNames();
std::vector valuevec = rw->GetDialValues();
// Set Bin Values from RW
for (UInt_t i = 0; i < namevec.size(); i++) {
// Loop over bins and check name matches
std::string syst = namevec.at(i);
double systval = valuevec.at(i);
std::vector allsyst = GeneralUtils::ParseToStr(syst, ",");
// Proper Reconf using RW
for (int j = 0; j < fMCHist->GetNbinsX(); j++) {
// Search for the name of this bin in the corrent dial
std::string binname = std::string(fMCHist->GetXaxis()->GetBinLabel(j + 1) );
// Check Full Name
if (!syst.compare(binname.c_str())) {
fMCHist->SetBinContent(j + 1, systval);
break;
}
std::vector splitbinname = GeneralUtils::ParseToStr(binname, ",");
for (size_t l = 0; l < splitbinname.size(); l++) {
std::string singlebinname = splitbinname[l];
for (size_t k = 0; k < allsyst.size(); k++) {
if (!allsyst[k].compare(singlebinname.c_str())) {
fMCHist->SetBinContent(j + 1, systval);
}
}
}
}
}
return;
};
//*******************************************************************************
void ParamPull::ResetToy(void) {
//*******************************************************************************
if (fDataHist) delete fDataHist;
LOG(DEB) << "Resetting toy" << std::endl;
LOG(DEB) << fDataTrue << std::endl;
fDataHist = (TH1D*)fDataTrue->Clone();
LOG(DEB) << "Setting name" << std::endl;
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
}
//*******************************************************************************
void ParamPull::SetFakeData(std::string fakeinput) {
//*******************************************************************************
// Set from MC Setting
if (!fakeinput.compare("MC")) {
// Copy MC into data
if (fDataHist) delete fDataHist;
fDataHist = (TH1D*)fMCHist->Clone();
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " fakedata" + fPlotTitles).c_str() );
// Copy original data errors
for (int i = 0; i < fDataOrig->GetNbinsX(); i++) {
fDataHist->SetBinError(i + 1, fDataOrig->GetBinError(i + 1) );
}
// Make True Toy Central Value Hist
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
(fName + " truedata" + fPlotTitles).c_str() );
} else {
ERR(FTL) << "Trying to set fake data for ParamPulls not from MC!" << std::endl;
ERR(FTL) << "Not currently implemented.." << std::endl;
throw;
}
}
//*******************************************************************************
void ParamPull::RemoveFakeData() {
//*******************************************************************************
delete fDataHist;
fDataHist = (TH1D*)fDataOrig->Clone();
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " data" + fPlotTitles).c_str() );
fDataTrue = (TH1D*) fDataHist->Clone();
fDataTrue->SetNameTitle( (fName + "_truedata").c_str(),
(fName + " truedata" + fPlotTitles).c_str() );
}
//*******************************************************************************
double ParamPull::GetLikelihood() {
//*******************************************************************************
double like = 0.0;
switch (fCalcType) {
// Gaussian Calculation with correlations
case kGausPull:
like = StatUtils::GetChi2FromCov(fDataHist, fMCHist, fInvCovar, NULL);
like *= 1E-76;
break;
// Default says this has no pull
case kNoThrow:
default:
like = 0.0;
break;
}
LOG(DEB) << "Likelihood = " << like << " " << fCalcType << std::endl;
return like;
};
//*******************************************************************************
int ParamPull::GetNDOF() {
//*******************************************************************************
int ndof = 0;
if (fCalcType != kNoThrow) {
ndof = fDataHist->GetNbinsX();
}
return ndof;
};
//*******************************************************************************
void ParamPull::ThrowCovariance() {
//*******************************************************************************
// Reset toy for throw
ResetToy();
LOG(FIT) << "Creating new toy dataset" << std::endl;
// Generate random Gaussian throws
std::vector randthrows;
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
double randtemp = 0.0;
switch (fThrowType) {
// Gaussian Throws
case kGausThrow:
randtemp = gRandom->Gaus(0.0, 1.0);
break;
// Uniform Throws
case kFlatThrow:
randtemp = gRandom->Uniform(0.0, 1.0);
if (fLimitHist) {
randtemp = fLimitHist->GetBinContent(i + 1) + fLimitHist->GetBinError(i + 1) * ( randtemp * 2 - 1 );
}
break;
// No Throws (DEFAULT)
default:
break;
}
randthrows.push_back(randtemp);
}
// Create Bin Modifications
double totalres = 0.0;
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
// Calc Bin Mod
double binmod = 0.0;
if (fThrowType == kGausThrow) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
binmod += (*fDecomp)(j, i) * randthrows.at(j);
}
} else if (fThrowType == kFlatThrow) {
binmod = randthrows.at(i) - fDataHist->GetBinContent(i + 1);
}
-
// Add up fraction dif
totalres += binmod;
// Add to current data
fDataHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) + binmod);
}
// Rename
fDataHist->SetNameTitle( (fName + "_data").c_str(),
(fName + " toydata" + fPlotTitles).c_str() );
// Check Limits
if (fLimitHist) {
for (int i = 0; i < fLimitHist->GetNbinsX(); i++) {
if (fLimitHist->GetBinError(i + 1) == 0.0) continue;
if (fDataHist->GetBinContent(i + 1) > fLimitHist->GetBinContent(i + 1) + fLimitHist->GetBinError(i + 1) ||
fDataHist->GetBinContent(i + 1) < fLimitHist->GetBinContent(i + 1) - fLimitHist->GetBinError(i + 1)) {
- this->ThrowCovariance();
+ LOG(FIT) << "Threw outside allowed region, rethrowing..." << std::endl;
+ ThrowCovariance();
}
}
}
-
- return;
};
//*******************************************************************************
TH2D ParamPull::GetCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fInvCovar);
for (int i = 0; i < tempCov.GetNbinsX(); i++) {
tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_INVCOV").c_str(),
(fName + " InvertedCovariance;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
TH2D ParamPull::GetFullCovar() {
//*******************************************************************************
TH2D tempCov = TH2D(*fCovar);
for (int i = 0; i < tempCov.GetNbinsX(); i++) {
tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_COV").c_str(),
(fName + " Covariance;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
TH2D ParamPull::GetDecompCovar() {
//*******************************************************************************
- TH2D tempCov = TH2D(*fCovar);
+ TH2D tempCov = TH2D(*fDecomp);
for (int i = 0; i < tempCov.GetNbinsX(); i++) {
tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1));
}
tempCov.SetNameTitle((fName + "_DEC").c_str(),
(fName + " Decomposition;Dials;Dials").c_str());
return tempCov;
}
//*******************************************************************************
void ParamPull::Write(std::string writeoptt) {
//*******************************************************************************
fDataHist->Write();
fMCHist->Write();
if (fLimitHist) {
fLimitHist->Write();
}
GetCovar().Write();
GetFullCovar().Write();
GetDecompCovar().Write();
return;
};
+void ParamPull::CheckHist(TH1D* hist) {
+ if (!hist) {
+ ERR(FTL) << "Can't find TH1D hist fit_dials in " << fName << std::endl;
+ ERR(FTL) << "File Entries:" << std::endl;
+ TFile *temp = new TFile(fInput.c_str(), "open");
+ temp->ls();
+ throw;
+ }
+}
diff --git a/src/FitBase/ParamPull.h b/src/FitBase/ParamPull.h
index aa99977..e62acb5 100644
--- a/src/FitBase/ParamPull.h
+++ b/src/FitBase/ParamPull.h
@@ -1,194 +1,195 @@
// 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 .
*******************************************************************************/
#ifndef PARAM_PULL_H_SEEN
#define PARAM_PULL_H_SEEN
/*!
* \addtogroup FitBase
* @{
*/
#include
#include
#include
#include
#include
#include
#include
#include
// ROOT includes
#include
#include
#include
#include
#include
#include
#include
#include "TRandom3.h"
// Fit Includes
#include "PlotUtils.h"
#include "StatUtils.h"
#include "FitWeight.h"
#include "FitLogger.h"
#include "EventManager.h"
#include "TVector.h"
using namespace std;
//! Enums to allow Non Guassian Pulls in Future
enum FitPullTypes {
kUnknownPull = -1,
kNoPull = 0,
kGausPull = 1
};
enum FitThrowTypes {
kUnknownThrow = -1,
kNoThrow = 0,
kGausThrow = 1,
kFlatThrow = 2
};
//! Used to produce gaussian penalty terms in the fit.
class ParamPull {
public:
//! Default Constructor
ParamPull(std::string name, std::string inputfile, std::string type, std::string dials="");
//! Default destructor
virtual ~ParamPull(void) {};
// Set dial types (DEFAULT,ABS,FRAC)
void SetType(std::string type);
// Setup Histogram inputs (from previous fit file, or ROOT file)
void SetupHistograms(std::string input);
// Read a previous NUISANCE file
void ReadFitFile(std::string input);
// Read a ROOT file with any histograms in
void ReadRootFile(std::string input);
// Read Text file
void ReadVectFile(std::string input);
// Read a single dial central value and error
void ReadDialInput(std::string input);
//! Reconfigure function reads in current RW engine dials and sets their value to MC
void Reconfigure(void);
//! Get likelihood given the current values
double GetLikelihood(void);
//! Get NDOF if used in likelihoods
int GetNDOF(void);
// Get Covariance Matrices as plots
TH2D GetCovar(void);
TH2D GetFullCovar(void);
TH2D GetDecompCovar(void);
// Get Covariance Matrices
inline TMatrixDSym GetCovarMatrix (void) const { return *fInvCovar; };
inline TMatrixDSym GetFullCovarMatrix (void) const { return *fCovar; };
inline TMatrixDSym GetDecompCovarMatrix (void) const { return *fDecomp; };
//! Save the histograms
void Write(std::string writeopt="");
//! Throw the dial values using the current covariance. Useful for parameter throws.
void ThrowCovariance(void);
//! Compare dials to RW
bool CheckDialsValid(void);
//! Reset toy data back to original data
void ResetToy(void);
//! Read fake data from MC
void SetFakeData(std::string fakeinput);
//! Reset fake data back to original data (before fake data or throws)
void RemoveFakeData();
// Get Functions
inline std::string GetName (void) const { return fName; };
inline std::string GetInput (void) const { return fInput; };
inline std::string GetType (void) const { return fType; };
inline std::string GetFileType (void) const { return fFileType; };
inline std::string GetDialOptions (void) const { return fDialOptions; };
std::map GetAllDials();
inline TH1D GetDataHist (void) const { return *fDataHist; };
inline TH1D GetDataTrue (void) const { return *fDataTrue; };
inline TH1D GetDataOrig (void) const { return *fDataOrig; };
inline TH1D GetMCHist (void) const { return *fMCHist; };
inline TH1D GetMaxHist (void) const { return *fMaxHist; };
inline TH1D GetMinHist (void) const { return *fMinHist; };
inline TH1I GetDialTypes (void) const { return *fTypeHist; };
inline TH1D GetLimitHist (void) const { return *fLimitHist; };
private:
+ void CheckHist(TH1D*);
TH1D RemoveBinsNotInString(TH1D hist, std::string mystr);
TH1I RemoveBinsNotInString(TH1I hist, std::string mystr);
std::string fName; //!< Pull Name
std::string fInput; //!< Pull input string
std::string fType; //!< Pull options type
std::string fFileType; //!< Pull input file types
std::string fPlotTitles; //! Axis format
std::string fDialOptions; //!< Dial handling options
std::string fDialSelection; //!< Dial Selection
TH1D* fMCHist; //!< Current MC Histogram
TH1D* fDataHist; //!< Current data Histogram
TH1D* fDataTrue; //!< True Data (before histogram throws)
TH1D* fDataOrig; //!< Orig Data (without toys or fake data)
TH1D* fMaxHist; //!< Maximum limit on MC/Data
TH1D* fMinHist; //!< Maximum limit on MC/Data
TH1I* fTypeHist; //!< Dial Types
int fCalcType; //!< Method to calculate likelihood
int fThrowType; //!< Method to calculate throws
TMatrixDSym* fCovar; //!< Covariance
TMatrixDSym* fInvCovar; //!< Inverted Covariance
TMatrixDSym* fDecomp; //!< Decomposition
TH1D* fLimitHist;
};
// Class TypeDefs
typedef std::list::const_iterator PullListConstIter;
typedef std::list::iterator PullListIter;
typedef std::vector::const_iterator PullVectConstIter;
typedef std::vector::iterator PullVectIter;
/*! @} */
#endif
diff --git a/src/Routines/ComparisonRoutines.cxx b/src/Routines/ComparisonRoutines.cxx
index 48b8abb..011e70d 100755
--- a/src/Routines/ComparisonRoutines.cxx
+++ b/src/Routines/ComparisonRoutines.cxx
@@ -1,536 +1,536 @@
// 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 .
*******************************************************************************/
#include "ComparisonRoutines.h"
/*
Constructor/Destructor
*/
//************************
void ComparisonRoutines::Init() {
//************************
fOutputFile = "";
fOutputRootFile = NULL;
fStrategy = "Compare";
fRoutines.clear();
fCardFile = "";
fFakeDataInput = "";
fSampleFCN = NULL;
fAllowedRoutines = ("Compare");
};
//*************************************
ComparisonRoutines::~ComparisonRoutines() {
//*************************************
delete fOutputRootFile;
};
/*
Input Functions
*/
//*************************************
ComparisonRoutines::ComparisonRoutines(int argc, char* argv[]) {
//*************************************
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector xmlcmds;
std::vector configargs;
// Make easier to handle arguments.
std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
nuiskey fCompKey;
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.Set("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.Set("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadSettings( fCompKey.GetS("cardfile"), "");
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")){
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::GetParI("VERBOSITY");
errorcount += Config::GetParI("ERROR");
bool trace = Config::GetParB("TRACE");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
SETVERBOSITY(verbocount);
SETTRACE(trace);
// Comparison Setup ========================================
// Proper Setup
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
SetupComparisonsFromXML();
SetupRWEngine();
SetupFCN();
return;
};
//*************************************
void ComparisonRoutines::SetupComparisonsFromXML() {
//*************************************
LOG(FIT) << "Setting up nuiscomp" << std::endl;
// Setup Parameters ------------------------------------------
std::vector parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
}
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
ERR(FTL) << "type='PARAMETER_TYPE'" << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
ERR(FTL) << "name='SAMPLE_NAME'" << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
ERR(FTL) << "nominal='NOMINAL_VALUE'" << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
double parlow = parnom - 1;
double parhigh = parnom + 1;
double parstep = 1;
// override if state not given
if (!key.Has("state")){
key.SetS("state","FIX");
}
std::string parstate = key.GetS("state");
// Check for incomplete limtis
int limdef = ((int)key.Has("low") +
(int)key.Has("high") +
(int)key.Has("step"));
if (limdef > 0 and limdef < 3){
ERR(FTL) << "Incomplete limit set given for parameter : " << parname << std::endl;
ERR(FTL) << "Requires: low='LOWER_LIMIT' high='UPPER_LIMIT' step='STEP_SIZE' " << std::endl;
throw;
}
// Extra limits
if (key.Has("low")) {
parlow = key.GetD("low");
parhigh = key.GetD("high");
parstep = key.GetD("step");
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parlow << " < p < " << parhigh
<< " : " << parstate << std::endl;
} else {
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parstate << std::endl;
}
// Convert if required
if (parstate.find("ABS") != std::string::npos) {
parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
} else if (parstate.find("FRAC") != std::string::npos) {
parnom = FitBase::RWFracToSigma( partype, parname, parnom );
parlow = FitBase::RWFracToSigma( partype, parname, parlow );
parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
parstep = FitBase::RWFracToSigma( partype, parname, parstep );
}
// Push into vectors
fParams.push_back(parname);
- fTypeVals[parname] = FitBase::ConvDialType(partype);;
+ fTypeVals[parname] = FitBase::ConvDialType(partype);
fCurVals[parname] = parnom;
fStateVals[parname] = parstate;
}
// Setup Samples ----------------------------------------------
std::vector samplekeys = Config::QueryKeys("sample");
if (!samplekeys.empty()) {
LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
}
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys.at(i);
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype =
key.Has("type") ? key.GetS("type") : "DEFAULT";
double samplenorm =
key.Has("norm") ? key.GetD("norm") : 1.0;
// Print out
LOG(FIT) << "Read Sample " << i << ". : "
<< samplename << " (" << sampletype << ") [Norm=" << samplenorm<<"]"<< std::endl
<< " -> input='" << samplefile << "'" << std::endl;
// If FREE add to parameters otherwise continue
if (sampletype.find("FREE") == std::string::npos) {
if (samplenorm != 1.0) {
ERR(FTL) << "You provided a sample normalisation but did not specify that the sample is free" << std::endl;
ERR(FTL) << "Change so sample contains type=\"FREE\" and re-run" << std::endl;
throw;
}
continue;
}
// Form norm dial from samplename + sampletype + "_norm";
std::string normname = samplename + "_norm";
// Check normname not already present
if (fTypeVals.find("normname") != fTypeVals.end()) {
continue;
}
// Add new norm dial to list if its passed above checks
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStateVals[normname] = sampletype;
fCurVals[normname] = samplenorm;
}
// Setup Fake Parameters -----------------------------
std::vector fakekeys = Config::QueryKeys("fakeparameter");
if (!fakekeys.empty()) {
LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
}
for (size_t i = 0; i < fakekeys.size(); i++) {
nuiskey key = fakekeys.at(i);
// Check for type,name,nom
if (!key.Has("name")) {
ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
throw;
}
// Get Inputs
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
// Push into vectors
fFakeVals[parname] = parnom;
}
}
//*************************************
void ComparisonRoutines::SetupRWEngine() {
//*************************************
LOG(FIT) << "Setting up FitWeight Engine" << std::endl;
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name));
}
return;
}
//*************************************
void ComparisonRoutines::SetupFCN() {
//*************************************
LOG(FIT) << "Building the SampleFCN" << std::endl;
if (fSampleFCN) delete fSampleFCN;
Config::Get().out = fOutputRootFile;
fOutputRootFile->cd();
fSampleFCN = new JointFCN(fOutputRootFile);
SetFakeData();
return;
}
//*************************************
void ComparisonRoutines::SetFakeData() {
//*************************************
if (fFakeDataInput.empty()) return;
if (fFakeDataInput.compare("MC") == 0) {
LOG(FIT) << "Setting fake data from MC starting prediction." << std::endl;
UpdateRWEngine(fFakeVals);
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->SetFakeData("MC");
UpdateRWEngine(fCurVals);
LOG(FIT) << "Set all data to fake MC predictions." << std::endl;
} else {
LOG(FIT) << "Setting fake data from: " << fFakeDataInput << std::endl;
fSampleFCN->SetFakeData(fFakeDataInput);
}
return;
}
/*
Fitting Functions
*/
//*************************************
void ComparisonRoutines::UpdateRWEngine(
std::map& updateVals) {
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
FitBase::GetRW()->SetDialValue(name, updateVals.at(name));
}
FitBase::GetRW()->Reconfigure();
return;
}
//*************************************
void ComparisonRoutines::Run() {
//*************************************
LOG(FIT) << "Running ComparisonRoutines : " << fStrategy << std::endl;
if (FitPar::Config().GetParB("save_nominal")) {
SaveNominal();
}
// Parse given routines
fRoutines = GeneralUtils::ParseToStr(fStrategy, ",");
if (fRoutines.empty()) {
ERR(FTL) << "Trying to run ComparisonRoutines with no routines given!" << std::endl;
throw;
}
for (UInt_t i = 0; i < fRoutines.size(); i++) {
std::string routine = fRoutines.at(i);
LOG(FIT) << "Routine: " << routine << std::endl;
if (!routine.compare("Compare")) {
UpdateRWEngine(fCurVals);
GenerateComparison();
PrintState();
SaveCurrentState();
}
}
return;
}
//*************************************
void ComparisonRoutines::GenerateComparison() {
//*************************************
LOG(FIT) << "Generating Comparison." << std::endl;
// Main Event Loop from event Manager
fSampleFCN->ReconfigureAllEvents();
return;
}
//*************************************
void ComparisonRoutines::PrintState() {
//*************************************
LOG(FIT) << "------------" << std::endl;
// Count max size
int maxcount = 0;
for (UInt_t i = 0; i < fParams.size(); i++) {
maxcount = max(int(fParams[i].size()), maxcount);
}
// Header
LOG(FIT) << " # " << left << setw(maxcount) << "Parameter "
<< " = " << setw(10) << "Value"
<< " +- " << setw(10) << "Error"
<< " " << setw(8) << "(Units)"
<< " " << setw(10) << "Conv. Val"
<< " +- " << setw(10) << "Conv. Err"
<< " " << setw(8) << "(Units)" << std::endl;
// Parameters
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
std::string curunits = "(sig.)";
double curval = fCurVals[syst];
double curerr = 0.0;
if (fStateVals[syst].find("ABS") != std::string::npos) {
curval = FitBase::RWSigmaToAbs(typestr, syst, curval);
curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
curunits = "(Abs.)";
} else if (fStateVals[syst].find("FRAC") != std::string::npos) {
curval = FitBase::RWSigmaToFrac(typestr, syst, curval);
curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) -
FitBase::RWSigmaToFrac(typestr, syst, 0.0));
curunits = "(Frac)";
}
std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")";
double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
std::ostringstream curparstring;
curparstring << " " << setw(3) << left << i << ". " << setw(maxcount)
<< syst << " = " << setw(10) << curval << " +- " << setw(10)
<< curerr << " " << setw(8) << curunits << " " << setw(10)
<< convval << " +- " << setw(10) << converr << " " << setw(8)
<< convunits;
LOG(FIT) << curparstring.str() << std::endl;
}
LOG(FIT) << "------------" << std::endl;
double like = fSampleFCN->GetLikelihood();
LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl;
LOG(FIT) << "------------" << std::endl;
}
/*
Write Functions
*/
//*************************************
void ComparisonRoutines::SaveCurrentState(std::string subdir) {
//*************************************
LOG(FIT) << "Saving current full FCN predictions" << std::endl;
// Setup DIRS
TDirectory* curdir = gDirectory;
if (!subdir.empty()) {
TDirectory* newdir = (TDirectory*)gDirectory->mkdir(subdir.c_str());
newdir->cd();
}
fSampleFCN->Write();
// Change back to current DIR
curdir->cd();
return;
}
//*************************************
void ComparisonRoutines::SaveNominal() {
//*************************************
fOutputRootFile->cd();
LOG(FIT) << "Saving Nominal Predictions (be cautious with this)" << std::endl;
FitBase::GetRW()->Reconfigure();
GenerateComparison();
SaveCurrentState("nominal");
};
diff --git a/src/Routines/SystematicRoutines.cxx b/src/Routines/SystematicRoutines.cxx
index 49988e6..04c98f6 100755
--- a/src/Routines/SystematicRoutines.cxx
+++ b/src/Routines/SystematicRoutines.cxx
@@ -1,1457 +1,1437 @@
// 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 .
*******************************************************************************/
#include "SystematicRoutines.h"
void SystematicRoutines::Init(){
fInputFile = "";
fInputRootFile = NULL;
fOutputFile = "";
fOutputRootFile = NULL;
fCovar = fCovarFree = NULL;
fCorrel = fCorrelFree = NULL;
fDecomp = fDecompFree = NULL;
fStrategy = "ErrorBands";
fRoutines.clear();
fRoutines.push_back("ErrorBands");
fCardFile = "";
fFakeDataInput = "";
fSampleFCN = NULL;
fAllowedRoutines = ("ErrorBands,PlotLimits");
};
SystematicRoutines::~SystematicRoutines(){
};
SystematicRoutines::SystematicRoutines(int argc, char* argv[]){
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector xmlcmds;
std::vector configargs;
fNThrows = 250;
fStartThrows = 0;
fThrowString = "";
// Make easier to handle arguments.
std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false);
ParserUtils::ParseArgument(args, "-s", fStartThrows, false, false);
ParserUtils::ParseArgument(args, "-t", fNThrows, false, false);
ParserUtils::ParseArgument(args, "-p", fThrowString, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.Set("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.Set("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadSettings( fCompKey.GetS("cardfile"), "");
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")){
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::GetParI("VERBOSITY");
errorcount += Config::GetParI("ERROR");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
SETVERBOSITY(verbocount);
// Proper Setup
if (fStrategy.find("ErrorBands") != std::string::npos ||
fStrategy.find("MergeErrors") != std::string::npos){
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
}
// fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
SetupSystematicsFromXML();
SetupCovariance();
SetupRWEngine();
SetupFCN();
GetCovarFromFCN();
// Run();
return;
};
void SystematicRoutines::SetupSystematicsFromXML(){
LOG(FIT) << "Setting up nuismin" << std::endl;
// Setup Parameters ------------------------------------------
std::vector parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
- }
+ }
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
double parlow = parnom - 1;
double parhigh = parnom + 1;
double parstep = 1;
-
// Override if state not given
if (!key.Has("state")){
key.SetS("state","FIX");
}
std::string parstate = key.GetS("state");
// Extra limits
if (key.Has("low")) {
parlow = key.GetD("low");
parhigh = key.GetD("high");
parstep = key.GetD("step");
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parlow << " < p < " << parhigh
<< " : " << parstate << std::endl;
} else {
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parstate << std::endl;
}
// Run Parameter Conversion if needed
if (parstate.find("ABS") != std::string::npos) {
parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
} else if (parstate.find("FRAC") != std::string::npos) {
parnom = FitBase::RWFracToSigma( partype, parname, parnom );
parlow = FitBase::RWFracToSigma( partype, parname, parlow );
parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
parstep = FitBase::RWFracToSigma( partype, parname, parstep );
}
// Push into vectors
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);;
fStartVals[parname] = parnom;
fCurVals[parname] = parnom;
fErrorVals[parname] = 0.0;
fStateVals[parname] = parstate;
bool fixstate = parstate.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = parlow;
fMaxVals[parname] = parhigh;
fStepVals[parname] = parstep;
}
// Setup Samples ----------------------------------------------
std::vector samplekeys = Config::QueryKeys("sample");
if (!samplekeys.empty()) {
LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
}
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys.at(i);
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype =
key.Has("type") ? key.GetS("type") : "DEFAULT";
double samplenorm =
key.Has("norm") ? key.GetD("norm") : 1.0;
// Print out
LOG(FIT) << "Read sample info " << i << " : "
<< samplename << std::endl
<< "\t\t input -> " << samplefile << std::endl
<< "\t\t state -> " << sampletype << std::endl
<< "\t\t norm -> " << samplenorm << std::endl;
// If FREE add to parameters otherwise continue
if (sampletype.find("FREE") == std::string::npos) {
continue;
}
// Form norm dial from samplename + sampletype + "_norm";
std::string normname = samplename + "_norm";
// Check normname not already present
if (fTypeVals.find(normname) != fTypeVals.end()) {
continue;
}
// Add new norm dial to list if its passed above checks
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStateVals[normname] = sampletype;
fCurVals[normname] = samplenorm;
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = sampletype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
}
// Setup Fake Parameters -----------------------------
std::vector fakekeys = Config::QueryKeys("fakeparameter");
if (!fakekeys.empty()) {
LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
}
for (size_t i = 0; i < fakekeys.size(); i++) {
nuiskey key = fakekeys.at(i);
// Check for type,name,nom
if (!key.Has("name")) {
ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
throw;
} else if (!key.Has("nom")) {
ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
throw;
}
// Get Inputs
std::string parname = key.GetS("name");
double parnom = key.GetD("nom");
// Push into vectors
fFakeVals[parname] = parnom;
}
+
}
/*
Setup Functions
*/
//*************************************
void SystematicRoutines::SetupRWEngine(){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
FitBase::GetRW() -> IncludeDial(name, fTypeVals.at(name) );
}
UpdateRWEngine(fStartVals);
return;
}
//*************************************
void SystematicRoutines::SetupFCN(){
//*************************************
LOG(FIT)<<"Making the jointFCN"<Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->SetFakeData("MC");
UpdateRWEngine(fCurVals);
LOG(FIT)<<"Set all data to fake MC predictions."<SetFakeData(fFakeDataInput);
}
return;
}
//*****************************************
void SystematicRoutines::GetCovarFromFCN(){
//*****************************************
- LOG(FIT) << "Loading ParamPull objects from FCN to build covar" << std::endl;
+ LOG(FIT) << "Loading ParamPull objects from FCN to build covariance..." << std::endl;
// Make helperstring
std::ostringstream helperstr;
// Keep track of what is being thrown
std::map dialthrowhandle;
// Get Covariance Objects from FCN
std::list inputpulls = fSampleFCN->GetPullList();
- for (PullListConstIter iter = inputpulls.begin();
- iter != inputpulls.end(); iter++){
+ for (PullListConstIter iter = inputpulls.begin(); iter != inputpulls.end(); iter++){
ParamPull* pull = (*iter);
-
- if (pull->GetType().find("THROW")){
+ if (pull->GetType().find("THROW") != std::string::npos){
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
LOG(FIT) << "Read ParamPull: " << pull->GetName() << " " << pull->GetType() << std::endl;
}
TH1D dialhist = pull->GetDataHist();
TH1D minhist = pull->GetMinHist();
TH1D maxhist = pull->GetMaxHist();
TH1I typehist = pull->GetDialTypes();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
dialthrowhandle[name] = pull->GetName();
if (fCurVals.find(name) == fCurVals.end()){
// Add to Containers
fParams.push_back(name);
fCurVals[name] = dialhist.GetBinContent(i+1);
fStartVals[name] = dialhist.GetBinContent(i+1);
fMinVals[name] = minhist.GetBinContent(i+1);
fMaxVals[name] = maxhist.GetBinContent(i+1);
fStepVals[name] = 1.0;
fFixVals[name] = false;
fStartFixVals[name] = false;
fTypeVals[name] = typehist.GetBinContent(i+1);
fStateVals[name] = "FREE" + pull->GetType();
// Maker Helper
helperstr << std::string(16, ' ' ) << FitBase::ConvDialType(fTypeVals[name]) << " "
<< name << " " << fMinVals[name] << " "
<< fMaxVals[name] << " " << fStepVals[name] << " " << fStateVals[name]
<< std::endl;
}
}
}
// Check if no throws given
if (fInputThrows.empty()){
ERR(WRN) << "No covariances given to nuissyst" << std::endl;
ERR(WRN) << "Pushing back an uncorrelated gaussian throw error for each free parameter using step size" << std::endl;
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (fFixVals[syst]) continue;
// Make Terms
std::string name = syst + "_pull";
std::ostringstream pullterm;
pullterm << "DIAL:" << syst << ";"
<< fStartVals[syst] << ";"
<< fStepVals[syst];
std::string type = "GAUSTHROW/NEUT";
// Push Back Pulls
ParamPull* pull = new ParamPull( name, pullterm.str(), type );
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
// Print Whats added
- ERR(WRN) << "Added ParamPull : " << name << " " << pullterm.str() << " " << type << std::endl;
+ LOG(FIT) << "Added ParamPull : " << name << " " << pullterm.str() << " " << type << std::endl;
// Add helper string for future fits
helperstr << std::string(16, ' ' ) << "covar " << name << " " << pullterm.str() << " " << type << std::endl;
// Keep Track of Throws
dialthrowhandle[syst] = pull->GetName();
}
}
// Print Helper String
if (!helperstr.str().empty()){
LOG(FIT) << "To remove these statements in future studies, add the lines below to your card:" << std::endl;
// Can't use the logger properly because this can be multi-line. Use cout and added spaces to look better!
std::cout << helperstr.str();
sleep(2);
}
// Print Throw State
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (dialthrowhandle.find(syst) != dialthrowhandle.end()){
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = THROWING with " << dialthrowhandle[syst] << std::endl;
} else {
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = FIXED" << std::endl;
}
}
// Pause anyway
sleep(1);
return;
}
/*
Fitting Functions
*/
//*************************************
void SystematicRoutines::UpdateRWEngine(std::map& updateVals){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
FitBase::GetRW()->SetDialValue(name,updateVals.at(name));
}
FitBase::GetRW()->Reconfigure();
return;
}
//*************************************
void SystematicRoutines::PrintState(){
//*************************************
LOG(FIT)<<"------------"<GetLikelihood();
LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl;
LOG(FIT)<<"------------"<cd();
SaveCurrentState();
}
//*************************************
void SystematicRoutines::SaveCurrentState(std::string subdir){
//*************************************
LOG(FIT)<<"Saving current full FCN predictions" <mkdir(subdir.c_str());
newdir->cd();
}
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
// Change back to current DIR
curdir->cd();
return;
}
//*************************************
void SystematicRoutines::SaveNominal(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Nominal Predictions (be cautious with this)" <Reconfigure();
SaveCurrentState("nominal");
};
//*************************************
void SystematicRoutines::SavePrefit(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Prefit Predictions"< 0){
fCovarFree = new TH2D("covariance_free",
"covariance_free",
NFREE,0,NFREE,
NFREE,0,NFREE);
}
// Set Bin Labels
int countall = 0;
int countfree = 0;
for (UInt_t i = 0; i < fParams.size(); i++){
fCovar->GetXaxis()->SetBinLabel(countall+1,fParams[i].c_str());
fCovar->GetYaxis()->SetBinLabel(countall+1,fParams[i].c_str());
countall++;
if (!fFixVals[fParams[i]] and NFREE > 0){
fCovarFree->GetXaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
fCovarFree->GetYaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
countfree++;
}
}
fCorrel = PlotUtils::GetCorrelationPlot(fCovar,"correlation");
fDecomp = PlotUtils::GetDecompPlot(fCovar,"decomposition");
if (NFREE > 0)fCorrelFree = PlotUtils::GetCorrelationPlot(fCovarFree, "correlation_free");
if (NFREE > 0)fDecompFree = PlotUtils::GetDecompPlot(fCovarFree,"decomposition_free");
return;
};
//*************************************
void SystematicRoutines::ThrowCovariance(bool uniformly){
//*************************************
// Set fThrownVals to all values in currentVals
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams.at(i);
fThrownVals[name] = fCurVals[name];
}
for (PullListConstIter iter = fInputThrows.begin();
iter != fInputThrows.end(); iter++){
ParamPull* pull = *iter;
pull->ThrowCovariance();
TH1D dialhist = pull->GetDataHist();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
if (fCurVals.find(name) != fCurVals.end()){
fThrownVals[name] = dialhist.GetBinContent(i+1);
}
}
- // Reset throw incase pulls are calculated.
+ // Reset throw in case pulls are calculated.
pull->ResetToy();
-
}
-
- return;
};
//*************************************
void SystematicRoutines::PlotLimits(){
//*************************************
std::cout << "Plotting Limits" << std::endl;
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
TDirectory* limfolder = (TDirectory*) fOutputRootFile->mkdir("Limits");
limfolder->cd();
// Set all parameters at their starting values
for (UInt_t i = 0; i < fParams.size(); i++){
fCurVals[fParams[i]] = fStartVals[fParams[i]];
}
TDirectory* nomfolder = (TDirectory*) limfolder->mkdir("nominal");
nomfolder->cd();
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
limfolder->cd();
std::vector allfolders;
// Loop through each parameter
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
std::cout << "Starting Param " << syst << std::endl;
if (fFixVals[syst]) continue;
// Loop Downwards
while (fCurVals[syst] > fMinVals[syst]){
fCurVals[syst] = fCurVals[syst] - fStepVals[syst];
// Check Limit
if (fCurVals[syst] < fMinVals[syst])
fCurVals[syst] = fMinVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder for variation
TDirectory* minfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
minfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to folder
fSampleFCN->Write();
}
// Reset before next loop
fCurVals[syst] = fStartVals[syst];
// Loop Upwards now
while (fCurVals[syst] < fMaxVals[syst]){
fCurVals[syst] = fCurVals[syst] + fStepVals[syst];
// Check Limit
if (fCurVals[syst] > fMaxVals[syst])
fCurVals[syst] = fMaxVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder
TDirectory* maxfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
maxfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to file
fSampleFCN->Write();
}
// Reset before leaving
fCurVals[syst] = fStartVals[syst];
UpdateRWEngine(fCurVals);
}
return;
}
//*************************************
void SystematicRoutines::Run(){
//*************************************
std::cout << "Running routines "<< std::endl;
fRoutines = GeneralUtils::ParseToStr(fStrategy,",");
for (UInt_t i = 0; i < fRoutines.size(); i++){
std::string routine = fRoutines.at(i);
int fitstate = kFitUnfinished;
LOG(FIT)<<"Running Routine: "<cd();
- int nthrows = fNThrows;
+ // For generating throws we check with the config
+ int nthrows = Config::GetParI("error_throws");
int startthrows = fStartThrows;
int endthrows = startthrows + nthrows;
if (nthrows < 0) nthrows = endthrows;
if (startthrows < 0) startthrows = 0;
if (endthrows < 0) endthrows = startthrows + nthrows;
// Setting Seed
// Matteo Mazzanti's Fix
struct timeval mytime;
gettimeofday(&mytime, NULL);
Double_t seed = time(NULL) + int(getpid())+ (mytime.tv_sec * 1000.) + (mytime.tv_usec / 1000.);
gRandom->SetSeed(seed);
// int seed = (gRandom->Uniform(0.0,1.0)*100000 + 100000000*(startthrows + endthrows) + time(NULL) + int(getpid()) );
// gRandom->SetSeed(seed);
LOG(FIT) << "Using Seed : " << seed << std::endl;
LOG(FIT) << "nthrows = " << nthrows << std::endl;
LOG(FIT) << "startthrows = " << startthrows << std::endl;
LOG(FIT) << "endthrows = " << endthrows << std::endl;
-
-
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
if (startthrows == 0){
LOG(FIT) << "Making nominal " << std::endl;
TDirectory* nominal = (TDirectory*) tempfile->mkdir("nominal");
nominal->cd();
fSampleFCN->Write();
}
- LOG(SAM) << "nthrows = " << nthrows << std::endl;
- LOG(SAM) << "startthrows = " << startthrows << std::endl;
- LOG(SAM) << "endthrows = " << endthrows << std::endl;
-
TTree* parameterTree = new TTree("throws","throws");
double chi2;
for (UInt_t i = 0; i < fParams.size(); i++)
parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]], (fParams[i] + "/D").c_str());
parameterTree->Branch("chi2",&chi2,"chi2/D");
fSampleFCN->CreateIterationTree("error_iterations", FitBase::GetRW());
// Would anybody actually want to do uniform throws of any parameter??
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Run Throws and save
for (Int_t i = 0; i < endthrows+1; i++){
- LOG(FIT) << "Loop " << i << std::endl;
ThrowCovariance(uniformly);
if (i < startthrows) continue;
if (i == 0) continue;
- LOG(FIT) << "Throw " << i << " ================================" << std::endl;
+ LOG(FIT) << "Throw " << i << "/" << endthrows << " ================================" << std::endl;
// Generate Random Parameter Throw
// ThrowCovariance(uniformly);
TDirectory* throwfolder = (TDirectory*)tempfile->mkdir(Form("throw_%i",i));
throwfolder->cd();
// Run Eval
double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals );
chi2 = fSampleFCN->DoEval( vals );
delete vals;
// Save the FCN
fSampleFCN->Write();
parameterTree->Fill();
}
tempfile->cd();
fSampleFCN->WriteIterationTree();
tempfile->Close();
}
-void SystematicRoutines::MergeThrows(){
-
-
- fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
- fOutputRootFile->cd();
+void SystematicRoutines::MergeThrows() {
+ fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
+ fOutputRootFile->cd();
// Make a container folder
TDirectory* errorDIR = (TDirectory*) fOutputRootFile->mkdir("error_bands");
errorDIR->cd();
TDirectory* outnominal = (TDirectory*) fOutputRootFile->mkdir("nominal_throw");
outnominal->cd();
// Split Input Files
if (!fThrowString.empty()) fThrowList = GeneralUtils::ParseToStr(fThrowString,",");
// Add default if no throwlist given
if (fThrowList.size() < 1) fThrowList.push_back( fOutputFile + ".throws.root" );
/// Save location of file containing nominal
std::string nominalfile;
bool nominalfound;
// Loop over files and check they exist.
for (uint i = 0; i < fThrowList.size(); i++){
std::string file = fThrowList[i];
bool found = false;
// normal
std::string newfile = file;
TFile* throwfile = new TFile(file.c_str(),"READ");
if (throwfile and !throwfile->IsZombie()){
- found = true;
+ found = true;
}
// normal.throws.root
if (!found){
newfile = file + ".throws.root";
throwfile = new TFile((file + ".throws.root").c_str(),"READ");
if (throwfile and !throwfile->IsZombie()) {
found = true;
}
}
// If its found save to throwlist, else save empty.
// Also search for nominal
if (found){
fThrowList[i] = newfile;
-
+
LOG(FIT) << "Throws File :" << newfile << std::endl;
// Find input which contains nominal
if (throwfile->Get("nominal")){
nominalfound = true;
nominalfile = newfile;
}
-
throwfile->Close();
-
} else {
fThrowList[i] = "";
}
-
delete throwfile;
}
// Make sure we have a nominal file
if (!nominalfound or nominalfile.empty()){
ERR(FTL) << "No nominal found when mergining! Exiting!" << std::endl;
throw;
}
-
-
- // Get the nominal throws file
+ // Get the nominal throws file
TFile* tempfile = new TFile((nominalfile).c_str(),"READ");
tempfile->cd();
TDirectory* nominal = (TDirectory*)tempfile->Get("nominal");
- // int nthrows = FitPar::Config().GetParI("error_throws");
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Check percentage of bad files is okay.
int badfilecount = 0;
for (uint i = 0; i < fThrowList.size(); i++){
if (!fThrowList[i].empty()){
LOG(FIT) << "Loading Throws From File " << i << " : "
- << fThrowList[i] << std::endl;
+ << fThrowList[i] << std::endl;
} else {
badfilecount++;
}
}
// Check we have at least one good file
if ((uint)badfilecount == fThrowList.size()){
ERR(FTL) << "Found no good throw files for MergeThrows" << std::endl;
throw;
} else if (badfilecount > fThrowList.size()*0.25){
ERR(WRN) << "Over 25% of your throw files are dodgy. Please check this is okay!" << std::endl;
ERR(WRN) << "Will continue for the time being..." << std::endl;
sleep(5);
}
// Now go through the keys in the temporary file and look for TH1D, and TH2D plots
TIter next(nominal->GetListOfKeys());
TKey *key;
while ((key = (TKey*)next())) {
TClass *cl = gROOT->GetClass(key->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
TH1* baseplot = (TH1D*)key->ReadObj();
std::string plotname = std::string(baseplot->GetName());
LOG(FIT) << "Creating error bands for " << plotname;
if (LOG_LEVEL(FIT)){
if (!uniformly) std::cout << " : Using COVARIANCE Throws! " << std::endl;
else std::cout << " : Using UNIFORM THROWS!!! " << std::endl;
}
int nbins = 0;
if (cl->InheritsFrom("TH1D")) nbins = ((TH1D*)baseplot)->GetNbinsX();
else nbins = ((TH1D*)baseplot)->GetNbinsX()* ((TH1D*)baseplot)->GetNbinsY();
// Setup TProfile with RMS option
TProfile* tprof = new TProfile((plotname + "_prof").c_str(),(plotname + "_prof").c_str(),nbins, 0, nbins, "S");
// Setup The TTREE
double* bincontents;
bincontents = new double[nbins];
double* binlowest;
binlowest = new double[nbins];
double* binhighest;
binhighest = new double[nbins];
errorDIR->cd();
TTree* bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str());
for (Int_t i = 0; i < nbins; i++){
bincontents[i] = 0.0;
binhighest[i] = 0.0;
binlowest[i] = 0.0;
bintree->Branch(Form("content_%i",i),&bincontents[i],Form("content_%i/D",i));
}
// Make new throw plot
TH1* newplot;
// Run Throw Merging.
for (UInt_t i = 0; i < fThrowList.size(); i++){
TFile* throwfile = new TFile(fThrowList[i].c_str(), "READ");
// Loop over all throws in a folder
TIter nextthrow(throwfile->GetListOfKeys());
TKey *throwkey;
while ((throwkey = (TKey*)nextthrow())) {
-
+
// Skip non throw folders
if (std::string(throwkey->GetName()).find("throw_") == std::string::npos) continue;
// Get Throw DIR
TDirectory* throwdir = (TDirectory*)throwkey->ReadObj();
// Get Plot From Throw
newplot = (TH1*)throwdir->Get(plotname.c_str());
if (!newplot) continue;
// Loop Over Plot
for (Int_t j = 0; j < nbins; j++){
tprof->Fill(j+0.5, newplot->GetBinContent(j+1));
bincontents[j] = newplot->GetBinContent(j+1);
if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j];
if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j];
}
errorDIR->cd();
bintree->Fill();
}
throwfile->Close();
delete throwfile;
}
errorDIR->cd();
if (uniformly){
LOG(FIT) << "Uniformly Calculating Plot Errors!" << std::endl;
}
TH1* statplot = (TH1*) baseplot->Clone();
for (Int_t j = 0; j < nbins; j++){
if (!uniformly){
-// if ((baseplot->GetBinError(j+1)/baseplot->GetBinContent(j+1)) < 1.0) {
- // baseplot->SetBinError(j+1,sqrt(pow(tprof->GetBinError(j+1),2) + pow(baseplot->GetBinError(j+1),2)));
- // } else {
- //baseplot->SetBinContent(j+1,tprof->GetBinContent(j+1));
- baseplot->SetBinError(j+1,tprof->GetBinError(j+1));
- // }
+ // if ((baseplot->GetBinError(j+1)/baseplot->GetBinContent(j+1)) < 1.0) {
+ // baseplot->SetBinError(j+1,sqrt(pow(tprof->GetBinError(j+1),2) + pow(baseplot->GetBinError(j+1),2)));
+ // } else {
+ //baseplot->SetBinContent(j+1,tprof->GetBinContent(j+1));
+ baseplot->SetBinError(j+1,tprof->GetBinError(j+1));
+ // }
} else {
- baseplot->SetBinContent(j+1, 0.0);//(binlowest[j] + binhighest[j]) / 2.0);
+ baseplot->SetBinContent(j+1, 0.0);//(binlowest[j] + binhighest[j]) / 2.0);
baseplot->SetBinError(j+1, 0.0); //(binhighest[j] - binlowest[j])/2.0);
}
}
errorDIR->cd();
baseplot->Write();
tprof->Write();
bintree->Write();
outnominal->cd();
for (int i = 0; i < nbins; i++){
baseplot->SetBinError(i+1, sqrt(pow(statplot->GetBinError(i+1),2) + pow(baseplot->GetBinError(i+1),2)));
}
baseplot->Write();
-
+
delete statplot;
delete baseplot;
delete tprof;
delete bintree;
delete [] bincontents;
}
return;
};
-void SystematicRoutines::EigenErrors(){
+void SystematicRoutines::EigenErrors() {
- fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
- fOutputRootFile->cd();
+ fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
+ fOutputRootFile->cd();
// Make Covariance
TMatrixDSym* fullcovar = new TMatrixDSym( fParams.size() );
// Extract covariance from all loaded ParamPulls
for (PullListConstIter iter = fInputThrows.begin();
- iter != fInputThrows.end(); iter++){
+ iter != fInputThrows.end(); iter++){
ParamPull* pull = *iter;
// Check pull is actualyl Gaussian
std::string pulltype = pull->GetType();
if (pulltype.find("GAUSTHROW") == std::string::npos){
THROW("Can only calculate EigenErrors for Gaussian pulls!");
}
// Get data and covariances
TH1D dialhist = pull->GetDataHist();
TH2D covhist = pull->GetFullCovar();
// Loop over all dials and compare names
for (size_t pari = 0; pari < fParams.size(); pari++){
for (size_t parj = 0; parj < fParams.size(); parj++){
-
- std::string name_pari = fParams[pari];
- std::string name_parj = fParams[parj];
-
- // Compare names to those in the pull
- for (int pulli = 0; pulli < dialhist.GetNbinsX(); pulli++){
- for (int pullj = 0; pullj < dialhist.GetNbinsX(); pullj++){
-
- std::string name_pulli = dialhist.GetXaxis()->GetBinLabel(pulli+1);
- std::string name_pullj = dialhist.GetXaxis()->GetBinLabel(pullj+1);
-
- if (name_pulli == name_pari && name_pullj == name_parj){
- (*fullcovar)[pari][parj] = covhist.GetBinContent(pulli+1, pullj+1);
- fCurVals[name_pari] = dialhist.GetBinContent(pulli+1);
- fCurVals[name_parj] = dialhist.GetBinContent(pullj+1);
- }
-
- }
- }
-
+
+ std::string name_pari = fParams[pari];
+ std::string name_parj = fParams[parj];
+
+ // Compare names to those in the pull
+ for (int pulli = 0; pulli < dialhist.GetNbinsX(); pulli++){
+ for (int pullj = 0; pullj < dialhist.GetNbinsX(); pullj++){
+
+ std::string name_pulli = dialhist.GetXaxis()->GetBinLabel(pulli+1);
+ std::string name_pullj = dialhist.GetXaxis()->GetBinLabel(pullj+1);
+
+ if (name_pulli == name_pari && name_pullj == name_parj){
+ (*fullcovar)[pari][parj] = covhist.GetBinContent(pulli+1, pullj+1);
+ fCurVals[name_pari] = dialhist.GetBinContent(pulli+1);
+ fCurVals[name_parj] = dialhist.GetBinContent(pullj+1);
+ }
+
+ }
+ }
+
}
}
}
/*
- TFile* test = new TFile("testingcovar.root","RECREATE");
- test->cd();
- TH2D* joinedcov = new TH2D("COVAR","COVAR",
- fullcovar->GetNrows(), 0.0, float(fullcovar->GetNrows()),
- fullcovar->GetNrows(), 0.0, float(fullcovar->GetNrows()));
- for (int i = 0; i < fullcovar->GetNrows(); i++){
- for (int j = 0; j < fullcovar->GetNcols(); j++){
- joinedcov->SetBinContent(i+1, j+1, (*fullcovar)[i][j]);
- }
- }
- joinedcov->Write("COVAR");
- test->Close();
- */
+ TFile* test = new TFile("testingcovar.root","RECREATE");
+ test->cd();
+ TH2D* joinedcov = new TH2D("COVAR","COVAR",
+ fullcovar->GetNrows(), 0.0, float(fullcovar->GetNrows()),
+ fullcovar->GetNrows(), 0.0, float(fullcovar->GetNrows()));
+ for (int i = 0; i < fullcovar->GetNrows(); i++){
+ for (int j = 0; j < fullcovar->GetNcols(); j++){
+ joinedcov->SetBinContent(i+1, j+1, (*fullcovar)[i][j]);
+ }
+ }
+ joinedcov->Write("COVAR");
+ test->Close();
+ */
// Calculator all EigenVectors and EigenValues
TMatrixDSymEigen* eigen = new TMatrixDSymEigen(*fullcovar);
const TVectorD eigenVals = eigen->GetEigenValues();
const TMatrixD eigenVect = eigen->GetEigenVectors();
eigenVals.Print();
eigenVect.Print();
TDirectory* outnominal = (TDirectory*) fOutputRootFile->mkdir("nominal");
outnominal->cd();
double *valst = FitUtils::GetArrayFromMap( fParams, fCurVals );
//double chi2 = fSampleFCN->DoEval( valst );
delete valst;
fSampleFCN->Write();
// Loop over all throws
TDirectory* throwsdir = (TDirectory*) fOutputRootFile->mkdir("throws");
throwsdir->cd();
int count = 0;
// Produce all error throws.
for (int i = 0; i < eigenVect.GetNrows(); i++){
TDirectory* throwfolder = (TDirectory*)throwsdir->mkdir(Form("throw_%i",count));
throwfolder->cd();
// Get New Parameter Vector
LOG(FIT) << "Parameter Set " << count << std::endl;
for (int j = 0; j < eigenVect.GetNrows(); j++){
std::string param = fParams[j];
LOG(FIT) << " " << j << ". " << param << " : " << fCurVals[param] + sqrt(eigenVals[i]) * eigenVect[j][i] << std::endl;
fThrownVals[param] = fCurVals[param] + sqrt(eigenVals[i]) * eigenVect[j][i];
}
// Run Eval
double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals );
double chi2 = fSampleFCN->DoEval( vals );
delete vals;
count++;
fSampleFCN->Write();
-
+
throwfolder = (TDirectory*)throwsdir->mkdir(Form("throw_%i",count));
throwfolder->cd();
// Get New Parameter Vector
LOG(FIT) << "Parameter Set " << count << std::endl;
for (int j = 0; j < eigenVect.GetNrows(); j++){
std::string param = fParams[j];
LOG(FIT) << " " << j << ". " << param << " : " <DoEval( vals2 );
delete vals2;
count++;
-
+
// Save the FCN
fSampleFCN->Write();
-
+
}
fOutputRootFile->Close();
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "UPDATE");
fOutputRootFile->cd();
throwsdir = (TDirectory*) fOutputRootFile->Get("throws");
outnominal = (TDirectory*) fOutputRootFile->Get("nominal");
// Loop through Error DIR
TDirectory* outerr = (TDirectory*) fOutputRootFile->mkdir("errors");
outerr->cd();
TIter next(outnominal->GetListOfKeys());
TKey *key;
while ((key = (TKey*)next())) {
-
+
TClass *cl = gROOT->GetClass(key->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
-
+
LOG(FIT) << "Creating error bands for " << key->GetName() << std::endl;
std::string plotname = std::string(key->GetName());
if (plotname.find("_EVT") != std::string::npos) continue;
if (plotname.find("_FLUX") != std::string::npos) continue;
if (plotname.find("_FLX") != std::string::npos) continue;
TH1* baseplot = (TH1D*)key->ReadObj()->Clone(Form("%s_ORIGINAL",key->GetName()));
TH1* errorplot_upper = (TH1D*)baseplot->Clone(Form("%s_ERROR_UPPER",key->GetName()));
TH1* errorplot_lower = (TH1D*)baseplot->Clone(Form("%s_ERROR_LOWER", key->GetName()));
TH1* meanplot = (TH1D*)baseplot->Clone(Form("%s_SET_MEAN", key->GetName()));
TH1* systplot = (TH1D*)baseplot->Clone(Form("%s_SYST", key->GetName()));
TH1* statplot = (TH1D*)baseplot->Clone(Form("%s_STAT", key->GetName()));
TH1* totlplot = (TH1D*)baseplot->Clone(Form("%s_TOTAL", key->GetName()));
int nbins = 0;
if (cl->InheritsFrom("TH1D")) nbins = ((TH1D*)baseplot)->GetNbinsX();
else nbins = ((TH1D*)baseplot)->GetNbinsX()* ((TH1D*)baseplot)->GetNbinsY();
meanplot->Reset();
errorplot_upper->Reset();
errorplot_lower->Reset();
for (int j = 0; j < nbins; j++){
errorplot_upper->SetBinError(j+1, 0.0);
errorplot_lower->SetBinError(j+1, 0.0);
}
// Loop over throws and calculate mean and error for +- throws
int addcount = 0;
// Add baseplot first to slightly bias to central value
meanplot->Add(baseplot);
addcount++;
for (int i = 0; i < count; i++){
TH1* newplot = (TH1D*) throwsdir->Get(Form("throw_%i/%s",i,plotname.c_str()));
if (!newplot){
- ERR(WRN) << "Cannot find new plot : " << Form("throw_%i/%s",i,plotname.c_str()) << std::endl;
- ERR(WRN) << "This plot will not have the correct errors!" << std::endl;
- continue;
+ ERR(WRN) << "Cannot find new plot : " << Form("throw_%i/%s",i,plotname.c_str()) << std::endl;
+ ERR(WRN) << "This plot will not have the correct errors!" << std::endl;
+ continue;
}
newplot->SetDirectory(0);
nbins = newplot->GetNbinsX();
-
+
for (int j = 0; j < nbins; j++){
- if (i % 2 == 0){
- // std::cout << plotname<< " : upper " << errorplot_upper->GetBinContent(j+1) << " adding " << pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2) << std::endl;
- // std::cout << " -> " << baseplot->GetBinContent(j+1) << " " <GetBinContent(j+1) << std::endl;
- errorplot_upper->SetBinContent(j+1, errorplot_upper->GetBinContent(j+1) +
- pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2));
- // newplot->Print();
- } else {
- // std::cout << plotname << " : lower " << errorplot_lower->GetBinContent(j+1) << " adding " << pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2) << std::endl;
- // std::cout << " -> " << baseplot->GetBinContent(j+1) << " " << newplot->GetBinContent(j+1) << std::endl;
- errorplot_lower->SetBinContent(j+1, errorplot_lower->GetBinContent(j+1) +
- pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2));
- // newplot->Print();
- }
- meanplot->SetBinContent(j+1, meanplot->GetBinContent(j+1) + baseplot->GetBinContent(j+1));
+ if (i % 2 == 0){
+ // std::cout << plotname<< " : upper " << errorplot_upper->GetBinContent(j+1) << " adding " << pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2) << std::endl;
+ // std::cout << " -> " << baseplot->GetBinContent(j+1) << " " <GetBinContent(j+1) << std::endl;
+ errorplot_upper->SetBinContent(j+1, errorplot_upper->GetBinContent(j+1) +
+ pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2));
+ // newplot->Print();
+ } else {
+ // std::cout << plotname << " : lower " << errorplot_lower->GetBinContent(j+1) << " adding " << pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2) << std::endl;
+ // std::cout << " -> " << baseplot->GetBinContent(j+1) << " " << newplot->GetBinContent(j+1) << std::endl;
+ errorplot_lower->SetBinContent(j+1, errorplot_lower->GetBinContent(j+1) +
+ pow(baseplot->GetBinContent(j+1) - newplot->GetBinContent(j+1),2));
+ // newplot->Print();
+ }
+ meanplot->SetBinContent(j+1, meanplot->GetBinContent(j+1) + baseplot->GetBinContent(j+1));
}
delete newplot;
addcount++;
}
-
+
// Get mean Average
for (int j = 0; j < nbins; j++){
meanplot->SetBinContent(j+1, meanplot->GetBinContent(j+1)/double(addcount));
}
for (int j = 0; j < nbins; j++){
errorplot_upper->SetBinContent(j+1, sqrt(errorplot_upper->GetBinContent(j+1)));
errorplot_lower->SetBinContent(j+1, sqrt(errorplot_lower->GetBinContent(j+1)));
-
+
statplot->SetBinError(j+1, baseplot->GetBinError(j+1) );
systplot->SetBinError(j+1, (errorplot_upper->GetBinContent(j+1) + errorplot_lower->GetBinContent(j+1))/2.0);
totlplot->SetBinError(j+1, sqrt( pow(statplot->GetBinError(j+1),2) + pow(systplot->GetBinError(j+1),2) ) );
meanplot->SetBinError(j+1, sqrt( pow(statplot->GetBinError(j+1),2) + pow(systplot->GetBinError(j+1),2) ) );
-
+
}
outerr->cd();
errorplot_upper->Write();
errorplot_lower->Write();
baseplot->Write();
meanplot->Write();
statplot->Write();
systplot->Write();
totlplot->Write();
delete errorplot_upper;
delete errorplot_lower;
delete baseplot;
delete meanplot;
delete statplot;
delete systplot;
delete totlplot;
}
}