diff --git a/.github/workflows/CDImage.yml b/.github/workflows/CDImage.yml index f336d1f4..1e49ad61 100644 --- a/.github/workflows/CDImage.yml +++ b/.github/workflows/CDImage.yml @@ -20,15 +20,12 @@ jobs: - os: Ubuntu22.04 file: Doc/MaCh3DockerFiles/Ubuntu22.04/Dockerfile tag_latest: ubuntu22.04latest - tag_release: ubuntu22.04v1.1.0 - os: Alma9 file: Doc/MaCh3DockerFiles/Alma9/Dockerfile tag_latest: alma9latest - tag_release: alma9v1.1.0 - os: Fedora32 file: Doc/MaCh3DockerFiles/Fedora32/Dockerfile tag_latest: fedora32latest - tag_release: fedora32v1.1.0 name: Image CD ${{ matrix.os }} @@ -43,6 +40,7 @@ jobs: run: | if [ "${{ github.event_name }}" == 'release' ]; then docker build . --file ${{ matrix.file }} --tag ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag_release }} --build-arg MACH3_VERSION=develop + docker build . --file ${{ matrix.file }} --tag ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.os }}v${{ github.event.release.tag_name }} --build-arg MACH3_VERSION=develop else docker build . --file ${{ matrix.file }} --tag ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag_latest }} --build-arg MACH3_VERSION=develop fi @@ -50,7 +48,7 @@ jobs: - name: Push Docker image run: | if [ "${{ github.event_name }}" == 'release' ]; then - docker push ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag_release }} + docker push ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.os }}v${{ github.event.release.tag_name }} else docker push ghcr.io/${{ github.repository_owner }}/mach3:${{ matrix.tag_latest }} fi diff --git a/Diagnostics/DiagMCMC.cpp b/Diagnostics/DiagMCMC.cpp index cbf256d4..f1245b73 100644 --- a/Diagnostics/DiagMCMC.cpp +++ b/Diagnostics/DiagMCMC.cpp @@ -2,26 +2,10 @@ #include "manager/manager.h" - -void DiagMCMC(std::string inputFile, std::string config); - -int main(int argc, char *argv[]) { - SetMaCh3LoggerFormat(); - if (argc != 3) - { - MACH3LOG_ERROR("How to use: DiagMCMC MCMC_Output.root config"); - throw MaCh3Exception(__FILE__ , __LINE__ ); - } - MACH3LOG_INFO("Producing single fit output"); - std::string filename = argv[1]; - std::string config = argv[2]; - DiagMCMC(filename, config); - - return 0; -} - - -void DiagMCMC(std::string inputFile, std::string config) +/// @brief Main function creating MCMCProcessor and calling MCMC Diagnostic +/// @param inputFile MCMC Chain +/// @param config Config file with settings +void DiagMCMC(const std::string& inputFile, const std::string& config) { MACH3LOG_INFO("File for study: {}", inputFile); @@ -47,4 +31,17 @@ void DiagMCMC(std::string inputFile, std::string config) delete Processor; } +int main(int argc, char *argv[]) { + SetMaCh3LoggerFormat(); + if (argc != 3) + { + MACH3LOG_ERROR("How to use: DiagMCMC MCMC_Output.root config"); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + MACH3LOG_INFO("Producing single fit output"); + std::string filename = argv[1]; + std::string config = argv[2]; + DiagMCMC(filename, config); + return 0; +} diff --git a/Diagnostics/ProcessMCMC.cpp b/Diagnostics/ProcessMCMC.cpp index 88d27b4e..80ef1e9c 100644 --- a/Diagnostics/ProcessMCMC.cpp +++ b/Diagnostics/ProcessMCMC.cpp @@ -2,15 +2,19 @@ #include "mcmc/MCMCProcessor.h" #include "manager/manager.h" -inline void ProcessMCMC(std::string inputFile); +/// @brief Main function processing MCMC and Producing plots +inline void ProcessMCMC(const std::string& inputFile); +/// @brief Function producing comparison of posterior and more betwen a few MCMC chains inline void MultipleProcessMCMC(); inline void CalcBayesFactor(MCMCProcessor* Processor); inline void CalcSavageDickey(MCMCProcessor* Processor); inline void CalcBipolarPlot(MCMCProcessor* Processor); inline void GetTrianglePlot(MCMCProcessor* Processor); -inline void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, std::string inputFile); +inline void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, const std::string& inputFile); inline void ReweightPrior(MCMCProcessor* Processor); -inline TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, std::string title); +/// @brief KS: Convert TMatrix to TH2D, mostly useful for making fancy plots +inline TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, const std::string& title); +/// @brief KS: Perform KS test to check if two posteriors for the same parameter came from the same distribution inline void KolmogorovSmirnovTest(MCMCProcessor** Processor, TCanvas* Posterior, TString canvasname); int nFiles; @@ -59,7 +63,7 @@ int main(int argc, char *argv[]) return 0; } -void ProcessMCMC(std::string inputFile) +void ProcessMCMC(const std::string& inputFile) { MACH3LOG_INFO("File for study: {} with config {}", inputFile, config); // Make the processor) @@ -391,8 +395,8 @@ void GetTrianglePlot(MCMCProcessor* Processor) { } } -//KS: You validate stability of posterior covariance matrix, you set burn calc cov matrix increase burn calc again and compare. By performing such operation several hundred times we can check when matrix becomes stable -void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, std::string inputFile) +/// @brief KS: You validate stability of posterior covariance matrix, you set burn calc cov matrix increase burn calc again and compare. By performing such operation several hundred times we can check when matrix becomes stable +void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, const std::string& inputFile) { //Turn of plots from Processor Processor->SetPrintToPDF(false); @@ -584,7 +588,7 @@ void ReweightPrior(MCMCProcessor* Processor) } //KS: Convert TMatrix to TH2D, mostly useful for making fancy plots -TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, std::string title) +TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, const std::string& title) { TH2D* hMatrix = new TH2D(title.c_str(), title.c_str(), Matrix->GetNrows(), 0.0, Matrix->GetNrows(), Matrix->GetNcols(), 0.0, Matrix->GetNcols()); for(int i = 0; i < Matrix->GetNrows(); i++) diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp index f8f10c6a..ae383cc0 100644 --- a/Diagnostics/RHat.cpp +++ b/Diagnostics/RHat.cpp @@ -30,8 +30,8 @@ #include "manager/manager.h" -//KS: This exe is meant to calculate R hat estimator. For a well converged this distribution should be centred at one. -//Based on Gelman et. al. arXiv:1903.08008v5 +/// KS: This exe is meant to calculate R hat estimator. For a well converged this distribution should be centred at one. +/// Based on Gelman et. al. arXiv:1903.08008v5 // ******************* int Ntoys; diff --git a/Doc/Doxyfile b/Doc/Doxyfile index 5ad4430e..83b4c161 100644 --- a/Doc/Doxyfile +++ b/Doc/Doxyfile @@ -115,7 +115,7 @@ ABBREVIATE_BRIEF = # description. # The default value is: NO. -ALWAYS_DETAILED_SEC = NO +ALWAYS_DETAILED_SEC = YES # If the INLINE_INHERITED_MEMB tag is set to YES, doxygen will show all # inherited members of a class in the documentation of that class as if those @@ -189,6 +189,14 @@ QT_AUTOBRIEF = NO MULTILINE_CPP_IS_BRIEF = NO +# By default Python docstrings are displayed as preformatted text and doxygen's +# special commands cannot be used. By setting PYTHON_DOCSTRING to NO the +# doxygen's special commands can be used and the contents of the docstring +# documentation blocks is shown as doxygen documentation. +# The default value is: YES. + +PYTHON_DOCSTRING = YES + # If the INHERIT_DOCS tag is set to YES then an undocumented member inherits the # documentation from any documented member that it re-implements. # The default value is: YES. @@ -1897,7 +1905,7 @@ ENABLE_PREPROCESSING = YES # The default value is: NO. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -MACRO_EXPANSION = NO +MACRO_EXPANSION = YES # If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES then # the macro expansion is limited to the macros specified with the PREDEFINED and @@ -2044,7 +2052,7 @@ HIDE_UNDOC_RELATIONS = YES # set to NO # The default value is: NO. -HAVE_DOT = NO +HAVE_DOT = YES # The DOT_NUM_THREADS specifies the number of dot invocations doxygen is allowed # to run in parallel. When set to 0 doxygen will base this on the number of diff --git a/cmake/Templates/Doxyfile.in b/cmake/Templates/Doxyfile.in index 59f5f822..10d1ba57 100644 --- a/cmake/Templates/Doxyfile.in +++ b/cmake/Templates/Doxyfile.in @@ -115,7 +115,7 @@ ABBREVIATE_BRIEF = # description. # The default value is: NO. -ALWAYS_DETAILED_SEC = NO +ALWAYS_DETAILED_SEC = YES # If the INLINE_INHERITED_MEMB tag is set to YES, doxygen will show all # inherited members of a class in the documentation of that class as if those @@ -189,6 +189,14 @@ QT_AUTOBRIEF = NO MULTILINE_CPP_IS_BRIEF = NO +# By default Python docstrings are displayed as preformatted text and doxygen's +# special commands cannot be used. By setting PYTHON_DOCSTRING to NO the +# doxygen's special commands can be used and the contents of the docstring +# documentation blocks is shown as doxygen documentation. +# The default value is: YES. + +PYTHON_DOCSTRING = YES + # If the INHERIT_DOCS tag is set to YES then an undocumented member inherits the # documentation from any documented member that it re-implements. # The default value is: YES. @@ -1897,7 +1905,7 @@ ENABLE_PREPROCESSING = YES # The default value is: NO. # This tag requires that the tag ENABLE_PREPROCESSING is set to YES. -MACRO_EXPANSION = NO +MACRO_EXPANSION = YES # If the EXPAND_ONLY_PREDEF and MACRO_EXPANSION tags are both set to YES then # the macro expansion is limited to the macros specified with the PREDEFINED and @@ -2044,7 +2052,7 @@ HIDE_UNDOC_RELATIONS = YES # set to NO # The default value is: NO. -HAVE_DOT = NO +HAVE_DOT = YES # The DOT_NUM_THREADS specifies the number of dot invocations doxygen is allowed # to run in parallel. When set to 0 doxygen will base this on the number of diff --git a/covariance/covarianceOsc.cpp b/covariance/covarianceOsc.cpp index e10b19a6..02142fe4 100644 --- a/covariance/covarianceOsc.cpp +++ b/covariance/covarianceOsc.cpp @@ -16,7 +16,6 @@ covarianceOsc::covarianceOsc(const char* name, const char *file) TVectorD* osc_baseline = (TVectorD*)infile->Get("osc_baseline"); TVectorD* osc_density = (TVectorD*)infile->Get("osc_density"); - double fScale = 1.0; //KS: Save all necessary information from covariance for(int io = 0; io < _fNumPar; io++) @@ -27,7 +26,7 @@ covarianceOsc::covarianceOsc(const char* name, const char *file) _fCurrVal[io] = _fPropVal[io] = _fPreFitValue[io]; _fError[io] = (*osc_sigma)(io); - _fIndivStepScale[io] = fScale * (*osc_stepscale)(io); + _fIndivStepScale[io] = (*osc_stepscale)(io); //KS: Set flat prior //HW: Might as well set it for everything in case default behaviour changes @@ -243,13 +242,13 @@ void covarianceOsc::CheckOrderOfParams() { std::vector wrongParam; bool wrongMatrix = false; - if(strcmp( _fNames[0].c_str(),"sin2th_12") != 0 ){wrongParam.push_back(0); wrongMatrix = true;}; - if(strcmp( _fNames[1].c_str(),"sin2th_23") != 0 ){wrongParam.push_back(1); wrongMatrix = true;}; - if(strcmp( _fNames[2].c_str(),"sin2th_13") != 0 ){wrongParam.push_back(2); wrongMatrix = true;}; - if(strcmp( _fNames[3].c_str(),"delm2_12") != 0 ){wrongParam.push_back(3); wrongMatrix = true;}; - if(strcmp( _fNames[4].c_str(),"delm2_23") != 0 ){wrongParam.push_back(4); wrongMatrix = true;}; - if(strcmp( _fNames[5].c_str(),"delta_cp") != 0 ){wrongParam.push_back(5); wrongMatrix = true;}; - if(PerformBetaStudy && strcmp( _fNames[6].c_str(),"beta") != 0 ){wrongParam.push_back(6); wrongMatrix = true;}; + if(_fNames[0] != "sin2th_12"){wrongParam.push_back(0); wrongMatrix = true;}; + if(_fNames[1] != "sin2th_23"){wrongParam.push_back(1); wrongMatrix = true;}; + if(_fNames[2] != "sin2th_13"){wrongParam.push_back(2); wrongMatrix = true;}; + if(_fNames[3] != "delm2_12"){wrongParam.push_back(3); wrongMatrix = true;}; + if(_fNames[4] != "delm2_23"){wrongParam.push_back(4); wrongMatrix = true;}; + if(_fNames[5] != "delta_cp"){wrongParam.push_back(5); wrongMatrix = true;}; + if(PerformBetaStudy && _fNames[6] != "beta"){wrongParam.push_back(6); wrongMatrix = true;}; if(wrongMatrix) { diff --git a/covariance/covarianceOsc.h b/covariance/covarianceOsc.h index 3a47130f..dc3bad71 100644 --- a/covariance/covarianceOsc.h +++ b/covariance/covarianceOsc.h @@ -67,4 +67,3 @@ class covarianceOsc : public covarianceBase /// Enum for special treatment of beta. int kBeta; }; - diff --git a/manager/MaCh3Exception.h b/manager/MaCh3Exception.h index 15b67912..637be3e3 100644 --- a/manager/MaCh3Exception.h +++ b/manager/MaCh3Exception.h @@ -21,8 +21,9 @@ class MaCh3Exception : public std::exception { std::string fileName = (lastSlashPos != std::string::npos) ? File.substr(lastSlashPos + 1) : File; errorMessage = ((Message.empty()) ? "Terminating MaCh3" : Message); - - MACH3LOG_ERROR("Find me here: {}::{}", fileName, Line); + // KS: Set logger format where we only have have "information type", line would be confusing + spdlog::set_pattern("[%^%l%$] %v"); + MACH3LOG_ERROR("Find me here: {}::{}", fileName, Line); } /// @brief Returns the error message associated with this exception. diff --git a/manager/Monitor.cpp b/manager/Monitor.cpp index b3da675d..1a3aeb6e 100644 --- a/manager/Monitor.cpp +++ b/manager/Monitor.cpp @@ -60,7 +60,7 @@ std::string GetMaCh3Version() { // Check if the file is opened successfully if (!versionFile.is_open()) { MACH3LOG_ERROR("Error: Couldn't open version.h {}", file); - throw; + throw MaCh3Exception(__FILE__, __LINE__); } std::string line; @@ -100,7 +100,7 @@ void GetOSInfo() { // ************************ //KS: Simple function retrieving CPU info -void GetCPUInfo(){ +void GetCPUInfo() { // ************************ //KS: Use -m 1 to limit to only one grep because in one computing node there is a lot of CPU which are the same @@ -162,7 +162,7 @@ std::string TerminalToString(const char* cmd) { std::string result; std::unique_ptr pipe(popen(cmd, "r"), pclose); if (!pipe) { - throw std::runtime_error("popen() failed!"); + throw MaCh3Exception(__FILE__, __LINE__, "popen() failed!"); } while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { result += buffer.data(); @@ -259,9 +259,8 @@ int getValue(std::string Type){ //Note: this value is in KB! } else { - std::cerr << "Not supported getValue: " << Type << std::endl; - std::cerr << __FILE__ << ":" << __LINE__ << std::endl; - throw; + MACH3LOG_ERROR("Not supported getValue: {}", Type); + throw MaCh3Exception(__FILE__, __LINE__); } return result; @@ -293,6 +292,4 @@ void PrintConfig(const YAML::Node& node){ } } - - } diff --git a/mcmc/FitterBase.cpp b/mcmc/FitterBase.cpp index 6672ab72..d792fa19 100644 --- a/mcmc/FitterBase.cpp +++ b/mcmc/FitterBase.cpp @@ -24,7 +24,7 @@ FitterBase::FitterBase(manager * const man) : fitMan(man) { stepClock = new TStopwatch; #ifdef DEBUG // Fit summary and debug info - debug = fitMan->raw()["General"]["Debug"].as(); + debug = GetFromManager(fitMan->raw()["General"]["Debug"], false); #endif std::string outfile = fitMan->raw()["General"]["OutputFile"].as(); @@ -73,7 +73,7 @@ FitterBase::FitterBase(manager * const man) : fitMan(man) { syst_llh = nullptr; TotalNSamples = 0; - fTestLikelihood = GetFromManager(fitMan->raw()["General"]["Fitter"]["FitTestLikelihood"], false);; + fTestLikelihood = GetFromManager(fitMan->raw()["General"]["Fitter"]["FitTestLikelihood"], false); } // ************************* @@ -113,12 +113,12 @@ void FitterBase::SaveSettings() { if (std::getenv("MaCh3_ROOT") == NULL) { MACH3LOG_ERROR("Need MaCh3_ROOT environment variable"); MACH3LOG_ERROR("Please remember about source bin/setup.MaCh3.sh"); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } if (std::getenv("MACH3") == NULL) { MACH3LOG_ERROR("Need MACH3 environment variable"); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } std::string header_path = std::string(std::getenv("MACH3")); @@ -176,8 +176,7 @@ void FitterBase::PrepareOutput() { // Check that we have added samples if (!samples.size()) { MACH3LOG_CRITICAL("No samples Found! If this is what you want find me here"); - MACH3LOG_CRITICAL("{}:{}", __FILE__, __LINE__); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } // Prepare the output trees diff --git a/mcmc/LikelihoodFit.cpp b/mcmc/LikelihoodFit.cpp index 96d18022..866dbc45 100644 --- a/mcmc/LikelihoodFit.cpp +++ b/mcmc/LikelihoodFit.cpp @@ -6,7 +6,7 @@ LikelihoodFit::LikelihoodFit(manager *man) : FitterBase(man) { // ******************* NPars = 0; NParsPCA = 0; - fMirroring = true; + fMirroring = GetFromManager(fitMan->raw()["General"]["Fitter"]["Mirroring"], false); } diff --git a/mcmc/MCMCProcessor.cpp b/mcmc/MCMCProcessor.cpp index 2e7e27c2..7540a09c 100644 --- a/mcmc/MCMCProcessor.cpp +++ b/mcmc/MCMCProcessor.cpp @@ -2,7 +2,6 @@ #include "TChain.h" - //Only if GPU is enabled #ifdef CUDA #include "mcmc/gpuMCMCProcessorUtils.cuh" @@ -98,7 +97,6 @@ MCMCProcessor::MCMCProcessor(const std::string &InputFile) : // The destructor MCMCProcessor::~MCMCProcessor() { // **************************** - // Close the pdf file MACH3LOG_INFO("Closing pdf in MCMCProcessor: {}", CanvasName); CanvasName += "]"; @@ -210,7 +208,7 @@ void MCMCProcessor::MakeOutputFile() { const int uniform = int(rand->Uniform(0, 10000)); // Open a TCanvas to write the posterior onto Posterior = new TCanvas(("Posterior" + std::to_string(uniform)).c_str(), ("Posterior" + std::to_string(uniform)).c_str(), 0, 0, 1024, 1024); - //KS: No idea why but ROOT changed treatment of viilin in R6. If you have non uniform binning this will results in very hard to see violin plots. + //KS: No idea why but ROOT changed treatment of violin in R6. If you have non uniform binning this will results in very hard to see violin plots. TCandle::SetScaledViolin(false); Posterior->SetGrid(); @@ -299,7 +297,7 @@ void MCMCProcessor::MakePostfit() { (*Means)(i) = Mean; (*Errors)(i) = Err; - GetGaussian(hpost[i], Mean, Err); + GetGaussian(hpost[i], Gauss, Mean, Err); (*Means_Gauss)(i) = Mean; (*Errors_Gauss)(i) = Err; @@ -679,7 +677,6 @@ void MCMCProcessor::DrawPostfit() { Start += NDSamplesBins[i]; } } - delete prefit; delete paramPlot; delete paramPlot_Gauss; @@ -694,7 +691,7 @@ void MCMCProcessor::DrawPostfit() { // Make fancy Credible Intervals plots void MCMCProcessor::MakeCredibleIntervals(const std::vector& CredibleIntervals, const std::vector& CredibleIntervalsColours, - bool CredibleInSigmas) { + const bool CredibleInSigmas) { // ********************* if(hpost[0] == nullptr) MakePostfit(); @@ -707,7 +704,7 @@ void MCMCProcessor::MakeCredibleIntervals(const std::vector& CredibleInt if(CredibleIntervals.size() != CredibleIntervalsColours.size()) { MACH3LOG_ERROR("Size of CredibleIntervals is not equal to size of CredibleIntervalsColours"); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } if(CredibleIntervals.size() > 1) @@ -1114,7 +1111,6 @@ void MCMCProcessor::MakeCovariance() { Correlation->Write("Correlation"); } - // *************** //KS: Cache all steps to allow multithreading, hit RAM quite a bit void MCMCProcessor::CacheSteps() { @@ -1128,7 +1124,7 @@ void MCMCProcessor::CacheSteps() { MACH3LOG_ERROR("It look like ParStep was already filled "); MACH3LOG_ERROR("Even though it is used for MakeCovariance_MP and for DiagMCMC "); MACH3LOG_ERROR("it has different structure in both for cache hits, sorry "); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } MACH3LOG_INFO("Caching input tree..."); @@ -1297,9 +1293,10 @@ void MCMCProcessor::MakeCovariance_MP(bool Mute) { }// End j loop }// End i loop - if(!Mute) clock.Stop(); - if(!Mute) MACH3LOG_INFO("Making Covariance took {:.2f}s to finish for {} steps", clock.RealTime(), nEntries); - + if(!Mute) { + clock.Stop(); + MACH3LOG_INFO("Making Covariance took {:.2f}s to finish for {} steps", clock.RealTime(), nEntries); + } OutputFile->cd(); if(printToPDF) { @@ -1327,15 +1324,17 @@ void MCMCProcessor::MakeCovariance_MP(bool Mute) { }// End j loop }// End i loop } //end if pdf - if(!Mute) Covariance->Write("Covariance"); - if(!Mute) Correlation->Write("Correlation"); + if(!Mute) { + Covariance->Write("Covariance"); + Correlation->Write("Correlation"); + } } // ********************* // Based on https://www.jstor.org/stable/25651249?seq=3, // all credits for finding and studying it goes to Henry -void MCMCProcessor::MakeSubOptimality(int NIntervals) { +void MCMCProcessor::MakeSubOptimality(const int NIntervals) { // ********************* //Save burn in cut, at the end of the loop we will return to default values @@ -1670,7 +1669,7 @@ void MCMCProcessor::DrawCorrelations1D() { void MCMCProcessor::MakeCredibleRegions(const std::vector& CredibleRegions, const std::vector& CredibleRegionStyle, const std::vector& CredibleRegionColor, - bool CredibleInSigmas) { + const bool CredibleInSigmas) { // ********************* if(hpost2D == nullptr) MakeCovariance_MP(); @@ -1813,7 +1812,6 @@ void MCMCProcessor::MakeCredibleRegions(const std::vector& CredibleRegio delete[] hpost_2D_cl; } - // ********************* // Make fancy triangle plot for selected parameters void MCMCProcessor::MakeTrianglePlot(const std::vector& ParNames, @@ -1825,7 +1823,7 @@ void MCMCProcessor::MakeTrianglePlot(const std::vector& ParNames, const std::vector& CredibleRegionStyle, const std::vector& CredibleRegionColor, // Other - bool CredibleInSigmas) { + const bool CredibleInSigmas) { // ********************* if(hpost2D == nullptr) MakeCovariance_MP(); @@ -2461,7 +2459,7 @@ void MCMCProcessor::FindInputFiles() { if (Config == nullptr) { MACH3LOG_ERROR("Didn't find MaCh3_Config tree in MCMC file! {}", MCMCFile); TempFile->ls(); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } //KS:Most inputs are in ${MACH3}/inputs/blarb.root if (std::getenv("MACH3") != nullptr) { @@ -2592,7 +2590,7 @@ void MCMCProcessor::ReadNDFile() { TFile *NDdetFile = new TFile(CovPos[kNDPar].back().c_str(), "open"); if (NDdetFile->IsZombie()) { MACH3LOG_ERROR("Couldn't find NDdetFile {}", CovPos[kNDPar].back()); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } NDdetFile->cd(); @@ -2747,243 +2745,6 @@ void MCMCProcessor::SetStepCut(const int Cuts) { BurnInCut = Cuts; } -// ************************** -//CW: Get the mean and RMS of a 1D posterior -void MCMCProcessor::GetArithmetic(TH1D * const hist, double& Mean, double& Error) { -// ************************** - Mean = hist->GetMean(); - Error = hist->GetRMS(); -} - -// ************************** -//CW: Get Gaussian characteristics -void MCMCProcessor::GetGaussian(TH1D *& hist, double& Mean, double& Error) { -// ************************** - const double meanval = hist->GetMean(); - const double err = hist->GetRMS(); - const double peakval = hist->GetBinCenter(hist->GetMaximumBin()); - - // Set the range for the Gaussian fit - Gauss->SetRange(meanval - 1.5*err , meanval + 1.5*err); - // Set the starting parameters close to RMS and peaks of the histograms - Gauss->SetParameters(hist->GetMaximum()*err*std::sqrt(2*3.14), peakval, err); - - // Perform the fit - hist->Fit(Gauss->GetName(),"Rq"); - hist->SetStats(0); - - Mean = Gauss->GetParameter(1); - Error = Gauss->GetParameter(2); -} - -// *************** -//CW: Get the highest posterior density from a TH1D -void MCMCProcessor::GetHPD(TH1D* const hist, double& Mean, double& Error, double& Error_p, double& Error_m, const double coverage) { -// *************** - // Get the bin which has the largest posterior density - const int MaxBin = hist->GetMaximumBin(); - // And it's value - const double peakval = hist->GetBinCenter(MaxBin); - - // The total integral of the posterior - const long double Integral = hist->Integral(); - //KS: and integral of left handed and right handed parts - const long double LowIntegral = hist->Integral(1, MaxBin-1) + hist->GetBinContent(MaxBin)/2.0; - const long double HighIntegral = hist->Integral(MaxBin+1, hist->GetNbinsX()) + hist->GetBinContent(MaxBin)/2.0; - - // Keep count of how much area we're covering - //KS: Take only half content of HPD bin as one half goes for right handed error and the other for left handed error - long double sum = hist->GetBinContent(MaxBin)/2.0; - - // Counter for current bin - int CurrBin = MaxBin; - while (sum/HighIntegral < coverage && CurrBin < hist->GetNbinsX()) { - CurrBin++; - sum += hist->GetBinContent(CurrBin); - } - const double sigma_p = std::fabs(hist->GetBinCenter(MaxBin)-hist->GetXaxis()->GetBinUpEdge(CurrBin)); - // Reset the sum - //KS: Take only half content of HPD bin as one half goes for right handed error and the other for left handed error - sum = hist->GetBinContent(MaxBin)/2.0; - - // Reset the bin counter - CurrBin = MaxBin; - // Counter for current bin - while (sum/LowIntegral < coverage && CurrBin > 1) { - CurrBin--; - sum += hist->GetBinContent(CurrBin); - } - const double sigma_m = std::fabs(hist->GetBinCenter(CurrBin)-hist->GetBinLowEdge(MaxBin)); - - // Now do the double sided HPD - //KS: Start sum from the HPD - sum = hist->GetBinContent(MaxBin); - int LowBin = MaxBin; - int HighBin = MaxBin; - long double LowCon = 0.0; - long double HighCon = 0.0; - - while (sum/Integral < coverage && (LowBin > 0 || HighBin < hist->GetNbinsX()+1)) - { - LowCon = 0.0; - HighCon = 0.0; - //KS:: Move further only if you haven't reached histogram end - if(LowBin > 1) - { - LowBin--; - LowCon = hist->GetBinContent(LowBin); - } - if(HighBin < hist->GetNbinsX()) - { - HighBin++; - HighCon = hist->GetBinContent(HighBin); - } - - // If we're on the last slice and the lower contour is larger than the upper - if ((sum+LowCon+HighCon)/Integral > coverage && LowCon > HighCon) { - sum += LowCon; - break; - // If we're on the last slice and the upper contour is larger than the lower - } else if ((sum+LowCon+HighCon)/Integral > coverage && HighCon >= LowCon) { - sum += HighCon; - break; - } else { - sum += LowCon + HighCon; - } - } - - double sigma_hpd = 0.0; - if (LowCon > HighCon) { - sigma_hpd = std::fabs(hist->GetBinLowEdge(LowBin)-hist->GetBinCenter(MaxBin)); - } else { - sigma_hpd = std::fabs(hist->GetXaxis()->GetBinUpEdge(HighBin)-hist->GetBinCenter(MaxBin)); - } - - Mean = peakval; - Error = sigma_hpd; - Error_p = sigma_p; - Error_m = sigma_m; -} - -// *************** -//KS: Get 1D histogram within credible interval, hpost_copy has to have the same binning, I don't do Copy() as this will lead to problems if this is used under multithreading -void MCMCProcessor::GetCredibleInterval(TH1D* const hist, TH1D* hpost_copy, const double coverage) { -// *************** - if(coverage > 1) - { - MACH3LOG_ERROR("Specified Credible Interval is greater that 1 and equal to {} Should be between 0 and 1", coverage); - throw MaCh3Exception(__FILE__ , __LINE__ ); - } - //KS: Reset first copy of histogram - hpost_copy->Reset(""); - hpost_copy->Fill(0.0, 0.0); - - //KS: Temporary structure to be thread save - double *hist_copy = new double[hist->GetXaxis()->GetNbins()+1]; - bool *hist_copy_fill = new bool[hist->GetXaxis()->GetNbins()+1]; - for (int i = 0; i <= hist->GetXaxis()->GetNbins(); ++i) - { - hist_copy[i] = hist->GetBinContent(i); - hist_copy_fill[i] = false; - } - - /// Loop over histogram bins with highest number of entries until covered 90 or 68.3% - const long double Integral = hist->Integral(); - long double sum = 0; - - while ((sum / Integral) < coverage) - { - /// Get bin of highest content and save the number of entries reached so far - int max_entry_bin = 0; - double max_entries = 0.; - for (int i = 0; i <= hist->GetXaxis()->GetNbins(); ++i) - { - if (hist_copy[i] > max_entries) - { - max_entries = hist_copy[i]; - max_entry_bin = i; - } - } - /// Replace bin value by -1 so it is not looped over as being maximum bin again - hist_copy[max_entry_bin] = -1.; - hist_copy_fill[max_entry_bin] = true; - - sum += max_entries; - } - //KS: Now fill our copy only for bins which got included in coverage region - for(int i = 0; i <= hist->GetXaxis()->GetNbins(); ++i) - { - if(hist_copy_fill[i]) hpost_copy->SetBinContent(i, hist->GetBinContent(i)); - } - - delete[] hist_copy; - delete[] hist_copy_fill; - - return; -} - -// *************** -//KS: Set 2D contour within some coverage -void MCMCProcessor::GetCredibleRegion(TH2D* const hist2D, const double coverage) { -// *************** - if(coverage > 1) - { - MACH3LOG_ERROR("Specified Credible Region is greater than 1 and equal to {:.2f} Should be between 0 and 1 {}", coverage); - throw MaCh3Exception(__FILE__ , __LINE__ ); - } - - //KS: Temporary structure to be thread save - double **hist_copy = new double*[hist2D->GetXaxis()->GetNbins()+1]; - for (int i = 0; i <= hist2D->GetXaxis()->GetNbins(); ++i) - { - hist_copy[i] = new double[hist2D->GetYaxis()->GetNbins()+1]; - for (int j = 0; j <= hist2D->GetYaxis()->GetNbins(); ++j) - { - hist_copy[i][j] = hist2D->GetBinContent(i,j); - } - } - - /// Loop over histogram bins with highest number of entries until covered 90 or 68.3% - const long double Integral = hist2D->Integral(); - long double sum = 0; - - //We need to as ROOT requires array to set to contour - double Contour[1]; - while ((sum / Integral) < coverage) - { - /// Get bin of highest content and save the number of entries reached so far - int max_entry_bin_x = 0; - int max_entry_bin_y = 0; - double max_entries = 0.; - for (int i = 0; i <= hist2D->GetXaxis()->GetNbins(); ++i) - { - for (int j = 0; j <= hist2D->GetYaxis()->GetNbins(); ++j) - { - if (hist_copy[i][j] > max_entries) - { - max_entries = hist_copy[i][j]; - max_entry_bin_x = i; - max_entry_bin_y = j; - } - } - } - /// Replace bin value by -1 so it is not looped over as being maximum bin again - hist_copy[max_entry_bin_x][max_entry_bin_y] = -1.; - - sum += max_entries; - Contour[0] = max_entries; - } - hist2D->SetContour(1, Contour); - - //Delete temporary arrays - for (int i = 0; i <= hist2D->GetXaxis()->GetNbins(); ++i) - { - delete[] hist_copy[i]; - } - delete[] hist_copy; - - return; -} // *************** // Pass central value @@ -3472,13 +3233,13 @@ void MCMCProcessor::PrepareDiagMCMC() { MACH3LOG_ERROR("It look like ParStep was already filled "); MACH3LOG_ERROR("Even though it is used for MakeCovariance_MP and for DiagMCMC"); MACH3LOG_ERROR("it has different structure in both for cache hits, sorry "); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } if(nBatches == 0) { MACH3LOG_ERROR("nBatches is equal to 0"); MACH3LOG_ERROR("please use SetnBatches to set other value fore example 20"); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } // Initialise ParStep @@ -4145,7 +3906,7 @@ void MCMCProcessor::BatchedAnalysis() { { MACH3LOG_ERROR("BatchedAverages haven't been initialises or have been deleted something is wrong"); MACH3LOG_ERROR("I need it and refuse to go further"); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } // Calculate variance estimator using batched means following https://arxiv.org/pdf/1911.00915.pdf see Eq. 1.2 diff --git a/mcmc/MCMCProcessor.h b/mcmc/MCMCProcessor.h index 0646cea3..a8858931 100644 --- a/mcmc/MCMCProcessor.h +++ b/mcmc/MCMCProcessor.h @@ -94,7 +94,7 @@ class MCMCProcessor { /// @param CredibleInSigmas Bool telling whether intervals are in percentage or in sigmas, then special conversions is used void MakeCredibleIntervals(const std::vector& CredibleIntervals = {0.99, 0.90, 0.68 }, const std::vector& CredibleIntervalsColours = {kCyan+4, kCyan-2, kCyan-10}, - bool CredibleInSigmas = false + const bool CredibleInSigmas = false ); /// @brief Draw the post-fit covariances void DrawCovariance(); @@ -106,7 +106,7 @@ class MCMCProcessor { void MakeCredibleRegions(const std::vector& CredibleRegions = {0.99, 0.90, 0.68}, const std::vector& CredibleRegionStyle = {kDashed, kSolid, kDotted}, const std::vector& CredibleRegionColor = {kGreen-3, kGreen-10, kGreen}, - bool CredibleInSigmas = false + const bool CredibleInSigmas = false ); /// @brief Make fancy triangle plot for selected parameters /// @param CredibleIntervals Vector with values of credible intervals, must be in descending order @@ -125,7 +125,7 @@ class MCMCProcessor { const std::vector& CredibleRegionStyle = {kDashed, kSolid, kDotted}, const std::vector& CredibleRegionColor = {kGreen-3, kGreen-10, kGreen}, // Other - bool CredibleInSigmas = false + const bool CredibleInSigmas = false ); /// @brief Make funny polar plot /// @param ParNames Vector with parameter names for which Polar Plot will be made @@ -266,34 +266,6 @@ class MCMCProcessor { /// @brief Prepare all objects used for output inline void SetupOutput(); - //Analyse posterior distribution - /// @brief Get Arithmetic mean from posterior - /// @param hist histograms from which we extract arithmetic mean - /// @param Mean Arithmetic Mean value - /// @param Error Arithmetic Error value - inline void GetArithmetic(TH1D * const hist, double& Mean, double& Error); - /// @brief Fit Gaussian to posterior - /// @param hist histograms to which we fit gaussian - /// @param Mean Gaussian Mean value - /// @param Error Gaussian Error value - inline void GetGaussian(TH1D *& hist, double& Mean, double& Error); - /// @brief Get Highest Posterior Density (HPD) - /// @param hist histograms from which we HPD - /// @param Mean HPD Mean value - /// @param Error HPD Error value - /// @param Error_p HPD Negative (left hand side) Error value - /// @param Error_m HPD Positive (right hand side) Error value - /// @param coverage What is defined coverage, by default 0.6827 (1 sigma) - inline void GetHPD(TH1D * const hist, double& Mean, double& Error, double& Error_p, double& Error_m, const double coverage = 0.6827); - /// @brief Get 1D Credible Interval - /// @param hist histograms based on which we calculate credible interval - /// @param coverage What is defined coverage, by default 0.6827 (1 sigma) - inline void GetCredibleInterval(TH1D* const hist, TH1D* hpost_copy, const double coverage = 0.6827); - /// @brief Get 2D Credible Region - /// @param hist2D histograms based on which we calculate credible regions - /// @param coverage What is defined coverage, by default 0.6827 (1 sigma) - inline void GetCredibleRegion(TH2D* hist2D, const double coverage = 0.6827); - // MCMC Diagnostic /// @brief CW: Prepare branches etc. for DiagMCMC inline void PrepareDiagMCMC(); diff --git a/mcmc/MinuitFit.cpp b/mcmc/MinuitFit.cpp index 93327ee0..a9cb8417 100644 --- a/mcmc/MinuitFit.cpp +++ b/mcmc/MinuitFit.cpp @@ -15,7 +15,7 @@ MinuitFit::MinuitFit(manager *man) : LikelihoodFit(man) { // Destructor: close the logger and output file MinuitFit::~MinuitFit() { // ************************* - if(minuit != NULL) delete minuit; + if(minuit != nullptr) delete minuit; } diff --git a/mcmc/StatisticalUtils.h b/mcmc/StatisticalUtils.h index beac09f7..5c19f67b 100644 --- a/mcmc/StatisticalUtils.h +++ b/mcmc/StatisticalUtils.h @@ -14,6 +14,7 @@ // ************************** /// @brief KS: Following H. Jeffreys. The theory of probability. UOP Oxford, 1998. DOI: 10.2307/3619118. +/// @param BayesFactor Obtained value of Bayes factor inline std::string GetJeffreysScale(const double BayesFactor){ // ************************** std::string JeffreysScale = ""; @@ -31,6 +32,7 @@ inline std::string GetJeffreysScale(const double BayesFactor){ // ************************** /// @brief KS: Based on Table 1 in https://www.t2k.org/docs/technotes/435 +/// @param BayesFactor Obtained value of Bayes factor inline std::string GetDunneKaboth(const double BayesFactor){ // ************************** std::string DunneKaboth = ""; @@ -73,7 +75,7 @@ inline double GetSigmaValue(const int sigma) { break; default: MACH3LOG_ERROR("{} is unsupported value of sigma", sigma); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); break; } return width; @@ -88,8 +90,7 @@ inline double GetBIC(const double llh, const int data, const int nPars){ if(nPars == 0) { MACH3LOG_ERROR("You haven't passed number of model parameters as it is still zero"); - std::cerr << __FILE__ << ":" << __LINE__ << std::endl; - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } const double BIC = nPars * logl(data) + llh; @@ -109,7 +110,12 @@ inline double GetNeffective(const int N1, const int N2) { // **************** /// @brief KS: For more see https://www.t2k.org/docs/technotes/429/TN429_v8#page=63 -inline void CheckBonferoniCorrectedpValue(std::vector SampleNameVec, std::vector PValVec, double Threshold = 0.05) { +/// @param SampleNameVec vector of sample names +/// @param PValVec pvalue for each sample +/// @param Threshold pvalue accepted threshold, usually 5% +inline void CheckBonferoniCorrectedpValue(const std::vector& SampleNameVec, + const std::vector& PValVec, + const double Threshold = 0.05) { // **************** std::cout<<""< SampleNameVec if(SampleNameVec.size() != PValVec.size()) { MACH3LOG_ERROR("Size of vectors do not match"); - throw; + throw MaCh3Exception(__FILE__ , __LINE__ ); } - const int NumberOfStatisticalTests = SampleNameVec.size(); //KS: 0.05 or 5% is value used by T2K. const double StatisticalSignificanceDown = Threshold / NumberOfStatisticalTests; @@ -142,20 +147,20 @@ inline void CheckBonferoniCorrectedpValue(std::vector SampleNameVec } } - if(Counter == 0) - { - std::cout<<" Every sample passed Bonferroni-corrected statistical significance level test"< GroupClasifier) { +inline int GetNumberOfRuns(const std::vector& GroupClasifier) { // **************** int NumberOfRuns = 0; @@ -185,6 +190,11 @@ inline int GetNumberOfRuns(std::vector GroupClasifier) { } // **************** +/// @brief KS: Calculate Beta parameter which will be different based on specified test statistic +/// @param data Number of data events in a bin +/// @param mc Number of MC events in a bin +/// @param w2 Value of weight squared in a bin +/// @param TestStat Test statistic based on which we calculate beta inline double GetBetaParameter(const double data, const double mc, const double w2, TestStatistic TestStat) { // **************** double Beta = 0.0; @@ -204,9 +214,8 @@ inline double GetBetaParameter(const double data, const double mc, const double // CW: b^2 - 4ac in quadratic equation const double temp2 = temp*temp + 4*data*fractional*fractional; if (temp2 < 0) { - std::cerr << "Negative square root in Barlow Beeston coefficient calculation!" << std::endl; - std::cerr << __FILE__ << ":" << __LINE__ << std::endl; - throw; + MACH3LOG_ERROR("Negative square root in Barlow Beeston coefficient calculation!"); + throw MaCh3Exception(__FILE__ , __LINE__ ); } // CW: Solve for the positive beta Beta = (-1*temp+std::sqrt(temp2))/2.; @@ -217,7 +226,8 @@ inline double GetBetaParameter(const double data, const double mc, const double // ********************* /// @brief Based on http://probability.ca/jeff/ftpdir/adaptex.pdf -inline double GetSubOptimality(std::vector EigenValues, const int TotalTarameters) { +/// @param EigenValues Eigen values of covariance matrix +inline double GetSubOptimality(const std::vector& EigenValues, const int TotalTarameters) { // ********************* double sum_eigenvalues_squared_inv = 0.0; double sum_eigenvalues_inv = 0.0; @@ -232,9 +242,268 @@ inline double GetSubOptimality(std::vector EigenValues, const int TotalT return SubOptimality; } + +// ************************** +/// @brief CW: Get Arithmetic mean from posterior +/// @param hist histograms from which we extract arithmetic mean +/// @param Mean Arithmetic Mean value +/// @param Error Arithmetic Error value +inline void GetArithmetic(TH1D * const hist, double& Mean, double& Error) { +// ************************** + Mean = hist->GetMean(); + Error = hist->GetRMS(); +} + +// ************************** +/// @brief CW: Fit Gaussian to posterior +/// @param hist histograms to which we fit gaussian +/// @param gauss tf1 with gaussian, we pass pointer to make things faster +/// @param Mean Gaussian Mean value +/// @param Error Gaussian Error value +inline void GetGaussian(TH1D*& hist, TF1* gauss, double& Mean, double& Error) { +// ************************** + const double meanval = hist->GetMean(); + const double err = hist->GetRMS(); + const double peakval = hist->GetBinCenter(hist->GetMaximumBin()); + + // Set the range for the Gaussian fit + gauss->SetRange(meanval - 1.5*err , meanval + 1.5*err); + // Set the starting parameters close to RMS and peaks of the histograms + gauss->SetParameters(hist->GetMaximum()*err*std::sqrt(2*3.14), peakval, err); + + // Perform the fit + hist->Fit(gauss->GetName(),"Rq"); + hist->SetStats(0); + + Mean = gauss->GetParameter(1); + Error = gauss->GetParameter(2); +} + +// *************** +/// @brief Get Highest Posterior Density (HPD) +/// @param hist histograms from which we HPD +/// @param Mean HPD Mean value +/// @param Error HPD Error value +/// @param Error_p HPD Negative (left hand side) Error value +/// @param Error_m HPD Positive (right hand side) Error value +/// @param coverage What is defined coverage, by default 0.6827 (1 sigma) +inline void GetHPD(TH1D* const hist, double& Mean, double& Error, double& Error_p, double& Error_m, const double coverage = 0.6827) { +// **************** + // Get the bin which has the largest posterior density + const int MaxBin = hist->GetMaximumBin(); + // And it's value + const double peakval = hist->GetBinCenter(MaxBin); + + // The total integral of the posterior + const long double Integral = hist->Integral(); + //KS: and integral of left handed and right handed parts + const long double LowIntegral = hist->Integral(1, MaxBin-1) + hist->GetBinContent(MaxBin)/2.0; + const long double HighIntegral = hist->Integral(MaxBin+1, hist->GetNbinsX()) + hist->GetBinContent(MaxBin)/2.0; + + // Keep count of how much area we're covering + //KS: Take only half content of HPD bin as one half goes for right handed error and the other for left handed error + long double sum = hist->GetBinContent(MaxBin)/2.0; + + // Counter for current bin + int CurrBin = MaxBin; + while (sum/HighIntegral < coverage && CurrBin < hist->GetNbinsX()) { + CurrBin++; + sum += hist->GetBinContent(CurrBin); + } + const double sigma_p = std::fabs(hist->GetBinCenter(MaxBin)-hist->GetXaxis()->GetBinUpEdge(CurrBin)); + // Reset the sum + //KS: Take only half content of HPD bin as one half goes for right handed error and the other for left handed error + sum = hist->GetBinContent(MaxBin)/2.0; + + // Reset the bin counter + CurrBin = MaxBin; + // Counter for current bin + while (sum/LowIntegral < coverage && CurrBin > 1) { + CurrBin--; + sum += hist->GetBinContent(CurrBin); + } + const double sigma_m = std::fabs(hist->GetBinCenter(CurrBin)-hist->GetBinLowEdge(MaxBin)); + + // Now do the double sided HPD + //KS: Start sum from the HPD + sum = hist->GetBinContent(MaxBin); + int LowBin = MaxBin; + int HighBin = MaxBin; + long double LowCon = 0.0; + long double HighCon = 0.0; + + while (sum/Integral < coverage && (LowBin > 0 || HighBin < hist->GetNbinsX()+1)) + { + LowCon = 0.0; + HighCon = 0.0; + //KS:: Move further only if you haven't reached histogram end + if(LowBin > 1) + { + LowBin--; + LowCon = hist->GetBinContent(LowBin); + } + if(HighBin < hist->GetNbinsX()) + { + HighBin++; + HighCon = hist->GetBinContent(HighBin); + } + + // If we're on the last slice and the lower contour is larger than the upper + if ((sum+LowCon+HighCon)/Integral > coverage && LowCon > HighCon) { + sum += LowCon; + break; + // If we're on the last slice and the upper contour is larger than the lower + } else if ((sum+LowCon+HighCon)/Integral > coverage && HighCon >= LowCon) { + sum += HighCon; + break; + } else { + sum += LowCon + HighCon; + } + } + + double sigma_hpd = 0.0; + if (LowCon > HighCon) { + sigma_hpd = std::fabs(hist->GetBinLowEdge(LowBin)-hist->GetBinCenter(MaxBin)); + } else { + sigma_hpd = std::fabs(hist->GetXaxis()->GetBinUpEdge(HighBin)-hist->GetBinCenter(MaxBin)); + } + + Mean = peakval; + Error = sigma_hpd; + Error_p = sigma_p; + Error_m = sigma_m; +} + + + +// *************** +/// @brief KS: Get 1D histogram within credible interval, hpost_copy has to have the same binning, I don't do Copy() as this will lead to problems if this is used under multithreading +/// @param hist histograms based on which we calculate credible interval +/// @param coverage What is defined coverage, by default 0.6827 (1 sigma) +inline void GetCredibleInterval(TH1D* const hist, TH1D* hpost_copy, const double coverage = 0.6827) { +// *************** + if(coverage > 1) + { + MACH3LOG_ERROR("Specified Credible Interval is greater that 1 and equal to {} Should be between 0 and 1", coverage); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + //KS: Reset first copy of histogram + hpost_copy->Reset(""); + hpost_copy->Fill(0.0, 0.0); + + //KS: Temporary structure to be thread save + double *hist_copy = new double[hist->GetXaxis()->GetNbins()+1]; + bool *hist_copy_fill = new bool[hist->GetXaxis()->GetNbins()+1]; + for (int i = 0; i <= hist->GetXaxis()->GetNbins(); ++i) + { + hist_copy[i] = hist->GetBinContent(i); + hist_copy_fill[i] = false; + } + + /// Loop over histogram bins with highest number of entries until covered 90 or 68.3% + const long double Integral = hist->Integral(); + long double sum = 0; + + while ((sum / Integral) < coverage) + { + /// Get bin of highest content and save the number of entries reached so far + int max_entry_bin = 0; + double max_entries = 0.; + for (int i = 0; i <= hist->GetXaxis()->GetNbins(); ++i) + { + if (hist_copy[i] > max_entries) + { + max_entries = hist_copy[i]; + max_entry_bin = i; + } + } + /// Replace bin value by -1 so it is not looped over as being maximum bin again + hist_copy[max_entry_bin] = -1.; + hist_copy_fill[max_entry_bin] = true; + + sum += max_entries; + } + //KS: Now fill our copy only for bins which got included in coverage region + for(int i = 0; i <= hist->GetXaxis()->GetNbins(); ++i) + { + if(hist_copy_fill[i]) hpost_copy->SetBinContent(i, hist->GetBinContent(i)); + } + + delete[] hist_copy; + delete[] hist_copy_fill; + + return; +} + +// *************** +/// @brief KS: Set 2D contour within some coverage +/// @param hist2D histograms based on which we calculate credible regions +/// @param coverage What is defined coverage, by default 0.6827 (1 sigma) +inline void GetCredibleRegion(TH2D* const hist2D, const double coverage = 0.6827) { +// *************** + if(coverage > 1) + { + MACH3LOG_ERROR("Specified Credible Region is greater than 1 and equal to {:.2f} Should be between 0 and 1 {}", coverage); + throw MaCh3Exception(__FILE__ , __LINE__ ); + } + + //KS: Temporary structure to be thread save + double **hist_copy = new double*[hist2D->GetXaxis()->GetNbins()+1]; + for (int i = 0; i <= hist2D->GetXaxis()->GetNbins(); ++i) + { + hist_copy[i] = new double[hist2D->GetYaxis()->GetNbins()+1]; + for (int j = 0; j <= hist2D->GetYaxis()->GetNbins(); ++j) + { + hist_copy[i][j] = hist2D->GetBinContent(i,j); + } + } + + /// Loop over histogram bins with highest number of entries until covered 90 or 68.3% + const long double Integral = hist2D->Integral(); + long double sum = 0; + + //We need to as ROOT requires array to set to contour + double Contour[1]; + while ((sum / Integral) < coverage) + { + /// Get bin of highest content and save the number of entries reached so far + int max_entry_bin_x = 0; + int max_entry_bin_y = 0; + double max_entries = 0.; + for (int i = 0; i <= hist2D->GetXaxis()->GetNbins(); ++i) + { + for (int j = 0; j <= hist2D->GetYaxis()->GetNbins(); ++j) + { + if (hist_copy[i][j] > max_entries) + { + max_entries = hist_copy[i][j]; + max_entry_bin_x = i; + max_entry_bin_y = j; + } + } + } + /// Replace bin value by -1 so it is not looped over as being maximum bin again + hist_copy[max_entry_bin_x][max_entry_bin_y] = -1.; + + sum += max_entries; + Contour[0] = max_entries; + } + hist2D->SetContour(1, Contour); + + //Delete temporary arrays + for (int i = 0; i <= hist2D->GetXaxis()->GetNbins(); ++i) + { + delete[] hist_copy[i]; + } + delete[] hist_copy; + + return; +} + // ********************* -// Interquartile Range (IQR) -inline double GetIQR(TH1D *Hist){ +/// @brief Interquartile Range (IQR) +/// @param hist histograms from which we IQR +inline double GetIQR(TH1D *Hist) { // ********************* if(Hist->Integral() == 0) return 0.0; diff --git a/samplePDF/Structs.cpp b/samplePDF/Structs.cpp index afc00156..84c90e26 100644 --- a/samplePDF/Structs.cpp +++ b/samplePDF/Structs.cpp @@ -434,8 +434,8 @@ double PolyIntegralWidth(TH2Poly *Histogram) { // ********************* //KS: Remove fitted TF1 from hist to make comparison easier -void RemoveFitter(TH1D* hist, std::string name) { - // ********************* +void RemoveFitter(TH1D* hist, const std::string& name) { +// ********************* TList *listOfFunctions = hist->GetListOfFunctions(); TF1 *fitter = dynamic_cast(listOfFunctions->FindObject(name.c_str())); diff --git a/samplePDF/Structs.h b/samplePDF/Structs.h index 4005d54c..b08eab8f 100644 --- a/samplePDF/Structs.h +++ b/samplePDF/Structs.h @@ -4,8 +4,11 @@ /// N.B. for 64 bit systems sizeof(float) == sizeof(double) so not a huge effect /// KS: Need more testing on FD #ifdef _LOW_MEMORY_STRUCTS_ +/// Custom floating point (float or double) #define _float_ float +/// Custom integer (int or short int) #define _int_ short int +/// Custom unsigned integer (unsigned short int or unsigned int) #define _unsigned_int_ unsigned short int #else #define _float_ double @@ -148,10 +151,10 @@ inline std::string SplineInterpolation_ToString(const SplineInterpolation i) { /// Make an enum of systematic type recognised by covariance class enum SystType { - kNorm, - kSpline, - kFunc, - kSystTypes //This only enumerates + kNorm, //!< For normalisation parameters + kSpline, //!< For splined parameters (1D) + kFunc, //!< For functional parameters + kSystTypes //!< This only enumerates }; @@ -192,7 +195,6 @@ inline std::string SystType_ToString(const SystType i) { return name; } - // *************************** // A handy namespace for variables extraction namespace MaCh3Utils { @@ -212,10 +214,10 @@ namespace MaCh3Utils { /// Enum to track the target material enum TargetMat { // ***************** - kTarget_H = 1, - kTarget_C = 12, + kTarget_H = 1, //!< Hydrogen + kTarget_C = 12, //!< Carbon 12 kTarget_N = 14, - kTarget_O = 16, + kTarget_O = 16, //!< Oxygen 16 kTarget_Al = 27, kTarget_Ar = 40, kTarget_Ti = 48, @@ -351,12 +353,12 @@ inline int ProbsToPDG(ProbNu NuType){ /// Make an enum of the test statistic that we're using enum TestStatistic { - kPoisson, - kBarlowBeeston, - kIceCube, - kPearson, - kDembinskiAbdelmottele, - kNTestStatistics //This only enumerates statistic + kPoisson, //!< Standard Poisson likelihood + kBarlowBeeston, //!< Barlow-Beeston following Conway https://cds.cern.ch/record/1333496? + kIceCube, //!< Based on https://arxiv.org/abs/1901.04645 + kPearson, //!< Standard Pearson likelihood + kDembinskiAbdelmottele, //!< Based on arXiv:2206.12346v2 + kNTestStatistics //!< This only enumerates statistic }; // ************************************************** @@ -391,10 +393,10 @@ inline std::string TestStatistic_ToString(TestStatistic i) { } /// @brief WP: Helper function for calculating unbinned Integral of TH2Poly i.e including overflow -double OverflowIntegral(TH2Poly*); +double OverflowIntegral(TH2Poly* poly); /// @brief WP: Helper function for calculating binned Integral of TH2Poly i.e not including overflow -double NoOverflowIntegral(TH2Poly*); +double NoOverflowIntegral(TH2Poly* poly); /// @brief WP: Poly Projectors TH1D* PolyProjectionX(TObject* poly, std::string TempName, std::vector xbins, bool computeErrors = false); @@ -414,13 +416,13 @@ TH2Poly* PolyScaleWidth(TH2Poly *Histogram, double scale); /// @brief WP: Helper to calc integral of th2poly analogous to th2d integra; with option "width" double PolyIntegralWidth(TH2Poly *Histogram); -/// @brief KS: Sanity check for TH2Poly +/// @brief 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 +/// @brief However if you use ROOT6 and have ROOT6 file exactly the same code will work. Something have changed with how TH2Poly bins are stored in TFile +/// @param file ROOT file that we will make version checks void CheckTH2PolyFileVersion(TFile *file); - - /// @brief KS: Remove fitted TF1 from hist to make comparison easier -void RemoveFitter(TH1D* hist, std::string name); +void RemoveFitter(TH1D* hist, const std::string& name); /// @brief Helper to check if files exist or not inline std::string file_exists(std::string filename) {