diff --git a/src/Config/NuisConfig.cxx b/src/Config/NuisConfig.cxx
index 7e78d66..38ca696 100644
--- a/src/Config/NuisConfig.cxx
+++ b/src/Config/NuisConfig.cxx
@@ -1,817 +1,817 @@
// 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 "NuisConfig.h"
#include "FitLogger.h"
#include "GeneralUtils.h"
#include "TXMLEngine.h"
namespace Config {
nuisconfig &Get() { return nuisconfig::GetConfig(); };
std::string GetPar(std::string const &name) { return Get().GetConfig(name); }
bool HasPar(std::string const &name) { return Get().HasConfig(name); }
std::string GetParS(std::string const &name) { return Get().GetConfigS(name); }
int GetParI(std::string const &name) { return Get().GetConfigI(name); }
bool GetParB(std::string const &name) { return Get().GetConfigB(name); }
float GetParF(std::string const &name) { return Get().GetConfigF(name); }
double GetParD(std::string const &name) { return Get().GetConfigD(name); }
void SetPar(std::string const &name, std::string const &val) {
Get().SetConfig(name, val);
}
void SetPar(std::string const &name, char const *val) {
Get().SetConfig(name, val);
}
void SetPar(std::string const &name, bool val) { Get().SetConfig(name, val); }
void SetPar(std::string const &name, int val) { Get().SetConfig(name, val); }
void SetPar(std::string const &name, float val) { Get().SetConfig(name, val); }
void SetPar(std::string const &name, double val) { Get().SetConfig(name, val); }
}
namespace FitPar {
std::string GetDataBase() { return GeneralUtils::GetTopLevelDir() + "/data/"; };
nuisconfig &Config() { return Config::Get(); };
}
nuisconfig *nuisconfig::m_nuisconfigInstance = NULL;
nuisconfig &nuisconfig::GetConfig(void) {
if (!m_nuisconfigInstance) m_nuisconfigInstance = new nuisconfig;
return *m_nuisconfigInstance;
};
// Main Class Definition
nuisconfig::nuisconfig() {
// Load default Parameters
std::string filename =
(GeneralUtils::GetTopLevelDir() + "/parameters/config.xml");
std::cout << "[ NUISANCE ]: Loading DEFAULT settings from : " << filename;
// Create XML Engine
fXML = new TXMLEngine;
fXML->SetSkipComments(true);
// Load in documents
fXMLDocs.push_back(fXML->ParseFile(filename.c_str(), 1000000));
if (!fXMLDocs[0]) {
THROW("Cannot Read Parameters File!");
}
// Setup Main XML Node to be the first file read
fMainNode = fXML->DocGetRootElement(fXMLDocs[0]);
// Print result
std::cout << " -> DONE." << std::endl;
}
nuisconfig::~nuisconfig() {
// Should really delete XML objects here but we don't
}
void nuisconfig::LoadSettings(std::string const &filename,
std::string const &state) {
// Open file and see if its XML
std::cout << "[ NUISANCE ]: Trying to parse file : " << filename;
StopTalking();
XMLDocPointer_t readdoc = fXML->ParseFile(filename.c_str(), 1000000);
StartTalking();
// If it is parse it as a nice XML config file
if (readdoc) {
std::cout << " -> Found XML file." << std::endl;
LoadXMLSettings(filename, state);
// Otherwise its an old simple card file
} else {
std::cout << " -> Assuming its a simple card file." << std::endl;
LoadCardSettings(filename, state);
}
}
void nuisconfig::LoadXMLSettings(std::string const &filename,
std::string const &state) {
std::cout << "[ NUISANCE ]: Loading XML settings from : " << filename
<< std::endl;
// Add new file to xml docs list
fXMLDocs.push_back(fXML->ParseFile(filename.c_str(), 1000000));
if (!fXMLDocs.back()) {
THROW("Failed to read: " << filename);
}
// Loop over children and add
XMLNodePointer_t child =
fXML->GetChild(fXML->DocGetRootElement(fXMLDocs.back()));
// // Here we manually load all the children from the card file into our root
// node
if (!child) {
THROW("CANNOT Find child inside settings file!");
}
while (child) {
// SPECIAL CONFIG CASE
// If its a config node, then remove previous attributes, overriding old
// value
if (!std::string(fXML->GetNodeName(child)).compare("config")) {
// Loop over attribues
XMLAttrPointer_t attr1 = fXML->GetFirstAttr(child);
while (attr1) {
// If a valid attribute name is given then compare
if (!GetConfigS(fXML->GetAttrName(attr1)).empty()) {
// Get full list of present configs
std::vector<XMLNodePointer_t> confignodes = GetNodes("config");
// Loop over present configs and compare
for (size_t i = 0; i < confignodes.size(); i++) {
// If we already have this config, free the old attribute
if (fXML->HasAttr(confignodes[i], fXML->GetAttrName(attr1))) {
fXML->FreeAttr(confignodes[i], fXML->GetAttrName(attr1));
break;
}
}
}
// Move onto next config attribute
attr1 = fXML->GetNextAttr(attr1);
}
}
TString nodeStr;
fXML->SaveSingleNode(child, &nodeStr);
XMLNodePointer_t copyNode = fXML->ReadSingleNode(nodeStr.Data());
// std::cout << "copying node..." << std::endl;
// PrintXML(copyNode);
// Add this child to the main config list
fXML->AddChild(fMainNode, copyNode);
// std::cout << "Done, was it added?" << std::endl;
// PrintXML(fMainNode);
// Get Next Child
child = fXML->GetNext(child);
}
std::cout << " -> DONE." << std::endl;
}
void nuisconfig::LoadCardSettings(std::string const &filename,
std::string const &state) {
std::cout << "[ NUISANCE ]: Loading simple config from : " << filename;
// Build XML Config from the card file by parsing each line
std::vector<std::string> cardlines =
GeneralUtils::ParseFileToStr(filename, "\n");
int linecount = 0;
// Loop over all input lines
for (std::vector<std::string>::iterator iter = cardlines.begin();
iter != cardlines.end(); iter++) {
std::string line = (*iter);
linecount++;
// Skip Comments
if (line.empty()) continue;
if (line.c_str()[0] == '#') continue;
// Parse whitespace
std::vector<std::string> strvct = GeneralUtils::ParseToStr(line, " ");
if (strvct.empty()) continue;
// Get Identifier
std::string id = strvct[0];
// Build backwards compatible xml configs
// Sample structure
if (!id.compare("sample")) {
CreateSampleNodeFromLine(line);
}
// Any parameter structure
if (id.find("_parameter") != std::string::npos) {
CreateParameterNodeFromLine(line);
}
// Any covar structure
if (!id.compare("covar") || !id.compare("pull") || !id.compare("throw")) {
CreatePullNodeFromLine(line);
}
// Any config structure
if (!id.compare("config")) {
CreateOldConfigNodeFromLine(line);
}
}
std::cout << " -> DONE." << std::endl;
}
XMLNodePointer_t nuisconfig::CreateSampleNodeFromLine(std::string const &line) {
// Create new node entry
XMLNodePointer_t samplenode = CreateNode("sample");
// Parse line
std::vector<std::string> strvct = GeneralUtils::ParseToStr(line, " ");
// Add line elements to the node
// name input type norm
if (strvct.size() > 1) Set(samplenode, "name", strvct[1]);
if (strvct.size() > 2) Set(samplenode, "input", strvct[2]);
if (strvct.size() > 3) Set(samplenode, "type", strvct[3]);
if (strvct.size() > 4) Set(samplenode, "norm", strvct[4]);
return samplenode;
}
XMLNodePointer_t nuisconfig::CreateParameterNodeFromLine(
std::string const &line) {
// Create new node entry
XMLNodePointer_t parnode = CreateNode("parameter");
// Parse line
std::vector<std::string> strvct = GeneralUtils::ParseToStr(line, " ");
// Add line elements to the node
// type name nominal [low] [high] [step] state
if (strvct.size() > 0) Set(parnode, "type", strvct[0]);
if (strvct.size() > 1) Set(parnode, "name", strvct[1]);
if (strvct.size() > 2) Set(parnode, "nominal", strvct[2]);
// If free structure
if (strvct.size() == 7) {
Set(parnode, "low", strvct[3]);
Set(parnode, "high", strvct[4]);
Set(parnode, "step", strvct[5]);
Set(parnode, "state", strvct[6]);
// Fixed param structure
} else if (strvct.size() == 3) {
Set(parnode, "state", "FIX");
} else if (strvct.size() == 4) {
Set(parnode, "state", strvct[3]);
}
return parnode;
}
XMLNodePointer_t nuisconfig::CreatePullNodeFromLine(std::string const &line) {
// Create new node entry
XMLNodePointer_t parnode = CreateNode("covar");
// Parse line
std::vector<std::string> strvct = GeneralUtils::ParseToStr(line, " ");
// Add line elements to the node
// name input type
if (strvct.size() > 1) Set(parnode, "name", strvct[1]);
if (strvct.size() > 2) Set(parnode, "input", strvct[2]);
if (strvct.size() > 3) Set(parnode, "type", strvct[3]);
return parnode;
}
XMLNodePointer_t nuisconfig::CreateOldConfigNodeFromLine(
std::string const &line) {
// Create new node entry
XMLNodePointer_t confignode = CreateNode("config");
// Parse line
std::vector<std::string> strvct = GeneralUtils::ParseToStr(line, " ");
// Add line elements to the node
// name value
if (strvct.size() > 2) Set(confignode, strvct[1], strvct[2]);
return confignode;
}
void nuisconfig::FinaliseSettings(std::string const &name) {
std::cout << "[ NUISANCE ]: Finalising run settings";
WriteSettings(name);
// Save full config to file
RemoveEmptyNodes();
RemoveIdenticalNodes();
std::cout << " -> DONE." << std::endl;
}
void nuisconfig::WriteSettings(std::string const &outputname) {
// Create a New XML Doc
XMLDocPointer_t newxmldoc = fXML->NewDoc();
fXML->DocSetRootElement(newxmldoc, fMainNode);
// Save document to file
if (GetConfigB("SaveParsedXMLFile")) {
fXML->SaveDoc(newxmldoc, outputname.c_str());
}
}
void nuisconfig::PrintXML(XMLNodePointer_t node, int indent) {
if (!node) {
node = fMainNode;
}
std::stringstream ss("");
for (int i = 0; i < indent; ++i) {
ss << " ";
}
std::cout << ss.str() << "<" << fXML->GetNodeName(node) << std::flush;
XMLAttrPointer_t attr = fXML->GetFirstAttr(node);
while (attr) {
std::cout << " " << fXML->GetAttrName(attr) << "=\""
<< fXML->GetAttrValue(attr) << "\"" << std::flush;
attr = fXML->GetNextAttr(attr);
}
if (!fXML->GetChild(node)) {
std::cout << " />" << std::endl;
return;
}
std::cout << " >" << std::endl;
XMLNodePointer_t child = fXML->GetChild(node);
while (child) {
PrintXML(child, indent + 1);
child = fXML->GetNext(child);
}
std::cout << ss.str() << "</" << fXML->GetNodeName(node) << ">" << std::endl;
}
XMLNodePointer_t nuisconfig::CreateNode(std::string const &name) {
return fXML->NewChild(fMainNode, 0, name.c_str());
}
XMLNodePointer_t nuisconfig::CreateNode(XMLNodePointer_t node,
std::string const &name) {
return fXML->NewChild(node, 0, name.c_str());
}
XMLNodePointer_t nuisconfig::GetNode(XMLNodePointer_t node,
std::string const &type) {
/// Loop over all children
XMLNodePointer_t child = fXML->GetChild(node);
while (child != 0) {
/// Get nodes for given type (if type empty return all)
if (std::string(fXML->GetNodeName(child)) == type.c_str() or type.empty()) {
return child;
}
// Next child
child = fXML->GetNext(child);
}
// Child not found
return 0;
}
void nuisconfig::RemoveEmptyNodes() {
std::vector<XMLNodePointer_t> nodelist = Config::Get().GetNodes();
for (size_t i = 0; i < nodelist.size(); i++) {
if (fXML->IsEmptyNode(nodelist[i])) {
std::cout << "Removing empty node: " << fXML->GetNodeName(nodelist[i])
<< ", with child ?" << bool(fXML->GetChild(nodelist[i]))
<< std::endl;
RemoveNode(nodelist[i]);
}
}
}
void nuisconfig::RemoveIdenticalNodes() {
std::vector<XMLNodePointer_t> removed;
// Loop over all nodes and check for identical nodes
std::vector<XMLNodePointer_t> nodelist = Config::Get().GetNodes();
for (size_t i = 0; i < nodelist.size(); i++) {
for (size_t j = 0; j < nodelist.size(); j++) {
if (i == j) continue;
XMLNodePointer_t node1 = nodelist[i];
XMLNodePointer_t node2 = nodelist[j];
// Check node already removed.
if (std::find(removed.begin(), removed.end(), node1) != removed.end()) {
continue;
}
if (std::find(removed.begin(), removed.end(), node2) != removed.end()) {
continue;
}
// Check matching
if (!MatchingNodes(node1, node2)) {
continue;
}
if (std::string(fXML->GetNodeName(node1)).compare("config") and
fXML->IsEmptyNode(node1)) {
// Matching so print out warning
std::cout << "Matching nodes given! Removing node1!" << std::endl
<< "Node 1" << std::endl;
PrintNode(node1);
std::cout << "Node 2" << std::endl;
PrintNode(node2);
}
// Remove node
removed.push_back(node1);
}
}
// Now go through and remove this node.
for (size_t i = 0; i < removed.size(); i++) {
RemoveNode(removed.at(i));
}
return;
}
void nuisconfig::RemoveNode(XMLNodePointer_t node) {
std::cout << "[INFO]: Removing node: " << fXML->GetNodeName(node)
<< std::endl;
fXML->FreeAllAttr(node);
fXML->CleanNode(node);
fXML->FreeNode(node);
fXML->UnlinkNode(node);
}
void nuisconfig::PrintNode(XMLNodePointer_t node) {
// Print Node Name
std::cout << fXML->GetNodeName(node) << std::endl;
// Loop and print all attributes
XMLAttrPointer_t attr = fXML->GetFirstAttr(node);
while (attr != 0) {
std::cout << " -> " << fXML->GetAttrName(attr) << " : "
<< fXML->GetAttrValue(attr) << std::endl;
attr = fXML->GetNextAttr(attr);
}
}
bool nuisconfig::MatchingNodes(XMLNodePointer_t node1, XMLNodePointer_t node2) {
bool matching = true;
XMLAttrPointer_t attr = fXML->GetFirstAttr(node1);
while (attr != 0) {
if (GetS(node2, fXML->GetAttrName(attr)) != fXML->GetAttrValue(attr))
matching = false;
attr = fXML->GetNextAttr(attr);
}
return matching;
}
XMLNodePointer_t nuisconfig::GetNode(std::string const &type) {
return GetNode(fMainNode, type);
}
std::vector<XMLNodePointer_t> nuisconfig::GetNodes(XMLNodePointer_t node,
std::string const &type) {
// Create new vector for nodes
std::vector<XMLNodePointer_t> nodelist;
/// Loop over all children
XMLNodePointer_t child = fXML->GetChild(node);
while (child != 0) {
/// Get nodes for given type (if type empty return all)
if (std::string(fXML->GetNodeName(child)) == type.c_str() or type.empty()) {
nodelist.push_back(child);
}
// Next child
child = fXML->GetNext(child);
}
// return list
return nodelist;
}
std::vector<XMLNodePointer_t> nuisconfig::GetNodes(std::string const &type) {
return GetNodes(fMainNode, type);
}
void nuisconfig::Set(XMLNodePointer_t node, std::string const &name,
std::string const &val) {
// Remove and re-add attribute
if (fXML->HasAttr(node, name.c_str())) {
fXML->FreeAttr(node, name.c_str());
}
fXML->NewAttr(node, 0, name.c_str(), val.c_str());
}
void nuisconfig::Set(XMLNodePointer_t node, std::string const &name,
char const *val) {
Set(node, name, std::string(val));
}
void nuisconfig::Set(XMLNodePointer_t node, std::string const &name, bool val) {
Set(node, name, GeneralUtils::BoolToStr(val));
}
void nuisconfig::Set(XMLNodePointer_t node, std::string const &name, int val) {
Set(node, name, GeneralUtils::IntToStr(val));
}
void nuisconfig::Set(XMLNodePointer_t node, std::string const &name,
float val) {
Set(node, name, GeneralUtils::DblToStr(val));
}
void nuisconfig::Set(XMLNodePointer_t node, std::string const &name,
double val) {
Set(node, name, GeneralUtils::DblToStr(val));
}
void nuisconfig::SetS(XMLNodePointer_t node, std::string const &name,
std::string const &val) {
Set(node, name, val);
}
void nuisconfig::SetB(XMLNodePointer_t node, std::string const &name,
bool val) {
Set(node, name, GeneralUtils::BoolToStr(val));
}
void nuisconfig::SetI(XMLNodePointer_t node, std::string const &name, int val) {
Set(node, name, GeneralUtils::IntToStr(val));
}
void nuisconfig::SetF(XMLNodePointer_t node, std::string const &name,
float val) {
Set(node, name, GeneralUtils::DblToStr(val));
}
void nuisconfig::SetD(XMLNodePointer_t node, std::string const &name,
double val) {
Set(node, name, GeneralUtils::DblToStr(val));
}
bool nuisconfig::Has(XMLNodePointer_t node, std::string const &name) {
// If node empty return empty
if (node == 0) return false;
// Search attributes
XMLAttrPointer_t attr = fXML->GetFirstAttr(node);
bool found = false;
// Loop over all attributes
while (attr != 0) {
// Find value of correct name
if (std::string(fXML->GetAttrName(attr)) == name.c_str()) {
found = true;
break;
}
// Next Attribute
attr = fXML->GetNextAttr(attr);
}
return found;
}
std::string nuisconfig::Get(XMLNodePointer_t node, std::string const &name) {
// If node empty return empty
if (node == 0) return "";
// Get Attribute from child with name
XMLAttrPointer_t attr = fXML->GetFirstAttr(node);
std::string temp = "";
// Loop over all attributes
while (attr != 0) {
// If valid match then save
if (std::string(fXML->GetAttrName(attr)) == name.c_str()) {
temp = fXML->GetAttrValue(attr);
}
// Next Attribute
attr = fXML->GetNextAttr(attr);
}
return temp;
}
std::string nuisconfig::GetS(XMLNodePointer_t node, std::string const &name) {
return Get(node, name);
}
bool nuisconfig::GetB(XMLNodePointer_t node, std::string const &name) {
std::string tempattr = Get(node, name);
return GeneralUtils::StrToBool(tempattr);
}
int nuisconfig::GetI(XMLNodePointer_t node, std::string const &name) {
std::string tempattr = Get(node, name);
return GeneralUtils::StrToInt(tempattr);
}
float nuisconfig::GetF(XMLNodePointer_t node, std::string const &name) {
std::string tempattr = Get(node, name);
return GeneralUtils::StrToDbl(tempattr);
}
double nuisconfig::GetD(XMLNodePointer_t node, std::string const &name) {
std::string tempattr = Get(node, name);
return GeneralUtils::StrToDbl(tempattr);
}
std::vector<std::string> nuisconfig::GetVS(XMLNodePointer_t node,
std::string const &name,
const char *del) {
std::string tempattr = Get(node, name);
return GeneralUtils::ParseToStr(tempattr, del);
}
// std::vector<int> nuisconfig::GetVB(XMLNodePointer_t node,
// std::string name,
// const char* del) {
// std::string tempattr = Get(node, name);
// return GeneralUtils::ParseToBool(tempattr, del);
// }
std::vector<int> nuisconfig::GetVI(XMLNodePointer_t node,
std::string const &name, const char *del) {
std::string tempattr = Get(node, name);
return GeneralUtils::ParseToInt(tempattr, del);
}
// std::vector<int> nuisconfig::GetVF(XMLNodePointer_t node,
// std::string name,
// const char* del) {
// std::string tempattr = Get(node, name);
// return GeneralUtils::ParseToDouble(tempattr, del);
// }
std::vector<double> nuisconfig::GetVD(XMLNodePointer_t node,
std::string const &name,
char const *del) {
std::string tempattr = Get(node, name);
return GeneralUtils::ParseToDbl(tempattr, del);
}
std::vector<std::string> nuisconfig::GetAllKeysForNode(XMLNodePointer_t node) {
- bool matching = true;
+ //bool matching = true;
XMLAttrPointer_t attr = fXML->GetFirstAttr(node);
std::vector<std::string> keys;
while (attr != 0) {
if (!std::string(fXML->GetAttrName(attr)).empty()) {
keys.push_back(std::string(fXML->GetAttrName(attr)));
}
attr = fXML->GetNextAttr(attr);
}
return keys;
}
XMLNodePointer_t nuisconfig::GetConfigNode(std::string const &name) {
// Loop over children and look for name
XMLNodePointer_t child = fXML->GetChild(fMainNode);
while (child != 0) {
// Select only config parameters
if (!std::string(fXML->GetNodeName(child)).compare("config")) {
// Loop over config attributes and search for name
XMLAttrPointer_t attr = fXML->GetFirstAttr(child);
while (attr != 0) {
// Save name value
if (std::string(fXML->GetAttrName(attr)) == name.c_str()) {
return child;
}
// Get Next Attribute
attr = fXML->GetNextAttr(attr);
}
}
// Next Child
child = fXML->GetNext(child);
}
return 0;
}
void nuisconfig::SetConfig(std::string const &name, std::string const &val) {
XMLNodePointer_t node = GetConfigNode(name);
if (!node) node = CreateNode("config");
Set(node, name, val);
}
void nuisconfig::SetConfig(std::string const &name, char const *val) {
SetConfig(name, std::string(val));
}
void nuisconfig::SetConfig(std::string const &name, bool val) {
XMLNodePointer_t node = GetConfigNode(name);
if (!node) node = CreateNode("config");
Set(node, name, val);
}
void nuisconfig::SetConfig(std::string const &name, int val) {
XMLNodePointer_t node = GetConfigNode(name);
if (!node) node = CreateNode("config");
Set(node, name, val);
}
void nuisconfig::SetConfig(std::string const &name, float val) {
XMLNodePointer_t node = GetConfigNode(name);
if (!node) node = CreateNode("config");
Set(node, name, val);
}
void nuisconfig::SetConfig(std::string const &name, double val) {
XMLNodePointer_t node = GetConfigNode(name);
if (!node) node = CreateNode("config");
Set(node, name, val);
}
void nuisconfig::OverrideConfig(std::string const &conf) {
std::vector<std::string> opts = GeneralUtils::ParseToStr(conf, "=");
SetConfig(opts[0], opts[1]);
}
std::string nuisconfig::GetConfig(std::string const &name) {
XMLNodePointer_t node = GetConfigNode(name);
if (!node) return "";
XMLAttrPointer_t attr = fXML->GetFirstAttr(node);
std::string temp = "";
// Loop config attributes
while (attr != 0) {
// Find match
if (std::string(fXML->GetAttrName(attr)) == name.c_str()) {
temp = fXML->GetAttrValue(attr);
}
// Get Next Attribute
attr = fXML->GetNextAttr(attr);
}
return temp;
}
bool nuisconfig::HasConfig(std::string const &name) {
return bool(GetConfigNode(name));
}
std::string nuisconfig::GetConfigS(std::string const &name) {
return GetConfig(name);
}
bool nuisconfig::GetConfigB(std::string const &name) {
std::string pars = GetConfig(name);
return GeneralUtils::StrToBool(pars);
}
int nuisconfig::GetConfigI(std::string const &name) {
std::string pars = GetConfig(name);
return GeneralUtils::StrToInt(pars);
}
float nuisconfig::GetConfigF(std::string const &name) {
std::string pars = GetConfig(name);
return GeneralUtils::StrToDbl(pars);
}
double nuisconfig::GetConfigD(std::string const &name) {
std::string pars = GetConfig(name);
return GeneralUtils::StrToDbl(pars);
}
std::string nuisconfig::GetParDIR(std::string const &parName) {
std::string outstr = this->GetParS(parName);
// Make replacements in the string
const int nfiletypes = 2;
const std::string filetypes[nfiletypes] = {"@data", "@nuisance"};
std::string filerepl[nfiletypes] = {FitPar::GetDataBase(),
FitPar::GetDataBase() + "/../"};
for (int i = 0; i < nfiletypes; i++) {
std::string findstring = filetypes[i];
std::string replstring = filerepl[i];
if (outstr.find(findstring) != std::string::npos) {
outstr.replace(outstr.find(findstring), findstring.size(), filerepl[i]);
break;
}
}
return outstr;
};
diff --git a/src/FitBase/SampleSettings.cxx b/src/FitBase/SampleSettings.cxx
index b94bb56..dc43b05 100644
--- a/src/FitBase/SampleSettings.cxx
+++ b/src/FitBase/SampleSettings.cxx
@@ -1,231 +1,231 @@
#include "SampleSettings.h"
SampleSettings::SampleSettings() {
};
SampleSettings::SampleSettings(nuiskey key) {
- fKeyValues = key;
- if (!key.Has("type")) key.Set("type", "DEFAULT");
+ fKeyValues = key;
+ if (!key.Has("type")) key.Set("type", "DEFAULT");
}
std::string SampleSettings::GetName() {
- return GetS("name");
+ return GetS("name");
}
void SampleSettings::SetS(std::string name, std::string val) {
- fKeyValues.SetS(name, val);
+ fKeyValues.SetS(name, val);
};
bool SampleSettings::Found(std::string name, std::string substr) {
- std::string compstring = fKeyValues.GetS(name);
- return compstring.find(substr) != std::string::npos;
+ std::string compstring = fKeyValues.GetS(name);
+ return compstring.find(substr) != std::string::npos;
}
void SampleSettings::SetXTitle(std::string name) {
- SetDefault("xtitle", name);
+ SetDefault("xtitle", name);
};
void SampleSettings::SetYTitle(std::string name) {
- SetDefault("ytitle", name);
+ SetDefault("ytitle", name);
};
void SampleSettings::SetZTitle(std::string name) {
- SetDefault("ztitle", name);
+ SetDefault("ztitle", name);
};
void SampleSettings::SetNormError(double norm) {
- SetDefault("norm_error", GeneralUtils::DblToStr(norm));
+ SetDefault("norm_error", GeneralUtils::DblToStr(norm));
};
double SampleSettings::GetNormError() {
return GetD("norm_error");
};
std::string SampleSettings::GetCovarInput() {
- return GetS("covar");
+ return GetS("covar");
}
void SampleSettings::SetAllowedTypes(std::string allowed, std::string defaulttype) {
- SetDefault("default_types", allowed);
- SetDefault("allowed_types", defaulttype);
+ SetDefault("default_types", allowed);
+ SetDefault("allowed_types", defaulttype);
};
void SampleSettings::SetEnuRangeFromFlux(TH1D* fluxhist) {
- double enu_min = fluxhist->GetXaxis()->GetXmin();
- double enu_max = fluxhist->GetXaxis()->GetXmax();
- SetEnuRange(enu_min, enu_max);
+ double enu_min = fluxhist->GetXaxis()->GetXmin();
+ double enu_max = fluxhist->GetXaxis()->GetXmax();
+ SetEnuRange(enu_min, enu_max);
};
void SampleSettings::SetEnuRange(double min, double max) {
- SetDefault("enu_min", min);
- SetDefault("enu_max", max);
+ SetDefault("enu_min", min);
+ SetDefault("enu_max", max);
};
void SampleSettings::Set(std::string name, double d) {
- SetDefault(name, d);
+ SetDefault(name, d);
}
void SampleSettings::Set(std::string name, int i) {
- SetDefault(name, i);
+ SetDefault(name, i);
}
void SampleSettings::Set(std::string name, std::string s) {
- SetDefault(name, s);
+ SetDefault(name, s);
}
void SampleSettings::DefineAllowedTargets(std::string targ) {
- // fAllowedTargets = TargetUtils::ParseTargetsToIntVect(targ);
+ // fAllowedTargets = TargetUtils::ParseTargetsToIntVect(targ);
};
void SampleSettings::FoundFill(std::string name, std::string substr, bool& cont, bool def) {
- std::string compstring = fKeyValues.GetS(name);
- if (compstring.find(substr) != std::string::npos) {
- cont = def;
- } else {
- cont = !def;
- }
+ std::string compstring = fKeyValues.GetS(name);
+ if (compstring.find(substr) != std::string::npos) {
+ cont = def;
+ } else {
+ cont = !def;
+ }
};
void SampleSettings::SetTitle(std::string val) {
- SetDefault("title", val);
+ SetDefault("title", val);
};
void SampleSettings::SetDataInput(std::string val) {
- SetDefault("data", val);
+ SetDefault("data", val);
};
std::string SampleSettings::GetMapInput() {
- return GetS("map");
+ return GetS("map");
}
void SampleSettings::SetMapInput(std::string val) {
- SetDefault("map", val);
+ SetDefault("map", val);
}
void SampleSettings::SetErrorInput(std::string val) {
- SetDefault("error", val);
+ SetDefault("error", val);
};
void SampleSettings::SetCovarInput(std::string val) {
- SetDefault("covar", val);
+ SetDefault("covar", val);
};
void SampleSettings::SetShapeCovarInput(std::string val) {
SetDefault("shapecovar", val);
};
void SampleSettings::SetDefault(std::string name, std::string val) {
- if (!fKeyValues.Has(name)) fKeyValues.Set(name, val);
+ if (!fKeyValues.Has(name)) fKeyValues.Set(name, val);
};
void SampleSettings::SetDefault(std::string name, double val ) {
- if (!fKeyValues.Has(name)) fKeyValues.Set(name, val);
+ if (!fKeyValues.Has(name)) fKeyValues.Set(name, val);
};
void SampleSettings::SetHasExtraHistograms(bool opt) {
- fHasExtraHistograms = opt;
+ fHasExtraHistograms = opt;
};
void SampleSettings::DefineAllowedSpecies(std::string species) {
- fAllowedTargets = BeamUtils::ParseSpeciesToIntVect(species);
+ fAllowedTargets = BeamUtils::ParseSpeciesToIntVect(species);
};
std::string SampleSettings::Title() {
- return GetS("title");
+ return GetS("title");
}
std::string SampleSettings::GetDataInput() {
- return GetS("data");
+ return GetS("data");
};
std::string SampleSettings::GetErrorInput() {
- return GetS("error");
+ return GetS("error");
};
std::string SampleSettings::PlotTitles() {
- return ";" + GetS("xtitle") + ";" + GetS("ytitle");
+ return ";" + GetS("xtitle") + ";" + GetS("ytitle");
};
std::string SampleSettings::GetFullTitles() {
- return Title() + PlotTitles();
+ return Title() + PlotTitles();
}
int SampleSettings::GetI(std::string name) {
- return fKeyValues.GetI(name);
+ return fKeyValues.GetI(name);
}
bool SampleSettings::GetB(std::string name){
- return fKeyValues.GetB(name);
+ return fKeyValues.GetB(name);
}
double SampleSettings::GetD(std::string name) {
- return fKeyValues.GetD(name);
+ return fKeyValues.GetD(name);
}
std::string SampleSettings::GetS(std::string name) {
- return fKeyValues.GetS(name);
+ return fKeyValues.GetS(name);
}
void SampleSettings::SetSuggestedFlux(std::string str) {
- SetS("suggested_flux", str);
+ SetS("suggested_flux", str);
}
void SampleSettings::SetDescription(std::string str) {
- SetS("description", str);
+ SetS("description", str);
}
void SampleSettings::Write(std::string name) {
- if (name.empty()) name = this->GetS("name") + "_settings";
+ if (name.empty()) name = this->GetS("name") + "_settings";
- // Make a new canvas
- TCanvas* c1 = new TCanvas( name.c_str(), name.c_str(), 800, 600 );
- c1->cd();
+ // Make a new canvas
+ TCanvas* c1 = new TCanvas( name.c_str(), name.c_str(), 800, 600 );
+ c1->cd();
- // Make new TPaveText
- TPaveText pttitle = TPaveText(0.05, 0.85, 0.95, 0.95);
- pttitle.AddText( name.c_str() );
- pttitle.SetTextAlign(11);
- pttitle.SetTextSize(15);
- pttitle.Draw();
+ // Make new TPaveText
+ TPaveText pttitle = TPaveText(0.05, 0.85, 0.95, 0.95);
+ pttitle.AddText( name.c_str() );
+ pttitle.SetTextAlign(11);
+ pttitle.SetTextSize(15);
+ pttitle.Draw();
- TPaveText pt = TPaveText(0.05, 0.05, 0.95, 0.80);
- std::vector<std::string> keys = fKeyValues.GetAllKeys();
- for (int i = 0; i < keys.size(); i++) {
+ TPaveText pt = TPaveText(0.05, 0.05, 0.95, 0.80);
+ std::vector<std::string> keys = fKeyValues.GetAllKeys();
+ for (size_t i = 0; i < keys.size(); i++) {
- std::string keyval = fKeyValues.GetS(keys[i]);
- std::vector<std::string> lines = GeneralUtils::ParseToStr(keyval, "\n");
+ std::string keyval = fKeyValues.GetS(keys[i]);
+ std::vector<std::string> lines = GeneralUtils::ParseToStr(keyval, "\n");
- if (lines.size() == 1) {
- pt.AddText( ("#bullet #bf{" + keys[i] + "} : " + fKeyValues.GetS(keys[i])).c_str() );
- } else {
- pt.AddText( ("#bullet #bf{" + keys[i] + "} : ").c_str() );
- for (int j = 0; j < lines.size(); j++) {
- pt.AddText((" |--> " + lines[j]).c_str() );
- }
- }
- }
+ if (lines.size() == 1) {
+ pt.AddText( ("#bullet #bf{" + keys[i] + "} : " + fKeyValues.GetS(keys[i])).c_str() );
+ } else {
+ pt.AddText( ("#bullet #bf{" + keys[i] + "} : ").c_str() );
+ for (size_t j = 0; j < lines.size(); j++) {
+ pt.AddText((" |--> " + lines[j]).c_str() );
+ }
+ }
+ }
- pt.SetTextAlign(11);
- pt.SetTextSize(14);
- pt.Draw("SAME");
+ pt.SetTextAlign(11);
+ pt.SetTextSize(14);
+ pt.Draw("SAME");
- c1->Write();
- delete c1;
+ c1->Write();
+ delete c1;
}
void SampleSettings::SetOnlyMC(bool state) {
- SetDefault("onlymc", state);
+ SetDefault("onlymc", state);
}
diff --git a/src/Logger/FitLogger.cxx b/src/Logger/FitLogger.cxx
index f675869..c94f8af 100644
--- a/src/Logger/FitLogger.cxx
+++ b/src/Logger/FitLogger.cxx
@@ -1,352 +1,352 @@
// 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 "FitLogger.h"
#include <fcntl.h>
#include <unistd.h>
namespace Logger {
// Logger Variables
int log_verb = 4;
bool use_colors = true;
bool showtrace = true;
std::ostream* __LOG_outstream(&std::cout);
std::ofstream __LOG_nullstream;
// Error Variables
int err_verb = 0;
std::ostream* __ERR_outstream(&std::cerr);
// Extra Variables
bool external_verb = false;
bool super_rainbow_mode = true; //!< For when fitting gets boring.
unsigned int super_rainbow_mode_colour = 0;
std::streambuf* default_cout = std::cout.rdbuf();
std::streambuf* default_cerr = std::cerr.rdbuf();
std::ofstream redirect_stream("/dev/null");
int silentfd = open("/dev/null", O_WRONLY);
int savedstdoutfd = dup(fileno(stdout));
int savedstderrfd = dup(fileno(stderr));
int nloggercalls = 0;
int timelastlog = 0;
}
// -------- Logging Functions --------- //
bool LOGGING(int level) {
// std::cout << "LOGGING : " << __FILENAME__ << " " << __FUNCTION__ <<
// std::endl;
return (Logger::log_verb >=
- (uint)__GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__));
+ (int)__GETLOG_LEVEL(level, __FILENAME__, __FUNCTION__));
};
int __GETLOG_LEVEL(int level, const char* filename, const char* funct) {
#ifdef __DEBUG__
int logfile = FitPar::Config().GetParI("logging." + std::string(filename));
if (logfile >= DEB and logfile <= EVT) {
level = logfile;
}
int logfunc = FitPar::Config().GetParI("logging." + std::string(funct));
if (logfunc >= DEB and logfunc <= EVT) {
level = logfunc;
}
#endif
return level;
};
std::ostream& __OUTLOG(int level, const char* filename, const char* funct,
int line) {
- if (Logger::log_verb < (unsigned int)level &&
- Logger::log_verb != (unsigned int)DEB) {
+ if (Logger::log_verb < (int)level &&
+ Logger::log_verb != (int)DEB) {
return (Logger::__LOG_nullstream);
} else {
if (Logger::use_colors) {
switch (level) {
case FIT:
std::cout << BOLDGREEN;
break;
case MIN:
std::cout << BOLDBLUE;
break;
case SAM:
std::cout << MAGENTA;
break;
case REC:
std::cout << BLUE;
break;
case SIG:
std::cout << GREEN;
break;
case DEB:
std::cout << CYAN;
break;
default:
break;
}
}
switch (level) {
case FIT:
std::cout << "[LOG Fitter]";
break;
case MIN:
std::cout << "[LOG Minmzr]";
break;
case SAM:
std::cout << "[LOG Sample]";
break;
case REC:
std::cout << "[LOG Reconf]";
break;
case SIG:
std::cout << "[LOG Signal]";
break;
case EVT:
std::cout << "[LOG Event ]";
break;
case DEB:
std::cout << "[LOG DEBUG ]";
break;
default:
std::cout << "[LOG INFO ]";
break;
}
// Apply indent
if (true) {
switch (level) {
case FIT:
std::cout << ": ";
break;
case MIN:
std::cout << ":- ";
break;
case SAM:
std::cout << ":-- ";
break;
case REC:
std::cout << ":--- ";
break;
case SIG:
std::cout << ":---- ";
break;
case EVT:
std::cout << ":----- ";
break;
case DEB:
std::cout << ":------ ";
break;
default:
std::cout << " ";
break;
}
}
if (Logger::use_colors) std::cout << RESET;
if (Logger::showtrace) {
std::cout << " : " << filename << "::" << funct << "[l. " << line
<< "] : ";
}
return *(Logger::__LOG_outstream);
}
}
void SETVERBOSITY(int level) { Logger::log_verb = level; }
void SETVERBOSITY(std::string verb) {
if (!verb.compare("DEB"))
Logger::log_verb = -1;
else if (!verb.compare("QUIET"))
Logger::log_verb = 0;
else if (!verb.compare("FIT"))
Logger::log_verb = 1;
else if (!verb.compare("MIN"))
Logger::log_verb = 2;
else if (!verb.compare("SAM"))
Logger::log_verb = 3;
else if (!verb.compare("REC"))
Logger::log_verb = 4;
else if (!verb.compare("SIG"))
Logger::log_verb = 5;
else if (!verb.compare("EVT"))
Logger::log_verb = 6;
else
Logger::log_verb = std::atoi(verb.c_str());
}
/// Set Trace Option
void SETTRACE(bool val) { Logger::showtrace = val; }
// ------ ERROR FUNCTIONS ---------- //
std::ostream& __OUTERR(int level, const char* filename, const char* funct,
int line) {
if (Logger::use_colors) std::cerr << RED;
switch (level) {
case FTL:
std::cerr << "[ERR FATAL ]: ";
break;
case WRN:
std::cerr << "[ERR WARN ]: ";
break;
}
if (Logger::use_colors) std::cerr << RESET;
// Allows enable error debugging trace
if (true or Logger::showtrace) {
std::cout << filename << "::" << funct << "[l. " << line << "] : ";
}
return *(Logger::__ERR_outstream);
}
// ----------- External Logging ----------- //
void SETEXTERNALVERBOSITY(int level) { Logger::external_verb = (level > 0); }
void StopTalking() {
// Check verbosity set correctly
if (!Logger::external_verb) return;
// Only redirect if we're not debugging
if (Logger::log_verb == (unsigned int)DEB) return;
std::cout.rdbuf(Logger::redirect_stream.rdbuf());
std::cerr.rdbuf(Logger::redirect_stream.rdbuf());
shhnuisancepythiaitokay_();
fflush(stdout);
fflush(stderr);
dup2(Logger::silentfd, fileno(stdout));
dup2(Logger::silentfd, fileno(stderr));
}
void StartTalking() {
// Check verbosity set correctly
if (!Logger::external_verb) return;
std::cout.rdbuf(Logger::default_cout);
std::cerr.rdbuf(Logger::default_cerr);
canihaznuisancepythia_();
fflush(stdout);
fflush(stderr);
dup2(Logger::savedstdoutfd, fileno(stdout));
dup2(Logger::savedstderrfd, fileno(stderr));
}
//******************************************
void LOG_VERB(std::string verb) {
//******************************************
if (!verb.compare("DEB"))
Logger::log_verb = -1;
else if (!verb.compare("QUIET"))
Logger::log_verb = 0;
else if (!verb.compare("FIT"))
Logger::log_verb = 1;
else if (!verb.compare("MIN"))
Logger::log_verb = 2;
else if (!verb.compare("SAM"))
Logger::log_verb = 3;
else if (!verb.compare("REC"))
Logger::log_verb = 4;
else if (!verb.compare("SIG"))
Logger::log_verb = 5;
else if (!verb.compare("EVT"))
Logger::log_verb = 6;
// else Logger::log_verb = GeneralUtils::StrToInt(verb);
std::cout << "Set logging verbosity to : " << Logger::log_verb << std::endl;
return;
}
//******************************************
void ERR_VERB(std::string verb) {
//******************************************
std::cout << "Setting ERROR VERB" << std::endl;
if (!verb.compare("ERRQUIET"))
Logger::err_verb = 0;
else if (!verb.compare("FTL"))
Logger::err_verb = 1;
else if (!verb.compare("WRN"))
Logger::err_verb = 2;
// else Logger::err_verb = GeneralUtils::StrToInt(verb);
std::cout << "Set error verbosity to : " << Logger::err_verb << std::endl;
return;
}
//******************************************
bool LOG_LEVEL(int level) {
//******************************************
if (Logger::log_verb == (unsigned int)DEB) {
return true;
}
- if (Logger::log_verb < (unsigned int)level) {
+ if (Logger::log_verb < (int)level) {
return false;
}
return true;
}
void SET_TRACE(bool val) { Logger::showtrace = val; }
//******************************************
std::ostream& _LOG(int level, const char* filename, const char* func, int line)
//******************************************
{
return __OUTLOG(level, filename, func, line);
}
//******************************************
std::ostream& _ERR(int level, const char* filename, const char* func, int line)
//******************************************
{
if (Logger::use_colors) std::cerr << RED;
if (Logger::showtrace) {
std::cout << filename << "::" << func << "[l. " << line << "] : ";
}
switch (level) {
case FTL:
std::cerr << "[ERR FATAL ]: ";
break;
case WRN:
std::cerr << "[ERR WARN ] : ";
break;
}
if (Logger::use_colors) std::cerr << RESET;
return *Logger::__ERR_outstream;
}
diff --git a/src/Reweight/MINERvAWeightCalcs.cxx b/src/Reweight/MINERvAWeightCalcs.cxx
index 5775b73..4afeecc 100644
--- a/src/Reweight/MINERvAWeightCalcs.cxx
+++ b/src/Reweight/MINERvAWeightCalcs.cxx
@@ -1,594 +1,594 @@
#ifdef __MINERVA_RW_ENABLED__
#ifdef __GENIE_ENABLED__
#include "MINERvAWeightCalcs.h"
#include "BaseFitEvt.h"
namespace nuisance {
namespace reweight {
//*******************************************************
MINERvAReWeight_QE::MINERvAReWeight_QE() {
//*******************************************************
fTwk_NormCCQE = 0.0;
fDef_NormCCQE = 1.0;
fCur_NormCCQE = fDef_NormCCQE;
}
//*******************************************************
MINERvAReWeight_QE::~MINERvAReWeight_QE(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_QE::CalcWeight(BaseFitEvt* evt) {
//*******************************************************
// Check GENIE
if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
const Interaction* interaction = ghep->Summary();
//const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
//const Target& tgt = init_state.Tgt();
// If the event is not QE this Calc doesn't handle it
if (!proc_info.IsQuasiElastic()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCQE Dial
if (!proc_info.IsWeakCC()) w *= fCur_NormCCQE;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_QE::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_QE::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
int curenum = rwenum % 1000;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCQE) {
fTwk_NormCCQE = val;
fCur_NormCCQE = fDef_NormCCQE + fTwk_NormCCQE;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCQE != 0.0));
}
//*******************************************************
bool MINERvAReWeight_QE::IsHandled(int rwenum) {
//*******************************************************
int curenum = rwenum % 1000;
switch (curenum) {
case Reweight::kMINERvARW_NormCCQE:
return true;
default:
return false;
}
}
//*******************************************************
MINERvAReWeight_MEC::MINERvAReWeight_MEC() {
//*******************************************************
fTwk_NormCCMEC = 0.0;
fDef_NormCCMEC = 1.0;
fCur_NormCCMEC = fDef_NormCCMEC;
}
//*******************************************************
MINERvAReWeight_MEC::~MINERvAReWeight_MEC(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_MEC::CalcWeight(BaseFitEvt* evt) {
//*******************************************************
// Check GENIE
if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
const Interaction* interaction = ghep->Summary();
//const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
//const Target& tgt = init_state.Tgt();
// If the event is not MEC this Calc doesn't handle it
if (!proc_info.IsMEC()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCMEC Dial
if (!proc_info.IsWeakCC()) w *= fCur_NormCCMEC;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_MEC::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_MEC::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
int curenum = rwenum % 1000;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCMEC) {
fTwk_NormCCMEC = val;
fCur_NormCCMEC = fDef_NormCCMEC + fTwk_NormCCMEC;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCMEC != 0.0));
}
//*******************************************************
bool MINERvAReWeight_MEC::IsHandled(int rwenum) {
//*******************************************************
int curenum = rwenum % 1000;
switch (curenum) {
case Reweight::kMINERvARW_NormCCMEC:
return true;
default:
return false;
}
}
//*******************************************************
MINERvAReWeight_RES::MINERvAReWeight_RES() {
//*******************************************************
fTwk_NormCCRES = 0.0;
fDef_NormCCRES = 1.0;
fCur_NormCCRES = fDef_NormCCRES;
}
//*******************************************************
MINERvAReWeight_RES::~MINERvAReWeight_RES(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_RES::CalcWeight(BaseFitEvt* evt) {
//*******************************************************
// std::cout << "Caculating RES" << std::endl;
// Check GENIE
if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
const Interaction* interaction = ghep->Summary();
//const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
//const Target& tgt = init_state.Tgt();
// If the event is not RES this Calc doesn't handle it
if (!proc_info.IsResonant()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCRES Dial
if (proc_info.IsWeakCC()) w *= fCur_NormCCRES;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_RES::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_RES::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
int curenum = rwenum % 1000;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCRES) {
fTwk_NormCCRES = val;
fCur_NormCCRES = fDef_NormCCRES + fTwk_NormCCRES;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCRES != 0.0));
}
//*******************************************************
bool MINERvAReWeight_RES::IsHandled(int rwenum) {
//*******************************************************
int curenum = rwenum % 1000;
switch (curenum) {
case Reweight::kMINERvARW_NormCCRES:
return true;
default:
return false;
}
}
//*******************************************************
RikRPA::RikRPA() {
//*******************************************************
// - Syst : kMINERvA_RikRPA_ApplyRPA
// - Type : Binary
// - Limits : 0.0 (false) -> 1.0 (true)
// - Default : 0.0
fApplyDial_RPACorrection = false;
// - Syst : kMINERvA_RikRPA_LowQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RPALowQ2 = 0.0;
fCurDial_RPALowQ2 = fDefDial_RPALowQ2;
fErrDial_RPALowQ2 = 0.0;
// - Syst : kMINERvA_RikRPA_HighQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RPAHighQ2 = 0.0;
fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2;
fErrDial_RPAHighQ2 = 1.0;
// - Syst : kMINERvA_RikRESRPA_ApplyRPA
// - Type : Binary
// - Limits : 0.0 (false) -> 1.0 (true)
// - Default : 0.0
fApplyDial_RESRPACorrection = false;
// - Syst : kMINERvA_RikRESRPA_LowQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RESRPALowQ2 = 0.0;
fCurDial_RESRPALowQ2 = fDefDial_RESRPALowQ2;
fErrDial_RESRPALowQ2 = 0.0;
// - Syst : kMINERvA_RikRESRPA_HighQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RESRPAHighQ2 = 0.0;
fCurDial_RESRPAHighQ2 = fDefDial_RESRPAHighQ2;
fErrDial_RESRPAHighQ2 = 1.0;
// Setup Calculators
fEventWeights = new double[5];
- for (size_t i = 0; i < kMaxCalculators; i++) {
+ for (int i = 0; i < kMaxCalculators; i++) {
fRPACalculators[i] = NULL;
}
fTweaked = false;
}
//*******************************************************
RikRPA::~RikRPA() {
//*******************************************************
// delete fEventWeights;
// for (size_t i = 0; i < kMaxCalculators; i++) {
// if (fRPACalculators[i]) delete fRPACalculators[i];
// fRPACalculators[i] = NULL;
// }
}
//*******************************************************
double RikRPA::CalcWeight(BaseFitEvt* evt) {
//*******************************************************
// LOG(FIT) << "Calculating RikRPA" << std::endl;
// Return 1.0 if not tweaked
if (!fTweaked) return 1.0;
double w = 1.0;
// Extract the GENIE Record
GHepRecord* ghep = static_cast<GHepRecord*>(evt->genie_event->event);
const Interaction* interaction = ghep->Summary();
const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
// const Kinematics & kine = interaction->Kine();
// const XclsTag & xcls = interaction->ExclTag();
const Target& tgt = init_state.Tgt();
// If not QE return 1.0
// LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() <<
// std::endl;
if (!tgt.IsNucleus()) return 1.0;
if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant()) return 1.0;
// Extract Beam and Target PDG
GHepParticle* neutrino = ghep->Probe();
int bpdg = neutrino->Pdg();
GHepParticle* target = ghep->Particle(1);
assert(target);
int tpdg = target->Pdg();
// Find the enum we need
int calcenum = GetRPACalcEnum(bpdg, tpdg);
if (calcenum == -1) return 1.0;
// Check we have the RPA Calc setup for this enum
// if not, set it up at that point
if (!fRPACalculators[calcenum]) SetupRPACalculator(calcenum);
weightRPA* rpacalc = fRPACalculators[calcenum];
if (!rpacalc) {
THROW("Failed to grab the RPA Calculator : " << calcenum);
}
// Extract Q0-Q3
GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
const TLorentzVector& k1 = *(neutrino->P4());
const TLorentzVector& k2 = *(fsl->P4());
double q0 = fabs((k1 - k2).E());
double q3 = fabs((k1 - k2).Vect().Mag());
double Q2 = fabs((k1 - k2).Mag2());
// Quasielastic
if (proc_info.IsQuasiElastic()){
// Now use q0-q3 and RPA Calculator to fill fWeights
rpacalc->getWeight(q0, q3, fEventWeights);
if (fApplyDial_RPACorrection) {
w *= fEventWeights[0]; // CV
}
// Syst Application : kMINERvA_RikRPA_LowQ2
if (fabs(fCurDial_RPALowQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RPALowQ2 > 0.0 && Q2 < 2.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[1]) *
fCurDial_RPALowQ2; // WLow+ } else if
} else if (fCurDial_RPALowQ2 < 0.0 && Q2 < 2.0) {
interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) *
fCurDial_RPALowQ2; // WLow-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
// Syst Application : kMINERvA_RikRPA_HighQ2
if (fabs(fCurDial_RPAHighQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RPAHighQ2 > 0.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
fCurDial_RPAHighQ2; // WHigh+
} else if (fCurDial_RPAHighQ2 < 0.0) {
interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
fCurDial_RPAHighQ2; // WHigh-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
}
// Resonant Events
if (proc_info.IsResonant()){
// Now use Q2 and RESRPA Calculator to fill fWeights
double CV = rpacalc->getWeight(Q2);
if (fApplyDial_RESRPACorrection) {
w *= CV; //fEventWeights[0]; // CVa
}
/*
// Syst Application : kMINERvA_RikRESRPA_LowQ2
if (fabs(fCurDial_RESRPAHighQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RESRPAHighQ2 > 0.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
fCurDial_RESRPAHighQ2; // WHigh+
} else if (fCurDial_RESRPAHighQ2 < 0.0) {
interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
fCurDial_RESRPAHighQ2; // WHigh-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
// Syst Application : kMINERvA_RikRESRPA_HighQ2
if (fabs(fCurDial_RESRPAHighQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RESRPAHighQ2 > 0.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
fCurDial_RESRPAHighQ2; // WHigh+
} else if (fCurDial_RESRPAHighQ2 < 0.0) {
interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
fCurDial_RESRPAHighQ2; // WHigh-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
*/
}
// LOG(FIT) << "RPA Weight = " << w << std::endl;
return w;
} // namespace reweight
//*******************************************************
void RikRPA::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void RikRPA::SetDialValue(int rwenum, double val) {
//*******************************************************
int curenum = rwenum % 1000;
// Check Handled
if (!IsHandled(curenum)) return;
if (curenum == Reweight::kMINERvARW_RikRPA_ApplyRPA) fApplyDial_RPACorrection = (val > 0.5);
if (curenum == Reweight::kMINERvARW_RikRPA_LowQ2) fCurDial_RPALowQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRPA_HighQ2) fCurDial_RPAHighQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRESRPA_ApplyRPA) fApplyDial_RESRPACorrection = (val > 0.5);
if (curenum == Reweight::kMINERvARW_RikRESRPA_LowQ2) fCurDial_RESRPALowQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRESRPA_HighQ2) fCurDial_RESRPAHighQ2 = val;
// Assign flag to say stuff has changed
fTweaked = (fApplyDial_RPACorrection ||
fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 ||
fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0 ||
fApplyDial_RESRPACorrection ||
fabs(fCurDial_RESRPAHighQ2 - fDefDial_RESRPAHighQ2) > 0.0 ||
fabs(fCurDial_RESRPALowQ2 - fDefDial_RESRPALowQ2) > 0.0);
}
//*******************************************************
bool RikRPA::IsHandled(int rwenum) {
//*******************************************************
int curenum = rwenum % 1000;
switch (curenum) {
case Reweight::kMINERvARW_RikRESRPA_ApplyRPA:
return true;
case Reweight::kMINERvARW_RikRESRPA_LowQ2:
return true;
case Reweight::kMINERvARW_RikRESRPA_HighQ2:
return true;
case Reweight::kMINERvARW_RikRPA_ApplyRPA:
return true;
case Reweight::kMINERvARW_RikRPA_LowQ2:
return true;
case Reweight::kMINERvARW_RikRPA_HighQ2:
return true;
default:
return false;
}
}
//*******************************************************
void RikRPA::SetupRPACalculator(int calcenum) {
//*******************************************************
std::string rwdir = FitPar::GetDataBase() + "/reweight/MINERvA/RikRPA/";
std::string fidir = "";
switch (calcenum) {
case kNuMuC12:
fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root";
break;
case kNuMuO16:
fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root";
break;
case kNuMuAr40:
fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root";
break;
case kNuMuCa40:
fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root";
break;
case kNuMuFe56:
fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root";
break;
case kNuMuBarC12:
fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root";
break;
case kNuMuBarO16:
fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root";
break;
case kNuMuBarAr40:
fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root";
break;
case kNuMuBarCa40:
fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root";
break;
case kNuMuBarFe56:
fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root";
break;
}
LOG(FIT) << "Loading RPA CALC : " << fidir << std::endl;
TDirectory* olddir = gDirectory;
std::cout << "***********************************************" << std::endl;
std::cout << "Loading a new weightRPA calculator" << std::endl;
std::cout << "Authors: Rik Gran, Heidi Schellman" << std::endl;
std::cout << "Citation: arXiv:1705.02932 [hep-ex]" << std::endl;
std::cout << "***********************************************" << std::endl;
fRPACalculators[calcenum] = new weightRPA(rwdir + "/" + fidir);
olddir->cd();
return;
}
//*******************************************************
int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) {
//*******************************************************
if (bpdg == 14 && tpdg == 1000060120)
return kNuMuC12;
else if (bpdg == 14 && tpdg == 1000080160)
return kNuMuO16;
else if (bpdg == 14 && tpdg == 1000180400)
return kNuMuAr40;
else if (bpdg == 14 && tpdg == 1000200400)
return kNuMuCa40;
else if (bpdg == 14 && tpdg == 1000280560)
return kNuMuFe56;
else if (bpdg == -14 && tpdg == 1000060120)
return kNuMuBarC12;
else if (bpdg == -14 && tpdg == 1000080160)
return kNuMuBarO16;
else if (bpdg == -14 && tpdg == 1000180400)
return kNuMuBarAr40;
else if (bpdg == -14 && tpdg == 1000200400)
return kNuMuBarCa40;
else if (bpdg == -14 && tpdg == 1000280560)
return kNuMuBarFe56;
else {
//ERROR(WRN, "Unknown beam and target combination for RPA Calcs! "
//<< bpdg << " " << tpdg);
}
return -1;
}
} // namespace reweight
} // namespace nuisance
#endif
#endif
diff --git a/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx b/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx
index 076d5a7..62aed1f 100644
--- a/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx
+++ b/src/T2K/T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx
@@ -1,198 +1,198 @@
// 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_CC0pinp_ifk_XSec_3Dinfip_nu.h"
//********************************************************************
T2K_CC0pinp_ifk_XSec_3Dinfip_nu::T2K_CC0pinp_ifk_XSec_3Dinfip_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pinp_ifk_XSec_3Dinfip_nu sample. \n" \
"Target: CH \n" \
"Flux: T2K 2.5 degree off-axis (ND280) \n" \
"Signal: CC0piNp (N>=1) with p_p>450MeV and cthp>0.4 \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("|#Delta p|");
fSettings.SetYTitle("p_#mu");
fSettings.SetZTitle("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");
fAnalysis = 1;
outOfBoundsMC = 0.0;
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pinp_ifk_XSec_3Dinfip_nu");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) / (TotalIntegratedFlux()));
//fScaleFactor = ((GetEventHistogram()->Integral("width")/(fNEvents+0.)) * 10 / (TotalIntegratedFlux()));
fSettings.SetDataInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root;result_tp" );
SetDataFromRootFile( fSettings.GetDataInput() );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root;cor_tp" );
SetCorrelationFromRootFile(fSettings.GetCovarInput() );
//SetCovarFromRootFile(FitPar::GetDataBase() + "/T2K/CC0pi/infkResults_origBin.root", "cov_tp" );
// Setup Histograms
SetHistograms();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
bool T2K_CC0pinp_ifk_XSec_3Dinfip_nu::isSignal(FitEvent *event){
return SignalDef::isT2K_CC0pi_ifk(event, EnuMin, EnuMax);
};
void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::FillEventVariables(FitEvent* event){
if (event->NumFSParticle(13) == 0 || event->NumFSParticle(2212) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector Pp = event->GetHMFSParticle(2212)->fP;
double pmu = Pmu.Vect().Mag()/1000.;
- double pp = Pp.Vect().Mag()/1000.;
+ //double pp = Pp.Vect().Mag()/1000.;
double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect()));
TVector3 tp_inf = FitUtils::tppInfK(Pmu, CosThetaMu, 25, true);
TVector3 Pp_mev(Pp.X()/1000,Pp.Y()/1000,Pp.Z()/1000);
TVector3 delta_tp = tp_inf-Pp_mev;
//std::cout << "Proton 3 mom is: " << std::endl;
//(Pp.Vect()).Print("all");
//std::cout << "Inferred Proton 3 mom is: " << std::endl;
//tp_inf.Print("all");
//std::cout << " " << std::endl;
fXVar = delta_tp.Mag();
fYVar = pmu;
fZVar = CosThetaMu;
return;
};
void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::FillHistograms(){
Measurement1D::FillHistograms();
if (Signal){
FillMCSlice( fXVar, fYVar, fZVar, Weight );
}
}
void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::ConvertEventRates(){
for (int i = 0; i < 7; i++){
fMCHist_Slices[i]->GetSumw2();
}
// Do standard conversion.
Measurement1D::ConvertEventRates();
// First scale MC slices also by their width in Y and Z
//MCHist_Slices[0]->Scale(1.0 / 1.00);
//MCHist_Slices[1]->Scale(1.0 / 0.60);
//MCHist_Slices[2]->Scale(1.0 / 0.10);
//MCHist_Slices[3]->Scale(1.0 / 0.10);
// Now Convert into 1D list
fMCHist->Reset();
//The first bin in the histogram in underflow, so set this and start bincount at 1
fMCHist->SetBinContent(1, outOfBoundsMC);
int bincount = 1;
for (int i = 0; i < 7; 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_CC0pinp_ifk_XSec_3Dinfip_nu::FillMCSlice(double x, double y, double z, double w){
// x is delta_tp
// y is pmu
// z is CosThetaMu
if (z <= -0.6) fMCHist_Slices[0]->Fill(x,w);
else if (z >= -0.6 and z < 0.0 and y < 0.25) fMCHist_Slices[1]->Fill(x,w);
else if (z >= -0.6 and z < 0.0 and y > 0.25) fMCHist_Slices[2]->Fill(x,w);
else if (z >= 0.0 and y < 0.25) fMCHist_Slices[3]->Fill(x,w);
else if (z >= 0.0 and z < 0.8 and y >= 0.25) fMCHist_Slices[4]->Fill(x,w);
else if (z >= 0.8 and z < 1.0 and y >= 0.25 and y < 0.75) fMCHist_Slices[5]->Fill(x,w);
else if (z >= 0.8 and z < 1.0 and y >= 0.75) fMCHist_Slices[6]->Fill(x,w);
else outOfBoundsMC += w;
}
void T2K_CC0pinp_ifk_XSec_3Dinfip_nu::SetHistograms(){
// Read in 1D Data Histograms
fInputFile = new TFile( (FitPar::GetDataBase() + "/T2K/CC0pi/STV/infkResults_origBin.root").c_str(),"READ");
//fInputFile->ls();
// Read in 1D Data
fDataHist = (TH1D*) fInputFile->Get("result_tp");
fDataHist->SetNameTitle("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data", "T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data");
SetAutoProcessTH1(fDataHist,kCMD_Write);
// Read in 2D Data Slices and Make MC Slices
for (int i = 0; i < 7; i++){ //both y and z slices
// Get Data Histogram
//fInputFile->ls();
fDataHist_Slices.push_back((TH1D*)fInputFile->Get(Form("resultBin%i_tp",i))->Clone());
fDataHist_Slices[i]->SetNameTitle(Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data_Slice%i",i), (Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_data_Slice%i",i)));
// Make MC Clones
fMCHist_Slices.push_back((TH1D*) fDataHist_Slices[i]->Clone());
fMCHist_Slices[i]->SetNameTitle(Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_MC_Slice%i",i), (Form("T2K_CC0pinp_ifk_XSec_3Dinfip_nu_MC_Slice%i",i)));
SetAutoProcessTH1(fDataHist_Slices[i],kCMD_Write);
SetAutoProcessTH1(fMCHist_Slices[i]);
fMCHist_Slices[i]->Reset();
}
return;
};
diff --git a/src/Tests/SmearceptanceTests.cxx b/src/Tests/SmearceptanceTests.cxx
index 8b1ed04..ad383fc 100644
--- a/src/Tests/SmearceptanceTests.cxx
+++ b/src/Tests/SmearceptanceTests.cxx
@@ -1,73 +1,73 @@
#include "SmearceptanceUtils.h"
#include "FitLogger.h"
#include "TDecompSVD.h"
#include <cassert>
#include <limits>
int main(int argc, char const *argv[]) {
double numerTol = 1.E2 * std::numeric_limits<double>::epsilon();
// a = 10, b = 4
// x = 5a + 3b (= 62)
// y = 6a - 2b (= 52)
// A = (10,4)
// M = ( (5, 3)
// (6, -2) )
// X = (62,52)
// A^T M = X
TVectorD A(2);
A[0] = 10;
A[1] = 4;
TVectorD X(2);
X[0] = 62;
X[1] = 52;
TMatrixD M(2, 2);
M[0][0] = 5;
M[0][1] = 3;
M[1][0] = 6;
M[1][1] = -2;
TVectorD ForwardSolve = M * A;
bool similar = true;
- for (size_t i = 0; i < A.GetNrows(); ++i) {
+ for (int i = 0; i < A.GetNrows(); ++i) {
if (fabs(X[i] - ForwardSolve[i]) > numerTol) {
ERROR(FTL, "Element " << i << " was not multiplied as expected: " << X[i]
<< " != " << ForwardSolve[i]
<< " (Tol: " << fabs(X[i] - ForwardSolve[i])
<< " > " << numerTol << ")");
similar = false;
}
}
assert(similar);
// TVectorD SVDSolve = SmearceptanceUtils::SVDInverseSolve(&A, &M);
TDecompSVD svd(M);
TMatrixD inv = svd.Invert();
TVectorD SVDSolve = inv * X;
similar = true;
- for (size_t i = 0; i < A.GetNrows(); ++i) {
+ for (int i = 0; i < A.GetNrows(); ++i) {
if (fabs(A[i] - SVDSolve[i]) > numerTol) {
ERROR(FTL, "Element " << i << " was not solved as expected: " << A[i]
<< " != " << SVDSolve[i]
<< " (Tol: " << fabs(A[i] - SVDSolve[i]) << " > "
<< numerTol << ")");
similar = false;
}
}
assert(similar);
}