diff --git a/CMakeLists.txt b/CMakeLists.txt
index aa07afe..590b423 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,223 +1,228 @@
# 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 .
################################################################################
cmake_minimum_required (VERSION 2.6 FATAL_ERROR)
project(NUISANCE)
include(ExternalProject)
enable_language(Fortran)
set (NUISANCE_VERSION_MAJOR 2)
set (NUISANCE_VERSION_MINOR 7)
set (NUISANCE_VERSION_REVISION 0)
set (NUISANCE_VERSION_STRING "v${NUISANCE_VERSION_MAJOR}r${NUISANCE_VERSION_MINOR}")
if(${NUISANCE_VERSION_REVISION} STRGREATER "0")
set (NUISANCE_VERSION_STRING "${NUISANCE_VERSION_STRING}p${NUISANCE_VERSION_REVISION}")
endif()
#Set this to TRUE to enable build debugging messages
set(BUILD_DEBUG_MSGS TRUE)
include(${CMAKE_SOURCE_DIR}/cmake/cmessage.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/cacheVariables.cmake)
cmessage(STATUS "CMAKE_INSTALL_PREFIX: \"${CMAKE_INSTALL_PREFIX}\"")
cmessage(STATUS "CMAKE_BUILD_TYPE: \"${CMAKE_BUILD_TYPE}\"")
################################################################################
# Check Dependencies
################################################################################
################################## ROOT ######################################
include(${CMAKE_SOURCE_DIR}/cmake/ROOTSetup.cmake)
################################# HEPMC ######################################
include(${CMAKE_SOURCE_DIR}/cmake/HepMC.cmake)
############################ Reweight Engines ################################
include(${CMAKE_SOURCE_DIR}/cmake/ReweightEnginesSetup.cmake)
############################ Other Generators ################################
include(${CMAKE_SOURCE_DIR}/cmake/GiBUUSetup.cmake)
if(USE_NUANCE)
LIST(APPEND EXTRA_CXX_FLAGS -D__NUANCE_ENABLED__)
endif()
################################# Pythia6/8 ####################################
include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/pythia8Setup.cmake)
################################# gperftools ###################################
include(${CMAKE_SOURCE_DIR}/cmake/gperfSetup.cmake)
if(NOT NOTEST)
enable_testing()
endif()
SET(GENERATOR_SUPPORT)
foreach(gen NEUT;NuWro;GENIE;GiBUU;NUANCE)
if(USE_${gen})
SET(GENERATOR_SUPPORT "${GENERATOR_SUPPORT}${gen} ")
endif()
endforeach(gen)
cmessage(STATUS "Generator Input Support: ${GENERATOR_SUPPORT}")
set(MINCODE
Routines
FCN)
set(CORE
MCStudies
Genie
FitBase
Config
Logger
InputHandler
Splines
Reweight
Utils
Statistical
#Devel
Smearceptance
)
LIST(APPEND ALLEXPERIMENTS
ANL
ArgoNeuT
BEBC
BNL
Electron
FNAL
GGM
K2K
MINERvA
MiniBooNE
SciBooNE
T2K)
foreach(exp ${ALLEXPERIMENTS})
if(NOT NO_${exp})
LIST(APPEND EXPERIMENTS_TO_BUILD ${exp})
else()
LIST(REVERSE EXTRA_CXX_FLAGS)
LIST(APPEND EXTRA_CXX_FLAGS -D__NO_${exp}__)
LIST(REVERSE EXTRA_CXX_FLAGS)
endif()
endforeach()
################################## COMPILER ####################################
include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake)
################################### doxygen ###################################
include(${CMAKE_SOURCE_DIR}/cmake/docsSetup.cmake)
################################################################################
set(MINIMUM_INCLUDE_DIRECTORIES)
LIST(APPEND MINIMUM_INCLUDE_DIRECTORIES
${RWENGINE_INCLUDE_DIRECTORIES}
${CMAKE_SOURCE_DIR}/src/FitBase
${CMAKE_SOURCE_DIR}/src/Reweight
${CMAKE_SOURCE_DIR}/src/InputHandler
${CMAKE_SOURCE_DIR}/src/Config
${CMAKE_SOURCE_DIR}/src/Logger
${CMAKE_SOURCE_DIR}/src/Statistical
${CMAKE_SOURCE_DIR}/src/Splines
${CMAKE_SOURCE_DIR}/src/Utils
${CMAKE_SOURCE_DIR}/src/Genie)
cmessage(DEBUG "Base include directories: ${MINIMUM_INCLUDE_DIRECTORIES}")
set(EXP_INCLUDE_DIRECTORIES)
foreach(edir ${EXPERIMENTS_TO_BUILD})
LIST(APPEND EXP_INCLUDE_DIRECTORIES ${CMAKE_SOURCE_DIR}/src/${edir})
endforeach()
cmessage(DEBUG "Included experiments: ${EXP_INCLUDE_DIRECTORIES}")
foreach(mdir ${MINCODE})
cmessage (DEBUG "Configuring directory: src/${mdir}")
add_subdirectory(src/${mdir})
endforeach()
foreach(edir ${EXPERIMENTS_TO_BUILD})
cmessage (DEBUG "Configuring directory: src/${edir}")
add_subdirectory(src/${edir})
endforeach()
foreach(cdir ${CORE})
cmessage (DEBUG "Configuring directory: src/${cdir}")
add_subdirectory(src/${cdir})
endforeach()
cmessage(DEBUG "Module targets: ${MODULETargets}")
add_subdirectory(app)
add_subdirectory(src/Tests)
configure_file(cmake/setup.sh.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" @ONLY)
install(FILES
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" DESTINATION
${CMAKE_INSTALL_PREFIX})
configure_file(cmake/MakeBinaryBlob.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/MakeBinaryBlob" @ONLY)
install(PROGRAMS
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/MakeBinaryBlob" DESTINATION
bin)
if(USE_DYNSAMPLES)
SET(ALL_INCLUDES ${MINIMUM_INCLUDE_DIRECTORIES})
LIST(APPEND ALL_INCLUDES ${CMAKE_SOURCE_DIR}/src/Smearceptance)
LIST(APPEND ALL_INCLUDES ${EXP_INCLUDE_DIRECTORIES})
string(REPLACE ";" " -I" ALL_INCLUDES_STR "${ALL_INCLUDES}")
+ cmessage(DEBUG ${CMAKE_DEPENDLIB_FLAGS})
+ string(REPLACE "-levent " "" CMAKE_DEPENDLIB_FLAGS_NEW ${CMAKE_DEPENDLIB_FLAGS})
+ set(CMAKE_DEPENDLIB_FLAGS ${CMAKE_DEPENDLIB_FLAGS_NEW})
+ cmessage(DEBUG ${CMAKE_DEPENDLIB_FLAGS})
+
string(REPLACE ";" " -l" ALL_MODULETARGETS_STR "${MODULETargets}")
configure_file(cmake/BuildDynamicSample.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSample" @ONLY)
install(PROGRAMS
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSample" DESTINATION
bin)
configure_file(cmake/BuildDynamicSmearcepter.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSmearcepter" @ONLY)
install(PROGRAMS
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSmearcepter" DESTINATION
bin)
endif()
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuiscardgen" DESTINATION
bin)
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuissamples" DESTINATION
bin)
diff --git a/src/Config/NuisConfig.cxx b/src/Config/NuisConfig.cxx
index f1f1033..7e78d66 100644
--- a/src/Config/NuisConfig.cxx
+++ b/src/Config/NuisConfig.cxx
@@ -1,802 +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 .
*******************************************************************************/
#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 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 cardlines =
GeneralUtils::ParseFileToStr(filename, "\n");
int linecount = 0;
// Loop over all input lines
for (std::vector::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 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 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 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 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 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() << "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 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 removed;
// Loop over all nodes and check for identical nodes
std::vector 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 nuisconfig::GetNodes(XMLNodePointer_t node,
std::string const &type) {
// Create new vector for nodes
std::vector 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 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 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 nuisconfig::GetVB(XMLNodePointer_t node,
// std::string name,
// const char* del) {
// std::string tempattr = Get(node, name);
// return GeneralUtils::ParseToBool(tempattr, del);
// }
std::vector 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 nuisconfig::GetVF(XMLNodePointer_t node,
// std::string name,
// const char* del) {
// std::string tempattr = Get(node, name);
// return GeneralUtils::ParseToDouble(tempattr, del);
// }
std::vector 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 nuisconfig::GetAllKeysForNode(XMLNodePointer_t node) {
bool matching = true;
XMLAttrPointer_t attr = fXML->GetFirstAttr(node);
std::vector 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 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/FCN/SampleList.cxx b/src/FCN/SampleList.cxx
index 857b9fb..36b1668 100644
--- a/src/FCN/SampleList.cxx
+++ b/src/FCN/SampleList.cxx
@@ -1,1307 +1,1307 @@
#include "SampleList.h"
#ifndef __NO_ANL__
#include "ANL_CCQE_Evt_1DQ2_nu.h"
#include "ANL_CCQE_XSec_1DEnu_nu.h"
// ANL CC1ppip
#include "ANL_CC1ppip_Evt_1DQ2_nu.h"
#include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1ppip_Evt_1DcosthAdler_nu.h"
#include "ANL_CC1ppip_Evt_1Dphi_nu.h"
#include "ANL_CC1ppip_Evt_1Dppi_nu.h"
#include "ANL_CC1ppip_Evt_1Dthpr_nu.h"
#include "ANL_CC1ppip_XSec_1DEnu_nu.h"
#include "ANL_CC1ppip_XSec_1DQ2_nu.h"
// ANL CC1npip
#include "ANL_CC1npip_Evt_1DQ2_nu.h"
#include "ANL_CC1npip_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1npip_Evt_1Dppi_nu.h"
#include "ANL_CC1npip_XSec_1DEnu_nu.h"
// ANL CC1pi0
#include "ANL_CC1pi0_Evt_1DQ2_nu.h"
#include "ANL_CC1pi0_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1pi0_XSec_1DEnu_nu.h"
// ANL NC1npip (mm, exotic!)
#include "ANL_NC1npip_Evt_1Dppi_nu.h"
// ANL NC1ppim (mm, exotic!)
#include "ANL_NC1ppim_Evt_1DcosmuStar_nu.h"
#include "ANL_NC1ppim_XSec_1DEnu_nu.h"
// ANL CC2pi 1pim1pip (mm, even more exotic!)
#include "ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu.h"
#include "ANL_CC2pi_1pim1pip_Evt_1Dppim_nu.h"
#include "ANL_CC2pi_1pim1pip_Evt_1Dppip_nu.h"
#include "ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu.h"
#include "ANL_CC2pi_1pim1pip_XSec_1DEnu_nu.h"
// ANL CC2pi 1pip1pip (mm, even more exotic!)
#include "ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu.h"
#include "ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu.h"
#include "ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu.h"
#include "ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu.h"
#include "ANL_CC2pi_1pip1pip_XSec_1DEnu_nu.h"
// ANL CC2pi 1pip1pi0 (mm, even more exotic!)
#include "ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu.h"
#include "ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu.h"
#include "ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu.h"
#include "ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu.h"
#include "ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu.h"
#endif
#ifndef __NO_ArgoNeuT__
// ArgoNeuT CC-inclusive
#include "ArgoNeuT_CCInc_XSec_1Dpmu_antinu.h"
#include "ArgoNeuT_CCInc_XSec_1Dpmu_nu.h"
#include "ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.h"
#include "ArgoNeuT_CCInc_XSec_1Dthetamu_nu.h"
#endif
#ifndef __NO_BNL__
// BNL CCQE
#include "BNL_CCQE_Evt_1DQ2_nu.h"
#include "BNL_CCQE_XSec_1DEnu_nu.h"
// BNL CC1ppip
#include "BNL_CC1ppip_Evt_1DQ2_nu.h"
#include "BNL_CC1ppip_Evt_1DQ2_nu.h"
#include "BNL_CC1ppip_Evt_1DcosthAdler_nu.h"
#include "BNL_CC1ppip_Evt_1Dphi_nu.h"
#include "BNL_CC1ppip_XSec_1DEnu_nu.h"
// BNL CC1npip
#include "BNL_CC1npip_Evt_1DQ2_nu.h"
#include "BNL_CC1npip_XSec_1DEnu_nu.h"
// BNL CC1pi0
#include "BNL_CC1pi0_Evt_1DQ2_nu.h"
#include "BNL_CC1pi0_XSec_1DEnu_nu.h"
#endif
#ifndef __NO_FNAL__
// FNAL CCQE
#include "FNAL_CCQE_Evt_1DQ2_nu.h"
// FNAL CC1ppip
#include "FNAL_CC1ppip_Evt_1DQ2_nu.h"
#include "FNAL_CC1ppip_XSec_1DEnu_nu.h"
#include "FNAL_CC1ppip_XSec_1DQ2_nu.h"
// FNAL CC1ppim
#include "FNAL_CC1ppim_XSec_1DEnu_antinu.h"
#endif
#ifndef __NO_BEBC__
// BEBC CCQE
#include "BEBC_CCQE_XSec_1DQ2_nu.h"
// BEBC CC1ppip
#include "BEBC_CC1ppip_XSec_1DEnu_nu.h"
#include "BEBC_CC1ppip_XSec_1DQ2_nu.h"
// BEBC CC1npip
#include "BEBC_CC1npip_XSec_1DEnu_nu.h"
#include "BEBC_CC1npip_XSec_1DQ2_nu.h"
// BEBC CC1pi0
#include "BEBC_CC1pi0_XSec_1DEnu_nu.h"
#include "BEBC_CC1pi0_XSec_1DQ2_nu.h"
// BEBC CC1npim
#include "BEBC_CC1npim_XSec_1DEnu_antinu.h"
#include "BEBC_CC1npim_XSec_1DQ2_antinu.h"
// BEBC CC1ppim
#include "BEBC_CC1ppim_XSec_1DEnu_antinu.h"
#include "BEBC_CC1ppim_XSec_1DQ2_antinu.h"
#endif
#ifndef __NO_GGM__
// GGM CC1ppip
#include "GGM_CC1ppip_Evt_1DQ2_nu.h"
#include "GGM_CC1ppip_XSec_1DEnu_nu.h"
#endif
#ifndef __NO_MiniBooNE__
// MiniBooNE CCQE
#include "MiniBooNE_CCQE_XSec_1DQ2_antinu.h"
#include "MiniBooNE_CCQE_XSec_1DQ2_nu.h"
#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h"
#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h"
#include "MiniBooNE_CCQE_XSec_2DTcos_nu.h"
// MiniBooNE CC1pi+ 1D
#include "MiniBooNE_CC1pip_XSec_1DEnu_nu.h"
#include "MiniBooNE_CC1pip_XSec_1DQ2_nu.h"
#include "MiniBooNE_CC1pip_XSec_1DTpi_nu.h"
#include "MiniBooNE_CC1pip_XSec_1DTu_nu.h"
// MiniBooNE CC1pi+ 2D
#include "MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTpiCospi_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTpiEnu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTuCosmu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTuEnu_nu.h"
// MiniBooNE CC1pi0
#include "MiniBooNE_CC1pi0_XSec_1DEnu_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1DQ2_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1DTu_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1Dppi0_nu.h"
#include "MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu.h"
#include "MiniBooNE_NC1pi0_XSec_1Dcospi0_nu.h"
#include "MiniBooNE_NC1pi0_XSec_1Dppi0_antinu.h"
#include "MiniBooNE_NC1pi0_XSec_1Dppi0_nu.h"
// MiniBooNE NC1pi0
//#include "MiniBooNE_NCpi0_XSec_1Dppi0_nu.h"
// MiniBooNE NCEL
#include "MiniBooNE_NCEL_XSec_Treco_nu.h"
#endif
#ifndef __NO_MINERvA__
// MINERvA CCQE
#include "MINERvA_CCQE_XSec_1DQ2_antinu.h"
#include "MINERvA_CCQE_XSec_1DQ2_joint.h"
#include "MINERvA_CCQE_XSec_1DQ2_nu.h"
// MINERvA CC0pi
#include "MINERvA_CC0pi_XSec_1DEe_nue.h"
#include "MINERvA_CC0pi_XSec_1DQ2_nu_proton.h"
#include "MINERvA_CC0pi_XSec_1DQ2_nue.h"
#include "MINERvA_CC0pi_XSec_1DThetae_nue.h"
// MINERvA CC1pi+
#include "MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h"
#include "MINERvA_CC1pip_XSec_1DTpi_nu.h"
#include "MINERvA_CC1pip_XSec_1Dth_20deg_nu.h"
#include "MINERvA_CC1pip_XSec_1Dth_nu.h"
// 2017 data update
#include "MINERvA_CC1pip_XSec_1D_2017Update.h"
// MINERvA CCNpi+
#include "MINERvA_CCNpip_XSec_1DEnu_nu.h"
#include "MINERvA_CCNpip_XSec_1DQ2_nu.h"
#include "MINERvA_CCNpip_XSec_1DTpi_nu.h"
#include "MINERvA_CCNpip_XSec_1Dpmu_nu.h"
#include "MINERvA_CCNpip_XSec_1Dth_nu.h"
#include "MINERvA_CCNpip_XSec_1Dthmu_nu.h"
// MINERvA CC1pi0
#include "MINERvA_CC1pi0_XSec_1DEnu_antinu.h"
#include "MINERvA_CC1pi0_XSec_1DQ2_antinu.h"
#include "MINERvA_CC1pi0_XSec_1DTpi0_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dpmu_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dppi0_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dth_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dthmu_antinu.h"
// MINERvA CC1pi0 neutrino
#include "MINERvA_CC1pi0_XSec_1D_nu.h"
// MINERvA CCINC
#include "MINERvA_CCinc_XSec_1DEnu_ratio.h"
#include "MINERvA_CCinc_XSec_1Dx_ratio.h"
#include "MINERvA_CCinc_XSec_2DEavq3_nu.h"
// MINERvA CCDIS
#include "MINERvA_CCDIS_XSec_1DEnu_ratio.h"
#include "MINERvA_CCDIS_XSec_1Dx_ratio.h"
// MINERvA CCCOH pion
#include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h"
#include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h"
#include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h"
#include "MINERvA_CCCOHPI_XSec_1DQ2_antinu.h"
#include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h"
#include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h"
#include "MINERvA_CCCOHPI_XSec_1Dth_nu.h"
#include "MINERvA_CCCOHPI_XSec_1Dth_nu.h"
#include "MINERvA_CCCOHPI_XSec_joint.h"
#include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h"
#include "MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h"
#include "MINERvA_CC0pi_XSec_2Dptpx_antinu.h"
#include "MINERvA_CC0pi_XSec_2Dptpx_nu.h"
#endif
#ifndef __NO_T2K__
// T2K CC0pi
#include "T2K_CC0pi_XSec_2DPcos_nu.h"
// T2K CC1pi+ on CH
#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h"
#include "T2K_CC1pip_CH_XSec_1DWrec_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dpmu_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dppi_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dq3_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h"
// T2K CC1pi+ on H2O
#include "T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h"
#include "T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dcospi_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dpmu_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dppi_nu.h"
// T2K STV CC0pi
#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h"
#include "T2K_CC0pi_XSec_2DPcos_nu_nonuniform.h"
#endif
#ifndef __NO_SciBooNE__
// SciBooNE COH studies
#include "SciBooNE_CCCOH_1TRK_1DQ2_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu.h"
#include "SciBooNE_CCCOH_MuPiVA_1DQ2_nu.h"
#include "SciBooNE_CCCOH_MuPr_1DQ2_nu.h"
#include "SciBooNE_CCCOH_STOPFINAL_1DQ2_nu.h"
#include "SciBooNE_CCCOH_STOP_NTrks_nu.h"
#endif
#ifndef __NO_K2K__
// K2K NC1pi0
#include "K2K_NC1pi0_Evt_1Dppi0_nu.h"
#endif
// MC Studies
#include "ExpMultDist_CCQE_XSec_1DVar_FakeStudy.h"
#include "ExpMultDist_CCQE_XSec_2DVar_FakeStudy.h"
#include "MCStudy_CCQEHistograms.h"
#include "GenericFlux_Tester.h"
#include "GenericFlux_Vectors.h"
#include "ElectronFlux_FlatTree.h"
#include "ElectronScattering_DurhamData.h"
#include "MCStudy_KaonPreSelection.h"
#include "MCStudy_MuonValidation.h"
#include "OfficialNIWGPlots.h"
#include "T2K2017_FakeData.h"
#include "Simple_Osc.h"
#include "Smear_SVDUnfold_Propagation_Osc.h"
#include "FitWeight.h"
#include "NuisConfig.h"
#include "NuisKey.h"
#ifdef __USE_DYNSAMPLES__
#include "TRegexp.h"
#include
// linux
#include
DynamicSampleFactory::DynamicSampleFactory() : NSamples(0), NManifests(0) {
LoadPlugins();
QLOG(FIT, "Loaded " << NSamples << " from " << NManifests
<< " shared object libraries.");
}
DynamicSampleFactory* DynamicSampleFactory::glblDSF = NULL;
DynamicSampleFactory::PluginManifest::~PluginManifest() {
for (size_t i_it = 0; i_it < Instances.size(); ++i_it) {
(*(DSF_DestroySample))(Instances[i_it]);
}
}
std::string EnsureTrailingSlash(std::string const& inp) {
if (!inp.length()) {
return "/";
}
if (inp[inp.length() - 1] == '/') {
return inp;
}
return inp + "/";
}
void DynamicSampleFactory::LoadPlugins() {
std::vector SearchDirectories;
if (Config::HasPar("dynamic_sample.path")) {
SearchDirectories =
GeneralUtils::ParseToStr(Config::GetParS("dynamic_sample.path"), ":");
}
char const* envPath = getenv("NUISANCE_DS_PATH");
if (envPath) {
std::vector envPaths = GeneralUtils::ParseToStr(envPath, ":");
for (size_t ep_it = 0; ep_it < envPaths.size(); ++ep_it) {
SearchDirectories.push_back(envPaths[ep_it]);
}
}
if (!SearchDirectories.size()) {
char const* pwdPath = getenv("PWD");
if (pwdPath) {
SearchDirectories.push_back(pwdPath);
}
}
for (size_t sp_it = 0; sp_it < SearchDirectories.size(); ++sp_it) {
std::string dirpath = EnsureTrailingSlash(SearchDirectories[sp_it]);
QLOG(FIT, "Searching for dynamic sample manifests in: " << dirpath);
Ssiz_t len = 0;
DIR* dir;
struct dirent* ent;
dir = opendir(dirpath.c_str());
if (dir != NULL) {
TRegexp matchExp("*.so", true);
while ((ent = readdir(dir)) != NULL) {
if (matchExp.Index(TString(ent->d_name), &len) != Ssiz_t(-1)) {
QLOG(FIT, "\tFound shared object: "
<< ent->d_name << " checking for relevant methods...");
void* dlobj =
dlopen((dirpath + ent->d_name).c_str(), RTLD_NOW | RTLD_GLOBAL);
char const* dlerr_cstr = dlerror();
std::string dlerr;
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN, "\tDL Load Error: " << dlerr);
continue;
}
PluginManifest plgManif;
plgManif.dllib = dlobj;
plgManif.soloc = (dirpath + ent->d_name);
plgManif.DSF_NSamples =
reinterpret_cast(dlsym(dlobj, "DSF_NSamples"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN, "\tFailed to load symbol \"DSF_NSamples\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.DSF_GetSampleName = reinterpret_cast(
dlsym(dlobj, "DSF_GetSampleName"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN, "\tFailed to load symbol \"DSF_GetSampleName\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.DSF_GetSample = reinterpret_cast(
dlsym(dlobj, "DSF_GetSample"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN, "\tFailed to load symbol \"DSF_GetSample\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.DSF_DestroySample = reinterpret_cast(
dlsym(dlobj, "DSF_DestroySample"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN, "Failed to load symbol \"DSF_DestroySample\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.NSamples = (*(plgManif.DSF_NSamples))();
QLOG(FIT, "\tSuccessfully loaded dynamic sample manifest: "
<< plgManif.soloc << ". Contains " << plgManif.NSamples
<< " samples.");
for (size_t smp_it = 0; smp_it < plgManif.NSamples; ++smp_it) {
char const* smp_name = (*(plgManif.DSF_GetSampleName))(smp_it);
if (!smp_name) {
THROW("Could not load sample " << smp_it << " / "
<< plgManif.NSamples << " from "
<< plgManif.soloc);
}
if (Samples.count(smp_name)) {
ERROR(WRN, "Already loaded a sample named: \""
<< smp_name << "\". cannot load duplciates. This "
"sample will be skipped.");
continue;
}
plgManif.SamplesProvided.push_back(smp_name);
Samples[smp_name] = std::make_pair(plgManif.soloc, smp_it);
QLOG(FIT, "\t\t" << smp_name);
}
if (plgManif.SamplesProvided.size()) {
Manifests[plgManif.soloc] = plgManif;
NSamples += plgManif.SamplesProvided.size();
NManifests++;
} else {
dlclose(dlobj);
}
}
}
closedir(dir);
} else {
ERROR(WRN, "Tried to open non-existant directory.");
}
}
}
DynamicSampleFactory& DynamicSampleFactory::Get() {
if (!glblDSF) {
glblDSF = new DynamicSampleFactory();
}
return *glblDSF;
}
void DynamicSampleFactory::Print() {
std::map > ManifestSamples;
for (std::map >::iterator smp_it =
Samples.begin();
smp_it != Samples.end(); ++smp_it) {
if (!ManifestSamples.count(smp_it->second.first)) {
ManifestSamples[smp_it->second.first] = std::vector();
}
ManifestSamples[smp_it->second.first].push_back(smp_it->first);
}
QLOG(FIT, "Dynamic sample manifest: ");
for (std::map >::iterator m_it =
ManifestSamples.begin();
m_it != ManifestSamples.end(); ++m_it) {
QLOG(FIT, "\tLibrary " << m_it->first << " contains: ");
for (size_t s_it = 0; s_it < m_it->second.size(); ++s_it) {
QLOG(FIT, "\t\t" << m_it->second[s_it]);
}
}
}
bool DynamicSampleFactory::HasSample(std::string const& name) {
return Samples.count(name);
}
bool DynamicSampleFactory::HasSample(nuiskey& samplekey) {
return HasSample(samplekey.GetS("name"));
}
MeasurementBase* DynamicSampleFactory::CreateSample(nuiskey& samplekey) {
if (!HasSample(samplekey)) {
ERROR(WRN, "Asked to load unknown sample: \"" << samplekey.GetS("name")
<< "\".");
return NULL;
}
std::pair sample = Samples[samplekey.GetS("name")];
QLOG(SAM, "\tLoading sample " << sample.second << " from " << sample.first);
return (*(Manifests[sample.first].DSF_GetSample))(sample.second, &samplekey);
}
DynamicSampleFactory::~DynamicSampleFactory() { Manifests.clear(); }
#endif
//! Functions to make it easier for samples to be created and handled.
namespace SampleUtils {
//! Create a given sample given its name, file, type, fakdata(fkdt) file and the
//! current rw engine and push it back into the list fChain.
MeasurementBase* CreateSample(std::string name, std::string file,
std::string type, std::string fkdt,
FitWeight* rw) {
nuiskey samplekey = Config::CreateKey("sample");
samplekey.Set("name", name);
samplekey.Set("input", file);
samplekey.Set("type", type);
return CreateSample(samplekey);
}
MeasurementBase* CreateSample(nuiskey samplekey) {
#ifdef __USE_DYNSAMPLES__
if (DynamicSampleFactory::Get().HasSample(samplekey)) {
QLOG(SAM, "Instantiating dynamic sample...");
MeasurementBase* ds = DynamicSampleFactory::Get().CreateSample(samplekey);
if (ds) {
QLOG(SAM, "Done.");
return ds;
}
THROW("Failed to instantiate dynamic sample.");
}
#endif
FitWeight* rw = FitBase::GetRW();
std::string name = samplekey.GetS("name");
std::string file = samplekey.GetS("input");
std::string type = samplekey.GetS("type");
std::string fkdt = "";
/*
ANL CCQE Samples
*/
#ifndef __NO_ANL__
if (!name.compare("ANL_CCQE_XSec_1DEnu_nu") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD26") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRL31") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD16")) {
return (new ANL_CCQE_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CCQE_Evt_1DQ2_nu") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRL31") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD26") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD16")) {
return (new ANL_CCQE_Evt_1DQ2_nu(samplekey));
/*
ANL CC1ppip samples
*/
} else if (!name.compare("ANL_CC1ppip_XSec_1DEnu_nu") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_Uncorr") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W16Cut_Uncorr")) {
return (new ANL_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_XSec_1DQ2_nu")) {
return (new ANL_CC1ppip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DQ2_nu") ||
!name.compare("ANL_CC1ppip_Evt_1DQ2_nu_W14Cut")) {
return (new ANL_CC1ppip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dppi_nu")) {
return (new ANL_CC1ppip_Evt_1Dppi_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dthpr_nu")) {
return (new ANL_CC1ppip_Evt_1Dthpr_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DcosmuStar_nu")) {
return (new ANL_CC1ppip_Evt_1DcosmuStar_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DcosthAdler_nu")) {
return (new ANL_CC1ppip_Evt_1DcosthAdler_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dphi_nu")) {
return (new ANL_CC1ppip_Evt_1Dphi_nu(samplekey));
/*
ANL CC1npip sample
*/
} else if (!name.compare("ANL_CC1npip_XSec_1DEnu_nu") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_Uncorr") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut_Uncorr") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_W16Cut_Uncorr")) {
return (new ANL_CC1npip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1DQ2_nu") ||
!name.compare("ANL_CC1npip_Evt_1DQ2_nu_W14Cut")) {
return (new ANL_CC1npip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1Dppi_nu")) {
return (new ANL_CC1npip_Evt_1Dppi_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1DcosmuStar_nu")) {
return (new ANL_CC1npip_Evt_1DcosmuStar_nu(samplekey));
/*
ANL CC1pi0 sample
*/
} else if (!name.compare("ANL_CC1pi0_XSec_1DEnu_nu") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_Uncorr") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut_Uncorr") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W16Cut_Uncorr")) {
return (new ANL_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC1pi0_Evt_1DQ2_nu") ||
!name.compare("ANL_CC1pi0_Evt_1DQ2_nu_W14Cut")) {
return (new ANL_CC1pi0_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1pi0_Evt_1DcosmuStar_nu")) {
return (new ANL_CC1pi0_Evt_1DcosmuStar_nu(samplekey));
/*
ANL NC1npip sample
*/
} else if (!name.compare("ANL_NC1npip_Evt_1Dppi_nu")) {
return (new ANL_NC1npip_Evt_1Dppi_nu(samplekey));
/*
ANL NC1ppim sample
*/
} else if (!name.compare("ANL_NC1ppim_XSec_1DEnu_nu")) {
return (new ANL_NC1ppim_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_NC1ppim_Evt_1DcosmuStar_nu")) {
return (new ANL_NC1ppim_Evt_1DcosmuStar_nu(samplekey));
/*
ANL CC2pi sample
*/
} else if (!name.compare("ANL_CC2pi_1pim1pip_XSec_1DEnu_nu")) {
return (new ANL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppip_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dppip_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppim_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dppim_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_XSec_1DEnu_nu")) {
return (new ANL_CC2pi_1pip1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu")) {
return (new ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu(samplekey));
/*
ArgoNeut Samples
*/
} else
#endif
#ifndef __NO_ArgoNeuT__
if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_antinu")) {
return (new ArgoNeuT_CCInc_XSec_1Dpmu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_nu")) {
return (new ArgoNeuT_CCInc_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_antinu")) {
return (new ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_nu")) {
return (new ArgoNeuT_CCInc_XSec_1Dthetamu_nu(samplekey));
/*
BNL Samples
*/
} else
#endif
#ifndef __NO_BNL__
if (!name.compare("BNL_CCQE_XSec_1DEnu_nu")) {
return (new BNL_CCQE_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CCQE_Evt_1DQ2_nu")) {
return (new BNL_CCQE_Evt_1DQ2_nu(samplekey));
/*
BNL CC1ppip samples
*/
} else if (!name.compare("BNL_CC1ppip_XSec_1DEnu_nu") ||
!name.compare("BNL_CC1ppip_XSec_1DEnu_nu_Uncorr") ||
!name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut") ||
!name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr")) {
return (new BNL_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1DQ2_nu") ||
!name.compare("BNL_CC1ppip_Evt_1DQ2_nu_W14Cut")) {
return (new BNL_CC1ppip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1DcosthAdler_nu")) {
return (new BNL_CC1ppip_Evt_1DcosthAdler_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1Dphi_nu")) {
return (new BNL_CC1ppip_Evt_1Dphi_nu(samplekey));
/*
BNL CC1npip samples
*/
} else if (!name.compare("BNL_CC1npip_XSec_1DEnu_nu") ||
!name.compare("BNL_CC1npip_XSec_1DEnu_nu_Uncorr")) {
return (new BNL_CC1npip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC1npip_Evt_1DQ2_nu")) {
return (new BNL_CC1npip_Evt_1DQ2_nu(samplekey));
/*
BNL CC1pi0 samples
*/
} else if (!name.compare("BNL_CC1pi0_XSec_1DEnu_nu")) {
return (new BNL_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC1pi0_Evt_1DQ2_nu")) {
return (new BNL_CC1pi0_Evt_1DQ2_nu(samplekey));
/*
FNAL Samples
*/
} else
#endif
#ifndef __NO_FNAL__
if (!name.compare("FNAL_CCQE_Evt_1DQ2_nu")) {
return (new FNAL_CCQE_Evt_1DQ2_nu(samplekey));
/*
FNAL CC1ppip
*/
} else if (!name.compare("FNAL_CC1ppip_XSec_1DEnu_nu")) {
return (new FNAL_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("FNAL_CC1ppip_XSec_1DQ2_nu")) {
return (new FNAL_CC1ppip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("FNAL_CC1ppip_Evt_1DQ2_nu")) {
return (new FNAL_CC1ppip_Evt_1DQ2_nu(samplekey));
/*
FNAL CC1ppim
*/
} else if (!name.compare("FNAL_CC1ppim_XSec_1DEnu_antinu")) {
return (new FNAL_CC1ppim_XSec_1DEnu_antinu(samplekey));
/*
BEBC Samples
*/
} else
#endif
#ifndef __NO_BEBC__
if (!name.compare("BEBC_CCQE_XSec_1DQ2_nu")) {
return (new BEBC_CCQE_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1ppip samples
*/
} else if (!name.compare("BEBC_CC1ppip_XSec_1DEnu_nu")) {
return (new BEBC_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BEBC_CC1ppip_XSec_1DQ2_nu")) {
return (new BEBC_CC1ppip_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1npip samples
*/
} else if (!name.compare("BEBC_CC1npip_XSec_1DEnu_nu")) {
return (new BEBC_CC1npip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BEBC_CC1npip_XSec_1DQ2_nu")) {
return (new BEBC_CC1npip_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1pi0 samples
*/
} else if (!name.compare("BEBC_CC1pi0_XSec_1DEnu_nu")) {
return (new BEBC_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BEBC_CC1pi0_XSec_1DQ2_nu")) {
return (new BEBC_CC1pi0_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1npim samples
*/
} else if (!name.compare("BEBC_CC1npim_XSec_1DEnu_antinu")) {
return (new BEBC_CC1npim_XSec_1DEnu_antinu(samplekey));
} else if (!name.compare("BEBC_CC1npim_XSec_1DQ2_antinu")) {
return (new BEBC_CC1npim_XSec_1DQ2_antinu(samplekey));
/*
BEBC CC1ppim samples
*/
} else if (!name.compare("BEBC_CC1ppim_XSec_1DEnu_antinu")) {
return (new BEBC_CC1ppim_XSec_1DEnu_antinu(samplekey));
} else if (!name.compare("BEBC_CC1ppim_XSec_1DQ2_antinu")) {
return (new BEBC_CC1ppim_XSec_1DQ2_antinu(samplekey));
/*
GGM CC1ppip samples
*/
} else
#endif
#ifndef __NO_GGM__
if (!name.compare("GGM_CC1ppip_XSec_1DEnu_nu")) {
return (new GGM_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("GGM_CC1ppip_Evt_1DQ2_nu")) {
return (new GGM_CC1ppip_Evt_1DQ2_nu(samplekey));
/*
MiniBooNE Samples
*/
/*
CCQE
*/
} else
#endif
#ifndef __NO_MiniBooNE__
if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_nu") ||
!name.compare("MiniBooNE_CCQELike_XSec_1DQ2_nu")) {
return (new MiniBooNE_CCQE_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_antinu") ||
!name.compare("MiniBooNE_CCQELike_XSec_1DQ2_antinu") ||
!name.compare("MiniBooNE_CCQE_CTarg_XSec_1DQ2_antinu")) {
return (new MiniBooNE_CCQE_XSec_1DQ2_antinu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_nu") ||
!name.compare("MiniBooNE_CCQELike_XSec_2DTcos_nu")) {
return (new MiniBooNE_CCQE_XSec_2DTcos_nu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_antinu") ||
!name.compare("MiniBooNE_CCQELike_XSec_2DTcos_antinu")) {
return (new MiniBooNE_CCQE_XSec_2DTcos_antinu(samplekey));
/*
MiniBooNE CC1pi+
*/
// 1D
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DEnu_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DQ2_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTpi_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DTpi_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTu_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DTu_nu(samplekey));
// 2D
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DQ2Enu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DQ2Enu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiCospi_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTpiCospi_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiEnu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTpiEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuCosmu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTuCosmu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuEnu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTuEnu_nu(samplekey));
/*
MiniBooNE CC1pi0
*/
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DEnu_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DQ2_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DTu_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1DTu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcosmu_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1Dcosmu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcospi0_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1Dcospi0_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dppi0_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1Dppi0_nu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_rhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_nu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_fhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_nu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_antinu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_rhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dppi0_antinu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_nu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_fhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dppi0_nu(samplekey));
/*
MiniBooNE NCEL
*/
} else if (!name.compare("MiniBooNE_NCEL_XSec_Treco_nu")) {
return (new MiniBooNE_NCEL_XSec_Treco_nu(samplekey));
/*
MINERvA Samples
*/
} else
#endif
#ifndef __NO_MINERvA__
if (!name.compare("MINERvA_CCQE_XSec_1DQ2_nu") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg_oldflux")) {
return (new MINERvA_CCQE_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg_oldflux")) {
return (new MINERvA_CCQE_XSec_1DQ2_antinu(samplekey));
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg")) {
return (new MINERvA_CCQE_XSec_1DQ2_joint(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DEe_nue")) {
return (new MINERvA_CC0pi_XSec_1DEe_nue(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nue")) {
return (new MINERvA_CC0pi_XSec_1DQ2_nue(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DThetae_nue")) {
return (new MINERvA_CC0pi_XSec_1DThetae_nue(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nu_proton")) {
return (new MINERvA_CC0pi_XSec_1DQ2_nu_proton(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtC_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtCH_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtFe_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtPb_nu")) {
return (new MINERvA_CC0pi_XSec_1DQ2_Tgt_nu(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioC_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioFe_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioPb_nu")) {
return (new MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpx_nu")) {
return (new MINERvA_CC0pi_XSec_2Dptpx_nu(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpx_antinu")) {
return (new MINERvA_CC0pi_XSec_2Dptpx_antinu(samplekey));
/*
CC1pi+
*/
// DONE
} else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu") ||
!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg") ||
!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_fluxcorr") ||
!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg_fluxcorr")) {
return (new MINERvA_CC1pip_XSec_1DTpi_nu(samplekey));
// DONE
} else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_nu") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_fluxcorr") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg_fluxcorr")) {
return (new MINERvA_CC1pip_XSec_1Dth_nu(samplekey));
} else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1Dpmu_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1Dthmu_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1DQ2_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1DEnu_nu_2017")) {
return (new MINERvA_CC1pip_XSec_1D_2017Update(samplekey));
/*
CCNpi+
*/
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dth_nu") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2016") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_fluxcorr") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg_fluxcorr")) {
return (new MINERvA_CCNpip_XSec_1Dth_nu(samplekey));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2016") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_fluxcorr") ||
!name.compare(
"MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg_fluxcorr")) {
return (new MINERvA_CCNpip_XSec_1DTpi_nu(samplekey));
// Done
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) {
return (new MINERvA_CCNpip_XSec_1Dthmu_nu(samplekey));
// Done
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) {
return (new MINERvA_CCNpip_XSec_1Dpmu_nu(samplekey));
// Done
} else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) {
return (new MINERvA_CCNpip_XSec_1DQ2_nu(samplekey));
// Done
} else if (!name.compare("MINERvA_CCNpip_XSec_1DEnu_nu")) {
return (new MINERvA_CCNpip_XSec_1DEnu_nu(samplekey));
/*
MINERvA CC1pi0 anti-nu
*/
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_fluxcorr") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015_fluxcorr") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016_fluxcorr")) {
return (new MINERvA_CC1pi0_XSec_1Dth_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu_fluxcorr")) {
return (new MINERvA_CC1pi0_XSec_1Dppi0_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi0_antinu")) {
return (new MINERvA_CC1pi0_XSec_1DTpi0_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DQ2_antinu")) {
return (new MINERvA_CC1pi0_XSec_1DQ2_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_antinu")) {
return (new MINERvA_CC1pi0_XSec_1Dthmu_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_antinu")) {
return (new MINERvA_CC1pi0_XSec_1Dpmu_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DEnu_antinu")) {
return (new MINERvA_CC1pi0_XSec_1DEnu_antinu(samplekey));
// MINERvA CC1pi0 nu
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DQ2_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DEnu_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DWexp_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DPPi0Mass_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DPPi0MassDelta_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DCosAdler_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DPhiAdler_nu")) {
return (new MINERvA_CC1pi0_XSec_1D_nu(samplekey));
/*
CCINC
*/
} else if (!name.compare("MINERvA_CCinc_XSec_2DEavq3_nu")) {
return (new MINERvA_CCinc_XSec_2DEavq3_nu(samplekey));
} else if (!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_C12_CH") ||
!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Pb208_CH")) {
return (new MINERvA_CCinc_XSec_1Dx_ratio(samplekey));
} else if (!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_C12_CH") ||
!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Pb208_CH")) {
return (new MINERvA_CCinc_XSec_1DEnu_ratio(samplekey));
/*
CCDIS
*/
} else if (!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_C12_CH") ||
!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Pb208_CH")) {
return (new MINERvA_CCDIS_XSec_1Dx_ratio(samplekey));
} else if (!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_C12_CH") ||
!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Pb208_CH")) {
return (new MINERvA_CCDIS_XSec_1DEnu_ratio(samplekey));
/*
CC-COH
*/
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_nu")) {
return (new MINERvA_CCCOHPI_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_nu")) {
return (new MINERvA_CCCOHPI_XSec_1DEpi_nu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_nu")) {
return (new MINERvA_CCCOHPI_XSec_1Dth_nu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_nu")) {
return (new MINERvA_CCCOHPI_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_antinu")) {
return (new MINERvA_CCCOHPI_XSec_1DEnu_antinu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_antinu")) {
return (new MINERvA_CCCOHPI_XSec_1DEpi_antinu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_antinu")) {
return (new MINERvA_CCCOHPI_XSec_1Dth_antinu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_antinu")) {
return (new MINERvA_CCCOHPI_XSec_1DQ2_antinu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_joint")) {
return (new MINERvA_CCCOHPI_XSec_joint(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_joint")) {
return (new MINERvA_CCCOHPI_XSec_joint(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_joint")) {
return (new MINERvA_CCCOHPI_XSec_joint(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_joint")) {
return (new MINERvA_CCCOHPI_XSec_joint(samplekey));
/*
T2K Samples
*/
} else
#endif
#ifndef __NO_T2K__
if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu") ||
!name.compare("T2K_CC0pi_XSec_2DPcos_nu_I") ||
!name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) {
return (new T2K_CC0pi_XSec_2DPcos_nu(samplekey));
} else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_nonuniform")) {
return (new T2K_CC0pi_XSec_2DPcos_nu_nonuniform(samplekey));
/*
T2K CC1pi+ CH samples
*/
// Comment these out for now because we don't have the proper data
-
+
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dpmu_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dppi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dppi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DQ2_nu")) {
return (new T2K_CC1pip_CH_XSec_1DQ2_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dq3_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dq3_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthpi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthq3pi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthq3pi_nu(file, rw, type, fkdt));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DWrec_nu")) {
return (new T2K_CC1pip_CH_XSec_1DWrec_nu(file, rw, type, fkdt));
/*
T2K CC1pi+ H2O samples
*/
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuDelta_nu")) {
return (new T2K_CC1pip_H2O_XSec_1DEnuDelta_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuMB_nu")) {
return (new T2K_CC1pip_H2O_XSec_1DEnuMB_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmu_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dcosmu_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmupi_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dcosmupi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcospi_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dcospi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dpmu_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dppi_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dppi_nu(samplekey));
/*
T2K CC0pi + np CH samples
*/
} else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddpt_nu")) {
return (new T2K_CC0pinp_STV_XSec_1Ddpt_nu(samplekey));
// SciBooNE COH studies
} else
#endif
#ifndef __NO_SciBooNE__
if (!name.compare("SciBooNE_CCCOH_STOP_NTrks_nu")) {
return (new SciBooNE_CCCOH_STOP_NTrks_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_1TRK_1DQ2_nu")) {
return (new SciBooNE_CCCOH_1TRK_1DQ2_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPr_1DQ2_nu")) {
return (new SciBooNE_CCCOH_MuPr_1DQ2_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1DQ2_nu")) {
return (new SciBooNE_CCCOH_MuPiVA_1DQ2_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_STOPFINAL_1DQ2_nu")) {
return (new SciBooNE_CCCOH_STOPFINAL_1DQ2_nu(samplekey));
/*
K2K Samples
*/
/*
NC1pi0
*/
} else
#endif
#ifndef __NO_K2K__
if (!name.compare("K2K_NC1pi0_Evt_1Dppi0_nu")) {
return (new K2K_NC1pi0_Evt_1Dppi0_nu(samplekey));
/*
Fake Studies
*/
} else
#endif
if (name.find("ExpMultDist_CCQE_XSec_1D") != std::string::npos &&
name.find("_FakeStudy") != std::string::npos) {
return (
new ExpMultDist_CCQE_XSec_1DVar_FakeStudy(name, file, rw, type, fkdt));
} else if (name.find("ExpMultDist_CCQE_XSec_2D") != std::string::npos &&
name.find("_FakeStudy") != std::string::npos) {
return (
new ExpMultDist_CCQE_XSec_2DVar_FakeStudy(name, file, rw, type, fkdt));
} else if (name.find("GenericFlux_") != std::string::npos) {
return (new GenericFlux_Tester(name, file, rw, type, fkdt));
} else if (name.find("GenericVectors_") != std::string::npos) {
return (new GenericFlux_Vectors(name, file, rw, type, fkdt));
} else if (!name.compare("T2K2017_FakeData")) {
return (new T2K2017_FakeData(samplekey));
} else if (!name.compare("MCStudy_CCQE")) {
return (new MCStudy_CCQEHistograms(name, file, rw, type, fkdt));
} else if (!name.compare("ElectronFlux_FlatTree")) {
return (new ElectronFlux_FlatTree(name, file, rw, type, fkdt));
} else if (name.find("ElectronData_") != std::string::npos) {
return new ElectronScattering_DurhamData(samplekey);
} else if (name.find("MuonValidation_") != std::string::npos) {
return (new MCStudy_MuonValidation(name, file, rw, type, fkdt));
} else if (!name.compare("NIWGOfficialPlots")) {
return (new OfficialNIWGPlots(samplekey));
} else if (!name.compare("Simple_Osc")) {
return (new Simple_Osc(samplekey));
} else if (!name.compare("Smear_SVDUnfold_Propagation_Osc")) {
return (new Smear_SVDUnfold_Propagation_Osc(samplekey));
} else {
ERR(FTL) << "Error: No such sample: " << name << std::endl;
exit(-1);
return NULL;
}
// Return NULL if no sample loaded.
return NULL;
}
}
diff --git a/src/MCStudies/Smearceptance_Tester.cxx b/src/MCStudies/Smearceptance_Tester.cxx
index 072f2bd..e19e54b 100644
--- a/src/MCStudies/Smearceptance_Tester.cxx
+++ b/src/MCStudies/Smearceptance_Tester.cxx
@@ -1,884 +1,881 @@
// 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 "Smearceptance_Tester.h"
#include "SmearceptanceUtils.h"
#include "Smearcepterton.h"
//#define DEBUG_SMEARTESTER
//********************************************************************
/// @brief Class to perform smearceptance MC Studies on a custom measurement
Smearceptance_Tester::Smearceptance_Tester(nuiskey samplekey) {
//********************************************************************
samplekey.Print();
// Sample overview ---------------------------------------------------
std::string descrip =
"Simple measurement class for producing an event summary tree of smeared "
"events.\n";
if (Config::HasPar("NPOT")) {
samplekey.SetS("NPOT", Config::GetParS("NPOT"));
}
if (Config::HasPar("FluxIntegralOverride")) {
samplekey.SetS("FluxIntegralOverride",
Config::GetParS("FluxIntegralOverride"));
}
if (Config::HasPar("TargetVolume")) {
samplekey.SetS("TargetVolume", Config::GetParS("TargetVolume"));
}
if (Config::HasPar("TargetMaterialDensity")) {
samplekey.SetS("TargetMaterialDensity",
Config::GetParS("TargetMaterialDensity"));
}
OutputSummaryTree = true;
if (Config::HasPar("smear.OutputSummaryTree")) {
OutputSummaryTree = Config::GetParI("smear.OutputSummaryTree");
}
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetTitle("Smearceptance Studies");
fSettings.SetDescription(descrip);
fSettings.SetXTitle("XXX");
fSettings.SetYTitle("Number of events");
fSettings.SetEnuRange(0.0, 1E5);
fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
fSettings.DefineAllowedTargets("*");
fSettings.DefineAllowedSpecies("*");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor =
(GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) /
TotalIntegratedFlux();
// Measurement Details
std::vector splitName = GeneralUtils::ParseToStr(fName, "_");
size_t firstUS = fName.find_first_of("_");
std::string smearceptorName = samplekey.GetS("smearceptor");
QLOG(SAM, "Using smearceptor: " << smearceptorName
<< " (parsed from: " << fName << ").");
fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
SetupDefaultHist();
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
covar = StatUtils::GetInvert(fFullCovar);
eventVariables = NULL;
QLOG(SAM, "Smearceptance Flux Scaling Factor = " << fScaleFactor);
if (fScaleFactor <= 0.0) {
ERROR(WRN, "SCALE FACTOR TOO LOW ");
sleep(20);
}
// Setup our TTrees
AddEventVariablesToTree();
smearceptor = &Smearcepterton::Get().GetSmearcepter(smearceptorName);
Int_t RecNBins = 20, TrueNBins = 20;
double RecBinL = 0xdeadbeef, TrueBinL = 0, RecBinH = 10, TrueBinH = 10;
if (Config::HasPar("smear.reconstructed.binning")) {
std::vector args = GeneralUtils::ParseToStr(
Config::GetParS("smear.reconstructed.binning"), ",");
RecNBins = GeneralUtils::StrToInt(args[0]);
RecBinL = GeneralUtils::StrToDbl(args[1]);
RecBinH = GeneralUtils::StrToDbl(args[2]);
TrueNBins = RecNBins;
TrueBinL = RecBinL;
TrueBinH = RecBinH;
}
if (Config::HasPar("smear.true.binning")) {
std::vector args =
GeneralUtils::ParseToStr(Config::GetParS("smear.true.binning"), ",");
TrueNBins = GeneralUtils::StrToInt(args[0]);
TrueBinL = GeneralUtils::StrToDbl(args[1]);
TrueBinH = GeneralUtils::StrToDbl(args[2]);
}
SVDTruncation = 0;
if (Config::HasPar("smear.true.binning")) {
SVDTruncation = Config::GetParI("smear.SVD.truncation");
QLOG(SAM, "Applying SVD truncation of: " << SVDTruncation)
}
ETrueDistrib = NULL;
ETrueDistrib_noweight = NULL;
ERecDistrib = NULL;
RecoSmear = NULL;
if (RecBinL != 0xdeadbeef) {
QLOG(SAM, "Using binning True: " << TrueNBins << ", [" << TrueBinL << " -- "
<< TrueBinH << "], Rec: " << RecNBins
<< ", [" << RecBinL << " -- " << RecBinH
<< "]");
ETrueDistrib = new TH1D("ELep_rate", ";True E_{#nu};Count", TrueNBins,
TrueBinL, TrueBinH);
ETrueDistrib_noweight =
new TH1D("ELep_rate_noweight", ";True E_{#nu};Count", TrueNBins,
TrueBinL, TrueBinH);
ERecDistrib = new TH1D("ELepRec_rate", ";Rec E_{#nu};Count", RecNBins,
RecBinL, RecBinH);
ETrueDistrib->Sumw2();
ERecDistrib->Sumw2();
RecoSmear =
new TH2D("ELepHadVis_Recon", ";True E_{#nu};Recon. E_{#nu}", RecNBins,
RecBinL, RecBinH, TrueNBins, TrueBinL, TrueBinH);
RecoSmear->Sumw2();
}
// Final setup ---------------------------------------------------
FinaliseMeasurement();
}
void Smearceptance_Tester::AddEventVariablesToTree() {
if (OutputSummaryTree) {
// Setup the TTree to save everything
if (!eventVariables) {
Config::Get().out->cd();
eventVariables =
new TTree((fName + "_VARS").c_str(), (fName + "_VARS").c_str());
}
LOG(SAM) << "Adding Event Variables" << std::endl;
eventVariables->Branch("Omega_true", &Omega_true, "Omega_true/F");
eventVariables->Branch("Q2_true", &Q2_true, "Q2_true/F");
eventVariables->Branch("Mode_true", &Mode_true, "Mode_true/I");
eventVariables->Branch("EISLep_true", &EISLep_true, "EISLep_true/F");
eventVariables->Branch("HMFS_clep_true", &HMFS_clep_true);
eventVariables->Branch("HMFS_pip_true", &HMFS_pip_true);
eventVariables->Branch("HMFS_pim_true", &HMFS_pim_true);
eventVariables->Branch("HMFS_cpi_true", &HMFS_cpi_true);
eventVariables->Branch("HMFS_pi0_true", &HMFS_pi0_true);
eventVariables->Branch("HMFS_cK_true", &HMFS_cK_true);
eventVariables->Branch("HMFS_K0_true", &HMFS_K0_true);
eventVariables->Branch("HMFS_p_true", &HMFS_p_true);
eventVariables->Branch("KEFSHad_cpip_true", &KEFSHad_cpip_true,
"KEFSHad_cpip_true/F");
eventVariables->Branch("KEFSHad_cpim_true", &KEFSHad_cpim_true,
"KEFSHad_cpim_true/F");
eventVariables->Branch("KEFSHad_cpi_true", &KEFSHad_cpi_true,
"KEFSHad_cpi_true/F");
eventVariables->Branch("TEFSHad_pi0_true", &TEFSHad_pi0_true,
"TEFSHad_pi0_true/F");
eventVariables->Branch("KEFSHad_cK_true", &KEFSHad_cK_true,
"KEFSHad_cK_true/F");
eventVariables->Branch("KEFSHad_K0_true", &KEFSHad_K0_true,
"KEFSHad_K0_true/F");
eventVariables->Branch("KEFSHad_p_true", &KEFSHad_p_true,
"KEFSHad_p_true/F");
eventVariables->Branch("KEFSHad_n_true", &KEFSHad_n_true,
"KEFSHad_n_true/F");
eventVariables->Branch("EFSHad_true", &EFSHad_true, "EFSHad_true/F");
eventVariables->Branch("EFSChargedEMHad_true", &EFSChargedEMHad_true,
"EFSChargedEMHad_true/F");
eventVariables->Branch("EFSLep_true", &EFSLep_true, "EFSLep_true/F");
eventVariables->Branch("EFSgamma_true", &EFSgamma_true, "EFSgamma_true/F");
eventVariables->Branch("PDGISLep_true", &PDGISLep_true, "PDGISLep_true/I");
eventVariables->Branch("PDGFSLep_true", &PDGFSLep_true, "PDGFSLep_true/I");
eventVariables->Branch("Nprotons_true", &Nprotons_true, "Nprotons_true/I");
eventVariables->Branch("Nneutrons_true", &Nneutrons_true,
"Nneutrons_true/I");
eventVariables->Branch("Ncpiplus_true", &Ncpiplus_true, "Ncpiplus_true/I");
eventVariables->Branch("Ncpiminus_true", &Ncpiminus_true,
"Ncpiminus_true/I");
eventVariables->Branch("Ncpi_true", &Ncpi_true, "Ncpi_true/I");
eventVariables->Branch("Npi0_true", &Npi0_true, "Npi0_true/I");
eventVariables->Branch("NcK_true", &NcK_true, "NcK_true/I");
eventVariables->Branch("NK0_true", &NK0_true, "NK0_true/I");
eventVariables->Branch("HMFS_clep_rec", &HMFS_clep_rec);
eventVariables->Branch("HMFS_pip_rec", &HMFS_pip_rec);
eventVariables->Branch("HMFS_pim_rec", &HMFS_pim_rec);
eventVariables->Branch("HMFS_cpi_rec", &HMFS_cpi_rec);
eventVariables->Branch("HMFS_pi0_rec", &HMFS_pi0_rec);
eventVariables->Branch("HMFS_cK_rec", &HMFS_cK_rec);
eventVariables->Branch("HMFS_K0_rec", &HMFS_K0_rec);
eventVariables->Branch("HMFS_p_rec", &HMFS_p_rec);
eventVariables->Branch("KEFSHad_cpip_rec", &KEFSHad_cpip_rec,
"KEFSHad_cpip_rec/F");
eventVariables->Branch("KEFSHad_cpim_rec", &KEFSHad_cpim_rec,
"KEFSHad_cpim_rec/F");
eventVariables->Branch("KEFSHad_cpi_rec", &KEFSHad_cpi_rec,
"KEFSHad_cpi_rec/F");
eventVariables->Branch("TEFSHad_pi0_rec", &TEFSHad_pi0_rec,
"TEFSHad_pi0_rec/F");
eventVariables->Branch("KEFSHad_cK_rec", &KEFSHad_cK_rec,
"KEFSHad_cK_rec/F");
eventVariables->Branch("KEFSHad_K0_rec", &KEFSHad_K0_rec,
"KEFSHad_K0_rec/F");
eventVariables->Branch("KEFSHad_p_rec", &KEFSHad_p_rec, "KEFSHad_p_rec/F");
eventVariables->Branch("KEFSHad_n_rec", &KEFSHad_n_rec, "KEFSHad_n_rec/F");
eventVariables->Branch("EFSHad_rec", &EFSHad_rec, "EFSHad_rec/F");
eventVariables->Branch("EFSLep_rec", &EFSLep_rec, "EFSLep_rec/F");
eventVariables->Branch("EFSVis_cpip", &EFSVis_cpip, "EFSVis_cpip/F");
eventVariables->Branch("EFSVis_cpim", &EFSVis_cpim, "EFSVis_cpim/F");
eventVariables->Branch("EFSVis_cpi", &EFSVis_cpi, "EFSVis_cpi/F");
eventVariables->Branch("EFSVis_pi0", &EFSVis_pi0, "EFSVis_pi0/F");
eventVariables->Branch("EFSVis_cK", &EFSVis_cK, "EFSVis_cK/F");
eventVariables->Branch("EFSVis_K0", &EFSVis_K0, "EFSVis_K0/F");
eventVariables->Branch("EFSVis_p", &EFSVis_p, "EFSVis_p/F");
eventVariables->Branch("EFSVis_n", &EFSVis_n, "EFSVis_n/F");
eventVariables->Branch("EFSVis_gamma", &EFSVis_gamma, "EFSVis_gamma/F");
eventVariables->Branch("EFSVis_other", &EFSVis_other, "EFSVis_other/F");
eventVariables->Branch("EFSVis", &EFSVis, "EFSVis/F");
eventVariables->Branch("FSCLep_seen", &FSCLep_seen, "FSCLep_seen/I");
eventVariables->Branch("Nprotons_seen", &Nprotons_seen, "Nprotons_seen/I");
eventVariables->Branch("Nneutrons_seen", &Nneutrons_seen,
"Nneutrons_seen/I");
eventVariables->Branch("Ncpip_seen", &Ncpip_seen, "Ncpip_seen/I");
eventVariables->Branch("Ncpim_seen", &Ncpim_seen, "Ncpim_seen/I");
eventVariables->Branch("Ncpi_seen", &Ncpi_seen, "Ncpi_seen/I");
eventVariables->Branch("Npi0_seen", &Npi0_seen, "Npi0_seen/I");
eventVariables->Branch("NcK_seen", &NcK_seen, "NcK_seen/I");
eventVariables->Branch("NK0_seen", &NK0_seen, "NK0_seen/I");
eventVariables->Branch("Nothers_seen", &Nothers_seen, "Nothers_seen/I");
eventVariables->Branch("EISLep_QE_rec", &EISLep_QE_rec, "EISLep_QE_rec/F");
eventVariables->Branch("EISLep_LepHad_rec", &EISLep_LepHad_rec,
"EISLep_LepHad_rec/F");
eventVariables->Branch("EISLep_LepHadVis_rec", &EISLep_LepHadVis_rec,
"EISLep_LepHadVis_rec/F");
eventVariables->Branch("Nprotons_contributed", &Nprotons_contributed,
"Nprotons_contributed/I");
eventVariables->Branch("Nneutrons_contributed", &Nneutrons_contributed,
"Nneutrons_contributed/I");
eventVariables->Branch("Ncpip_contributed", &Ncpip_contributed,
"Ncpip_contributed/I");
eventVariables->Branch("Ncpim_contributed", &Ncpim_contributed,
"Ncpim_contributed/I");
eventVariables->Branch("Ncpi_contributed", &Ncpi_contributed,
"Ncpi_contributed/I");
eventVariables->Branch("Npi0_contributed", &Npi0_contributed,
"Npi0_contributed/I");
eventVariables->Branch("NcK_contributed", &NcK_contributed,
"NcK_contributed/I");
eventVariables->Branch("NK0_contributed", &NK0_contributed,
"NK0_contributed/I");
eventVariables->Branch("Ngamma_contributed", &Ngamma_contributed,
"Ngamma_contributed/I");
eventVariables->Branch("Nothers_contibuted", &Nothers_contibuted,
"Nothers_contibuted/I");
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("FluxWeight", &FluxWeight, "FluxWeight/F");
eventVariables->Branch("EffWeight", &EffWeight, "EffWeight/F");
xsecScaling = fScaleFactor;
eventVariables->Branch("xsecScaling", &xsecScaling, "xsecScaling/F");
eventVariables->Branch("flagCCINC_true", &flagCCINC_true,
"flagCCINC_true/O");
- eventVariables->Branch("flagCC0K_true", &flagCC0K_true,
- "flagCC0K_true/O");
+ eventVariables->Branch("flagCC0K_true", &flagCC0K_true, "flagCC0K_true/O");
eventVariables->Branch("flagCC0Pi_true", &flagCC0Pi_true,
"flagCC0Pi_true/O");
eventVariables->Branch("flagCC1Pi_true", &flagCC1Pi_true,
"flagCC1Pi_true/O");
eventVariables->Branch("flagCCINC_rec", &flagCCINC_rec, "flagCCINC_rec/O");
eventVariables->Branch("flagCC0K_rec", &flagCC0K_rec, "flagCC0K_rec/O");
eventVariables->Branch("flagCC0Pi_rec", &flagCC0Pi_rec, "flagCC0Pi_rec/O");
eventVariables->Branch("flagCC1Pi_rec", &flagCC1Pi_rec, "flagCC1Pi_rec/O");
}
PredEvtRateWeight = 1;
if (fEvtRateScaleFactor != 0xdeadbeef) {
if (OutputSummaryTree) {
eventVariables->Branch("PredEvtRateWeight", &PredEvtRateWeight,
"PredEvtRateWeight/F");
}
PredEvtRateWeight = fScaleFactor * fEvtRateScaleFactor;
}
}
template
int CountNPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum +=
std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it]);
}
return sum;
}
template
int CountNNotPdgsSeen(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
- for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
- sum +=
- (std::count(ri.RecObjClass.begin(), ri.RecObjClass.end(), pdgs[pdg_it])
- ? 0
- : 1);
+ for (size_t p_it = 0; p_it < ri.RecObjClass.size(); ++p_it) {
+ if (!std::count(pdgs, pdgs + N, ri.RecObjClass[p_it])) {
+ sum++;
+ }
}
return sum;
}
template
int CountNPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
sum += std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
pdgs[pdg_it]);
}
return sum;
}
template
int CountNNotPdgsContributed(RecoInfo ri, int const (&pdgs)[N]) {
int sum = 0;
- for (size_t pdg_it = 0; pdg_it < N; ++pdg_it) {
- sum += (std::count(ri.TrueContribPDGs.begin(), ri.TrueContribPDGs.end(),
- pdgs[pdg_it])
- ? 0
- : 1);
+ for (size_t p_it = 0; p_it < ri.TrueContribPDGs.size(); ++p_it) {
+ if (!std::count(pdgs, pdgs + N, ri.TrueContribPDGs[p_it])) {
+ sum++;
+ }
}
return sum;
}
TLorentzVector GetHMFSRecParticles(RecoInfo ri, int pdg) {
TLorentzVector mom(0, 0, 0, 0);
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if ((ri.RecObjClass[p_it] == pdg) &&
(mom.Mag() < ri.RecObjMom[p_it].Mag())) {
mom.SetXYZM(ri.RecObjMom[p_it].X(), ri.RecObjMom[p_it].Y(),
ri.RecObjMom[p_it].Z(),
PhysConst::GetMass(ri.RecObjClass[p_it]) * 1.0E3);
}
}
return mom;
}
template
double SumKE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.RecObjClass[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass) - mass;
}
return sum;
}
template
double SumTE_RecoInfo(RecoInfo ri, int const (&pdgs)[N], double mass) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecObjMom.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.RecObjClass[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += sqrt(ri.RecObjMom[p_it].Mag2() + mass * mass);
}
return sum;
}
template
double SumVisE_RecoInfo(RecoInfo ri, int const (&pdgs)[N]) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
if (!std::count(pdgs, pdgs + N,
ri.TrueContribPDGs[p_it])) { // If we don't care about this
// particle type.
continue;
}
sum += ri.RecVisibleEnergy[p_it];
}
return sum;
}
template
double SumVisE_RecoInfo_NotPdgs(RecoInfo ri, int const (&pdgs)[N]) {
double sum = 0;
for (size_t p_it = 0; p_it < ri.RecVisibleEnergy.size(); ++p_it) {
if (std::count(pdgs, pdgs + N,
ri.TrueContribPDGs[p_it])) { // If we know about this
// particle type.
continue;
}
sum += ri.RecVisibleEnergy[p_it];
}
return sum;
}
//********************************************************************
void Smearceptance_Tester::FillEventVariables(FitEvent *event) {
//********************************************************************
static int const cpipPDG[] = {211};
static int const cpimPDG[] = {-211};
static int const pi0PDG[] = {111};
static int const cKPDG[] = {321, -321};
static int const K0PDG[] = {311, 310, 130};
static int const ProtonPDG[] = {2212};
static int const NeutronPDG[] = {2112};
static int const GammaPDG[] = {22};
static int const CLeptonPDGs[] = {11, 13, 15, -11, -13, -15};
static int const ExplicitPDGs[] = {211, -211, 111, 321, -321, 311, 310, 130,
2212, 2112, 22, 11, 13, 15, 12, 14,
16, -11, -13, -15, -12, -14, -16};
RecoInfo *ri = smearceptor->Smearcept(event);
//** START Pions
HMFS_clep_true = TLorentzVector(0, 0, 0, 0);
HMFS_clep_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsCLep = event->GetHMFSParticle(CLeptonPDGs);
if (fsCLep) {
HMFS_clep_true = fsCLep->P4();
HMFS_clep_rec = GetHMFSRecParticles(*ri, fsCLep->PDG());
}
//** END Charged leptons
//** START Pions
HMFS_pip_true = TLorentzVector(0, 0, 0, 0);
HMFS_pip_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPip = event->GetHMFSPiPlus();
if (fsPip) {
HMFS_pip_true = fsPip->P4();
HMFS_pip_rec = GetHMFSRecParticles(*ri, fsPip->PDG());
}
HMFS_pim_true = TLorentzVector(0, 0, 0, 0);
HMFS_pim_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPim = event->GetHMFSPiMinus();
if (fsPim) {
HMFS_pim_true = fsPim->P4();
HMFS_pim_rec = GetHMFSRecParticles(*ri, fsPim->PDG());
}
HMFS_cpi_true = TLorentzVector(0, 0, 0, 0);
HMFS_cpi_rec = TLorentzVector(0, 0, 0, 0);
if (fsPip || fsPim) {
if (!fsPip) {
HMFS_cpi_true = HMFS_pim_true;
HMFS_cpi_rec = HMFS_pim_rec;
} else if (!fsPim) {
HMFS_cpi_true = HMFS_pip_true;
HMFS_cpi_rec = HMFS_pip_rec;
} else {
HMFS_cpi_true =
(fsPip->p2() > fsPim->p2()) ? HMFS_pip_true : HMFS_pim_true;
HMFS_cpi_rec = (fsPip->p2() > fsPim->p2()) ? HMFS_pip_rec : HMFS_pim_rec;
}
}
HMFS_pi0_true = TLorentzVector(0, 0, 0, 0);
HMFS_pi0_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsPi0 = event->GetHMFSPiZero();
if (fsPi0) {
HMFS_pi0_true = fsPi0->P4();
HMFS_pi0_rec = GetHMFSRecParticles(*ri, fsPi0->PDG());
}
//** END Pions
//** START Kaons
HMFS_cK_true = TLorentzVector(0, 0, 0, 0);
HMFS_cK_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fscK = event->GetHMFSParticle(cKPDG);
if (fscK) {
HMFS_cK_true = fscK->P4();
HMFS_cK_rec = GetHMFSRecParticles(*ri, fscK->PDG());
}
HMFS_K0_true = TLorentzVector(0, 0, 0, 0);
HMFS_K0_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsK0 = event->GetHMFSParticle(K0PDG);
if (fsK0) {
HMFS_K0_true = fsK0->P4();
HMFS_K0_rec = GetHMFSRecParticles(*ri, fsK0->PDG());
}
//** END Kaons
//** START Nucleons
HMFS_p_true = TLorentzVector(0, 0, 0, 0);
HMFS_p_rec = TLorentzVector(0, 0, 0, 0);
FitParticle *fsP = event->GetHMFSProton();
if (fsP) {
HMFS_p_true = fsP->P4();
HMFS_p_rec = GetHMFSRecParticles(*ri, fsP->PDG());
}
TLorentzVector FourMomentumTransfer =
(event->GetHMISAnyLeptons()->P4() - event->GetHMFSAnyLeptons()->P4());
Omega_true = FourMomentumTransfer.E();
Q2_true = -1 * FourMomentumTransfer.Mag2();
Mode_true = event->Mode;
EISLep_true = event->GetHMISAnyLeptons()->E();
KEFSHad_cpip_true = FitUtils::SumTE_PartVect(event->GetAllFSPiPlus());
KEFSHad_cpim_true = FitUtils::SumTE_PartVect(event->GetAllFSPiMinus());
KEFSHad_cpi_true = KEFSHad_cpip_true + KEFSHad_cpim_true;
TEFSHad_pi0_true = FitUtils::SumTE_PartVect(event->GetAllFSPiZero());
KEFSHad_cK_true = FitUtils::SumTE_PartVect(event->GetAllFSParticle(cKPDG));
KEFSHad_K0_true = FitUtils::SumTE_PartVect(event->GetAllFSParticle(K0PDG));
KEFSHad_p_true = FitUtils::SumKE_PartVect(event->GetAllFSProton());
KEFSHad_n_true = FitUtils::SumKE_PartVect(event->GetAllFSNeutron());
EFSHad_true =
KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true + KEFSHad_n_true;
EFSChargedEMHad_true = KEFSHad_cpi_true + TEFSHad_pi0_true + KEFSHad_p_true +
KEFSHad_cK_true + KEFSHad_K0_true;
EFSLep_true = event->GetHMFSAnyLeptons()->E();
EFSgamma_true = FitUtils::SumTE_PartVect(event->GetAllFSPhoton());
PDGISLep_true = event->GetHMISAnyLeptons()->PDG();
PDGFSLep_true = event->GetHMFSAnyLeptons()->PDG();
Nprotons_true = event->GetAllFSProton().size();
Nneutrons_true = event->GetAllFSNeutron().size();
Ncpiplus_true = event->GetAllFSPiPlus().size();
Ncpiminus_true = event->GetAllFSPiMinus().size();
Ncpi_true = Ncpiplus_true + Ncpiminus_true;
Npi0_true = event->GetAllFSPiZero().size();
NcK_true = event->GetAllFSParticle(cKPDG).size();
NK0_true = event->GetAllFSParticle(K0PDG).size();
KEFSHad_cpip_rec =
SumKE_RecoInfo(*ri, cpipPDG, PhysConst::mass_cpi * PhysConst::mass_MeV);
KEFSHad_cpim_rec =
SumKE_RecoInfo(*ri, cpimPDG, PhysConst::mass_cpi * PhysConst::mass_MeV);
KEFSHad_cpi_rec = KEFSHad_cpip_rec + KEFSHad_cpim_rec;
TEFSHad_pi0_rec =
SumTE_RecoInfo(*ri, pi0PDG, PhysConst::mass_pi0 * PhysConst::mass_MeV);
KEFSHad_cK_rec =
SumKE_RecoInfo(*ri, cKPDG, PhysConst::mass_cK * PhysConst::mass_MeV);
KEFSHad_K0_rec =
SumKE_RecoInfo(*ri, K0PDG, PhysConst::mass_K0 * PhysConst::mass_MeV);
KEFSHad_p_rec = SumKE_RecoInfo(*ri, ProtonPDG,
PhysConst::mass_proton * PhysConst::mass_MeV);
KEFSHad_n_rec = SumKE_RecoInfo(*ri, NeutronPDG,
PhysConst::mass_neutron * PhysConst::mass_MeV);
EFSHad_rec =
KEFSHad_cpi_rec + TEFSHad_pi0_rec + KEFSHad_p_rec + KEFSHad_n_rec;
TLorentzVector FSLepMom_rec(0, 0, 0, 0);
if (event->GetHMFSAnyLeptons()) {
FSLepMom_rec = GetHMFSRecParticles(*ri, event->GetHMFSAnyLeptons()->PDG());
EFSLep_rec = FSLepMom_rec.E();
} else {
EFSLep_rec = 0;
}
EFSVis_cpip = SumVisE_RecoInfo(*ri, cpipPDG);
EFSVis_cpim = SumVisE_RecoInfo(*ri, cpimPDG);
EFSVis_cpi = EFSVis_cpip + EFSVis_cpim;
EFSVis_pi0 = SumVisE_RecoInfo(*ri, pi0PDG);
EFSVis_cK = SumVisE_RecoInfo(*ri, cKPDG);
EFSVis_K0 = SumVisE_RecoInfo(*ri, K0PDG);
EFSVis_p = SumVisE_RecoInfo(*ri, ProtonPDG);
EFSVis_n = SumVisE_RecoInfo(*ri, NeutronPDG);
EFSVis_gamma = SumVisE_RecoInfo(*ri, GammaPDG);
EFSVis_other = SumVisE_RecoInfo_NotPdgs(*ri, ExplicitPDGs);
EFSVis = EFSVis_cpi + EFSVis_pi0 + EFSVis_p + EFSVis_n + EFSVis_gamma +
EFSVis_cK + EFSVis_K0;
FSCLep_seen = CountNPdgsSeen(*ri, CLeptonPDGs);
Nprotons_seen = CountNPdgsSeen(*ri, ProtonPDG);
Nneutrons_seen = CountNPdgsSeen(*ri, NeutronPDG);
Ncpip_seen = CountNPdgsSeen(*ri, cpipPDG);
Ncpim_seen = CountNPdgsSeen(*ri, cpimPDG);
Ncpi_seen = Ncpip_seen + Ncpim_seen;
Npi0_seen = CountNPdgsSeen(*ri, pi0PDG);
NcK_seen = CountNPdgsSeen(*ri, cKPDG);
NK0_seen = CountNPdgsSeen(*ri, K0PDG);
Nothers_seen = CountNNotPdgsSeen(*ri, ExplicitPDGs);
if (FSCLep_seen && (FSLepMom_rec.Mag() > 1E-8)) {
EISLep_QE_rec =
FitUtils::EnuQErec(FSLepMom_rec.Mag() / 1000.0, FSLepMom_rec.CosTheta(),
34, PDGFSLep_true > 0) *
1000.0;
} else {
EISLep_QE_rec = 0;
}
EISLep_LepHad_rec = EFSLep_rec + EFSHad_rec;
EISLep_LepHadVis_rec = EFSLep_rec + EFSHad_rec + EFSVis;
Nprotons_contributed = CountNPdgsContributed(*ri, ProtonPDG);
Nneutrons_contributed = CountNPdgsContributed(*ri, NeutronPDG);
Ncpip_contributed = CountNPdgsContributed(*ri, cpipPDG);
Ncpim_contributed = CountNPdgsContributed(*ri, cpimPDG);
Ncpi_contributed = Ncpip_contributed + Ncpim_contributed;
Npi0_contributed = CountNPdgsContributed(*ri, pi0PDG);
NcK_contributed = CountNPdgsContributed(*ri, cKPDG);
NK0_contributed = CountNPdgsContributed(*ri, K0PDG);
Ngamma_contributed = CountNPdgsContributed(*ri, GammaPDG);
Nothers_contibuted = CountNNotPdgsContributed(*ri, ExplicitPDGs);
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
FluxWeight = GetFluxHistogram()->GetBinContent(
GetFluxHistogram()->FindBin(EISLep_true)) /
GetFluxHistogram()->Integral();
EffWeight = ri->Weight;
flagCCINC_true = PDGFSLep_true & 1;
flagCC0K_true = (NcK_true + NK0_true) == 0;
flagCC0Pi_true =
flagCCINC_true && flagCC0K_true && ((Ncpi_true + Npi0_true) == 0);
flagCC1Pi_true =
flagCCINC_true && flagCC0K_true && ((Ncpi_true + Npi0_true) == 1);
flagCCINC_rec = FSCLep_seen && flagCCINC_true;
flagCC0K_rec = (NcK_seen + NK0_seen) == 0;
flagCC0Pi_rec =
flagCCINC_rec && flagCC0K_rec && ((Ncpi_seen + Npi0_seen) == 0);
flagCC1Pi_rec =
flagCCINC_rec && flagCC0K_rec && ((Ncpi_seen + Npi0_seen) == 1);
if (OutputSummaryTree) {
// Fill the eventVariables Tree
eventVariables->Fill();
}
if (RecoSmear) {
RecoSmear->Fill(EISLep_true / 1000.0,
flagCCINC_rec ? EISLep_LepHadVis_rec / 1000.0 : -1, Weight);
ETrueDistrib_noweight->Fill(EISLep_true / 1000.0,
flagCCINC_true ? Weight : 0);
ETrueDistrib->Fill(EISLep_true / 1000.0,
flagCCINC_true ? Weight * PredEvtRateWeight : 0);
ERecDistrib->Fill(EISLep_LepHadVis_rec / 1000.0,
flagCCINC_rec ? Weight * PredEvtRateWeight : 0);
}
};
//********************************************************************
void Smearceptance_Tester::Write(std::string drawOpt) {
//********************************************************************
if (OutputSummaryTree) {
// First save the TTree
eventVariables->Write();
}
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
if (!RecoSmear) {
return;
}
TH2D *SmearMatrix_ev =
static_cast(RecoSmear->Clone("ELepHadVis_Smear_ev"));
for (Int_t trueAxis_it = 1;
trueAxis_it < RecoSmear->GetXaxis()->GetNbins() + 1; ++trueAxis_it) {
double NEISLep = ETrueDistrib_noweight->GetBinContent(trueAxis_it);
for (Int_t recoAxis_it = 1;
recoAxis_it < RecoSmear->GetYaxis()->GetNbins() + 1; ++recoAxis_it) {
if (NEISLep > std::numeric_limits::epsilon()) {
SmearMatrix_ev->SetBinContent(
trueAxis_it, recoAxis_it,
SmearMatrix_ev->GetBinContent(trueAxis_it, recoAxis_it) / NEISLep);
}
}
}
ETrueDistrib_noweight->Write();
ETrueDistrib->Write();
ERecDistrib->Write();
RecoSmear->Write();
SmearMatrix_ev->Write();
TH2D *ResponseMatrix_ev =
SmearceptanceUtils::SVDGetInverse(SmearMatrix_ev, SVDTruncation);
ResponseMatrix_ev->SetName("ResponseMatrix_ev");
ResponseMatrix_ev->Write();
#ifdef DEBUG_SMEARTESTER
TMatrixD SmearMatrix_ev_md = SmearceptanceUtils::GetMatrix(SmearMatrix_ev);
TH1D *SmearedEvt = static_cast(ERecDistrib->Clone());
SmearedEvt->SetNameTitle("SmearedEvt", ";Rec E_{#nu}; count");
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
ETrueDistrib, SmearedEvt, SmearMatrix_ev_md, 5000, false);
SmearedEvt->Write();
SmearedEvt->Scale(1, "width");
SmearedEvt->SetName("SmearedEvt_bw");
SmearedEvt->Write();
#endif
#ifdef __PROB3PP_ENABLED__
FitWeight *fw = FitBase::GetRW();
if (fw->HasRWEngine(kOSCILLATION)) {
OscWeightEngine *oscWE =
dynamic_cast(fw->GetRWEngine(kOSCILLATION));
TGraph POsc;
POsc.Set(1E4 - 1);
double min = ETrueDistrib->GetXaxis()->GetBinLowEdge(1);
double step = (ETrueDistrib->GetXaxis()->GetBinUpEdge(
ETrueDistrib->GetXaxis()->GetNbins()) -
ETrueDistrib->GetXaxis()->GetBinLowEdge(1)) /
double(1E4);
for (size_t i = 1; i < 1E4; ++i) {
double enu = min + i * step;
double ow = oscWE->CalcWeight(enu, 14);
if (ow != ow) {
std::cout << "Bad osc weight for ENu: " << enu << std::endl;
}
POsc.SetPoint(i - 1, enu, ow);
}
POsc.Write("POsc", TObject::kOverwrite);
}
#endif
TMatrixD ResponseMatrix_evt_md =
SmearceptanceUtils::GetMatrix(ResponseMatrix_ev);
TH1D *Unfolded_enu_obs = static_cast(ETrueDistrib->Clone());
Unfolded_enu_obs->SetNameTitle("UnfoldedENu_evt", ";True E_{#nu};count");
SmearceptanceUtils::PushTH1ThroughMatrixWithErrors(
ERecDistrib, Unfolded_enu_obs, ResponseMatrix_evt_md, 5000, false);
Unfolded_enu_obs->Write();
Unfolded_enu_obs->Scale(1, "width");
Unfolded_enu_obs->SetName("UnfoldedENu_evt_bw");
Unfolded_enu_obs->Write();
ETrueDistrib->Scale(1, "width");
ETrueDistrib->SetName("ELep_rate_bw");
ETrueDistrib->Write();
ERecDistrib->Scale(1, "width");
ERecDistrib->SetName("ELepRec_rate_bw");
ERecDistrib->Write();
}
// -------------------------------------------------------------------
// Purely MC Plot
// Following functions are just overrides to handle this
// -------------------------------------------------------------------
//********************************************************************
/// Everything is classed as signal...
bool Smearceptance_Tester::isSignal(FitEvent *event) {
//********************************************************************
(void)event;
return true;
};
//********************************************************************
void Smearceptance_Tester::ScaleEvents() {
//********************************************************************
// Saving everything to a TTree so no scaling required
return;
}
//********************************************************************
void Smearceptance_Tester::ApplyNormScale(float norm) {
//********************************************************************
// Saving everything to a TTree so no scaling required
fCurrentNorm = norm;
return;
}
//********************************************************************
void Smearceptance_Tester::FillHistograms() {
//********************************************************************
// No Histograms need filling........
return;
}
//********************************************************************
void Smearceptance_Tester::ResetAll() {
//********************************************************************
if (OutputSummaryTree) {
eventVariables->Reset();
}
return;
}
//********************************************************************
float Smearceptance_Tester::GetChi2() {
//********************************************************************
// No Likelihood to test, purely MC
return 0.0;
}
diff --git a/src/Smearceptance/Smearcepterton.cxx b/src/Smearceptance/Smearcepterton.cxx
index 192a1e9..005644e 100644
--- a/src/Smearceptance/Smearcepterton.cxx
+++ b/src/Smearceptance/Smearcepterton.cxx
@@ -1,345 +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 .
*******************************************************************************/
#include "Smearcepterton.h"
#include "EfficiencyApplicator.h"
#include "GaussianSmearer.h"
#include "MetaSimpleSmearcepter.h"
#include "ThresholdAccepter.h"
#include "TrackedMomentumMatrixSmearer.h"
#include "VisECoalescer.h"
#include
#ifdef __USE_DYNSAMPLES__
#include "TRegexp.h"
#include
// linux
#include
DynamicSmearceptorFactory::DynamicSmearceptorFactory()
: NSmearceptors(0), NManifests(0) {
LoadPlugins();
QLOG(FIT, "Loaded " << NSmearceptors << " from " << NManifests
<< " shared object libraries.");
}
DynamicSmearceptorFactory* DynamicSmearceptorFactory::glblDSF = NULL;
DynamicSmearceptorFactory::PluginManifest::~PluginManifest() {
for (size_t i_it = 0; i_it < Instances.size(); ++i_it) {
(*(DSF_DestroySmearceptor))(Instances[i_it]);
}
}
std::string EnsureTrailingSlash(std::string const& inp) {
if (!inp.length()) {
return "/";
}
if (inp[inp.length() - 1] == '/') {
return inp;
}
return inp + "/";
}
void DynamicSmearceptorFactory::LoadPlugins() {
std::vector SearchDirectories;
if (Config::HasPar("dynamic_smearceptor.path")) {
SearchDirectories = GeneralUtils::ParseToStr(
Config::GetParS("dynamic_smearceptor.path"), ":");
}
char const* envPath = getenv("NUISANCE_DS_PATH");
if (envPath) {
std::vector envPaths = GeneralUtils::ParseToStr(envPath, ":");
for (size_t ep_it = 0; ep_it < envPaths.size(); ++ep_it) {
SearchDirectories.push_back(envPaths[ep_it]);
}
}
if (!SearchDirectories.size()) {
char const* pwdPath = getenv("PWD");
if (pwdPath) {
SearchDirectories.push_back(pwdPath);
}
}
for (size_t sp_it = 0; sp_it < SearchDirectories.size(); ++sp_it) {
std::string dirpath = EnsureTrailingSlash(SearchDirectories[sp_it]);
QLOG(FIT, "Searching for dynamic smearceptor manifests in: " << dirpath);
Ssiz_t len = 0;
DIR* dir;
struct dirent* ent;
dir = opendir(dirpath.c_str());
if (dir != NULL) {
TRegexp matchExp("*.so", true);
while ((ent = readdir(dir)) != NULL) {
if (matchExp.Index(TString(ent->d_name), &len) != Ssiz_t(-1)) {
QLOG(FIT, "\tFound shared object: "
<< ent->d_name << " checking for relevant methods...");
void* dlobj =
dlopen((dirpath + ent->d_name).c_str(), RTLD_NOW | RTLD_GLOBAL);
char const* dlerr_cstr = dlerror();
std::string dlerr;
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN, "\tDL Load Error: " << dlerr);
continue;
}
PluginManifest plgManif;
plgManif.dllib = dlobj;
plgManif.soloc = (dirpath + ent->d_name);
plgManif.DSF_NSmearceptors = reinterpret_cast(
dlsym(dlobj, "DSF_NSmearceptors"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN, "\tFailed to load symbol \"DSF_NSmearceptors\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.DSF_GetSmearceptorName =
reinterpret_cast(
dlsym(dlobj, "DSF_GetSmearceptorName"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN,
"\tFailed to load symbol \"DSF_GetSmearceptorName\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.DSF_GetSmearceptor =
reinterpret_cast(
dlsym(dlobj, "DSF_GetSmearceptor"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN, "\tFailed to load symbol \"DSF_GetSmearceptor\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.DSF_DestroySmearceptor =
reinterpret_cast(
dlsym(dlobj, "DSF_DestroySmearceptor"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
ERROR(WRN, "Failed to load symbol \"DSF_DestroySmearceptor\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.NSmearceptors = (*(plgManif.DSF_NSmearceptors))();
QLOG(FIT, "\tSuccessfully loaded dynamic smearceptor manifest: "
<< plgManif.soloc << ". Contains "
<< plgManif.NSmearceptors << " smearceptors.");
for (size_t smp_it = 0; smp_it < plgManif.NSmearceptors; ++smp_it) {
char const* smp_name = (*(plgManif.DSF_GetSmearceptorName))(smp_it);
if (!smp_name) {
THROW("Could not load smearceptor "
<< smp_it << " / " << plgManif.NSmearceptors << " from "
<< plgManif.soloc);
}
if (Smearceptors.count(smp_name)) {
ERROR(WRN, "Already loaded a smearceptor named: \""
<< smp_name << "\". cannot load duplciates. This "
"smearceptor will be skipped.");
continue;
}
plgManif.SmearceptorsProvided.push_back(smp_name);
Smearceptors[smp_name] = std::make_pair(plgManif.soloc, smp_it);
QLOG(FIT, "\t\t" << smp_name);
}
if (plgManif.SmearceptorsProvided.size()) {
Manifests[plgManif.soloc] = plgManif;
NSmearceptors += plgManif.SmearceptorsProvided.size();
NManifests++;
} else {
dlclose(dlobj);
}
}
}
closedir(dir);
} else {
ERROR(WRN, "Tried to open non-existant directory.");
}
}
}
DynamicSmearceptorFactory& DynamicSmearceptorFactory::Get() {
if (!glblDSF) {
glblDSF = new DynamicSmearceptorFactory();
}
return *glblDSF;
}
void DynamicSmearceptorFactory::Print() {
std::map > ManifestSmearceptors;
for (std::map >::iterator smp_it =
Smearceptors.begin();
smp_it != Smearceptors.end(); ++smp_it) {
if (!ManifestSmearceptors.count(smp_it->second.first)) {
ManifestSmearceptors[smp_it->second.first] = std::vector();
}
ManifestSmearceptors[smp_it->second.first].push_back(smp_it->first);
}
QLOG(FIT, "Dynamic smearceptor manifest: ");
for (std::map >::iterator m_it =
ManifestSmearceptors.begin();
m_it != ManifestSmearceptors.end(); ++m_it) {
QLOG(FIT, "\tLibrary " << m_it->first << " contains: ");
for (size_t s_it = 0; s_it < m_it->second.size(); ++s_it) {
QLOG(FIT, "\t\t" << m_it->second[s_it]);
}
}
}
bool DynamicSmearceptorFactory::HasSmearceptor(std::string const& name) {
return Smearceptors.count(name);
}
bool DynamicSmearceptorFactory::HasSmearceptor(nuiskey& smearceptorkey) {
return HasSmearceptor(smearceptorkey.GetElementName());
}
ISmearcepter* DynamicSmearceptorFactory::CreateSmearceptor(
nuiskey& smearceptorkey) {
if (!HasSmearceptor(smearceptorkey)) {
ERROR(WRN, "Asked to load unknown smearceptor: \""
<< smearceptorkey.GetElementName() << "\".");
return NULL;
}
std::pair smearceptor =
Smearceptors[smearceptorkey.GetElementName()];
QLOG(SAM, "\tLoading smearceptor " << smearceptor.second << " from "
<< smearceptor.first);
ISmearcepter* smear = (*(Manifests[smearceptor.first].DSF_GetSmearceptor))(
smearceptor.second, &smearceptorkey);
return smear;
}
DynamicSmearceptorFactory::~DynamicSmearceptorFactory() { Manifests.clear(); }
#endif
Smearcepterton* Smearcepterton::_inst = NULL;
Smearcepterton& Smearcepterton::Get() {
if (!_inst) {
_inst = new Smearcepterton();
}
return *_inst;
}
Smearcepterton::Smearcepterton() { InitialiserSmearcepters(); }
void Smearcepterton::InitialiserSmearcepters() {
// hard coded list of tag name -> smearcepter factories, add here to add your
// own.
std::map factories;
factories["ThresholdAccepter"] = &BuildSmearcepter;
factories["EfficiencyApplicator"] = &BuildSmearcepter;
factories["GaussianSmearer"] = &BuildSmearcepter;
factories["TrackedMomentumMatrixSmearer"] =
&BuildSmearcepter;
factories["VisECoalescer"] = &BuildSmearcepter;
factories["MetaSimpleSmearcepter"] = &BuildSmearcepter;
+ Config::Get().PrintXML(NULL);
+
std::vector smearcepterBlocks = Config::QueryKeys("smearcepters");
+ // std::cout << "[INFO]: " << nodelist.size() << " smearcepter nodes." << std::endl;
+ // for(size_t i = 0; i < nodelist.size(); ++i){
+ // Config::Get().PrintXML(nodelist[i]);
+ // }
+
for (size_t smearB_it = 0; smearB_it < smearcepterBlocks.size();
++smearB_it) {
std::vector smearcepters =
smearcepterBlocks[smearB_it].GetListOfChildNodes();
for (size_t smear_it = 0; smear_it < smearcepters.size(); ++smear_it) {
std::string const& smearType = smearcepters[smear_it].GetElementName();
ISmearcepter* smearer = NULL;
#ifdef __USE_DYNSAMPLES__
if (DynamicSmearceptorFactory::Get().HasSmearceptor(smearType)) {
smearer = DynamicSmearceptorFactory::Get().CreateSmearceptor(
smearcepters[smear_it]);
} else
#endif
{
if (!factories.count(smearType)) {
ERROR(WRN, "No known smearer accepts elements named: \"" << smearType
<< "\"");
continue;
}
smearer = factories[smearType](smearcepters[smear_it]);
}
if (!smearer) {
THROW("Failed to load smearceptor.");
}
if (!smearer->GetName().length()) {
THROW("Smearcepter type " << smearer->GetElementName()
<< " had no instance name.");
}
Smearcepters[smearer->GetName()] = smearer;
QLOG(FIT, "Configured smearer named: " << smearer->GetName()
<< " of type: "
<< smearer->GetElementName());
}
}
}
diff --git a/src/Smearceptance/ThresholdAccepter.cxx b/src/Smearceptance/ThresholdAccepter.cxx
index adc7239..7cb716d 100644
--- a/src/Smearceptance/ThresholdAccepter.cxx
+++ b/src/Smearceptance/ThresholdAccepter.cxx
@@ -1,330 +1,330 @@
// 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 "ThresholdAccepter.h"
namespace {
ThresholdAccepter::KineVar GetKineType(nuiskey &nk) {
if (nk.Has("RecoThresholdMomentum_MeV")) {
return ThresholdAccepter::kMomentum;
} else if (nk.Has("RecoThresholdKE_MeV")) {
return ThresholdAccepter::kKE;
} else if (nk.Has("RecoThresholdCosTheta_Max")) {
return ThresholdAccepter::kCosTheta_Max;
} else if (nk.Has("RecoThresholdCosTheta_Min")) {
return ThresholdAccepter::kCosTheta_Min;
} else if (nk.Has("RecoThresholdAbsCosTheta_Max")) {
return ThresholdAccepter::kAbsCosTheta_Max;
} else if (nk.Has("RecoThresholdAbsCosTheta_Min")) {
return ThresholdAccepter::kAbsCosTheta_Min;
} else {
THROW("Cannot determine the threshold type for Smearcepter element.");
}
return ThresholdAccepter::kNoVar;
}
std::string GetKineTypeName(ThresholdAccepter::KineVar kv) {
switch (kv) {
case ThresholdAccepter::kMomentum:
return "Momentum";
case ThresholdAccepter::kKE:
return "KE";
case ThresholdAccepter::kCosTheta_Max:
return "CosTheta_Max";
case ThresholdAccepter::kCosTheta_Min:
return "CosTheta_Min";
case ThresholdAccepter::kAbsCosTheta_Max:
return "AbsCosTheta_Max";
case ThresholdAccepter::kAbsCosTheta_Min:
return "CosTheta_Min";
default:
return "NoVar";
}
}
double GetKineThreshold(nuiskey &nk, ThresholdAccepter::KineVar kv) {
switch (kv) {
case ThresholdAccepter::kMomentum:
return nk.GetD("RecoThresholdMomentum_MeV");
case ThresholdAccepter::kKE:
return nk.GetD("RecoThresholdKE_MeV");
case ThresholdAccepter::kCosTheta_Max:
return nk.GetD("RecoThresholdCosTheta_Max");
case ThresholdAccepter::kCosTheta_Min:
return nk.GetD("RecoThresholdCosTheta_Min");
case ThresholdAccepter::kAbsCosTheta_Max:
return nk.GetD("RecoThresholdAbsCosTheta_Max");
case ThresholdAccepter::kAbsCosTheta_Min:
return nk.GetD("RecoThresholdAbsCosTheta_Min");
default:
return 0;
}
}
double GetKineVal(FitParticle *fp, ThresholdAccepter::Thresh &rt) {
switch (rt.ThresholdType) {
case ThresholdAccepter::kMomentum:
return fp->P3().Mag();
case ThresholdAccepter::kKE:
return fp->KE();
case ThresholdAccepter::kCosTheta_Max:
return fp->P3().CosTheta();
case ThresholdAccepter::kCosTheta_Min:
return fp->P3().CosTheta();
case ThresholdAccepter::kAbsCosTheta_Max:
return fabs(fp->P3().CosTheta());
case ThresholdAccepter::kAbsCosTheta_Min:
return fabs(fp->P3().CosTheta());
default:
return 0;
}
}
bool PassesThreshold(FitParticle *fp, ThresholdAccepter::Thresh &rt) {
switch (rt.ThresholdType) {
case ThresholdAccepter::kMomentum:
return (fp->P3().Mag() > rt.ThresholdVal);
case ThresholdAccepter::kKE:
return (fp->KE() > rt.ThresholdVal);
case ThresholdAccepter::kCosTheta_Max:
return (fp->P3().CosTheta() < rt.ThresholdVal);
case ThresholdAccepter::kCosTheta_Min:
return (fp->P3().CosTheta() > rt.ThresholdVal);
case ThresholdAccepter::kAbsCosTheta_Max:
return (fabs(fp->P3().CosTheta()) < rt.ThresholdVal);
case ThresholdAccepter::kAbsCosTheta_Min:
return (fabs(fp->P3().CosTheta()) > rt.ThresholdVal);
default:
return 0;
}
}
}
/// Reads particle threshold nodes
///
/// Nodes look like:
///
///
///
///
///
///
///
///
void ThresholdAccepter::SpecifcSetup(nuiskey &nk) {
std::vector recoThresholdDescriptors =
nk.GetListOfChildNodes("RecoThreshold");
for (size_t t_it = 0; t_it < recoThresholdDescriptors.size(); ++t_it) {
std::string pdgs_s = recoThresholdDescriptors[t_it].GetS("PDG");
std::vector pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");
for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
Thresh t;
t.ThresholdType = GetKineType(recoThresholdDescriptors[t_it]);
t.ThresholdVal =
GetKineThreshold(recoThresholdDescriptors[t_it], t.ThresholdType);
ReconThresholds[pdgs_i[pdg_it]].push_back(t);
QLOG(FIT, "Added reconstruction threshold of type: "
<< ReconThresholds[pdgs_i[pdg_it]].back().ThresholdVal
<< " "
<< GetKineTypeName(
ReconThresholds[pdgs_i[pdg_it]].back().ThresholdType)
<< ", for PDG: " << pdgs_i[pdg_it]);
}
}
std::vector visThresholdDescriptors =
nk.GetListOfChildNodes("VisThreshold");
for (size_t t_it = 0; t_it < visThresholdDescriptors.size(); ++t_it) {
std::string pdgs_s = visThresholdDescriptors[t_it].GetS("PDG");
std::vector pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ",");
for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) {
if (VisThresholds.count(pdgs_i[pdg_it])) {
ERROR(WRN, "Smearceptor " << ElementName << ":" << InstanceName
<< " already has a threshold for PDG: "
<< pdgs_i[pdg_it]);
}
VisThresh vt;
vt.UseKE = visThresholdDescriptors[t_it].Has("Contrib")
? (visThresholdDescriptors[t_it].GetS("Contrib") == "K")
: false;
vt.Fraction = visThresholdDescriptors[t_it].Has("Fraction")
? visThresholdDescriptors[t_it].GetD("Fraction")
: 1;
if (visThresholdDescriptors[t_it].Has("VisThresholdKE_MeV")) {
vt.ThresholdType = ThresholdAccepter::kKE;
vt.ThresholdVal =
visThresholdDescriptors[t_it].GetD("VisThresholdKE_MeV");
} else if (visThresholdDescriptors[t_it].Has(
"VisThresholdMomentum_MeV")) {
vt.ThresholdType = ThresholdAccepter::kMomentum;
vt.ThresholdVal =
visThresholdDescriptors[t_it].GetD("VisThresholdMomentum_MeV");
;
} else {
ERROR(WRN, "Smearceptor "
<< ElementName << ":" << InstanceName
<< " cannot find threshold information for PDG: "
<< pdgs_i[pdg_it]);
continue;
}
VisThresholds[pdgs_i[pdg_it]] = vt;
QLOG(FIT,
"Added visibility threshold of MeV "
<< VisThresholds[pdgs_i[pdg_it]].ThresholdVal << " "
<< GetKineTypeName(VisThresholds[pdgs_i[pdg_it]].ThresholdType)
<< ", for PDG: " << pdgs_i[pdg_it]
<< ". If visible, particle deposits: "
<< (VisThresholds[pdgs_i[pdg_it]].UseKE ? "KE" : "TE"));
}
}
}
void ThresholdAccepter::SmearceptOneParticle(RecoInfo *ri, FitParticle *fp
#ifdef DEBUG_THRESACCEPT
,
size_t p_it
#endif
) {
#ifdef DEBUG_THRESACCEPT
std::cout << std::endl;
std::cout << "[" << p_it << " = " << fp << "]: " << fp->PDG() << ", "
<< fp->Status() << ", " << fp->E() << " -- KE:" << fp->KE()
<< " Mom: " << fp->P3().Mag() << std::flush;
#endif
if (fp->Status() != kFinalState) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Not final state." << std::flush;
#endif
return;
}
if ((ReconThresholds.count(fp->PDG()) + VisThresholds.count(fp->PDG())) ==
0) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Undetectable." << std::flush;
#endif
return;
}
// If no reco thresholds it should fall through to EVis
bool Passes = ReconThresholds[fp->PDG()].size();
bool FailEnergyThresh = !ReconThresholds[fp->PDG()].size();
for (size_t rt_it = 0; rt_it < ReconThresholds[fp->PDG()].size(); ++rt_it) {
bool Passed = PassesThreshold(fp, ReconThresholds[fp->PDG()][rt_it]);
if (!Passed) {
#ifdef DEBUG_THRESACCEPT
std::cout << "\n\t -- Rejected. ("
<< GetKineTypeName(
ReconThresholds[fp->PDG()][rt_it].ThresholdType)
<< " Threshold: "
<< ReconThresholds[fp->PDG()][rt_it].ThresholdVal << " | "
<< GetKineVal(fp, ReconThresholds[fp->PDG()][rt_it]) << ")"
<< std::flush;
+#endif
if ((ReconThresholds[fp->PDG()][rt_it].ThresholdType ==
ThresholdAccepter::kMomentum) ||
(ReconThresholds[fp->PDG()][rt_it].ThresholdType ==
ThresholdAccepter::kKE)) {
FailEnergyThresh = true;
}
-#endif
} else {
#ifdef DEBUG_THRESACCEPT
std::cout << "\n\t -- Accepted. ("
<< GetKineTypeName(
ReconThresholds[fp->PDG()][rt_it].ThresholdType)
<< " Threshold: "
<< ReconThresholds[fp->PDG()][rt_it].ThresholdVal << " | "
<< GetKineVal(fp, ReconThresholds[fp->PDG()][rt_it]) << ")"
<< std::flush;
#endif
}
Passes = Passes && Passed;
}
if (Passes) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Reconstructed." << std::flush;
#endif
ri->RecObjMom.push_back(fp->P3());
ri->RecObjClass.push_back(fp->PDG());
return;
} else if (!FailEnergyThresh) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Failed non-Energy threshold, no chance for EVis."
<< std::flush;
#endif
return;
}
if (((VisThresholds[fp->PDG()].ThresholdType == ThresholdAccepter::kKE) &&
(VisThresholds[fp->PDG()].ThresholdVal <
fp->KE())) // Above KE-style threshold
|| ((VisThresholds[fp->PDG()].ThresholdType ==
ThresholdAccepter::kMomentum) &&
(VisThresholds[fp->PDG()].ThresholdVal <
fp->P3().Mag())) // Above mom-style threshold
) {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Contributed to VisE. ("
<< GetKineTypeName(VisThresholds[fp->PDG()].ThresholdType) << ": "
<< VisThresholds[fp->PDG()].ThresholdVal << ")" << std::flush;
#endif
ri->RecVisibleEnergy.push_back(
VisThresholds[fp->PDG()].Fraction *
(VisThresholds[fp->PDG()].UseKE ? fp->KE() : fp->E()));
ri->TrueContribPDGs.push_back(fp->PDG());
return;
} else {
#ifdef DEBUG_THRESACCEPT
std::cout << " -- Rejected. "
<< " Vis: ("
<< GetKineTypeName(VisThresholds[fp->PDG()].ThresholdType) << ": "
<< VisThresholds[fp->PDG()].ThresholdVal << ")" << std::flush;
#endif
}
}
RecoInfo *ThresholdAccepter::Smearcept(FitEvent *fe) {
RecoInfo *ri = new RecoInfo();
for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) {
FitParticle *fp = fe->GetParticle(p_it);
SmearceptOneParticle(ri, fp
#ifdef DEBUG_THRESACCEPT
,
p_it
#endif
);
}
#ifdef DEBUG_THRESACCEPT
std::cout << std::endl;
#endif
return ri;
}