Merge pull request #136 from mach3-software/feature_FixSegFault

Fix seg fault
mach3-software · Oct 1, 2024 · e94eac3 · e94eac3
2 parents 34bb829 + fa7d227
commit e94eac3
Show file tree

Hide file tree

Showing 5 changed files with 127 additions and 128 deletions.
diff --git a/.github/workflows/CIValidations.yml b/.github/workflows/CIValidations.yml
@@ -22,37 +22,35 @@ jobs:
       matrix:
         include:
           - name: Spline Validations
-            test_1: Apps/SplineValidations
+            test_1: ./Apps/SplineValidations
             test_2: empty
             test_3: empty
             test_4: empty
             test_5: empty
             test_6: empty
 
           - name: Covariance Validations
-            test_1: Apps/CovarianceValidations
-            test_2: Apps/MaCh3ModeValidations
+            test_1: ./Apps/CovarianceValidations
+            test_2: ./Apps/MaCh3ModeValidations
             test_3: empty
             test_4: empty
             test_5: empty
             test_6: empty
 
           - name: Fitter Validations
-            test_1: Apps/FitterValidations
+            test_1: ./Apps/FitterValidations
             test_2: ./bin/ProcessMCMC bin/TutorialDiagConfig.yaml MCMC_Test.root
             test_3: ./bin/DiagMCMC MCMC_Test.root bin/TutorialDiagConfig.yaml
             test_4: ./bin/RHat 10 MCMC_Test.root MCMC_Test.root MCMC_Test.root MCMC_Test.root
             test_5: ./bin/CombineMaCh3Chains -o MergedChain.root MCMC_Test.root MCMC_Test.root MCMC_Test.root MCMC_Test.root
-            test_6: empty
+            test_6: ./bin/GetPenaltyTerm MCMC_Test.root bin/TutorialDiagConfig.yaml
             # TODO
-            # Stupid no idea why ROOT segfault when cleaning program...
-            #test_4: ./bin/GetPenaltyTerm MCMC_Test.root bin/TutorialDiagConfig.yaml
             #need fixing config
             #test_7: bin/GetPostfitParamPlots -o MCMC_Test_Process.root -c Inputs/PlottingConfig.yaml
             #test_8: bin/PlotLLH -o MCMC_Test.root -c Inputs/PlottingConfig.yaml
 
           - name: NuMCMC Tools Validations
-            test_1: Apps/NuMCMCvalidations.sh
+            test_1: ./Apps/NuMCMCvalidations.sh
             test_2: empty
             test_3: empty
             test_4: empty

diff --git a/Diagnostics/CombineMaCh3Chains.cpp b/Diagnostics/CombineMaCh3Chains.cpp
@@ -1,3 +1,4 @@
+// ROOT includes
 #include "TList.h"
 #include "TFile.h"
 #include "TMacro.h"
@@ -7,6 +8,7 @@
 #include "TKey.h"
 #include "TROOT.h"
 
+// MaCh3 includes
 #include "manager/manager.h"
 
 std::string OutFileName = "";
@@ -15,7 +17,7 @@ std::vector<std::string> inpFileList;
 bool forceOverwrite = false;
 
 /// EM: Will compare the version header contained in the two provided files and shout if they don't match
-bool checkSoftwareVersions(TFile *file, TFile *prevFile, std::string ConfigName)
+bool checkSoftwareVersions(TFile *file, TFile *prevFile, const std::string& ConfigName)
 {
   bool weirdFile = false;
 
@@ -269,6 +271,7 @@ void ParseArg(int argc, char *argv[]){
 int main(int argc, char *argv[])
 {
   SetMaCh3LoggerFormat();
+  MaCh3Utils::MaCh3Welcome();
   ParseArg(argc, argv);
   CombineChain();
   return 0;

diff --git a/Diagnostics/GetPenaltyTerm.cpp b/Diagnostics/GetPenaltyTerm.cpp
@@ -1,10 +1,3 @@
-// C++ includes
-#include <iostream>
-#include <sstream>
-#include <iomanip>
-#include <vector>
-#include <algorithm>
-
 // ROOT includes
 #include "TFile.h"
 #include "TBranch.h"
@@ -29,11 +22,14 @@
 
 #include "manager/manager.h"
 
-//KS: Simple script to get whatever penalty term you like, since flux and xsec are on systematic we cannot just directly take it from the chain
-//g++ `root-config --cflags` -std=c++11 -g -o GetPenaltyTerm GetPenaltyTerm.cpp -I`root-config --incdir` `root-config --glibs --libs`
 
-void GetPenaltyTerm(std::string inputFile, std::string configFile);
-void ReadXSecFile(std::string inputFile);
+/// @file GetPenaltyTerm.cpp
+/// @brief KS: This file contains the implementation of the function to extract specific penalty terms from systematic chains.
+///
+/// This script is designed to retrieve penalty terms from various sources, such as flux and cross-section systematic chains.
+/// Since flux and cross-section uncertainties are handled systematically, the penalty term cannot be taken directly from the chain.
+///
+/// @todo KS: This should really be moved to MCMC Processor
 
 std::vector <int> nominal;
 std::vector <bool> isFlat;
@@ -43,24 +39,92 @@ int size;
 
 double** invCovMatrix;
 
-int main(int argc, char *argv[]) 
+void ReadXSecFile(const std::string& inputFile)
 {
-  SetMaCh3LoggerFormat();
-  if (argc != 3 )
+  // Now read the MCMC file
+  TFile *TempFile = new TFile(inputFile.c_str(), "open");
+
+  // Get the matrix
+  TMatrixDSym *XSecMatrix = (TMatrixDSym*)(TempFile->Get("CovarianceFolder/xsec_cov"));
+
+  // Get the settings for the MCMC
+  TMacro *Config = (TMacro*)(TempFile->Get("MaCh3_Config"));
+  if (Config == nullptr) {
+    MACH3LOG_ERROR("Didn't find MaCh3_Config tree in MCMC file! {}", inputFile);
+    TempFile->ls();
+    throw MaCh3Exception(__FILE__ , __LINE__ );
+  }
+
+  YAML::Node Settings = TMacroToYAML(*Config);
+
+  //CW: Get the xsec Covariance matrix
+  std::vector<std::string> XsecCovPos = GetFromManager<std::vector<std::string>>(Settings["General"]["Systematics"]["XsecCovFile"], {"none"});
+  if(XsecCovPos.back() == "none")
   {
-    MACH3LOG_WARN("Something went wrong ");
-    MACH3LOG_WARN("./GetPenaltyTerm root_file_to_analyse.root ");
+    MACH3LOG_WARN("Couldn't find XsecCov branch in output");
+    MaCh3Utils::PrintConfig(Settings);
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
-  std::string filename = argv[1];
-  std::string config = argv[2];
-  GetPenaltyTerm(filename, config);
 
-  return 0;
+  //KS:Most inputs are in ${MACH3}/inputs/blarb.root
+  if (std::getenv("MACH3") != nullptr) {
+    MACH3LOG_INFO("Found MACH3 environment variable: {}", std::getenv("MACH3"));
+    for(unsigned int i = 0; i < XsecCovPos.size(); i++)
+      XsecCovPos[i].insert(0, std::string(std::getenv("MACH3"))+"/");
+  }
+
+  YAML::Node XSecFile;
+  XSecFile["Systematics"] = YAML::Node(YAML::NodeType::Sequence);
+  for(unsigned int i = 0; i < XsecCovPos.size(); i++)
+  {
+    YAML::Node YAMLDocTemp = YAML::LoadFile(XsecCovPos[i]);
+    for (const auto& item : YAMLDocTemp["Systematics"]) {
+      XSecFile["Systematics"].push_back(item);
+    }
+  }
+
+  size = XSecMatrix->GetNrows();
+
+  auto systematics = XSecFile["Systematics"];
+  for (auto it = systematics.begin(); it != systematics.end(); ++it)
+  {
+    auto const &param = *it;
+
+    ParamNames.push_back(param["Systematic"]["Names"]["FancyName"].as<std::string>());
+    nominal.push_back( param["Systematic"]["ParameterValues"]["PreFitValue"].as<double>() );
+
+    bool flat = false;
+    if (param["Systematic"]["FlatPrior"]) { flat = param["Systematic"]["FlatPrior"].as<bool>(); }
+    isFlat.push_back( flat );
+  }
+
+  XSecMatrix->Invert();
+  //KS: Let's use double as it is faster than TMatrix
+  invCovMatrix = new double*[size];
+  for (int i = 0; i < size; i++)
+  {
+    invCovMatrix[i] = new double[size];
+    for (int j = 0; j < size; ++j)
+    {
+      invCovMatrix[i][j] = -999;
+    }
+  }
+  #ifdef MULTITHREAD
+  #pragma omp parallel for
+  #endif
+  for (int i = 0; i < size; i++)
+  {
+    for (int j = 0; j < size; ++j)
+    {
+      invCovMatrix[i][j] = (*XSecMatrix)(i,j);
+    }
+  }
+
+  TempFile->Close();
+  delete TempFile;
 }
 
-//KS: This should really be moved to MCMC Processor
-void GetPenaltyTerm(std::string inputFile, std::string configFile)
+void GetPenaltyTerm(const std::string& inputFile, const std::string& configFile)
 {
   TCanvas* canvas = new TCanvas("canvas", "canvas", 0, 0, 1024, 1024);
   canvas->SetGrid();
@@ -76,11 +140,6 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile)
   gStyle->SetPalette(51);
 
   ReadXSecFile(inputFile);
-
-  // KS: This can reduce time necessary for caching even by half
-  #ifdef MULTITHREAD
-  //ROOT::EnableImplicitMT();
-  #endif
 
   // Open the Chain
   TChain* Chain = new TChain("posteriors","");
@@ -194,7 +253,7 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile)
     hLogL[i]->SetLineColor(kBlue);
   }
   double* logL = new double[NSets]();
-  for(int n = 0; n < AllEvents; n++)
+  for(int n = 0; n < AllEvents; ++n)
   {
     if(n%10000 == 0) MaCh3Utils::PrintProgressBar(n, AllEvents);
 
@@ -282,7 +341,7 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile)
   canvas->Print(Form("%s_PenaltyTerm.pdf[",inputFile.c_str()), "pdf");
   for(int i = 0; i < NSets; i++)
   {
-    double Maximum = hLogL[i]->GetMaximum();
+    const double Maximum = hLogL[i]->GetMaximum();
     hLogL[i]->GetYaxis()->SetRangeUser(0., Maximum*1.2);
     hLogL[i]->SetTitle(FancyTittle[i].c_str());
     hLogL[i]->GetXaxis()->SetTitle("Step");
@@ -308,89 +367,23 @@ void GetPenaltyTerm(std::string inputFile, std::string configFile)
     delete[] invCovMatrix[i];
   }
   delete[] invCovMatrix;
+  OutputFile->Close();
+  delete OutputFile;
 }
 
-void ReadXSecFile(std::string inputFile)
+int main(int argc, char *argv[])
 {
-  // Now read the MCMC file
-  TFile *TempFile = new TFile(inputFile.c_str(), "open");
-
-  // Get the matrix
-  TMatrixDSym *XSecMatrix = (TMatrixDSym*)(TempFile->Get("CovarianceFolder/xsec_cov"));
-
-  // Get the settings for the MCMC
-  TMacro *Config = (TMacro*)(TempFile->Get("MaCh3_Config"));
-  if (Config == nullptr) {
-    MACH3LOG_ERROR("Didn't find MaCh3_Config tree in MCMC file! {}", inputFile);
-    TempFile->ls();
-    throw MaCh3Exception(__FILE__ , __LINE__ );
-  }
-
-  YAML::Node Settings = TMacroToYAML(*Config);
-
-  //CW: Get the xsec Covariance matrix
-  std::vector<std::string> XsecCovPos = GetFromManager<std::vector<std::string>>(Settings["General"]["Systematics"]["XsecCovFile"], {"none"});
-  if(XsecCovPos.back() == "none")
+  SetMaCh3LoggerFormat();
+  MaCh3Utils::MaCh3Welcome();
+  if (argc != 3 )
   {
-    MACH3LOG_WARN("Couldn't find XsecCov branch in output");
-    MaCh3Utils::PrintConfig(Settings);
+    MACH3LOG_WARN("Something went wrong ");
+    MACH3LOG_WARN("./GetPenaltyTerm root_file_to_analyse.root ");
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
+  std::string filename = argv[1];
+  std::string config = argv[2];
+  GetPenaltyTerm(filename, config);
 
-  //KS:Most inputs are in ${MACH3}/inputs/blarb.root
-  if (std::getenv("MACH3") != nullptr) {
-    MACH3LOG_INFO("Found MACH3 environment variable: {}", std::getenv("MACH3"));
-    for(unsigned int i = 0; i < XsecCovPos.size(); i++)
-      XsecCovPos[i].insert(0, std::string(std::getenv("MACH3"))+"/");
-  }
-
-  YAML::Node XSecFile;
-  XSecFile["Systematics"] = YAML::Node(YAML::NodeType::Sequence);
-  for(unsigned int i = 0; i < XsecCovPos.size(); i++)
-  {
-    YAML::Node YAMLDocTemp = YAML::LoadFile(XsecCovPos[i]);
-    for (const auto& item : YAMLDocTemp["Systematics"]) {
-      XSecFile["Systematics"].push_back(item);
-    }
-  }
-
-  size = XSecMatrix->GetNrows();
-
-  auto systematics = XSecFile["Systematics"];
-  for (auto it = systematics.begin(); it != systematics.end(); ++it)
-  {
-    auto const &param = *it;
-
-    ParamNames.push_back(param["Systematic"]["Names"]["FancyName"].as<std::string>());
-    nominal.push_back( param["Systematic"]["ParameterValues"]["PreFitValue"].as<double>() );
-
-    bool flat = false;
-    if (param["Systematic"]["FlatPrior"]) { flat = param["Systematic"]["FlatPrior"].as<bool>(); }
-    isFlat.push_back( flat );
-  }
-
-  XSecMatrix->Invert();
-  //KS: Let's use double as it is faster than TMatrix
-  invCovMatrix = new double*[size];
-  for (int i = 0; i < size; i++) 
-  {
-    invCovMatrix[i] = new double[size];
-    for (int j = 0; j < size; ++j) 
-    {
-      invCovMatrix[i][j] = -999;
-    }
-  }
-#ifdef MULTITHREAD
-#pragma omp parallel for
-#endif
-  for (int i = 0; i < size; i++)
-  {
-    for (int j = 0; j < size; ++j) 
-    {
-      invCovMatrix[i][j] = (*XSecMatrix)(i,j);
-    }
-  }  
-
-  TempFile->Close();
-  delete TempFile;
+  return 0;
 }
diff --git a/Diagnostics/RHat.cpp b/Diagnostics/RHat.cpp
@@ -1,12 +1,3 @@
-// C++ includes
-#include <iostream>
-#include <sstream>
-#include <iomanip>
-#include <vector>
-#include <algorithm>
-#include <cmath>
-#include <math.h>
-
 // ROOT includes
 #include "TObjArray.h"
 #include "TChain.h"
@@ -27,11 +18,17 @@
 #include "omp.h"
 #endif
 
+// MaCh3 includes
 #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
+/// @file RHat.cpp
+/// @brief This executable calculates the \f$ \hat{R} \f$ estimator for Markov Chain Monte Carlo (MCMC) convergence.
+///
+/// KS: This exe is meant to calculate the \f$ \hat{R} \f$ estimator. For a well-converged chain, this distribution
+/// should be centered at one. The \f$ \hat{R} \f$ statistic is used to assess the convergence of MCMC simulations
+/// and helps determine whether the chains have reached a stable distribution.
+///
+/// @cite gelman2019.
 
 // *******************
 int Ntoys;

diff --git a/Doc/bibliography.bib b/Doc/bibliography.bib
@@ -160,3 +160,11 @@ @misc{rossetti2024batch
     chapter = "5.4",
     url = "https://rossetti.github.io/RossettiArenaBook/ch5-BatchMeansMethod.html",
 }
+
+@article{gelman2019,
+    author       = "Aki Vehtari and Andrew Gelman and Daniel Simpson and Bob Carpenter and Paul-Christian Bürkner",
+    title        = "{Rank-normalization, folding, and localization: An improved \(\widehat{R}\) for assessing convergence of MCMC}",
+    year         = "2019",
+    doi          = "10.48550/arXiv.1903.08008",
+    url          = "https://doi.org/10.48550/arXiv.1903.08008"
+}