diff --git a/.github/CIscripts/RAM.sh b/.github/CIscripts/RAM.sh index a0385c3cb..bd254bf72 100755 --- a/.github/CIscripts/RAM.sh +++ b/.github/CIscripts/RAM.sh @@ -16,7 +16,7 @@ source bin/setup.NuOscillator.sh NUM_CORES=$(nproc) # Start the executable in the background -./bin/MCMCTutorial Inputs/FitterConfig.yaml & +./bin/MCMCTutorial Inputs/FitterConfig.yaml General:MCMC:NSteps:100000 & # Get the PID of the process PID=$! diff --git a/.github/workflows/CIBuild.yml b/.github/workflows/CIBuild.yml index 81109e679..809da7196 100644 --- a/.github/workflows/CIBuild.yml +++ b/.github/workflows/CIBuild.yml @@ -34,6 +34,10 @@ jobs: file: Doc/MaCh3DockerFiles/Rocky9/Dockerfile tag: rocky9latest cmakeoptions: -DMaCh3_PYTHON_ENABLED=ON + - os: Alma9 NoMultithreading + file: Doc/MaCh3DockerFiles/Alma9/Dockerfile + tag: alma9latest + cmakeoptions: -DMaCh3_PYTHON_ENABLED=ON -DMaCh3_MULTITHREAD_ENABLED=OFF name: Build CI ${{ matrix.os }} @@ -44,4 +48,4 @@ jobs: run: echo "${{ secrets.GITHUB_TOKEN }}" | docker login ghcr.io -u ${{ github.actor }} --password-stdin - name: Build the Docker image - run: docker build . --file ${{ matrix.file }} --tag ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag }} --build-arg MACH3_VERSION=${{ github.head_ref }} --build-arg CMAKE_OPTIONS="${{ matrix.cmakeoptions }}" + run: docker build . --file ${{ matrix.file }} --tag ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag }} --build-arg MACH3_VERSION="${{ github.head_ref }}" --build-arg CMAKE_OPTIONS="${{ matrix.cmakeoptions }}" diff --git a/.github/workflows/README.md b/.github/workflows/README.md new file mode 100644 index 000000000..85cd53fa2 --- /dev/null +++ b/.github/workflows/README.md @@ -0,0 +1,51 @@ +# Actions + +Summary of currently implemented bots + +## Benchmark +* Performs several action like reweight and save output in nice chart + +## CDImage +* Creates container after each commit to develop or after tag + +## CIBuild +* Checks compilation on several resources + +## CIPythonValidations +* Testin of pyMaCh3 + +## CIValidations +* Unit and integration test of MaCh3 + +## CodeQL +* Code scanning of MaCh3 + +## Doxygen +* Creates doxygen and spinhx documentation for MaCh3 + +## Label +* Label each PR based on modified files + +## Linter +* Lint code base (currently C++ turned off) + +## MakeRelase +* Make realse after tag is made + +## Meme +* Nothing to add, cool meme for PR + +## Newsletter +* Make weekly summary of + +## PRtittleChecker +* Ensures each PR starts with suffix like tidy: for example + +## Stale +* Tags stale PRs and issues + +## TagTutorial +* Make tag of tutorial after MaCh3 is tagged to ensure both have corresponding tags + +## Telemetry +* Add telemetry like RAM and CPU usage for each PR diff --git a/CMakeLists.txt b/CMakeLists.txt index 49ecf45ac..4618e1343 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ set(CMAKE_VERBOSE_MAKEFILE ON) # CMake version check cmake_minimum_required(VERSION 3.14 FATAL_ERROR) -project(MaCh3 VERSION 1.3.4 LANGUAGES CXX) +project(MaCh3 VERSION 1.4.0 LANGUAGES CXX) set(MaCh3_VERSION ${PROJECT_VERSION}) #LP - This option name is confusing, but I wont change it now. diff --git a/Diagnostics/CombineMaCh3Chains.cpp b/Diagnostics/CombineMaCh3Chains.cpp index b104bbda1..571a3c30c 100644 --- a/Diagnostics/CombineMaCh3Chains.cpp +++ b/Diagnostics/CombineMaCh3Chains.cpp @@ -1,13 +1,11 @@ -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" // C++ includes #include +// MaCh3 includes +#include "manager/manager.h" +#include "samplePDF/Structs.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT includes #include "TList.h" #include "TFile.h" @@ -17,10 +15,7 @@ #include "TFileMerger.h" #include "TKey.h" #include "TROOT.h" -#pragma GCC diagnostic pop - -// MaCh3 includes -#include "manager/manager.h" +_MaCh3_Safe_Include_End_ //} /// @file CombineMaCh3Chains /// @author Ewan Miller diff --git a/Diagnostics/DiagMCMC.cpp b/Diagnostics/DiagMCMC.cpp index 0b34744e6..4ce1c4d52 100644 --- a/Diagnostics/DiagMCMC.cpp +++ b/Diagnostics/DiagMCMC.cpp @@ -9,7 +9,7 @@ void DiagMCMC(const std::string& inputFile, const std::string& config) { MACH3LOG_INFO("File for study: {}", inputFile); - YAML::Node Settings = YAML::LoadFile(config); + YAML::Node Settings = M3OpenConfig(config); // Make the processor auto Processor = std::make_unique(inputFile); diff --git a/Diagnostics/GetPenaltyTerm.cpp b/Diagnostics/GetPenaltyTerm.cpp index 85b60f839..547b764b6 100644 --- a/Diagnostics/GetPenaltyTerm.cpp +++ b/Diagnostics/GetPenaltyTerm.cpp @@ -1,10 +1,8 @@ -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +// MaCh3 includes +#include "manager/manager.h" +#include "samplePDF/Structs.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT includes #include "TFile.h" #include "TBranch.h" @@ -22,14 +20,7 @@ #include "TColor.h" #include "TObjString.h" #include "TROOT.h" -#pragma GCC diagnostic pop - -#ifdef MULTITHREAD -#include "omp.h" -#endif - -#include "manager/manager.h" - +_MaCh3_Safe_Include_End_ //} /// @file GetPenaltyTerm.cpp /// @brief KS: This file contains the implementation of the function to extract specific penalty terms from systematic chains. @@ -85,7 +76,7 @@ void ReadXSecFile(const std::string& inputFile) XSecFile["Systematics"] = YAML::Node(YAML::NodeType::Sequence); for(unsigned int i = 0; i < XsecCovPos.size(); i++) { - YAML::Node YAMLDocTemp = YAML::LoadFile(XsecCovPos[i]); + YAML::Node YAMLDocTemp = M3OpenConfig(XsecCovPos[i]); for (const auto& item : YAMLDocTemp["Systematics"]) { XSecFile["Systematics"].push_back(item); } diff --git a/Diagnostics/ProcessMCMC.cpp b/Diagnostics/ProcessMCMC.cpp index 0fe320a95..1b544ca5b 100644 --- a/Diagnostics/ProcessMCMC.cpp +++ b/Diagnostics/ProcessMCMC.cpp @@ -74,7 +74,7 @@ void ProcessMCMC(const std::string& inputFile) // Make the processor) auto Processor = std::make_unique(inputFile); - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); + YAML::Node card_yaml = M3OpenConfig(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; const bool PlotCorr = GetFromManager(Settings["PlotCorr"], false); @@ -152,7 +152,7 @@ void ProcessMCMC(const std::string& inputFile) void MultipleProcessMCMC() { - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); + YAML::Node card_yaml = M3OpenConfig(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; constexpr Color_t PosteriorColor[] = {kBlue-1, kRed, kGreen+2}; @@ -318,7 +318,7 @@ void MultipleProcessMCMC() // KS: Calculate Bayes factor for a given hypothesis, most informative are those related to osc params. However, it make relative easy interpretation for switch dials void CalcBayesFactor(const std::unique_ptr& Processor) { - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); + YAML::Node card_yaml = M3OpenConfig(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; std::vector ParNames; @@ -338,7 +338,7 @@ void CalcBayesFactor(const std::unique_ptr& Processor) void CalcSavageDickey(const std::unique_ptr& Processor) { - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); + YAML::Node card_yaml = M3OpenConfig(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; std::vector ParNames; @@ -356,7 +356,7 @@ void CalcSavageDickey(const std::unique_ptr& Processor) void CalcParameterEvolution(const std::unique_ptr& Processor) { - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); + YAML::Node card_yaml = M3OpenConfig(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; std::vector ParNames; @@ -371,7 +371,7 @@ void CalcParameterEvolution(const std::unique_ptr& Processor) void CalcBipolarPlot(const std::unique_ptr& Processor) { - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); + YAML::Node card_yaml = M3OpenConfig(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; std::vector ParNames; @@ -383,7 +383,7 @@ void CalcBipolarPlot(const std::unique_ptr& Processor) } void GetTrianglePlot(const std::unique_ptr& Processor) { - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); + YAML::Node card_yaml = M3OpenConfig(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; for (const auto& dg : Settings["TrianglePlot"]) @@ -433,7 +433,7 @@ void DiagnoseCovarianceMatrix(const std::unique_ptr& Processor, c Canvas->Print(Form("Correlation_%s.pdf[", OutName.c_str()), "pdf"); Canvas->Print(Form("Covariance_%s.pdf[", OutName.c_str()), "pdf"); - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); + YAML::Node card_yaml = M3OpenConfig(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; const int entries = int(Processor->GetnSteps()); @@ -579,7 +579,7 @@ void DiagnoseCovarianceMatrix(const std::unique_ptr& Processor, c void ReweightPrior(const std::unique_ptr& Processor) { - YAML::Node card_yaml = YAML::LoadFile(config.c_str()); + YAML::Node card_yaml = M3OpenConfig(config.c_str()); YAML::Node Settings = card_yaml["ProcessMCMC"]; const auto& Prior = Settings["PriorReweighting"]; diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp index 8103b56e2..7ee4f98f8 100644 --- a/Diagnostics/RHat.cpp +++ b/Diagnostics/RHat.cpp @@ -1,10 +1,8 @@ -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +// MaCh3 includes +#include "manager/manager.h" +#include "samplePDF/Structs.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT includes #include "TObjArray.h" #include "TChain.h" @@ -20,14 +18,7 @@ #include "TColor.h" #include "TStyle.h" #include "TROOT.h" -#pragma GCC diagnostic pop - -#ifdef MULTITHREAD -#include "omp.h" -#endif - -// MaCh3 includes -#include "manager/manager.h" +_MaCh3_Safe_Include_End_ //} /// @file RHat.cpp /// @brief This executable calculates the \f$ \hat{R} \f$ estimator for Markov Chain Monte Carlo (MCMC) convergence. diff --git a/Diagnostics/RHat_HighMem.cpp b/Diagnostics/RHat_HighMem.cpp index f27e7c8ec..48ffed571 100644 --- a/Diagnostics/RHat_HighMem.cpp +++ b/Diagnostics/RHat_HighMem.cpp @@ -1,10 +1,8 @@ -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +// MaCh3 includes +#include "manager/manager.h" +#include "samplePDF/Structs.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT includes #include "TObjArray.h" #include "TChain.h" @@ -20,14 +18,7 @@ #include "TColor.h" #include "TStyle.h" #include "TROOT.h" -#pragma GCC diagnostic pop - -#ifdef MULTITHREAD -#include "omp.h" -#endif - -// MaCh3 includes -#include "manager/manager.h" +_MaCh3_Safe_Include_End_ //} /// @file RHat.cpp /// @brief This executable calculates the \f$ \hat{R} \f$ estimator for Markov Chain Monte Carlo (MCMC) convergence. diff --git a/Doc/Doxyfile b/Doc/Doxyfile index 4d41d8071..b15bc2855 100644 --- a/Doc/Doxyfile +++ b/Doc/Doxyfile @@ -38,7 +38,7 @@ PROJECT_NAME = "MaCh3" # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = 1.3.4 +PROJECT_NUMBER = 1.4.0 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/Doc/MaCh3DockerFiles/Rocky9/Dockerfile b/Doc/MaCh3DockerFiles/Rocky9/Dockerfile index ae38be7c7..880b18f6c 100644 --- a/Doc/MaCh3DockerFiles/Rocky9/Dockerfile +++ b/Doc/MaCh3DockerFiles/Rocky9/Dockerfile @@ -5,6 +5,8 @@ FROM kamilskwarczynski/nukamil:latest AS mach3_build # Add a label for the author LABEL maintainer="The MaCh3 Collaboration" LABEL website="https://mach3-software.github.io/MaCh3/" +LABEL compiler="GNU 11.4.1" +LABEL root_version="v6.30.04" LABEL org.opencontainers.image.description="Official MaCh3 container" # Declare the build argument diff --git a/README.md b/README.md index e5204d848..e54623d14 100644 --- a/README.md +++ b/README.md @@ -128,7 +128,6 @@ The following fitting algorithms are available: | MINUIT2 | [Ref](https://cds.cern.ch/record/2296388/) | Yes | | PSO | [Ref](https://doi.org/10.1162/EVCO_r_00180) | No | - ## Debug Several debugging options are available which are heavy for RAM and performance and, therefore not used by default. To enable it: ```bash @@ -160,6 +159,7 @@ Based on several test here are recommended version: | Fedora32 | ✅ | | CentOS7 | ❔ | | Windows | ❌ | +| MacOS | ❌ | ✅ - Part of CI/CD
❔ - Not part of CI/CD but used by some users/developers so it might work
diff --git a/cmake/Templates/setup.MaCh3.sh.in b/cmake/Templates/setup.MaCh3.sh.in index 48a8cf9c0..a37d78b3d 100644 --- a/cmake/Templates/setup.MaCh3.sh.in +++ b/cmake/Templates/setup.MaCh3.sh.in @@ -59,13 +59,19 @@ function add_to_PYTHONPATH () { d=$(cd -- "$d" && { pwd -P || pwd; }) 2>/dev/null # canonicalize symbolic links if [ -z "$d" ]; then continue; fi # skip nonexistent directory - case ":$PYTHONPATH:" in - *":$d:"*) :;; - *) export PYTHONPATH=$d:$PYTHONPATH;; - esac + if [ "$d" == "/usr/bin" ] || [ "$d" == "/usr/bin64" ] || [ "$d" == "/usr/local/bin" ] || [ "$d" == "/usr/local/bin64" ]; then + case ":$PYTHONPATH:" in + *":$d:"*) :;; + *) export PYTHONPATH=$PYTHONPATH:$d;; + esac + else + case ":$PYTHONPATH:" in + *":$d:"*) :;; + *) export PYTHONPATH=$d:$PYTHONPATH;; + esac + fi done } - fi diff --git a/covariance/CovarianceUtils.h b/covariance/CovarianceUtils.h index 19d5c8ca7..6b5bbf86e 100644 --- a/covariance/CovarianceUtils.h +++ b/covariance/CovarianceUtils.h @@ -1,12 +1,9 @@ #pragma once -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +// MaCh3 includes +#include "samplePDF/Structs.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT includes #include "TMatrixT.h" #include "TMatrixDSym.h" @@ -24,26 +21,10 @@ #include "TMatrixDSymEigen.h" #include "TMatrixDEigen.h" #include "TDecompSVD.h" -#pragma GCC diagnostic pop - -#ifdef MULTITHREAD -#include "omp.h" -#endif +_MaCh3_Safe_Include_End_ //} -namespace MaCh3Utils +namespace M3 { - /// @brief number of threads which we need for example for TRandom3 - inline int GetNThreads() - { - #ifdef MULTITHREAD - int nThreads = omp_get_max_threads(); - #else - int nThreads = 1; - #endif - - return nThreads; - } - /// @brief CW: Multi-threaded matrix multiplication inline double* MatrixMult(double *A, double *B, int n) { //CW: First transpose to increse cache hits diff --git a/covariance/covarianceBase.cpp b/covariance/covarianceBase.cpp index 4756ffb83..d43f0393a 100644 --- a/covariance/covarianceBase.cpp +++ b/covariance/covarianceBase.cpp @@ -1,4 +1,5 @@ #include "covariance/covarianceBase.h" +#include "regex" // ******************************************** covarianceBase::covarianceBase(std::string name, std::string file, double threshold, int FirstPCA, int LastPCA) : inputFile(file), pca(false), @@ -116,7 +117,7 @@ void covarianceBase::init(std::string name, std::string file) { PrintLength = 35; - const int nThreads = MaCh3Utils::GetNThreads(); + const int nThreads = M3::GetNThreads(); //KS: set Random numbers for each thread so each thread has different seed //or for one thread if without MULTITHREAD random_number.reserve(nThreads); @@ -171,7 +172,7 @@ void covarianceBase::init(const std::vector& YAMLFile) { } } - const int nThreads = MaCh3Utils::GetNThreads(); + const int nThreads = M3::GetNThreads(); //KS: set Random numbers for each thread so each thread has different seed //or for one thread if without MULTITHREAD random_number.reserve(nThreads); @@ -210,6 +211,7 @@ void covarianceBase::init(const std::vector& YAMLFile) { _fGenerated[i] = Get(param["Systematic"]["ParameterValues"]["Generated"], __FILE__ , __LINE__); _fIndivStepScale[i] = Get(param["Systematic"]["StepScale"]["MCMC"], __FILE__ , __LINE__); _fError[i] = Get(param["Systematic"]["Error"], __FILE__ , __LINE__); + _fDetID[i] = GetFromManager>(param["Systematic"]["DetID"], {}, __FILE__, __LINE__); if(_fError[i] <= 0) { MACH3LOG_ERROR("Error for param {}({}) is negative and eqaul to {}", _fFancyNames[i], i, _fError[i]); throw MaCh3Exception(__FILE__ , __LINE__ ); @@ -333,6 +335,7 @@ void covarianceBase::ReserveMemory(const int SizeVec) { _fUpBound = std::vector(SizeVec); _fFlatPrior = std::vector(SizeVec); _fIndivStepScale = std::vector(SizeVec); + _fDetID = std::vector>(_fNumPar); corr_throw = new double[SizeVec](); // set random parameter vector (for correlated steps) @@ -1302,3 +1305,35 @@ void covarianceBase::SaveUpdatedMatrixConfig() { fout << copyNode; fout.close(); } + +// ******************************************** +bool covarianceBase::AppliesToDetID(const int SystIndex, const std::string& DetID) const { +// ******************************************** + // Empty means apply to all + if (_fDetID[SystIndex].size() == 0) return true; + + // Make a copy and to lower case to not be case sensitive + std::string DetIDCopy = DetID; + std::transform(DetIDCopy.begin(), DetIDCopy.end(), DetIDCopy.begin(), ::tolower); + bool Applies = false; + + for (size_t i = 0; i < _fDetID[SystIndex].size(); i++) { + // Convert to low case to not be case sensitive + std::string pattern = _fDetID[SystIndex][i]; + std::transform(pattern.begin(), pattern.end(), pattern.begin(), ::tolower); + + // Replace '*' in the pattern with '.*' for regex matching + std::string regexPattern = "^" + std::regex_replace(pattern, std::regex("\\*"), ".*") + "$"; + try { + std::regex regex(regexPattern); + if (std::regex_match(DetIDCopy, regex)) { + Applies = true; + break; + } + } catch (const std::regex_error& e) { + // Handle regex error (for invalid patterns) + MACH3LOG_ERROR("Regex error: {}", e.what()); + } + } + return Applies; +} diff --git a/covariance/covarianceBase.h b/covariance/covarianceBase.h index 1a2d2614b..0c789a820 100644 --- a/covariance/covarianceBase.h +++ b/covariance/covarianceBase.h @@ -425,6 +425,11 @@ class covarianceBase { /// @cite haario2001adaptive void updateAdaptiveCovariance(); + /// @brief Check if parameter is affecting given det ID + /// @param SystIndex number of parameter + /// @param DetID The Detector ID used to filter parameters. + bool AppliesToDetID(const int SystIndex, const std::string& DetID) const; + /// The input root file we read in const std::string inputFile; @@ -476,6 +481,8 @@ class covarianceBase { std::vector _fIndivStepScale; /// Whether to apply flat prior or not std::vector _fFlatPrior; + /// Tells to which samples object param should be applied + std::vector> _fDetID; /// perform PCA or not bool pca; diff --git a/covariance/covarianceOsc.cpp b/covariance/covarianceOsc.cpp index f5070f830..d4d902dd4 100644 --- a/covariance/covarianceOsc.cpp +++ b/covariance/covarianceOsc.cpp @@ -81,11 +81,34 @@ void covarianceOsc::Print() { MACH3LOG_INFO("Number of pars: {}", _fNumPar); MACH3LOG_INFO("Current: {} parameters:", matrixName); - MACH3LOG_INFO("{:<5} | {:<25} | {:<10} | {:<15} | {:<15} | {:<10}", - "#", "Name", "Nom.", "IndivStepScale", "_fError", "FlatPrior"); + MACH3LOG_INFO("================================================================================================================================="); + MACH3LOG_INFO("{:<5} {:2} {:<25} {:2} {:<10} {:2} {:<15} {:2} {:<15} {:2} {:<10} {:2} {:<10}", + "#", "|", "Name", "|", "Prior", "|", "IndivStepScale", "|", "Error", "|", "FlatPrior", "|", "DetID"); + MACH3LOG_INFO("---------------------------------------------------------------------------------------------------------------------------------"); + for (int i = 0; i < _fNumPar; i++) { + std::string detIdString = ""; + for (const auto& detID : _fDetID[i]) { + if (!detIdString.empty()) { + detIdString += ", "; + } + detIdString += detID; + } - for(int i = 0; i < _fNumPar; i++) { - MACH3LOG_INFO("{:<5} | {:<25} | {:<10.4f} | {:<15.2f} | {:<15.4f} | {:<10}", - i, _fNames[i].c_str(), _fPreFitValue[i], _fIndivStepScale[i], _fError[i], _fFlatPrior[i]); + MACH3LOG_INFO("{:<5} {:2} {:<25} {:2} {:<10.4f} {:2} {:<15.2f} {:2} {:<15.4f} {:2} {:<10} {:2} {:<10}", + i, "|", _fNames[i].c_str(), "|", _fPreFitValue[i], "|", _fIndivStepScale[i], "|", _fError[i], "|", _fFlatPrior[i], "|", detIdString); + } + MACH3LOG_INFO("================================================================================================================================="); +} + +// ******************************************** +// DB Grab the Normalisation parameters for the relevant DetID +std::vector covarianceOsc::GetOscParsFromDetID(const std::string& DetID) { +// ******************************************** + std::vector returnVec; + for (int i = 0; i < _fNumPar; ++i) { + if (AppliesToDetID(i, DetID)) { + returnVec.push_back(retPointer(i)); + } } + return returnVec; } diff --git a/covariance/covarianceOsc.h b/covariance/covarianceOsc.h index ac4330c13..d56ac39af 100644 --- a/covariance/covarianceOsc.h +++ b/covariance/covarianceOsc.h @@ -17,7 +17,8 @@ class covarianceOsc : public covarianceBase void proposeStep() override; /// @brief Sets whether to flip delta M23. void setFlipDeltaM23(bool flip){flipdelM = flip;} - + /// @brief Get pointers to Osc params from detId + std::vector GetOscParsFromDetID(const std::string& DetID); /// @brief KS: Print all useful information's after initialization void Print(); diff --git a/covariance/covarianceXsec.cpp b/covariance/covarianceXsec.cpp index ca52db737..77ec410ed 100644 --- a/covariance/covarianceXsec.cpp +++ b/covariance/covarianceXsec.cpp @@ -28,7 +28,6 @@ void covarianceXsec::InitXsecFromConfig() { // ******************************************** _fSystToGlobalSystIndexMap.resize(SystType::kSystTypes); - _fDetID = std::vector(_fNumPar); _fParamType = std::vector(_fNumPar); _ParameterGroup = std::vector(_fNumPar); @@ -43,7 +42,6 @@ void covarianceXsec::InitXsecFromConfig() { //PreFitValues etc etc. for (auto const ¶m : _fYAMLDoc["Systematics"]) { - _fDetID[i] = Get(param["Systematic"]["DetID"], __FILE__ , __LINE__); _ParameterGroup[i] = Get(param["Systematic"]["ParameterGroup"], __FILE__ , __LINE__); //Fill the map to get the correlations later as well @@ -105,7 +103,7 @@ covarianceXsec::~covarianceXsec() { // ******************************************** // DB Grab the Spline Names for the relevant DetID -const std::vector covarianceXsec::GetSplineParsNamesFromDetID(const int DetID) { +const std::vector covarianceXsec::GetSplineParsNamesFromDetID(const std::string& DetID) { // ******************************************** std::vector returnVec; for (auto &pair : _fSystToGlobalSystIndexMap[SystType::kSpline]) { @@ -120,7 +118,7 @@ const std::vector covarianceXsec::GetSplineParsNamesFromDetID(const } // ******************************************** -const std::vector covarianceXsec::GetSplineInterpolationFromDetID(const int DetID) { +const std::vector covarianceXsec::GetSplineInterpolationFromDetID(const std::string& DetID) { // ******************************************** std::vector returnVec; for (auto &pair : _fSystToGlobalSystIndexMap[SystType::kSpline]) { @@ -136,7 +134,7 @@ const std::vector covarianceXsec::GetSplineInterpolationFro // ******************************************** // DB Grab the Spline Modes for the relevant DetID -const std::vector< std::vector > covarianceXsec::GetSplineModeVecFromDetID(const int DetID) { +const std::vector< std::vector > covarianceXsec::GetSplineModeVecFromDetID(const std::string& DetID) { // ******************************************** std::vector< std::vector > returnVec; //Need a counter or something to correctly get the index in _fSplineModes since it's not of length nPars @@ -144,7 +142,7 @@ const std::vector< std::vector > covarianceXsec::GetSplineModeVecFromDetID( for (auto &pair : _fSystToGlobalSystIndexMap[SystType::kSpline]) { auto &SplineIndex = pair.first; auto &SystIndex = pair.second; - if ((GetParDetID(SystIndex) & DetID)) { //If parameter applies to required DetID + if (AppliesToDetID(SystIndex, DetID)) { //If parameter applies to required DetID returnVec.push_back(SplineParams.at(SplineIndex)._fSplineModes); } } @@ -212,7 +210,7 @@ XsecNorms4 covarianceXsec::GetXsecNorm(const YAML::Node& param, const int Index) // ******************************************** // Grab the global syst index for the relevant DetID // i.e. get a vector of size nSplines where each entry is filled with the global syst number -const std::vector covarianceXsec::GetGlobalSystIndexFromDetID(const int DetID, const SystType Type) { +const std::vector covarianceXsec::GetGlobalSystIndexFromDetID(const std::string& DetID, const SystType Type) { // ******************************************** std::vector returnVec; for (auto &pair : _fSystToGlobalSystIndexMap[Type]) { @@ -227,7 +225,7 @@ const std::vector covarianceXsec::GetGlobalSystIndexFromDetID(const int Det // ******************************************** // Grab the global syst index for the relevant DetID // i.e. get a vector of size nSplines where each entry is filled with the global syst number -const std::vector covarianceXsec::GetSystIndexFromDetID(int DetID, const SystType Type) { +const std::vector covarianceXsec::GetSystIndexFromDetID(const std::string& DetID, const SystType Type) { // ******************************************** std::vector returnVec; for (auto &pair : _fSystToGlobalSystIndexMap[Type]) { @@ -266,7 +264,7 @@ XsecSplines1 covarianceXsec::GetXsecSpline(const YAML::Node& param) { // ******************************************** // DB Grab the Normalisation parameters for the relevant DetID -const std::vector covarianceXsec::GetNormParsFromDetID(const int DetID) { +const std::vector covarianceXsec::GetNormParsFromDetID(const std::string& DetID) { // ******************************************** std::vector returnVec; for (auto &pair : _fSystToGlobalSystIndexMap[SystType::kNorm]) { @@ -281,7 +279,7 @@ const std::vector covarianceXsec::GetNormParsFromDetID(const int Det // ******************************************** // DB Grab the number of parameters for the relevant DetID -int covarianceXsec::GetNumParamsFromDetID(const int DetID, const SystType Type) { +int covarianceXsec::GetNumParamsFromDetID(const std::string& DetID, const SystType Type) { // ******************************************** int returnVal = 0; IterateOverParams(DetID, @@ -293,7 +291,7 @@ int covarianceXsec::GetNumParamsFromDetID(const int DetID, const SystType Type) // ******************************************** // DB Grab the parameter names for the relevant DetID -const std::vector covarianceXsec::GetParsNamesFromDetID(const int DetID, const SystType Type) { +const std::vector covarianceXsec::GetParsNamesFromDetID(const std::string& DetID, const SystType Type) { // ******************************************** std::vector returnVec; IterateOverParams(DetID, @@ -305,7 +303,7 @@ const std::vector covarianceXsec::GetParsNamesFromDetID(const int D // ******************************************** // DB DB Grab the parameter indices for the relevant DetID -const std::vector covarianceXsec::GetParsIndexFromDetID(const int DetID, const SystType Type) { +const std::vector covarianceXsec::GetParsIndexFromDetID(const std::string& DetID, const SystType Type) { // ******************************************** std::vector returnVec; IterateOverParams(DetID, @@ -317,7 +315,7 @@ const std::vector covarianceXsec::GetParsIndexFromDetID(const int DetID, co // ******************************************** template -void covarianceXsec::IterateOverParams(const int DetID, FilterFunc filter, ActionFunc action) { +void covarianceXsec::IterateOverParams(const std::string& DetID, FilterFunc filter, ActionFunc action) { // ******************************************** for (int i = 0; i < _fNumPar; ++i) { if ((AppliesToDetID(i, DetID)) && filter(i)) { // Common filter logic @@ -382,11 +380,18 @@ void covarianceXsec::Print() { void covarianceXsec::PrintGlobablInfo() { // ******************************************** MACH3LOG_INFO("============================================================================================================================================================"); - MACH3LOG_INFO("{:<5} {:2} {:<40} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<5} {:2} {:<10}", "#", "|", "Name", "|", "Gen.", "|", "Prior", "|", "Error", "|", "Lower", "|", "Upper", "|", "StepScale", "|", "DetID", "|", "Type"); + MACH3LOG_INFO("{:<5} {:2} {:<40} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10}", "#", "|", "Name", "|", "Gen.", "|", "Prior", "|", "Error", "|", "Lower", "|", "Upper", "|", "StepScale", "|", "DetID", "|", "Type"); MACH3LOG_INFO("------------------------------------------------------------------------------------------------------------------------------------------------------------"); for (int i = 0; i < GetNumParams(); i++) { std::string ErrString = fmt::format("{:.2f}", _fError[i]); - MACH3LOG_INFO("{:<5} {:2} {:<40} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<5} {:2} {:<10}", i, "|", GetParFancyName(i), "|", _fGenerated[i], "|", _fPreFitValue[i], "|", "+/- " + ErrString, "|", _fLowBound[i], "|", _fUpBound[i], "|", _fIndivStepScale[i], "|", _fDetID[i], "|", SystType_ToString(_fParamType[i])); + std::string detIdString = ""; + for (const auto& detID : _fDetID[i]) { + if (!detIdString.empty()) { + detIdString += ", "; + } + detIdString += detID; + } + MACH3LOG_INFO("{:<5} {:2} {:<40} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10} {:2} {:<10}", i, "|", GetParFancyName(i), "|", _fGenerated[i], "|", _fPreFitValue[i], "|", "+/- " + ErrString, "|", _fLowBound[i], "|", _fUpBound[i], "|", _fIndivStepScale[i], "|", detIdString, "|", SystType_ToString(_fParamType[i])); } MACH3LOG_INFO("============================================================================================================================================================"); } @@ -589,9 +594,7 @@ void covarianceXsec::DumpMatrixToFile(const std::string& Name) { TVectorD* xsec_param_knot_weight_lb = new TVectorD(_fNumPar); TVectorD* xsec_param_knot_weight_ub = new TVectorD(_fNumPar); - TVectorD* xsec_error = new TVectorD(_fNumPar); - TVectorD* xsec_param_id = new TVectorD(_fNumPar); for(int i = 0; i < _fNumPar; ++i) { @@ -609,7 +612,6 @@ void covarianceXsec::DumpMatrixToFile(const std::string& Name) { (*xsec_flat_prior)[i] = _fFlatPrior[i]; (*xsec_stepscale)[i] = _fIndivStepScale[i]; (*xsec_error)[i] = _fError[i]; - (*xsec_param_id)[i] = _fDetID[i]; (*xsec_param_lb)[i] = _fLowBound[i]; (*xsec_param_ub)[i] = _fUpBound[i]; @@ -656,9 +658,6 @@ void covarianceXsec::DumpMatrixToFile(const std::string& Name) { delete xsec_param_knot_weight_lb; xsec_param_knot_weight_ub->Write("xsec_param_knot_weight_ub"); delete xsec_param_knot_weight_ub; - - xsec_param_id->Write("xsec_param_id"); - delete xsec_param_id; xsec_error->Write("xsec_error"); delete xsec_error; diff --git a/covariance/covarianceXsec.h b/covariance/covarianceXsec.h index 47fbf8531..999313aff 100644 --- a/covariance/covarianceXsec.h +++ b/covariance/covarianceXsec.h @@ -24,7 +24,7 @@ class covarianceXsec : public covarianceBase { // General Getter functions not split by detector /// @brief ETA - just return the int of the DetID, this can be removed to do a string comp at some point. /// @param i parameter index - inline int GetParDetID(const int i) const { return _fDetID[i];}; + inline std::vector GetParDetID(const int i) const { return _fDetID[i];}; /// @brief ETA - just return a string of "spline", "norm" or "functional" /// @param i parameter index inline std::string GetParamTypeString(const int i) const { return SystType_ToString(_fParamType[i]); } @@ -36,13 +36,13 @@ class covarianceXsec : public covarianceBase { /// @param i spline parameter index, not confuse with global index inline SplineInterpolation GetParSplineInterpolation(const int i) {return SplineParams.at(i)._SplineInterpolationType;} /// @brief Get the interpolation types for splines affecting a particular DetID - const std::vector GetSplineInterpolationFromDetID(int DetID); + const std::vector GetSplineInterpolationFromDetID(const std::string& DetID); /// @brief Get the name of the spline associated with the spline at index i /// @param i spline parameter index, not to be confused with global index std::string GetParSplineName(const int i) {return _fSplineNames[i];} - //DB Get spline parameters depending on given DetID - const std::vector GetGlobalSystIndexFromDetID(const int DetID, const SystType Type); + /// @brief DB Get spline parameters depending on given DetID + const std::vector GetGlobalSystIndexFromDetID(const std::string& DetID, const SystType Type); /// @brief EM: value at which we cap spline knot weight /// @param i spline parameter index, not confuse with global index inline double GetParSplineKnotUpperBound(const int i) {return SplineParams.at(i)._SplineKnotUpBound;} @@ -52,26 +52,26 @@ class covarianceXsec : public covarianceBase { /// @brief DB Grab the number of parameters for the relevant DetID /// @param Type Type of syst, for example kNorm, kSpline etc - int GetNumParamsFromDetID(const int DetID, const SystType Type); + int GetNumParamsFromDetID(const std::string& DetID, const SystType Type); /// @brief DB Grab the parameter names for the relevant DetID /// @param Type Type of syst, for example kNorm, kSpline etc - const std::vector GetParsNamesFromDetID(const int DetID, const SystType Type); + const std::vector GetParsNamesFromDetID(const std::string& DetID, const SystType Type); /// @brief DB Grab the parameter indices for the relevant DetID /// @param Type Type of syst, for example kNorm, kSpline etc - const std::vector GetParsIndexFromDetID(const int DetID, const SystType Type); + const std::vector GetParsIndexFromDetID(const std::string& DetID, const SystType Type); /// @brief DB Get spline parameters depending on given DetID - const std::vector GetSplineParsNamesFromDetID(const int DetID); + const std::vector GetSplineParsNamesFromDetID(const std::string& DetID); /// @brief DB Get spline parameters depending on given DetID - const std::vector GetSplineFileParsNamesFromDetID(const int DetID); + const std::vector GetSplineFileParsNamesFromDetID(const std::string& DetID); /// @brief DB Grab the Spline Modes for the relevant DetID - const std::vector< std::vector > GetSplineModeVecFromDetID(const int DetID); + const std::vector< std::vector > GetSplineModeVecFromDetID(const std::string& DetID); /// @brief Grab the index of the syst relative to global numbering. /// @param Type Type of syst, for example kNorm, kSpline etc - const std::vector GetSystIndexFromDetID(const int DetID, const SystType Type); + const std::vector GetSystIndexFromDetID(const std::string& DetID, const SystType Type); /// @brief DB Get norm/func parameters depending on given DetID - const std::vector GetNormParsFromDetID(const int DetID); + const std::vector GetNormParsFromDetID(const std::string& DetID); /// @brief KS: For most covariances prior and fparInit (prior) are the same, however for Xsec those can be different std::vector getNominalArray() override @@ -126,13 +126,8 @@ class covarianceXsec : public covarianceBase { /// parameter. /// @param DetID The Detector ID used to filter parameters. template - void IterateOverParams(const int DetID, FilterFunc filter, ActionFunc action); - /// @brief Check if parameter is affecting given det ID - /// @param SystIndex number of parameter - /// @param DetID The Detector ID used to filter parameters. - bool AppliesToDetID(const int SystIndex, const int DetID) const { - return (GetParDetID(SystIndex) & DetID) != 0; - } + void IterateOverParams(const std::string& DetID, FilterFunc filter, ActionFunc action); + /// @brief Initializes the systematic parameters from the configuration file. /// This function loads parameters like normalizations and splines from the provided YAML file. /// @note This is used internally during the object's initialization process. @@ -148,8 +143,6 @@ class covarianceXsec : public covarianceBase { /// @param param Yaml node describing param inline XsecSplines1 GetXsecSpline(const YAML::Node& param); - /// Tells to which samples object param should be applied - std::vector _fDetID; /// Type of parameter like norm, spline etc. std::vector _fParamType; diff --git a/manager/CMakeLists.txt b/manager/CMakeLists.txt index 7e4058af5..765bd3268 100644 --- a/manager/CMakeLists.txt +++ b/manager/CMakeLists.txt @@ -6,6 +6,7 @@ set(HEADERS Monitor.h MaCh3Exception.h gpuUtils.cuh + Core.h ) add_library(Manager SHARED diff --git a/manager/Core.h b/manager/Core.h new file mode 100755 index 000000000..7b38677e9 --- /dev/null +++ b/manager/Core.h @@ -0,0 +1,90 @@ +#pragma once + +/// @file Core.h +/// @author Clarence Wret +/// @author Kamil Skwarczynski +/// @author Luke Pickering + +/// Run low or high memory versions of structs +/// N.B. for 64 bit systems sizeof(float) == sizeof(double) so not a huge effect +namespace M3 { + #ifdef _LOW_MEMORY_STRUCTS_ + /// Custom floating point (float or double) + using float_t = float; + /// Custom integer (int or short int) + using int_t = short; + /// Custom unsigned integer (unsigned short int or unsigned int) + using uint_t = unsigned short; + #else + using float_t = double; + using int_t = int; + using uint_t = unsigned; + #endif + + /// Function template for fused multiply-add + template + constexpr T fmaf_t(T x, T y, T z) { + #ifdef _LOW_MEMORY_STRUCTS_ + return std::fmaf(x, y, z); + #else + return std::fma(x, y, z); + #endif + } +} + +/// KS: noexcept can help with performance but is terrible for debugging, this is meant to help easy way of of turning it on or off. In near future move this to struct or other central class. +//#define SafeException +#ifndef SafeException +#define _noexcept_ noexcept +#else +#define _noexcept_ +#endif + +/// KS: Using restrict limits the effects of pointer aliasing, aiding optimizations. While reading I found that there might be some compilers which don't have __restrict__. As always we use _restrict_ to more easily turn off restrict in the code +#ifndef DEBUG +#define _restrict_ __restrict__ +#else +#define _restrict_ +#endif + +/// Number of overflow bins in TH2Poly, +#define _TH2PolyOverflowBins_ 9 + +constexpr static const double _BAD_DOUBLE_ = -999.99; +constexpr static const int _BAD_INT_ = -999; + +constexpr static const double _DEFAULT_RETURN_VAL_ = -999999.123456; + +// Some commonly used variables to which we set pointers to +constexpr static const double Unity = 1.; +constexpr static const double Zero = 0.; +constexpr static const float Unity_F = 1.; +constexpr static const float Zero_F = 0.; +constexpr static const int Unity_Int = 1; + +// C++ includes +#include +#include +#include +#include +#include + +#ifdef MULTITHREAD +#include "omp.h" +#endif + +/// @brief KS: Avoiding warning checking for headers +/// @details Many external files don't strictly adhere to rigorous C++ standards. +/// Inline functions in such headers may cause errors when compiled with strict MaCh3 compiler flags. +/// @warning Use this for any external header file to avoid warnings. +#define _MaCh3_Safe_Include_Start_ \ +_Pragma("GCC diagnostic push") \ +_Pragma("GCC diagnostic ignored \"-Wuseless-cast\"") \ +_Pragma("GCC diagnostic ignored \"-Wfloat-conversion\"") \ +_Pragma("GCC diagnostic ignored \"-Wold-style-cast\"") \ +_Pragma("GCC diagnostic ignored \"-Wformat-nonliteral\"") \ +_Pragma("GCC diagnostic ignored \"-Wconversion\"") + +/// @brief KS: Restore warning checking after including external headers +#define _MaCh3_Safe_Include_End_ \ +_Pragma("GCC diagnostic pop") diff --git a/manager/MaCh3Exception.h b/manager/MaCh3Exception.h index 637be3e35..effb75591 100644 --- a/manager/MaCh3Exception.h +++ b/manager/MaCh3Exception.h @@ -7,6 +7,7 @@ // MaCh3 Includes #include "manager/MaCh3Logger.h" +#include "manager/Core.h" /// @brief Custom exception class for MaCh3 errors. class MaCh3Exception : public std::exception { diff --git a/manager/MaCh3Modes.cpp b/manager/MaCh3Modes.cpp index fef23a63b..863596f85 100644 --- a/manager/MaCh3Modes.cpp +++ b/manager/MaCh3Modes.cpp @@ -6,7 +6,7 @@ MaCh3Modes::MaCh3Modes(std::string const &filename) { // ******************* // Load config - YAML::Node config = YAML::LoadFile(filename); + YAML::Node config = M3OpenConfig(filename); std::string GetMaCh3ModeName(const int Index); NModes = 0; diff --git a/manager/Monitor.cpp b/manager/Monitor.cpp index fc2dbf34e..1b19e7f62 100644 --- a/manager/Monitor.cpp +++ b/manager/Monitor.cpp @@ -110,6 +110,7 @@ void GetCPUInfo() { MACH3LOG_INFO("{}", TerminalToString("lscpu | grep -m 1 -E 'L3 |L3:'")); MACH3LOG_INFO("{}", TerminalToString("lscpu | grep -m 1 -E 'Thread.* per core:|Wątków na rdzeń:'")); MACH3LOG_INFO("{}", TerminalToString("lscpu | grep -m 1 -E '^CPU(:|\\(s\\)):?\\s+[0-9]+'")); + MACH3LOG_INFO("With available threads {}", M3::GetNThreads()); //KS: /proc/cpuinfo and lscpu holds much more info I have limited it but one can expand it if needed } @@ -287,5 +288,19 @@ void MaCh3Usage(int argc, char **argv){ throw MaCh3Exception(__FILE__, __LINE__); } } +} //end namespace +namespace M3 { +// *************************************************************************** +int GetNThreads() { +// *************************************************************************** + #ifdef MULTITHREAD + int nThreads = omp_get_max_threads(); + #else + int nThreads = 1; + #endif + + return nThreads; +} } //end namespace + diff --git a/manager/Monitor.h b/manager/Monitor.h index c01b62e15..df93a9348 100644 --- a/manager/Monitor.h +++ b/manager/Monitor.h @@ -7,25 +7,20 @@ #include #include -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" +// MaCh3 includes +#include "manager/MaCh3Logger.h" +#include "manager/MaCh3Exception.h" +#include "manager/YamlHelper.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT include #include "TTree.h" #include "TBranch.h" #include "TMacro.h" #include "TChain.h" #include "TStopwatch.h" -#pragma GCC diagnostic pop +_MaCh3_Safe_Include_End_ //} -// MaCh3 includes -#include "manager/MaCh3Logger.h" -#include "samplePDF/Structs.h" -#include "manager/MaCh3Exception.h" -#include "manager/YamlHelper.h" /// @file Monitor.h /// @brief System and monitoring utilities for printing system information and status updates. @@ -76,3 +71,8 @@ namespace MaCh3Utils { /// @details This function prints a simple usage guide for MaCh3 executables, typically called when incorrect arguments are passed. void MaCh3Usage(int argc, char **argv); } + +namespace M3 { + /// @brief number of threads which we need for example for TRandom3 + int GetNThreads(); +} diff --git a/manager/YamlHelper.h b/manager/YamlHelper.h index 51c445378..dc7275474 100644 --- a/manager/YamlHelper.h +++ b/manager/YamlHelper.h @@ -6,24 +6,20 @@ #include #include -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" +// MaCh3 Includes +#include "manager/MaCh3Exception.h" +#include "manager/Core.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT Includes #include "TMacro.h" #include "TList.h" #include "TObjString.h" -#pragma GCC diagnostic pop +_MaCh3_Safe_Include_End_ //} // yaml Includes #include "yaml-cpp/yaml.h" -// MaCh3 Includes -#include "manager/MaCh3Exception.h" - /// @file YamlHelper.h /// @brief Utility functions for handling YAML nodes /// @author Kamil Skwarczynski @@ -312,6 +308,13 @@ Type GetFromManager(const YAML::Node& node, Type defval, const std::string File /// @param Line number where function is called inline YAML::Node LoadYamlConfig(const std::string& filename, const std::string& File, const int Line) { // ********************** + // KS: YAML can be dumb and not throw error if you pass toml for example... + if (!(filename.length() >= 5 && filename.compare(filename.length() - 5, 5, ".yaml") == 0) && + !(filename.length() >= 4 && filename.compare(filename.length() - 4, 4, ".yml") == 0)) { + MACH3LOG_ERROR("Invalid file extension: {}", filename); + throw MaCh3Exception(File, Line); + } + try { return YAML::LoadFile(filename); } catch (const std::exception& e) { diff --git a/manager/manager.cpp b/manager/manager.cpp index 3a02d786f..7610c15b3 100644 --- a/manager/manager.cpp +++ b/manager/manager.cpp @@ -1,14 +1,9 @@ -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" +#include "manager/manager.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT include #include "TFile.h" -#pragma GCC diagnostic pop - -#include "manager/manager.h" +_MaCh3_Safe_Include_End_ //} // ************************* manager::manager(std::string const &filename) diff --git a/mcmc/FitterBase.cpp b/mcmc/FitterBase.cpp index 1918db645..9036179de 100644 --- a/mcmc/FitterBase.cpp +++ b/mcmc/FitterBase.cpp @@ -1,9 +1,11 @@ #include "FitterBase.h" +_MaCh3_Safe_Include_Start_ //{ #include "TRandom.h" #include "TStopwatch.h" #include "TTree.h" #include "TGraphAsymmErrors.h" +_MaCh3_Safe_Include_End_ //} #pragma GCC diagnostic ignored "-Wuseless-cast" @@ -406,6 +408,8 @@ void FitterBase::ProcessMCMC() { void FitterBase::DragRace(const int NLaps) { // ************************* MACH3LOG_INFO("Let the Race Begin!"); + MACH3LOG_INFO("All tests will be performed with {} threads", M3::GetNThreads()); + // Reweight the MC for(unsigned int ivs = 0; ivs < samples.size(); ivs++ ) { diff --git a/mcmc/MCMCProcessor.cpp b/mcmc/MCMCProcessor.cpp index d93827199..fa34df61f 100644 --- a/mcmc/MCMCProcessor.cpp +++ b/mcmc/MCMCProcessor.cpp @@ -1,15 +1,9 @@ #include "MCMCProcessor.h" -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +_MaCh3_Safe_Include_Start_ //{ #include "TChain.h" #include "TF1.h" -#pragma GCC diagnostic pop +_MaCh3_Safe_Include_End_ //} //Only if GPU is enabled #ifdef CUDA @@ -577,7 +571,7 @@ void MCMCProcessor::DrawPostfit() { CompLeg->SetLineStyle(0); CompLeg->SetBorderSize(0); - const double BottomMargin = Posterior->GetBottomMargin(); + const std::vector Margins = GetMargins(Posterior); Posterior->SetBottomMargin(0.2); OutputFile->cd(); @@ -692,7 +686,7 @@ void MCMCProcessor::DrawPostfit() { delete paramPlot_HPD; //KS: Return Margin to default one - Posterior->SetBottomMargin(BottomMargin); + SetMargins(Posterior, Margins); } // ********************* @@ -1445,8 +1439,7 @@ void MCMCProcessor::DrawCovariance() { void MCMCProcessor::DrawCorrelations1D() { // ********************* //KS: Store it as we go back to them at the end - const double TopMargin = Posterior->GetTopMargin(); - const double BottomMargin = Posterior->GetBottomMargin(); + const std::vector Margins = GetMargins(Posterior); const int OptTitle = gStyle->GetOptTitle(); Posterior->SetTopMargin(0.1); @@ -1577,8 +1570,7 @@ void MCMCProcessor::DrawCorrelations1D() { delete CorrDir; OutputFile->cd(); - Posterior->SetTopMargin(TopMargin); - Posterior->SetBottomMargin(BottomMargin); + SetMargins(Posterior, Margins); gStyle->SetOptTitle(OptTitle); } @@ -1736,10 +1728,7 @@ void MCMCProcessor::MakeTrianglePlot(const std::vector& ParNames, } //KS: Store it as we go back to them at the end - const double TopMargin = Posterior->GetTopMargin(); - const double BottomMargin = Posterior->GetBottomMargin(); - const double LeftMargin = Posterior->GetLeftMargin(); - const double RighMargin = Posterior->GetRightMargin(); + const std::vector Margins = GetMargins(Posterior); Posterior->SetTopMargin(0.001); Posterior->SetBottomMargin(0.001); Posterior->SetLeftMargin(0.001); @@ -1970,25 +1959,20 @@ void MCMCProcessor::MakeTrianglePlot(const std::vector& ParNames, for(int i = 0; i < nParamPlot; i++) { delete hpost_copy[i]; - for (int j = 0; j < nCredibleIntervals; ++j) - { + for (int j = 0; j < nCredibleIntervals; ++j) { delete hpost_cl[i][j]; } } for(int i = 0; i < Npad - nParamPlot; i++) { delete hpost_2D_copy[i]; - for (int j = 0; j < nCredibleRegions; ++j) - { + for (int j = 0; j < nCredibleRegions; ++j) { delete hpost_2D_cl[i][j]; } } //KS: Restore margin - Posterior->SetTopMargin(TopMargin); - Posterior->SetLeftMargin(BottomMargin); - Posterior->SetLeftMargin(LeftMargin); - Posterior->SetRightMargin(RighMargin); + SetMargins(Posterior, Margins); } // ************************** @@ -2601,10 +2585,7 @@ void MCMCProcessor::GetPolarPlot(const std::vector& ParNames){ // ************************** if(hpost[0] == nullptr) MakePostfit(); - const double TopMargin = Posterior->GetTopMargin(); - const double BottomMargin = Posterior->GetBottomMargin(); - const double LeftMargin = Posterior->GetLeftMargin(); - const double RightMargin = Posterior->GetRightMargin(); + std::vector Margins = GetMargins(Posterior); Posterior->SetTopMargin(0.1); Posterior->SetBottomMargin(0.1); @@ -2665,10 +2646,7 @@ void MCMCProcessor::GetPolarPlot(const std::vector& ParNames){ OutputFile->cd(); - Posterior->SetTopMargin(TopMargin); - Posterior->SetBottomMargin(BottomMargin); - Posterior->SetLeftMargin(LeftMargin); - Posterior->SetRightMargin(RightMargin); + SetMargins(Posterior, Margins); } // ************************** @@ -2724,7 +2702,7 @@ void MCMCProcessor::GetBayesFactor(const std::vector& ParNames, MACH3LOG_INFO("{} for {}", Name, ParNames[k]); MACH3LOG_INFO("Following Jeffreys Scale = {}", JeffreysScale); MACH3LOG_INFO("Following Dunne-Kaboth Scale = {}", DunneKabothScale); - std::cout<& Names, double old_chi = -1; double old_prior = -1; - if(FlatPrior[j]) - { + if(FlatPrior[j]) { old_prior = 1.0; - } - else - { + } else { old_chi = (ParameterPos[j] - OldCentral[j])/OldError[j]; old_prior = std::exp(-0.5 * old_chi * old_chi); } @@ -2991,7 +2966,6 @@ void MCMCProcessor::ParameterEvolution(const std::vector& Names, } const int IntervalsSize = nSteps/NIntervals[k]; - // ROOT won't overwrite gifs so we need to delete the file if it's there already int ret = system(fmt::format("rm {}.gif",Names[k]).c_str()); if (ret != 0){ @@ -3047,7 +3021,6 @@ void MCMCProcessor::ParameterEvolution(const std::vector& Names, // Diagnose the MCMC void MCMCProcessor::DiagMCMC() { // ************************** - // Prepare branches etc for DiagMCMC PrepareDiagMCMC(); @@ -3349,12 +3322,12 @@ void MCMCProcessor::AutoCorrelation_FFT() { std::vector> LagL(nDraw); // Arrays needed to perform FFT using ROOT - double* ACFFT = new double[nEntries](); // Main autocorrelation array - double* ParVals = new double[nEntries](); // Param values for full chain - double* ParValsFFTR = new double[nEntries](); // FFT Real part - double* ParValsFFTI = new double[nEntries](); // FFT Imaginary part - double* ParValsFFTSquare = new double[nEntries](); // FFT Absolute square - double* ParValsComplex = new double[nEntries](); // Input Imaginary values (0) + std::vector ACFFT(nEntries, 0.0); // Main autocorrelation array + std::vector ParVals(nEntries, 0.0); // Param values for full chain + std::vector ParValsFFTR(nEntries, 0.0); // FFT Real part + std::vector ParValsFFTI(nEntries, 0.0); // FFT Imaginary part + std::vector ParValsFFTSquare(nEntries, 0.0); // FFT Absolute square + std::vector ParValsComplex(nEntries, 0.0); // Input Imaginary values (0) // Create forward and reverse FFT objects. I don't love using ROOT here, // but it works so I can't complain @@ -3371,9 +3344,9 @@ void MCMCProcessor::AutoCorrelation_FFT() { } // Transform - fftf->SetPointsComplex(ParVals, ParValsComplex); + fftf->SetPointsComplex(ParVals.data(), ParValsComplex.data()); fftf->Transform(); - fftf->GetPointsComplex(ParValsFFTR, ParValsFFTI); + fftf->GetPointsComplex(ParValsFFTR.data(), ParValsFFTI.data()); // Square the results to get the power spectrum for (int i = 0; i < nEntries; ++i) { @@ -3381,9 +3354,9 @@ void MCMCProcessor::AutoCorrelation_FFT() { } // Transforming back gives the autocovariance - fftb->SetPointsComplex(ParValsFFTSquare, ParValsComplex); + fftb->SetPointsComplex(ParValsFFTSquare.data(), ParValsComplex.data()); fftb->Transform(); - fftb->GetPointsComplex(ACFFT, ParValsComplex); + fftb->GetPointsComplex(ACFFT.data(), ParValsComplex.data()); // Divide by norm to get autocorrelation double normAC = ACFFT[0]; @@ -3418,14 +3391,6 @@ void MCMCProcessor::AutoCorrelation_FFT() { //KS: This is different diagnostic however it relies on calculated Lag, thus we call it before we delete LagKPlots CalculateESS(nLags, LagL); - // Clean up - delete[] ACFFT; - delete[] ParVals; - delete[] ParValsFFTR; - delete[] ParValsFFTI; - delete[] ParValsFFTSquare; - delete[] ParValsComplex; - AutoCorrDir->Close(); delete AutoCorrDir; @@ -4317,3 +4282,28 @@ void MCMCProcessor::PrintInfo() const { MACH3LOG_INFO("# Osc params: \033[1;32m {} starting at {} \033[0m ", nParam[kOSCPar], ParamTypeStartPos[kOSCPar]); MACH3LOG_INFO("************************************************"); } + + +// ************************** +std::vector MCMCProcessor::GetMargins(const std::unique_ptr& Canv) const { +// ************************** + return std::vector{Canv->GetTopMargin(), Canv->GetBottomMargin(), + Canv->GetLeftMargin(), Canv->GetRightMargin()}; +} + +// ************************** +void MCMCProcessor::SetMargins(std::unique_ptr& Canv, const std::vector& margins) { +// ************************** + if (!Canv) { + MACH3LOG_ERROR("Canv is nullptr"); + throw MaCh3Exception(__FILE__, __LINE__); + } + if (margins.size() != 4) { + MACH3LOG_ERROR("Margin vector must have exactly 4 elements"); + throw MaCh3Exception(__FILE__, __LINE__); + } + Canv->SetTopMargin(margins[0]); + Canv->SetBottomMargin(margins[1]); + Canv->SetLeftMargin(margins[2]); + Canv->SetRightMargin(margins[3]); +} diff --git a/mcmc/MCMCProcessor.h b/mcmc/MCMCProcessor.h index 3a243f78a..3a4d767f7 100644 --- a/mcmc/MCMCProcessor.h +++ b/mcmc/MCMCProcessor.h @@ -8,13 +8,11 @@ #include #include -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +// MaCh3 includes +#include "mcmc/StatisticalUtils.h" +#include "samplePDF/HistogramUtils.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT includes #include "TFile.h" #include "TBranch.h" @@ -38,13 +36,9 @@ #include "TMath.h" #include "TMatrixDSymEigen.h" #include "TVirtualFFT.h" -#pragma GCC diagnostic pop +_MaCh3_Safe_Include_End_ //} -// MaCh3 includes -#include "mcmc/StatisticalUtils.h" -#include "samplePDF/HistogramUtils.h" - //KS: Joy of forward declaration https://gieseanw.wordpress.com/2018/02/25/the-joys-of-forward-declarations-results-from-the-real-world/ class TChain; class TF1; @@ -349,6 +343,11 @@ class MCMCProcessor { /// @author Richard Calland inline void PowerSpectrumAnalysis(); + /// @brief Get TCanvas margins, to be able to reset them if particular function need different margins + std::vector GetMargins(const std::unique_ptr& Canv) const; + /// @brief Set TCanvas margins to specified values + void SetMargins(std::unique_ptr& Canv, const std::vector& margins); + /// Name of MCMC file std::string MCMCFile; /// Output file suffix useful when running over same file with different settings diff --git a/mcmc/SampleSummary.cpp b/mcmc/SampleSummary.cpp index d9b3cf096..19a14f6c4 100644 --- a/mcmc/SampleSummary.cpp +++ b/mcmc/SampleSummary.cpp @@ -764,7 +764,7 @@ void SampleSummary::Write() { TempHistogram->GetXaxis()->SetTitle("-2LLH(Draw Fluc, Draw)"); TempHistogram->GetYaxis()->SetTitle("-2LLH(Data, Draw)"); TempHistogram->SetNameTitle((SampleNames[i]+"_drawfluc_draw").c_str(), (SampleNames[i]+"_drawfluc_draw").c_str()); - MakeCutLLH2D(TempHistogram); + Get2DBayesianpValue(TempHistogram); TempHistogram->Write(); delete TempHistogram; @@ -774,7 +774,7 @@ void SampleSummary::Write() { TempHistogram2->GetXaxis()->SetTitle("-2LLH(Pred Fluc, Draw)"); TempHistogram2->GetYaxis()->SetTitle("-2LLH(Data, Draw)"); TempHistogram2->SetNameTitle((SampleNames[i]+"_predfluc_draw").c_str(), (SampleNames[i]+"_predfluc_draw").c_str()); - MakeCutLLH2D(TempHistogram2); + Get2DBayesianpValue(TempHistogram2); TempHistogram2->Write(); delete TempHistogram2; @@ -784,7 +784,7 @@ void SampleSummary::Write() { TempHistogram3->GetXaxis()->SetTitle("-2LLH(Pred Fluc, Draw)"); TempHistogram3->GetYaxis()->SetTitle("-2LLH(Data, Draw)"); TempHistogram3->SetNameTitle((SampleNames[i]+"_rate_predfluc_draw").c_str(), (SampleNames[i]+"_rate_predfluc_draw").c_str()); - MakeCutLLH2D(TempHistogram3); + Get2DBayesianpValue(TempHistogram3); TempHistogram3->Write(); delete TempHistogram3; @@ -794,7 +794,7 @@ void SampleSummary::Write() { TempHistogram4->GetXaxis()->SetTitle(("-2LLH_{Draw Fluc, Draw} for " + SamplePDF->GetKinVarLabel(i, 0)).c_str()); TempHistogram4->GetYaxis()->SetTitle(("-2LLH_{Data, Draw} for " + SamplePDF->GetKinVarLabel(i, 0)).c_str()); TempHistogram4->SetNameTitle((SampleNames[i]+"_drawfluc_draw_ProjectX").c_str(), (SampleNames[i]+"_drawfluc_draw_ProjectX").c_str()); - MakeCutLLH2D(TempHistogram4); + Get2DBayesianpValue(TempHistogram4); TempHistogram4->Write(); delete TempHistogram4; @@ -1438,10 +1438,10 @@ void SampleSummary::MakeCutLLH() { MakeCutLLH1D(lnLHist_drawdata.get()); MakeCutLLH1D(lnLHist_drawflucdraw.get()); - MakeCutLLH2D(lnLDrawHist.get()); - MakeCutLLH2D(lnLFlucHist.get()); - MakeCutLLH2D(lnLDrawHistRate.get()); - MakeCutLLH2D(lnLFlucHist_ProjectX.get()); + Get2DBayesianpValue(lnLDrawHist.get()); + Get2DBayesianpValue(lnLFlucHist.get()); + Get2DBayesianpValue(lnLDrawHistRate.get()); + Get2DBayesianpValue(lnLFlucHist_ProjectX.get()); } // **************** @@ -1506,60 +1506,6 @@ void SampleSummary::MakeCutLLH1D(TH1D *Histogram, double llh_ref) { delete TempCanvas; } -// **************** -// Make the 2D cut distribution and give the 2D p-value -void SampleSummary::MakeCutLLH2D(TH2D *Histogram) { -// **************** - const double TotalIntegral = Histogram->Integral(); - // Count how many fills are above y=x axis - // This is the 2D p-value - double Above = 0.0; - for (int i = 0; i < Histogram->GetXaxis()->GetNbins(); ++i) { - const double xvalue = Histogram->GetXaxis()->GetBinCenter(i+1); - for (int j = 0; j < Histogram->GetYaxis()->GetNbins(); ++j) { - const double yvalue = Histogram->GetYaxis()->GetBinCenter(j+1); - // We're only interested in being _ABOVE_ the y=x axis - if (xvalue >= yvalue) { - Above += Histogram->GetBinContent(i+1, j+1); - } - } - } - const double pvalue = Above/TotalIntegral; - std::stringstream ss; - ss << int(Above) << "/" << int(TotalIntegral) << "=" << pvalue; - Histogram->SetTitle((std::string(Histogram->GetTitle())+"_"+ss.str()).c_str()); - - // Now add the TLine going diagonally - double minimum = Histogram->GetXaxis()->GetBinLowEdge(1); - if (Histogram->GetYaxis()->GetBinLowEdge(1) > minimum) { - minimum = Histogram->GetYaxis()->GetBinLowEdge(1); - } - double maximum = Histogram->GetXaxis()->GetBinLowEdge(Histogram->GetXaxis()->GetNbins()); - if (Histogram->GetYaxis()->GetBinLowEdge(Histogram->GetYaxis()->GetNbins()) < maximum) { - maximum = Histogram->GetYaxis()->GetBinLowEdge(Histogram->GetYaxis()->GetNbins()); - //KS: Extend by bin width to perfectly fit canvas - maximum += Histogram->GetYaxis()->GetBinWidth(Histogram->GetYaxis()->GetNbins()); - } - else maximum += Histogram->GetXaxis()->GetBinWidth(Histogram->GetXaxis()->GetNbins()); - auto TempLine = std::make_unique(minimum, minimum, maximum, maximum); - TempLine->SetLineColor(kRed); - TempLine->SetLineWidth(2); - - std::string Title = Histogram->GetName(); - Title += "_canv"; - TCanvas *TempCanvas = new TCanvas(Title.c_str(), Title.c_str(), 1024, 1024); - TempCanvas->SetGridx(); - TempCanvas->SetGridy(); - TempCanvas->cd(); - gStyle->SetPalette(81); - Histogram->SetMinimum(-0.01); - Histogram->Draw("colz"); - TempLine->Draw("same"); - - TempCanvas->Write(); - delete TempCanvas; -} - // **************** // Make the 1D Event Rate Hist void SampleSummary::MakeCutEventRate(TH1D *Histogram, const double DataRate) { diff --git a/mcmc/SampleSummary.h b/mcmc/SampleSummary.h index b357e20f8..be1c88253 100644 --- a/mcmc/SampleSummary.h +++ b/mcmc/SampleSummary.h @@ -88,8 +88,6 @@ class SampleSummary { inline void MakeCutLLH(); // Make the 1D cut distribution and give the 1D p-value inline void MakeCutLLH1D(TH1D *Histogram, double llh_ref = -999); - /// @brief Make the 2D cut distribution and give the 2D p-value - inline void MakeCutLLH2D(TH2D *Histogram); /// @brief Make the 1D Event Rate Hist inline void MakeCutEventRate(TH1D *Histogram, const double DataRate); /// @brief Make the fluctuated histograms (2D and 1D) for the chi2s diff --git a/mcmc/StatisticalUtils.cpp b/mcmc/StatisticalUtils.cpp index c6a48166c..9458db3f6 100644 --- a/mcmc/StatisticalUtils.cpp +++ b/mcmc/StatisticalUtils.cpp @@ -595,3 +595,57 @@ double GetModeError(TH1D* hpost) { return Mode_Error; } + +// **************** +// Make the 2D cut distribution and give the 2D p-value +void Get2DBayesianpValue(TH2D *Histogram) { +// **************** + const double TotalIntegral = Histogram->Integral(); + // Count how many fills are above y=x axis + // This is the 2D p-value + double Above = 0.0; + for (int i = 0; i < Histogram->GetXaxis()->GetNbins(); ++i) { + const double xvalue = Histogram->GetXaxis()->GetBinCenter(i+1); + for (int j = 0; j < Histogram->GetYaxis()->GetNbins(); ++j) { + const double yvalue = Histogram->GetYaxis()->GetBinCenter(j+1); + // We're only interested in being _ABOVE_ the y=x axis + if (xvalue >= yvalue) { + Above += Histogram->GetBinContent(i+1, j+1); + } + } + } + const double pvalue = Above/TotalIntegral; + std::stringstream ss; + ss << int(Above) << "/" << int(TotalIntegral) << "=" << pvalue; + Histogram->SetTitle((std::string(Histogram->GetTitle())+"_"+ss.str()).c_str()); + + // Now add the TLine going diagonally + double minimum = Histogram->GetXaxis()->GetBinLowEdge(1); + if (Histogram->GetYaxis()->GetBinLowEdge(1) > minimum) { + minimum = Histogram->GetYaxis()->GetBinLowEdge(1); + } + double maximum = Histogram->GetXaxis()->GetBinLowEdge(Histogram->GetXaxis()->GetNbins()); + if (Histogram->GetYaxis()->GetBinLowEdge(Histogram->GetYaxis()->GetNbins()) < maximum) { + maximum = Histogram->GetYaxis()->GetBinLowEdge(Histogram->GetYaxis()->GetNbins()); + //KS: Extend by bin width to perfectly fit canvas + maximum += Histogram->GetYaxis()->GetBinWidth(Histogram->GetYaxis()->GetNbins()); + } + else maximum += Histogram->GetXaxis()->GetBinWidth(Histogram->GetXaxis()->GetNbins()); + auto TempLine = std::make_unique(minimum, minimum, maximum, maximum); + TempLine->SetLineColor(kRed); + TempLine->SetLineWidth(2); + + std::string Title = Histogram->GetName(); + Title += "_canv"; + TCanvas *TempCanvas = new TCanvas(Title.c_str(), Title.c_str(), 1024, 1024); + TempCanvas->SetGridx(); + TempCanvas->SetGridy(); + TempCanvas->cd(); + gStyle->SetPalette(81); + Histogram->SetMinimum(-0.01); + Histogram->Draw("colz"); + TempLine->Draw("same"); + + TempCanvas->Write(); + delete TempCanvas; +} diff --git a/mcmc/StatisticalUtils.h b/mcmc/StatisticalUtils.h index e610c93db..5adc786dc 100644 --- a/mcmc/StatisticalUtils.h +++ b/mcmc/StatisticalUtils.h @@ -13,6 +13,20 @@ #include "manager/manager.h" #include "samplePDF/HistogramUtils.h" +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wuseless-cast" +#pragma GCC diagnostic ignored "-Wfloat-conversion" +#pragma GCC diagnostic ignored "-Wfloat-conversion" +#pragma GCC diagnostic ignored "-Wold-style-cast" +#pragma GCC diagnostic ignored "-Wconversion" +#pragma GCC diagnostic ignored "-Wformat-nonliteral" +// ROOT includes +#include "TCanvas.h" +#include "TLine.h" +#include "TROOT.h" +#include "TStyle.h" +#pragma GCC diagnostic pop + /// @file StatisticalUtils.h /// @brief Utility functions for statistical interpretations in MaCh3 /// @author Kamil Skwarczynski @@ -156,3 +170,8 @@ double GetPValueFromZScore(const double zScore); /// @brief Get the mode error from a TH1D /// @param hpost hist from which we extract mode error double GetModeError(TH1D* hpost); + +/// @brief Calculates the 2D Bayesian p-value and generates a visualization. +/// @param Histogram A pointer to a TH2D histogram object. The function modifies the histogram's title to include the p-value information. +/// @warning The canvas is saved to the current ROOT file using `TempCanvas->Write()`. +void Get2DBayesianpValue(TH2D *Histogram); diff --git a/plotting/GetPostfitParamPlots.cpp b/plotting/GetPostfitParamPlots.cpp index 882f464a1..89dfd1add 100644 --- a/plotting/GetPostfitParamPlots.cpp +++ b/plotting/GetPostfitParamPlots.cpp @@ -3,13 +3,11 @@ #include #include -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +#include "plottingUtils/plottingUtils.h" +#include "plottingUtils/plottingManager.h" +#include "samplePDF/Structs.h" + +_MaCh3_Safe_Include_Start_ //{ #include "TROOT.h" #include "TGaxis.h" #include "TString.h" @@ -26,10 +24,7 @@ #include "TCandle.h" #include "TFrame.h" #include "TGraphAsymmErrors.h" -#pragma GCC diagnostic pop - -#include "plottingUtils/plottingUtils.h" -#include "plottingUtils/plottingManager.h" +_MaCh3_Safe_Include_End_ //} /// @file GetPostfitParamPlots /// This script generates post-fit parameter plots. The central postfit value is diff --git a/plotting/MatrixPlotter.cpp b/plotting/MatrixPlotter.cpp index 62a760097..264b02b97 100755 --- a/plotting/MatrixPlotter.cpp +++ b/plotting/MatrixPlotter.cpp @@ -53,7 +53,7 @@ TH2D* GetSubMatrix(TH2D *MatrixFull, const std::string& Title, const std::vector void SetupInfo(const std::string& Config, std::vector& Title, std::vector>& Params) { // Load the YAML file - YAML::Node config = YAML::LoadFile(Config); + YAML::Node config = M3OpenConfig(Config); // Access the "MatrixPlotter" section YAML::Node settings = config["MatrixPlotter"]; diff --git a/plotting/plottingUtils/inputManager.cpp b/plotting/plottingUtils/inputManager.cpp index 01ce016a0..39865e678 100644 --- a/plotting/plottingUtils/inputManager.cpp +++ b/plotting/plottingUtils/inputManager.cpp @@ -5,7 +5,7 @@ namespace MaCh3Plotting { // this is the constructor with user specified translation config file InputManager::InputManager(const std::string &translationConfigName) { // read the config file - _translatorConfig = YAML::LoadFile(translationConfigName); + _translatorConfig = M3OpenConfig(translationConfigName); MACH3LOG_DEBUG("InputManager: have loaded translation config file"); diff --git a/plotting/plottingUtils/plottingManager.cpp b/plotting/plottingUtils/plottingManager.cpp index e99d6e4cd..47bde186b 100644 --- a/plotting/plottingUtils/plottingManager.cpp +++ b/plotting/plottingUtils/plottingManager.cpp @@ -27,7 +27,7 @@ void PlottingManager::initialise() { /// config files and make sure that all provided options are valid and all necessary options are provided /// as it can be pretty annoying and difficult to identify what's going wrong when yaml just fails /// to find an option at runtime - _plottingConfig = YAML::LoadFile(_configFileName); + _plottingConfig = M3OpenConfig(_configFileName); MACH3LOG_DEBUG("Initialising PlottingManager with plotting congif {}", _configFileName); // read options from the config diff --git a/samplePDF/FarDetectorCoreInfoStruct.h b/samplePDF/FarDetectorCoreInfoStruct.h index 4919ba824..b50315831 100644 --- a/samplePDF/FarDetectorCoreInfoStruct.h +++ b/samplePDF/FarDetectorCoreInfoStruct.h @@ -11,7 +11,7 @@ struct FarDetectorCoreInfo { FarDetectorCoreInfo& operator=(FarDetectorCoreInfo const &other) = delete; FarDetectorCoreInfo& operator=(FarDetectorCoreInfo &&other) = delete; - ~FarDetectorCoreInfo(){ delete [] isNC; } + ~FarDetectorCoreInfo(){if(isNC != nullptr) delete [] isNC;} int nEvents; ///< how many MC events are there double ChannelIndex; @@ -21,8 +21,6 @@ struct FarDetectorCoreInfo { std::vector nupdg; std::vector nupdgUnosc; - int SampleDetID; - //THe x_var and y_vars that you're binning in std::vector x_var; std::vector y_var; diff --git a/samplePDF/HistogramUtils.cpp b/samplePDF/HistogramUtils.cpp index 75b147ed7..5e6ff193e 100644 --- a/samplePDF/HistogramUtils.cpp +++ b/samplePDF/HistogramUtils.cpp @@ -1,15 +1,9 @@ -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +#include "samplePDF/HistogramUtils.h" + +_MaCh3_Safe_Include_Start_ //{ #include "TList.h" #include "TObjArray.h" -#pragma GCC diagnostic pop - -#include "samplePDF/HistogramUtils.h" +_MaCh3_Safe_Include_End_ //} // ************************************************** //KS: ROOT changes something with binning when moving from ROOT 5 to ROOT 6. If you open ROOT5 produced file with ROOT6 you will be missing 9 last bins diff --git a/samplePDF/HistogramUtils.h b/samplePDF/HistogramUtils.h index 6e3b34076..872c16b7f 100644 --- a/samplePDF/HistogramUtils.h +++ b/samplePDF/HistogramUtils.h @@ -1,20 +1,14 @@ #pragma once -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +// MaCh3 includes +#include "samplePDF/Structs.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT include #include "TGraphAsymmErrors.h" #include "TObjString.h" #include "TRandom3.h" -#pragma GCC diagnostic pop - -// MaCh3 inlcudes -#include "samplePDF/Structs.h" +_MaCh3_Safe_Include_End_ //} /// @file HistogramUtils.h /// @author Will Parker diff --git a/samplePDF/Structs.h b/samplePDF/Structs.h index dbb7151d9..084e5823f 100644 --- a/samplePDF/Structs.h +++ b/samplePDF/Structs.h @@ -1,76 +1,16 @@ #pragma once -/// Run low or high memory versions of structs -/// N.B. for 64 bit systems sizeof(float) == sizeof(double) so not a huge effect -/// KS: Need more testing on FD -namespace M3 { -#ifdef _LOW_MEMORY_STRUCTS_ -/// Custom floating point (float or double) -using float_t = float; -/// Custom integer (int or short int) -using int_t = short; -/// Custom unsigned integer (unsigned short int or unsigned int) -using uint_t = unsigned short; -#else -using float_t = double; -using int_t = int; -using uint_t = unsigned; -#endif -} - -/// KS: noexcept can help with performance but is terrible for debugging, this is meant to help easy way of of turning it on or off. In near future move this to struct or other central class. -//#define SafeException -#ifndef SafeException -#define _noexcept_ noexcept -#else -#define _noexcept_ -#endif - -/// KS: Using restrict limits the effects of pointer aliasing, aiding optimizations. While reading I found that there might be some compilers which don't have __restrict__. As always we use _restrict_ to more easily turn off restrict in the code -#ifndef DEBUG -#define _restrict_ __restrict__ -#else -#define _restrict_ -#endif - -/// Number of overflow bins in TH2Poly, -#define _TH2PolyOverflowBins_ 9 - -constexpr static const double _BAD_DOUBLE_ = -999.99; -constexpr static const int _BAD_INT_ = -999; - -constexpr static const double _DEFAULT_RETURN_VAL_ = -999999.123456; - -// Some commonly used variables to which we set pointers to -constexpr static const double Unity = 1.; -constexpr static const double Zero = 0.; -constexpr static const float Unity_F = 1.; -constexpr static const float Zero_F = 0.; -constexpr static const int Unity_Int = 1; - // C++ includes -#include -#include -#include -#include -#include #include #include #include -#ifdef MULTITHREAD -#include "omp.h" -#endif - +// MaCh3 includes #include "manager/MaCh3Exception.h" #include "manager/MaCh3Logger.h" +#include "manager/Core.h" -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" +_MaCh3_Safe_Include_Start_ //{ // ROOT include #include "TSpline.h" #include "TObjString.h" @@ -80,8 +20,7 @@ constexpr static const int Unity_Int = 1; #include "TH1.h" // NuOscillator includes #include "Constants/OscillatorConstants.h" -#pragma GCC diagnostic pop - +_MaCh3_Safe_Include_End_ //} /// @file Structs.h /// @author Asher Kaboth @@ -99,6 +38,45 @@ std::vector MakeVector( const T (&data)[N] ) { return std::vector(data, data+N); } +// ******************* +/// @brief Generic cleanup function +template +void CleanVector(std::vector& vec) { +// ******************* + vec.clear(); + vec.shrink_to_fit(); +} + +// ******************* +/// @brief Generic cleanup function +template +void CleanVector(std::vector>& vec) { +// ******************* + for (auto& innerVec : vec) { + innerVec.clear(); + innerVec.shrink_to_fit(); + } + vec.clear(); + vec.shrink_to_fit(); +} + +// ******************* +/// @brief Generic cleanup function +template +void CleanVector(std::vector>>& vec) { +// ******************* + for (auto& inner2DVec : vec) { + for (auto& innerVec : inner2DVec) { + innerVec.clear(); + innerVec.shrink_to_fit(); + } + inner2DVec.clear(); + inner2DVec.shrink_to_fit(); + } + vec.clear(); + vec.shrink_to_fit(); +} + // ******************* /// @brief KS: This is mad way of converting string to int. Why? To be able to use string with switch constexpr unsigned int str2int(const char* str, int h = 0) { @@ -567,19 +545,4 @@ namespace MaCh3Utils { return NuOscillatorFlavour; } // *************************** - - /// @brief DB Anything added here must be of the form 2^X, where X is an integer - /// @warning DB Used to contain which DetIDs are supported - static const std::unordered_mapKnownDetIDsMap({ - {0,1}, //ND - {1,8}, //FD - {2,16}, //SK1Rmu - {3,32}, //Nova - {4,64}, //Atm SubGeV e-like - {5,128}, //Atm SubGeV mu-like - {6,256}, //Atm MultiGeV e-like - {7,512}, //Atm MultiGeV mu-like - }); - static const int nKnownDetIDs = int(KnownDetIDsMap.size()); - } // end MaCh3Utils namespace diff --git a/samplePDF/samplePDFBase.h b/samplePDF/samplePDFBase.h index bc1accd0b..d706bb511 100644 --- a/samplePDF/samplePDFBase.h +++ b/samplePDF/samplePDFBase.h @@ -3,13 +3,12 @@ //C++ includes #include -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +//MaCh3 includes +#include "samplePDF/Structs.h" +#include "samplePDF/HistogramUtils.h" +#include "manager/manager.h" + +_MaCh3_Safe_Include_Start_ //{ //ROOT includes #include "TTree.h" #include "TH1D.h" @@ -19,12 +18,7 @@ #include "TROOT.h" #include "TRandom3.h" #include "TString.h" -#pragma GCC diagnostic pop - -//MaCh3 includes -#include "samplePDF/Structs.h" -#include "samplePDF/HistogramUtils.h" -#include "manager/manager.h" +_MaCh3_Safe_Include_End_ //} /// @brief Class responsible for handling implementation of samples used in analysis, reweighting and returning LLH class samplePDFBase diff --git a/samplePDF/samplePDFFDBase.cpp b/samplePDF/samplePDFFDBase.cpp index 2559a9cd0..515ed82a6 100644 --- a/samplePDF/samplePDFFDBase.cpp +++ b/samplePDF/samplePDFFDBase.cpp @@ -29,7 +29,7 @@ samplePDFFDBase::samplePDFFDBase(std::string ConfigFileName, covarianceXsec* xse samplePDFFD_array = nullptr; samplePDFFD_data = nullptr; - + SampleDetID = ""; SampleManager = std::unique_ptr(new manager(ConfigFileName.c_str())); } @@ -72,7 +72,7 @@ void samplePDFFDBase::ReadSampleConfig() MACH3LOG_ERROR("ID not defined in {}, please add this!", SampleManager->GetFileName()); throw MaCh3Exception(__FILE__, __LINE__); } - SampleDetID = SampleManager->raw()["DetID"].as(); + SampleDetID = SampleManager->raw()["DetID"].as(); if (!CheckNodeExists(SampleManager->raw(), "NuOsc", "NuOscConfigFile")) { MACH3LOG_ERROR("NuOsc::NuOscConfigFile is not defined in {}, please add this!", SampleManager->GetFileName()); @@ -228,7 +228,7 @@ void samplePDFFDBase::fill2DHist() // ************************************************ /// @function samplePDFFDBase::SetupSampleBinning() -/// @brief Function to setup the binning of your sample histograms and the underlying +/// @brief Function to setup the binning of your sample histograms and the underlying /// arrays that get handled in fillArray() and fillArray_MP(). /// The SampleXBins are filled in the daughter class from the sample config file. /// This "passing" can be removed. @@ -297,44 +297,6 @@ bool samplePDFFDBase::IsEventSelected(const int iSample, const int iEvent) { return false; } } - - //DB To avoid unnecessary checks, now return false rather than setting bool to true and continuing to check - return true; -} - -// ************************************************ -bool samplePDFFDBase::IsEventSelected(const std::vector< std::string >& ParameterStr, const int iSample, const int iEvent) { -// ************************************************ - double Val; - for (unsigned int iSelection=0;iSelection=SelectionBounds[iSelection][1])) { - return false; - } - } - - //DB To avoid unnecessary checks, now return false rather than setting bool to true and continuing to check - return true; -} - -// ************************************************ -//Same as the function above but just acts on the vector and the event -bool samplePDFFDBase::IsEventSelected(const std::vector< std::string >& ParameterStr, - const std::vector< std::vector > &SelectionCuts, - const int iSample, const int iEvent) { -// ************************************************ - - double Val; - for (unsigned int iSelection=0;iSelection= SelectionCuts[iSelection][1] && SelectionCuts[iSelection][0] != -999){ - return false; - } - else if(Val < SelectionCuts[iSelection][0] && SelectionCuts[iSelection][1] != -999){ - return false; - } - } //DB To avoid unnecessary checks, now return false rather than setting bool to true and continuing to check return true; } @@ -349,14 +311,14 @@ void samplePDFFDBase::reweight() { //You only need to do these things if OscCov has been initialised //if not then you're not considering oscillations if (OscCov) { - std::vector OscVec(OscCov->GetNumParams()); - for (int iPar=0;iParGetNumParams();iPar++) { + std::vector OscVec(OscParams.size()); + for (size_t iPar = 0; iPar < OscParams.size(); ++iPar) { #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wuseless-cast" - OscVec[iPar] = M3::float_t(OscCov->getParProp(iPar)); + OscVec[iPar] = M3::float_t(*OscParams[iPar]); #pragma GCC diagnostic pop } - for (int iSample=0;iSampleCalculateProbabilities(OscVec); } } @@ -383,7 +345,7 @@ void samplePDFFDBase::fillArray() { #else //ETA we should probably store this in samplePDFFDBase size_t nXBins = int(XBinEdges.size()-1); - size_t nYBins = int(YBinEdges.size()-1); + //size_t nYBins = int(YBinEdges.size()-1); PrepFunctionalParameters(); if(SplineHandler){ @@ -398,14 +360,12 @@ void samplePDFFDBase::fillArray() { continue; } - std::cout << "Event passed selection, here we go!!" << std::endl; - double splineweight = 1.0; double normweight = 1.0; double totalweight = 1.0; if(SplineHandler){ - splineweight *= CalcXsecWeightSpline(iSample, iEvent); + splineweight = CalcWeightSpline(iSample, iEvent); } //DB Catch negative spline weights and skip any event with a negative event. Previously we would set weight to zero and continue but that is inefficient. Do this on a spline-by-spline basis if (splineweight <= 0.){ @@ -414,7 +374,7 @@ void samplePDFFDBase::fillArray() { } //Loop over stored normalisation and function pointers - normweight *= CalcXsecWeightNorm(iSample, iEvent); + normweight = CalcWeightNorm(iSample, iEvent); //DB Catch negative norm weights and skip any event with a negative event. Previously we would set weight to zere and continue but that is inefficient if (normweight <= 0.){ @@ -551,7 +511,7 @@ void samplePDFFDBase::fillArray_MP() { //As weights were skdet::fParProp, and we use the non-shifted erec, we might as well cache the corresponding fParProp index for each event and the pointer to it if(SplineHandler){ - splineweight *= CalcXsecWeightSpline(iSample, iEvent); + splineweight = CalcWeightSpline(iSample, iEvent); } //DB Catch negative spline weights and skip any event with a negative event. Previously we would set weight to zero and continue but that is inefficient if (splineweight <= 0.){ @@ -559,7 +519,7 @@ void samplePDFFDBase::fillArray_MP() { continue; } - normweight *= CalcXsecWeightNorm(iSample, iEvent); + normweight = CalcWeightNorm(iSample, iEvent); //DB Catch negative norm weights and skip any event with a negative event. Previously we would set weight to zero and continue but that is inefficient if (normweight <= 0.){ MCSamples[iSample].xsec_w[iEvent] = 0.; @@ -593,7 +553,7 @@ void samplePDFFDBase::fillArray_MP() { int YBinToFill = MCSamples[iSample].NomYBin[iEvent]; //DB Check to see if momentum shift has moved bins - //DB - First, check to see if the event is still in the nominal bin + //DB - First, check to see if the event is still in the nominal bin if (XVar < MCSamples[iSample].rw_upper_xbinedge[iEvent] && XVar >= MCSamples[iSample].rw_lower_xbinedge[iEvent]) { XBinToFill = MCSamples[iSample].NomXBin[iEvent]; } @@ -651,7 +611,6 @@ void samplePDFFDBase::fillArray_MP() { } #endif - // ************************************************** // Helper function to reset the data and MC histograms void samplePDFFDBase::ResetHistograms() { @@ -670,7 +629,7 @@ void samplePDFFDBase::ResetHistograms() { // *************************************************************************** // Calculate the spline weight for one event -M3::float_t samplePDFFDBase::CalcXsecWeightSpline(const int iSample, const int iEvent) const { +M3::float_t samplePDFFDBase::CalcWeightSpline(const int iSample, const int iEvent) const { // *************************************************************************** M3::float_t xsecw = 1.0; //DB Xsec syst @@ -683,7 +642,7 @@ M3::float_t samplePDFFDBase::CalcXsecWeightSpline(const int iSample, const int i // *************************************************************************** // Calculate the normalisation weight for one event -M3::float_t samplePDFFDBase::CalcXsecWeightNorm(const int iSample, const int iEvent) const { +M3::float_t samplePDFFDBase::CalcWeightNorm(const int iSample, const int iEvent) const { // *************************************************************************** M3::float_t xsecw = 1.0; //Loop over stored normalisation and function pointers @@ -730,6 +689,13 @@ void samplePDFFDBase::SetupNormParameters() { } } } + // If not debugging let's clear memory + #ifndef DEBUG + for (size_t iSample = 0; iSample < MCSamples.size(); ++iSample) { + MCSamples[iSample].xsec_norms_bins.clear(); + MCSamples[iSample].xsec_norms_bins.shrink_to_fit(); + } + #endif } // ************************************************ @@ -806,17 +772,15 @@ void samplePDFFDBase::CalcXsecNormsBins(int iSample) { // All normalisations are just 1 bin for 2015, so bin = index (where index is just the bin for that normalisation) int bin = (*it).index; - //If syst on applies to a particular detector - if ((XsecCov->GetParDetID(bin) & SampleDetID)==SampleDetID) { - XsecBins.push_back(bin); - MACH3LOG_TRACE("Event {}, will be affected by dial {}", iEvent, (*it).name); - #ifdef DEBUG - VerboseCounter[std::distance(xsec_norms.begin(), it)]++; - #endif - } + XsecBins.push_back(bin); + MACH3LOG_TRACE("Event {}, will be affected by dial {}", iEvent, (*it).name); + #ifdef DEBUG + VerboseCounter[std::distance(xsec_norms.begin(), it)]++; + #endif + //} } // end iteration over xsec_norms } // end if (xsecCov) - fdobj->xsec_norms_bins[iEvent]=XsecBins; + fdobj->xsec_norms_bins[iEvent] = XsecBins; }//end loop over events #ifdef DEBUG MACH3LOG_DEBUG("Channel {}", iSample); @@ -1337,6 +1301,8 @@ void samplePDFFDBase::SetupNuOscillator() { } // end loop over events }// end loop over channels delete OscillFactory; + + OscParams = OscCov->GetOscParsFromDetID(SampleDetID); } M3::float_t samplePDFFDBase::GetEventWeight(const int iSample, const int iEntry) const { @@ -1377,9 +1343,9 @@ void samplePDFFDBase::fillSplineBins() { throw MaCh3Exception(__FILE__, __LINE__); } MCSamples[i].xsec_spline_pointers[j].resize(MCSamples[i].nxsec_spline_pointers[j]); - for(int spline=0; splineretPointer(EventSplines[spline][0], EventSplines[spline][1], EventSplines[spline][2], + MCSamples[i].xsec_spline_pointers[j][spline] = SplineHandler->retPointer(EventSplines[spline][0], EventSplines[spline][1], EventSplines[spline][2], EventSplines[spline][3], EventSplines[spline][4], EventSplines[spline][5], EventSplines[spline][6]); } } @@ -1401,7 +1367,7 @@ double samplePDFFDBase::GetLikelihood() { double negLogL = 0.; #ifdef MULTITHREAD - #pragma omp parallel for reduction(+:negLogL) + #pragma omp parallel for collapse(2) reduction(+:negLogL) #endif for (int xBin = 0; xBin < nXBins; ++xBin) { @@ -1461,7 +1427,6 @@ void samplePDFFDBase::InitialiseSingleFDMCObject(int iSample, int nEvents_) { #else fdobj->osc_w_pointer.resize(nEvents, &Unity); #endif - fdobj->SampleDetID = -1; for(int iEvent = 0 ; iEvent < fdobj->nEvents ; ++iEvent){ fdobj->isNC[iEvent] = false; @@ -1484,7 +1449,6 @@ void samplePDFFDBase::InitialiseSplineObject() { } SplineHandler->AddSample(samplename, SampleDetID, spline_filepaths, SplineVarNames); - SplineHandler->PrintArrayDimension(); SplineHandler->CountNumberOfLoadedSplines(false, 1); SplineHandler->TransferToMonolith(); diff --git a/samplePDF/samplePDFFDBase.h b/samplePDF/samplePDFFDBase.h index 8d07a25a4..2aace0efe 100644 --- a/samplePDF/samplePDFFDBase.h +++ b/samplePDF/samplePDFFDBase.h @@ -21,6 +21,7 @@ class samplePDFFDBase : public samplePDFBase //######################################### Functions ######################################### /// @param ConfigFileName Name of config to initialise the sample object samplePDFFDBase(std::string ConfigFileName, covarianceXsec* xsec_cov, covarianceOsc* osc_cov = nullptr); + /// @brief destructor virtual ~samplePDFFDBase(); int GetNDim(){return nDimensions;} //DB Function to differentiate 1D or 2D binning @@ -82,7 +83,7 @@ class samplePDFFDBase : public samplePDFBase virtual void SetupSplines() = 0; //DB Require all objects to have a function which reads in the MC - // @brief Initialise any variables that your experiment specific samplePDF needs + /// @brief Initialise any variables that your experiment specific samplePDF needs virtual void Init() = 0; /// @brief Experiment specific setup, returns the number of events which were loaded @@ -126,15 +127,13 @@ class samplePDFFDBase : public samplePDFBase virtual void applyShifts(int iSample, int iEvent){(void) iSample; (void) iEvent;}; /// @brief DB Function which determines if an event is selected, where Selection double looks like {{ND280KinematicTypes Var1, douuble LowBound} bool IsEventSelected(const int iSample, const int iEvent); - bool IsEventSelected(const std::vector& ParameterStr, const int iSample, const int iEvent); - bool IsEventSelected(const std::vector& ParameterStr, const std::vector> &SelectionCuts, const int iSample, const int iEvent); /// @brief Check whether a normalisation systematic affects an event or not void CalcXsecNormsBins(int iSample); /// @brief Calculate the spline weight for a given event - M3::float_t CalcXsecWeightSpline(const int iSample, const int iEvent) const; + M3::float_t CalcWeightSpline(const int iSample, const int iEvent) const; /// @brief Calculate the norm weight for a given event - M3::float_t CalcXsecWeightNorm(const int iSample, const int iEvent) const; + M3::float_t CalcWeightNorm(const int iSample, const int iEvent) const; /// @brief Calculate weights for function parameters /// @@ -206,7 +205,7 @@ class samplePDFFDBase : public samplePDFBase /// @brief Keep track of the dimensions of the sample binning int nDimensions = _BAD_INT_; /// @brief A unique ID for each sample based on powers of two for quick binary operator comparisons - int SampleDetID = _BAD_INT_; + std::string SampleDetID; /// holds "TrueNeutrinoEnergy" and the strings used for the sample binning. std::vector SplineBinnedVars; @@ -215,7 +214,8 @@ class samplePDFFDBase : public samplePDFBase /// @brief Information to store for normalisation pars std::vector xsec_norms; - + /// pointer to osc params, since not all params affect every sample, we perform some operations before hand for speed + std::vector OscParams; //=========================================================================== //DB Vectors to store which kinematic cuts we apply //like in XsecNorms but for events in sample. Read in from sample yaml file diff --git a/splines/SplineBase.h b/splines/SplineBase.h index 7dc3e9bc1..62494626f 100644 --- a/splines/SplineBase.h +++ b/splines/SplineBase.h @@ -4,13 +4,12 @@ #include #include -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" +// MaCh3 includes +#include "manager/MaCh3Logger.h" +#include "splines/SplineStructs.h" +#include "splines/SplineCommon.h" + +_MaCh3_Safe_Include_Start_ //{ // ROOT include #include "TFile.h" #include "TH1F.h" @@ -21,16 +20,7 @@ #include "TCanvas.h" #include "TStyle.h" #include "TTree.h" -#pragma GCC diagnostic pop - -#ifdef MULTITHREAD -#include "omp.h" -#endif - -// MaCh3 includes -#include "manager/MaCh3Logger.h" -#include "splines/SplineStructs.h" -#include "splines/SplineCommon.h" +_MaCh3_Safe_Include_End_ //} /// @brief Base class for calculating weight from spline class SplineBase { diff --git a/splines/SplineStructs.h b/splines/SplineStructs.h index 6a7a239b6..ff56a0fb4 100644 --- a/splines/SplineStructs.h +++ b/splines/SplineStructs.h @@ -11,7 +11,6 @@ #pragma GCC diagnostic ignored "-Wfloat-conversion" #pragma GCC diagnostic ignored "-Wnull-dereference" - /// @file SplineStructs.h /// @brief Contains structures and helper functions for handling spline representations of systematic parameters in the MaCh3. /// @author Clarence Wret diff --git a/splines/splineFDBase.cpp b/splines/splineFDBase.cpp index 32c9a8de7..16023eb6a 100644 --- a/splines/splineFDBase.cpp +++ b/splines/splineFDBase.cpp @@ -1,6 +1,5 @@ #include "splineFDBase.h" #include -#include "samplePDF/Structs.h" #pragma GCC diagnostic ignored "-Wuseless-cast" #pragma GCC diagnostic ignored "-Wfloat-conversion" @@ -29,21 +28,32 @@ splineFDBase::~splineFDBase(){ //**************************************** void splineFDBase::cleanUpMemory() { //**************************************** - //Call once everything's been allocated in samplePDFSKBase, cleans up junk from memory! //Not a huge saving but it's better than leaving everything up to the compiler MACH3LOG_INFO("Cleaning up spline memory"); + CleanVector(indexvec); + CleanVector(SplineFileParPrefixNames); + CleanVector(GlobalSystIndex); + CleanVector(SplineModeVecs); + CleanVector(UniqueSystNames); + CleanVector(SplineInterpolationTypes); + - indexvec.clear(); - indexvec.shrink_to_fit(); - SplineFileParPrefixNames.clear(); - SplineFileParPrefixNames.shrink_to_fit(); + + for (auto& Binning2D : SplineBinning) { + for (auto& Binning1D : Binning2D) { + for (TAxis* axis : Binning1D) { + delete axis; + } + Binning1D.clear(); + Binning1D.shrink_to_fit(); + } + Binning2D.clear(); + Binning2D.shrink_to_fit(); + } SplineBinning.clear(); SplineBinning.shrink_to_fit(); - GlobalSystIndex.clear(); - GlobalSystIndex.shrink_to_fit(); - UniqueSystNames.clear(); - UniqueSystNames.shrink_to_fit(); + //Really make sure all the memory is cleared for(auto Spline : splinevec_Monolith){ if(Spline){delete Spline;} @@ -54,7 +64,10 @@ void splineFDBase::cleanUpMemory() { } //**************************************** -bool splineFDBase::AddSample(std::string SampleName, int DetID, std::vector OscChanFileNames, std::vector SplineVarNames) +bool splineFDBase::AddSample(const std::string& SampleName, + const std::string& DetID, + const std::vector& OscChanFileNames, + const std::vector& SplineVarNames) //Adds samples to the large array //**************************************** { @@ -75,8 +88,6 @@ bool splineFDBase::AddSample(std::string SampleName, int DetID, std::vector SplineParsIndex_Sample_temp = xsec->GetSplineParsIndexFromDetID(DetID); - std::vector SplineFileParPrefixNames_Sample = xsec->GetSplineParsNamesFromDetID(DetID); SplineFileParPrefixNames.push_back(SplineFileParPrefixNames_Sample); @@ -189,7 +200,7 @@ void splineFDBase::TransferToMonolith() xcoeff_arr[iCoeff+i]=tmpXCoeffArr[i]; for(int j=0; j<4; j++){ - manycoeff_arr[(iCoeff+i)*4+j]=tmpManyCoeffArr[i*4+j]; + manycoeff_arr[(iCoeff+i)*4+j]=tmpManyCoeffArr[i*4+j]; } } delete[] tmpXCoeffArr; @@ -225,23 +236,19 @@ void splineFDBase::Evaluate() { // ETA - find the spline segment that the current parameter // value is in. This is now extremely similar to the // function in SplineMonolith.cpp -void splineFDBase::FindSplineSegment() +void splineFDBase::FindSplineSegment() { //**************************************** -{ //HW okay let's try this, we delete+refill a new array which we'll fill with x-s for our segment - #ifdef MULTITHREAD - #pragma omp parallel //for schedule(dynamic) - #endif - for (int iSyst = 0; iSyst < nUniqueSysts; iSyst++) { - int nPoints = UniqueSystNKnots[iSyst]; - std::vector xArray = UniqueSystXPts[iSyst]; + for (int iSyst = 0; iSyst < nUniqueSysts; ++iSyst) { + const int nPoints = UniqueSystNKnots[iSyst]; + // KS: By reference to not copy memory like a fool + const std::vector& xArray = UniqueSystXPts[iSyst]; // Get the variation for this reconfigure for the ith parameter - int GlobalIndex = UniqueSystIndices[iSyst]; - - M3::float_t xvar = M3::float_t(xsec->getParProp(GlobalIndex)); + const int GlobalIndex = UniqueSystIndices[iSyst]; + const M3::float_t xvar = M3::float_t(xsec->getParProp(GlobalIndex)); - xVarArray[iSyst]=xvar; + xVarArray[iSyst] = xvar; M3::int_t segment = 0; M3::int_t kHigh = M3::int_t(nPoints - 1); @@ -279,23 +286,6 @@ void splineFDBase::FindSplineSegment() if (segment >= nPoints-1 && nPoints > 1){segment = M3::int_t(nPoints-2);} UniqueSystCurrSegment[iSyst] = segment; - - //#ifdef DEBUG - // if (SplineInfoArray[i].xPts[segment] > xvar && segment != 0) { - // std::cerr << "Found a segment which is _ABOVE_ the variation!" << std::endl; - // std::cerr << "IT SHOULD ALWAYS BE BELOW! (except when segment 0)" << std::endl; - // std::cerr << "Spline: "<< i << std::endl; - // - // std::cerr << "Found segment = " << segment << std::endl; - // std::cerr << "Doing variation = " << xvar << std::endl; - // std::cerr << "x in spline = " << SplineInfoArray[i].xPts[segment] << std::endl; - // for (_M3::int_t_ j = 0; j < SplineInfoArray[j].nPts; ++j) { - // std::cerr << " " << j << " = " << SplineInfoArray[i].xPts[j] << std::endl; - // } - // std::cerr << __FILE__ << ":" << __LINE__ << std::endl; - // throw; - // } - //#endif } //end loop over params } @@ -306,51 +296,40 @@ void splineFDBase::CalcSplineWeights() #ifdef MULTITHREAD #pragma omp parallel for simd #endif - for (unsigned int iCoeff = 0; iCoeff < uniquecoeffindices.size(); iCoeff++) + for (size_t iCoeff = 0; iCoeff < uniquecoeffindices.size(); ++iCoeff) { + const int iSpline = uniquecoeffindices[iCoeff]; + const short int uniqueIndex = short(uniquesplinevec_Monolith[iSpline]); + const short int currentsegment = short(UniqueSystCurrSegment[uniqueIndex]); - int iSpline = uniquecoeffindices[iCoeff]; - short int uniqueIndex=short(uniquesplinevec_Monolith[iSpline]); - short int currentsegment=short(UniqueSystCurrSegment[uniqueIndex]); - - int segCoeff = coeffindexvec[iSpline]+currentsegment; - + const int segCoeff = coeffindexvec[iSpline]+currentsegment; + const int coeffOffset = segCoeff * 4; // These are what we can extract from the TSpline3 - M3::float_t x = xcoeff_arr[segCoeff]; - M3::float_t y = manycoeff_arr[(segCoeff)*4+kCoeffY]; - M3::float_t b = manycoeff_arr[(segCoeff)*4+kCoeffB]; - M3::float_t c = manycoeff_arr[(segCoeff)*4+kCoeffC]; - M3::float_t d = manycoeff_arr[(segCoeff)*4+kCoeffD]; + const M3::float_t x = xcoeff_arr[segCoeff]; + const M3::float_t y = manycoeff_arr[coeffOffset+kCoeffY]; + const M3::float_t b = manycoeff_arr[coeffOffset+kCoeffB]; + const M3::float_t c = manycoeff_arr[coeffOffset+kCoeffC]; + const M3::float_t d = manycoeff_arr[coeffOffset+kCoeffD]; // Get the variation for this reconfigure for the ith parameter - M3::float_t xvar = xVarArray[uniqueIndex]; + const M3::float_t xvar = xVarArray[uniqueIndex]; // The Delta(x) - M3::float_t dx = xvar - x; + const M3::float_t dx = xvar - x; //Speedy 1% time boost https://en.cppreference.com/w/c/numeric/math/fma (see ND code!) - M3::float_t weight = 0; -#ifdef _LOW_MEMORY_STRUCTS_ - weight = std::fmaf(dx, std::fmaf(dx, std::fmaf(dx, d, c), b), y); -#else - weight = std::fma(dx, std::fma(dx, std::fma(dx, d, c), b), y); -#endif + M3::float_t weight = M3::fmaf_t(dx, M3::fmaf_t(dx, M3::fmaf_t(dx, d, c), b), y); //This is the speedy version of writing dx^3+b*dx^2+c*dx+d - //ETA - do we need this? We check later for negative weights and I wonder if this is even //possible with the fmaf line above? - if(weight<0){weight=0;} //Stops is getting negative weights + if(weight < 0){weight = 0.;} //Stops is getting negative weights -// LP - ignore the diagnostic here as it is only useless if M3::float_t = double -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" - weightvec_Monolith[iSpline] = double(weight); -#pragma GCC diagnostic pop + weightvec_Monolith[iSpline] = weight; } } //**************************************** -void splineFDBase::BuildSampleIndexingArray(std::string SampleName) +void splineFDBase::BuildSampleIndexingArray(const std::string& SampleName) //Creates an array to be filled with monolith indexes for each sample (allows for indexing between 7D binning and 1D Vector) //Only need 1 indexing array everything else interfaces with this to get binning properties //**************************************** @@ -525,7 +504,7 @@ int splineFDBase::CountNumberOfLoadedSplines(bool NonFlat, int Verbosity) { // Loop over oscillation channels for (unsigned int iSyst = 0; iSyst < indexvec[iSample][iOscChan].size(); iSyst++) { // Loop over systematics - for (unsigned int iMode = 0; iMode < indexvec[iSample][iOscChan][iSyst].size(); iMode++) + for (unsigned int iMode = 0; iMode < indexvec[iSample][iOscChan][iSyst].size(); iMode++) { // Loop over modes for (unsigned int iVar1 = 0; iVar1 < indexvec[iSample][iOscChan][iSyst][iMode].size(); iVar1++) { // Loop over first dimension @@ -585,8 +564,8 @@ void splineFDBase::PrepForReweight() { std::string SystName = SplineFileParPrefixNames[iSample][iSyst]; bool FoundSyst = false; - //ETA - this always seems to be empty to begin with?? - //so this loop never gets used? + //ETA - this always seems to be empty to begin with?? + //so this loop never gets used? for (unsigned int iFoundSyst = 0; iFoundSyst < UniqueSystNames.size(); iFoundSyst++) { if (SystName == UniqueSystNames[iFoundSyst]) @@ -641,10 +620,10 @@ void splineFDBase::PrepForReweight() { break; } }//osc loop end - //ETA - only push back unique name if a non-flat response has been found - if(FoundNonFlatSpline){ - UniqueSystNames.push_back(SystName); - } + //ETA - only push back unique name if a non-flat response has been found + if(FoundNonFlatSpline){ + UniqueSystNames.push_back(SystName); + } if (!FoundNonFlatSpline) { @@ -675,7 +654,7 @@ void splineFDBase::PrepForReweight() { UniqueSystSplines[iSpline]->GetKnot(iKnot, xPoint, yPoint); UniqueSystXPts[iSpline][iKnot] = xPoint; } - //ETA - let this just be set as the first segment by default + //ETA - let this just be set as the first segment by default UniqueSystCurrSegment[iSpline] = 0; xVarArray[iSpline]=0; } @@ -742,10 +721,10 @@ void splineFDBase::getSplineCoeff_SepMany(int splineindex, M3::float_t* &xArray, } } - for(int i=0; iGetCoeff(i, x, y, b, c, d); - // Store the coefficients for each knot contiguously in memory - // 4 because manyArray stores y,b,c,d + // Store the coefficients for each knot contiguously in memory + // 4 because manyArray stores y,b,c,d xArray[i] = x; manyArray[i*4] = y; manyArray[i*4+1] = b; @@ -782,7 +761,7 @@ std::string splineFDBase::getDimLabel(int iSample, unsigned int Axis) } //Returns sample index in -int splineFDBase::getSampleIndex(std::string SampleName){ +int splineFDBase::getSampleIndex(const std::string& SampleName){ int SampleIndex = -1; for (unsigned int iSample = 0; iSample < SampleNames.size(); iSample++) { @@ -800,7 +779,7 @@ int splineFDBase::getSampleIndex(std::string SampleName){ } //**************************************** -void splineFDBase::PrintSampleDetails(std::string SampleName) +void splineFDBase::PrintSampleDetails(const std::string& SampleName) //**************************************** { int iSample = getSampleIndex(SampleName); @@ -813,59 +792,33 @@ void splineFDBase::PrintSampleDetails(std::string SampleName) } //**************************************** -void splineFDBase::PrintArrayDetails(std::string SampleName) +void splineFDBase::PrintArrayDetails(const std::string& SampleName) //**************************************** { int iSample = getSampleIndex(SampleName); int nOscChannels = int(indexvec[iSample].size()); - MACH3LOG_INFO("Sample {} has {} oscillation channels", iSample, nOscChannels); + MACH3LOG_INFO("Sample {} has {} oscillation channels", SampleName, nOscChannels); for (int iOscChan = 0; iOscChan < nOscChannels; iOscChan++) { int nSysts = int(indexvec[iSample][iOscChan].size()); MACH3LOG_INFO("Oscillation channel {} has {} systematics", iOscChan, nSysts); } -} -//**************************************** -void splineFDBase::PrintArrayDimension() { -//**************************************** MACH3LOG_INFO("#----------------------------------------------------------------------------------------------------------------------------------#"); - MACH3LOG_INFO("Array dimensions.."); - MACH3LOG_INFO("{:<20}{}", "nSamples:", indexvec.size()); - - MACH3LOG_INFO("{:<20}", "nOscChans:"); - std::string oscChans; - for (unsigned int iSample = 0; iSample < indexvec.size(); iSample++) { - oscChans += fmt::format("{} ", indexvec[iSample].size()); - } - MACH3LOG_INFO("{}", oscChans); - - MACH3LOG_INFO("{:<20}", "nSysts:"); - for (unsigned int iSample = 0; iSample < indexvec.size(); iSample++) { - std::string systs = fmt::format("\tSample: {}\t", iSample); - for (unsigned int iOscChan = 0; iOscChan < indexvec[iSample].size(); iOscChan++) { - systs += fmt::format("{} ", indexvec[iSample][iOscChan].size()); + MACH3LOG_INFO("Printing no. of modes affected by each systematic for each oscillation channel"); + for (unsigned int iOscChan = 0; iOscChan < indexvec[iSample].size(); iOscChan++) { + std::string modes = fmt::format("OscChan: {}\t", iOscChan); + for (unsigned int iSyst = 0; iSyst < indexvec[iSample][iOscChan].size(); iSyst++) { + modes += fmt::format("{} ", indexvec[iSample][iOscChan][iSyst].size()); } - MACH3LOG_INFO("{}", systs); - } - - MACH3LOG_INFO("{:<20}", "nModes:"); - for (unsigned int iSample = 0; iSample < indexvec.size(); iSample++) { - MACH3LOG_INFO("\tSample: {}\t--------------------------", iSample); - for (unsigned int iOscChan = 0; iOscChan < indexvec[iSample].size(); iOscChan++) { - std::string modes = fmt::format("\t\tOscChan: {}\t", iOscChan); - for (unsigned int iSyst = 0; iSyst < indexvec[iSample][iOscChan].size(); iSyst++) { - modes += fmt::format("{} ", indexvec[iSample][iOscChan][iSyst].size()); - } - MACH3LOG_INFO("{}", modes); - } - MACH3LOG_INFO(""); // Empty line for spacing + MACH3LOG_INFO("{}", modes); } MACH3LOG_INFO("#----------------------------------------------------------------------------------------------------------------------------------#"); } + //**************************************** -bool splineFDBase::isValidSplineIndex(std::string SampleName, int iOscChan, int iSyst, int iMode, int iVar1, int iVar2, int iVar3) +bool splineFDBase::isValidSplineIndex(const std::string& SampleName, int iOscChan, int iSyst, int iMode, int iVar1, int iVar2, int iVar3) //**************************************** { int iSample=getSampleIndex(SampleName); diff --git a/splines/splineFDBase.h b/splines/splineFDBase.h index faee085bb..0fddff4ba 100644 --- a/splines/splineFDBase.h +++ b/splines/splineFDBase.h @@ -1,20 +1,14 @@ #pragma once -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wuseless-cast" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wfloat-conversion" -#pragma GCC diagnostic ignored "-Wold-style-cast" -#pragma GCC diagnostic ignored "-Wconversion" -#pragma GCC diagnostic ignored "-Wformat-nonliteral" -// ROOT includes -#include "TH3F.h" -#pragma GCC diagnostic pop - -//MaCh3 +//MaCh3 includes #include "samplePDF/Structs.h" #include "splines/SplineBase.h" +_MaCh3_Safe_Include_Start_ //{ +// ROOT includes +#include "TH3F.h" +_MaCh3_Safe_Include_End_ //} + /// @brief Bin-by-bin class calculating response for spline parameters. /// @see For more details, visit the [Wiki](https://github.com/mach3-software/MaCh3/wiki/05.-Splines). /// @author Dan Barrow @@ -24,23 +18,23 @@ class splineFDBase : public SplineBase { /// @todo ETA - do all of these functions and members actually need to be public? public: /// @brief Constructor - splineFDBase(covarianceXsec *xsec_ = NULL); + splineFDBase(covarianceXsec *xsec_ = nullptr); /// @brief Destructor /// @todo it need some love virtual ~splineFDBase(); - void SetupSplines(); - void SetupSplines(int BinningOpt); /// @brief CW: This Eval should be used when using two separate x,{y,a,b,c,d} arrays /// to store the weights; probably the best one here! Same thing but pass parameter /// spline segments instead of variations void Evaluate(); - //Spline Monolith things - //Essential methods used externally - //Move these to splineFDBase in core - bool AddSample(std::string SampleName, int DetID, std::vector OscChanFileNames, std::vector SplineVarNames); - void TransferToMonolith(); + /// @brief add oscillation channel to spline monolith + bool AddSample(const std::string& SampleName, + const std::string& DetID, + const std::vector& OscChanFileNames, + const std::vector& SplineVarNames); + void TransferToMonolith(); + /// @brief Remove setup variables not needed for spline evaluations void cleanUpMemory(); //Have to define this in your own class @@ -57,18 +51,18 @@ class splineFDBase : public SplineBase { int CountNumberOfLoadedSplines(bool NonFlat=false, int Verbosity=0); int getNDim(int BinningOpt); std::string getDimLabel(int BinningOpt, unsigned int Axis); - int getSampleIndex(std::string SampleName); - bool isValidSplineIndex(std::string SampleName, int iSyst, int iOscChan, int iMode, int iVar1, int iVar2, int iVar3); + int getSampleIndex(const std::string& SampleName); + bool isValidSplineIndex(const std::string& SampleName, int iSyst, int iOscChan, int iMode, int iVar1, int iVar2, int iVar3); - void BuildSampleIndexingArray(std::string SampleName); + void BuildSampleIndexingArray(const std::string& SampleName); void PrepForReweight(); void getSplineCoeff_SepMany(int splineindex, M3::float_t *& xArray, M3::float_t *&manyArray); void PrintBinning(TAxis* Axis); - void PrintSampleDetails(std::string SampleName); - void PrintArrayDetails(std::string SampleName); - void PrintArrayDimension(); - - const M3::float_t* retPointer(int sample, int oscchan, int syst, int mode, int var1bin, int var2bin, int var3bin){ + void PrintSampleDetails(const std::string& SampleName); + void PrintArrayDetails(const std::string& SampleName); + + const M3::float_t* retPointer(const int sample, const int oscchan, const int syst, const int mode, + const int var1bin, const int var2bin, const int var3bin) const{ int index = indexvec[sample][oscchan][syst][mode][var1bin][var2bin][var3bin]; return &weightvec_Monolith[index]; } @@ -85,19 +79,17 @@ class splineFDBase : public SplineBase { //And now the actual member variables std::vector SampleNames; - std::vector BinningOpts; std::vector Dimensions; std::vector> DimensionLabels; - std::vector DetIDs; + std::vector DetIDs; std::vector nSplineParams; std::vector nOscChans; - std::vector< std::vector > SplineParsIndex; std::vector< std::vector< std::vector > > SplineBinning; std::vector< std::vector > SplineFileParPrefixNames; - //A vector of vectors of the spline modes that a systematic applies to - //This gets compared against the event mode to figure out if a syst should - //apply to an event or not + /// A vector of vectors of the spline modes that a systematic applies to + /// This gets compared against the event mode to figure out if a syst should + /// apply to an event or not std::vector< std::vector< std::vector > > SplineModeVecs; /// @brief This holds the global spline index and is used to grab the current parameter value /// to evaluate splines at. Each internal vector will be of size of the number of spline @@ -117,7 +109,7 @@ class splineFDBase : public SplineBase { /// @brief Variables related to determined which modes have splines and which piggy-back of other modes std::vector< std::vector< std::vector< std::vector< std::vector< std::vector< std::vector< int > > > > > > > indexvec; std::vector coeffindexvec; - std::vectoruniquecoeffindices; //Unique coefficient indices + std::vector uniquecoeffindices; //Unique coefficient indices std::vector< TSpline3_red* > splinevec_Monolith; @@ -125,11 +117,14 @@ class splineFDBase : public SplineBase { int MonolithIndex; int CoeffIndex; - //Probably need to clear these arrays up at some point + /// Probably need to clear these arrays up at some point M3::float_t *xVarArray; - bool *isflatarray; // Need to keep track of which splines are flat and which aren't - M3::float_t *xcoeff_arr; //x coefficients for each spline - M3::float_t *manycoeff_arr; //ybcd coefficients for each spline + /// Need to keep track of which splines are flat and which aren't + bool *isflatarray; + /// x coefficients for each spline + M3::float_t *xcoeff_arr; + /// ybcd coefficients for each spline + M3::float_t *manycoeff_arr; std::vector weightvec_Monolith; std::vector uniquesplinevec_Monolith;