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() << "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; }