diff --git a/data/T2K/CCinc/nue_2019/FHC_nue_unfold_with_neut.txt b/data/T2K/CCinc/nue_2019/FHC_nue_unfold_with_neut.txt
new file mode 100644
index 0000000..61d7536
--- /dev/null
+++ b/data/T2K/CCinc/nue_2019/FHC_nue_unfold_with_neut.txt
@@ -0,0 +1,4 @@
+0.3 3.10407E-39 0.793005E-39
+1.6 2.41606E-39 0.759596E-39
+3.2 1.09267E-39 0.510421E-39
+30.0 0.00 0.00
diff --git a/data/T2K/CCinc/nue_2019/RHC_nue_unfold_with_neut.txt b/data/T2K/CCinc/nue_2019/RHC_nue_unfold_with_neut.txt
new file mode 100644
index 0000000..34d8b89
--- /dev/null
+++ b/data/T2K/CCinc/nue_2019/RHC_nue_unfold_with_neut.txt
@@ -0,0 +1,4 @@
+0.3 3.97726E-39 1.88965E-39
+1.6 10.5856E-39 3.08965E-39
+30.0 0.00 0.00
+
diff --git a/data/T2K/CCinc/nue_2019/RHC_nuebar_unfold_with_neut.txt b/data/T2K/CCinc/nue_2019/RHC_nuebar_unfold_with_neut.txt
new file mode 100644
index 0000000..fa4ba4e
--- /dev/null
+++ b/data/T2K/CCinc/nue_2019/RHC_nuebar_unfold_with_neut.txt
@@ -0,0 +1,3 @@
+0.3 1.01288E-39 0.551926E-39
+1.6 1.99471E-39 0.768926E-39
+30.0 0.00 0.00
diff --git a/data/T2K/CCinc/nue_2019/fract_covar_unfold_with_neut.txt b/data/T2K/CCinc/nue_2019/fract_covar_unfold_with_neut.txt
new file mode 100644
index 0000000..e90f33a
--- /dev/null
+++ b/data/T2K/CCinc/nue_2019/fract_covar_unfold_with_neut.txt
@@ -0,0 +1,11 @@
+0.0716895 0.0166114 0.065507 0.0232243 0.0315353 0.0398808 0.0347689
+0.0166114 0.0935784 0.0262206 0.0337246 0.0293716 0.0431473 0.0384238
+0.065507 0.0262206 0.230568 0.0585712 0.0650908 0.0850014 0.0798983
+0.0232243 0.0337246 0.0585712 0.242144 0.00368026 0.0536218 0.0440962
+0.0315353 0.0293716 0.0650908 0.00368026 0.0832293 0.0450185 0.0400094
+0.0398808 0.0431473 0.0850014 0.0536218 0.0450185 0.292844 0.0440944
+0.0347689 0.0384238 0.0798983 0.0440962 0.0400094 0.0440944 0.151416
+
+
+
+
diff --git a/parameters/config.xml b/parameters/config.xml
index 05044e1..24db849 100644
--- a/parameters/config.xml
+++ b/parameters/config.xml
@@ -1,233 +1,233 @@
-
-
+
diff --git a/parameters/fitter.config.dat b/parameters/fitter.config.dat
index e8ab1ba..a0c71d0 100644
--- a/parameters/fitter.config.dat
+++ b/parameters/fitter.config.dat
@@ -1,158 +1,158 @@
# ######################################################
# NUISANCE CONFIGURATION OPTIONS
# This file is read in by default at runtime
# If you want to override on a case by case bases use -q at runtime
# ######################################################
# MAIN Configs
# ######################################################
# Logger goes from
# 1 Quiet
# 2 Fitter
# 3 Samples
# 4 Reconfigure Loops
# 5 Every Event print out (SHOUT)
# -1 DEBUGGING
config verbosity 2
config VERBOSITY 2
# ERROR goes from
# 0 NONE
# 1 FATAL
# 2 WARN
config ERROR 2
# Input Configs
# ######################################################
# Default Requirements file for the externalDataFitter Package
# MAX Events : -1 is no cut. Events will be scaled automatically to give good xsec predictions.
config input.maxevents -1
config MAXEVENTS -1
config input.MAXEVENTS -1
# Turn on/off event manager
# EventManager enables us to only loop number of events once for multiple projections of the same measurements
# e.g. MiniBooNE CC1pi+ Q2 and MiniBooNE CC1pi+ Tmu would ordinarily require 2 reconfigures, but with this enabled it requires only one
config input.eventmanager 0
config EventManager 0
# Event Directories
# Can setup default directories and use @EVENT_DIR/path to link to it
config EVENT_DIR /data2/stowell/NIWG/
config NEUT_EVENT_DIR /data2/stowell/NIWG/neut/fit_samples_neut5.3.3/
config GENIE_EVENT_DIR /data2/stowell/NIWG/genie/fit_samples_R.2.10.0/
config NUWRO_EVENT_DIR /data2/stowell/NIWG/nuwro/fit_samples/
config GIBUU_EVENT_DIR /data/GIBUU/DIR/
# In PrepareGENIE the reconstructed splines can be saved into the file
config save_genie_splines 1
# In InputHandler the option to regenerate NuWro flux/xsec plots is available
# Going to move this to its own app soon
config input.regen_nuwro_plots 0
# DEVEL CONFIG OPTION, don't touch!
config cachesize 0
# ReWeighting Configuration Options
# ######################################################
# Set absolute twkdial for parameters
config params.setabstwk 0
# Convert Dials in output statements using dial conversion card
config convert_dials 0
# Make RW Calculations be quiet
condif params.silentweighting 0
# Vetos can be used to specify RW dials NOT to be loaded in
# Useful if one specific one has an issue
config FitWeight.fNIWGRW_veto ''
config FitWeight.fNuwroRW_veto ''
config FitWeight.fNeutRW_veto ''
config FitWeight.fGenieRW_veto ''
# Output Options
# ######################################################
# Save Nominal prediction with all rw engines at default
config savenominal 0
# Save prefit with values at starting values
config saveprefit 0
# Here's the full list of drawing options
# See src/FitBase/Measurement1D::Write for more info
#config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/RESIDUAL/MATRIX/FLUX/MASK/MAP
#config drawopts DATA/MC
-config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/FLUX/XSEC/MASK/COV/INCOV/DECOMP/CANVPDG/CANVMC
+config drawopts DATA/MC/EVT/FINE/RATIO/MODES/SHAPE/FLUX/XSEC/MASK/COV/INVCOV/DECOMP/CANVPDG/CANVMC
# Save the shape scaling applied with option SHAPE into the main MC hist
config saveshapescaling 0
# Set style of 1D output histograms
config linecolour 1
config linestyle 1
config linewidth 1
# For GenericFlux
config isLiteMode 0
# Statistical Options
# ######################################################
# Add MC Statistical error to likelihoods
config statutils.addmcerror 0
# NUISMIN Configurations
# ######################################################
config minimizer.maxcalls 1000000
config minimizer.maxiterations 1000000
config minimizer.tolerance 0.001
# Number of events required in low stats routines
config minimizer.lowstatevents 25000
# Error band generator configs
# ######################################################
# For -f ErrorBands creates error bands for given measurements
# How many throws do we want (The higher the better precision)
config error_throws 250
# Are we throwing uniform or according to Gaussian?
# Only use uniform if wanting to study the limits of a dial.
config error_uniform 0
# Other Individual Case Configs
# ######################################################
# Covariance throw options for fake data studies with MiniBooNE data.
config thrown_covariance FULL
config throw_mc_stat 0.0
config throw_diag_syst 0
config throw_corr_syst 0
config throw_mc_stat 0
# Apply a shift to the muon momentum before calculation of Q2
config muon_momentum_shift 0.0
config muon_momentum_throw 0
# MINERvA Specific Configs
config MINERvA_XSec_CCinc_2DEavq3_nu.hadron_cut 0
config MINERvA_CCinc_XSec_2DEavq3_nu.useq3true 0
config Modes.split_PN_NN 0
# SciBooNE specific
config SciBarDensity 1.04
#config SciBarDensity 2.00
config SciBarRecoDist 10.0
config PenetratingMuonEnergy 1.4
config NumRangeSteps 50
diff --git a/src/Config/NuisConfig.cxx b/src/Config/NuisConfig.cxx
index c1d3915..4a33b29 100644
--- a/src/Config/NuisConfig.cxx
+++ b/src/Config/NuisConfig.cxx
@@ -1,817 +1,817 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#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;
+ std::cout << "[ NUISANCE ]: Loading DEFAULT settings from : " << filename << std::endl;
// Create XML Engine
fXML = new TXMLEngine;
fXML->SetSkipComments(true);
// Load in documents
fXMLDocs.push_back(fXML->ParseFile(filename.c_str(), 1000000));
if (!fXMLDocs[0]) {
NUIS_ABORT("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;
+ std::cout << "[ NUISANCE ]: Finished nuisconfig." << 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()) {
NUIS_ABORT("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) {
NUIS_ABORT("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;
+ std::cout << "[ NUISANCE ]: Finished loading XML settings" << std::endl;
}
void nuisconfig::LoadCardSettings(std::string const &filename,
std::string const &state) {
- std::cout << "[ NUISANCE ]: Loading simple config from : " << filename;
+ std::cout << "[ NUISANCE ]: Loading simple config from : " << filename << std::endl;
// 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;
+ std::cout << "[ NUISANCE ]: Finished loading simple config" << 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";
+ std::cout << "[ NUISANCE ]: Finalising run settings" << std::endl;
WriteSettings(name);
// Save full config to file
RemoveEmptyNodes();
RemoveIdenticalNodes();
- std::cout << " -> DONE." << std::endl;
+ std::cout << "[ NUISANCE ]: Finished finalising run settings" << std::endl;
}
void nuisconfig::WriteSettings(std::string const &outputname) {
// Create a New XML Doc
XMLDocPointer_t newxmldoc = fXML->NewDoc();
fXML->DocSetRootElement(newxmldoc, fMainNode);
// Save document to file
if (GetConfigB("SaveParsedXMLFile")) {
fXML->SaveDoc(newxmldoc, outputname.c_str());
}
}
void nuisconfig::PrintXML(XMLNodePointer_t node, int indent) {
if (!node) {
node = fMainNode;
}
std::stringstream ss("");
for (int i = 0; i < indent; ++i) {
ss << " ";
}
std::cout << ss.str() << "<" << fXML->GetNodeName(node) << std::flush;
XMLAttrPointer_t attr = fXML->GetFirstAttr(node);
while (attr) {
std::cout << " " << fXML->GetAttrName(attr) << "=\""
<< fXML->GetAttrValue(attr) << "\"" << std::flush;
attr = fXML->GetNextAttr(attr);
}
if (!fXML->GetChild(node)) {
std::cout << " />" << std::endl;
return;
}
std::cout << " >" << std::endl;
XMLNodePointer_t child = fXML->GetChild(node);
while (child) {
PrintXML(child, indent + 1);
child = fXML->GetNext(child);
}
std::cout << ss.str() << "" << fXML->GetNodeName(node) << ">" << std::endl;
}
XMLNodePointer_t nuisconfig::CreateNode(std::string const &name) {
return fXML->NewChild(fMainNode, 0, name.c_str());
}
XMLNodePointer_t nuisconfig::CreateNode(XMLNodePointer_t node,
std::string const &name) {
return fXML->NewChild(node, 0, name.c_str());
}
XMLNodePointer_t nuisconfig::GetNode(XMLNodePointer_t node,
std::string const &type) {
/// Loop over all children
XMLNodePointer_t child = fXML->GetChild(node);
while (child != 0) {
/// Get nodes for given type (if type empty return all)
if (std::string(fXML->GetNodeName(child)) == type.c_str() or type.empty()) {
return child;
}
// Next child
child = fXML->GetNext(child);
}
// Child not found
return 0;
}
void nuisconfig::RemoveEmptyNodes() {
std::vector 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 93aff01..fb63bef 100644
--- a/src/FCN/SampleList.cxx
+++ b/src/FCN/SampleList.cxx
@@ -1,1435 +1,1447 @@
#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_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 CC1Pi
#include "ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dpmu_nu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu.h"
// 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_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_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_MicroBooNE__
#include "MicroBooNE_CCInc_XSec_2DPcos_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"
// 2018 MINERvA CC0pi STV
#include "MINERvA_CC0pinp_STV_XSec_1D_nu.h"
// 2018 MINERvA CC0pi 2D
#include "MINERvA_CC0pi_XSec_1D_2018_nu.h"
#include "MINERvA_CC0pi_XSec_2D_nu.h"
// #include "MINERvA_CC0pi_XSec_3DptpzTp_nu.h"
// 2018 MINERvA CC0pi 2D antinu
#include "MINERvA_CC0pi_XSec_2D_antinu.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_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_joint.h"
#include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h"
#include "MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h"
#endif
#ifndef __NO_T2K__
// T2K CC0pi 2016
#include "T2K_CC0pi_XSec_2DPcos_nu_I.h"
#include "T2K_CC0pi_XSec_2DPcos_nu_II.h"
// T2K CC-inclusive with full acceptance 2018
#include "T2K_CCinc_XSec_2DPcos_nu_nonuniform.h"
+// T2K nue CC-inclusive 2019
+#include "T2K_nueCCinc_XSec_1Dpe.h"
+#include "T2K_nueCCinc_XSec_joint.h"
+
// T2K STV CC0pi 2018
#include "T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h"
#include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h"
#include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h"
#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h"
#include "T2K_CC0pinp_ifk_XSec_3Dinfa_nu.h"
#include "T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h"
#include "T2K_CC0pinp_ifk_XSec_3Dinfp_nu.h"
// T2K CC1pi+ on CH
#include "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h"
#include "T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h"
#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dppi_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h"
#include "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h"
//#include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h"
//#include "T2K_CC1pip_CH_XSec_1DWrec_nu.h"
//#include "T2K_CC1pip_CH_XSec_1Dq3_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"
#endif
#ifndef __NO_SciBooNE__
// SciBooNE COH studies
#include "SciBooNE_CCCOH_1TRK_1DQ2_nu.h"
#include "SciBooNE_CCCOH_1TRK_1Dpmu_nu.h"
#include "SciBooNE_CCCOH_1TRK_1Dthetamu_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_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_MuPiVA_1Dpmu_nu.h"
#include "SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu.h"
#include "SciBooNE_CCCOH_MuPr_1DQ2_nu.h"
#include "SciBooNE_CCCOH_MuPr_1Dpmu_nu.h"
#include "SciBooNE_CCCOH_MuPr_1Dthetamu_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();
NUIS_LOG(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]);
NUIS_LOG(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)) {
NUIS_LOG(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()) {
NUIS_ERR(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()) {
NUIS_ERR(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()) {
NUIS_ERR(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()) {
NUIS_ERR(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()) {
NUIS_ERR(WRN, "Failed to load symbol \"DSF_DestroySample\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.NSamples = (*(plgManif.DSF_NSamples))();
NUIS_LOG(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) {
NUIS_ABORT("Could not load sample "
<< smp_it << " / " << plgManif.NSamples << " from "
<< plgManif.soloc);
}
if (Samples.count(smp_name)) {
NUIS_ERR(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);
NUIS_LOG(FIT, "\t\t" << smp_name);
}
if (plgManif.SamplesProvided.size()) {
Manifests[plgManif.soloc] = plgManif;
NSamples += plgManif.SamplesProvided.size();
NManifests++;
} else {
dlclose(dlobj);
}
}
}
closedir(dir);
} else {
NUIS_ERR(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);
}
NUIS_LOG(FIT, "Dynamic sample manifest: ");
for (std::map >::iterator m_it =
ManifestSamples.begin();
m_it != ManifestSamples.end(); ++m_it) {
NUIS_LOG(FIT, "\tLibrary " << m_it->first << " contains: ");
for (size_t s_it = 0; s_it < m_it->second.size(); ++s_it) {
NUIS_LOG(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)) {
NUIS_ERR(WRN, "Asked to load unknown sample: \"" << samplekey.GetS("name")
<< "\".");
return NULL;
}
std::pair sample = Samples[samplekey.GetS("name")];
NUIS_LOG(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)) {
NUIS_LOG(SAM, "Instantiating dynamic sample...");
MeasurementBase *ds = DynamicSampleFactory::Get().CreateSample(samplekey);
if (ds) {
NUIS_LOG(SAM, "Done.");
return ds;
}
NUIS_ABORT("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));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_nu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu(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));
} else
#endif
#ifndef __NO_MicroBooNE__
/*
MicroBooNE Samples
*/
/*
MicroBooNE CCinclusive
*/
if (!name.compare("MicroBooNE_CCInc_XSec_2DPcos_nu")) {
return (new MicroBooNE_CCInc_XSec_2DPcos_nu(samplekey));
} else
#endif
#ifndef __NO_MINERvA__
/*
MINERvA Samples
*/
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_CC0pinp_STV_XSec_1Dpmu_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu")) {
return (new MINERvA_CC0pinp_STV_XSec_1D_nu(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));
// Dan Ruterbories measurements of late 2018
} else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu")) {
return (new MINERvA_CC0pi_XSec_2D_nu(samplekey));
// } else if (!name.compare("MINERvA_CC0pi_XSec_3DptpzTp_nu")) {
// return (new MINERvA_CC0pi_XSec_3DptpzTp_nu(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1Dpt_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1Dpz_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2QE_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DEnuQE_nu")) {
return (new MINERvA_CC0pi_XSec_1D_2018_nu(samplekey));
// C. Patrick's early 2018 measurements
} else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_antinu") ||
!name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuQE_antinu") ||
!name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuTrue_antinu")) {
return (new MINERvA_CC0pi_XSec_2D_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));
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) {
return (new MINERvA_CCNpip_XSec_1Dthmu_nu(samplekey));
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) {
return (new MINERvA_CCNpip_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) {
return (new MINERvA_CCNpip_XSec_1DQ2_nu(samplekey));
} 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_I")) {
return (new T2K_CC0pi_XSec_2DPcos_nu_I(samplekey));
} else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) {
return (new T2K_CC0pi_XSec_2DPcos_nu_II(samplekey));
} else if (!name.compare("T2K_CCinc_XSec_2DPcos_nu_nonuniform")) {
return (new T2K_CCinc_XSec_2DPcos_nu_nonuniform(samplekey));
+ } else if (!name.compare("T2K_nueCCinc_XSec_1Dpe_FHC") ||
+ !name.compare("T2K_nueCCinc_XSec_1Dpe_RHC") ||
+ !name.compare("T2K_nuebarCCinc_XSec_1Dpe_RHC")) {
+ return (new T2K_nueCCinc_XSec_1Dpe(samplekey));
+
+ } else if (!name.compare("T2K_nueCCinc_XSec_joint")) {
+ return (new T2K_nueCCinc_XSec_joint(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_2Dpmucosmu_nu")) {
return (new T2K_CC1pip_CH_XSec_2Dpmucosmu_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_1Dthpi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthpi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DQ2_nu")) {
return (new T2K_CC1pip_CH_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DAdlerPhi_nu")) {
return (new T2K_CC1pip_CH_XSec_1DAdlerPhi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DCosThAdler_nu")) {
return (new T2K_CC1pip_CH_XSec_1DCosThAdler_nu(samplekey));
// Maybe something for the future: were in Raquel's thesis
//} 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_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));
} else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddphit_nu")) {
return (new T2K_CC0pinp_STV_XSec_1Ddphit_nu(samplekey));
} else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddat_nu")) {
return (new T2K_CC0pinp_STV_XSec_1Ddat_nu(samplekey));
} else if (!name.compare("T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform")) {
return (new T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform(samplekey));
} else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfp_nu")) {
return (new T2K_CC0pinp_ifk_XSec_3Dinfp_nu(samplekey));
} else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfa_nu")) {
return (new T2K_CC0pinp_ifk_XSec_3Dinfa_nu(samplekey));
} else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfip_nu")) {
return (new T2K_CC0pinp_ifk_XSec_3Dinfip_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_1TRK_1Dpmu_nu")) {
return (new SciBooNE_CCCOH_1TRK_1Dpmu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_1TRK_1Dthetamu_nu")) {
return (new SciBooNE_CCCOH_1TRK_1Dthetamu_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_MuPr_1Dpmu_nu")) {
return (new SciBooNE_CCCOH_MuPr_1Dpmu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPr_1Dthetamu_nu")) {
return (new SciBooNE_CCCOH_MuPr_1Dthetamu_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_MuPiVA_1Dpmu_nu")) {
return (new SciBooNE_CCCOH_MuPiVA_1Dpmu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu")) {
return (new SciBooNE_CCCOH_MuPiVA_1Dthetamu_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_MuPiNoVA_1Dthetamu_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1Dpmu_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 {
NUIS_ABORT("Error: No such sample: " << name << std::endl);
}
// Return NULL if no sample loaded.
return NULL;
}
} // namespace SampleUtils
diff --git a/src/FitBase/JointMeas1D.cxx b/src/FitBase/JointMeas1D.cxx
index 7d7afec..93c4abd 100644
--- a/src/FitBase/JointMeas1D.cxx
+++ b/src/FitBase/JointMeas1D.cxx
@@ -1,2222 +1,2222 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This ile is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "JointMeas1D.h"
//********************************************************************
JointMeas1D::JointMeas1D(void) {
//********************************************************************
// XSec Scalings
fScaleFactor = -1.0;
fCurrentNorm = 1.0;
// Histograms
fDataHist = NULL;
fDataTrue = NULL;
fMCHist = NULL;
fMCFine = NULL;
fMCWeighted = NULL;
fMaskHist = NULL;
// Covar
covar = NULL;
fFullCovar = NULL;
fCovar = NULL;
fInvert = NULL;
fDecomp = NULL;
// Fake Data
fFakeDataInput = "";
fFakeDataFile = NULL;
// Options
fDefaultTypes = "FIX/FULL/CHI2";
fAllowedTypes =
"FIX,FREE,SHAPE/FULL,DIAG/CHI2/NORM/ENUCORR/Q2CORR/ENU1D/MASK";
fIsFix = false;
fIsShape = false;
fIsFree = false;
fIsDiag = false;
fIsFull = false;
fAddNormPen = false;
fIsMask = false;
fIsChi2SVD = false;
fIsRawEvents = false;
fIsDifXSec = false;
fIsEnu1D = false;
// Inputs
fInput = NULL;
fRW = NULL;
// Extra Histograms
fMCHist_Modes = NULL;
for (std::vector::const_iterator iter = fSubChain.begin();
iter != fSubChain.end(); iter++) {
MeasurementBase *exp = *iter;
if (exp)
delete exp;
}
fSubChain.clear();
// Flags for Joint Measurements
fIsRatio = false;
fIsSummed = false;
fSaveSubMeas = false;
fIsJoint = true;
}
//********************************************************************
void JointMeas1D::SetupDefaultHist() {
//********************************************************************
// Setup fMCHist
fMCHist = (TH1D *)fDataHist->Clone();
fMCHist->SetNameTitle((fName + "_MC").c_str(),
(fName + "_MC" + fPlotTitles).c_str());
// Setup fMCFine
Int_t nBins = fMCHist->GetNbinsX();
fMCFine = new TH1D(
(fName + "_MC_FINE").c_str(), (fName + "_MC_FINE" + fPlotTitles).c_str(),
nBins * 6, fMCHist->GetBinLowEdge(1), fMCHist->GetBinLowEdge(nBins + 1));
fMCStat = (TH1D *)fMCHist->Clone();
fMCStat->Reset();
fMCHist->Reset();
fMCFine->Reset();
// Setup the NEUT Mode Array
PlotUtils::CreateNeutModeArray((TH1D *)fMCHist, (TH1 **)fMCHist_PDG);
PlotUtils::ResetNeutModeArray((TH1 **)fMCHist_PDG);
// Setup bin masks using sample name
if (fIsMask) {
std::string maskloc = FitPar::Config().GetParDIR(fName + ".mask");
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + fName + ".mask";
}
SetBinMask(maskloc);
}
fMCHist_Modes =
new TrueModeStack((fName + "_MODES").c_str(), ("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes);
return;
}
//********************************************************************
JointMeas1D::~JointMeas1D(void) {
//********************************************************************
if (fDataHist)
delete fDataHist;
if (fDataTrue)
delete fDataTrue;
if (fMCHist)
delete fMCHist;
if (fMCFine)
delete fMCFine;
if (fMCWeighted)
delete fMCWeighted;
if (fMaskHist)
delete fMaskHist;
if (covar)
delete covar;
if (fFullCovar)
delete fFullCovar;
if (fCovar)
delete fCovar;
if (fInvert)
delete fInvert;
if (fDecomp)
delete fDecomp;
for (std::vector::const_iterator iter = fSubChain.begin();
iter != fSubChain.end(); iter++) {
MeasurementBase *exp = *iter;
if (exp)
delete exp;
}
fSubChain.clear();
}
//********************************************************************
SampleSettings JointMeas1D::LoadSampleSettings(nuiskey samplekey) {
//********************************************************************
SampleSettings s = MeasurementBase::LoadSampleSettings(samplekey);
// Parse Inputs
fSubInFiles.clear();
std::vector entries =
GeneralUtils::ParseToStr(s.GetS("input"), ";");
if (entries.size() < 2) {
NUIS_ABORT("Joint measurement expected to recieve at least two semi-colon "
"separated input files, but recieved: \""
<< s.GetS("input") << "\"");
}
std::vector first_file_descriptor =
GeneralUtils::ParseToStr(entries.front(), ":");
if (first_file_descriptor.size() != 2) {
NUIS_ABORT("Found Joint measurement where the input file had no type: \""
<< s.GetS("input")
<< "\", expected \"INPUTTYPE:File.root;File2.root\".");
}
std::string inpType = first_file_descriptor[0];
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
if (GeneralUtils::ParseToStr(*iter, ":").size() != 2) {
std::stringstream ss("");
ss << inpType << ":" << (*iter);
fSubInFiles.push_back(ss.str());
} else {
fSubInFiles.push_back(*iter);
}
}
return s;
}
//********************************************************************
void JointMeas1D::FinaliseSampleSettings() {
//********************************************************************
// Setup naming + renaming
fName = fSettings.GetName();
fSettings.SetS("originalname", fName);
if (fSettings.Has("rename")) {
fName = fSettings.GetS("rename");
fSettings.SetS("name", fName);
}
// Setup all other options
NUIS_LOG(SAM, "Finalising Sample Settings: " << fName);
if ((fSettings.GetS("originalname").find("Evt") != std::string::npos)) {
fIsRawEvents = true;
NUIS_LOG(SAM, "Found event rate measurement but using poisson likelihoods.");
}
if (fSettings.GetS("originalname").find("XSec_1DEnu") != std::string::npos) {
fIsEnu1D = true;
NUIS_LOG(SAM, "::" << fName << "::");
NUIS_LOG(SAM, "Found XSec Enu measurement, applying flux integrated scaling, "
<< "not flux averaged!");
}
if (fIsEnu1D && fIsRawEvents) {
NUIS_LOG(SAM, "Found 1D Enu XSec distribution AND fIsRawEvents, is this "
"really correct?!");
NUIS_LOG(SAM,
"Check experiment constructor for " << fName << " and correct this!");
NUIS_LOG(SAM, "I live in " << __FILE__ << ":" << __LINE__);
exit(-1);
}
// Parse Inputs
fSubInFiles.clear();
std::vector entries =
GeneralUtils::ParseToStr(fSettings.GetS("input"), ";");
if (entries.size() < 2) {
NUIS_ABORT("Joint measurement expected to recieve at least two semi-colon "
"separated input files, but recieved: \""
<< fSettings.GetS("input") << "\"");
}
std::vector first_file_descriptor =
GeneralUtils::ParseToStr(entries.front(), ":");
if (first_file_descriptor.size() != 2) {
NUIS_ABORT("Found Joint measurement where the input file had no type: \""
<< fSettings.GetS("input")
<< "\", expected \"INPUTTYPE:File.root;File2.root\".");
}
std::string inpType = first_file_descriptor[0];
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
if (GeneralUtils::ParseToStr(*iter, ":").size() != 2) {
std::stringstream ss("");
ss << inpType << ":" << (*iter);
fSubInFiles.push_back(ss.str());
} else {
fSubInFiles.push_back(*iter);
}
}
// Setup options
SetFitOptions(fDefaultTypes); // defaults
SetFitOptions(fSettings.GetS("type")); // user specified
EnuMin = GeneralUtils::StrToDbl(fSettings.GetS("enu_min"));
EnuMax = GeneralUtils::StrToDbl(fSettings.GetS("enu_max"));
if (fAddNormPen) {
if (fNormError <= 0.0) {
NUIS_ERR(FTL, "Norm error for class " << fName << " is 0.0!");
NUIS_ERR(FTL, "If you want to use it please add fNormError=VAL");
throw;
}
}
if (!fRW)
fRW = FitBase::GetRW();
NUIS_LOG(SAM, "Finalised Sample Settings");
}
//********************************************************************
void JointMeas1D::SetDataFromTextFile(std::string datafile) {
//********************************************************************
NUIS_LOG(SAM, "Reading data from text file: " << datafile);
fDataHist = PlotUtils::GetTH1DFromFile(
datafile, fSettings.GetName() + "_data", fSettings.GetFullTitles());
}
//********************************************************************
void JointMeas1D::SetDataFromRootFile(std::string datafile,
std::string histname) {
//********************************************************************
NUIS_LOG(SAM, "Reading data from root file: " << datafile << ";" << histname);
fDataHist = PlotUtils::GetTH1DFromRootFile(datafile, histname);
fDataHist->SetNameTitle((fSettings.GetName() + "_data").c_str(),
(fSettings.GetFullTitles()).c_str());
return;
};
//********************************************************************
void JointMeas1D::SetPoissonErrors() {
//********************************************************************
if (!fDataHist) {
NUIS_ERR(FTL, "Need a data hist to setup possion errors! ");
NUIS_ABORT("Setup Data First!");
}
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, sqrt(fDataHist->GetBinContent(i + 1)));
}
}
//********************************************************************
void JointMeas1D::SetCovarFromDiagonal(TH1D *data) {
//********************************************************************
if (!data and fDataHist) {
data = fDataHist;
}
if (data) {
NUIS_LOG(SAM, "Setting diagonal covariance for: " << data->GetName());
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(data);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
} else {
NUIS_ERR(FTL, "No data input provided to set diagonal covar from!");
}
if (!fIsDiag) {
NUIS_ERR(FTL, "SetCovarMatrixFromDiag called for measurement "
<< "that is not set as diagonal.");
throw;
}
}
//********************************************************************
void JointMeas1D::SetCovarFromTextFile(std::string covfile, int dim) {
//********************************************************************
NUIS_LOG(SAM, "Reading covariance from text file: " << covfile);
fFullCovar = StatUtils::GetCovarFromTextFile(covfile, dim);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCovarFromMultipleTextFiles(std::string covfiles, int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
std::vector covList = GeneralUtils::ParseToStr(covfiles, ";");
fFullCovar = new TMatrixDSym(dim);
for (uint i = 0; i < covList.size(); ++i) {
NUIS_LOG(SAM, "Reading covariance from text file: " << covList[i]);
TMatrixDSym *temp_cov = StatUtils::GetCovarFromTextFile(covList[i], dim);
(*fFullCovar) += (*temp_cov);
delete temp_cov;
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCovarFromRootFile(std::string covfile,
std::string histname) {
//********************************************************************
NUIS_LOG(SAM,
"Reading covariance from text file: " << covfile << ";" << histname);
fFullCovar = StatUtils::GetCovarFromRootFile(covfile, histname);
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCovarInvertFromTextFile(std::string covfile, int dim) {
//********************************************************************
NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile);
covar = StatUtils::GetCovarFromTextFile(covfile, dim);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCovarInvertFromRootFile(std::string covfile,
std::string histname) {
//********************************************************************
NUIS_LOG(SAM, "Reading inverted covariance from text file: " << covfile << ";"
<< histname);
covar = StatUtils::GetCovarFromRootFile(covfile, histname);
fFullCovar = StatUtils::GetInvert(covar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCorrelationFromTextFile(std::string covfile, int dim) {
//********************************************************************
if (dim == -1)
dim = fDataHist->GetNbinsX();
NUIS_LOG(SAM,
"Reading data correlations from text file: " << covfile << ";" << dim);
TMatrixDSym *correlation = StatUtils::GetCovarFromTextFile(covfile, dim);
if (!fDataHist) {
NUIS_ABORT("Trying to set correlations from text file but there is no "
"data to build it from. \n"
<< "In constructor make sure data is set before "
"SetCorrelationFromTextFile is called. \n");
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) *
fDataHist->GetBinError(i + 1) *
fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
//********************************************************************
void JointMeas1D::SetCorrelationFromMultipleTextFiles(std::string corrfiles,
int dim) {
//********************************************************************
if (dim == -1) {
dim = fDataHist->GetNbinsX();
}
std::vector corrList = GeneralUtils::ParseToStr(corrfiles, ";");
fFullCovar = new TMatrixDSym(dim);
for (uint i = 0; i < corrList.size(); ++i) {
NUIS_LOG(SAM, "Reading covariance from text file: " << corrList[i]);
TMatrixDSym *temp_cov = StatUtils::GetCovarFromTextFile(corrList[i], dim);
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
// Note that there is a factor of 1E76 here. It is very silly indeed.
// However, if you remove it, you also need to fix the factors of 1E38
// added to the chi2 calculations!
(*temp_cov)(i, j) = (*temp_cov)(i, j) * fDataHist->GetBinError(i + 1) *
fDataHist->GetBinError(j + 1) * 1E76;
}
}
(*fFullCovar) += (*temp_cov);
delete temp_cov;
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
//********************************************************************
void JointMeas1D::SetCorrelationFromRootFile(std::string covfile,
std::string histname) {
//********************************************************************
NUIS_LOG(SAM, "Reading data correlations from text file: " << covfile << ";"
<< histname);
TMatrixDSym *correlation = StatUtils::GetCovarFromRootFile(covfile, histname);
if (!fDataHist) {
NUIS_ABORT("Trying to set correlations from text file but there is no "
"data to build it from. \n"
<< "In constructor make sure data is set before "
"SetCorrelationFromTextFile is called. \n");
}
// Fill covar from data errors and correlations
fFullCovar = new TMatrixDSym(fDataHist->GetNbinsX());
for (int i = 0; i < fDataHist->GetNbinsX(); i++) {
for (int j = 0; j < fDataHist->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = (*correlation)(i, j) *
fDataHist->GetBinError(i + 1) *
fDataHist->GetBinError(j + 1) * 1.E76;
}
}
// Fill other covars.
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete correlation;
}
void JointMeas1D::SetShapeCovar() {
// Return if this is missing any pre-requisites
if (!fFullCovar)
return;
if (!fDataHist)
return;
// Also return if it's bloody stupid under the circumstances
if (fIsDiag)
return;
fShapeCovar = StatUtils::ExtractShapeOnlyCovar(fFullCovar, fDataHist);
return;
}
//********************************************************************
void JointMeas1D::SetCholDecompFromTextFile(std::string covfile, int dim) {
//********************************************************************
NUIS_LOG(SAM, "Reading cholesky from text file: " << covfile);
TMatrixD *temp = StatUtils::GetMatrixFromTextFile(covfile, dim, dim);
TMatrixD *trans = (TMatrixD *)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(dim, trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void JointMeas1D::SetCholDecompFromRootFile(std::string covfile,
std::string histname) {
//********************************************************************
NUIS_LOG(SAM, "Reading cholesky decomp from root file: " << covfile << ";"
<< histname);
TMatrixD *temp = StatUtils::GetMatrixFromRootFile(covfile, histname);
TMatrixD *trans = (TMatrixD *)temp->Clone();
trans->T();
(*trans) *= (*temp);
fFullCovar = new TMatrixDSym(temp->GetNrows(), trans->GetMatrixArray(), "");
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
delete temp;
delete trans;
}
//********************************************************************
void JointMeas1D::ScaleData(double scale) {
//********************************************************************
fDataHist->Scale(scale);
}
//********************************************************************
void JointMeas1D::ScaleCovar(double scale) {
//********************************************************************
(*fFullCovar) *= scale;
(*covar) *= 1.0 / scale;
(*fDecomp) *= sqrt(scale);
}
//********************************************************************
void JointMeas1D::SetBinMask(std::string maskfile) {
//********************************************************************
if (!fIsMask)
return;
NUIS_LOG(SAM, "Reading bin mask from file: " << maskfile);
// Create a mask histogram with dim of data
int nbins = fDataHist->GetNbinsX();
fMaskHist = new TH1I((fSettings.GetName() + "_BINMASK").c_str(),
(fSettings.GetName() + "_BINMASK; Bin; Mask?").c_str(),
nbins, 0, nbins);
std::string line;
std::ifstream mask(maskfile.c_str(), std::ifstream::in);
if (!mask.is_open()) {
NUIS_ABORT(" Cannot find mask file: " << maskfile);
}
while (std::getline(mask >> std::ws, line, '\n')) {
std::vector entries = GeneralUtils::ParseToInt(line, " ");
// Skip lines with poorly formatted lines
if (entries.size() < 2) {
NUIS_LOG(WRN, "JointMeas1D::SetBinMask(), couldn't parse line: " << line);
continue;
}
// The first index should be the bin number, the second should be the mask
// value.
int val = 0;
if (entries[1] > 0)
val = 1;
fMaskHist->SetBinContent(entries[0], val);
}
// Apply masking by setting masked data bins to zero
PlotUtils::MaskBins(fDataHist, fMaskHist);
return;
}
//********************************************************************
void JointMeas1D::FinaliseMeasurement() {
//********************************************************************
NUIS_LOG(SAM, "Finalising Measurement: " << fName);
// Make sure data is setup
if (!fDataHist) {
NUIS_ABORT("No data has been setup inside " << fName << " constructor!");
}
// Make sure covariances are setup
if (!fFullCovar) {
fIsDiag = true;
SetCovarFromDiagonal(fDataHist);
}
if (!covar) {
covar = StatUtils::GetInvert(fFullCovar);
}
if (!fDecomp) {
fDecomp = StatUtils::GetDecomp(fFullCovar);
}
// Push the diagonals of fFullCovar onto the data histogram
// Comment out until scaling is used consistently...
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar, 1E-38);
// Setup fMCHist from data
fMCHist = (TH1D *)fDataHist->Clone();
fMCHist->SetNameTitle((fSettings.GetName() + "_MC").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCHist->Reset();
// Setup fMCFine
fMCFine = new TH1D("mcfine", "mcfine", fDataHist->GetNbinsX(),
fMCHist->GetBinLowEdge(1),
fMCHist->GetBinLowEdge(fDataHist->GetNbinsX() + 1));
fMCFine->SetNameTitle((fSettings.GetName() + "_MC_FINE").c_str(),
(fSettings.GetFullTitles()).c_str());
fMCFine->Reset();
// Setup MC Stat
fMCStat = (TH1D *)fMCHist->Clone();
fMCStat->Reset();
// Search drawopts for possible types to include by default
std::string drawopts = FitPar::Config().GetParS("drawopts");
if (drawopts.find("MODES") != std::string::npos) {
fMCHist_Modes = new TrueModeStack((fSettings.GetName() + "_MODES").c_str(),
("True Channels"), fMCHist);
SetAutoProcessTH1(fMCHist_Modes);
}
// Setup bin masks using sample name
if (fIsMask) {
std::string curname = fName;
std::string origname = fSettings.GetS("originalname");
// Check rename.mask
std::string maskloc = FitPar::Config().GetParDIR(curname + ".mask");
// Check origname.mask
if (maskloc.empty())
maskloc = FitPar::Config().GetParDIR(origname + ".mask");
// Check database
if (maskloc.empty()) {
maskloc = FitPar::GetDataBase() + "/masks/" + origname + ".mask";
}
// Setup Bin Mask
SetBinMask(maskloc);
}
/*
if (fScaleFactor < 0) {
ERR(FTL) << "I found a negative fScaleFactor in " << __FILE__ << ":" <<
__LINE__ << std::endl; ERR(FTL) << "fScaleFactor = " << fScaleFactor <<
std::endl; ERR(FTL) << "EXITING" << std::endl; throw;
}
*/
// Create and fill Weighted Histogram
if (!fMCWeighted) {
fMCWeighted = (TH1D *)fMCHist->Clone();
fMCWeighted->SetNameTitle((fName + "_MCWGHTS").c_str(),
(fName + "_MCWGHTS" + fPlotTitles).c_str());
fMCWeighted->GetYaxis()->SetTitle("Weighted Events");
}
}
//********************************************************************
void JointMeas1D::SetFitOptions(std::string opt) {
//********************************************************************
// Do nothing if default given
if (opt == "DEFAULT")
return;
// CHECK Conflicting Fit Options
std::vector fit_option_allow =
GeneralUtils::ParseToStr(fAllowedTypes, "/");
for (UInt_t i = 0; i < fit_option_allow.size(); i++) {
std::vector fit_option_section =
GeneralUtils::ParseToStr(fit_option_allow.at(i), ",");
bool found_option = false;
for (UInt_t j = 0; j < fit_option_section.size(); j++) {
std::string av_opt = fit_option_section.at(j);
if (!found_option and opt.find(av_opt) != std::string::npos) {
found_option = true;
} else if (found_option and opt.find(av_opt) != std::string::npos) {
NUIS_ABORT("ERROR: Conflicting fit options provided: "
<< opt << std::endl
<< "Conflicting group = " << fit_option_section.at(i)
<< std::endl
<< "You should only supply one of these options in card file.");
}
}
}
// Check all options are allowed
std::vector fit_options_input =
GeneralUtils::ParseToStr(opt, "/");
for (UInt_t i = 0; i < fit_options_input.size(); i++) {
if (fAllowedTypes.find(fit_options_input.at(i)) == std::string::npos) {
NUIS_ERR(FTL, "ERROR: Fit Option '"
<< fit_options_input.at(i)
<< "' Provided is not allowed for this measurement.");
NUIS_ERR(FTL, "Fit Options should be provided as a '/' seperated list "
"(e.g. FREE/DIAG/NORM)");
NUIS_ERR(FTL, "Available options for " << fName << " are '" << fAllowedTypes
<< "'");
throw;
}
}
// Set TYPE
fFitType = opt;
// FIX,SHAPE,FREE
if (opt.find("FIX") != std::string::npos) {
fIsFree = fIsShape = false;
fIsFix = true;
} else if (opt.find("SHAPE") != std::string::npos) {
fIsFree = fIsFix = false;
fIsShape = true;
} else if (opt.find("FREE") != std::string::npos) {
fIsFix = fIsShape = false;
fIsFree = true;
}
// DIAG,FULL (or default to full)
if (opt.find("DIAG") != std::string::npos) {
fIsDiag = true;
fIsFull = false;
} else if (opt.find("FULL") != std::string::npos) {
fIsDiag = false;
fIsFull = true;
}
// CHI2/LL (OTHERS?)
if (opt.find("LOG") != std::string::npos) {
fIsChi2 = false;
NUIS_ERR(FTL, "No other LIKELIHOODS properly supported!");
NUIS_ERR(FTL, "Try to use a chi2!");
throw;
} else {
fIsChi2 = true;
}
// EXTRAS
if (opt.find("RAW") != std::string::npos)
fIsRawEvents = true;
if (opt.find("DIF") != std::string::npos)
fIsDifXSec = true;
if (opt.find("ENU1D") != std::string::npos)
fIsEnu1D = true;
if (opt.find("NORM") != std::string::npos)
fAddNormPen = true;
if (opt.find("MASK") != std::string::npos)
fIsMask = true;
return;
};
//********************************************************************
void JointMeas1D::SetSmearingMatrix(std::string smearfile, int truedim,
int recodim) {
//********************************************************************
// The smearing matrix describes the migration from true bins (rows) to reco
// bins (columns)
// Counter over the true bins!
int row = 0;
std::string line;
std::ifstream smear(smearfile.c_str(), std::ifstream::in);
// Note that the smearing matrix may be rectangular.
fSmearMatrix = new TMatrixD(truedim, recodim);
if (smear.is_open()) {
NUIS_LOG(SAM, "Reading smearing matrix from file: " << smearfile);
} else {
NUIS_ABORT("Smearing matrix provided is incorrect: " << smearfile);
}
while (std::getline(smear >> std::ws, line, '\n')) {
int column = 0;
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*fSmearMatrix)(row, column) = (*iter) / 100.; // Convert to fraction from
// percentage (this may not be
// general enough)
column++;
}
row++;
}
return;
}
//********************************************************************
void JointMeas1D::ApplySmearingMatrix() {
//********************************************************************
if (!fSmearMatrix) {
NUIS_ERR(WRN,
fName << ": attempted to apply smearing matrix, but none was set");
return;
}
TH1D *unsmeared = (TH1D *)fMCHist->Clone();
TH1D *smeared = (TH1D *)fMCHist->Clone();
smeared->Reset();
// Loop over reconstructed bins
// true = row; reco = column
for (int rbin = 0; rbin < fSmearMatrix->GetNcols(); ++rbin) {
// Sum up the constributions from all true bins
double rBinVal = 0;
// Loop over true bins
for (int tbin = 0; tbin < fSmearMatrix->GetNrows(); ++tbin) {
rBinVal +=
(*fSmearMatrix)(tbin, rbin) * unsmeared->GetBinContent(tbin + 1);
}
smeared->SetBinContent(rbin + 1, rBinVal);
}
fMCHist = (TH1D *)smeared->Clone();
return;
}
/*
Reconfigure LOOP
*/
//********************************************************************
void JointMeas1D::ResetAll() {
//********************************************************************
fMCHist->Reset();
fMCFine->Reset();
fMCStat->Reset();
return;
};
//********************************************************************
void JointMeas1D::FillHistograms() {
//********************************************************************
if (Signal) {
fMCHist->Fill(fXVar, Weight);
fMCFine->Fill(fXVar, Weight);
fMCStat->Fill(fXVar, 1.0);
if (fMCHist_Modes)
fMCHist_Modes->Fill(Mode, fXVar, Weight);
}
return;
};
//********************************************************************
void JointMeas1D::ScaleEvents() {
//********************************************************************
NUIS_LOG(FIT, "Scaling JointMeas1D");
// Fill MCWeighted;
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
fMCWeighted->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1));
fMCWeighted->SetBinError(i + 1, fMCHist->GetBinError(i + 1));
}
// Setup Stat ratios for MC and MC Fine
double *statratio = new double[fMCHist->GetNbinsX()];
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
if (fMCHist->GetBinContent(i + 1) != 0) {
statratio[i] =
fMCHist->GetBinError(i + 1) / fMCHist->GetBinContent(i + 1);
} else {
statratio[i] = 0.0;
}
}
double *statratiofine = new double[fMCFine->GetNbinsX()];
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
if (fMCFine->GetBinContent(i + 1) != 0) {
statratiofine[i] =
fMCFine->GetBinError(i + 1) / fMCFine->GetBinContent(i + 1);
} else {
statratiofine[i] = 0.0;
}
}
// Scaling for raw event rates
if (fIsRawEvents) {
double datamcratio = fDataHist->Integral() / fMCHist->Integral();
fMCHist->Scale(datamcratio);
fMCFine->Scale(datamcratio);
if (fMCHist_Modes)
fMCHist_Modes->Scale(datamcratio);
// Scaling for XSec as function of Enu
} else if (fIsEnu1D) {
PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor, fNEvents);
PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram(),
GetEventHistogram(), fScaleFactor, fNEvents);
// if (fMCHist_Modes) {
// PlotUtils::FluxUnfoldedScaling(fMCHist_Modes, GetFluxHistogram(),
// GetEventHistogram(), fScaleFactor,
// fNEvents);
// }
// Any other differential scaling
} else {
fMCHist->Scale(fScaleFactor, "width");
fMCFine->Scale(fScaleFactor, "width");
if (fMCHist_Modes) {
fMCHist_Modes->Scale(fScaleFactor, "width");
}
}
// Proper error scaling - ROOT Freaks out with xsec weights sometimes
for (int i = 0; i < fMCStat->GetNbinsX(); i++) {
fMCHist->SetBinError(i + 1, fMCHist->GetBinContent(i + 1) * statratio[i]);
}
for (int i = 0; i < fMCFine->GetNbinsX(); i++) {
fMCFine->SetBinError(i + 1,
fMCFine->GetBinContent(i + 1) * statratiofine[i]);
}
// Clean up
delete statratio;
delete statratiofine;
return;
};
//********************************************************************
void JointMeas1D::ApplyNormScale(double norm) {
//********************************************************************
fCurrentNorm = norm;
fMCHist->Scale(1.0 / norm);
fMCFine->Scale(1.0 / norm);
return;
};
/*
Statistic Functions - Outsources to StatUtils
*/
//********************************************************************
int JointMeas1D::GetNDOF() {
//********************************************************************
int ndof = fDataHist->GetNbinsX();
if (fMaskHist)
ndof -= fMaskHist->Integral();
return ndof;
}
//********************************************************************
double JointMeas1D::GetLikelihood() {
//********************************************************************
// If this is for a ratio, there is no data histogram to compare to!
if (fNoData || !fDataHist)
return 0.;
// Apply Masking to MC if Required.
if (fIsMask and fMaskHist) {
PlotUtils::MaskBins(fMCHist, fMaskHist);
}
// Sort Shape Scaling
double scaleF = 0.0;
if (fIsShape) {
if (fMCHist->Integral(1, fMCHist->GetNbinsX(), "width")) {
scaleF = fDataHist->Integral(1, fDataHist->GetNbinsX(), "width") /
fMCHist->Integral(1, fMCHist->GetNbinsX(), "width");
fMCHist->Scale(scaleF);
fMCFine->Scale(scaleF);
}
}
// Likelihood Calculation
double stat = 0.;
if (fIsChi2) {
if (fIsRawEvents) {
stat = StatUtils::GetChi2FromEventRate(fDataHist, fMCHist, fMaskHist);
} else if (fIsDiag) {
stat = StatUtils::GetChi2FromDiag(fDataHist, fMCHist, fMaskHist);
} else if (!fIsDiag and !fIsRawEvents) {
stat = StatUtils::GetChi2FromCov(fDataHist, fMCHist, covar, fMaskHist);
}
}
// Sort Penalty Terms
if (fAddNormPen) {
double penalty =
(1. - fCurrentNorm) * (1. - fCurrentNorm) / (fNormError * fNormError);
stat += penalty;
}
// Return to normal scaling
if (fIsShape and !FitPar::Config().GetParB("saveshapescaling")) {
fMCHist->Scale(1. / scaleF);
fMCFine->Scale(1. / scaleF);
}
fLikelihood = stat;
return stat;
}
/*
Fake Data Functions
*/
//********************************************************************
void JointMeas1D::SetFakeDataValues(std::string fakeOption) {
//********************************************************************
// Setup original/datatrue
TH1D *tempdata = (TH1D *)fDataHist->Clone();
if (!fIsFakeData) {
fIsFakeData = true;
// Make a copy of the original data histogram.
if (!fDataOrig)
fDataOrig = (TH1D *)fDataHist->Clone((fName + "_data_original").c_str());
} else {
ResetFakeData();
}
// Setup Inputs
fFakeDataInput = fakeOption;
NUIS_LOG(SAM, "Setting fake data from : " << fFakeDataInput);
// From MC
if (fFakeDataInput.compare("MC") == 0) {
fDataHist = (TH1D *)fMCHist->Clone((fName + "_MC").c_str());
// Fake File
} else {
if (!fFakeDataFile)
fFakeDataFile = new TFile(fFakeDataInput.c_str(), "READ");
fDataHist = (TH1D *)fFakeDataFile->Get((fName + "_MC").c_str());
}
// Setup Data Hist
fDataHist->SetNameTitle((fName + "_FAKE").c_str(),
(fName + fPlotTitles).c_str());
// Replace Data True
if (fDataTrue)
delete fDataTrue;
fDataTrue = (TH1D *)fDataHist->Clone();
fDataTrue->SetNameTitle((fName + "_FAKE_TRUE").c_str(),
(fName + fPlotTitles).c_str());
// Make a new covariance for fake data hist.
int nbins = fDataHist->GetNbinsX();
double alpha_i = 0.0;
double alpha_j = 0.0;
for (int i = 0; i < nbins; i++) {
for (int j = 0; j < nbins; j++) {
alpha_i =
fDataHist->GetBinContent(i + 1) / tempdata->GetBinContent(i + 1);
alpha_j =
fDataHist->GetBinContent(j + 1) / tempdata->GetBinContent(j + 1);
(*fFullCovar)(i, j) = alpha_i * alpha_j * (*fFullCovar)(i, j);
}
}
// Setup Covariances
if (covar)
delete covar;
covar = StatUtils::GetInvert(fFullCovar);
if (fDecomp)
delete fDecomp;
fDecomp = StatUtils::GetInvert(fFullCovar);
delete tempdata;
return;
};
//********************************************************************
void JointMeas1D::ResetFakeData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist)
delete fDataHist;
fDataHist =
(TH1D *)fDataTrue->Clone((fSettings.GetName() + "_FKDAT").c_str());
}
}
//********************************************************************
void JointMeas1D::ResetData() {
//********************************************************************
if (fIsFakeData) {
if (fDataHist)
delete fDataHist;
fDataHist =
(TH1D *)fDataOrig->Clone((fSettings.GetName() + "_data").c_str());
}
fIsFakeData = false;
}
//********************************************************************
void JointMeas1D::ThrowCovariance() {
//********************************************************************
// Take a fDecomposition and use it to throw the current dataset.
// Requires fDataTrue also be set incase used repeatedly.
if (fDataHist)
delete fDataHist;
fDataHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
return;
};
//********************************************************************
void JointMeas1D::ThrowDataToy() {
//********************************************************************
if (!fDataTrue)
fDataTrue = (TH1D *)fDataHist->Clone();
if (fMCHist)
delete fMCHist;
fMCHist = StatUtils::ThrowHistogram(fDataTrue, fFullCovar);
}
/*
Access Functions
*/
//********************************************************************
TH1D *JointMeas1D::GetMCHistogram() {
//********************************************************************
if (!fMCHist)
return fMCHist;
std::ostringstream chi2;
chi2 << "#chi^{2}=" << std::setprecision(5) << this->GetLikelihood();
int linecolor = kRed;
int linestyle = 1;
int linewidth = 1;
int fillcolor = 0;
int fillstyle = 1001;
if (fSettings.Has("linecolor"))
linecolor = fSettings.GetI("linecolor");
if (fSettings.Has("linestyle"))
linestyle = fSettings.GetI("linestyle");
if (fSettings.Has("linewidth"))
linewidth = fSettings.GetI("linewidth");
if (fSettings.Has("fillcolor"))
fillcolor = fSettings.GetI("fillcolor");
if (fSettings.Has("fillstyle"))
fillstyle = fSettings.GetI("fillstyle");
fMCHist->SetTitle(chi2.str().c_str());
fMCHist->SetLineColor(linecolor);
fMCHist->SetLineStyle(linestyle);
fMCHist->SetLineWidth(linewidth);
fMCHist->SetFillColor(fillcolor);
fMCHist->SetFillStyle(fillstyle);
return fMCHist;
};
//********************************************************************
TH1D *JointMeas1D::GetDataHistogram() {
//********************************************************************
if (!fDataHist)
return fDataHist;
int datacolor = kBlack;
int datastyle = 1;
int datawidth = 1;
if (fSettings.Has("datacolor"))
datacolor = fSettings.GetI("datacolor");
if (fSettings.Has("datastyle"))
datastyle = fSettings.GetI("datastyle");
if (fSettings.Has("datawidth"))
datawidth = fSettings.GetI("datawidth");
fDataHist->SetLineColor(datacolor);
fDataHist->SetLineWidth(datawidth);
fDataHist->SetMarkerStyle(datastyle);
return fDataHist;
};
/*
Write Functions
*/
// Save all the histograms at once
//********************************************************************
void JointMeas1D::Write(std::string drawOpt) {
//********************************************************************
// Get Draw Options
drawOpt = FitPar::Config().GetParS("drawopts");
// Write Settigns
if (drawOpt.find("SETTINGS") != std::string::npos) {
fSettings.Set("#chi^{2}", fLikelihood);
fSettings.Set("NDOF", this->GetNDOF());
fSettings.Set("#chi^{2}/NDOF", fLikelihood / this->GetNDOF());
fSettings.Write();
}
// Write Data/MC
GetDataHistogram()->Write();
GetMCHistogram()->Write();
// Write Fine Histogram
if (drawOpt.find("FINE") != std::string::npos)
GetFineList().at(0)->Write();
// Write Weighted Histogram
if (drawOpt.find("WEIGHTS") != std::string::npos && fMCWeighted)
fMCWeighted->Write();
// Save Flux/Evt if no event manager
if (!FitPar::Config().GetParB("EventManager")) {
if (drawOpt.find("FLUX") != std::string::npos && GetFluxHistogram())
GetFluxHistogram()->Write();
if (drawOpt.find("EVT") != std::string::npos && GetEventHistogram())
GetEventHistogram()->Write();
if (drawOpt.find("XSEC") != std::string::npos && GetEventHistogram())
GetEventHistogram()->Write();
}
// Write Mask
if (fIsMask && (drawOpt.find("MASK") != std::string::npos)) {
fMaskHist->Write();
}
// Write Covariances
if (drawOpt.find("COV") != std::string::npos && fFullCovar) {
- PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName());
+ PlotUtils::GetFullCovarPlot(fFullCovar, fSettings.GetName())->Write();
}
if (drawOpt.find("INVCOV") != std::string::npos && covar) {
- PlotUtils::GetInvCovarPlot(covar, fSettings.GetName());
+ PlotUtils::GetInvCovarPlot(covar, fSettings.GetName())->Write();
}
if (drawOpt.find("DECOMP") != std::string::npos && fDecomp) {
- PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName());
+ PlotUtils::GetDecompCovarPlot(fDecomp, fSettings.GetName())->Write();
}
// // Likelihood residual plots
// if (drawOpt.find("RESIDUAL") != std::string::npos) {
// WriteResidualPlots();
// }
// Ratio and Shape Plots
if (drawOpt.find("RATIO") != std::string::npos) {
WriteRatioPlot();
}
if (drawOpt.find("SHAPE") != std::string::npos) {
WriteShapePlot();
if (drawOpt.find("RATIO") != std::string::npos)
WriteShapeRatioPlot();
}
// // RATIO
// if (drawOpt.find("CANVMC") != std::string::npos) {
// TCanvas* c1 = WriteMCCanvas(fDataHist, fMCHist);
// c1->Write();
// delete c1;
// }
// // PDG
// if (drawOpt.find("CANVPDG") != std::string::npos && fMCHist_Modes) {
// TCanvas* c2 = WritePDGCanvas(fDataHist, fMCHist, fMCHist_Modes);
// c2->Write();
// delete c2;
// }
// Write Extra Histograms
AutoWriteExtraTH1();
WriteExtraHistograms();
if (fSaveSubMeas) {
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
exp->Write(drawOpt);
}
}
// Returning
NUIS_LOG(SAM, "Written Histograms: " << fName);
return;
}
//********************************************************************
void JointMeas1D::WriteRatioPlot() {
//********************************************************************
// Setup mc data ratios
TH1D *dataRatio = (TH1D *)fDataHist->Clone((fName + "_data_RATIO").c_str());
TH1D *mcRatio = (TH1D *)fMCHist->Clone((fName + "_MC_RATIO").c_str());
// Extra MC Data Ratios
for (int i = 0; i < mcRatio->GetNbinsX(); i++) {
dataRatio->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) /
fMCHist->GetBinContent(i + 1));
dataRatio->SetBinError(i + 1, fDataHist->GetBinError(i + 1) /
fMCHist->GetBinContent(i + 1));
mcRatio->SetBinContent(i + 1, fMCHist->GetBinContent(i + 1) /
fMCHist->GetBinContent(i + 1));
mcRatio->SetBinError(i + 1, fMCHist->GetBinError(i + 1) /
fMCHist->GetBinContent(i + 1));
}
// Write ratios
mcRatio->Write();
dataRatio->Write();
delete mcRatio;
delete dataRatio;
}
//********************************************************************
void JointMeas1D::WriteShapePlot() {
//********************************************************************
TH1D *mcShape = (TH1D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
std::stringstream ss;
ss << "Scale=" << shapeScale;
mcShape->SetTitle(ss.str().c_str());
mcShape->SetLineWidth(3);
mcShape->SetLineStyle(7);
mcShape->Write();
delete mcShape;
}
//********************************************************************
void JointMeas1D::WriteShapeRatioPlot() {
//********************************************************************
// Get a mcshape histogram
TH1D *mcShape = (TH1D *)fMCHist->Clone((fName + "_MC_SHAPE").c_str());
double shapeScale = 1.0;
if (fIsRawEvents) {
shapeScale = fDataHist->Integral() / fMCHist->Integral();
} else {
shapeScale = fDataHist->Integral("width") / fMCHist->Integral("width");
}
mcShape->Scale(shapeScale);
// Create shape ratio histograms
TH1D *mcShapeRatio =
(TH1D *)mcShape->Clone((fName + "_MC_SHAPE_RATIO").c_str());
TH1D *dataShapeRatio =
(TH1D *)fDataHist->Clone((fName + "_data_SHAPE_RATIO").c_str());
// Divide the histograms
mcShapeRatio->Divide(mcShape);
dataShapeRatio->Divide(mcShape);
// Colour the shape ratio plots
mcShapeRatio->SetLineWidth(3);
mcShapeRatio->SetLineStyle(7);
mcShapeRatio->Write();
dataShapeRatio->Write();
delete mcShapeRatio;
delete dataShapeRatio;
}
/*
Setup Functions
*/
//********************************************************************
void JointMeas1D::SetupMeasurement(std::string input, std::string type,
FitWeight *rw, std::string fkdt) {
//********************************************************************
// For joint samples, input files are given as a semi-colon seperated list.
// Parse this list and save it for later, and set up the types etc.
if (FitPar::Config().GetParB("EventManager")) {
NUIS_ERR(FTL, "Event Manager does not yet work with JointMeas1D Samples");
NUIS_ERR(FTL, "If you want good predictions for "
<< fName
<< " then run with it turned off! (-q EventManager=0)");
}
fSubInFiles.clear();
std::vector entries = GeneralUtils::ParseToStr(input, ";");
if (entries.size() < 2) {
NUIS_ABORT("Joint measurement expected to recieve at least two semi-colon "
"separated input files, but recieved: \""
<< input << "\"");
}
std::vector first_file_descriptor =
GeneralUtils::ParseToStr(entries.front(), ":");
if (first_file_descriptor.size() != 2) {
NUIS_ABORT("Found Joint measurement where the input file had no type: \""
<< input << "\", expected \"INPUTTYPE:File.root;File2.root\".");
}
std::string inpType = first_file_descriptor[0];
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
if (GeneralUtils::ParseToStr(*iter, ":").size() != 2) {
std::stringstream ss("");
ss << inpType << ":" << (*iter);
fSubInFiles.push_back(ss.str());
} else {
fSubInFiles.push_back(*iter);
}
}
// Set Engine and Fake Data
fRW = rw;
fFakeDataInput = fkdt;
// Set Fit Options
SetFitOptions(type);
return;
}
/*
XSec Functions
*/
//********************************************************************
double JointMeas1D::TotalIntegratedFlux(std::string intOpt, double low,
double high) {
//********************************************************************
double totalflux = 0.0;
// Destroy the job for sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
double expflux = exp->TotalIntegratedFlux(intOpt, low, high);
// Fill flux options
if (fIsRatio) {
totalflux = expflux;
break;
}
if (fIsSummed) {
totalflux += expflux;
}
}
return totalflux;
}
/*
Reconfigure Functions
*/
//********************************************************************
void JointMeas1D::Reconfigure() {
//********************************************************************
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
exp->Reconfigure();
}
ConvertEventRates();
return;
}
//********************************************************************
void JointMeas1D::ConvertEventRates() {
//********************************************************************
// Apply Event Scaling
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = static_cast(*expIter);
exp->ScaleEvents();
}
// Joint function called by top level class
MakePlots();
// Do Final Normalisation
ApplyNormScale(fRW->GetSampleNorm(this->fName));
}
//********************************************************************
void JointMeas1D::MakePlots() {
//********************************************************************
// Reset the 1D histograms but not the subClasses
ResetAll();
// If Summed
if (fIsSummed) {
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = static_cast(*expIter);
this->fMCHist->Add(exp->GetMCList().at(0));
this->fMCFine->Add(exp->GetFineList().at(0));
}
return;
}
// If Ratio
if (fIsRatio) {
int sample = 0;
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
if (sample == 0) {
this->fMCHist->Add(exp->GetMCList().at(0));
this->fMCFine->Add(exp->GetFineList().at(0));
} else if (sample == 1) {
this->fMCHist->Divide(exp->GetMCList().at(0));
this->fMCFine->Divide(exp->GetFineList().at(0));
} else {
break;
}
sample++;
}
return;
}
return;
}
/*
Access Functions
*/
//********************************************************************
std::vector JointMeas1D::GetMCList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(this->fMCHist);
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
std::vector subTempVect = exp->GetMCList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetDataList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(this->fDataHist);
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
std::vector subTempVect = exp->GetDataList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetFineList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(this->fMCFine);
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
std::vector subTempVect = exp->GetFineList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetMaskList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(this->fMaskHist);
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
std::vector subTempVect = exp->GetMaskList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetFluxList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(MeasurementBase::GetFluxHistogram());
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
std::vector subTempVect = exp->GetFluxList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetEventRateList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(MeasurementBase::GetEventHistogram());
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
std::vector subTempVect = exp->GetEventRateList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
std::vector JointMeas1D::GetXSecList() {
//********************************************************************
// Make Default Vector
std::vector tempVect;
tempVect.push_back(MeasurementBase::GetXSecHistogram());
// Return vector from all sub samples
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
std::vector subTempVect = exp->GetXSecList();
for (UInt_t i = 0; i < subTempVect.size(); i++) {
tempVect.push_back(subTempVect.at(i));
}
}
return tempVect;
}
//********************************************************************
TH1D *JointMeas1D::GetCombinedFlux() {
//********************************************************************
TH1D *newflux = NULL;
int sample = 0;
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
// Get flux from experiment
std::vector fluxVect = exp->GetFluxList();
// Setup newflux
if (sample == 0) {
newflux = (TH1D *)fluxVect.at(0);
newflux->Reset();
}
// Add all fluxes
for (UInt_t i = 0; i < fluxVect.size(); i++) {
newflux->Add((TH1D *)fluxVect.at(i));
sample++;
}
}
if (!newflux) {
NUIS_ABORT("No combined flux setup in JointMeas1D");
}
return newflux;
}
//********************************************************************
TH1D *JointMeas1D::GetCombinedEventRate() {
//********************************************************************
TH1D *newflux = NULL;
int sample = 0;
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
MeasurementBase *exp = *expIter;
// Get flux from experiment
std::vector fluxVect = exp->GetFluxList();
// Setup newflux
if (sample == 0) {
newflux = (TH1D *)fluxVect.at(0);
newflux->Reset();
}
// Add all fluxes
for (UInt_t i = 0; i < fluxVect.size(); i++) {
newflux->Add(fluxVect.at(i));
sample++;
}
}
if (!newflux) {
NUIS_ABORT("No combined event rate setup in JointMeas1D");
}
return newflux;
}
//********************************************************************
std::vector JointMeas1D::GetSubSamples() {
//********************************************************************
std::vector exps;
for (std::vector::const_iterator expIter =
fSubChain.begin();
expIter != fSubChain.end(); expIter++) {
exps.push_back(*expIter);
}
return exps;
}
//// CRAP TO BE REMOVED
//********************************************************************
void JointMeas1D::SetDataValues(std::string dataFile) {
//********************************************************************
// Override this function if the input file isn't in a suitable format
NUIS_LOG(SAM, "Reading data from: " << dataFile.c_str());
fDataHist =
PlotUtils::GetTH1DFromFile(dataFile, (fName + "_data"), fPlotTitles);
fDataTrue = (TH1D *)fDataHist->Clone();
// Number of data points is number of bins
fNDataPointsX = fDataHist->GetXaxis()->GetNbins();
return;
};
//********************************************************************
void JointMeas1D::SetDataFromDatabase(std::string inhistfile,
std::string histname) {
//********************************************************************
NUIS_LOG(SAM, "Filling histogram from " << inhistfile << "->" << histname);
fDataHist = PlotUtils::GetTH1DFromRootFile(
(GeneralUtils::GetTopLevelDir() + "/data/" + inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void JointMeas1D::SetDataFromFile(std::string inhistfile,
std::string histname) {
//********************************************************************
NUIS_LOG(SAM, "Filling histogram from " << inhistfile << "->" << histname);
fDataHist = PlotUtils::GetTH1DFromRootFile((inhistfile), histname);
fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_data").c_str());
return;
};
//********************************************************************
void JointMeas1D::SetCovarMatrix(std::string covarFile) {
//********************************************************************
// Covariance function, only really used when reading in the MB Covariances.
TFile *tempFile = new TFile(covarFile.c_str(), "READ");
TH2D *covarPlot = new TH2D();
TH2D *fFullCovarPlot = new TH2D();
std::string covName = "";
std::string covOption = FitPar::Config().GetParS("thrown_covariance");
if (fIsShape || fIsFree)
covName = "shp_";
if (fIsDiag)
covName += "diag";
else
covName += "full";
covarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str());
if (!covOption.compare("SUB"))
fFullCovarPlot = (TH2D *)tempFile->Get((covName + "cov").c_str());
else if (!covOption.compare("FULL"))
fFullCovarPlot = (TH2D *)tempFile->Get("fullcov");
else {
NUIS_ERR(WRN, "Incorrect thrown_covariance option in parameters.");
}
int dim = int(fDataHist->GetNbinsX()); //-this->masked->Integral());
int covdim = int(fDataHist->GetNbinsX());
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
fDecomp = new TMatrixDSym(dim);
int row, column = 0;
row = 0;
column = 0;
for (Int_t i = 0; i < covdim; i++) {
// if (this->masked->GetBinContent(i+1) > 0) continue;
for (Int_t j = 0; j < covdim; j++) {
// if (this->masked->GetBinContent(j+1) > 0) continue;
(*this->covar)(row, column) = covarPlot->GetBinContent(i + 1, j + 1);
(*fFullCovar)(row, column) = fFullCovarPlot->GetBinContent(i + 1, j + 1);
column++;
}
column = 0;
row++;
}
// Set bin errors on data
if (!fIsDiag) {
StatUtils::SetDataErrorFromCov(fDataHist, fFullCovar);
}
// Get Deteriminant and inverse matrix
// fCovDet = this->covar->Determinant();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
// Sets the covariance matrix from a provided file in a text format
// scale is a multiplicative pre-factor to apply in the case where the
// covariance is given in some unit (e.g. 1E-38)
void JointMeas1D::SetCovarMatrixFromText(std::string covarFile, int dim,
double scale) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
std::ifstream covarread(covarFile.c_str(), std::ifstream::in);
this->covar = new TMatrixDSym(dim);
fFullCovar = new TMatrixDSym(dim);
if (covarread.is_open()) {
NUIS_LOG(SAM, "Reading covariance matrix from file: " << covarFile);
} else {
NUIS_ERR(FTL, "Covariance matrix provided is incorrect: " << covarFile);
}
// Loop over the lines in the file
while (std::getline(covarread >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
if (entries.size() <= 1) {
NUIS_ERR(WRN, "SetCovarMatrixFromText -> Covariance matrix only has <= 1 "
"entries on this line: "
<< row);
}
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
(*covar)(row, column) = *iter;
(*fFullCovar)(row, column) = *iter;
column++;
}
row++;
}
covarread.close();
// Scale the actualy covariance matrix by some multiplicative factor
(*fFullCovar) *= scale;
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
// THIS IS ACTUALLY THE INVERSE COVARIANCE MATRIXA AAAAARGH
delete this->covar;
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
// Now need to multiply by the scaling factor
// If the covariance
(*this->covar) *= 1. / (scale);
return;
};
//********************************************************************
void JointMeas1D::SetCovarMatrixFromCorrText(std::string corrFile, int dim) {
//********************************************************************
// Make a counter to track the line number
int row = 0;
std::string line;
std::ifstream corr(corrFile.c_str(), std::ifstream::in);
this->covar = new TMatrixDSym(dim);
this->fFullCovar = new TMatrixDSym(dim);
if (corr.is_open()) {
NUIS_LOG(SAM,
"Reading and converting correlation matrix from file: " << corrFile);
} else {
NUIS_ERR(FTL, "Correlation matrix provided is incorrect: " << corrFile);
exit(-1);
}
while (std::getline(corr >> std::ws, line, '\n')) {
int column = 0;
// Loop over entries and insert them into matrix
// Multiply by the errors to get the covariance, rather than the correlation
// matrix
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
for (std::vector::iterator iter = entries.begin();
iter != entries.end(); iter++) {
double val = (*iter) * this->fDataHist->GetBinError(row + 1) * 1E38 *
this->fDataHist->GetBinError(column + 1) * 1E38;
if (val == 0) {
NUIS_ERR(FTL, "Found a zero value in the covariance matrix, assuming "
"this is an error!");
exit(-1);
}
(*this->covar)(row, column) = val;
(*this->fFullCovar)(row, column) = val;
column++;
}
row++;
}
// Robust matrix inversion method
TDecompSVD LU = TDecompSVD(*this->covar);
delete this->covar;
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
return;
};
//********************************************************************
// FullUnits refers to if we have "real" unscaled units in the covariance
// matrix, e.g. 1E-76. If this is the case we need to scale it so that the chi2
// contribution is correct NUISANCE internally assumes the covariance matrix has
// units of 1E76
void JointMeas1D::SetCovarFromDataFile(std::string covarFile,
std::string covName, bool FullUnits) {
//********************************************************************
NUIS_LOG(SAM, "Getting covariance from " << covarFile << "->" << covName);
TFile *tempFile = new TFile(covarFile.c_str(), "READ");
TH2D *covPlot = (TH2D *)tempFile->Get(covName.c_str());
covPlot->SetDirectory(0);
// Scale the covariance matrix if it comes in normal units
if (FullUnits) {
covPlot->Scale(1.E76);
}
int dim = covPlot->GetNbinsX();
fFullCovar = new TMatrixDSym(dim);
for (int i = 0; i < dim; i++) {
for (int j = 0; j < dim; j++) {
(*fFullCovar)(i, j) = covPlot->GetBinContent(i + 1, j + 1);
}
}
this->covar = (TMatrixDSym *)fFullCovar->Clone();
fDecomp = (TMatrixDSym *)fFullCovar->Clone();
TDecompSVD LU = TDecompSVD(*this->covar);
this->covar = new TMatrixDSym(dim, LU.Invert().GetMatrixArray(), "");
TDecompChol LUChol = TDecompChol(*fDecomp);
LUChol.Decompose();
fDecomp = new TMatrixDSym(dim, LU.GetU().GetMatrixArray(), "");
return;
};
// std::vector JointMeas1D::GetMCList(void){
// std::vector temp;
// return temp;
// }
// std::vector JointMeas1D::GetDataList(void){
// std::vector temp;
// return temp;
// }
// std::vector JointMeas1D::GetMaskList(void){
// std::vector temp;
// return temp;
// }
// std::vector JointMeas1D::GetFineList(void){
// std::vector temp;
// return temp;
// }
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 24dad08..cb1b9fe 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,599 +1,603 @@
// 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 .
*******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#ifdef GENIE_PRE_R3
#include "Messenger/Messenger.h"
#else
#include "Framework/Messenger/Messenger.h"
#endif
#include "InputUtils.h"
GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
void GENIEGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
}
void GENIEGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
}
void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
fGenieParticlePDGs = new int[stacksize];
}
void GENIEGeneratorInfo::DeallocateParticleStack() {
delete fGenieParticlePDGs;
}
void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord *ntpl) {
Reset();
// Check for GENIE Event
if (!ntpl)
return;
if (!ntpl->event)
return;
// Cast Event Record
GHepRecord *ghep = static_cast(ntpl->event);
if (!ghep)
return;
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(ghep);
// Loop over all particles
int i = 0;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p)
continue;
// Get PDG
fGenieParticlePDGs[i] = p->Pdg();
i++;
}
}
void GENIEGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fGenieParticlePDGs[i] = 0;
}
}
GENIEInputHandler::GENIEInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating GENIEInputHandler : " << handle);
+ // Plz no shouting
+ StopTalking();
genie::Messenger::Instance()->SetPriorityLevel("GHepUtils", pFATAL);
+ StartTalking();
+ // Shout all you want
// Run a joint input handling
fName = handle;
// Setup the TChain
fGENIETree = new TChain("gtree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Are we running with NOvA weights
fNOvAWeights = FitPar::Config().GetParB("NOvA_Weights");
MAQEw = 1.0;
NonResw = 1.0;
RPAQEw = 1.0;
RPARESw = 1.0;
MECw = 1.0;
DISw = 1.0;
NOVAw = 1.0;
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile *inp_file = new TFile(
InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
NUIS_ABORT(
"GENIE File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D *fluxhist = (TH1D *)inp_file->Get("nuisance_flux");
TH1D *eventhist = (TH1D *)inp_file->Get("nuisance_events");
if (!fluxhist or !eventhist) {
NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
NUIS_ABORT("GENIE FILE doesn't contain flux/xsec info."
<< std::endl
<< "Try running the app PrepareGENIE first on :"
<< inputs[inp_it] << std::endl
<< "$ PrepareGENIE -h");
}
// Get N Events
TTree *genietree = (TTree *)inp_file->Get("gtree");
if (!genietree) {
NUIS_ERR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
NUIS_ABORT(
"Check your inputs, they may need to be completely regenerated!");
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
NUIS_ABORT("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Check for precomputed weights
TTree *weighttree = (TTree *)inp_file->Get("nova_wgts");
if (fNOvAWeights) {
if (!weighttree) {
NUIS_ABORT("Did not find nova_wgts tree in file "
<< inputs[inp_it] << " but you specified it" << std::endl);
} else {
NUIS_LOG(FIT, "Found nova_wgts tree in file " << inputs[inp_it]);
}
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fGENIETree->AddFile(inputs[inp_it].c_str());
if (weighttree != NULL)
fGENIETree->AddFriend(weighttree);
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
// Set up the custom weights
if (fNOvAWeights) {
fGENIETree->SetBranchAddress("MAQEwgt", &MAQEw);
fGENIETree->SetBranchAddress("nonResNormWgt", &NonResw);
fGENIETree->SetBranchAddress("RPAQEWgt", &RPAQEw);
fGENIETree->SetBranchAddress("RPARESWgt", &RPARESw);
fGENIETree->SetBranchAddress("MECWgt", &MECw);
fGENIETree->SetBranchAddress("DISWgt", &DISw);
fGENIETree->SetBranchAddress("nova2018CVWgt", &NOVAw);
}
// Libraries should be seen but not heard...
StopTalking();
fGENIETree->GetEntry(0);
StartTalking();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
};
GENIEInputHandler::~GENIEInputHandler() {
// if (fGenieGHep) delete fGenieGHep;
// if (fGenieNtpl) delete fGenieNtpl;
// if (fGENIETree) delete fGENIETree;
// if (fGenieInfo) delete fGenieInfo;
}
void GENIEInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 1);
fGENIETree->SetCacheSize(fCacheSize);
}
}
void GENIEInputHandler::RemoveCache() {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 0);
fGENIETree->SetCacheSize(0);
}
FitEvent *GENIEInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
UInt_t entry = ent + fSkip;
if (entry >= (UInt_t)fNEvents)
return NULL;
// Clear the previous event (See Note 1 in ROOT TClonesArray documentation)
if (fGenieNtpl) {
fGenieNtpl->Clear();
}
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGENIETree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
// Check for GENIE Event
if (!fGenieNtpl)
return NULL;
if (!fGenieNtpl->event)
return NULL;
// Cast Event Record
fGenieGHep = static_cast(fGenieNtpl->event);
if (!fGenieGHep)
return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle *p;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
if (state != genie::kIStInitialState) {
continue;
}
fNUISANCEEvent->probe_E = p->E() * 1.E3;
fNUISANCEEvent->probe_pdg = p->Pdg();
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle *p,
int mode) {
/*
kIStUndefined = -1,
kIStInitialState = 0, / generator-level initial state /
kIStStableFinalState = 1, / generator-level final state:
particles to be tracked by detector-level MC /
kIStIntermediateState = 2,
kIStDecayedState = 3,
kIStCorrelatedNucleon = 10,
kIStNucleonTarget = 11,
kIStDISPreFragmHadronicState = 12,
kIStPreDecayResonantState = 13,
kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked
for hadron transport modules to act on /
kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments
entering the record collectively as a 'hadronic blob' pseudo-particle /
kIStNucleonClusterTarget = 16, // for composite nucleons before
phase space decay
*/
int state = kUndefinedState;
switch (p->Status()) {
case genie::kIStNucleonTarget:
case genie::kIStInitialState:
case genie::kIStCorrelatedNucleon:
case genie::kIStNucleonClusterTarget:
state = kInitialState;
break;
case genie::kIStStableFinalState:
state = kFinalState;
break;
case genie::kIStHadronInTheNucleus:
if (abs(mode) == 2)
state = kInitialState;
else
state = kFSIState;
break;
case genie::kIStPreDecayResonantState:
case genie::kIStDISPreFragmHadronicState:
case genie::kIStIntermediateState:
state = kFSIState;
break;
case genie::kIStFinalStateNuclearRemnant:
case genie::kIStUndefined:
case genie::kIStDecayedState:
default:
break;
}
// Flag to remove nuclear part in genie
if (p->Pdg() > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord *gheprec) {
// Electron Scattering
if (gheprec->Summary()->ProcInfo().IsEM()) {
if (gheprec->Summary()->InitState().ProbePdg() == 11) {
if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
return 1;
else if (gheprec->Summary()->ProcInfo().IsMEC())
return 2;
else if (gheprec->Summary()->ProcInfo().IsResonant())
return 13;
else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
return 26;
else {
NUIS_ERR(WRN,
"Unknown GENIE Electron Scattering Mode!"
<< std::endl
<< "ScatteringTypeId = "
<< gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
<< "InteractionTypeId = "
<< gheprec->Summary()->ProcInfo().InteractionTypeId()
<< std::endl
<< genie::ScatteringType::AsString(
gheprec->Summary()->ProcInfo().ScatteringTypeId())
<< " "
<< genie::InteractionType::AsString(
gheprec->Summary()->ProcInfo().InteractionTypeId())
<< " " << gheprec->Summary()->ProcInfo().IsMEC());
return 0;
}
}
// Weak CC
} else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 2;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -2;
// CC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
// Weak NC
} else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 32;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -32;
// NC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
}
return 0;
}
void GENIEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Check for GENIE Event
if (!fGenieNtpl)
return;
if (!fGenieNtpl->event)
return;
// Cast Event Record
fGenieGHep = static_cast(fGenieNtpl->event);
if (!fGenieGHep)
return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
// Have a bool storing if interaction happened on free or bound nucleon
bool IsFree = false;
// Set the TargetPDG
if (fGenieGHep->TargetNucleus() != NULL) {
fNUISANCEEvent->fTargetPDG = fGenieGHep->TargetNucleus()->Pdg();
IsFree = false;
// Sometimes GENIE scatters off free nucleons, electrons, photons
// In which TargetNucleus is NULL and we need to find the initial state
// particle
} else {
// Check the particle is an initial state particle
// Follows GHepRecord::TargetNucleusPosition but doesn't do check on
// pdg::IsIon
GHepParticle *p = fGenieGHep->Particle(1);
// Check that particle 1 actually exists
if (!p) {
NUIS_ABORT("Can't find particle 1 for GHepRecord");
}
// If not an ion but is an initial state particle
if (!pdg::IsIon(p->Pdg()) && p->Status() == kIStInitialState) {
IsFree = true;
fNUISANCEEvent->fTargetPDG = p->Pdg();
// Catch if something strange happens:
// Here particle 1 is not an initial state particle OR
// particle 1 is an ion OR
// both
} else {
if (pdg::IsIon(p->Pdg())) {
NUIS_ABORT(
"Particle 1 in GHepRecord stack is an ion but isn't an initial "
"state particle");
} else {
NUIS_ABORT(
"Particle 1 in GHepRecord stack is not an ion but is an initial "
"state particle");
}
}
}
// Set the A and Z and H from the target PDG
// Depends on if we scattered off a free or bound nucleon
if (!IsFree) {
fNUISANCEEvent->fTargetA =
TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetZ =
TargetUtils::GetTargetZFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetH = 0;
} else {
// If free proton scattering
if (fNUISANCEEvent->fTargetPDG == 2212) {
fNUISANCEEvent->fTargetA = 1;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 1;
// If free neutron scattering
} else if (fNUISANCEEvent->fTargetPDG == 2112) {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 0;
// If neither
} else {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 0;
fNUISANCEEvent->fTargetH = 0;
}
}
fNUISANCEEvent->fBound = !IsFree;
fNUISANCEEvent->InputWeight =
1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
// And the custom weights
if (fNOvAWeights) {
fNUISANCEEvent->CustomWeight = NOVAw;
fNUISANCEEvent->CustomWeightArray[0] = MAQEw;
fNUISANCEEvent->CustomWeightArray[1] = NonResw;
fNUISANCEEvent->CustomWeightArray[2] = RPAQEw;
fNUISANCEEvent->CustomWeightArray[3] = RPARESw;
fNUISANCEEvent->CustomWeightArray[4] = MECw;
fNUISANCEEvent->CustomWeightArray[5] = NOVAw;
} else {
fNUISANCEEvent->CustomWeight = 1.0;
fNUISANCEEvent->CustomWeightArray[0] = 1.0;
fNUISANCEEvent->CustomWeightArray[1] = 1.0;
fNUISANCEEvent->CustomWeightArray[2] = 1.0;
fNUISANCEEvent->CustomWeightArray[3] = 1.0;
fNUISANCEEvent->CustomWeightArray[4] = 1.0;
fNUISANCEEvent->CustomWeightArray[5] = 1.0;
}
// Get N Particle Stack
unsigned int npart = fGenieGHep->GetEntries();
unsigned int kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
NUIS_ERR(WRN, "GENIE has too many particles, expanding stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(fGenieGHep);
fNUISANCEEvent->fNParticles = 0;
// Loop over all particles
while ((p = (dynamic_cast((iter).Next())))) {
if (!p)
continue;
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState)
continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState)
continue;
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// Fill Vectors
int curpart = fNUISANCEEvent->fNParticles;
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
// Set if the particle was on the fundamental vertex
fNUISANCEEvent->fPrimaryVertex[curpart] = (p->FirstMother() < 2);
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
NUIS_ERR(WRN, "Number of GENIE Particles exceeds maximum!");
NUIS_ERR(WRN, "Extend kMax, or run without including FSI particles!");
break;
}
}
// Fill Extra Stack
if (fSaveExtra)
fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle *ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
if (ISAnyLepton) {
fNUISANCEEvent->probe_E = ISAnyLepton->E();
fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/T2K/CMakeLists.txt b/src/T2K/CMakeLists.txt
index 6246dfe..b706840 100644
--- a/src/T2K/CMakeLists.txt
+++ b/src/T2K/CMakeLists.txt
@@ -1,102 +1,106 @@
# 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 .
################################################################################
set(IMPLFILES
T2K_CC0pi_XSec_2DPcos_nu_I.cxx
T2K_CC0pi_XSec_2DPcos_nu_II.cxx
T2K_CCinc_XSec_2DPcos_nu_nonuniform.cxx
+T2K_nueCCinc_XSec_1Dpe.cxx
+T2K_nueCCinc_XSec_joint.cxx
T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.cxx
T2K_CC1pip_CH_XSec_1Dppi_nu.cxx
T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx
T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx
T2K_CC1pip_CH_XSec_1DQ2_nu.cxx
T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.cxx
T2K_CC1pip_CH_XSec_1DCosThAdler_nu.cxx
T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.cxx
T2K_CC1pip_H2O_XSec_1DEnuMB_nu.cxx
T2K_CC1pip_H2O_XSec_1Dcosmu_nu.cxx
T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.cxx
T2K_CC1pip_H2O_XSec_1Dcospi_nu.cxx
T2K_CC1pip_H2O_XSec_1Dpmu_nu.cxx
T2K_CC1pip_H2O_XSec_1Dppi_nu.cxx
T2K_CC0pinp_STV_XSec_1Ddpt_nu.cxx
T2K_CC0pinp_STV_XSec_1Ddphit_nu.cxx
T2K_CC0pinp_STV_XSec_1Ddat_nu.cxx
T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.cxx
T2K_CC0pinp_ifk_XSec_3Dinfp_nu.cxx
T2K_CC0pinp_ifk_XSec_3Dinfa_nu.cxx
T2K_CC0pinp_ifk_XSec_3Dinfip_nu.cxx
T2K_SignalDef.cxx
)
set(HEADERFILES
T2K_CC0pi_XSec_2DPcos_nu_I.h
T2K_CC0pi_XSec_2DPcos_nu_II.h
T2K_CCinc_XSec_2DPcos_nu_nonuniform.h
+T2K_nueCCinc_XSec_1Dpe.h
+T2K_nueCCinc_XSec_joint.h
T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h
T2K_CC1pip_CH_XSec_1Dppi_nu.h
T2K_CC1pip_CH_XSec_1Dthpi_nu.h
T2K_CC1pip_CH_XSec_1Dthmupi_nu.h
T2K_CC1pip_CH_XSec_1DQ2_nu.h
T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h
T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h
T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h
T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h
T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h
T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h
T2K_CC1pip_H2O_XSec_1Dcospi_nu.h
T2K_CC1pip_H2O_XSec_1Dpmu_nu.h
T2K_CC1pip_H2O_XSec_1Dppi_nu.h
T2K_CC0pinp_STV_XSec_1Ddpt_nu.h
T2K_CC0pinp_STV_XSec_1Ddphit_nu.h
T2K_CC0pinp_STV_XSec_1Ddat_nu.h
T2K_CC0pi1p_XSec_3DPcoscos_nu_nonuniform.h
T2K_CC0pinp_ifk_XSec_3Dinfp_nu.h
T2K_CC0pinp_ifk_XSec_3Dinfa_nu.h
T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h
T2K_SignalDef.h
)
set(LIBNAME expT2K)
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
add_library(${LIBNAME} STATIC ${IMPLFILES})
else(CMAKE_BUILD_TYPE MATCHES RELEASE)
add_library(${LIBNAME} SHARED ${IMPLFILES})
endif()
include_directories(${MINIMUM_INCLUDE_DIRECTORIES})
set_target_properties(${LIBNAME} PROPERTIES VERSION
"${NUISANCE_VERSION_MAJOR}.${NUISANCE_VERSION_MINOR}.${NUISANCE_VERSION_REVISION}")
#set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS})
if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES)
add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES})
endif()
install(TARGETS ${LIBNAME} DESTINATION lib)
#Can uncomment this to install the headers... but is it really neccessary?
install(FILES ${HEADERFILES} DESTINATION include/T2K)
set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE)
diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
index 882cc20..b359e64 100644
--- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
+++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
@@ -1,328 +1,328 @@
// 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 "T2K_SignalDef.h"
#include "T2K_CC0pi_XSec_2DPcos_nu_I.h"
static size_t nangbins = 9;
static double angular_binning_costheta[] = {-1, 0, 0.6, 0.7, 0.8,
0.85, 0.9, 0.94, 0.98, 1};
//********************************************************************
T2K_CC0pi_XSec_2DPcos_nu_I::T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu_I sample. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0pi\n"
"https://journals.aps.org/prd/abstract/10.1103/"
"PhysRevD.93.112012 Analysis I";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("P_{#mu} (GeV)");
fSettings.SetYTitle("cos#theta_{#mu}");
fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV",
"FIX/FULL");
fSettings.DefineAllowedTargets("C,H");
fSettings.SetEnuRangeFromFlux(fFluxHist);
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pi_XSec_2DPcos_nu_I");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
fMaskMomOverflow = false;
if (samplekey.Has("mask_mom_overflow")) {
fMaskMomOverflow = samplekey.GetB("mask_mom_overflow");
}
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) *
1E-38 / (TotalIntegratedFlux()));
assert(std::isnormal(fScaleFactor));
// Setup Histograms
SetHistograms();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
bool T2K_CC0pi_XSec_2DPcos_nu_I::isSignal(FitEvent *event) {
return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_I);
};
void T2K_CC0pi_XSec_2DPcos_nu_I::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double pmu = Pmu.Vect().Mag() / 1000.;
double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect()));
fXVar = pmu;
fYVar = CosThetaMu;
return;
};
void T2K_CC0pi_XSec_2DPcos_nu_I::FillHistograms() {
Measurement1D::FillHistograms();
if (Signal) {
fMCHist_Fine2D->Fill(fXVar, fYVar, Weight);
FillMCSlice(fXVar, fYVar, Weight);
}
}
// Modification is needed after the full reconfigure to move bins around
// Otherwise this would need to be replaced by a TH2Poly which is too awkward.
void T2K_CC0pi_XSec_2DPcos_nu_I::ConvertEventRates() {
// Do standard conversion.
Measurement1D::ConvertEventRates();
// Scale MC slices by their bin area
for (size_t i = 0; i < nangbins; ++i) {
fMCHist_Slices[i]->Scale(fScaleFactor / (angular_binning_costheta[i + 1] -
angular_binning_costheta[i]),
"width");
}
// Now Convert into 1D list
fMCHist->Reset();
int bincount = 0;
for (size_t i = 0; i < nangbins; i++) {
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
fMCHist->SetBinContent(bincount + 1,
fMCHist_Slices[i]->GetBinContent(j + 1));
fMCHist->SetBinError(bincount + 1, fMCHist_Slices[i]->GetBinError(j + 1));
bincount++;
}
}
return;
}
void T2K_CC0pi_XSec_2DPcos_nu_I::FillMCSlice(double x, double y, double w) {
for (size_t i = 0; i < nangbins; ++i) {
if ((y >= angular_binning_costheta[i]) &&
(y < angular_binning_costheta[i + 1])) {
fMCHist_Slices[i]->Fill(x, w);
}
}
}
void T2K_CC0pi_XSec_2DPcos_nu_I::SetHistograms() {
// Read in 1D Data Histograms
TFile input(
(FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root")
.c_str(),
"READ");
fMCHist_Fine2D = new TH2D("T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D",
"T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", 400, 0.0, 30.0,
100, -1.0, 1.0);
fMCHist_Fine2D->SetDirectory(NULL);
SetAutoProcessTH1(fMCHist_Fine2D);
TH2D *tempcov = (TH2D *)input.Get("analysis1_totcov");
fFullCovar = new TMatrixDSym(tempcov->GetNbinsX());
for (int i = 0; i < tempcov->GetNbinsX(); i++) {
for (int j = 0; j < tempcov->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
// Read in 2D Data Slices and Make MC Slices
int bincount = 0;
for (size_t i = 0; i < nangbins; i++) {
fDataHist_Slices.push_back(
(TH1D *)input.Get(Form("dataslice_%lu", i))->Clone());
fDataHist_Slices[i]->SetNameTitle(
Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu", i),
(Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu, cos(#theta) [%f,%f] ",
i, angular_binning_costheta[i],
angular_binning_costheta[i + 1])));
fDataHist_Slices.back()->SetDirectory(NULL);
// Loop over nbins and set errors from covar
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
fDataHist_Slices[i]->SetBinError(
j + 1, sqrt((*fFullCovar)(bincount, bincount)) * 1E-38);
bincount++;
}
}
assert(bincount == tempcov->GetNbinsX());
if (fMaskMomOverflow) {
MaskMomOverflow();
bincount = fFullCovar->GetNcols();
}
- std::vector> data_slice_bcbes;
+ std::vector > data_slice_bcbes;
for (size_t i = 0; i < nangbins; i++) {
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
data_slice_bcbes.push_back(
std::make_pair(fDataHist_Slices[i]->GetBinContent(j + 1),
fDataHist_Slices[i]->GetBinError(j + 1)));
}
}
for (size_t i = 0; i < nangbins; i++) {
fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone());
fMCHist_Slices.back()->SetDirectory(NULL);
fMCHist_Slices.back()->Reset();
fMCHist_Slices.back()->SetLineColor(kRed);
fMCHist_Slices[i]->SetNameTitle(
Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu", i),
(Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu, cos(#theta) [%f,%f] ", i,
angular_binning_costheta[i], angular_binning_costheta[i + 1])));
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
fDataHist =
new TH1D("T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D",
"T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D", bincount, 0, bincount);
fDataHist->SetDirectory(NULL);
for (size_t i = 0; i < data_slice_bcbes.size(); ++i) {
fDataHist->SetBinContent(i + 1, data_slice_bcbes[i].first);
fDataHist->SetBinError(i + 1, data_slice_bcbes[i].second);
}
fMCHist = (TH1D *)fDataHist->Clone("T2K_CC0pi_XSec_2DPcos_nu_I_MC_1D");
fMCHist->SetDirectory(NULL);
return;
}
void T2K_CC0pi_XSec_2DPcos_nu_I::MaskMomOverflow() {
std::vector bins_to_cut;
size_t nallbins = 0;
for (size_t i = 0; i < nangbins; i++) {
std::vector slice_bin_edges;
slice_bin_edges.push_back(
fDataHist_Slices[i]->GetXaxis()->GetBinLowEdge(1));
for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) {
slice_bin_edges.push_back(
fDataHist_Slices[i]->GetXaxis()->GetBinUpEdge(j + 1));
nallbins++;
}
bins_to_cut.push_back(nallbins++);
TH1D *tmp = new TH1D(fDataHist_Slices[i]->GetName(),
fDataHist_Slices[i]->GetTitle(),
slice_bin_edges.size() - 1, slice_bin_edges.data());
tmp->SetDirectory(NULL);
for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) {
tmp->SetBinContent(j + 1, fDataHist_Slices[i]->GetBinContent(j + 1));
tmp->SetBinError(j + 1, fDataHist_Slices[i]->GetBinError(j + 1));
}
delete fDataHist_Slices[i];
fDataHist_Slices[i] = tmp;
}
TMatrixDSym *tmpcovar = new TMatrixDSym(nallbins - bins_to_cut.size());
int icut = 0;
for (int ifull = 0; ifull < fFullCovar->GetNcols(); ifull++) {
if (std::find(bins_to_cut.begin(), bins_to_cut.end(), ifull) !=
bins_to_cut.end()) {
continue;
}
int jcut = 0;
for (int jfull = 0; jfull < fFullCovar->GetNcols(); jfull++) {
if (std::find(bins_to_cut.begin(), bins_to_cut.end(), jfull) !=
bins_to_cut.end()) {
continue;
}
(*tmpcovar)[icut][jcut] = (*fFullCovar)[ifull][jfull];
jcut++;
}
icut++;
}
delete fFullCovar;
fFullCovar = tmpcovar;
}
void T2K_CC0pi_XSec_2DPcos_nu_I::Write(std::string drawopt) {
this->Measurement1D::Write(drawopt);
for (size_t i = 0; i < nangbins; i++) {
fMCHist_Slices[i]->Write();
fDataHist_Slices[i]->Write();
}
if (fResidualHist) {
std::vector MCResidual_Slices;
size_t tb_it = 0;
for (size_t i = 0; i < fMCHist_Slices.size(); ++i) {
std::string name =
Form("T2K_CC0pi_XSec_2DPcos_nu_I_RESIDUAL_Slice%lu", i);
MCResidual_Slices.push_back(
static_cast(fMCHist_Slices[i]->Clone(name.c_str())));
MCResidual_Slices.back()->Reset();
for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) {
double bc = fResidualHist->GetBinContent(tb_it + 1);
MCResidual_Slices.back()->SetBinContent(j + 1, bc);
tb_it++;
}
MCResidual_Slices.back()->Write();
}
}
if (fChi2LessBinHist) {
std::vector MCChi2LessBin_Slices;
size_t tb_it = 0;
for (size_t i = 0; i < fMCHist_Slices.size(); ++i) {
std::string name =
Form("T2K_CC0pi_XSec_2DPcos_nu_I_Chi2NMinusOne_Slice%lu", i);
MCChi2LessBin_Slices.push_back(
static_cast(fMCHist_Slices[i]->Clone(name.c_str())));
MCChi2LessBin_Slices.back()->Reset();
for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) {
double bc = fChi2LessBinHist->GetBinContent(tb_it + 1);
MCChi2LessBin_Slices.back()->SetBinContent(j + 1, bc);
tb_it++;
}
MCChi2LessBin_Slices.back()->Write();
}
}
}
diff --git a/src/T2K/T2K_nueCCinc_XSec_1Dpe.cxx b/src/T2K/T2K_nueCCinc_XSec_1Dpe.cxx
new file mode 100644
index 0000000..18bb0f6
--- /dev/null
+++ b/src/T2K/T2K_nueCCinc_XSec_1Dpe.cxx
@@ -0,0 +1,85 @@
+#include "T2K_nueCCinc_XSec_1Dpe.h"
+#include "T2K_SignalDef.h"
+
+//********************************************************************
+T2K_nueCCinc_XSec_1Dpe::T2K_nueCCinc_XSec_1Dpe(nuiskey samplekey) {
+//********************************************************************
+
+ fSettings = LoadSampleSettings(samplekey);
+ std::string name = fSettings.GetS("name");
+ std::string descrip = "";
+ // This has to deal with nue FHC, and nue/nuebar RHC
+ if (!name.compare("T2K_nueCCinc_XSec_1Dpe_FHC")){
+ descrip = name +". \n"
+ "Target: CH \n"
+ "Flux: T2K FHC nue \n"
+ "Signal: CC-inclusive \n";
+
+ fSettings.SetTitle("T2K FHC #nu_{e}-CC-inclusive p_{e}");
+ fSettings.DefineAllowedSpecies("nue");
+ fSettings.SetDataInput(FitPar::GetDataBase() + "/T2K/CCinc/nue_2019/FHC_nue_unfold_with_neut.txt");
+ nuPDG = 12;
+ lepPDG = 11;
+ } else if (!name.compare("T2K_nueCCinc_XSec_1Dpe_RHC")){
+ descrip = name +". \n"
+ "Target: CH \n"
+ "Flux: T2K RHC nue \n"
+ "Signal: CC-inclusive \n";
+
+ fSettings.SetTitle("T2K RHC #nu_{e}-CC-inclusive p_{e}");
+ fSettings.DefineAllowedSpecies("nue");
+ fSettings.SetDataInput(FitPar::GetDataBase() + "/T2K/CCinc/nue_2019/RHC_nue_unfold_with_neut.txt");
+ nuPDG = 12;
+ lepPDG = 11;
+ } else if (!name.compare("T2K_nuebarCCinc_XSec_1Dpe_RHC")){
+ descrip = name +". \n"
+ "Target: CH \n"
+ "Flux: T2K RHC nuebar \n"
+ "Signal: CC-inclusive \n";
+
+ fSettings.SetTitle("T2K RHC #bar{#nu}_{e}-CC-inclusive p_{e}");
+ fSettings.DefineAllowedSpecies("nueb");
+ fSettings.SetDataInput(FitPar::GetDataBase() + "/T2K/CCinc/nue_2019/RHC_nuebar_unfold_with_neut.txt");
+ nuPDG = -12;
+ lepPDG = -11;
+ }
+
+ // Setup common settings
+ fSettings.SetDescription(descrip);
+ fSettings.SetXTitle("p_{e} (GeV)");
+ fSettings.SetYTitle("#frac{d#sigma}{dp_{e}} (cm^{2}/nucleon)");
+ fSettings.SetEnuRange(0.0, 30.0);
+ fSettings.DefineAllowedTargets("C,H");
+ FinaliseSampleSettings();
+
+ SetDataFromTextFile( fSettings.GetDataInput() );
+
+ // Scaling Setup ---------------------------------------------------
+ // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
+ fScaleFactor = GetEventHistogram()->Integral("width") * 1E-38 / fNEvents /
+ TotalIntegratedFlux();
+
+ // Don't scale by width
+ fIsNoWidth = true;
+
+ FinaliseMeasurement();
+};
+
+// Signal is simply a CC inclusive without any angular/momentum cuts
+bool T2K_nueCCinc_XSec_1Dpe::isSignal(FitEvent *event) {
+ if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false;
+
+ // Phase space restriction
+ TLorentzVector Pnu = event->GetHMISParticle(nuPDG)->fP;
+ TLorentzVector Pe = event->GetHMFSParticle(lepPDG)->fP;
+ static double max_angle = 45*TMath::Pi()/180.;
+ if (FitUtils::th(Pnu, Pe) > max_angle) return false;
+ return true;
+};
+
+void T2K_nueCCinc_XSec_1Dpe::FillEventVariables(FitEvent *event) {
+
+ if (event->NumFSParticle(lepPDG) == 0) return;
+ TLorentzVector Pe = event->GetHMFSParticle(lepPDG)->fP;
+ fXVar = Pe.Vect().Mag() / 1000.;
+};
diff --git a/src/T2K/T2K_nueCCinc_XSec_1Dpe.h b/src/T2K/T2K_nueCCinc_XSec_1Dpe.h
new file mode 100644
index 0000000..e76736b
--- /dev/null
+++ b/src/T2K/T2K_nueCCinc_XSec_1Dpe.h
@@ -0,0 +1,22 @@
+#ifndef T2K_NUECCINC_XSEC_1DPE_H_SEEN
+#define T2K_NUECCINC_XSEC_1DPE_H_SEEN
+
+#include "Measurement1D.h"
+
+class T2K_nueCCinc_XSec_1Dpe : public Measurement1D {
+public:
+
+ T2K_nueCCinc_XSec_1Dpe(nuiskey samplekey);
+
+ ~T2K_nueCCinc_XSec_1Dpe() {};
+
+ bool isSignal(FitEvent *event);
+
+ void FillEventVariables(FitEvent* event);
+
+ private:
+ int nuPDG;
+ int lepPDG;
+};
+
+#endif
diff --git a/src/T2K/T2K_nueCCinc_XSec_joint.cxx b/src/T2K/T2K_nueCCinc_XSec_joint.cxx
new file mode 100644
index 0000000..a11ebf2
--- /dev/null
+++ b/src/T2K/T2K_nueCCinc_XSec_joint.cxx
@@ -0,0 +1,166 @@
+#include "T2K_nueCCinc_XSec_joint.h"
+
+//********************************************************************
+T2K_nueCCinc_XSec_joint::T2K_nueCCinc_XSec_joint(nuiskey samplekey){
+//********************************************************************
+
+ fSettings = LoadSampleSettings(samplekey);
+ std::string descrip = "T2K_nueCCinc_XSec_joint. \n"
+ "Target: CH \n"
+ "Flux: T2K FHC nue \n"
+ "Signal: CC-inclusive \n";
+ fSettings.SetTitle("T2K #nu_{e}-CC-inclusive p_{e} joint");
+ fSettings.DefineAllowedSpecies("nue, nueb");
+ fSettings.SetCovarInput(FitPar::GetDataBase() + "/T2K/CCinc/nue_2019/fract_covar_unfold_with_neut.txt");
+ fSettings.SetDescription(descrip);
+ fSettings.SetXTitle("p_{e} (GeV)");
+ fSettings.SetYTitle("#frac{d#sigma}{dp_{e}} (cm^{2}/nucleon)");
+ fSettings.SetEnuRange(0.0, 30.0);
+ fSettings.DefineAllowedTargets("C,H");
+ FinaliseSampleSettings();
+
+ if (fSubInFiles.size() != 3) {
+ NUIS_ABORT("T2K nue joint requires input files in format: FHC nue; RHC nue; RHC nuebar");
+ }
+
+ std::string inFileFHCNue = fSubInFiles.at(0);
+ std::string inFileRHCNue = fSubInFiles.at(1);
+ std::string inFileRHCNuebar = fSubInFiles.at(2);
+
+ // Create some config keys
+ nuiskey FHCNueKey = Config::CreateKey("sample");
+ FHCNueKey.SetS("input", inFileFHCNue);
+ FHCNueKey.SetS("type", fSettings.GetS("type"));
+ FHCNueKey.SetS("name", "T2K_nueCCinc_XSec_1Dpe_FHC");
+ FHC_nue = new T2K_nueCCinc_XSec_1Dpe(FHCNueKey);
+
+ nuiskey RHCNueKey = Config::CreateKey("sample");
+ RHCNueKey.SetS("input", inFileRHCNue);
+ RHCNueKey.SetS("type", fSettings.GetS("type"));
+ RHCNueKey.SetS("name", "T2K_nueCCinc_XSec_1Dpe_RHC");
+ RHC_nue = new T2K_nueCCinc_XSec_1Dpe(RHCNueKey);
+
+ nuiskey RHCNuebarKey = Config::CreateKey("sample");
+ RHCNuebarKey.SetS("input", inFileRHCNuebar);
+ RHCNuebarKey.SetS("type", fSettings.GetS("type"));
+ RHCNuebarKey.SetS("name", "T2K_nuebarCCinc_XSec_1Dpe_RHC");
+ RHC_nuebar = new T2K_nueCCinc_XSec_1Dpe(RHCNuebarKey);
+
+ // Sort out the data hist
+ this->CombineDataHists();
+
+ // This is a fractional covariance. Need to account for that
+ SetFractCovarFromTextFile(fSettings.GetCovarInput());
+ ScaleCovar(1E76);
+ SetShapeCovar();
+
+ // Add to chain for processing
+ fSubChain.clear();
+ fSubChain.push_back(FHC_nue);
+ fSubChain.push_back(RHC_nue);
+ fSubChain.push_back(RHC_nuebar);
+
+ // This saves information from the sub-measurements
+ fSaveSubMeas = true;
+ FinaliseMeasurement();
+};
+
+//********************************************************************
+void T2K_nueCCinc_XSec_joint::SetFractCovarFromTextFile(std::string covfile){
+//********************************************************************
+
+ if (!fDataHist){
+ NUIS_ERR(FTL, "fDataHist is required to call T2K_nueCCinc_XSec_joint::SetFractCovarFromTextFile");
+ throw;
+ }
+
+ NUIS_LOG(SAM, "Reading fractional covariance from text file: " << covfile);
+
+ int dim = fDataHist->GetNbinsX();
+ TMatrixD *tempmat = StatUtils::GetMatrixFromTextFile(covfile, dim, dim);
+
+ // Make a symmetric covariance
+ fFullCovar = new TMatrixDSym(tempmat->GetNrows());
+ for (int i = 0; i < tempmat->GetNrows(); i++) {
+ for (int j = 0; j < tempmat->GetNrows(); j++) {
+ (*fFullCovar)(i, j) = fDataHist->GetBinContent(i+1)
+ *(*tempmat)(i, j)*fDataHist->GetBinContent(j+1);
+ }
+ }
+
+ covar = StatUtils::GetInvert(fFullCovar);
+ fDecomp = StatUtils::GetDecomp(fFullCovar);
+}
+
+
+//********************************************************************
+void T2K_nueCCinc_XSec_joint::CombineDataHists(){
+//********************************************************************
+
+ TH1D *FHC_nue_data = (TH1D*)FHC_nue->GetDataHistogram();
+ TH1D *RHC_nue_data = (TH1D*)RHC_nue->GetDataHistogram();
+ TH1D *RHC_nuebar_data = (TH1D*)RHC_nuebar->GetDataHistogram();
+
+ int nbins = FHC_nue_data->GetNbinsX() +
+ RHC_nue_data->GetNbinsX() +
+ RHC_nuebar_data->GetNbinsX();
+
+ fDataHist = new TH1D((fSettings.GetName() + "_data").c_str(),
+ (fSettings.GetFullTitles()).c_str(), nbins, 0, nbins);
+ fDataHist->SetDirectory(0);
+
+ // Bit ugly, but...
+ int count = 0;
+ for (int x=0; xGetNbinsX(); ++x){
+ fDataHist->SetBinContent(count+1, FHC_nue_data->GetBinContent(x+1));
+ fDataHist->SetBinError(count+1, FHC_nue_data->GetBinError(x+1));
+ fDataHist->GetXaxis()->SetBinLabel(count+1, Form("FHC #nu_{e} %.1f-%.1f",
+ FHC_nue_data->GetXaxis()->GetBinLowEdge(x+1),
+ FHC_nue_data->GetXaxis()->GetBinUpEdge(x+1)));
+ count++;
+ }
+ for (int x=0; xGetNbinsX(); ++x){
+ fDataHist->SetBinContent(count+1, RHC_nue_data->GetBinContent(x+1));
+ fDataHist->SetBinError(count+1, RHC_nue_data->GetBinError(x+1));
+ fDataHist->GetXaxis()->SetBinLabel(count+1, Form("RHC #nu_{e} %.1f-%.1f",
+ RHC_nue_data->GetXaxis()->GetBinLowEdge(x+1),
+ RHC_nue_data->GetXaxis()->GetBinUpEdge(x+1)));
+ count++;
+ }
+ for (int x=0; xGetNbinsX(); ++x){
+ fDataHist->SetBinContent(count+1, RHC_nuebar_data->GetBinContent(x+1));
+ fDataHist->SetBinError(count+1, RHC_nuebar_data->GetBinError(x+1));
+ fDataHist->GetXaxis()->SetBinLabel(count+1, Form("RHC #bar{#nu}_{e} %.1f-%.1f",
+ RHC_nuebar_data->GetXaxis()->GetBinLowEdge(x+1),
+ RHC_nuebar_data->GetXaxis()->GetBinUpEdge(x+1)));
+ count++;
+ }
+}
+
+//********************************************************************
+void T2K_nueCCinc_XSec_joint::MakePlots() {
+//********************************************************************
+
+ TH1D *FHC_nue_mc = (TH1D*)FHC_nue->GetMCHistogram();
+ TH1D *RHC_nue_mc = (TH1D*)RHC_nue->GetMCHistogram();
+ TH1D *RHC_nuebar_mc = (TH1D*)RHC_nuebar->GetMCHistogram();
+
+ int count = 0;
+ for (int i = 0; i < FHC_nue_mc->GetNbinsX(); ++i) {
+ fMCHist->SetBinContent(count + 1, FHC_nue_mc->GetBinContent(i + 1));
+ fMCHist->SetBinError(count + 1, FHC_nue_mc->GetBinError(i + 1));
+ count++;
+ }
+ for (int i = 0; i < RHC_nue_mc->GetNbinsX(); ++i) {
+ fMCHist->SetBinContent(count + 1, RHC_nue_mc->GetBinContent(i + 1));
+ fMCHist->SetBinError(count + 1, RHC_nue_mc->GetBinError(i + 1));
+ count++;
+ }
+ for (int i = 0; i < RHC_nuebar_mc->GetNbinsX(); ++i) {
+ fMCHist->SetBinContent(count + 1, RHC_nuebar_mc->GetBinContent(i + 1));
+ fMCHist->SetBinError(count + 1, RHC_nuebar_mc->GetBinError(i + 1));
+ count++;
+ }
+
+ return;
+}
diff --git a/src/T2K/T2K_nueCCinc_XSec_joint.h b/src/T2K/T2K_nueCCinc_XSec_joint.h
new file mode 100644
index 0000000..3cd1210
--- /dev/null
+++ b/src/T2K/T2K_nueCCinc_XSec_joint.h
@@ -0,0 +1,34 @@
+#ifndef T2K_NUECCINC_XSEC_JOINT_H_SEEN
+#define T2K_NUECCINC_XSEC_JOINT_H_SEEN
+
+#include "JointMeas1D.h"
+#include "Measurement1D.h"
+#include "T2K_nueCCinc_XSec_1Dpe.h"
+
+class T2K_nueCCinc_XSec_joint : public JointMeas1D {
+public:
+
+ T2K_nueCCinc_XSec_joint(nuiskey samplekey);
+ virtual ~T2K_nueCCinc_XSec_joint() {};
+
+ void MakePlots();
+
+ // Makes a data hist from all the smaller ones
+ void CombineDataHists();
+
+ // Interprets the fractional covariance provided
+ void SetFractCovarFromTextFile(std::string covfile);
+
+ // No need for these
+ void ScaleEvents(){return;};
+ bool isSignal(){return false;};
+ void FillEventVariables(){return;};
+
+ private:
+ // The separate measurements that go into this
+ Measurement1D* FHC_nue;
+ Measurement1D* RHC_nue;
+ Measurement1D* RHC_nuebar;
+};
+
+#endif
diff --git a/src/Utils/PlotUtils.cxx b/src/Utils/PlotUtils.cxx
index 612c493..1d641ff 100644
--- a/src/Utils/PlotUtils.cxx
+++ b/src/Utils/PlotUtils.cxx
@@ -1,1206 +1,1206 @@
// 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 "PlotUtils.h"
#include "FitEvent.h"
#include "StatUtils.h"
// MOVE TO MB UTILS!
// This function is intended to be modified to enforce a consistent masking for
// all models.
TH2D *PlotUtils::SetMaskHist(std::string type, TH2D *data) {
TH2D *fMaskHist = (TH2D *)data->Clone("fMaskHist");
for (int xBin = 0; xBin < fMaskHist->GetNbinsX(); ++xBin) {
for (int yBin = 0; yBin < fMaskHist->GetNbinsY(); ++yBin) {
if (data->GetBinContent(xBin + 1, yBin + 1) == 0)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 0);
else
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 0.5);
if (!type.compare("MB_numu_2D")) {
if (yBin == 19 && xBin < 8)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
} else {
if (yBin == 19 && xBin < 11)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
}
if (yBin == 18 && xBin < 3)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
if (yBin == 17 && xBin < 2)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
if (yBin == 16 && xBin < 1)
fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
}
}
return fMaskHist;
};
// MOVE TO GENERAL UTILS?
bool PlotUtils::CheckObjectWithName(TFile *inFile, std::string substring) {
TIter nextkey(inFile->GetListOfKeys());
TKey *key;
while ((key = (TKey *)nextkey())) {
std::string test(key->GetName());
if (test.find(substring) != std::string::npos)
return true;
}
return false;
};
// MOVE TO GENERAL UTILS?
std::string PlotUtils::GetObjectWithName(TFile *inFile, std::string substring) {
TIter nextkey(inFile->GetListOfKeys());
TKey *key;
std::string output = "";
while ((key = (TKey *)nextkey())) {
std::string test(key->GetName());
if (test.find(substring) != std::string::npos)
output = test;
}
return output;
};
void PlotUtils::CreateNeutModeArray(TH1 *hist, TH1 *neutarray[]) {
for (int i = 0; i < 60; i++) {
neutarray[i] = (TH1 *)hist->Clone(Form("%s_NMODE_%i", hist->GetName(), i));
}
return;
};
void PlotUtils::DeleteNeutModeArray(TH1 *neutarray[]) {
for (int i = 0; i < 60; i++) {
delete neutarray[i];
}
return;
};
void PlotUtils::FillNeutModeArray(TH1D *hist[], int mode, double xval,
double weight) {
if (abs(mode) > 60)
return;
hist[abs(mode)]->Fill(xval, weight);
return;
};
void PlotUtils::FillNeutModeArray(TH2D *hist[], int mode, double xval,
double yval, double weight) {
if (abs(mode) > 60)
return;
hist[abs(mode)]->Fill(xval, yval, weight);
return;
};
THStack PlotUtils::GetNeutModeStack(std::string title, TH1 *ModeStack[],
int option) {
(void)option;
THStack allmodes = THStack(title.c_str(), title.c_str());
for (int i = 0; i < 60; i++) {
allmodes.Add(ModeStack[i]);
}
// Credit to Clarence for copying all this out.
// CC
ModeStack[1]->SetTitle("CCQE");
ModeStack[1]->SetFillColor(kBlue);
// ModeStack[1]->SetFillStyle(3444);
ModeStack[1]->SetLineColor(kBlue);
ModeStack[2]->SetTitle("2p/2h Nieves");
ModeStack[2]->SetFillColor(kRed);
// ModeStack[2]->SetFillStyle(3344);
ModeStack[2]->SetLineColor(kRed);
// ModeStack[11]->SetTitle("#it{#nu + p #rightarrow l^{-} + p + #pi^{+}}");
ModeStack[11]->SetTitle("CC1#pi^{+} on p");
ModeStack[11]->SetFillColor(kGreen);
// ModeStack[11]->SetFillStyle(3004);
ModeStack[11]->SetLineColor(kGreen);
// ModeStack[12]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #pi^{0}}");
ModeStack[12]->SetTitle("CC1#pi^{0} on n");
ModeStack[12]->SetFillColor(kGreen + 3);
// ModeStack[12]->SetFillStyle(3005);
ModeStack[12]->SetLineColor(kGreen);
// ModeStack[13]->SetTitle("#it{#nu + n #rightarrow l^{-} + n + #pi^{+}}");
ModeStack[13]->SetTitle("CC1#pi^{+} on n");
ModeStack[13]->SetFillColor(kGreen - 2);
// ModeStack[13]->SetFillStyle(3004);
ModeStack[13]->SetLineColor(kGreen);
ModeStack[16]->SetTitle("CC coherent");
ModeStack[16]->SetFillColor(kBlue);
// ModeStack[16]->SetFillStyle(3644);
ModeStack[16]->SetLineColor(kBlue);
// ModeStack[17]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #gamma}");
ModeStack[17]->SetTitle("CC1#gamma");
ModeStack[17]->SetFillColor(kMagenta);
// ModeStack[17]->SetFillStyle(3001);
ModeStack[17]->SetLineColor(kMagenta);
ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)");
ModeStack[21]->SetFillColor(kYellow);
// ModeStack[21]->SetFillStyle(3005);
ModeStack[21]->SetLineColor(kYellow);
// ModeStack[22]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #eta^{0}}");
ModeStack[22]->SetTitle("CC1#eta^{0} on n");
ModeStack[22]->SetFillColor(kYellow - 2);
// ModeStack[22]->SetFillStyle(3013);
ModeStack[22]->SetLineColor(kYellow - 2);
// ModeStack[23]->SetTitle("#it{#nu + n #rightarrow l^{-} + #Lambda +
// K^{+}}");
ModeStack[23]->SetTitle("CC1#Labda1K^{+}");
ModeStack[23]->SetFillColor(kYellow - 6);
// ModeStack[23]->SetFillStyle(3013);
ModeStack[23]->SetLineColor(kYellow - 6);
ModeStack[26]->SetTitle("DIS (W > 2.0)");
ModeStack[26]->SetFillColor(kRed);
// ModeStack[26]->SetFillStyle(3006);
ModeStack[26]->SetLineColor(kRed);
// NC
// ModeStack[31]->SetTitle("#it{#nu + n #rightarrow #nu + n + #pi^{0}}");
ModeStack[31]->SetTitle("NC1#pi^{0} on n");
ModeStack[31]->SetFillColor(kBlue);
// ModeStack[31]->SetFillStyle(3004);
ModeStack[31]->SetLineColor(kBlue);
// ModeStack[32]->SetTitle("#it{#nu + p #rightarrow #nu + p + #pi^{0}}");
ModeStack[32]->SetTitle("NC1#pi^{0} on p");
ModeStack[32]->SetFillColor(kBlue + 3);
// ModeStack[32]->SetFillStyle(3004);
ModeStack[32]->SetLineColor(kBlue + 3);
// ModeStack[33]->SetTitle("#it{#nu + n #rightarrow #nu + p + #pi^{-}}");
ModeStack[33]->SetTitle("NC1#pi^{-} on n");
ModeStack[33]->SetFillColor(kBlue - 2);
// ModeStack[33]->SetFillStyle(3005);
ModeStack[33]->SetLineColor(kBlue - 2);
// ModeStack[34]->SetTitle("#it{#nu + p #rightarrow #nu + n + #pi^{+}}");
ModeStack[34]->SetTitle("NC1#pi^{+} on p");
ModeStack[34]->SetFillColor(kBlue - 8);
// ModeStack[34]->SetFillStyle(3005);
ModeStack[34]->SetLineColor(kBlue - 8);
ModeStack[36]->SetTitle("NC Coherent");
ModeStack[36]->SetFillColor(kBlue + 8);
// ModeStack[36]->SetFillStyle(3644);
ModeStack[36]->SetLineColor(kBlue + 8);
// ModeStack[38]->SetTitle("#it{#nu + n #rightarrow #nu + n + #gamma}");
ModeStack[38]->SetTitle("NC1#gamma on n");
ModeStack[38]->SetFillColor(kMagenta);
// ModeStack[38]->SetFillStyle(3001);
ModeStack[38]->SetLineColor(kMagenta);
// ModeStack[39]->SetTitle("#it{#nu + p #rightarrow #nu + p + #gamma}");
ModeStack[39]->SetTitle("NC1#gamma on p");
ModeStack[39]->SetFillColor(kMagenta - 10);
// ModeStack[39]->SetFillStyle(3001);
ModeStack[39]->SetLineColor(kMagenta - 10);
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
ModeStack[41]->SetFillColor(kBlue - 10);
// ModeStack[41]->SetFillStyle(3005);
ModeStack[41]->SetLineColor(kBlue - 10);
// ModeStack[42]->SetTitle("#it{#nu + n #rightarrow #nu + n + #eta^{0}}");
ModeStack[42]->SetTitle("NC1#eta^{0} on n");
ModeStack[42]->SetFillColor(kYellow - 2);
// ModeStack[42]->SetFillStyle(3013);
ModeStack[42]->SetLineColor(kYellow - 2);
// ModeStack[43]->SetTitle("#it{#nu + p #rightarrow #nu + p + #eta^{0}}");
ModeStack[43]->SetTitle("NC1#eta^{0} on p");
ModeStack[43]->SetFillColor(kYellow - 4);
// ModeStack[43]->SetFillStyle(3013);
ModeStack[43]->SetLineColor(kYellow - 4);
// ModeStack[44]->SetTitle("#it{#nu + n #rightarrow #nu + #Lambda + K^{0}}");
ModeStack[44]->SetTitle("NC1#Lambda1K^{0} on n");
ModeStack[44]->SetFillColor(kYellow - 6);
// ModeStack[44]->SetFillStyle(3014);
ModeStack[44]->SetLineColor(kYellow - 6);
// ModeStack[45]->SetTitle("#it{#nu + p #rightarrow #nu + #Lambda + K^{+}}");
ModeStack[45]->SetTitle("NC1#Lambda1K^{+}");
ModeStack[45]->SetFillColor(kYellow - 10);
// ModeStack[45]->SetFillStyle(3014);
ModeStack[45]->SetLineColor(kYellow - 10);
ModeStack[46]->SetTitle("DIS (W > 2.0)");
ModeStack[46]->SetFillColor(kRed);
// ModeStack[46]->SetFillStyle(3006);
ModeStack[46]->SetLineColor(kRed);
// ModeStack[51]->SetTitle("#it{#nu + p #rightarrow #nu + p}");
ModeStack[51]->SetTitle("NC on p");
ModeStack[51]->SetFillColor(kBlack);
// ModeStack[51]->SetFillStyle(3444);
ModeStack[51]->SetLineColor(kBlack);
// ModeStack[52]->SetTitle("#it{#nu + n #rightarrow #nu + n}");
ModeStack[52]->SetTitle("NC on n");
ModeStack[52]->SetFillColor(kGray);
// ModeStack[52]->SetFillStyle(3444);
ModeStack[52]->SetLineColor(kGray);
return allmodes;
};
TLegend PlotUtils::GenerateStackLegend(THStack stack, int xlow, int ylow,
int xhigh, int yhigh) {
TLegend leg = TLegend(xlow, ylow, xhigh, yhigh);
TObjArray *histarray = stack.GetStack();
int nhist = histarray->GetEntries();
for (int i = 0; i < nhist; i++) {
TH1 *hist = (TH1 *)(histarray->At(i));
leg.AddEntry((hist), ((TH1 *)histarray->At(i))->GetTitle(), "fl");
}
leg.SetName(Form("%s_LEG", stack.GetName()));
return leg;
};
void PlotUtils::ScaleNeutModeArray(TH1 *hist[], double factor,
std::string option) {
for (int i = 0; i < 60; i++) {
if (hist[i])
hist[i]->Scale(factor, option.c_str());
}
return;
};
void PlotUtils::ResetNeutModeArray(TH1 *hist[]) {
for (int i = 0; i < 60; i++) {
if (hist[i])
hist[i]->Reset();
}
return;
};
//********************************************************************
// This assumes the Enu axis is the x axis, as is the case for MiniBooNE 2D
// distributions
void PlotUtils::FluxUnfoldedScaling(TH2D *fMCHist, TH1D *fhist, TH1D *ehist,
double scalefactor) {
//********************************************************************
// Make clones to avoid changing stuff
TH1D *eventhist = (TH1D *)ehist->Clone();
TH1D *fFluxHist = (TH1D *)fhist->Clone();
// Undo width integral in SF
fMCHist->Scale(scalefactor /
eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width"));
// Standardise The Flux
eventhist->Scale(1.0 / fFluxHist->Integral());
fFluxHist->Scale(1.0 / fFluxHist->Integral());
// Do interpolation for 2D plots?
// fFluxHist = PlotUtils::InterpolateFineHistogram(fFluxHist,100,"width");
// eventhist = PlotUtils::InterpolateFineHistogram(eventhist,100,"width");
// eventhist->Scale(1.0/fFluxHist->Integral());
// fFluxHist->Scale(1.0/fFluxHist->Integral());
// Scale fMCHist by eventhist integral
fMCHist->Scale(eventhist->Integral(1, eventhist->GetNbinsX() + 1));
// Find which axis is the Enu axis
bool EnuOnXaxis = false;
std::string xaxis = fMCHist->GetXaxis()->GetTitle();
if (xaxis.find("E") != std::string::npos &&
xaxis.find("nu") != std::string::npos)
EnuOnXaxis = true;
std::string yaxis = fMCHist->GetYaxis()->GetTitle();
if (yaxis.find("E") != std::string::npos &&
xaxis.find("nu") != std::string::npos) {
// First check that xaxis didn't also find Enu
if (EnuOnXaxis) {
NUIS_ERR(FTL, fMCHist->GetTitle() << " error:");
NUIS_ERR(FTL, "Found Enu in xaxis title: " << xaxis);
NUIS_ERR(FTL, "AND");
NUIS_ERR(FTL, "Found Enu in yaxis title: " << yaxis);
NUIS_ABORT("Enu on x and Enu on y flux unfolded scaling isn't "
"implemented, please modify "
<< __FILE__ << ":" << __LINE__);
}
EnuOnXaxis = false;
}
// Now Get a flux PDF assuming X axis is Enu
TH1D *pdfflux = NULL;
// If xaxis is Enu
if (EnuOnXaxis)
pdfflux = (TH1D *)fMCHist->ProjectionX()->Clone();
// If yaxis is Enu
else
pdfflux = (TH1D *)fMCHist->ProjectionY()->Clone();
// pdfflux->Write( (std::string(fMCHist->GetName()) + "_PROJX").c_str());
pdfflux->Reset();
// Awful MiniBooNE Check for the time being
// Needed because the flux is in GeV whereas the measurement is in MeV
bool ismb =
std::string(fMCHist->GetName()).find("MiniBooNE") != std::string::npos;
for (int i = 0; i < pdfflux->GetNbinsX(); i++) {
double Ml = pdfflux->GetXaxis()->GetBinLowEdge(i + 1);
double Mh = pdfflux->GetXaxis()->GetBinLowEdge(i + 2);
// double Mc = pdfflux->GetXaxis()->GetBinCenter(i+1);
// double Mw = pdfflux->GetBinWidth(i+1);
double fluxint = 0.0;
// Scaling to match flux for MB
if (ismb) {
Ml /= 1.E3;
Mh /= 1.E3;
// Mc /= 1.E3;
// Mw /= 1.E3;
}
for (int j = 0; j < fFluxHist->GetNbinsX(); j++) {
// double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1);
double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1);
double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2);
double Fe = fFluxHist->GetBinContent(j + 1);
double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1);
if (Fl >= Ml and Fh <= Mh) {
fluxint += Fe;
} else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) {
fluxint += Fe * (Fh - Ml) / Fw;
} else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) {
fluxint += Fe * (Mh - Fl) / Fw;
} else if (Ml >= Fl and Mh <= Fh) {
fluxint += Fe * (Mh - Ml) / Fw;
} else {
continue;
}
}
pdfflux->SetBinContent(i + 1, fluxint);
}
// Then finally divide by the bin-width in
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
for (int j = 0; j < fMCHist->GetNbinsY(); j++) {
if (pdfflux->GetBinContent(i + 1) == 0.0)
continue;
// Different scaling depending on if Enu is on x or y axis
double scaling = 1.0;
// If Enu is on the x-axis, we want the ith entry of the flux
// And to divide by the bin width of the jth bin
if (EnuOnXaxis) {
double binWidth = fMCHist->GetYaxis()->GetBinLowEdge(j + 2) -
fMCHist->GetYaxis()->GetBinLowEdge(j + 1);
scaling = pdfflux->GetBinContent(i + 1) * binWidth;
} else {
double binWidth = fMCHist->GetXaxis()->GetBinLowEdge(i + 2) -
fMCHist->GetXaxis()->GetBinLowEdge(i + 1);
scaling = pdfflux->GetBinContent(j + 1) * binWidth;
}
// fMCHist->SetBinContent(i + 1, j + 1,
// fMCHist->GetBinContent(i + 1, j + 1) /
// pdfflux->GetBinContent(i + 1) / binWidth);
// fMCHist->SetBinError(i + 1, j + 1, fMCHist->GetBinError(i + 1, j + 1) /
// pdfflux->GetBinContent(i + 1) /
// binWidth);
fMCHist->SetBinContent(i + 1, j + 1,
fMCHist->GetBinContent(i + 1, j + 1) / scaling);
fMCHist->SetBinError(i + 1, j + 1,
fMCHist->GetBinError(i + 1, j + 1) / scaling);
}
}
delete eventhist;
delete fFluxHist;
};
TH1D *PlotUtils::InterpolateFineHistogram(TH1D *hist, int res,
std::string opt) {
int nbins = hist->GetNbinsX();
double elow = hist->GetXaxis()->GetBinLowEdge(1);
double ehigh = hist->GetXaxis()->GetBinLowEdge(nbins + 1);
bool width = true; // opt.find("width") != std::string::npos;
TH1D *fine = new TH1D("fine", "fine", nbins * res, elow, ehigh);
TGraph *temp = new TGraph();
for (int i = 0; i < nbins; i++) {
double E = hist->GetXaxis()->GetBinCenter(i + 1);
double C = hist->GetBinContent(i + 1);
double W = hist->GetXaxis()->GetBinWidth(i + 1);
if (!width)
W = 1.0;
if (W != 0.0)
temp->SetPoint(temp->GetN(), E, C / W);
}
for (int i = 0; i < fine->GetNbinsX(); i++) {
double E = fine->GetXaxis()->GetBinCenter(i + 1);
double W = fine->GetBinWidth(i + 1);
if (!width)
W = 1.0;
fine->SetBinContent(i + 1, temp->Eval(E, 0, "S") * W);
}
fine->Scale(hist->Integral(1, hist->GetNbinsX() + 1) /
fine->Integral(1, fine->GetNbinsX() + 1));
// std::cout << "Interpolation Difference = "
//<< fine->Integral(1, fine->GetNbinsX() + 1) << "/"
//<< hist->Integral(1, hist->GetNbinsX() + 1) << std::endl;
return fine;
}
//********************************************************************
// This interpolates the flux by a TGraph instead of requiring the flux and MC
// flux to have the same binning
void PlotUtils::FluxUnfoldedScaling(TH1D *mcHist, TH1D *fhist, TH1D *ehist,
double scalefactor, int nevents) {
//********************************************************************
TH1D *eventhist = (TH1D *)ehist->Clone();
TH1D *fFluxHist = (TH1D *)fhist->Clone();
if (FitPar::Config().GetParB("save_flux_debug")) {
std::string name = std::string(mcHist->GetName());
mcHist->Write((name + "_UNF_MC").c_str());
fFluxHist->Write((name + "_UNF_FLUX").c_str());
eventhist->Write((name + "_UNF_EVT").c_str());
TH1D *scalehist = new TH1D("scalehist", "scalehist", 1, 0.0, 1.0);
scalehist->SetBinContent(1, scalefactor);
scalehist->SetBinContent(2, nevents);
scalehist->Write((name + "_UNF_SCALE").c_str());
}
// Undo width integral in SF
mcHist->Scale(scalefactor /
eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width"));
// Standardise The Flux
eventhist->Scale(1.0 / fFluxHist->Integral());
fFluxHist->Scale(1.0 / fFluxHist->Integral());
// Scale mcHist by eventhist integral
mcHist->Scale(eventhist->Integral(1, eventhist->GetNbinsX() + 1));
// Now Get a flux PDF
TH1D *pdfflux = (TH1D *)mcHist->Clone();
pdfflux->Reset();
for (int i = 0; i < mcHist->GetNbinsX(); i++) {
double Ml = mcHist->GetXaxis()->GetBinLowEdge(i + 1);
double Mh = mcHist->GetXaxis()->GetBinLowEdge(i + 2);
// double Mc = mcHist->GetXaxis()->GetBinCenter(i+1);
// double Me = mcHist->GetBinContent(i+1);
// double Mw = mcHist->GetBinWidth(i+1);
double fluxint = 0.0;
for (int j = 0; j < fFluxHist->GetNbinsX(); j++) {
// double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1);
double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1);
double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2);
double Fe = fFluxHist->GetBinContent(j + 1);
double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1);
if (Fl >= Ml and Fh <= Mh) {
fluxint += Fe;
} else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) {
fluxint += Fe * (Fh - Ml) / Fw;
} else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) {
fluxint += Fe * (Mh - Fl) / Fw;
} else if (Ml >= Fl and Mh <= Fh) {
fluxint += Fe * (Mh - Ml) / Fw;
} else {
continue;
}
}
pdfflux->SetBinContent(i + 1, fluxint);
}
// Scale MC hist by pdfflux
for (int i = 0; i < mcHist->GetNbinsX(); i++) {
if (pdfflux->GetBinContent(i + 1) == 0.0)
continue;
mcHist->SetBinContent(i + 1, mcHist->GetBinContent(i + 1) /
pdfflux->GetBinContent(i + 1));
mcHist->SetBinError(i + 1, mcHist->GetBinError(i + 1) /
pdfflux->GetBinContent(i + 1));
}
delete eventhist;
delete fFluxHist;
};
// MOVE TO GENERAL UTILS
//********************************************************************
void PlotUtils::Set2DHistFromText(std::string dataFile, TH2 *hist, double norm,
bool skipbins) {
//********************************************************************
std::string line;
std::ifstream data(dataFile.c_str(), std::ifstream::in);
int yBin = 0;
while (std::getline(data >> std::ws, line, '\n')) {
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
// Loop over entries and insert them into the histogram
for (uint xBin = 0; xBin < entries.size(); xBin++) {
if (!skipbins || entries[xBin] != -1.0)
hist->SetBinContent(xBin + 1, yBin + 1, entries[xBin] * norm);
}
yBin++;
}
return;
}
// MOVE TO GENERAL UTILS
TH1D *PlotUtils::GetTH1DFromFile(std::string dataFile, std::string title,
std::string fPlotTitles,
std::string alt_name) {
TH1D *tempPlot;
// If format is a root file
if (dataFile.find(".root") != std::string::npos) {
TFile *temp_infile = new TFile(dataFile.c_str(), "READ");
tempPlot = (TH1D *)temp_infile->Get(title.c_str());
tempPlot->SetDirectory(0);
temp_infile->Close();
delete temp_infile;
// Else its a space seperated txt file
} else {
// Make a TGraph Errors
TGraphErrors *gr = new TGraphErrors(dataFile.c_str(), "%lg %lg %lg");
if (gr->IsZombie()) {
NUIS_ABORT(
dataFile
<< " is a zombie and could not be read. Are you sure it exists?"
<< std::endl);
}
double *bins = gr->GetX();
double *values = gr->GetY();
double *errors = gr->GetEY();
int npoints = gr->GetN();
// Fill the histogram from it
tempPlot = new TH1D(title.c_str(), title.c_str(), npoints - 1, bins);
for (int i = 0; i < npoints; ++i) {
tempPlot->SetBinContent(i + 1, values[i]);
// If only two columns are present in the input file, use the sqrt(values)
// as the error equivalent to assuming that the error is statistical. Also
// check that we're looking at an event rate rather than a cross section
if (!errors[i] && values[i] > 1E-30) {
tempPlot->SetBinError(i + 1, sqrt(values[i]));
} else {
tempPlot->SetBinError(i + 1, errors[i]);
}
}
delete gr;
}
// Allow alternate naming for root files
if (!alt_name.empty()) {
tempPlot->SetNameTitle(alt_name.c_str(), alt_name.c_str());
}
// Allow alternate axis titles
if (!fPlotTitles.empty()) {
tempPlot->SetNameTitle(
tempPlot->GetName(),
(std::string(tempPlot->GetTitle()) + fPlotTitles).c_str());
}
return tempPlot;
};
TH1D *PlotUtils::GetRatioPlot(TH1D *hist1, TH1D *hist2) {
// make copy of first hist
TH1D *new_hist = (TH1D *)hist1->Clone();
// Do bins and errors ourselves as scales can go awkward
for (int i = 0; i < new_hist->GetNbinsX(); i++) {
if (hist2->GetBinContent(i + 1) == 0.0) {
new_hist->SetBinContent(i + 1, 0.0);
}
new_hist->SetBinContent(i + 1, hist1->GetBinContent(i + 1) /
hist2->GetBinContent(i + 1));
new_hist->SetBinError(i + 1, hist1->GetBinError(i + 1) /
hist2->GetBinContent(i + 1));
}
return new_hist;
};
TH1D *PlotUtils::GetRenormalisedPlot(TH1D *hist1, TH1D *hist2) {
// make copy of first hist
TH1D *new_hist = (TH1D *)hist1->Clone();
if (hist1->Integral("width") == 0 or hist2->Integral("width") == 0) {
new_hist->Reset();
return new_hist;
}
Double_t scaleF = hist2->Integral("width") / hist1->Integral("width");
new_hist->Scale(scaleF);
return new_hist;
};
TH1D *PlotUtils::GetShapePlot(TH1D *hist1) {
// make copy of first hist
TH1D *new_hist = (TH1D *)hist1->Clone();
if (hist1->Integral("width") == 0) {
new_hist->Reset();
return new_hist;
}
Double_t scaleF1 = 1.0 / hist1->Integral("width");
new_hist->Scale(scaleF1);
return new_hist;
};
TH1D *PlotUtils::GetShapeRatio(TH1D *hist1, TH1D *hist2) {
TH1D *new_hist1 = GetShapePlot(hist1);
TH1D *new_hist2 = GetShapePlot(hist2);
// Do bins and errors ourselves as scales can go awkward
for (int i = 0; i < new_hist1->GetNbinsX(); i++) {
if (hist2->GetBinContent(i + 1) == 0) {
new_hist1->SetBinContent(i + 1, 0.0);
}
new_hist1->SetBinContent(i + 1, new_hist1->GetBinContent(i + 1) /
new_hist2->GetBinContent(i + 1));
new_hist1->SetBinError(i + 1, new_hist1->GetBinError(i + 1) /
new_hist2->GetBinContent(i + 1));
}
delete new_hist2;
return new_hist1;
};
TH2D *PlotUtils::GetCovarPlot(TMatrixDSym *cov, std::string name,
std::string title) {
TH2D *CovarPlot;
if (cov)
CovarPlot = new TH2D((*cov));
else
CovarPlot = new TH2D(name.c_str(), title.c_str(), 1, 0, 1, 1, 0, 1);
CovarPlot->SetName(name.c_str());
CovarPlot->SetTitle(title.c_str());
return CovarPlot;
}
TH2D *PlotUtils::GetFullCovarPlot(TMatrixDSym *cov, std::string name) {
return PlotUtils::GetCovarPlot(
cov, name + "_COV", name + "_COV;Bins;Bins;Covariance (#times10^{-76})");
}
TH2D *PlotUtils::GetInvCovarPlot(TMatrixDSym *cov, std::string name) {
return PlotUtils::GetCovarPlot(
cov, name + "_INVCOV",
name + "_INVCOV;Bins;Bins;Inv. Covariance (#times10^{-76})");
}
TH2D *PlotUtils::GetDecompCovarPlot(TMatrixDSym *cov, std::string name) {
return PlotUtils::GetCovarPlot(
cov, name + "_DECCOV",
name + "_DECCOV;Bins;Bins;Decomp Covariance (#times10^{-76})");
}
TH1D *PlotUtils::GetTH1DFromRootFile(std::string file, std::string name) {
if (name.empty()) {
std::vector tempfile = GeneralUtils::ParseToStr(file, ";");
file = tempfile[0];
name = tempfile[1];
}
TFile *rootHistFile = new TFile(file.c_str(), "READ");
TH1D *tempHist = (TH1D *)rootHistFile->Get(name.c_str())->Clone();
if (tempHist == NULL) {
NUIS_ABORT("Could not find distribution " << name << " in file " << file);
}
tempHist->SetDirectory(0);
rootHistFile->Close();
return tempHist;
}
TH2D *PlotUtils::GetTH2DFromRootFile(std::string file, std::string name) {
if (name.empty()) {
std::vector tempfile = GeneralUtils::ParseToStr(file, ";");
file = tempfile[0];
name = tempfile[1];
}
TFile *rootHistFile = new TFile(file.c_str(), "READ");
TH2D *tempHist = (TH2D *)rootHistFile->Get(name.c_str())->Clone();
tempHist->SetDirectory(0);
rootHistFile->Close();
delete rootHistFile;
return tempHist;
}
TH1 *PlotUtils::GetTH1FromRootFile(std::string file, std::string name) {
if (name.empty()) {
std::vector tempfile = GeneralUtils::ParseToStr(file, ";");
file = tempfile[0];
name = tempfile[1];
}
TFile *rootHistFile = new TFile(file.c_str(), "READ");
if (!rootHistFile || rootHistFile->IsZombie()) {
NUIS_ABORT("Couldn't open root file: \"" << file << "\".");
}
TH1 *tempHist = dynamic_cast(rootHistFile->Get(name.c_str())->Clone());
if (!tempHist) {
NUIS_ABORT("Couldn't retrieve: \"" << name << "\" from root file: \""
<< file << "\".");
}
tempHist->SetDirectory(0);
rootHistFile->Close();
delete rootHistFile;
return tempHist;
}
TGraph *PlotUtils::GetTGraphFromRootFile(std::string file, std::string name) {
if (name.empty()) {
std::vector tempfile = GeneralUtils::ParseToStr(file, ";");
file = tempfile[0];
name = tempfile[1];
}
TDirectory *olddir = gDirectory;
TFile *rootHistFile = new TFile(file.c_str(), "READ");
if (!rootHistFile || rootHistFile->IsZombie()) {
NUIS_ABORT("Couldn't open root file: \"" << file << "\".");
}
TDirectory *newdir = gDirectory;
TGraph *temp =
dynamic_cast(rootHistFile->Get(name.c_str())->Clone());
if (!temp) {
NUIS_ABORT("Couldn't retrieve: \"" << name << "\" from root file: \""
<< file << "\".");
}
newdir->Remove(temp);
olddir->Append(temp);
rootHistFile->Close();
olddir->cd();
return temp;
}
/// Returns a vector of named TH1*s found in a single input file.
///
/// Expects a descriptor like: file.root[hist1|hist2|...]
std::vector
PlotUtils::GetTH1sFromRootFile(std::string const &descriptor) {
std::vector descriptors =
GeneralUtils::ParseToStr(descriptor, ",");
std::vector hists;
for (size_t d_it = 0; d_it < descriptors.size(); ++d_it) {
std::string &d = descriptors[d_it];
std::vector fname = GeneralUtils::ParseToStr(d, "[");
if (!fname.size() || !fname[0].length()) {
NUIS_ABORT("Couldn't find input file when attempting to parse : \""
<< d << "\". Expected input.root[hist1|hist2|...].");
}
if (fname[1][fname[1].length() - 1] == ']') {
fname[1] = fname[1].substr(0, fname[1].length() - 1);
}
std::vector histnames =
GeneralUtils::ParseToStr(fname[1], "|");
if (!histnames.size()) {
NUIS_ABORT(
"Couldn't find any histogram name specifiers when attempting to "
"parse "
": \""
<< fname[1] << "\". Expected hist1|hist2|...");
}
TFile *rootHistFile = new TFile(fname[0].c_str(), "READ");
if (!rootHistFile || rootHistFile->IsZombie()) {
NUIS_ABORT("Couldn't open root file: \"" << fname[0] << "\".");
}
for (size_t i = 0; i < histnames.size(); ++i) {
TH1 *tempHist =
dynamic_cast(rootHistFile->Get(histnames[i].c_str())->Clone());
if (!tempHist) {
NUIS_ABORT("Couldn't retrieve: \""
<< histnames[i] << "\" from root file: \"" << fname[0]
<< "\".");
}
tempHist->SetDirectory(0);
hists.push_back(tempHist);
}
rootHistFile->Close();
}
return hists;
}
// Create an array from an input file
std::vector PlotUtils::GetArrayFromTextFile(std::string DataFile) {
std::string line;
std::ifstream data(DataFile.c_str(), std::ifstream::in);
// Get first line
std::getline(data >> std::ws, line, '\n');
// Convert from a string into a vector of double
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
return entries;
}
// Get a 2D array from a text file
-std::vector>
+std::vector >
PlotUtils::Get2DArrayFromTextFile(std::string DataFile) {
std::string line;
- std::vector> DataArray;
+ std::vector > DataArray;
std::ifstream data(DataFile.c_str(), std::ifstream::in);
while (std::getline(data >> std::ws, line, '\n')) {
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
DataArray.push_back(entries);
}
return DataArray;
}
TH2D *PlotUtils::GetTH2DFromTextFile(std::string data, std::string binx,
std::string biny) {
// First read in the binning
// Array of x binning
std::vector xbins = GetArrayFromTextFile(binx);
// Array of y binning
std::vector ybins = GetArrayFromTextFile(biny);
// Read in the data
- std::vector> Data = Get2DArrayFromTextFile(data);
+ std::vector > Data = Get2DArrayFromTextFile(data);
// And finally fill the data
TH2D *DataPlot = new TH2D("TempHist", "TempHist", xbins.size() - 1, &xbins[0],
ybins.size() - 1, &ybins[0]);
int nBinsX = 0;
int nBinsY = 0;
- for (std::vector>::iterator it = Data.begin();
+ for (std::vector >::iterator it = Data.begin();
it != Data.end(); ++it) {
nBinsX++;
// Get the inner vector
std::vector temp = *it;
// Save the previous number[of bins to make sure it's uniform binning
int oldBinsY = nBinsY;
// Reset the counter
nBinsY = 0;
for (std::vector::iterator jt = temp.begin(); jt != temp.end();
++jt) {
nBinsY++;
DataPlot->SetBinContent(nBinsX, nBinsY, *jt);
DataPlot->SetBinError(nBinsX, nBinsY, 0.0);
}
if (oldBinsY > 0 && oldBinsY != nBinsY) {
NUIS_ERR(FTL, "Found non-uniform y-binning in " << data);
NUIS_ERR(FTL, "Previous slice: " << oldBinsY);
NUIS_ERR(FTL, "Current slice: " << nBinsY);
NUIS_ABORT("Non-uniform binning is not supported in "
"PlotUtils::GetTH2DFromTextFile");
}
}
// Check x bins
if (size_t(nBinsX + 1) != xbins.size()) {
NUIS_ERR(FTL,
"Number of x bins in data histogram does not match the binning "
"histogram!");
NUIS_ERR(
FTL,
"Are they the wrong way around (i.e. xbinning should be ybinning)?");
NUIS_ERR(FTL, "Data: " << nBinsX);
NUIS_ABORT("From " << binx << " binning: " << xbins.size());
}
// Check y bins
if (size_t(nBinsY + 1) != ybins.size()) {
NUIS_ERR(FTL,
"Number of y bins in data histogram does not match the binning "
"histogram!");
NUIS_ERR(
FTL,
"Are they the wrong way around (i.e. xbinning should be ybinning)?");
NUIS_ERR(FTL, "Data: " << nBinsY);
NUIS_ABORT("From " << biny << " binning: " << ybins.size());
}
return DataPlot;
}
TH1D *PlotUtils::GetSliceY(TH2D *Hist, int SliceNo) {
TH1D *Slice = Hist->ProjectionX(Form("%s_SLICEY%i", Hist->GetName(), SliceNo),
SliceNo, SliceNo, "e");
Slice->SetTitle(Form("%s, %.2f-%.2f", Hist->GetYaxis()->GetTitle(),
Hist->GetYaxis()->GetBinLowEdge(SliceNo),
Hist->GetYaxis()->GetBinLowEdge(SliceNo + 1)));
Slice->GetYaxis()->SetTitle(Hist->GetZaxis()->GetTitle());
return Slice;
}
TH1D *PlotUtils::GetSliceX(TH2D *Hist, int SliceNo) {
TH1D *Slice = Hist->ProjectionY(Form("%s_SLICEX%i", Hist->GetName(), SliceNo),
SliceNo, SliceNo, "e");
Slice->SetTitle(Form("%s, %.2f-%.2f", Hist->GetXaxis()->GetTitle(),
Hist->GetXaxis()->GetBinLowEdge(SliceNo),
Hist->GetXaxis()->GetBinLowEdge(SliceNo + 1)));
Slice->GetYaxis()->SetTitle(Hist->GetZaxis()->GetTitle());
return Slice;
}
void PlotUtils::AddNeutModeArray(TH1D *hist1[], TH1D *hist2[], double scaling) {
for (int i = 0; i < 60; i++) {
if (!hist2[i])
continue;
if (!hist1[i])
continue;
hist1[i]->Add(hist2[i], scaling);
}
return;
}
void PlotUtils::ScaleToData(TH1D *data, TH1D *mc, TH1I *mask) {
double scaleF = GetDataMCRatio(data, mc, mask);
mc->Scale(scaleF);
return;
}
void PlotUtils::MaskBins(TH1D *hist, TH1I *mask) {
for (int i = 0; i < hist->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1) <= 0.5)
continue;
hist->SetBinContent(i + 1, 0.0);
hist->SetBinError(i + 1, 0.0);
NUIS_LOG(DEB, "MaskBins: Set " << hist->GetName() << " Bin " << i + 1
<< " to 0.0 +- 0.0");
}
return;
}
void PlotUtils::MaskBins(TH2D *hist, TH2I *mask) {
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (mask->GetBinContent(i + 1, j + 1) <= 0.5)
continue;
hist->SetBinContent(i + 1, j + 1, 0.0);
hist->SetBinError(i + 1, j + 1, 0.0);
NUIS_LOG(DEB, "MaskBins: Set " << hist->GetName() << " Bin " << i + 1
<< " " << j + 1 << " to 0.0 +- 0.0");
}
}
return;
}
double PlotUtils::GetDataMCRatio(TH1D *data, TH1D *mc, TH1I *mask) {
double rat = 1.0;
TH1D *newmc = (TH1D *)mc->Clone();
TH1D *newdt = (TH1D *)data->Clone();
if (mask) {
MaskBins(newmc, mask);
MaskBins(newdt, mask);
}
rat = newdt->Integral() / newmc->Integral();
return rat;
}
TH2D *PlotUtils::GetCorrelationPlot(TH2D *cov, std::string name) {
TH2D *cor = (TH2D *)cov->Clone();
cor->Reset();
for (int i = 0; i < cov->GetNbinsX(); i++) {
for (int j = 0; j < cov->GetNbinsY(); j++) {
if (cov->GetBinContent(i + 1, i + 1) != 0.0 and
cov->GetBinContent(j + 1, j + 1) != 0.0)
cor->SetBinContent(i + 1, j + 1,
cov->GetBinContent(i + 1, j + 1) /
(sqrt(cov->GetBinContent(i + 1, i + 1) *
cov->GetBinContent(j + 1, j + 1))));
}
}
if (!name.empty()) {
cor->SetNameTitle(name.c_str(), (name + ";;correlation").c_str());
}
cor->SetMinimum(-1);
cor->SetMaximum(1);
return cor;
}
TH2D *PlotUtils::GetDecompPlot(TH2D *cov, std::string name) {
TMatrixDSym *covarmat = new TMatrixDSym(cov->GetNbinsX());
for (int i = 0; i < cov->GetNbinsX(); i++)
for (int j = 0; j < cov->GetNbinsY(); j++)
(*covarmat)(i, j) = cov->GetBinContent(i + 1, j + 1);
TMatrixDSym *decompmat = StatUtils::GetDecomp(covarmat);
TH2D *dec = (TH2D *)cov->Clone();
for (int i = 0; i < cov->GetNbinsX(); i++)
for (int j = 0; j < cov->GetNbinsY(); j++)
dec->SetBinContent(i + 1, j + 1, (*decompmat)(i, j));
delete covarmat;
delete decompmat;
dec->SetNameTitle(name.c_str(), (name + ";;;decomposition").c_str());
return dec;
}
TH2D *PlotUtils::MergeIntoTH2D(TH1D *xhist, TH1D *yhist, std::string zname) {
std::vector xedges, yedges;
for (int i = 0; i < xhist->GetNbinsX() + 2; i++) {
xedges.push_back(xhist->GetXaxis()->GetBinLowEdge(i + 1));
}
for (int i = 0; i < yhist->GetNbinsX() + 2; i++) {
yedges.push_back(yhist->GetXaxis()->GetBinLowEdge(i + 1));
}
int nbinsx = xhist->GetNbinsX();
int nbinsy = yhist->GetNbinsX();
std::string name =
std::string(xhist->GetName()) + "_vs_" + std::string(yhist->GetName());
std::string titles = ";" + std::string(xhist->GetXaxis()->GetTitle()) + ";" +
std::string(yhist->GetXaxis()->GetTitle()) + ";" + zname;
TH2D *newplot = new TH2D(name.c_str(), (name + titles).c_str(), nbinsx,
&xedges[0], nbinsy, &yedges[0]);
return newplot;
}
//***************************************************
void PlotUtils::MatchEmptyBins(TH1D *data, TH1D *mc) {
//**************************************************
for (int i = 0; i < data->GetNbinsX(); i++) {
if (data->GetBinContent(i + 1) == 0.0 or data->GetBinError(i + 1) == 0.0)
mc->SetBinContent(i + 1, 0.0);
}
return;
}
//***************************************************
void PlotUtils::MatchEmptyBins(TH2D *data, TH2D *mc) {
//**************************************************
for (int i = 0; i < data->GetNbinsX(); i++) {
for (int j = 0; j < data->GetNbinsY(); j++) {
if (data->GetBinContent(i + 1, j + 1) == 0.0 or
data->GetBinError(i + 1, j + 1) == 0.0)
mc->SetBinContent(i + 1, j + 1, 0.0);
}
}
return;
}
//***************************************************
TH1D *PlotUtils::GetProjectionX(TH2D *hist, TH2I *mask) {
//***************************************************
TH2D *maskedhist = StatUtils::ApplyHistogramMasking(hist, mask);
// This includes the underflow/overflow
TH1D *hist_X =
maskedhist->ProjectionX("_px", 1, maskedhist->GetXaxis()->GetNbins());
hist_X->SetTitle(Form("%s x no under/overflow", hist_X->GetTitle()));
delete maskedhist;
return hist_X;
}
//***************************************************
TH1D *PlotUtils::GetProjectionY(TH2D *hist, TH2I *mask) {
//***************************************************
TH2D *maskedhist = StatUtils::ApplyHistogramMasking(hist, mask);
// This includes the underflow/overflow
TH1D *hist_Y =
maskedhist->ProjectionY("_py", 1, maskedhist->GetYaxis()->GetNbins());
hist_Y->SetTitle(Form("%s y no under/overflow", hist_Y->GetTitle()));
delete maskedhist;
return hist_Y;
}