diff --git a/.github/workflows/CIValidations.yml b/.github/workflows/CIValidations.yml
index 6ea11dce..aae7f4a5 100644
--- a/.github/workflows/CIValidations.yml
+++ b/.github/workflows/CIValidations.yml
@@ -28,6 +28,7 @@ jobs:
             test_4: empty
             test_5: empty
             test_6: empty
+            test_7: empty
 
           - name: Covariance Validations
             test_1: ./CIValidations/CovarianceValidations
@@ -36,6 +37,7 @@ jobs:
             test_4: empty
             test_5: empty
             test_6: empty
+            test_7: empty
 
           - name: Fitter Validations
             test_1: ./CIValidations/FitterValidations
@@ -44,10 +46,11 @@ jobs:
             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: ./bin/GetPenaltyTerm MCMC_Test.root bin/TutorialDiagConfig.yaml
+            test_7: ./bin/MatrixPlotter bin/TutorialDiagConfig.yaml MCMC_Test_drawCorr.root
             # TODO
             #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
+            #test_8: bin/GetPostfitParamPlots -o MCMC_Test_drawCorr.root -c Inputs/PlottingConfig.yaml
+            #test_9: bin/PlotLLH -o MCMC_Test.root -c Inputs/PlottingConfig.yaml
 
           - name: NuMCMC Tools Validations
             test_1: ./CIValidations/NuMCMCvalidations.sh
@@ -56,6 +59,7 @@ jobs:
             test_4: empty
             test_5: empty
             test_6: empty
+            test_7: empty
 
     container:
       image: ghcr.io/mach3-software/mach3:alma9latest
@@ -118,3 +122,9 @@ jobs:
             echo "Performing  test 6"
             ${{ matrix.test_6 }}
         fi
+
+        if [[ "${{ matrix.test_7 }}" != "empty" ]]; then
+            echo " "
+            echo "Performing  test 7"
+            ${{ matrix.test_7 }}
+        fi
diff --git a/CMakeLists.txt b/CMakeLists.txt
index bd14b63d..ef0e3956 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -131,6 +131,7 @@ else()
   target_compile_options(MaCh3CompilerOptions INTERFACE
     -O3                                   # Optimize code for maximum speed
     -finline-limit=100000000              # Increase the limit for inlining functions to improve performance
+    #-fstrict-aliasing                    # Assume that pointers of different types do not point to the same memory location
     # KS: After benchmarking below didn't in fact worse performance, leave it for future tests and documentation
     #-funroll-loops                       # Unroll loops where possible for performance
     #--param=max-vartrack-size=100000000  # Set maximum size of variable tracking data to avoid excessive memory usage
diff --git a/covariance/covarianceBase.cpp b/covariance/covarianceBase.cpp
index 633f5d9a..c86f4a15 100644
--- a/covariance/covarianceBase.cpp
+++ b/covariance/covarianceBase.cpp
@@ -117,7 +117,7 @@ void covarianceBase::init(std::string name, std::string file) {
   }
 
   // Should put in a 
-  TMatrixDSym *CovMat = (TMatrixDSym*)(infile->Get(name.c_str()));
+  TMatrixDSym *CovMat = static_cast<TMatrixDSym*>(infile->Get(name.c_str()));
   if (CovMat == nullptr) {
     MACH3LOG_ERROR("Could not find covariance matrix name {} in file {}", name, file);
     MACH3LOG_ERROR("Are you really sure {} exists in the file?", name);
@@ -348,7 +348,7 @@ void covarianceBase::setCovMatrix(TMatrixDSym *cov) {
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   covMatrix = cov;
-  invCovMatrix = (TMatrixDSym*)cov->Clone();
+  invCovMatrix = static_cast<TMatrixDSym*>(cov->Clone());
   invCovMatrix->Invert();
   //KS: ROOT has bad memory management, using standard double means we can decrease most operation by factor 2 simply due to cache hits
   for (int i = 0; i < _fNumPar; i++)
@@ -899,7 +899,7 @@ void covarianceBase::SetBranches(TTree &tree, bool SaveProposal) {
   // When running PCA, also save PCA parameters
   if (pca) {
     for (int i = 0; i < _fNumParPCA; ++i) {
-      tree.Branch(Form("%s_PCA", _fNames[i].c_str()), (double*)&(fParCurr_PCA.GetMatrixArray()[i]), Form("%s_PCA/D", _fNames[i].c_str()));
+      tree.Branch(Form("%s_PCA", _fNames[i].c_str()), static_cast<double*>(&fParCurr_PCA.GetMatrixArray()[i]), Form("%s_PCA/D", _fNames[i].c_str()));
     }
   }
 
@@ -912,7 +912,7 @@ void covarianceBase::SetBranches(TTree &tree, bool SaveProposal) {
     // When running PCA, also save PCA parameters
     if (pca) {
       for (int i = 0; i < _fNumParPCA; ++i) {
-        tree.Branch(Form("%s_PCA_Prop", _fNames[i].c_str()), (double*)&(fParProp_PCA.GetMatrixArray()[i]), Form("%s_PCA_Prop/D", _fNames[i].c_str()));
+        tree.Branch(Form("%s_PCA_Prop", _fNames[i].c_str()), static_cast<double*>(&fParProp_PCA.GetMatrixArray()[i]), Form("%s_PCA_Prop/D", _fNames[i].c_str()));
       }
     }
   }
@@ -1155,7 +1155,7 @@ void covarianceBase::setThrowMatrix(TMatrixDSym *cov){
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
 
-  throwMatrix = (TMatrixDSym*)cov->Clone();
+  throwMatrix = static_cast<TMatrixDSym*>(cov->Clone());
   if(use_adaptive && AdaptiveHandler.AdaptionUpdate()) makeClosestPosDef(throwMatrix);
   else MakePosDef(throwMatrix);
   
@@ -1274,11 +1274,11 @@ void covarianceBase::makeClosestPosDef(TMatrixDSym *cov) {
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   
-  TMatrixD cov_sym_v = (TMatrixD)cov_sym_svd.GetV();
+  TMatrixD cov_sym_v = static_cast<TMatrixD>(cov_sym_svd.GetV());
   TMatrixD cov_sym_vt = cov_sym_v;
   cov_sym_vt.T();
   //SVD returns as vector (grrr) so need to get into matrix form for multiplying!
-  TVectorD cov_sym_sigvect = (TVectorD)cov_sym_svd.GetSig();
+  TVectorD cov_sym_sigvect = static_cast<TVectorD>(cov_sym_svd.GetSig());
   const Int_t nCols = cov_sym_v.GetNcols(); //square so only need rows hence lack of cols
   TMatrixDSym cov_sym_sig(nCols);
   TMatrixDDiag cov_sym_sig_diag(cov_sym_sig);
diff --git a/mcmc/FitterBase.cpp b/mcmc/FitterBase.cpp
index b5f10514..7875505c 100644
--- a/mcmc/FitterBase.cpp
+++ b/mcmc/FitterBase.cpp
@@ -13,21 +13,20 @@ FitterBase::FitterBase(manager * const man) : fitMan(man) {
 
   //Get mach3 modes from manager
   Modes = fitMan->GetMaCh3Modes();
-
-  random = new TRandom3(fitMan->raw()["General"]["Seed"].as<int>());
+  random = std::make_unique<TRandom3>(fitMan->raw()["General"]["Seed"].as<int>());
 
   // Counter of the accepted # of steps
   accCount = 0;
   step = 0;
 
-  clock = new TStopwatch;
-  stepClock = new TStopwatch;
+  clock = std::make_unique<TStopwatch>();
+  stepClock = std::make_unique<TStopwatch>();
   #ifdef DEBUG
   // Fit summary and debug info
   debug = GetFromManager<bool>(fitMan->raw()["General"]["Debug"], false);
   #endif
 
-  std::string outfile = fitMan->raw()["General"]["OutputFile"].as<std::string>();
+  auto outfile = fitMan->raw()["General"]["OutputFile"].as<std::string>();
   // Save output every auto_save steps
   //you don't want this too often https://root.cern/root/html606/TTree_8cxx_source.html#l01229
   auto_save = fitMan->raw()["General"]["MCMC"]["AutoSave"].as<int>();
@@ -63,13 +62,6 @@ FitterBase::FitterBase(manager * const man) : fitMan(man) {
   if (debug) debugFile.open((outfile+".log").c_str());
   #endif
 
-  // Clear the samples and systematics
-  samples.clear();
-  systematics.clear();
-
-  sample_llh = nullptr;
-  syst_llh = nullptr;
-
   TotalNSamples = 0;
   fTestLikelihood = GetFromManager<bool>(fitMan->raw()["General"]["Fitter"]["FitTestLikelihood"], false);
 }
@@ -80,18 +72,8 @@ FitterBase::~FitterBase() {
 // *************************
   SaveOutput();
 
-  if(random != nullptr) delete random;
-  random = nullptr;
-  if(sample_llh != nullptr) delete[] sample_llh;
-  sample_llh = nullptr;
-  if(syst_llh != nullptr) delete[] syst_llh;
-  syst_llh = nullptr;
   if(outputFile != nullptr) delete outputFile;
   outputFile = nullptr;
-  if(clock != nullptr) delete clock;
-  clock = nullptr;
-  if(stepClock != nullptr) delete stepClock;
-  stepClock = nullptr;
   MACH3LOG_DEBUG("Closing MaCh3 Fitter Engine");
 }
 
@@ -99,7 +81,6 @@ FitterBase::~FitterBase() {
 // Prepare the output tree
 void FitterBase::SaveSettings() {
 // *******************
-
   if(SettingsSaved) return;
 
   outputFile->cd();
@@ -188,8 +169,8 @@ void FitterBase::PrepareOutput() {
 
     // Store individual likelihood components
     // Using double means double as large output file!
-    sample_llh = new double[samples.size()];
-    syst_llh = new double[systematics.size()];
+    sample_llh.resize(samples.size());
+    syst_llh.resize(systematics.size());
 
     for (size_t i = 0; i < samples.size(); ++i) {
       std::stringstream oss, oss2;
@@ -231,7 +212,6 @@ void FitterBase::PrepareOutput() {
 // *******************
 void FitterBase::SaveOutput() {
 // *******************
-
   if(FileSaved) return;
   //Stop Clock
   clock->Stop();
@@ -268,18 +248,16 @@ void FitterBase::addSamplePDF(samplePDFBase * const sample) {
 // Add flux systematics, cross-section systematics, ND systematics to the chain
 void FitterBase::addSystObj(covarianceBase * const cov) {
 // *************************
-
   MACH3LOG_INFO("Adding systematic object {}, with {} params", cov->getName(), cov->GetNumParams());
   systematics.push_back(cov);
 
   CovFolder->cd();
-  double *n_vec = new double[cov->GetNumParams()];
+  std::vector<double> n_vec(cov->GetNumParams());
   for (int i = 0; i < cov->GetNumParams(); ++i)
     n_vec[i] = cov->getParInit(i);
 
-  TVectorT<double> t_vec(cov->GetNumParams(), n_vec);
+  TVectorT<double> t_vec(cov->GetNumParams(), n_vec.data());
   t_vec.Write((std::string(cov->getName()) + "_prior").c_str());
-  delete[] n_vec;
 
   cov->getCovMatrix()->Write(cov->getName().c_str());
 
@@ -374,7 +352,7 @@ void FitterBase::StartFromPreviousFit(const std::string& FitName) {
 // Process the MCMC output to get postfit etc
 void FitterBase::ProcessMCMC() {
 // *******************
-  if (fitMan == NULL) return;
+  if (fitMan == nullptr) return;
 
   // Process the MCMC
   if (fitMan->raw()["General"]["ProcessMCMC"].as<bool>()) {
@@ -522,10 +500,10 @@ void FitterBase::RunLLHScan() {
   TDirectory *Sample_LLH = outputFile->mkdir("Sample_LLH");
   TDirectory *Total_LLH = outputFile->mkdir("Total_LLH");
 
-  TDirectory **SampleSplit_LLH = nullptr;
+  std::vector<TDirectory *>SampleSplit_LLH;
   if(PlotAllNDsamplesLLH)
   {
-    SampleSplit_LLH = new TDirectory*[TotalNSamples];
+    SampleSplit_LLH.resize(TotalNSamples);
     int SampleIterator = 0;
     for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
     {
@@ -590,12 +568,11 @@ void FitterBase::RunLLHScan() {
       if (IsPCA) {
         lower = prior - nSigma*sqrt(((*it)->getEigenValues())(i));
         upper = prior + nSigma*sqrt(((*it)->getEigenValues())(i));
-
-        std::cout << "eval " << i << " = " << (*it)->getEigenValues()(i) << std::endl;
-        std::cout << "prior " << i << " = " << prior << std::endl;
-        std::cout << "lower " << i << " = " << lower << std::endl;
-        std::cout << "upper " << i << " = " << upper << std::endl;
-        std::cout << "nSigma " << nSigma << std::endl;
+        MACH3LOG_INFO("eval {} = {:.2f}", i, (*it)->getEigenValues()(i));
+        MACH3LOG_INFO("prior {} = {:.2f}", i, prior);
+        MACH3LOG_INFO("lower {} = {:.2f}", i, lower);
+        MACH3LOG_INFO("upper {} = {:.2f}", i, upper);
+        MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
       }
 
       // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables
@@ -608,15 +585,14 @@ void FitterBase::RunLLHScan() {
       MACH3LOG_INFO("Scanning {} with {} steps, from {:.2f} - {:.2f}, prior = {:.2f}", name, n_points, lower, upper, prior);
 
       // Make the TH1D
-      TH1D *hScan = new TH1D((name+"_full").c_str(), (name+"_full").c_str(), n_points, lower, upper);
+      auto hScan = std::make_unique<TH1D>((name + "_full").c_str(), (name + "_full").c_str(), n_points, lower, upper);
       hScan->SetTitle(std::string(std::string("2LLH_full, ") + name + ";" + name + "; -2(ln L_{sample} + ln L_{xsec+flux} + ln L_{det})").c_str());
 
-      TH1D *hScanSam = new TH1D((name+"_sam").c_str(), (name+"_sam").c_str(), n_points, lower, upper);
+      auto hScanSam = std::make_unique<TH1D>((name + "_sam").c_str(), (name + "_sam").c_str(), n_points, lower, upper);
       hScanSam->SetTitle(std::string(std::string("2LLH_sam, ") + name + ";" + name + "; -2(ln L_{sample})").c_str());
 
-
-      TH1D **hScanSample = new TH1D*[samples.size()];
-      double *nSamLLH = new double[samples.size()];
+      std::vector<TH1D*> hScanSample(samples.size());
+      std::vector<double> nSamLLH(samples.size());
       for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
       {
         std::string NameTemp = samples[ivs]->GetName();
@@ -625,8 +601,8 @@ void FitterBase::RunLLHScan() {
         nSamLLH[ivs] = 0.;
       }
 
-      TH1D **hScanCov = new TH1D*[systematics.size()];
-      double *nCovLLH = new double[systematics.size()];
+      std::vector<TH1D*> hScanCov(systematics.size());
+      std::vector<double> nCovLLH(systematics.size());
       for(unsigned int ivc = 0; ivc < systematics.size(); ivc++ )
       {
         std::string NameTemp = systematics[ivc]->getName();
@@ -637,18 +613,18 @@ void FitterBase::RunLLHScan() {
         nCovLLH[ivc] = 0.;
       }
 
-      TH1D **hScanSamSplit = nullptr;
-      double *sampleSplitllh = nullptr;
+      std::vector<TH1D *> hScanSamSplit;
+      std::vector<double> sampleSplitllh;
       if(PlotAllNDsamplesLLH)
       {
         int SampleIterator = 0;
         for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
         {
-          hScanSamSplit = new TH1D*[TotalNSamples];
-          sampleSplitllh = new double[TotalNSamples];
+          hScanSamSplit.resize(TotalNSamples);
+          sampleSplitllh.resize(TotalNSamples);
           for(_int_ is = 0; is < samples[ivs]->GetNsamples(); is++ )
           {
-            hScanSamSplit[SampleIterator] = new TH1D( (name+samples[ivs]->GetSampleName(is)).c_str(), (name+samples[ivs]->GetSampleName(is)).c_str(), n_points, lower, upper );
+            hScanSamSplit[SampleIterator] = new TH1D((name+samples[ivs]->GetSampleName(is)).c_str(), (name+samples[ivs]->GetSampleName(is)).c_str(), n_points, lower, upper);
             hScanSamSplit[SampleIterator]->SetTitle(std::string(std::string("2LLH_sam, ") + name + ";" + name + "; -2(ln L_{sample})").c_str());
             SampleIterator++;
           }
@@ -745,26 +721,11 @@ void FitterBase::RunLLHScan() {
         hScanSample[ivs]->Write();
         delete hScanSample[ivs];
       }
-
       Sample_LLH->cd();
       hScanSam->Write();
       Total_LLH->cd();
       hScan->Write();
 
-      delete[] hScanCov;
-      delete[] nCovLLH;
-      delete[] hScanSample;
-      delete[] nSamLLH;
-      delete hScanSam;
-      delete hScan;
-
-      hScanCov = nullptr;
-      nCovLLH = nullptr;
-      hScanSample = nullptr;
-      nSamLLH = nullptr;
-      hScanSam = nullptr;
-      hScan = nullptr;
-
       if(PlotAllNDsamplesLLH)
       {
         int SampleIterator = 0;
@@ -778,7 +739,6 @@ void FitterBase::RunLLHScan() {
             SampleIterator++;
           }
         }
-        delete[] hScanSamSplit;
       }
 
       // Reset the parameters to their prior central values
@@ -827,8 +787,7 @@ void FitterBase::RunLLHScan() {
 //LLH scan is good first estimate of step scale
 void FitterBase::GetStepScaleBasedOnLLHScan() {
 // *************************
-  TDirectory *Sample_LLH = (TDirectory*) outputFile->Get("Sample_LLH");
-
+  TDirectory *Sample_LLH = outputFile->Get<TDirectory>("Sample_LLH");
   MACH3LOG_INFO("Starting Get Step Scale Based On LLHScan");
 
   if(Sample_LLH->IsZombie())
@@ -850,8 +809,7 @@ void FitterBase::GetStepScaleBasedOnLLHScan() {
       if (isxsec) name = (*it)->GetParFancyName(i);
 
       StepScale[i] = (*it)->GetIndivStepScale(i);
-      TH1D* LLHScan = (TH1D*) Sample_LLH->Get((name+"_sam").c_str());
-
+      TH1D* LLHScan = static_cast<TH1D*>(Sample_LLH->Get((name + "_sam").c_str()));
       if(LLHScan == nullptr)
       {
         MACH3LOG_WARN("Couldn't find LLH scan, for {}, skipping", name);
@@ -936,14 +894,13 @@ void FitterBase::Run2DLLHScan() {
       double upper_x = prior_x + nSigma*(*it)->getDiagonalError(i);
       // If PCA, transform these parameter values to the PCA basis
       if (IsPCA) {
-        lower_x = prior_x - nSigma*sqrt(((*it)->getEigenValues())(i));
-        upper_x = prior_x + nSigma*sqrt(((*it)->getEigenValues())(i));
-
-        std::cout << "eval " << i << " = " << (*it)->getEigenValues()(i) << std::endl;
-        std::cout << "prior " << i << " = " << prior_x << std::endl;
-        std::cout << "lower " << i << " = " << lower_x << std::endl;
-        std::cout << "upper " << i << " = " << upper_x << std::endl;
-        std::cout << "nSigma " << nSigma << std::endl;
+        lower_x = prior_x - nSigma*std::sqrt(((*it)->getEigenValues())(i));
+        upper_x = prior_x + nSigma*std::sqrt(((*it)->getEigenValues())(i));
+        MACH3LOG_INFO("eval {} = {:.2f}", i, (*it)->getEigenValues()(i));
+        MACH3LOG_INFO("prior {} = {:.2f}", i, prior_x);
+        MACH3LOG_INFO("lower {} = {:.2f}", i, lower_x);
+        MACH3LOG_INFO("upper {} = {:.2f}", i, upper_x);
+        MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
       }
 
       // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables
@@ -994,12 +951,11 @@ void FitterBase::Run2DLLHScan() {
         if (IsPCA) {
           lower_y = prior_y - nSigma*sqrt(((*it)->getEigenValues())(j));
           upper_y = prior_y + nSigma*sqrt(((*it)->getEigenValues())(j));
-
-          std::cout << "eval " << j << " = " << (*it)->getEigenValues()(j) << std::endl;
-          std::cout << "prior " << j << " = " << prior_y << std::endl;
-          std::cout << "lower " << j << " = " << lower_y << std::endl;
-          std::cout << "upper " << j << " = " << upper_y << std::endl;
-          std::cout << "nSigma " << nSigma << std::endl;
+          MACH3LOG_INFO("eval {} = {:.2f}", i, (*it)->getEigenValues()(j));
+          MACH3LOG_INFO("prior {} = {:.2f}", i, prior_y);
+          MACH3LOG_INFO("lower {} = {:.2f}", i, lower_y);
+          MACH3LOG_INFO("upper {} = {:.2f}", i, upper_y);
+          MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
         }
 
         // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables
@@ -1012,7 +968,8 @@ void FitterBase::Run2DLLHScan() {
         MACH3LOG_INFO("Scanning X {} with {} steps, from {} - {}, prior = {}", name_x, n_points, lower_x, upper_x, prior_x);
         MACH3LOG_INFO("Scanning Y {} with {} steps, from {} - {}, prior = {}", name_y, n_points, lower_y, upper_y, prior_y);
 
-        TH2D *hScanSam = new TH2D((name_x + "_" + name_y + "_sam").c_str(), (name_x + "_" + name_y + "_sam").c_str(), n_points, lower_x, upper_x, n_points, lower_y, upper_y);
+        auto hScanSam = std::make_unique<TH2D>((name_x + "_" + name_y + "_sam").c_str(), (name_x + "_" + name_y + "_sam").c_str(),
+                                                n_points, lower_x, upper_x, n_points, lower_y, upper_y);
         hScanSam->GetXaxis()->SetTitle(name_x.c_str());
         hScanSam->GetYaxis()->SetTitle(name_y.c_str());
         hScanSam->GetZaxis()->SetTitle("2LLH_sam");
@@ -1052,10 +1009,6 @@ void FitterBase::Run2DLLHScan() {
 
         Sample_2DLLH->cd();
         hScanSam->Write();
-
-        delete hScanSam;
-        hScanSam = nullptr;
-
         // Reset the parameters to their prior central values
         if (IsPCA) {
           (*it)->setParProp_PCA(i, prior_x);
@@ -1103,14 +1056,13 @@ void FitterBase::RunSigmaVar() {
   }
 
   bool isxsec = false;
-
   for (std::vector<covarianceBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
   {
     TMatrixDSym *Cov = (*it)->getCovMatrix();
 
     if((*it)->IsPCA())
     {
-      MACH3LOG_ERROR("Usin PCAed matrix not implemented within simga var code, I am sorry :(");
+      MACH3LOG_ERROR("Using PCAed matrix not implemented within sigma var code, I am sorry :(");
       throw MaCh3Exception(__FILE__ , __LINE__ );
     }
 
@@ -1218,7 +1170,7 @@ void FitterBase::RunSigmaVar() {
             std::string parVarTitle = std::string(name) + "_" + ss.str();
 
             // This is a TH2D
-            TH2Poly* currSamp = (TH2Poly*)(samples[ivs]->getPDF(k))->Clone();
+            TH2Poly* currSamp = static_cast<TH2Poly*>(samples[ivs]->getPDF(k)->Clone());
             currSamp->SetDirectory(0);
             // Set a descriptiv-ish title
             std::string title_long = std::string(currSamp->GetName())+"_"+parVarTitle;
@@ -1235,11 +1187,11 @@ void FitterBase::RunSigmaVar() {
               sigmaArray_mode_x[j][SampleIterator] = new TH1D*[nRelevantModes]();
               sigmaArray_mode_y[j][SampleIterator] = new TH1D*[nRelevantModes]();
               // Now get the TH2D mode variations
-              TH2Poly** currSampMode; currSampMode = new TH2Poly*[nRelevantModes]();
+              TH2Poly** currSampMode = new TH2Poly*[nRelevantModes]();
               std::string mode_title_long;
               for(int ir = 0; ir < nRelevantModes; ir++)
               {
-                currSampMode[ir] = (TH2Poly*)(samples[ivs]->getPDFMode(k, RelevantModes[ir]))->Clone();
+                currSampMode[ir] = static_cast<TH2Poly*>(samples[ivs]->getPDFMode(k, RelevantModes[ir])->Clone());
                 currSampMode[ir]->SetDirectory(0);
 
                 mode_title_long = title_long + "_" + Modes->GetMaCh3ModeName(RelevantModes[ir]);
@@ -1256,18 +1208,17 @@ void FitterBase::RunSigmaVar() {
               delete[] currSampMode;
             }
 
-            //KS: This will give different results depending if data or asimov, both have their uses
+            //KS: This will give different results depending if data or Asimov, both have their uses
             if (PlotLLHperBin)
             {
-              TH2Poly* currLLHSamp = (TH2Poly*)(samples[ivs]->getPDF(k))->Clone();
+              TH2Poly* currLLHSamp = static_cast<TH2Poly*>(samples[ivs]->getPDF(k)->Clone());
               currLLHSamp->SetDirectory(0);
               currLLHSamp->Reset("");
               currLLHSamp->Fill(0.0, 0.0, 0.0);
 
-              TH2Poly* MCpdf = (TH2Poly*)(samples[ivs]->getPDF(k));
-              TH2Poly* Datapdf = (TH2Poly*)(samples[ivs]->getData(k));
-              TH2Poly* W2pdf = (TH2Poly*)(samples[ivs]->getW2(k));
-
+              TH2Poly* MCpdf = static_cast<TH2Poly*>(samples[ivs]->getPDF(k));
+              TH2Poly* Datapdf = static_cast<TH2Poly*>(samples[ivs]->getData(k));
+              TH2Poly* W2pdf = static_cast<TH2Poly*>(samples[ivs]->getW2(k));
               for(int bin = 1; bin < currLLHSamp->GetNumberOfBins()+1; bin++)
               {
                 const double mc = MCpdf->GetBinContent(bin);
@@ -1289,10 +1240,10 @@ void FitterBase::RunSigmaVar() {
             sigmaArray_y[j][SampleIterator]->SetDirectory(0);
             sigmaArray_y[j][SampleIterator]->GetXaxis()->SetTitle(currSamp->GetYaxis()->GetTitle());
 
-            sigmaArray_x_norm[j][SampleIterator] = (TH1D*)sigmaArray_x[j][SampleIterator]->Clone();
+            sigmaArray_x_norm[j][SampleIterator] = static_cast<TH1D*>(sigmaArray_x[j][SampleIterator]->Clone());
             sigmaArray_x_norm[j][SampleIterator]->SetDirectory(0);
             sigmaArray_x_norm[j][SampleIterator]->Scale(1., "width");
-            sigmaArray_y_norm[j][SampleIterator] = (TH1D*)sigmaArray_y[j][SampleIterator]->Clone();
+            sigmaArray_y_norm[j][SampleIterator] = static_cast<TH1D*>(sigmaArray_y[j][SampleIterator]->Clone());
             sigmaArray_y_norm[j][SampleIterator]->SetDirectory(0);
             sigmaArray_y_norm[j][SampleIterator]->Scale(1., "width");
 
@@ -1419,35 +1370,3 @@ void FitterBase::RunSigmaVar() {
     } // end looping over xsec parameters (i)
   } // end looping over covarianceBase objects
 }
-
-// *************************
-TGraphAsymmErrors* FitterBase::MakeAsymGraph(TH1D* sigmaArrayLeft, TH1D* sigmaArrayCentr, TH1D* sigmaArrayRight, const std::string& title) {
-// *************************
-
-  TGraphAsymmErrors *var = new TGraphAsymmErrors(sigmaArrayCentr);
-  var->SetNameTitle((title).c_str(), (title).c_str());
-
-  // Need to draw TGraphs to set axes labels
-  var->Draw("AP");
-  var->GetXaxis()->SetTitle(sigmaArrayCentr->GetXaxis()->GetTitle());
-  var->GetYaxis()->SetTitle("Number of events/bin");
-
-  for (int m = 0; m < var->GetN(); ++m)
-  {
-    double xlow = sigmaArrayLeft->GetBinContent(m+1);
-    double xhigh = sigmaArrayRight->GetBinContent(m+1);
-    double xtemp;
-
-    // Figure out which variation is larger so we set the error correctly
-    if (xlow > xhigh)
-    {
-      xtemp = xlow;
-      xlow = xhigh;
-      xhigh = xtemp;
-    }
-
-    var->SetPointEYhigh(m, xhigh - var->GetY()[m]);
-    var->SetPointEYlow(m, var->GetY()[m] - xlow);
-  }
-  return var;
-}
diff --git a/mcmc/FitterBase.h b/mcmc/FitterBase.h
index 71d33aec..ccca5a41 100644
--- a/mcmc/FitterBase.h
+++ b/mcmc/FitterBase.h
@@ -77,14 +77,6 @@ class FitterBase {
   /// @brief Save the settings that the MCMC was run with.
   void SaveSettings();
 
-  /// @brief Used by sigma variation, check how 1 sigma changes spectra
-  /// @param sigmaArrayLeft sigma var hist at -1 or -3 sigma shift
-  /// @param sigmaArrayCentr sigma var hist at prior values
-  /// @param sigmaArrayRight sigma var hist at +1 or +3 sigma shift
-  /// @param title A tittle for returned object
-  /// @return A `TGraphAsymmErrors` object that visualizes the sigma variation of spectra, showing confidence intervals between different sigma shifts.
-  inline TGraphAsymmErrors* MakeAsymGraph(TH1D* sigmaArrayLeft, TH1D* sigmaArrayCentr, TH1D* sigmaArrayRight, const std::string& title);
-
   /// The manager
   manager *fitMan;
 
@@ -103,9 +95,9 @@ class FitterBase {
   int accCount;
 
   /// store the llh breakdowns
-  double *sample_llh;
+  std::vector<double> sample_llh;
   /// systematic llh breakdowns
-  double *syst_llh;
+  std::vector<double> syst_llh;
 
   /// Sample holder
   std::vector<samplePDFBase*> samples;
@@ -116,14 +108,14 @@ class FitterBase {
   std::vector<covarianceBase*> systematics;
 
   /// tells global time how long fit took
-  TStopwatch* clock;
+  std::unique_ptr<TStopwatch> clock;
   /// tells how long single step/fit iteration took
-  TStopwatch* stepClock;
+  std::unique_ptr<TStopwatch> stepClock;
   /// Time of single step
   double stepTime;
 
   /// Random number
-  TRandom3* random;
+  std::unique_ptr<TRandom3> random;
 
   /// Output
   TFile *outputFile;
diff --git a/mcmc/LikelihoodFit.cpp b/mcmc/LikelihoodFit.cpp
index cdc8696b..afe501e2 100644
--- a/mcmc/LikelihoodFit.cpp
+++ b/mcmc/LikelihoodFit.cpp
@@ -54,10 +54,10 @@ double LikelihoodFit::CalcChi2(const double* x) {
     if(!(*it)->IsPCA())
     {
       std::vector<double> pars;
-      const int Size = (*it)->GetNumParams();
+      const int NumPar = (*it)->GetNumParams();
       //KS: Avoid push back as they are slow
-      pars.resize(Size);
-      for(int i = 0; i < Size; ++i, ++ParCounter)
+      pars.resize(NumPar);
+      for(int i = 0; i < NumPar; ++i, ++ParCounter)
       {
         double ParVal = x[ParCounter];
         //KS: Basically apply mirroring for parameters out of bounds
@@ -79,10 +79,10 @@ double LikelihoodFit::CalcChi2(const double* x) {
     else
     {
       std::vector<double> pars;
-      const int Size = (*it)->getNpars();
+      const int NumPar = (*it)->getNpars();
       //KS: Avoid push back as they are slow
-      pars.resize(Size);
-      for(int i = 0; i < Size; ++i, ++ParCounter)
+      pars.resize(NumPar);
+      for(int i = 0; i < NumPar; ++i, ++ParCounter)
       {
         double ParVal = x[ParCounter];
         //KS: Basically apply mirroring for parameters out of bounds
diff --git a/mcmc/MCMCProcessor.cpp b/mcmc/MCMCProcessor.cpp
index 2ea3955b..badf70c4 100644
--- a/mcmc/MCMCProcessor.cpp
+++ b/mcmc/MCMCProcessor.cpp
@@ -197,8 +197,8 @@ void MCMCProcessor::GetCovariance(TMatrixDSym *&Cov, TMatrixDSym *&Corr) {
 // ***************
   if (CacheMCMC) MakeCovariance_MP();
   else MakeCovariance();
-  Cov = (TMatrixDSym*)Covariance->Clone();
-  Corr = (TMatrixDSym*)Correlation->Clone();
+  Cov = static_cast<TMatrixDSym*>(Covariance->Clone());
+  Corr = static_cast<TMatrixDSym*>(Correlation->Clone());
 }
 
 // ***************
@@ -720,18 +720,16 @@ void MCMCProcessor::MakeCredibleIntervals(const std::vector<double>& CredibleInt
   }
 
   const int nCredible = CredibleIntervals.size();
-  TH1D** hpost_copy = new TH1D*[nDraw];
-  TH1D*** hpost_cl = new TH1D**[nDraw];
-
+  std::vector<TH1D*> hpost_copy(nDraw);
+  std::vector<std::vector<TH1D*>> hpost_cl(nDraw);
   //KS: Copy all histograms to be thread safe
   for (int i = 0; i < nDraw; ++i)
   {
-    hpost_copy[i] = (TH1D*) hpost[i]->Clone(Form("hpost_copy_%i", i));
-    hpost_cl[i] = new TH1D*[nCredible];
-
+    hpost_copy[i] = static_cast<TH1D*>(hpost[i]->Clone(Form("hpost_copy_%i", i)));
+    hpost_cl[i].resize(nCredible);
     for (int j = 0; j < nCredible; ++j)
     {
-      hpost_cl[i][j] = (TH1D*) hpost[i]->Clone( Form("hpost_copy_%i_CL_%f", i, CredibleIntervals[j]));
+      hpost_cl[i][j] = static_cast<TH1D*>(hpost[i]->Clone(Form("hpost_copy_%i_CL_%f", i, CredibleIntervals[j])));
 
       //KS: Reset to get rid to TF1 otherwise we run into segfault :(
       hpost_cl[i][j]->Reset("");
@@ -826,17 +824,12 @@ void MCMCProcessor::MakeCredibleIntervals(const std::vector<double>& CredibleInt
   }
 
   //KS: Remove histograms
-  for (int i = 0; i < nDraw; ++i)
-  {
+  for (int i = 0; i < nDraw; ++i) {
     delete hpost_copy[i];
-    for (int j = 0; j < nCredible; ++j)
-    {
+    for (int j = 0; j < nCredible; ++j) {
       delete hpost_cl[i][j];
     }
-    delete[] hpost_cl[i];
   }
-  delete[] hpost_copy;
-  delete[] hpost_cl;
 
   CredibleDir->Close();
   delete CredibleDir;
@@ -1054,8 +1047,7 @@ void MCMCProcessor::MakeCovariance() {
       const double min_j = Chain->GetMinimum(BranchNames[j]);
 
       // TH2F to hold the Correlation 
-      TH2D *hpost_2D = new TH2D(DrawMe, DrawMe, nBins, min_i, max_i, nBins, min_j, max_j);
-
+      std::unique_ptr<TH2D> hpost_2D = std::make_unique<TH2D>(DrawMe, DrawMe, nBins, min_i, max_i, nBins, min_j, max_j);
       hpost_2D->SetMinimum(0);
       hpost_2D->GetXaxis()->SetTitle(Title_i);
       hpost_2D->GetYaxis()->SetTitle(Title_j);
@@ -1090,8 +1082,6 @@ void MCMCProcessor::MakeCovariance() {
       // Write it to root file
       //OutputFile->cd();
       //if( std::fabs((*Correlation)(i,j)) > Post2DPlotThreshold ) hpost_2D->Write();
-
-      delete hpost_2D;
     } // End j loop
   } // End i loop
   OutputFile->cd();
@@ -1213,7 +1203,6 @@ void MCMCProcessor::CacheSteps() {
 // Make the post-fit covariance matrix in all dimensions
 void MCMCProcessor::MakeCovariance_MP(bool Mute) {
 // *********************
-    
   if (OutputFile == nullptr) MakeOutputFile();
     
   if(!CacheMCMC) CacheSteps();
@@ -1340,7 +1329,7 @@ void MCMCProcessor::MakeSubOptimality(const int NIntervals) {
   TStopwatch clock;
   clock.Start();
 
-  TH1D* SubOptimality = new TH1D("Suboptimality", "Suboptimality", NIntervals, MinStep, MaxStep);
+  std::unique_ptr<TH1D> SubOptimality = std::make_unique<TH1D>("Suboptimality", "Suboptimality", NIntervals, MinStep, MaxStep);
   SubOptimality->GetXaxis()->SetTitle("Step");
   SubOptimality->GetYaxis()->SetTitle("Suboptimality");
   SubOptimality->SetLineWidth(2);
@@ -1386,8 +1375,6 @@ void MCMCProcessor::MakeSubOptimality(const int NIntervals) {
   // Write it to root file
   OutputFile->cd();
   Posterior->Write();
-
-  delete SubOptimality;
 }
 
 // *********************
@@ -1398,13 +1385,13 @@ void MCMCProcessor::DrawCovariance() {
   Posterior->SetRightMargin(0.15);
 
   // The Covariance matrix from the fit
-  TH2D* hCov = new TH2D("hCov", "hCov", nDraw, 0, nDraw, nDraw, 0, nDraw);
+  std::unique_ptr<TH2D> hCov = std::make_unique<TH2D>("hCov", "hCov", nDraw, 0, nDraw, nDraw, 0, nDraw);
   hCov->GetZaxis()->SetTitle("Covariance");
   // The Covariance matrix square root, with correct sign
-  TH2D* hCovSq = new TH2D("hCovSq", "hCovSq", nDraw, 0, nDraw, nDraw, 0, nDraw);
+  std::unique_ptr<TH2D> hCovSq = std::make_unique<TH2D>("hCovSq", "hCovSq", nDraw, 0, nDraw, nDraw, 0, nDraw);
   hCovSq->GetZaxis()->SetTitle("Covariance");
   // The Correlation
-  TH2D* hCorr = new TH2D("hCorr", "hCorr", nDraw, 0, nDraw, nDraw, 0, nDraw);
+  std::unique_ptr<TH2D> hCorr = std::make_unique<TH2D>("hCorr", "hCorr", nDraw, 0, nDraw, nDraw, 0, nDraw);
   hCorr->GetZaxis()->SetTitle("Correlation");
   hCorr->SetMinimum(-1);
   hCorr->SetMaximum(1);
@@ -1480,10 +1467,6 @@ void MCMCProcessor::DrawCovariance() {
   
   //Back to normal
   Posterior->SetRightMargin(RightMargin);
-  delete hCov;
-  delete hCovSq;
-  delete hCorr;
-
   DrawCorrelations1D();
 }
 
@@ -1491,7 +1474,6 @@ void MCMCProcessor::DrawCovariance() {
 //KS: Make the 1D projections of Correlations inspired by Henry's slides (page 28) https://www.t2k.org/asg/oagroup/meeting/2023/2023-07-10-oa-pre-meeting/MaCh3FDUpdate
 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();
@@ -1512,7 +1494,7 @@ void MCMCProcessor::DrawCorrelations1D() {
   std::vector<std::vector<std::string>> NameCorrOfInterest;
   NameCorrOfInterest.resize(nDraw);
 
-  TH1D ***Corr1DHist = new TH1D**[nDraw]();
+  std::vector<std::vector<std::unique_ptr<TH1D>>> Corr1DHist(nDraw);
   //KS: Initialising ROOT objects is never safe in MP loop
   for(int i = 0; i < nDraw; ++i)
   {
@@ -1521,10 +1503,10 @@ void MCMCProcessor::DrawCorrelations1D() {
     double PriorError = 1.0;
     GetNthParameter(i, Prior, PriorError, Title);
 
-    Corr1DHist[i] = new TH1D*[Nhists]();
+    Corr1DHist[i].resize(Nhists);
     for(int j = 0; j < Nhists; ++j)
     {
-      Corr1DHist[i][j] = new TH1D(Form("Corr1DHist_%i_%i", i, j), Form("Corr1DHist_%i_%i", i, j), nDraw, 0, nDraw);
+      Corr1DHist[i][j] = std::make_unique<TH1D>(Form("Corr1DHist_%i_%i", i, j), Form("Corr1DHist_%i_%i", i, j), nDraw, 0, nDraw);
       Corr1DHist[i][j]->SetTitle(Form("%s",Title.Data()));
       Corr1DHist[i][j]->GetYaxis()->SetTitle("Correlation");
       Corr1DHist[i][j]->SetFillColor(CorrColours[j]);
@@ -1583,7 +1565,7 @@ void MCMCProcessor::DrawCorrelations1D() {
     leg->SetTextSize(0.02);
     for(int k = 0; k < Nhists; k++)
     {
-      leg->AddEntry(Corr1DHist[i][k], Form("%.2f > |Corr| >= %.2f", Thresholds[k+1], Thresholds[k]), "f");
+      leg->AddEntry(Corr1DHist[i][k].get(), Form("%.2f > |Corr| >= %.2f", Thresholds[k+1], Thresholds[k]), "f");
     }
     leg->SetLineColor(0);
     leg->SetLineStyle(0);
@@ -1625,16 +1607,6 @@ void MCMCProcessor::DrawCorrelations1D() {
     delete Corr1DHist_Reduced;
   }
 
-  for(int i = 0; i < nDraw; i++)
-  {
-    for(int k = 1; k < Nhists; k++)
-    {
-      delete Corr1DHist[i][k];
-    }
-    delete[] Corr1DHist[i];
-  }
-  delete[] Corr1DHist;
-
   CorrDir->Close();
   delete CorrDir;
 
@@ -1662,22 +1634,21 @@ void MCMCProcessor::MakeCredibleRegions(const std::vector<double>& CredibleRegio
     throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   const int nCredible = CredibleRegions.size();
-  TH2D*** hpost_2D_copy = new TH2D**[nDraw];
-  TH2D**** hpost_2D_cl = new TH2D***[nDraw];
 
+  std::vector<std::vector<TH2D*>> hpost_2D_copy(nDraw);
+  std::vector<std::vector<std::vector<TH2D*>>> hpost_2D_cl(nDraw);
   //KS: Copy all histograms to be thread safe
   for (int i = 0; i < nDraw; ++i)
   {
-    hpost_2D_copy[i] = new TH2D*[nDraw];
-    hpost_2D_cl[i] = new TH2D**[nDraw];
+    hpost_2D_copy[i].resize(nDraw);
+    hpost_2D_cl[i].resize(nDraw);
     for (int j = 0; j <= i; ++j)
     {
-      hpost_2D_copy[i][j] = (TH2D*) hpost2D[i][j]->Clone( Form("hpost_copy_%i_%i", i, j));
-
-      hpost_2D_cl[i][j] = new TH2D*[nCredible];
+      hpost_2D_copy[i][j] = static_cast<TH2D*>(hpost2D[i][j]->Clone(Form("hpost_copy_%i_%i", i, j)));
+      hpost_2D_cl[i][j].resize(nCredible);
       for (int k = 0; k < nCredible; ++k)
       {
-        hpost_2D_cl[i][j][k] = (TH2D*)hpost2D[i][j]->Clone( Form("hpost_copy_%i_%i_CL_%f", i, j, CredibleRegions[k]));;
+        hpost_2D_cl[i][j][k] = static_cast<TH2D*>(hpost2D[i][j]->Clone(Form("hpost_copy_%i_%i_CL_%f", i, j, CredibleRegions[k])));
       }
     }
   }
@@ -1728,7 +1699,7 @@ void MCMCProcessor::MakeCredibleRegions(const std::vector<double>& CredibleRegio
       legend->SetBorderSize(0);
 
       //Get Best point
-      TGraph *bestfitM = new TGraph(1);
+      auto bestfitM = std::make_unique<TGraph>(1);
       const int MaxBin = hpost_2D_copy[i][j]->GetMaximumBin();
       int Mbx, Mby, Mbz;
       hpost_2D_copy[i][j]->GetBinXYZ(MaxBin, Mbx, Mby, Mbz);
@@ -1739,7 +1710,7 @@ void MCMCProcessor::MakeCredibleRegions(const std::vector<double>& CredibleRegio
       bestfitM->SetMarkerStyle(22);
       bestfitM->SetMarkerSize(1);
       bestfitM->SetMarkerColor(kMagenta);
-      legend->AddEntry(bestfitM,"Best Fit","p");
+      legend->AddEntry(bestfitM.get(),"Best Fit","p");
 
       //Plot default 2D posterior
       hpost_2D_copy[i][j]->Draw("COLZ");
@@ -1754,7 +1725,6 @@ void MCMCProcessor::MakeCredibleRegions(const std::vector<double>& CredibleRegio
         else
           legend->AddEntry(hpost_2D_cl[i][j][k], Form("%.0f%% Credible Region", CredibleRegions[k]*100), "l");
       }
-
       legend->Draw("SAME");
       bestfitM->Draw("SAME.P");
 
@@ -1767,8 +1737,6 @@ void MCMCProcessor::MakeCredibleRegions(const std::vector<double>& CredibleRegio
       // Write it to root file
       //OutputFile->cd();
       //if( std::fabs((*Correlation)(i,j)) > Post2DPlotThreshold ) Posterior->Write();
-
-      delete bestfitM;
     }
   }
 
@@ -1783,13 +1751,8 @@ void MCMCProcessor::MakeCredibleRegions(const std::vector<double>& CredibleRegio
       {
         delete hpost_2D_cl[i][j][k];
       }
-      delete[] hpost_2D_cl[i][j];
     }
-    delete[] hpost_2D_copy[i];
-    delete[] hpost_2D_cl[i];
   }
-  delete[] hpost_2D_copy;
-  delete[] hpost_2D_cl;
 }
 
 // *********************
@@ -1964,13 +1927,13 @@ void MCMCProcessor::MakeTrianglePlot(const std::vector<std::string>& ParNames,
       //KS:if diagonal plot main posterior
       if(x == y)
       {
-        hpost_copy[counterPost] = (TH1D*) hpost[ParamNumber[x]]->Clone(Form("hpost_copy_%i", ParamNumber[x]));
+        hpost_copy[counterPost] = static_cast<TH1D*>(hpost[ParamNumber[x]]->Clone(Form("hpost_copy_%i", ParamNumber[x])));
         hpost_cl[counterPost] = new TH1D*[nCredibleIntervals];
         /// Scale the histograms so it shows the posterior probability
         hpost_copy[counterPost]->Scale(1. / hpost_copy[counterPost]->Integral());
         for (int j = 0; j < nCredibleIntervals; ++j)
         {
-          hpost_cl[counterPost][j] = (TH1D*) hpost[ParamNumber[x]]->Clone( Form("hpost_copy_%i_CL_%f", ParamNumber[x], CredibleIntervals[j]));
+          hpost_cl[counterPost][j] = static_cast<TH1D*>(hpost[ParamNumber[x]]->Clone(Form("hpost_copy_%i_CL_%f", ParamNumber[x], CredibleIntervals[j])));
           //KS: Reset to get rid to TF1 otherwise we run into segfault :(
           hpost_cl[counterPost][j]->Reset("");
           hpost_cl[counterPost][j]->Fill(0.0, 0.0);
@@ -2016,12 +1979,14 @@ void MCMCProcessor::MakeTrianglePlot(const std::vector<std::string>& ParNames,
       //KS: Here we plot 2D credible regions
       else
       {
-        hpost_2D_copy[counter2DPost] = (TH2D*) hpost2D[ParamNumber[x]][ParamNumber[y]]->Clone( Form("hpost_copy_%i_%i", ParamNumber[x], ParamNumber[y]));
+        hpost_2D_copy[counter2DPost] = static_cast<TH2D*>(hpost2D[ParamNumber[x]][ParamNumber[y]]->Clone(
+          Form("hpost_copy_%i_%i", ParamNumber[x], ParamNumber[y])));
         hpost_2D_cl[counter2DPost] = new TH2D*[nCredibleRegions];
         //KS: Now copy for every credible region
         for (int k = 0; k < nCredibleRegions; ++k)
         {
-          hpost_2D_cl[counter2DPost][k] = (TH2D*)hpost2D[ParamNumber[x]][ParamNumber[y]]->Clone( Form("hpost_copy_%i_%i_CL_%f", ParamNumber[x], ParamNumber[y], CredibleRegions[k]));
+          hpost_2D_cl[counter2DPost][k] = static_cast<TH2D*>(hpost2D[ParamNumber[x]][ParamNumber[y]]->Clone(
+            Form("hpost_copy_%i_%i_CL_%f", ParamNumber[x], ParamNumber[y], CredibleRegions[k])));
 
           if(CredibleInSigmas)
           {
@@ -2141,7 +2106,6 @@ void MCMCProcessor::MakeTrianglePlot(const std::vector<std::string>& ParNames,
     }
     delete[] hpost_2D_cl[i];
   }
-
   delete[] hpost_copy;
   delete[] hpost_cl;
   delete[] hpost_2D_copy;
@@ -2175,7 +2139,7 @@ void MCMCProcessor::ScanInput() {
   UpperCut = nEntries+1;
 
   // Get the list of branches
-  TObjArray* brlis = (TObjArray*)(Chain->GetListOfBranches());
+  TObjArray* brlis = static_cast<TObjArray*>(Chain->GetListOfBranches());
 
   // Get the number of branches
   nBranches = brlis->GetEntries();
@@ -2192,7 +2156,7 @@ void MCMCProcessor::ScanInput() {
   for (int i = 0; i < nBranches; i++)
   {
     // Get the TBranch and its name
-    TBranch* br = (TBranch*)brlis->At(i);
+    TBranch* br = static_cast<TBranch*>(brlis->At(i));
     TString bname = br->GetName();
 
     //KS: Exclude parameter types
@@ -2290,7 +2254,7 @@ void MCMCProcessor::SetupOutput() {
   CanvasName.ReplaceAll("[","");
 
   // We fit with this Gaussian
-  Gauss = new TF1("gauss","[0]/sqrt(2.0*3.14159)/[2]*TMath::Exp(-0.5*pow(x-[1],2)/[2]/[2])",   -5, 5);
+  Gauss = new TF1("gauss","[0]/sqrt(2.0*3.14159)/[2]*TMath::Exp(-0.5*pow(x-[1],2)/[2]/[2])", -5, 5);
 
   // Declare the TVectors
   Covariance = new TMatrixDSym(nDraw);
@@ -2432,10 +2396,10 @@ void MCMCProcessor::FindInputFiles() {
   // Now read the MCMC file
   TFile *TempFile = new TFile(MCMCFile.c_str(), "open");
 
-  TDirectory* CovarianceFolder = (TDirectory*)TempFile->Get("CovarianceFolder");
+  TDirectory* CovarianceFolder = static_cast<TDirectory*>(TempFile->Get("CovarianceFolder"));
 
   // Get the settings for the MCMC
-  TMacro *Config = (TMacro*)(TempFile->Get("MaCh3_Config"));
+  TMacro* Config = static_cast<TMacro*>(TempFile->Get("MaCh3_Config"));
   if (Config == nullptr) {
     MACH3LOG_ERROR("Didn't find MaCh3_Config tree in MCMC file! {}", MCMCFile);
     TempFile->ls();
@@ -2459,7 +2423,7 @@ void MCMCProcessor::FindInputFiles() {
     InputNotFound = true;
   }
 
-  TMacro *XsecConfig = (TMacro*)(CovarianceFolder->Get("Config_xsec_cov"));
+  TMacro *XsecConfig = static_cast<TMacro*>(CovarianceFolder->Get("Config_xsec_cov"));
   if (XsecConfig == nullptr) {
     MACH3LOG_WARN("Didn't find Config_xsec_cov tree in MCMC file! {}", MCMCFile);
   } else {
@@ -2580,10 +2544,9 @@ void MCMCProcessor::ReadNDFile() {
   }
   NDdetFile->cd();
 
-  TMatrixDSym *NDdetMatrix = (TMatrixDSym*)(NDdetFile->Get("nddet_cov"));
-  TVectorD *NDdetNominal = (TVectorD*)(NDdetFile->Get("det_weights"));
-  TDirectory *BinningDirectory = (TDirectory*)NDdetFile->Get("Binning")->Clone();
-
+  TMatrixDSym *NDdetMatrix = static_cast<TMatrixDSym*>(NDdetFile->Get("nddet_cov"));
+  TVectorD *NDdetNominal = static_cast<TVectorD*>(NDdetFile->Get("det_weights"));
+  TDirectory *BinningDirectory = static_cast<TDirectory*>(NDdetFile->Get("Binning")->Clone());
   for (int i = 0; i < NDdetNominal->GetNrows(); ++i)
   {
     ParamNom[kNDPar].push_back( (*NDdetNominal)(i) );
@@ -2597,7 +2560,6 @@ void MCMCProcessor::ReadNDFile() {
 
   TIter next(BinningDirectory->GetListOfKeys());
   TKey *key = nullptr;
-
   // Loop through all entries
   while ((key = (TKey*)next()))
   {
@@ -2829,7 +2791,7 @@ void MCMCProcessor::GetPolarPlot(const std::vector<std::string>& ParNames){
       y_val[ipt] = hpost[ParamNo]->GetBinContent(ipt+1)/Integral;
     }
 
-    TGraphPolar* PolarGraph = new TGraphPolar(nBins, x_val.data(), y_val.data());
+    auto PolarGraph = std::make_unique<TGraphPolar>(nBins, x_val.data(), y_val.data());
     PolarGraph->SetLineWidth(2);
     PolarGraph->SetFillStyle(3001);
     PolarGraph->SetLineColor(kRed);
@@ -2843,8 +2805,6 @@ void MCMCProcessor::GetPolarPlot(const std::vector<std::string>& ParNames){
 
     Posterior->Print(CanvasName);
     Posterior->Write(Title);
-
-    delete PolarGraph;
   } //End loop over parameters
 
   PolarDir->Close();
@@ -3025,13 +2985,13 @@ void MCMCProcessor::GetSavageDickey(const std::vector<std::string>& ParNames,
     
     std::string DunneKabothScale = GetDunneKaboth(SavageDickey);
     //Get Best point
-    TGraph *PostPoint = new TGraph(1);
+    std::unique_ptr<TGraph> PostPoint(new TGraph(1));
     PostPoint->SetPoint(0, EvaluationPoint[k], ProbPosterior);
     PostPoint->SetMarkerStyle(20);
     PostPoint->SetMarkerSize(1);
     PostPoint->Draw("P same");
     
-    TGraph *PriorPoint = new TGraph(1);
+    std::unique_ptr<TGraph> PriorPoint(new TGraph(1));
     PriorPoint->SetPoint(0, EvaluationPoint[k], ProbPrior);
     PriorPoint->SetMarkerStyle(20);
     PriorPoint->SetMarkerSize(1);
@@ -3041,7 +3001,7 @@ void MCMCProcessor::GetSavageDickey(const std::vector<std::string>& ParNames,
     legend->SetTextSize(0.04);
     legend->AddEntry(PriorHist, "Prior", "l");
     legend->AddEntry(PosteriorHist, "Posterior", "l");
-    legend->AddEntry(PostPoint, Form("SavageDickey = %.2f, (%s)", SavageDickey, DunneKabothScale.c_str()),"");
+    legend->AddEntry(PostPoint.get(), Form("SavageDickey = %.2f, (%s)", SavageDickey, DunneKabothScale.c_str()),"");
     legend->SetLineColor(0);
     legend->SetLineStyle(0);
     legend->SetFillColor(0);
@@ -3054,8 +3014,6 @@ void MCMCProcessor::GetSavageDickey(const std::vector<std::string>& ParNames,
     
     delete PosteriorHist;
     delete PriorHist;
-    delete PostPoint;
-    delete PriorPoint;
   } //End loop over parameters
 
   SavageDickeyDir->Close();
@@ -3934,7 +3892,6 @@ void MCMCProcessor::BatchedMeans() {
 // Get the batched means variance estimation and variable indicating if number of batches is sensible
 void MCMCProcessor::BatchedAnalysis() {
 // **************************
-
   if(BatchedAverages == nullptr)
   {
     MACH3LOG_ERROR("BatchedAverages haven't been initialises or have been deleted something is wrong");
@@ -4193,7 +4150,7 @@ void MCMCProcessor::GewekeDiagnostic() {
   constexpr int NChecks = 100;
   constexpr double Division = (UpperThreshold - LowerThreshold)/NChecks;
 
-  TH1D** GewekePlots = new TH1D*[nDraw];
+  std::vector<std::unique_ptr<TH1D>> GewekePlots(nDraw);
   for (int j = 0; j < nDraw; ++j)
   {
     TString Title = "";
@@ -4201,7 +4158,7 @@ void MCMCProcessor::GewekeDiagnostic() {
     double PriorError = 1.0;
     GetNthParameter(j, Prior, PriorError, Title);
     std::string HistName = Form("%s_%s_Geweke", Title.Data(), BranchNames[j].Data());
-    GewekePlots[j] = new TH1D(HistName.c_str(), HistName.c_str(), NChecks, 0.0, 100*UpperThreshold);
+    GewekePlots[j] = std::make_unique<TH1D>(HistName.c_str(), HistName.c_str(), NChecks, 0.0, 100 * UpperThreshold);
     GewekePlots[j]->GetXaxis()->SetTitle("Burn-In (%)");
     GewekePlots[j]->GetYaxis()->SetTitle("Geweke T score");
   }
@@ -4283,7 +4240,6 @@ void MCMCProcessor::GewekeDiagnostic() {
     for (int j = 0; j < nDraw; ++j)
     {
       double T_score = std::fabs((MeanDown[j] - MeanUp[j])/std::sqrt(SpectralVarianceDown[j]/DenomCounterDown[j] + SpectralVarianceUp[j]/DenomCounterUp[j]));
-
       GewekePlots[j]->SetBinContent(k, T_score);
     }
   } //end loop over intervals
@@ -4298,9 +4254,7 @@ void MCMCProcessor::GewekeDiagnostic() {
   {
     GewekeDir->cd();
     GewekePlots[j]->Write();
-    delete GewekePlots[j];
   }
-  delete[] GewekePlots;
   for (int i = 0; i < nDraw; ++i) {
     delete[] ParStep[i];
   }
diff --git a/mcmc/PSO.cpp b/mcmc/PSO.cpp
index 479f3323..1fe496c5 100644
--- a/mcmc/PSO.cpp
+++ b/mcmc/PSO.cpp
@@ -2,512 +2,507 @@
 
 PSO::PSO(manager *man) : LikelihoodFit(man) {
 
-    fConstriction = fitMan->raw()["General"]["PSO"]["Constriction"].as<double>();
-    fInertia = fitMan->raw()["General"]["PSO"]["Inertia"].as<double>()*fConstriction;
-    fOne = fitMan->raw()["General"]["PSO"]["One"].as<double>()*fConstriction;
-    fTwo = fitMan->raw()["General"]["PSO"]["Two"].as<double>()*fConstriction;
-    fParticles = fitMan->raw()["General"]["PSO"]["Particles"].as<int>();
-    fIterations = fitMan->raw()["General"]["PSO"]["Iterations"].as<int>();
-    fConvergence = fitMan->raw()["General"]["PSO"]["Convergence"].as<double>();
+  fConstriction = fitMan->raw()["General"]["PSO"]["Constriction"].as<double>();
+  fInertia = fitMan->raw()["General"]["PSO"]["Inertia"].as<double>()*fConstriction;
+  fOne = fitMan->raw()["General"]["PSO"]["One"].as<double>()*fConstriction;
+  fTwo = fitMan->raw()["General"]["PSO"]["Two"].as<double>()*fConstriction;
+  fParticles = fitMan->raw()["General"]["PSO"]["Particles"].as<int>();
+  fIterations = fitMan->raw()["General"]["PSO"]["Iterations"].as<int>();
+  fConvergence = fitMan->raw()["General"]["PSO"]["Convergence"].as<double>();
 
-    fDim = 0;
+  fDim = 0;
 
-    if(fTestLikelihood)
-    {
-      fDim = fitMan->raw()["General"]["PSO"]["TestLikelihoodDim"].as<int>();
-    }
+  if(fTestLikelihood)
+  {
+    fDim = fitMan->raw()["General"]["PSO"]["TestLikelihoodDim"].as<int>();
+  }
 }
 
 void PSO::runMCMC(){
 
-    PrepareFit();
+  PrepareFit();
 
-    if(fTestLikelihood){ 
-        outTree->Branch("nParts", &fParticles, "nParts/I");
-        for(int i = 0; i < fDim; ++i){
-            par = new double[fParticles];
-            paramlist.push_back(par);
-            outTree->Branch(Form("Parameter_%d", i), paramlist[i], Form("Parameter_%d[nParts]/D",i));
-        }
-//        vel = new double[fParticles];
-        outTree->Branch("vel", vel, "vel[nParts]/D");
+  if(fTestLikelihood){
+    outTree->Branch("nParts", &fParticles, "nParts/I");
+    for(int i = 0; i < fDim; ++i){
+      par = new double[fParticles];
+      paramlist.push_back(par);
+      outTree->Branch(Form("Parameter_%d", i), paramlist[i], Form("Parameter_%d[nParts]/D",i));
     }
+//  vel = new double[fParticles];
+    outTree->Branch("vel", vel, "vel[nParts]/D");
+  }
 
-    init();
-    run();
-    WriteOutput();
-    return;
+  init();
+  run();
+  WriteOutput();
+  return;
 }
 
 
 void PSO::init(){
 
-    fBestValue = 1234567890.0;
+  fBestValue = 1234567890.0;
 
-    //KS: For none PCA this will be eqaul to normal parameters
-    //const int NparsPSOFull = NPars;
-    //const int NparsPSO = NParsPCA;
+  //KS: For none PCA this will be eqaul to normal parameters
+  //const int NparsPSOFull = NPars;
+  //const int NparsPSO = NParsPCA;
 
-    std::cout << "Preparing PSO" << std::endl;
+  std::cout << "Preparing PSO" << std::endl;
 
-    // Initialise bounds on parameters
-    if(fTestLikelihood){
-        for (int i = 0; i < fDim; i++){
-            // Test function ranges
-            ranges_min.push_back(-5);
-            ranges_max.push_back(5);
+  // Initialise bounds on parameters
+  if(fTestLikelihood){
+    for (int i = 0; i < fDim; i++){
+      // Test function ranges
+      ranges_min.push_back(-5);
+      ranges_max.push_back(5);
+      fixed.push_back(0);
+    }
+  }
+  else{
+    for (std::vector<covarianceBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it){
+      if(!(*it)->IsPCA())
+      {
+        fDim += (*it)->GetNumParams();
+        for(int i = 0; i < (*it)->GetNumParams(); ++i)
+        {
+          double curr = (*it)->getParInit(i);
+          double lim = 10.0*(*it)->getDiagonalError(i);
+          double low = (*it)->GetLowerBound(i);
+          double high = (*it)->GetUpperBound(i);
+          if(low > curr - lim) ranges_min.push_back(low);
+          else ranges_min.push_back(curr - lim);
+          if(high < curr + lim) ranges_min.push_back(high);
+          else ranges_min.push_back(curr + lim);
+          prior.push_back(curr);
+
+          if((*it)->isParameterFixed(i)){
+            fixed.push_back(1);
+          }
+          else{
             fixed.push_back(0);
+          }
         }
-    }
-    else{
-        for (std::vector<covarianceBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it){
-            if(!(*it)->IsPCA())
-            {
-                fDim += (*it)->GetNumParams();
-                for(int i = 0; i < (*it)->GetNumParams(); ++i)
-                {
-                    double curr = (*it)->getParInit(i);
-                    double lim = 10.0*(*it)->getDiagonalError(i);
-                    double low = (*it)->GetLowerBound(i);
-                    double high = (*it)->GetUpperBound(i);
-                    if(low > curr - lim) ranges_min.push_back(low);
-                    else ranges_min.push_back(curr - lim);
-                    if(high < curr + lim) ranges_min.push_back(high);
-                    else ranges_min.push_back(curr + lim);
-                    prior.push_back(curr);
-
-                    if((*it)->isParameterFixed(i)){
-                        fixed.push_back(1);
-                    }
-                    else{
-                        fixed.push_back(0);
-                    }
-                }
-            }
-            else
-            {
-                fDim += (*it)->getNpars();
-                for(int i = 0; i < (*it)->getNpars(); ++i)
-                {
-                    ranges_min.push_back(-100.0);
-                    ranges_max.push_back(100.0);
-                    prior.push_back((*it)->getParInit(i));
-                    if((*it)->isParameterFixedPCA(i)){
-                        fixed.push_back(1);
-                    }
-                    else{
-                        fixed.push_back(0);
-                    }
-                }
-            }
+      }
+      else
+      {
+        fDim += (*it)->getNpars();
+        for(int i = 0; i < (*it)->getNpars(); ++i)
+        {
+          ranges_min.push_back(-100.0);
+          ranges_max.push_back(100.0);
+          prior.push_back((*it)->getParInit(i));
+          if((*it)->isParameterFixedPCA(i)){
+            fixed.push_back(1);
+          }
+          else{
+            fixed.push_back(0);
+          }
         }
+      }
     }
+  }
 
-    std::cout << "Printing Minimums and Maximums of Variables to be minimised" << std::endl;
-    for (int i =0; i<fDim; i++){
-        std::cout << "Variable "<< i<<" :" << ranges_min[i] << ", "<< ranges_max[i] << std::endl;
-    }
+  std::cout << "Printing Minimums and Maximums of Variables to be minimised" << std::endl;
+  for (int i =0; i<fDim; i++){
+    std::cout << "Variable "<< i<<" :" << ranges_min[i] << ", "<< ranges_max[i] << std::endl;
+  }
 
-    // Initialise particle positions
-    for (int i = 0; i < fParticles; ++i){
-        std::vector<double> init_position;
-        std::vector<double> init_velocity;
-
-        //Initialising in +/- 5sigma of prior value from BANFF interface
-        for (int j=0; j<fDim; ++j){
-            if(fixed[j]){
-                init_position.push_back(prior[j]);
-                init_velocity.push_back(0.0);
-            }
-            else{
-                double dist = fabs(ranges_max[j]-ranges_min[j]);
-                //Initialise to random position uniform in space
-                init_position.push_back(ranges_min[j] + random->Rndm()*dist);
-                //Initialise velocity to random position uniform in space
-                init_velocity.push_back((2.0*random->Rndm()-1.0));//*dist);
-            }
-        }
+  // Initialise particle positions
+  for (int i = 0; i < fParticles; ++i){
+    std::vector<double> init_position;
+    std::vector<double> init_velocity;
+
+    //Initialising in +/- 5sigma of prior value from BANFF interface
+    for (int j=0; j<fDim; ++j){
+      if(fixed[j]){
+        init_position.push_back(prior[j]);
+        init_velocity.push_back(0.0);
+      }
+      else{
+        double dist = fabs(ranges_max[j]-ranges_min[j]);
+        //Initialise to random position uniform in space
+        init_position.push_back(ranges_min[j] + random->Rndm()*dist);
+        //Initialise velocity to random position uniform in space
+        init_velocity.push_back((2.0*random->Rndm()-1.0));//*dist);
+      }
+    }
 
-        particle* new_particle = new particle(init_position,init_velocity);
-        new_particle->set_personal_best_position(init_position);
-        double new_value = CalcChi(init_position);
-        new_particle->set_personal_best_value(new_value);
-        new_particle->set_value(new_value);
-        system.push_back(new_particle);
-        if(new_value < fBestValue){
-            fBestValue = new_value;
-            set_best_particle(new_particle);
-        }
+    particle* new_particle = new particle(init_position,init_velocity);
+    new_particle->set_personal_best_position(init_position);
+    double new_value = CalcChi(init_position);
+    new_particle->set_personal_best_value(new_value);
+    new_particle->set_value(new_value);
+    system.push_back(new_particle);
+    if(new_value < fBestValue){
+      fBestValue = new_value;
+      set_best_particle(new_particle);
     }
+  }
 }
 
 std::vector<std::vector<double> > PSO::bisection(std::vector<double>position,double minimum, double range, double precision){
-    std::vector<std::vector<double>> uncertainties_list;
-    for (unsigned int i = 0; i< position.size(); ++i){
-        std::cout << i << std::endl;
-        //std::vector<double> uncertainties;
-        std::vector<double> new_position = position; new_position[i] = position[i]-range;
-        double val_1 = CalcChi(new_position)-minimum-1.0;
-        while (val_1*-1.0> 0.0){
-            new_position[i] -= range;
-            val_1 = CalcChi(new_position)-minimum-1.0;
-        }
-        std::vector<double> bisect_position = position; bisect_position[i] = bisect_position[i] - (position[i]-new_position[i])/2;
-        std::vector<std::vector<double>> position_list{new_position,bisect_position,position};
-        double val_2 = CalcChi(bisect_position)-minimum-1.0;
-        std::vector<double> value_list{val_1,val_2, -1.0};
-        double res = 1.0;
-        while (res > precision){
-            if (value_list[0] * value_list[1] < 0){
-                std::vector<double> new_bisect_position = position_list[0];new_bisect_position[i] =new_bisect_position[i]+ (position_list[1][i]-position_list[0][i])/2;
-                double new_val = CalcChi(new_bisect_position)-minimum-1.0;
-                position_list[2] = position_list[1];
-                value_list[2] = value_list[1];
-                value_list[1] = new_val;
-                position_list[1] = new_bisect_position;
-                res = abs(position[2]-position[0]);
-            }
-            else{
-                std::vector<double> new_bisect_position = position_list[1];new_bisect_position[i] += (position_list[2][i]-position_list[1][i])/2;
-                double new_val = CalcChi(new_bisect_position)-minimum-1.0;
-                position_list[0] = position_list[1];
-                value_list[0] = value_list[1];
-                value_list[1] = new_val;
-                position_list[1] = new_bisect_position;
-                res = abs(position_list[2][i]-position_list[1][i]);
-            }
-        }
-        //do the same thing for position uncertainty
-        std::vector<double> new_position_p = position; new_position_p[i] = position[i]+range;
-        double val_1_p = CalcChi(new_position_p)-minimum-1.0;
-        while (val_1_p * -1.0 > 0.0){
-            new_position_p[i] += range;
-            val_1_p = CalcChi(new_position_p)-minimum-1.0;
-        }
-        std::vector<double> bisect_position_p = position; bisect_position_p[i] = bisect_position_p[i] += (new_position_p[i]-position[i])/2;
-        std::vector<std::vector<double>> position_list_p{position,bisect_position_p,new_position_p};
-        double val_2_p = CalcChi(bisect_position_p)-minimum-1.0;
-        std::vector<double> value_list_p{-1.0,val_2_p, val_1_p};
-        double res_p = 1.0;
-        while (res_p > precision){
-            if (value_list_p[0] * value_list_p[1] < 0){
-                std::vector<double> new_bisect_position_p = position_list_p[0];new_bisect_position_p[i] += (position_list_p[1][i]-position_list_p[0][i])/2;
-                double new_val_p = CalcChi(new_bisect_position_p)-minimum-1.0;
-                position_list_p[2] = position_list_p[1];
-                value_list_p[2] = value_list_p[1];
-                value_list_p[1] = new_val_p;
-                position_list_p[1] = new_bisect_position_p;
-                res = abs(position[2]-position[0]);
-                res_p = abs(position_list_p[1][i]-position_list_p[0][i]);
-                //std::cout << "Pos midpoint is " << position_list_p[1][i] << std::endl;
-            }
-            else{
-                std::vector<double> new_bisect_position_p = position_list_p[1];new_bisect_position_p[i] += (position_list_p[2][i]-position_list_p[1][i])/2;
-                double new_val_p = CalcChi(new_bisect_position_p)-minimum-1.0;
-                position_list_p[0] = position_list_p[1];
-                value_list_p[0] = value_list_p[1];
-                value_list_p[1] = new_val_p;
-                position_list_p[1] = new_bisect_position_p;
-                res_p = abs(position_list_p[2][i]-position_list_p[1][i]);
-                //std::cout << "Pos midpoint is " << position_list_p[1][i] << std::endl;
-            }
-        }
-        uncertainties_list.push_back({abs(position[i]-position_list[1][i]),abs(position[i]-position_list_p[1][i])});
-        std::cout << "Uncertainty finished for d = "<< i << std::endl;
-        std::cout << std::setprecision(10)<< "LLR values for ± positive and negative uncertainties are " << CalcChi(position_list[1]) << " and " << CalcChi(position_list_p[1]) << std::endl;
+  std::vector<std::vector<double>> uncertainties_list;
+  for (unsigned int i = 0; i< position.size(); ++i){
+    std::cout << i << std::endl;
+    //std::vector<double> uncertainties;
+    std::vector<double> new_position = position; new_position[i] = position[i]-range;
+    double val_1 = CalcChi(new_position)-minimum-1.0;
+    while (val_1*-1.0> 0.0){
+      new_position[i] -= range;
+      val_1 = CalcChi(new_position)-minimum-1.0;
     }
-    return uncertainties_list;
+    std::vector<double> bisect_position = position; bisect_position[i] = bisect_position[i] - (position[i]-new_position[i])/2;
+    std::vector<std::vector<double>> position_list{new_position,bisect_position,position};
+    double val_2 = CalcChi(bisect_position)-minimum-1.0;
+    std::vector<double> value_list{val_1,val_2, -1.0};
+    double res = 1.0;
+    while (res > precision){
+      if (value_list[0] * value_list[1] < 0){
+        std::vector<double> new_bisect_position = position_list[0];new_bisect_position[i] =new_bisect_position[i]+ (position_list[1][i]-position_list[0][i])/2;
+        double new_val = CalcChi(new_bisect_position)-minimum-1.0;
+        position_list[2] = position_list[1];
+        value_list[2] = value_list[1];
+        value_list[1] = new_val;
+        position_list[1] = new_bisect_position;
+        res = abs(position[2]-position[0]);
+      }
+      else{
+        std::vector<double> new_bisect_position = position_list[1];new_bisect_position[i] += (position_list[2][i]-position_list[1][i])/2;
+        double new_val = CalcChi(new_bisect_position)-minimum-1.0;
+        position_list[0] = position_list[1];
+        value_list[0] = value_list[1];
+        value_list[1] = new_val;
+        position_list[1] = new_bisect_position;
+        res = abs(position_list[2][i]-position_list[1][i]);
+      }
+    }
+    //do the same thing for position uncertainty
+    std::vector<double> new_position_p = position; new_position_p[i] = position[i]+range;
+    double val_1_p = CalcChi(new_position_p)-minimum-1.0;
+    while (val_1_p * -1.0 > 0.0){
+      new_position_p[i] += range;
+      val_1_p = CalcChi(new_position_p)-minimum-1.0;
+    }
+    std::vector<double> bisect_position_p = position; bisect_position_p[i] = bisect_position_p[i] += (new_position_p[i]-position[i])/2;
+    std::vector<std::vector<double>> position_list_p{position,bisect_position_p,new_position_p};
+    double val_2_p = CalcChi(bisect_position_p)-minimum-1.0;
+    std::vector<double> value_list_p{-1.0,val_2_p, val_1_p};
+    double res_p = 1.0;
+    while (res_p > precision){
+      if (value_list_p[0] * value_list_p[1] < 0){
+        std::vector<double> new_bisect_position_p = position_list_p[0];new_bisect_position_p[i] += (position_list_p[1][i]-position_list_p[0][i])/2;
+        double new_val_p = CalcChi(new_bisect_position_p)-minimum-1.0;
+        position_list_p[2] = position_list_p[1];
+        value_list_p[2] = value_list_p[1];
+        value_list_p[1] = new_val_p;
+        position_list_p[1] = new_bisect_position_p;
+        res = abs(position[2]-position[0]);
+        res_p = abs(position_list_p[1][i]-position_list_p[0][i]);
+        //std::cout << "Pos midpoint is " << position_list_p[1][i] << std::endl;
+      }
+      else{
+        std::vector<double> new_bisect_position_p = position_list_p[1];new_bisect_position_p[i] += (position_list_p[2][i]-position_list_p[1][i])/2;
+        double new_val_p = CalcChi(new_bisect_position_p)-minimum-1.0;
+        position_list_p[0] = position_list_p[1];
+        value_list_p[0] = value_list_p[1];
+        value_list_p[1] = new_val_p;
+        position_list_p[1] = new_bisect_position_p;
+        res_p = abs(position_list_p[2][i]-position_list_p[1][i]);
+        //std::cout << "Pos midpoint is " << position_list_p[1][i] << std::endl;
+      }
+    }
+    uncertainties_list.push_back({abs(position[i]-position_list[1][i]),abs(position[i]-position_list_p[1][i])});
+    std::cout << "Uncertainty finished for d = "<< i << std::endl;
+    std::cout << std::setprecision(10)<< "LLR values for ± positive and negative uncertainties are " << CalcChi(position_list[1]) << " and " << CalcChi(position_list_p[1]) << std::endl;
+  }
+  return uncertainties_list;
 }
 
 std::vector<std::vector<double>> PSO::calc_uncertainty(std::vector<double>position,double minimum) {
-    std::vector<double> pos_uncertainty(position.size());
-    std::vector<double> neg_uncertainty(position.size());
-    int num = 200;
-    std::vector<double> pos = position;
-    for (unsigned int i = 0; i < position.size(); ++i) {
-
-        double curr_ival = pos[i];
-
-        double neg_stop = position[i] - 5e-2;
-        double pos_stop = position[i] + 5e-2;
-        double start = position[i];
-        std::vector<double> x(num);
-        std::vector<double> y(num);
-        double StepPoint = (start-neg_stop) / (num - 1);
-        double value = start;
-        for (int j = 0; j < num; ++j) {
-          pos[i] = value;
-          double LLR = CalcChi(position) - minimum - 1.0;
-          x[j] = value;
-          y[j] = LLR;
-          value -= StepPoint;
-        }
-        pos[i] = curr_ival;
-
-        int closest_index = 0;
-        double closest_value = abs(y[0]); // Initialize with the first element
-        for (unsigned int ii = 1; ii < y.size(); ++ii) {
-          double abs_y = abs(y[ii]);
-          if (abs_y < closest_value) {
-            closest_index = ii;
-            closest_value = abs_y;
-          }
-        }
-        neg_uncertainty[i] = x[closest_index];
-        std::cout << "Neg" << std::endl;
-        x.assign(num, 0);
-        y.assign(num, 0);
-        StepPoint = (pos_stop-start) / (num - 1);
-        value = start;
-        for (int j = 0; j < num; ++j) {
-          pos[i] = value;
-          double LLR = CalcChi(position) - minimum - 1.0;
-          x[j] = value;
-          y[j] = LLR;
-          value += StepPoint;
-        }
-        pos[i] = curr_ival;
-        closest_index = 0;
-        closest_value = abs(y[0]); // Initialize with the first element
-        for (unsigned int ii = 1; ii < y.size(); ++ii) {
-          double abs_y = abs(y[ii]);
-          if (abs_y < closest_value) {
-            closest_index = ii;
-            closest_value = abs_y;
-          }
-        } 
-        pos_uncertainty[i] = x[closest_index];
-    } 
-    std::vector<std::vector<double>> res{neg_uncertainty,pos_uncertainty};
-    return res;
+  std::vector<double> pos_uncertainty(position.size());
+  std::vector<double> neg_uncertainty(position.size());
+  int num = 200;
+  std::vector<double> pos = position;
+  for (unsigned int i = 0; i < position.size(); ++i) {
+    double curr_ival = pos[i];
+
+    double neg_stop = position[i] - 5e-2;
+    double pos_stop = position[i] + 5e-2;
+    double start = position[i];
+    std::vector<double> x(num);
+    std::vector<double> y(num);
+    double StepPoint = (start-neg_stop) / (num - 1);
+    double value = start;
+    for (int j = 0; j < num; ++j) {
+      pos[i] = value;
+      double LLR = CalcChi(position) - minimum - 1.0;
+      x[j] = value;
+      y[j] = LLR;
+      value -= StepPoint;
+    }
+    pos[i] = curr_ival;
+
+    int closest_index = 0;
+    double closest_value = abs(y[0]); // Initialize with the first element
+    for (unsigned int ii = 1; ii < y.size(); ++ii) {
+      double abs_y = abs(y[ii]);
+      if (abs_y < closest_value) {
+        closest_index = ii;
+        closest_value = abs_y;
+      }
+    }
+    neg_uncertainty[i] = x[closest_index];
+    std::cout << "Neg" << std::endl;
+    x.assign(num, 0);
+    y.assign(num, 0);
+    StepPoint = (pos_stop-start) / (num - 1);
+    value = start;
+    for (int j = 0; j < num; ++j) {
+      pos[i] = value;
+      double LLR = CalcChi(position) - minimum - 1.0;
+      x[j] = value;
+      y[j] = LLR;
+      value += StepPoint;
+    }
+    pos[i] = curr_ival;
+    closest_index = 0;
+    closest_value = abs(y[0]); // Initialize with the first element
+    for (unsigned int ii = 1; ii < y.size(); ++ii) {
+      double abs_y = abs(y[ii]);
+      if (abs_y < closest_value) {
+        closest_index = ii;
+        closest_value = abs_y;
+      }
+    }
+    pos_uncertainty[i] = x[closest_index];
+  }
+  std::vector<std::vector<double>> res{neg_uncertainty,pos_uncertainty};
+  return res;
 }
 
 void PSO::uncertainty_check(std::vector<double> previous_pos){
-    std::vector<std::vector<double >> x_list;
-    std::vector<std::vector<double >> y_list;
-    std::vector<double> position = previous_pos;
-    int num = 5000;
-    for (unsigned int i = 0;i<previous_pos.size();++i){
-        double curr_ival = position[i];
-        double start = previous_pos[i] - 1e-1;
-        double stop = previous_pos[i] + 1e-1;
-        std::vector<double> x(num);
-        std::vector<double> y(num);
-        double StepPoint = (stop - start) / (num - 1);
-        double value = start;
-        // std::cout << "result for fDim " << 1 << std::endl;
-        for (int j =0;j< num; ++j){
-            position[i] = value;
-            double LLR = CalcChi(position);
-            x[j] = value;
-            y[j] = LLR;
-            value += StepPoint;
-        }
-        position[i] = curr_ival;
-        std::cout << " " << std::endl;
-        std::cout << "For fDim" << i+1  << " x list is " ;
-        for (unsigned int k= 0;k<x.size(); ++k){
-            std::cout << x[k] << " , " ;
-        } std::cout << " " << std::endl;
-        std::cout << "  " << std::endl;
-        std::cout << "For fDim" << i+1 << " y list is " ;
-        for (unsigned int k= 0;k<x.size(); ++k){
-            std::cout << y[k] << " , " ;
-        } std::cout <<  " " << std::endl;
-        std::cout << " " <<std::endl;
+  std::vector<std::vector<double >> x_list;
+  std::vector<std::vector<double >> y_list;
+  std::vector<double> position = previous_pos;
+  int num = 5000;
+  for (unsigned int i = 0;i<previous_pos.size();++i){
+    double curr_ival = position[i];
+    double start = previous_pos[i] - 1e-1;
+    double stop = previous_pos[i] + 1e-1;
+    std::vector<double> x(num);
+    std::vector<double> y(num);
+    double StepPoint = (stop - start) / (num - 1);
+    double value = start;
+    // std::cout << "result for fDim " << 1 << std::endl;
+    for (int j =0;j< num; ++j){
+      position[i] = value;
+      double LLR = CalcChi(position);
+      x[j] = value;
+      y[j] = LLR;
+      value += StepPoint;
     }
+    position[i] = curr_ival;
+    std::cout << " " << std::endl;
+    std::cout << "For fDim" << i+1  << " x list is " ;
+    for (unsigned int k= 0;k<x.size(); ++k){
+      std::cout << x[k] << " , " ;
+    } std::cout << " " << std::endl;
+    std::cout << "  " << std::endl;
+    std::cout << "For fDim" << i+1 << " y list is " ;
+    for (unsigned int k= 0;k<x.size(); ++k){
+      std::cout << y[k] << " , " ;
+    } std::cout <<  " " << std::endl;
+    std::cout << " " <<std::endl;
+  }
 }
 
 double PSO::swarmIterate(){
+  std::vector<double> total_pos(fDim,0.0);
+
+  for (int i = 0; i < fParticles; ++i) {
+    std::vector<double> part1 = vector_multiply(system[i]->get_velocity(), fInertia);
+    std::vector<double> part2 = vector_multiply(vector_subtract(system[i]->get_personal_best_position(), system[i]->get_position()), (fOne * random->Rndm()));
+    std::vector<double> part3 = vector_multiply(vector_subtract(get_best_particle()->get_personal_best_position(), system[i]->get_position()),(fTwo * random->Rndm()));
+    std::vector<double> new_velocity = three_vector_addition(part1, part2, part3);
+    std::vector<double> new_pos = vector_add(system[i]->get_position(), new_velocity);
+    transform(total_pos.begin(), total_pos.end(), new_pos.begin(), total_pos.begin(),[](double x, double y) {return x+y;});
+
+    for (int j = 0; j < fDim; ++j) {
+      // Check if out of bounds and reflect if so
+      if(ranges_min[j] > new_pos[j]){
+          new_pos[j] = ranges_min[j];
+      }
+      else if(new_pos[j] > ranges_max[j]) {
+          new_pos[j] = ranges_max[j];
+      }
+      // If parameter fixed don't update it
+      if(fixed[j]) new_pos[j] = system[i]->get_position()[j];
+    }
 
-    std::vector<double> total_pos(fDim,0.0);
-
-    for (int i = 0; i < fParticles; ++i) {
-
-        std::vector<double> part1 = vector_multiply(system[i]->get_velocity(), fInertia);
-        std::vector<double> part2 = vector_multiply(vector_subtract(system[i]->get_personal_best_position(), system[i]->get_position()), (fOne * random->Rndm()));
-        std::vector<double> part3 = vector_multiply(vector_subtract(get_best_particle()->get_personal_best_position(), system[i]->get_position()),(fTwo * random->Rndm()));
-        std::vector<double> new_velocity = three_vector_addition(part1, part2, part3);
-        std::vector<double> new_pos = vector_add(system[i]->get_position(), new_velocity);
-        transform(total_pos.begin(), total_pos.end(), new_pos.begin(), total_pos.begin(),[](double x, double y) {return x+y;});
-
-        for (int j = 0; j < fDim; ++j) {
-            // Check if out of bounds and reflect if so
-            if(ranges_min[j] > new_pos[j]){
-                new_pos[j] = ranges_min[j];
-            }
-            else if(new_pos[j] > ranges_max[j]) {
-                new_pos[j] = ranges_max[j];
-            }
-            // If parameter fixed don't update it
-            if(fixed[j]) new_pos[j] = system[i]->get_position()[j];
-        }
-
-        if(fTestLikelihood){
-            double velo = 0.0;
-            for (int j = 0; j < fDim; ++j) {
-                paramlist[j][i] = new_pos[j];
-                velo += new_velocity[j]*new_velocity[j];
-            }
-            vel[i] = sqrt(velo);
-        }
+    if(fTestLikelihood){
+      double velo = 0.0;
+      for (int j = 0; j < fDim; ++j) {
+        paramlist[j][i] = new_pos[j];
+        velo += new_velocity[j]*new_velocity[j];
+      }
+      vel[i] = sqrt(velo);
+    }
 
-        system[i]->set_velocity(new_velocity);
-        system[i]->set_position(new_pos);
-        double new_value = CalcChi(new_pos);
-        if(new_value <= system[i]->get_personal_best_value()) {
-            system[i]->set_personal_best_value(new_value);
-            system[i]->set_personal_best_position(new_pos);
-            if(new_value < fBestValue){
-                fBestValue = new_value;
-                set_best_particle(system[i]);
-            }
-        }
+    system[i]->set_velocity(new_velocity);
+    system[i]->set_position(new_pos);
+    double new_value = CalcChi(new_pos);
+    if(new_value <= system[i]->get_personal_best_value()) {
+      system[i]->set_personal_best_value(new_value);
+      system[i]->set_personal_best_position(new_pos);
+      if(new_value < fBestValue){
+        fBestValue = new_value;
+        set_best_particle(system[i]);
+      }
     }
+  }
 
-    std::vector<double> best_pos = get_best_particle()->get_personal_best_position();
-    std::vector<double> result(best_pos.size(), 0.0);
-    transform(total_pos.begin(), total_pos.end(), total_pos.begin(), [=](double val){return val/fParticles;});
-    transform(total_pos.begin(),total_pos.end(),best_pos.begin(),result.begin(),[](double x, double y) {return x-y;});
+  std::vector<double> best_pos = get_best_particle()->get_personal_best_position();
+  std::vector<double> result(best_pos.size(), 0.0);
+  transform(total_pos.begin(), total_pos.end(), total_pos.begin(), [=](double val){return val/fParticles;});
+  transform(total_pos.begin(),total_pos.end(),best_pos.begin(),result.begin(),[](double x, double y) {return x-y;});
 
-    double mean_dist_sq = 0;
-    for (int i = 0; i<fDim; i++){
-        mean_dist_sq += result[i]*result[i];
-    }
+  double mean_dist_sq = 0;
+  for (int i = 0; i<fDim; i++){
+    mean_dist_sq += result[i]*result[i];
+  }
 
-    return mean_dist_sq;
+  return mean_dist_sq;
 }
 
 void PSO::run() {
 
-    double mean_dist_sq = 0;
+  double mean_dist_sq = 0;
 
-    int iter = 0;
-    for(int i = 0; i < fIterations; ++i, ++iter){
-        mean_dist_sq = swarmIterate();
-        //double meanVel = std::accumulate(vel, vel + fParticles, 0) / fParticles;
+  int iter = 0;
+  for(int i = 0; i < fIterations; ++i, ++iter){
+    mean_dist_sq = swarmIterate();
+    //double meanVel = std::accumulate(vel, vel + fParticles, 0) / fParticles;
 
-        // Weight inertia randomly but scaled by total distance of swarm from global minimum - proxy for total velocity
-        // fWeight = ((random->Rndm()+1.0)*0.5)*(10.0/meanVel);
+    // Weight inertia randomly but scaled by total distance of swarm from global minimum - proxy for total velocity
+    // fWeight = ((random->Rndm()+1.0)*0.5)*(10.0/meanVel);
 
-        logLCurr = fBestValue;
+    logLCurr = fBestValue;
 
-        outTree->Fill();
-        // Auto save the output
-        if (step % auto_save == 0) outTree->AutoSave();
-        step++;
-        accCount++;
+    outTree->Fill();
+    // Auto save the output
+    if (step % auto_save == 0) outTree->AutoSave();
+    step++;
+    accCount++;
 
-        if (i%100 == 0){
-            std::cout << "Mean Dist Sq = " << mean_dist_sq <<std::endl;
-            std::cout << "Current LLR = " << fBestValue << std::endl;
-            std::cout << "Position = " <<std::endl;
-            for (int j = 0; j< fDim; ++j){
-                std::cout << "    Dim " << j << " = " << std::setprecision(10) << get_best_particle()->get_personal_best_position()[j] << std::endl;
-            }
-            
-        }
-        if(fConvergence > 0.0){
-            if(mean_dist_sq < fConvergence){
-                break;
-            }
-        }
+    if (i%100 == 0){
+      std::cout << "Mean Dist Sq = " << mean_dist_sq <<std::endl;
+      std::cout << "Current LLR = " << fBestValue << std::endl;
+      std::cout << "Position = " <<std::endl;
+      for (int j = 0; j< fDim; ++j){
+        std::cout << "    Dim " << j << " = " << std::setprecision(10) << get_best_particle()->get_personal_best_position()[j] << std::endl;
+      }
     }
-    std::cout << "Finished after " << iter <<" runs out of "<< fIterations << std::endl;
-    std::cout << "Mean Dist " << mean_dist_sq <<std::endl;
-    std::cout << "Best LLR " << get_best_particle()->get_personal_best_value() << std::endl;
-
-    uncertainties = bisection(get_best_particle()->get_personal_best_position(),get_best_particle()->get_personal_best_value(),0.5,0.005);
-    std::cout << "Position for Global Minimum = "<<std::endl;
-    for (int i = 0; i< fDim; ++i){
-        std::cout << "    Dim " << i << " = " << std::setprecision(10) << get_best_particle()->get_personal_best_position()[i]  << " +" << uncertainties[i][1] << ", -" <<  uncertainties[i][0] << std::endl;
+    if(fConvergence > 0.0){
+      if(mean_dist_sq < fConvergence){
+        break;
+      }
     }
+  }
+  std::cout << "Finished after " << iter <<" runs out of "<< fIterations << std::endl;
+  std::cout << "Mean Dist " << mean_dist_sq <<std::endl;
+  std::cout << "Best LLR " << get_best_particle()->get_personal_best_value() << std::endl;
+
+  uncertainties = bisection(get_best_particle()->get_personal_best_position(),get_best_particle()->get_personal_best_value(),0.5,0.005);
+  std::cout << "Position for Global Minimum = "<<std::endl;
+  for (int i = 0; i< fDim; ++i){
+    std::cout << "    Dim " << i << " = " << std::setprecision(10) << get_best_particle()->get_personal_best_position()[i]  << " +" << uncertainties[i][1] << ", -" <<  uncertainties[i][0] << std::endl;
+  }
 }
 
 void PSO::WriteOutput(){
+  outputFile->cd();
 
-    outputFile->cd();
-
-    TVectorD* PSOParValue = new TVectorD(fDim);
-    TVectorD* PSOParError = new TVectorD(fDim);
+  TVectorD* PSOParValue = new TVectorD(fDim);
+  TVectorD* PSOParError = new TVectorD(fDim);
 
-    for(int i = 0; i < fDim; ++i)
-    {
-        (*PSOParValue)(i) = 0;
-        (*PSOParError)(i) = 0;
-    }
+  for(int i = 0; i < fDim; ++i)
+  {
+    (*PSOParValue)(i) = 0;
+    (*PSOParError)(i) = 0;
+  }
 
-    std::vector<double> minimum = get_best_particle()->get_personal_best_position();
+  std::vector<double> minimum = get_best_particle()->get_personal_best_position();
 
-    int ParCounter = 0;
+  int ParCounter = 0;
 
-    if(fTestLikelihood){
-        for(int i = 0; i < fDim; ++i){
-            (*PSOParValue)(i) = minimum[i];
-            (*PSOParError)(i) = (uncertainties[i][0]+uncertainties[i][1])/2.0;
-        }
+  if(fTestLikelihood){
+    for(int i = 0; i < fDim; ++i){
+      (*PSOParValue)(i) = minimum[i];
+      (*PSOParError)(i) = (uncertainties[i][0]+uncertainties[i][1])/2.0;
     }
-    else{
-        for (std::vector<covarianceBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
+  }
+  else{
+    for (std::vector<covarianceBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
+    {
+      if(!(*it)->IsPCA())
+      {
+        for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
         {
-            if(!(*it)->IsPCA())
-            {
-                for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
-                {
-                    double ParVal = minimum[ParCounter];
-                    //KS: Basically apply mirroring for parameters out of bounds
-                    (*PSOParValue)(ParCounter) = ParVal;
-                    (*PSOParError)(ParCounter) = (uncertainties[ParCounter][0]+uncertainties[ParCounter][1])/2.0;
-                    //KS: For fixed params HESS will not calcuate error so we need to pass prior error
-                    if((*it)->isParameterFixed(i))
-                    {
-                        (*PSOParError)(ParCounter) = (*it)->getDiagonalError(i);
-                    }
-                }
-            }
-            else
-            {
-                //KS: We need to convert parameters from PCA to normal base
-                TVectorD ParVals((*it)->GetNumParams());
-                TVectorD ParVals_PCA((*it)->getNpars());
-
-                TVectorD ErrorVals((*it)->GetNumParams());
-                TVectorD ErrorVals_PCA((*it)->getNpars());
-
-                //First save them
-                //KS: This code is super convoluted as MaCh3 can store separate matrices while PSO has one matrix. In future this will be simplified, keep it like this for now.
-                const int StartVal = ParCounter;
-                for(int i = 0; i < (*it)->getNpars(); ++i, ++ParCounter)
-                {
-                    ParVals_PCA(i) = minimum[ParCounter];
-                    ErrorVals_PCA(i) = (uncertainties[ParCounter][0]+uncertainties[ParCounter][1])/2.0;
-                }
-                ParVals = ((*it)->getTransferMatrix())*ParVals_PCA;
-                ErrorVals = ((*it)->getTransferMatrix())*ErrorVals_PCA;
-
-                ParCounter = StartVal;
-                //KS: Now after going from PCA to normal let';s save it
-                for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
-                {
-                    (*PSOParValue)(ParCounter) = ParVals(i);
-                    (*PSOParError)(ParCounter) = std::fabs(ErrorVals(i));
-                    //int ParCounterMatrix = StartVal;
-                    //If fixed take prior
-                    if((*it)->isParameterFixedPCA(i))
-                    {
-                        (*PSOParError)(ParCounter) = (*it)->getDiagonalError(i);
-                    }
-                }
-            }
+          double ParVal = minimum[ParCounter];
+          //KS: Basically apply mirroring for parameters out of bounds
+          (*PSOParValue)(ParCounter) = ParVal;
+          (*PSOParError)(ParCounter) = (uncertainties[ParCounter][0]+uncertainties[ParCounter][1])/2.0;
+          //KS: For fixed params HESS will not calcuate error so we need to pass prior error
+          if((*it)->isParameterFixed(i))
+          {
+            (*PSOParError)(ParCounter) = (*it)->getDiagonalError(i);
+          }
         }
+      }
+      else
+      {
+        //KS: We need to convert parameters from PCA to normal base
+        TVectorD ParVals((*it)->GetNumParams());
+        TVectorD ParVals_PCA((*it)->getNpars());
+
+        TVectorD ErrorVals((*it)->GetNumParams());
+        TVectorD ErrorVals_PCA((*it)->getNpars());
+
+        //First save them
+        //KS: This code is super convoluted as MaCh3 can store separate matrices while PSO has one matrix. In future this will be simplified, keep it like this for now.
+        const int StartVal = ParCounter;
+        for(int i = 0; i < (*it)->getNpars(); ++i, ++ParCounter)
+        {
+          ParVals_PCA(i) = minimum[ParCounter];
+          ErrorVals_PCA(i) = (uncertainties[ParCounter][0]+uncertainties[ParCounter][1])/2.0;
+        }
+        ParVals = ((*it)->getTransferMatrix())*ParVals_PCA;
+        ErrorVals = ((*it)->getTransferMatrix())*ErrorVals_PCA;
+
+        ParCounter = StartVal;
+        //KS: Now after going from PCA to normal let';s save it
+        for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
+        {
+          (*PSOParValue)(ParCounter) = ParVals(i);
+          (*PSOParError)(ParCounter) = std::fabs(ErrorVals(i));
+          //int ParCounterMatrix = StartVal;
+          //If fixed take prior
+          if((*it)->isParameterFixedPCA(i))
+          {
+            (*PSOParError)(ParCounter) = (*it)->getDiagonalError(i);
+          }
+        }
+      }
     }
+  }
 
-    PSOParValue->Write("PSOParValue");
-    PSOParError->Write("PSOParError");
-    delete PSOParValue;
-    delete PSOParError;
-    // Save all the output
-    SaveOutput();
+  PSOParValue->Write("PSOParValue");
+  PSOParError->Write("PSOParError");
+  delete PSOParValue;
+  delete PSOParError;
+  // Save all the output
+  SaveOutput();
 }
 
 
diff --git a/mcmc/PSO.h b/mcmc/PSO.h
index 472cd5bf..3596a45a 100644
--- a/mcmc/PSO.h
+++ b/mcmc/PSO.h
@@ -16,144 +16,144 @@
 /// @brief Class particle - stores the position, velocity and personal best
 /// With functions which move particle and update velocity
 class particle{
-    public:
-        particle(){};
-        particle(std::vector<double> pos, std::vector<double> vel) : position(pos), velocity(vel){};
-        /// @brief Destructor
-        virtual ~particle() {};
-
-        void set_position(std::vector<double> new_position) {
-            position = new_position;
-        };
-
-        std::vector<double> get_position() {
-            return position;
-        };
-
-        void set_personal_best_position(std::vector<double> new_pos){
-            personal_best_position = new_pos;
-        };
-
-        std::vector<double> get_personal_best_position(){
-            return personal_best_position;
-        };
-
-        void set_personal_best_value(double new_val){
-            personal_best_value = new_val;
-        };
-
-        double get_personal_best_value(){
-            return personal_best_value;
-        };
-
-        std::vector<double> get_velocity(){
-            return velocity;
-        };
-
-        void set_velocity(std::vector<double> new_velocity){
-            velocity = new_velocity;
-        };
-
-        double get_value(){
-            return curr_value;
-        };
-        void set_value(double val){
-            curr_value = val;
-        };
-
-    private:
-        std::vector<double> position;
-        std::vector<double> velocity;
-        double personal_best_value;
-        double curr_value;
-        std::vector<double> personal_best_position;
+  public:
+      particle(){};
+      particle(std::vector<double> pos, std::vector<double> vel) : position(pos), velocity(vel){};
+      /// @brief Destructor
+      virtual ~particle() {};
+
+      void set_position(std::vector<double> new_position) {
+        position = new_position;
+      };
+
+      std::vector<double> get_position() {
+        return position;
+      };
+
+      void set_personal_best_position(std::vector<double> new_pos){
+        personal_best_position = new_pos;
+      };
+
+      std::vector<double> get_personal_best_position(){
+        return personal_best_position;
+      };
+
+      void set_personal_best_value(double new_val){
+        personal_best_value = new_val;
+      };
+
+      double get_personal_best_value(){
+        return personal_best_value;
+      };
+
+      std::vector<double> get_velocity(){
+        return velocity;
+      };
+
+      void set_velocity(std::vector<double> new_velocity){
+        velocity = new_velocity;
+      };
+
+      double get_value(){
+        return curr_value;
+      };
+      void set_value(double val){
+        curr_value = val;
+      };
+
+  private:
+    std::vector<double> position;
+    std::vector<double> velocity;
+    double personal_best_value;
+    double curr_value;
+    std::vector<double> personal_best_position;
 };
 
  /// @brief Class PSO, consist of a vector of object Class Particle and global best
  /// Takes in the size (number of particle) and number of iteration
  /// functions includes: finding global best, updating velocity, actual minimisation function
 class PSO : public LikelihoodFit {
-    public:
-        /// @brief constructor
-        PSO(manager * const fitMan);
-        /// @brief Destructor
-        virtual ~PSO() {};
-
-        particle* get_best_particle(){
-            return best_particle;
-        }
-        void set_best_particle(particle* n){
-            best_particle = n;
-        }
-
-        std::vector<std::vector<double>> bisection(std::vector<double>position,double minimum, double range, double precision);
-        std::vector<std::vector<double>> calc_uncertainty(std::vector<double>position,double minimum);
-        void init();
-        void uncertainty_check(std::vector<double> previous_pos);
-        void run();
-        void WriteOutput();
-        void runMCMC() override;
-        double CalcChi2(const double* x);
-        double rastriginFunc(const double* x);
-        double swarmIterate();
-
-        std::vector<double> vector_multiply(std::vector<double> velocity, double mul){
-            transform(velocity.begin(),velocity.end(),velocity.begin(),std::bind1st(std::multiplies<double>(),mul));
-            return velocity;
-        };
-
-        std::vector<double> vector_add(std::vector<double> v1, std::vector<double> v2){
-            std::vector<double> v3;
-            transform(v1.begin(), v1.end(), v2.begin(), back_inserter(v3), std::plus<double>());
-            return v3;
-        };
-        std::vector<double> vector_subtract(std::vector<double> v1, std::vector<double> v2){
-            std::vector<double> v3 ;
-            transform(v1.begin(), v1.end(), v2.begin(), back_inserter(v3), std::minus<double>());
-            return v3;
-        };
-        std::vector<double> three_vector_addition(std::vector<double> vec1, std::vector<double> vec2,std::vector<double> vec3){
-            for (size_t i = 0; i < vec1.size(); ++i) {
-                vec1[i] += vec2[i] + vec3[i];
-            }
-            return vec1;
-        };
-        std::vector<double> four_vector_addition(std::vector<double> vec1, std::vector<double> vec2,std::vector<double> vec3,std::vector<double> vec4){
-            for (size_t i = 0; i < vec1.size(); ++i) {
-                vec1[i] += vec2[i] + vec3[i] + vec4[i];
-            }
-            return vec1;
-        };
-
-        double CalcChi(std::vector<double> x){
-            double* a = &x[0];
-            return CalcChi2(a);
-        };
-
-    inline std::string GetName()const {return "PSO";};
-
-    private:
-        particle* best_particle;
-        double fBestValue;
-        std::vector<double> prior;
-        std::vector<bool> fixed;
-        std::vector<double> ranges_max;
-        std::vector<double> ranges_min;
-        std::vector<particle*> system;
-        double fInertia;
-        double fOne;
-        double fTwo;
-        double fConvergence;
-        int fIterations;
-        double fConstriction;
-        std::vector<std::vector<double> > uncertainties;
-
-        int fParticles;
-        static const int kMaxParticles = 10000;
-        std::vector<double*> paramlist;
-        double vel[kMaxParticles];
-        double* par;
-
-        int fDim;
+  public:
+    /// @brief constructor
+    PSO(manager * const fitMan);
+    /// @brief Destructor
+    virtual ~PSO() {};
+
+    particle* get_best_particle(){
+        return best_particle;
+    }
+    void set_best_particle(particle* n){
+      best_particle = n;
+    }
+
+    std::vector<std::vector<double>> bisection(std::vector<double>position,double minimum, double range, double precision);
+    std::vector<std::vector<double>> calc_uncertainty(std::vector<double>position,double minimum);
+    void init();
+    void uncertainty_check(std::vector<double> previous_pos);
+    void run();
+    void WriteOutput();
+    void runMCMC() override;
+    double CalcChi2(const double* x);
+    double rastriginFunc(const double* x);
+    double swarmIterate();
+
+    std::vector<double> vector_multiply(std::vector<double> velocity, double mul){
+      transform(velocity.begin(),velocity.end(),velocity.begin(),std::bind1st(std::multiplies<double>(),mul));
+      return velocity;
+    };
+
+    std::vector<double> vector_add(std::vector<double> v1, std::vector<double> v2){
+      std::vector<double> v3;
+      transform(v1.begin(), v1.end(), v2.begin(), back_inserter(v3), std::plus<double>());
+      return v3;
+    };
+    std::vector<double> vector_subtract(std::vector<double> v1, std::vector<double> v2){
+      std::vector<double> v3 ;
+      transform(v1.begin(), v1.end(), v2.begin(), back_inserter(v3), std::minus<double>());
+      return v3;
+    };
+    std::vector<double> three_vector_addition(std::vector<double> vec1, std::vector<double> vec2,std::vector<double> vec3){
+      for (size_t i = 0; i < vec1.size(); ++i) {
+        vec1[i] += vec2[i] + vec3[i];
+      }
+      return vec1;
+    };
+    std::vector<double> four_vector_addition(std::vector<double> vec1, std::vector<double> vec2,std::vector<double> vec3,std::vector<double> vec4){
+      for (size_t i = 0; i < vec1.size(); ++i) {
+        vec1[i] += vec2[i] + vec3[i] + vec4[i];
+      }
+      return vec1;
+    };
+
+    double CalcChi(std::vector<double> x){
+      double* a = &x[0];
+      return CalcChi2(a);
+    };
+
+  inline std::string GetName()const {return "PSO";};
+
+  private:
+    particle* best_particle;
+    double fBestValue;
+    std::vector<double> prior;
+    std::vector<bool> fixed;
+    std::vector<double> ranges_max;
+    std::vector<double> ranges_min;
+    std::vector<particle*> system;
+    double fInertia;
+    double fOne;
+    double fTwo;
+    double fConvergence;
+    int fIterations;
+    double fConstriction;
+    std::vector<std::vector<double> > uncertainties;
+
+    int fParticles;
+    static const int kMaxParticles = 10000;
+    std::vector<double*> paramlist;
+    double vel[kMaxParticles];
+    double* par;
+
+    int fDim;
 };
 
diff --git a/mcmc/StatisticalUtils.h b/mcmc/StatisticalUtils.h
index d0604773..4680ef0d 100644
--- a/mcmc/StatisticalUtils.h
+++ b/mcmc/StatisticalUtils.h
@@ -491,3 +491,44 @@ inline double GetIQR(TH1D *Hist) {
 
   return quartiles[2] - quartiles[0];
 }
+
+// ********************
+/// @brief Compute the Kullback-Leibler divergence between two TH2Poly histograms.
+///
+/// @param DataPoly Pointer to the data histogram (TH2Poly).
+/// @param PolyMC Pointer to the Monte Carlo histogram (TH2Poly).
+/// @return The Kullback-Leibler divergence value. Returns 0 if the data or MC integral is zero.
+inline double ComputeKLDivergence(TH2Poly* DataPoly, TH2Poly* PolyMC) {
+// *********************
+  double klDivergence = 0.0;
+  double DataIntegral = NoOverflowIntegral(DataPoly);
+  double MCIntegral = NoOverflowIntegral(PolyMC);
+  for (int i = 1; i < DataPoly->GetNumberOfBins()+1; ++i)
+  {
+    if (DataPoly->GetBinContent(i) > 0 && PolyMC->GetBinContent(i) > 0) {
+      klDivergence += DataPoly->GetBinContent(i) / DataIntegral *
+      std::log((DataPoly->GetBinContent(i) / DataIntegral) / ( PolyMC->GetBinContent(i) / MCIntegral));
+    }
+  }
+  return klDivergence;
+}
+// ********************
+/// @brief KS: Combine p-values using Fisher's method.
+///
+/// @param pvalues A vector of individual p-values to combine.
+/// @return The combined p-value, representing the overall significance.
+inline double FisherCombinedPValue(const std::vector<double>& pvalues) {
+// ********************
+
+  double testStatistic = 0;
+  for(size_t i = 0; i < pvalues.size(); i++)
+  {
+    const double pval = std::max(0.00001, pvalues[i]);
+    testStatistic += -2.0 * std::log(pval);
+  }
+  // Degrees of freedom is twice the number of p-values
+  int degreesOfFreedom = 2 * pvalues.size();
+  double pValue = TMath::Prob(testStatistic, degreesOfFreedom);
+
+  return pValue;
+}
diff --git a/plotting/CMakeLists.txt b/plotting/CMakeLists.txt
index 285d439e..2323f166 100644
--- a/plotting/CMakeLists.txt
+++ b/plotting/CMakeLists.txt
@@ -4,6 +4,7 @@ add_subdirectory(plottingUtils)
 foreach(app
     GetPostfitParamPlots
     PlotLLH
+    MatrixPlotter
     )
     add_executable( ${app} ${app}.cpp )
     target_link_libraries( ${app} MaCh3::Plotting )
diff --git a/plotting/MatrixPlotter.cpp b/plotting/MatrixPlotter.cpp
new file mode 100755
index 00000000..9556feef
--- /dev/null
+++ b/plotting/MatrixPlotter.cpp
@@ -0,0 +1,255 @@
+#include <iostream>
+#include <algorithm>
+
+#include "TFile.h"
+#include "TH2D.h"
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TError.h"
+#include "TPaletteAxis.h"
+#include "TColor.h"
+
+#include "plottingUtils/plottingUtils.h"
+#include "plottingUtils/plottingManager.h"
+
+TH2D* GetSubMatrix(TH2D *MatrixFull, const std::string& Title, const std::vector<std::string>& Params)
+{
+  std::vector<int> ParamIndex(Params.size(), -999);
+
+  for(size_t i = 0; i < Params.size(); i++)
+  {
+    for(int j = 0; j < MatrixFull->GetNbinsX(); j++)
+    {
+      if(MatrixFull->GetXaxis()->GetBinLabel(j+1) == Params[i])
+      {
+        ParamIndex[i] = j;
+        break;
+      }
+    }
+  }
+
+  for(size_t i = 0; i < Params.size(); i++)
+  {
+    if(ParamIndex[i] == -999 )
+      MACH3LOG_ERROR("Didn't find param {} in matrix within {} sub-block", Params[i], Title);
+  }
+
+  auto new_end = std::remove(ParamIndex.begin(), ParamIndex.end(), -999);
+  ParamIndex.erase(new_end, ParamIndex.end());
+
+  TH2D* Hist = new TH2D(Title.c_str(), Title.c_str(), ParamIndex.size(), 0, ParamIndex.size(), ParamIndex.size(), 0, ParamIndex.size());
+  Hist->GetZaxis()->SetTitle("Correlation");
+  Hist->SetMinimum(-1);
+  Hist->SetMaximum(1);
+  Hist->GetXaxis()->SetLabelSize(0.015);
+  Hist->GetYaxis()->SetLabelSize(0.015);
+
+  for(size_t x = 0; x < ParamIndex.size(); x++)
+  {
+    for(size_t y = 0; y < ParamIndex.size(); y++)
+    {
+      Hist->SetBinContent(x+1, y+1, MatrixFull->GetBinContent(ParamIndex[x]+1, ParamIndex[y]+1));
+    }
+    Hist->GetXaxis()->SetBinLabel(x+1, MatrixFull->GetXaxis()->GetBinLabel(ParamIndex[x]+1));
+    Hist->GetYaxis()->SetBinLabel(x+1, MatrixFull->GetXaxis()->GetBinLabel(ParamIndex[x]+1));
+  }
+  return Hist;
+}
+
+void SetupInfo(const std::string& Config, std::vector<std::string>& Title, std::vector<std::vector<std::string>>& Params)
+{
+  // Load the YAML file
+  YAML::Node config = YAML::LoadFile(Config);
+
+  // Access the "MatrixPlotter" section
+  YAML::Node settings = config["MatrixPlotter"];
+
+  // Retrieve the list of Titles
+  Title = settings["Titles"].as<std::vector<std::string>>();
+  Params.resize(Title.size());
+
+  // Retrieve parameters for each title (using the titles as keys)
+  for(size_t it = 0; it < Title.size(); it++)
+  {
+    if (settings[Title[it]]) {
+      Params[it] = settings[Title[it]].as<std::vector<std::string>>();
+    } else {
+      MACH3LOG_ERROR("Missing key in YAML for: {}", Title[it]);
+      throw MaCh3Exception(__FILE__ , __LINE__ );
+    }
+  }
+
+  // Check if sizes match
+  if(Title.size() != Params.size())
+  {
+    MACH3LOG_ERROR("Size doesn't match");
+    throw MaCh3Exception(__FILE__ , __LINE__ );
+  }
+}
+
+void PlotMatrix(std::string Config, std::string File)
+{
+  // Open the ROOT file
+  TFile *file = TFile::Open(File.c_str());
+  if (!file || file->IsZombie()) {
+    MACH3LOG_ERROR("Error opening file");
+    throw MaCh3Exception(__FILE__ , __LINE__ );
+  }
+
+  TH2D *MatrixFull = nullptr;
+  file->GetObject("Correlation_plot", MatrixFull);
+
+  if (!MatrixFull) {
+    MACH3LOG_ERROR("Error: Could not retrieve histogram");
+    throw MaCh3Exception(__FILE__ , __LINE__ );
+  }
+
+  auto MatrixPlot = std::make_unique<TCanvas>("MatrixPlot", "MatrixPlot", 0, 0, 1024, 1024);
+  MatrixPlot->SetGrid();
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+  MatrixPlot->SetTickx();
+  MatrixPlot->SetTicky();
+  MatrixPlot->SetBottomMargin(0.2);
+  MatrixPlot->SetTopMargin(0.1);
+  MatrixPlot->SetRightMargin(0.15);
+  MatrixPlot->SetLeftMargin(0.15);
+
+  gStyle->SetOptTitle(1);
+  gStyle->SetPaintTextFormat("4.1f");
+
+  // Make pretty Correlation colors (red to blue)
+  const int NRGBs = 5;
+  TColor::InitializeColors();
+  Double_t stops[NRGBs] = { 0.00, 0.25, 0.50, 0.75, 1.00 };
+  Double_t red[NRGBs]   = { 0.00, 0.25, 1.00, 1.00, 0.50 };
+  Double_t green[NRGBs] = { 0.00, 0.25, 1.00, 0.25, 0.00 };
+  Double_t blue[NRGBs]  = { 0.50, 1.00, 1.00, 0.25, 0.00 };
+  TColor::CreateGradientColorTable(5, stops, red, green, blue, 255);
+  gStyle->SetNumberContours(255);
+
+  //To avoid TCanvas::Print> messages
+  gErrorIgnoreLevel = kWarning;
+
+  MatrixPlot->Print("MatrixPlot.pdf[");
+
+  MatrixFull->Draw("COLZ");
+  MatrixPlot->Print("MatrixPlot.pdf");
+  std::vector<std::string> Title;
+  std::vector<std::vector<std::string>> Params;
+
+  SetupInfo(Config, Title, Params);
+
+  for(size_t it = 0; it < Title.size(); it++)
+  {
+    TH2D* Hist = GetSubMatrix(MatrixFull, Title[it], Params[it]);
+    Hist->GetXaxis()->LabelsOption("v");
+    Hist->Draw("COLZ TEXT");
+    MatrixPlot->Print("MatrixPlot.pdf");
+    delete Hist;
+  }
+  MatrixPlot->Print("MatrixPlot.pdf]");
+
+  delete MatrixFull;
+  file->Close();
+  delete file;
+}
+
+void CompareMatrices(std::string Config, std::string File1, std::string Title1, std::string File2, std::string Title2)
+{
+  // Open the ROOT file
+  const int NFiles = 2;
+  TFile *file[NFiles];
+  file[0] = TFile::Open(File1.c_str());
+  file[1] = TFile::Open(File2.c_str());
+
+  TH2D *MatrixFull[2] = {nullptr};
+  for(int i = 0; i < NFiles; i++) {
+    file[i]->GetObject("Correlation_plot", MatrixFull[i]);
+  }
+  auto MatrixPlot = std::make_unique<TCanvas>("MatrixPlot", "MatrixPlot", 0, 0, 1024, 1024);
+  MatrixPlot->SetGrid();
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+  MatrixPlot->SetTickx();
+  MatrixPlot->SetTicky();
+  MatrixPlot->SetBottomMargin(0.2);
+  MatrixPlot->SetTopMargin(0.1);
+  MatrixPlot->SetRightMargin(0.15);
+  MatrixPlot->SetLeftMargin(0.15);
+
+  gStyle->SetOptTitle(1);
+
+  //KS: Fancy colots
+  const int NRGBs = 10;
+  TColor::InitializeColors();
+  Double_t stops[NRGBs] = { 0.00, 0.10, 0.25, 0.35, 0.50, 0.60, 0.65, 0.75, 0.90, 1.00 };
+  Double_t red[NRGBs]   = { 0.50, 1.00, 1.00, 0.25, 0.00, 0.10, 0.50, 1.00, 0.75, 0.55 };
+  Double_t green[NRGBs] = { 0.00, 0.25, 1.00, 0.25, 0.00, 0.60, 0.90, 1.00, 0.75, 0.75 };
+  Double_t blue[NRGBs]  = { 0.00, 0.25, 1.00, 1.00, 0.50, 0.60, 0.90, 1.00, 0.05, 0.05 };
+  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, 255);
+  gStyle->SetNumberContours(255);
+
+  //To avoid TCanvas::Print> messages
+  gErrorIgnoreLevel = kWarning;
+
+  MatrixPlot->Print("MatrixComparePlot.pdf[");
+
+  std::vector<std::string> Title;
+  std::vector<std::vector<std::string>> Params;
+
+  SetupInfo(Config, Title, Params);
+
+  for(size_t it = 0; it < Title.size(); it++)
+  {
+    TH2D* Hist[2] = {nullptr};
+    for(int i = 0; i < NFiles; i++)
+    {
+      Hist[i] = GetSubMatrix(MatrixFull[i], Title[it], Params[it]);
+    }
+    Hist[0]->GetZaxis()->SetTitle( (Title1 + "/" +  Title2).c_str());
+    Hist[0]->GetXaxis()->LabelsOption("v");
+
+    Hist[0]->Draw("COLZ");
+    MatrixPlot->Print("MatrixComparePlot.pdf");
+
+    for(int i = 0; i < NFiles; i++)
+    {
+      delete Hist[i];
+    }
+  }
+  MatrixPlot->Print("MatrixComparePlot.pdf]");
+
+  for(int i = 0; i < NFiles; i++)
+  {
+    delete MatrixFull[i];
+    file[i]->Close();
+    delete file[i];
+  }
+}
+
+int main(int argc, char *argv[]) 
+{
+  SetMaCh3LoggerFormat();
+
+  if (argc != 3 && argc != 6)
+  {
+    MACH3LOG_INFO("How to use: {} MCMC_Processor_Output.root config", argv[0]);
+    throw MaCh3Exception(__FILE__ , __LINE__ );
+  }
+
+  if (argc == 3)
+  {
+    PlotMatrix(std::string(argv[1]), std::string(argv[2]));
+  }
+
+  if (argc == 6)
+  {
+    MACH3LOG_INFO("Comparing matrices");
+    CompareMatrices(std::string(argv[1]), std::string(argv[2]), std::string(argv[3]), std::string(argv[4]), std::string(argv[5]));
+  }
+
+  MACH3LOG_INFO("Finished plotting matrices");
+  return 0;
+}
+
diff --git a/plotting/README.md b/plotting/README.md
index 6cef87b8..e819506f 100644
--- a/plotting/README.md
+++ b/plotting/README.md
@@ -86,3 +86,6 @@ options:
     -o the name of the output pdf
 
     -d a string specifying additional drawing options to pass to the histogram draw calls, e.g. `-d "C"` will plot smooth curves through the histogram bins. See https://root.cern/doc/master/classTHistPainter.html#HP01a for possible options.
+
+
+**MatrixPlotter** - As input you need output from `ProcessMCMC`, keep in mind you need to run it with `PlotCorr`. The executable allows to plot submatrices and whatever combination of parameters you like.
diff --git a/samplePDF/Structs.cpp b/samplePDF/Structs.cpp
index b2015586..3253ac5d 100644
--- a/samplePDF/Structs.cpp
+++ b/samplePDF/Structs.cpp
@@ -391,7 +391,6 @@ TH2Poly* PolyScaleWidth(TH2Poly *Histogram, double scale) {
     area = (xup-xlow)*(yup-ylow);
     HistCopy->SetBinContent(i, Histogram->GetBinContent(i)/(area*scale));
   }
-
   return HistCopy;
 }
 
@@ -412,7 +411,6 @@ double PolyIntegralWidth(TH2Poly *Histogram) {
     area = (xup-xlow)*(yup-ylow);
     integral += Histogram->GetBinContent(i)*area;
   }
-
   return integral;
 }
 
@@ -428,6 +426,37 @@ void RemoveFitter(TH1D* hist, const std::string& name) {
   delete fitter;
 }
 
+// *************************
+TGraphAsymmErrors* MakeAsymGraph(TH1D* sigmaArrayLeft, TH1D* sigmaArrayCentr, TH1D* sigmaArrayRight, const std::string& title) {
+// *************************
+  TGraphAsymmErrors *var = new TGraphAsymmErrors(sigmaArrayCentr);
+  var->SetNameTitle((title).c_str(), (title).c_str());
+
+  // Need to draw TGraphs to set axes labels
+  var->Draw("AP");
+  var->GetXaxis()->SetTitle(sigmaArrayCentr->GetXaxis()->GetTitle());
+  var->GetYaxis()->SetTitle("Number of events/bin");
+
+  for (int m = 0; m < var->GetN(); ++m)
+  {
+    double xlow = sigmaArrayLeft->GetBinContent(m+1);
+    double xhigh = sigmaArrayRight->GetBinContent(m+1);
+    double xtemp;
+
+    // Figure out which variation is larger so we set the error correctly
+    if (xlow > xhigh)
+    {
+      xtemp = xlow;
+      xlow = xhigh;
+      xhigh = xtemp;
+    }
+
+    var->SetPointEYhigh(m, xhigh - var->GetY()[m]);
+    var->SetPointEYlow(m, var->GetY()[m] - xlow);
+  }
+  return var;
+}
+
 // ****************
 //DB Get the Cernekov momentum threshold in MeV
 double returnCherenkovThresholdMomentum(int PDG) {
@@ -442,7 +471,6 @@ double returnCherenkovThresholdMomentum(int PDG) {
 // Recalculate Q^2 after Eb shift. Takes in shifted lepton momentum, lepton angle, and true neutrino energy
 double CalculateQ2(double PLep, double PUpd, double EnuTrue, double InitialQ2){
 // ***************************************************************************
-
   const double MLep = 0.10565837;
 
   // Caluclate muon energy
@@ -458,12 +486,10 @@ double CalculateQ2(double PLep, double PUpd, double EnuTrue, double InitialQ2){
   return Q2Upd - InitialQ2;
 }
 
-
 // **************************************************************************
 // Recalculate Enu after Eb shift. Takes in shifted lepton momentum, lepton angle, and binding energy change, and if nu/anu
 double CalculateEnu(double PLep, double costh, double Eb, bool neutrino){
 // ***************************************************************************
-
   double mNeff = 0.93956536 - Eb / 1000.;
   double mNoth = 0.93827203;
 
@@ -478,5 +504,4 @@ double CalculateEnu(double PLep, double costh, double Eb, bool neutrino){
   double Enu = (2 * mNeff * eLep - mLep * mLep + mNoth * mNoth - mNeff * mNeff) /(2 * (mNeff - eLep + PLep * costh));
 
   return Enu;
-
 }
diff --git a/samplePDF/Structs.h b/samplePDF/Structs.h
index d4e075eb..02efe365 100644
--- a/samplePDF/Structs.h
+++ b/samplePDF/Structs.h
@@ -57,6 +57,7 @@
 #include "TObjString.h"
 #include "TH2Poly.h"
 #include "TFile.h"
+#include "TGraphAsymmErrors.h"
 
 #ifdef MULTITHREAD
 #include "omp.h"
@@ -141,10 +142,9 @@ inline std::string GetTF1(const SplineInterpolation i) {
       Func = "([1]+[0]*x)";
       break;
     default:
-      std::cerr << "UNKNOWN SPECIFIED!" << std::endl;
-      std::cerr << "You gave  " << static_cast<int>(i) << std::endl;
-      std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
-      throw;
+      MACH3LOG_ERROR("UNKNOWN SPLINE INTERPOLATION SPECIFIED!");
+      MACH3LOG_ERROR("You gave {}", static_cast<int>(i));
+      throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   return Func;
 }
@@ -174,10 +174,9 @@ inline RespFuncType SplineInterpolation_ToRespFuncType(const SplineInterpolation
       Type = RespFuncType::kTF1_red;
       break;
     default:
-      std::cerr << "UNKNOWN SPLINE INTERPOLATION SPECIFIED!" << std::endl;
-      std::cerr << "You gave  " << static_cast<int>(i) << std::endl;
-      std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
-      throw;
+      MACH3LOG_ERROR("UNKNOWN SPLINE INTERPOLATION SPECIFIED!");
+      MACH3LOG_ERROR("You gave {}", static_cast<int>(i));
+      throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   return Type;
 }
@@ -206,10 +205,9 @@ inline std::string SplineInterpolation_ToString(const SplineInterpolation i) {
       name = "LinearFunc";
       break;
     default:
-      std::cerr << "UNKNOWN SPLINE INTERPOLATION SPECIFIED!" << std::endl;
-      std::cerr << "You gave  " << static_cast<int>(i) << std::endl;
-      std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
-      throw;
+      MACH3LOG_ERROR("UNKNOWN SPLINE INTERPOLATION SPECIFIED!");
+      MACH3LOG_ERROR("You gave {}", static_cast<int>(i));
+      throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   return name;
 }
@@ -255,10 +253,9 @@ inline std::string SystType_ToString(const SystType i) {
       name = "Functional";
       break;
     default:
-      std::cerr << "UNKNOWN SYST TYPE SPECIFIED!" << std::endl;
-      std::cerr << "You gave  " << static_cast<int>(i) << std::endl;
-      std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
-      throw;
+      MACH3LOG_ERROR("UNKNOWN SYST TYPE SPECIFIED!");
+      MACH3LOG_ERROR("You gave {}", static_cast<int>(i));
+      throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   return name;
 }
@@ -450,10 +447,9 @@ inline std::string TestStatistic_ToString(TestStatistic i) {
     name = "DembinskiAbdelmottele";
     break;
     default:
-      std::cerr << "UNKNOWN LIKELHOOD SPECIFIED!" << std::endl;
-      std::cerr << "You gave test-statistic " << i << std::endl;
-      std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
-      throw;
+      MACH3LOG_ERROR("UNKNOWN LIKELIHOOD SPECIFIED!");
+      MACH3LOG_ERROR("You gave test-statistic {}", static_cast<int>(i));
+      throw MaCh3Exception(__FILE__ , __LINE__ );
   }
   return name;
 }
@@ -490,20 +486,28 @@ void CheckTH2PolyFileVersion(TFile *file);
 /// @brief KS: Remove fitted TF1 from hist to make comparison easier
 void RemoveFitter(TH1D* hist, const std::string& name);
 
+/// @brief Used by sigma variation, check how 1 sigma changes spectra
+/// @param sigmaArrayLeft sigma var hist at -1 or -3 sigma shift
+/// @param sigmaArrayCentr sigma var hist at prior values
+/// @param sigmaArrayRight sigma var hist at +1 or +3 sigma shift
+/// @param title A tittle for returned object
+/// @return A `TGraphAsymmErrors` object that visualizes the sigma variation of spectra, showing confidence intervals between different sigma shifts.
+TGraphAsymmErrors* MakeAsymGraph(TH1D* sigmaArrayLeft, TH1D* sigmaArrayCentr, TH1D* sigmaArrayRight, const std::string& title);
+
 /// @brief Helper to check if files exist or not
 inline std::string file_exists(std::string filename) {
   std::ifstream infile(filename.c_str());
   if (!infile.good()) {
-    std::cerr << "*** ERROR ***" << std::endl;
-    std::cerr << "File " << filename << " does not exist" << std::endl;
-    std::cerr << "Please try again" << std::endl;
-    std::cerr << "*************" << std::endl;
-    throw;
+    MACH3LOG_ERROR("*** ERROR ***");
+    MACH3LOG_ERROR("File {} does not exist", filename);
+    MACH3LOG_ERROR("Please try again");
+    MACH3LOG_ERROR("*************");
+    throw MaCh3Exception(__FILE__ , __LINE__ );
   }
-
   return filename;
 }
-/// @brief DB Get the Cernekov momentum threshold in MeV
+
+/// @brief DB Get the Cherenkov momentum threshold in MeV
 double returnCherenkovThresholdMomentum(int PDG);
 
 double CalculateQ2(double PLep, double PUpd, double EnuTrue, double InitialQ2 = 0.0);
diff --git a/samplePDF/samplePDFBase.cpp b/samplePDF/samplePDFBase.cpp
index 1a755ba3..d6a2fc2b 100644
--- a/samplePDF/samplePDFBase.cpp
+++ b/samplePDF/samplePDFBase.cpp
@@ -217,7 +217,7 @@ double samplePDFBase::getTestStatLLH(const double data, const double mc, const d
   // Need some MC
   if (mc == 0) return 0.0;
 
-  // The MC used in the likeliihood calculation
+  // The MC used in the likelihood calculation
   // Is allowed to be changed by Barlow Beeston beta parameters
   double newmc = mc;
 
@@ -226,7 +226,7 @@ double samplePDFBase::getTestStatLLH(const double data, const double mc, const d
     //CW: Not full Barlow-Beeston or what is referred to as "light": we're not introducing any more parameters
     // Assume the MC has a Gaussian distribution around generated
     // As in https://arxiv.org/abs/1103.0354 eq 10, 11
-    //CW: Calculate the Barlow-Beeston likelhood contribution from MC statistics
+    //CW: Calculate the Barlow-Beeston likelihood contribution from MC statistics
     // Assumes the beta scaling parameters are Gaussian distributed
     // Follows arXiv:1103.0354 section 5 and equation 8, 9, 10, 11 on page 4/5
     // Essentially solves equation 11
@@ -241,9 +241,8 @@ double samplePDFBase::getTestStatLLH(const double data, const double mc, const d
       // 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__ );
       }
       // Solve for the positive beta
       const double beta = (-1*temp+sqrt(temp2))/2.;
@@ -312,7 +311,7 @@ double samplePDFBase::getTestStatLLH(const double data, const double mc, const d
 
         return stat;
       }
-      // Auxillary variables
+      // Auxiliary variables
       const long double b = mc/w2;
       const long double a = mc*b+1;
       const long double k = data;
@@ -335,16 +334,15 @@ double samplePDFBase::getTestStatLLH(const double data, const double mc, const d
     break;
     case (kPoisson):
     {
-	  //Just call getTestStatLLH which doesn't take in weights
-	  //and is a poisson likelihood comparison.
+      //Just call getTestStatLLH which doesn't take in weights
+      //and is a Poisson likelihood comparison.
       return getTestStatLLH(data, mc);//stat;
     }
     break;
 
     default:
     std::cerr << "Couldn't find TestStatistic " << fTestStatistic << " exiting!" << std::endl;
-    std::cerr << __FILE__ << ":" << __LINE__ << std::endl;
-    throw;
+    throw MaCh3Exception(__FILE__ , __LINE__ );
   } // end switch
 }