ProcessMCMC.cpp File Reference

Main exectable responsible for different types of MCMC processing like drawing posteriors, triangle plots etc. Actual implantation of methods is in MCMCProcessor. More...

#include "Fitters/OscProcessor.h"
#include "Manager/Manager.h"

Include dependency graph for ProcessMCMC.cpp:

Go to the source code of this file.

Functions
void	ProcessMCMC (const std::string &inputFile)
	Main function processing MCMC and Producing plots. More...
void	MultipleProcessMCMC ()
	Function producing comparison of posterior and more betwen a few MCMC chains. More...
void	CalcBayesFactor (MCMCProcessor *Processor)
	KS: Calculate Bayes factor for a given hypothesis, most informative are those related to osc params. However, it make relative easy interpretation for switch dials. More...
void	CalcSavageDickey (MCMCProcessor *Processor)
void	CalcBipolarPlot (MCMCProcessor *Processor)
void	CalcParameterEvolution (MCMCProcessor *Processor)
void	GetTrianglePlot (MCMCProcessor *Processor)
void	DiagnoseCovarianceMatrix (MCMCProcessor *Processor, const std::string &inputFile)
	KS: You validate stability of posterior covariance matrix, you set burn calc cov matrix increase burn calc again and compare. By performing such operation several hundred times we can check when matrix becomes stable. More...
void	ReweightPrior (MCMCProcessor *Processor)
TH2D *	TMatrixIntoTH2D (TMatrixDSym *Matrix, const std::string &title)
	KS: Convert TMatrix to TH2D, mostly useful for making fancy plots. More...
void	KolmogorovSmirnovTest (const std::vector< std::unique_ptr< MCMCProcessor >> &Processor, const std::unique_ptr< TCanvas > &Posterior, const TString &canvasname)
	KS: Perform KS test to check if two posteriors for the same parameter came from the same distribution. More...
int	main (int argc, char *argv[])
std::map< std::string, std::pair< double, double > >	GetCustomBinning (const YAML::Node &Settings)
	Parse custom binning edges from a YAML configuration node. More...

Variables
int	nFiles
std::vector< std::string >	FileNames
std::vector< std::string >	TitleNames
std::string	config

Detailed Description

Main exectable responsible for different types of MCMC processing like drawing posteriors, triangle plots etc. Actual implantation of methods is in MCMCProcessor.

Author

Kamil Skwarczynski

Definition in file ProcessMCMC.cpp.

Function Documentation

CalcBayesFactor()

void CalcBayesFactor ( MCMCProcessor *  Processor )

inline

KS: Calculate Bayes factor for a given hypothesis, most informative are those related to osc params. However, it make relative easy interpretation for switch dials.

Definition at line 420 of file ProcessMCMC.cpp.

{
  YAML::Node card_yaml = M3OpenConfig(config.c_str());
  YAML::Node Settings = card_yaml["ProcessMCMC"];
 
  std::vector<std::string> ParNames;
  std::vector<std::vector<double>> Model1Bounds;
  std::vector<std::vector<double>> Model2Bounds;
  std::vector<std::vector<std::string>> ModelNames;
  for (const auto& dg : Settings["BayesFactor"])
  {
    ParNames.push_back(dg[0].as<std::string>());
    ModelNames.push_back(dg[1].as<std::vector<std::string>>());
    Model1Bounds.push_back(dg[2].as<std::vector<double>>());
    Model2Bounds.push_back(dg[3].as<std::vector<double>>());
  }
 
  Processor->GetBayesFactor(ParNames, Model1Bounds, Model2Bounds, ModelNames);
}

CalcBipolarPlot()

void CalcBipolarPlot ( MCMCProcessor *  Processor )

inline

Definition at line 473 of file ProcessMCMC.cpp.

{
  YAML::Node card_yaml = M3OpenConfig(config.c_str());
  YAML::Node Settings = card_yaml["ProcessMCMC"];
 
  std::vector<std::string> ParNames;
  for (const auto& d : Settings["BipolarPlot"])
  {
    ParNames.push_back(d[0].as<std::string>());
  }
  Processor->GetPolarPlot(ParNames);
}

CalcParameterEvolution()

void CalcParameterEvolution ( MCMCProcessor *  Processor )

inline

Definition at line 458 of file ProcessMCMC.cpp.

{
  YAML::Node card_yaml = M3OpenConfig(config.c_str());
  YAML::Node Settings = card_yaml["ProcessMCMC"];
 
  std::vector<std::string> ParNames;
  std::vector<int> Intervals;
  for (const auto& d : Settings["ParameterEvolution"])
  {
    ParNames.push_back(d[0].as<std::string>());
    Intervals.push_back(d[1].as<int>());
  }
  Processor->ParameterEvolution(ParNames, Intervals);
}

CalcSavageDickey()

void CalcSavageDickey ( MCMCProcessor *  Processor )

inline

Definition at line 440 of file ProcessMCMC.cpp.

{
  YAML::Node card_yaml = M3OpenConfig(config.c_str());
  YAML::Node Settings = card_yaml["ProcessMCMC"];
 
  std::vector<std::string> ParNames;
  std::vector<double> EvaluationPoint;
  std::vector<std::vector<double>> Bounds;
  
  for (const auto& d : Settings["SavageDickey"])
  {
    ParNames.push_back(d[0].as<std::string>());
    EvaluationPoint.push_back(d[1].as<double>());
    Bounds.push_back(d[2].as<std::vector<double>>());
  }
  Processor->GetSavageDickey(ParNames, EvaluationPoint, Bounds);
}

DiagnoseCovarianceMatrix()

void DiagnoseCovarianceMatrix ( MCMCProcessor *  Processor, const std::string &  inputFile  )

inline

KS: You validate stability of posterior covariance matrix, you set burn calc cov matrix increase burn calc again and compare. By performing such operation several hundred times we can check when matrix becomes stable.

Definition at line 506 of file ProcessMCMC.cpp.

{
  //Turn of plots from Processor
  Processor->SetPrintToPDF(false);
  // Open a TCanvas to write the posterior onto
  auto Canvas = std::make_unique<TCanvas>("Canvas", "Canvas", 0, 0, 1024, 1024);
  Canvas->SetGrid();
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  Canvas->SetTickx();
  Canvas->SetTicky();
  Canvas->SetBottomMargin(0.1f);
  Canvas->SetTopMargin(0.05f);
  Canvas->SetRightMargin(0.15f);
  Canvas->SetLeftMargin(0.10f);
  
  //KS: Fancy colours
  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);
  
  std::string OutName = inputFile;
  OutName = OutName.substr(0, OutName.find(".root"));
  Canvas->Print(Form("Correlation_%s.pdf[", OutName.c_str()), "pdf");
  Canvas->Print(Form("Covariance_%s.pdf[", OutName.c_str()), "pdf");
  
  YAML::Node card_yaml = M3OpenConfig(config.c_str());
  YAML::Node Settings = card_yaml["ProcessMCMC"];
 
  const int entries = int(Processor->GetnSteps());
  const int NIntervals = GetFromManager<int>(Settings["NIntervals"], 5);
  const int IntervalsSize = entries/NIntervals;
  //We start with burn from 0 (no burn in at all)
  int BurnIn = 0;
  MACH3LOG_INFO("Diagnosing matrices with entries={}, NIntervals={} and IntervalsSize={}", entries, NIntervals, IntervalsSize);
 
  TMatrixDSym *Covariance = nullptr;
  TMatrixDSym *Correlation = nullptr;
  
  TH2D *CovariancePreviousHist = nullptr;
  TH2D *CorrelationPreviousHist = nullptr;
 
  TH2D *CovarianceHist = nullptr;
  TH2D *CorrelationHist = nullptr;
 
  //KS: Get first covariances, we need two for comparison...
  Processor->SetStepCut(BurnIn);
  Processor->GetCovariance(Covariance, Correlation);
  
  CovariancePreviousHist = TMatrixIntoTH2D(Covariance, "Covariance"); 
  CorrelationPreviousHist = TMatrixIntoTH2D(Correlation, "Correlation");
      
  delete Covariance;
  Covariance = nullptr;
  delete Correlation;
  Correlation = nullptr;
  
  //KS: Loop over all desired cuts
  for(int k = 1; k < NIntervals; ++k)
  {
    BurnIn = k*IntervalsSize;
    Processor->SetStepCut(BurnIn);
    Processor->GetCovariance(Covariance, Correlation);
    Processor->ResetHistograms();
    
    CovarianceHist = TMatrixIntoTH2D(Covariance, "Covariance"); 
    CorrelationHist = TMatrixIntoTH2D(Correlation, "Correlation");
            
    TH2D *CovarianceDiff = static_cast<TH2D*>(CovarianceHist->Clone("Covariance_Ratio"));
    TH2D *CorrelationDiff = static_cast<TH2D*>(CorrelationHist->Clone("Correlation_Ratio"));
    
    //KS: Bit messy but quite often covariance is 0 is divided by 0 is problematic so
    #ifdef MULTITHREAD
    #pragma omp parallel for
    #endif
    for (int j = 1; j < CovarianceDiff->GetXaxis()->GetNbins()+1; ++j)
    {
      for (int i = 1; i < CovarianceDiff->GetYaxis()->GetNbins()+1; ++i)
      {
        if( std::fabs (CovarianceDiff->GetBinContent(j, i)) < 1.e-5 && std::fabs (CovariancePreviousHist->GetBinContent(j, i)) < 1.e-5)
        {
          CovarianceDiff->SetBinContent(j, i, M3::_BAD_DOUBLE_);
          CovariancePreviousHist->SetBinContent(j, i, M3::_BAD_DOUBLE_);
        }
        if( std::fabs (CorrelationDiff->GetBinContent(j, i)) < 1.e-5 && std::fabs (CorrelationPreviousHist->GetBinContent(j, i)) < 1.e-5)
        {
          CorrelationDiff->SetBinContent(j, i, M3::_BAD_DOUBLE_);
          CorrelationPreviousHist->SetBinContent(j, i, M3::_BAD_DOUBLE_);
        }
      }
    }
    //Divide matrices
    CovarianceDiff->Divide(CovariancePreviousHist);
    CorrelationDiff->Divide(CorrelationPreviousHist);
    
    //Now it is time for fancy names etc.
    for (int j = 0; j < CovarianceDiff->GetXaxis()->GetNbins(); ++j)
    {
      TString Title = "";
      double Prior = 1.0;
      double PriorError = 1.0;
  
      Processor->GetNthParameter(j, Prior, PriorError, Title);
      
      CovarianceDiff->GetXaxis()->SetBinLabel(j+1, Title);
      CovarianceDiff->GetYaxis()->SetBinLabel(j+1, Title);
      CorrelationDiff->GetXaxis()->SetBinLabel(j+1, Title);
      CorrelationDiff->GetYaxis()->SetBinLabel(j+1, Title);
    }
    CovarianceDiff->GetXaxis()->SetLabelSize(0.015f);
    CovarianceDiff->GetYaxis()->SetLabelSize(0.015f);
    CorrelationDiff->GetXaxis()->SetLabelSize(0.015f);
    CorrelationDiff->GetYaxis()->SetLabelSize(0.015f);
    
    std::stringstream ss;
    ss << "BCut_";
    ss << BurnIn;
    ss << "/";
    ss << "BCut_";
    ss << (k-1)*IntervalsSize;
    std::string str = ss.str();
    
    TString Title = "Cov " + str;
    CovarianceDiff->GetZaxis()->SetTitle( Title );
    Title = "Corr " + str;
    CorrelationDiff->GetZaxis()->SetTitle(Title);
    
    CovarianceDiff->SetMinimum(-2);
    CovarianceDiff->SetMaximum(2);
    CorrelationDiff->SetMinimum(-2);
    CorrelationDiff->SetMaximum(2);
 
    Canvas->cd();
    CovarianceDiff->Draw("colz");
    Canvas->Print(Form("Covariance_%s.pdf", OutName.c_str()), "pdf");
 
    CorrelationDiff->Draw("colz");
    Canvas->Print(Form("Correlation_%s.pdf", OutName.c_str()), "pdf");
        
    //KS: Current hist become previous as we need it for further comparison
    delete CovariancePreviousHist;
    CovariancePreviousHist = static_cast<TH2D*>(CovarianceHist->Clone());
    delete CorrelationPreviousHist;
    CorrelationPreviousHist = static_cast<TH2D*>(CorrelationHist->Clone());
    
    delete CovarianceHist;
    CovarianceHist = nullptr;
    delete CorrelationHist;
    CorrelationHist = nullptr;
    
    delete CovarianceDiff;
    delete CorrelationDiff;
    delete Covariance;
    Covariance = nullptr;
    delete Correlation;
    Correlation = nullptr;
  }
  Canvas->cd();
  Canvas->Print(Form("Covariance_%s.pdf]", OutName.c_str()), "pdf");
  Canvas->Print(Form("Correlation_%s.pdf]", OutName.c_str()), "pdf");
  
  Processor->SetPrintToPDF(true);
  if(Covariance != nullptr)              delete Covariance;
  if(Correlation != nullptr)             delete Correlation;
  if(CovariancePreviousHist != nullptr)  delete CovariancePreviousHist;
  if(CorrelationPreviousHist != nullptr) delete CorrelationPreviousHist;
  if(CovarianceHist != nullptr)          delete CovarianceHist;
  if(CorrelationHist != nullptr)         delete CorrelationHist;
}

GetCustomBinning()

std::map<std::string, std::pair<double, double> > GetCustomBinning

(

const YAML::Node &

Settings

)

Parse custom binning edges from a YAML configuration node.

This function reads the CustomBinEdges section of the YAML config and returns a mapping of parameter names to their lower and upper edges.

The expected YAML syntax is:

CustomBinEdges:
  delta_cp: [-3.141592, 3.141592]
  another_param: [min, max]

Parameters

Settings

YAML configuration node containing the optional CustomBinEdges section.

Definition at line 92 of file ProcessMCMC.cpp.

{
  std::map<std::string, std::pair<double, double>> CustomBinning;
  if (Settings["CustomBinEdges"]) {
    const YAML::Node& edges = Settings["CustomBinEdges"];
 
    for (const auto& node : edges) {
      std::string key = node.first.as<std::string>();
      auto values = node.second.as<std::vector<double>>();
 
      if (values.size() == 2) {
        CustomBinning[key] = std::make_pair(values[0], values[1]);
        MACH3LOG_DEBUG("Adding custom binning {} with {:.4f}, {:.4f}", key, values[0], values[1]);
      } else {
        MACH3LOG_ERROR("Invalid number of values for key: {}", key);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
    }
  }
  return CustomBinning;
}

GetTrianglePlot()

void GetTrianglePlot ( MCMCProcessor *  Processor )

inline

Definition at line 486 of file ProcessMCMC.cpp.

                                               {
  YAML::Node card_yaml = M3OpenConfig(config.c_str());
  YAML::Node Settings = card_yaml["ProcessMCMC"];
 
  for (const auto& dg : Settings["TrianglePlot"])
  {
    std::string ParName = dg[0].as<std::string>();
 
    std::vector<std::string> NameVec = dg[1].as<std::vector<std::string>>();
    Processor->MakeTrianglePlot(NameVec,
                                GetFromManager<std::vector<double>>(Settings["CredibleIntervals"], {0.99, 0.90, 0.68}),
                                GetFromManager<std::vector<short int>>(Settings["CredibleIntervalsColours"], {436, 430, 422}),
                                GetFromManager<std::vector<double>>(Settings["CredibleRegions"], {0.99, 0.90, 0.68}),
                                GetFromManager<std::vector<short int>>(Settings["CredibleRegionStyle"], {2, 1, 3}),
                                GetFromManager<std::vector<short int>>(Settings["CredibleRegionColor"], {413, 406, 416}),
                                GetFromManager<bool>(Settings["CredibleInSigmas"], false));
  }
}

KolmogorovSmirnovTest()

void KolmogorovSmirnovTest ( const std::vector< std::unique_ptr< MCMCProcessor >> &  Processor, const std::unique_ptr< TCanvas > &  Posterior, const TString &  canvasname  )

inline

KS: Perform KS test to check if two posteriors for the same parameter came from the same distribution.

Definition at line 710 of file ProcessMCMC.cpp.

{
  constexpr Color_t CumulativeColor[] = {kBlue-1, kRed, kGreen+2};
  constexpr Style_t CumulativeStyle[] = {kSolid, kDashed, kDotted};
 
  for(int i = 0; i < Processor[0]->GetNParams(); ++i) 
  {
    // This holds the posterior density
    std::vector<std::unique_ptr<TH1D>> hpost(nFiles);
    std::vector<std::unique_ptr<TH1D>> CumulativeDistribution(nFiles);
       
    TString Title;
    double Prior = 1.0;
    double PriorError = 1.0;
 
    Processor[0]->GetNthParameter(i, Prior, PriorError, Title);
    bool Skip = false;
    for (int ik = 0 ; ik < nFiles;  ik++)
    {
      int Index = 0;
      if(ik == 0 ) Index = i;
      else
      {
        // KS: If somehow this chain doesn't given params we skip it
        Index = Processor[ik]->GetParamIndexFromName(hpost[0]->GetTitle());
        if(Index == M3::_BAD_INT_)
        {
          Skip = true;
          break;
        }
      }
      hpost[ik] = M3::Clone(Processor[ik]->GetHpost(Index));
      CumulativeDistribution[ik] = M3::Clone(Processor[ik]->GetHpost(Index));
      CumulativeDistribution[ik]->Fill(0., 0.);
      CumulativeDistribution[ik]->Reset();
      CumulativeDistribution[ik]->SetMaximum(1.);
      TString TempTitle = Title+" Kolmogorov Smirnov";
      CumulativeDistribution[ik]->SetTitle(TempTitle);
      
      TempTitle = Title+" Value";
      CumulativeDistribution[ik]->GetXaxis()->SetTitle(TempTitle);
      CumulativeDistribution[ik]->GetYaxis()->SetTitle("Cumulative Probability");
      
      CumulativeDistribution[ik]->SetLineWidth(2);
      CumulativeDistribution[ik]->SetLineColor(CumulativeColor[ik]);
      CumulativeDistribution[ik]->SetLineStyle(CumulativeStyle[ik]);
    }
 
    // Don't plot if this is a fixed histogram (i.e. the peak is the whole integral)
    if(hpost[0]->GetMaximum() == hpost[0]->Integral()*1.5 || Skip) {
      continue;
    }
    
    for (int ik = 0 ; ik < nFiles;  ik++)
    {
      const int NumberOfBins = hpost[ik]->GetXaxis()->GetNbins();
      double Cumulative = 0;
      const double Integral = hpost[ik]->Integral();
      for (int j = 1; j < NumberOfBins+1; ++j)
      {
        Cumulative += hpost[ik]->GetBinContent(j)/Integral;
        CumulativeDistribution[ik]->SetBinContent(j, Cumulative);
      }
      //KS: Set overflow to 1 just in case
      CumulativeDistribution[ik]->SetBinContent(NumberOfBins+1, 1.);
    }
    
    std::vector<int> TestStatBin(nFiles, 0);
    std::vector<double> TestStatD(nFiles, -999);
    std::vector<std::unique_ptr<TLine>> LineD(nFiles);
    //Find KS statistic
    for (int ik = 1 ; ik < nFiles;  ik++)
    {
      const int NumberOfBins = CumulativeDistribution[0]->GetXaxis()->GetNbins();
      for (int j = 1; j < NumberOfBins+1; ++j)
      {
        const double BinValue = CumulativeDistribution[0]->GetBinCenter(j);
        const int BinNumber = CumulativeDistribution[ik]->FindBin(BinValue);
        //KS: Calculate D statistic for this bin, only save it if it's bigger than previously found value
        double TempDstat = std::fabs(CumulativeDistribution[0]->GetBinContent(j) - CumulativeDistribution[ik]->GetBinContent(BinNumber));
        if(TempDstat > TestStatD[ik])
        {
          TestStatD[ik] = TempDstat;
          TestStatBin[ik] = j;
        }
      }
    }
 
    for (int ik = 0 ; ik < nFiles;  ik++)
    {
      LineD[ik] = std::make_unique<TLine>(CumulativeDistribution[0]->GetBinCenter(TestStatBin[ik]), 0, CumulativeDistribution[0]->GetBinCenter(TestStatBin[ik]), CumulativeDistribution[0]->GetBinContent(TestStatBin[ik]));
      LineD[ik]->SetLineColor(CumulativeColor[ik]);
      LineD[ik]->SetLineWidth(2.0);
    }
    CumulativeDistribution[0]->Draw();
    for (int ik = 0 ; ik < nFiles;  ik++)
      CumulativeDistribution[ik]->Draw("SAME");
    
    auto leg = std::make_unique<TLegend>(0.15, 0.7, 0.5, 0.90);
    leg->SetTextSize(0.04f);
    for (int ik = 0; ik < nFiles;  ik++)
      leg->AddEntry(CumulativeDistribution[ik].get(), TitleNames[ik].c_str(), "l");
    for (int ik = 1; ik < nFiles;  ik++)
      leg->AddEntry(LineD[ik].get(), Form("#Delta D = %.4f", TestStatD[ik]), "l");
    
    leg->SetLineColor(0);
    leg->SetLineStyle(0);
    leg->SetFillColor(0);
    leg->SetFillStyle(0);
    leg->Draw("SAME");
 
    for (int ik = 1; ik < nFiles;  ik++)
      LineD[ik]->Draw("sam");
      
    Posterior->cd();
    Posterior->Print(canvasname);
  } //End loop over parameter
}

main()

int main	(	int	argc,
		char *	argv[]
	)

Definition at line 31 of file ProcessMCMC.cpp.

{
  SetMaCh3LoggerFormat();
  nFiles = 0;
  if (argc != 3 && argc !=6 && argc != 8)
  {
    MACH3LOG_ERROR("How to use: ");
    MACH3LOG_ERROR("  single chain: {} <Config> <MCMM_ND_Output.root>", argv[0]);
    MACH3LOG_ERROR("  two chain:    {} <Config> <MCMM_ND_Output_1.root> <Title 1> <MCMC_ND_Output_2.root> <Title 2>", argv[0]);
    MACH3LOG_ERROR("  three chain:  {} <Config> <MCMM_ND_Output_1.root> <Title 1> <MCMC_ND_Output_2.root> <Title 2> <MCMC_ND_Output_3.root> <Title 3>", argv[0]);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  config = argv[1];
  YAML::Node card_yaml = M3OpenConfig(config);
  if (!CheckNodeExists(card_yaml, "ProcessMCMC")) {
    MACH3LOG_ERROR("The 'ProcessMCMC' node is not defined in the YAML configuration.");
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  if (argc == 3)
  {
    MACH3LOG_INFO("Producing single fit output");
    std::string filename = argv[2];
    ProcessMCMC(filename);
  }
  // If we want to compare two or more fits (e.g. binning changes or introducing new params/priors)
  else if (argc == 6 || argc == 8)
  {
    MACH3LOG_INFO("Producing two fit comparison");
    FileNames.push_back(argv[2]);
    TitleNames.push_back(argv[3]);
 
    FileNames.push_back(argv[4]);
    TitleNames.push_back(argv[5]);
    //KS: If there is third file add it
    if(argc == 8)
    {
      FileNames.push_back(argv[6]);
      TitleNames.push_back(argv[7]);
    }
 
    MultipleProcessMCMC();
  }
  return 0;
}

MultipleProcessMCMC()

void MultipleProcessMCMC ( )

inline

Function producing comparison of posterior and more betwen a few MCMC chains.

Definition at line 211 of file ProcessMCMC.cpp.

{
  YAML::Node card_yaml = M3OpenConfig(config.c_str());
  YAML::Node Settings = card_yaml["ProcessMCMC"];
 
  constexpr Color_t PosteriorColor[] = {kBlue-1, kRed, kGreen+2};
  //constexpr Style_t PosteriorStyle[] = {kSolid, kDashed, kDotted};
  nFiles = int(FileNames.size());
  std::vector<std::unique_ptr<MCMCProcessor>> Processor(nFiles);
 
  if(!Settings["BurnInSteps"]) {
    MACH3LOG_WARN("BurnInSteps not set, defaulting to 20%");
  }
 
  for (int ik = 0; ik < nFiles; ik++)
  {
    MACH3LOG_INFO("File for study: {}", FileNames[ik]);
    // Make the processor
    Processor[ik] = std::make_unique<MCMCProcessor>(FileNames[ik]);
    Processor[ik]->SetOutputSuffix(("_" + std::to_string(ik)).c_str());
 
    Processor[ik]->SetExcludedTypes(GetFromManager<std::vector<std::string>>(Settings["ExcludedTypes"], {}));
    Processor[ik]->SetExcludedNames(GetFromManager<std::vector<std::string>>(Settings["ExcludedNames"], {}));
    Processor[ik]->SetExcludedGroups(GetFromManager<std::vector<std::string>>(Settings["ExcludedGroups"], {}));
 
    //Apply additional cuts to 1D posterior
    Processor[ik]->SetPosterior1DCut(GetFromManager<std::string>(Settings["Posterior1DCut"], ""));
 
    Processor[ik]->SetPlotRelativeToPrior(GetFromManager<bool>(Settings["PlotRelativeToPrior"], false));
    Processor[ik]->SetFancyNames(GetFromManager<bool>(Settings["FancyNames"], true));
    Processor[ik]->Initialise();
 
    if(Settings["BurnInSteps"]) {
      Processor[ik]->SetStepCut(Settings["BurnInSteps"].as<int>());
    }else {
      Processor[ik]->SetStepCut(static_cast<int>(Processor[ik]->GetnSteps()/5));
    }
 
    if(Settings["MaxEntries"]) {
      Processor[ik]->SetEntries(Get<int>(Settings["MaxEntries"], __FILE__, __LINE__));
    }
    if(Settings["NBins"]) {
      Processor[ik]->SetNBins(Get<int>(Settings["NBins"], __FILE__, __LINE__));
    }
  }
 
  Processor[0]->MakePostfit(GetCustomBinning(Settings));
  Processor[0]->DrawPostfit();
  // Get edges from first histogram to ensure all params use same binning
  std::map<std::string, std::pair<double, double>> ParamEdges;
  for(int i = 0; i < Processor[0]->GetNParams(); ++i) {
    // Get the histogram for the i-th parameter
    TH1D* hist = Processor[0]->GetHpost(i);
    if (!hist) {
      MACH3LOG_DEBUG("Histogram for parameter {} is null.", i);
      continue;
    }
 
    // Get the parameter name (title of the histogram)
    std::string paramName = hist->GetTitle();
 
    // Get the axis limits (edges)
    TAxis* axis = hist->GetXaxis();
    double xmin = axis->GetXmin();
    double xmax = axis->GetXmax();
 
    MACH3LOG_DEBUG("Adding bin edges for {} equal to {:.4f}, {:.4f}",paramName, xmin, xmax);
    // Insert into the map
    ParamEdges[paramName] = std::make_pair(xmin, xmax);
  }
 
  //KS: Multithreading here is very tempting but there are some issues with root that need to be resovled :(
  for (int ik = 1; ik < nFiles; ik++)
  {
    // Make the postfit
    Processor[ik]->MakePostfit(ParamEdges);
    Processor[ik]->DrawPostfit();
  }
 
  // Open a TCanvas to write the posterior onto
  auto Posterior = std::make_unique<TCanvas>("PosteriorMulti", "PosteriorMulti", 0, 0, 1024, 1024);
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  Posterior->SetGrid();
  Posterior->SetBottomMargin(0.1f);
  Posterior->SetTopMargin(0.05f);
  Posterior->SetRightMargin(0.03f);
  Posterior->SetLeftMargin(0.15f);
  
  // First filename: keep path, just remove ".root"
  // Would be nice to specify outpath in a later update
  size_t pos = FileNames[0].rfind(".root");
  std::string base = (pos == std::string::npos) ? FileNames[0] : FileNames[0].substr(0, pos);
  TString canvasname = base;
 
  // Remaining filenames: strip path and ".root"
  // So if you have /path/to/file1.root and /path/to/file2.root or /another/path/to/file2.root, canvasname = /path/to/file1_file2.root
  for (int ik = 1; ik < nFiles; ik++) {
    pos = FileNames[ik].find_last_of('/');
    base = (pos == std::string::npos) ? FileNames[ik] : FileNames[ik].substr(pos + 1);
    pos = base.rfind(".root");
    if (pos != std::string::npos) base = base.substr(0, pos);
    canvasname += "_" + TString(base);
  }
  
  canvasname = canvasname +".pdf[";
 
  Posterior->Print(canvasname);
  // Once the pdf file is open no longer need to bracket
  canvasname.ReplaceAll("[","");
 
  for(int i = 0; i < Processor[0]->GetNParams(); ++i) 
  {
    // This holds the posterior density
    std::vector<std::unique_ptr<TH1D>> hpost(nFiles);
    std::vector<std::unique_ptr<TLine>> hpd(nFiles);
    hpost[0] = M3::Clone(Processor[0]->GetHpost(i));
    hpost[0]->GetYaxis()->SetTitle("Posterior Density");
    bool Skip = false;
    for (int ik = 1 ; ik < nFiles;  ik++)
    {
      // KS: If somehow this chain doesn't given params we skip it
      const int Index = Processor[ik]->GetParamIndexFromName(hpost[0]->GetTitle());
      if(Index == M3::_BAD_INT_)
      {
        Skip = true;
        break;
      }
      hpost[ik] = M3::Clone(Processor[ik]->GetHpost(Index));
    }
 
    // Don't plot if this is a fixed histogram (i.e. the peak is the whole integral)
    if(hpost[0]->GetMaximum() == hpost[0]->Integral()*1.5 || Skip) {
      continue;
    }
    for (int ik = 0; ik < nFiles;  ik++)
    {
      RemoveFitter(hpost[ik].get(), "Gauss");
 
      // Set some nice colours
      hpost[ik]->SetLineColor(PosteriorColor[ik]);
      //hpost[ik]->SetLineStyle(PosteriorStyle[ik]);
      hpost[ik]->SetLineWidth(2);
 
      // Area normalise the distributions
      hpost[ik]->Scale(1./hpost[ik]->Integral());
    }
    TString Title;
    double Prior = 1.0;
    double PriorError = 1.0;
 
    Processor[0]->GetNthParameter(i, Prior, PriorError, Title);
 
    // Now make the TLine for the Asimov
    auto Asimov = std::make_unique<TLine>(Prior, hpost[0]->GetMinimum(), Prior, hpost[0]->GetMaximum());
    Asimov->SetLineColor(kRed-3);
    Asimov->SetLineWidth(2);
    Asimov->SetLineStyle(kDashed);
 
    // Make a nice little TLegend
    auto leg = std::make_unique<TLegend>(0.20, 0.7, 0.6, 0.97);
    leg->SetTextSize(0.03f);
    leg->SetFillColor(0);
    leg->SetFillStyle(0);
    leg->SetLineColor(0);
    leg->SetLineStyle(0);
    TString asimovLeg = Form("#splitline{Prior}{x = %.2f , #sigma = %.2f}", Prior, PriorError);
    leg->AddEntry(Asimov.get(), asimovLeg, "l");
 
    for (int ik = 0; ik < nFiles; ik++)
    {
      TString rebinLeg = Form("#splitline{%s}{#mu = %.2f, #sigma = %.2f}", TitleNames[ik].c_str(), hpost[ik]->GetMean(), hpost[ik]->GetRMS());
      leg->AddEntry(hpost[ik].get(),  rebinLeg, "l");
 
      hpd[ik] = std::make_unique<TLine>(hpost[ik]->GetBinCenter(hpost[ik]->GetMaximumBin()), hpost[ik]->GetMinimum(),
                                        hpost[ik]->GetBinCenter(hpost[ik]->GetMaximumBin()), hpost[ik]->GetMaximum());
      hpd[ik]->SetLineColor(hpost[ik]->GetLineColor());
      hpd[ik]->SetLineWidth(2);
      hpd[ik]->SetLineStyle(kSolid);
    }
 
    // Find the maximum value to nicely resize hist
    double maximum = 0;
    for (int ik = 0; ik < nFiles; ik++) maximum = std::max(maximum, hpost[ik]->GetMaximum());
    for (int ik = 0; ik < nFiles; ik++) hpost[ik]->SetMaximum(1.3*maximum);
 
    hpost[0]->Draw("hist");
    for (int ik = 1; ik < nFiles; ik++) hpost[ik]->Draw("hist same");
    Asimov->Draw("same");
    for (int ik = 0; ik < nFiles; ik++) hpd[ik]->Draw("same");
    leg->Draw("same");
    Posterior->cd();
    Posterior->Print(canvasname);
  }//End loop over parameters
    
  // Finally draw the parameter plot onto the PDF
  // Close the .pdf file with all the posteriors
  Posterior->cd();
  Posterior->Clear();
 
  if(GetFromManager<bool>(Settings["PerformKStest"], true)) KolmogorovSmirnovTest(Processor, Posterior, canvasname);
  
  // Close the pdf file
  MACH3LOG_INFO("Closing pdf {}", canvasname);
  canvasname+="]";
  Posterior->Print(canvasname);
}

ProcessMCMC()

void ProcessMCMC ( const std::string &  inputFile )

inline

Main function processing MCMC and Producing plots.

Definition at line 114 of file ProcessMCMC.cpp.

{
  MACH3LOG_INFO("File for study: {} with config  {}", inputFile, config);
  // Make the processor)
  auto Processor = std::make_unique<OscProcessor>(inputFile);
 
  YAML::Node card_yaml = M3OpenConfig(config.c_str());
  YAML::Node Settings = card_yaml["ProcessMCMC"];
 
  const bool PlotCorr = GetFromManager<bool>(Settings["PlotCorr"], false);
 
  Processor->SetExcludedTypes(GetFromManager<std::vector<std::string>>(Settings["ExcludedTypes"], {}));
  Processor->SetExcludedNames(GetFromManager<std::vector<std::string>>(Settings["ExcludedNames"], {}));
  Processor->SetExcludedGroups(GetFromManager<std::vector<std::string>>(Settings["ExcludedGroups"], {}));
 
  //Apply additional cuts to 1D posterior
  Processor->SetPosterior1DCut(GetFromManager<std::string>(Settings["Posterior1DCut"], ""));
 
  if(PlotCorr) Processor->SetOutputSuffix("_drawCorr");
  //KS:Turn off plotting detector and some other setting, should be via some config
  Processor->SetPlotRelativeToPrior(GetFromManager<bool>(Settings["PlotRelativeToPrior"], false));
  Processor->SetPrintToPDF(GetFromManager<bool>(Settings["PrintToPDF"], true));
 
  //KS: Whether you want prior error bands for parameters with flat prior or not
  Processor->SetPlotErrorForFlatPrior(GetFromManager<bool>(Settings["PlotErrorForFlatPrior"], true));
  Processor->SetFancyNames(GetFromManager<bool>(Settings["FancyNames"], true));
  Processor->SetPlotBinValue(GetFromManager<bool>(Settings["PlotBinValue"], false));
  //KS: Plot only 2D posteriors with correlations greater than 0.2
  Processor->SetPost2DPlotThreshold(GetFromManager<double>(Settings["Post2DPlotThreshold"], 0.2));
 
  Processor->Initialise();
 
  if(Settings["BurnInSteps"]) {
    Processor->SetStepCut(Settings["BurnInSteps"].as<int>());
  } else {
    MACH3LOG_WARN("BurnInSteps not set, defaulting to 20%");
    Processor->SetStepCut(static_cast<int>(Processor->GetnSteps()/5));
  }
  if(Settings["MaxEntries"]) {
    Processor->SetEntries(Get<int>(Settings["MaxEntries"], __FILE__, __LINE__));
  }
  if(Settings["NBins"]) {
    Processor->SetNBins(Get<int>(Settings["NBins"], __FILE__, __LINE__));
  }
  if(Settings["Thinning"])
  {
    if(Settings["Thinning"][0].as<bool>()){
      Processor->ThinMCMC(Settings["Thinning"][1].as<int>());
    }
  }
  // Make the postfit
  Processor->MakePostfit(GetCustomBinning(Settings));
  Processor->DrawPostfit();
  //KS: Should set via config whether you want below or not
  if(GetFromManager<bool>(Settings["MakeCredibleIntervals"], true)) {
    Processor->MakeCredibleIntervals(GetFromManager<std::vector<double>>(Settings["CredibleIntervals"], {0.99, 0.90, 0.68}),
                                     GetFromManager<std::vector<short int>>(Settings["CredibleIntervalsColours"], {436, 430, 422}),
                                     GetFromManager<bool>(Settings["CredibleInSigmas"], false));
  }
  if(GetFromManager<bool>(Settings["CalcBayesFactor"], true))  CalcBayesFactor(Processor.get());
  if(GetFromManager<bool>(Settings["CalcSavageDickey"], true)) CalcSavageDickey(Processor.get());
  if(GetFromManager<bool>(Settings["CalcBipolarPlot"], false)) CalcBipolarPlot(Processor.get());
  if(GetFromManager<bool>(Settings["CalcParameterEvolution"], false)) CalcParameterEvolution(Processor.get());
 
  if(PlotCorr)
  {
    Processor->SetSmoothing(GetFromManager<bool>(Settings["Smoothing"], true));
    // Make the covariance matrix
    //We have different treatment for multithread
    Processor->CacheSteps();
    //KS: Since we cached let's make fancy violins :)
    if(GetFromManager<bool>(Settings["MakeViolin"], true)) Processor->MakeViolin();
    Processor->MakeCovariance_MP();
 
    Processor->DrawCovariance();
    if(GetFromManager<bool>(Settings["MakeCovarianceYAML"], true)) Processor->MakeCovarianceYAML(GetFromManager<std::string>(Settings["CovarianceYAMLOutName"], "UpdatedCorrelationMatrix.yaml"), GetFromManager<std::string>(Settings["CovarianceYAMLMeansMethod"], "HPD"));
 
    auto const &MakeSubOptimality = Settings["MakeSubOptimality"];
    if(MakeSubOptimality[0].as<bool>()) Processor->MakeSubOptimality(MakeSubOptimality[1].as<int>());
 
    if(GetFromManager<bool>(Settings["MakeCredibleRegions"], false)) {
      Processor->MakeCredibleRegions(GetFromManager<std::vector<double>>(Settings["CredibleRegions"], {0.99, 0.90, 0.68}),
                                     GetFromManager<std::vector<short int>>(Settings["CredibleRegionStyle"], {2, 1, 3}),
                                     GetFromManager<std::vector<short int>>(Settings["CredibleRegionColor"], {413, 406, 416}),
                                     GetFromManager<bool>(Settings["CredibleInSigmas"], false), 
                                     GetFromManager<bool>(Settings["Draw2DPosterior"], true),
                                     GetFromManager<bool>(Settings["DrawBestFit"], true));
    }
    if(GetFromManager<bool>(Settings["GetTrianglePlot"], true)) GetTrianglePlot(Processor.get());
 
    //KS: When creating covariance matrix longest time is spend on caching every step, since we already cached we can run some fancy covariance stability diagnostic
    if(GetFromManager<bool>(Settings["DiagnoseCovarianceMatrix"], false)) DiagnoseCovarianceMatrix(Processor.get(), inputFile);
  }
  if(GetFromManager<bool>(Settings["ReweightPrior"], false)) ReweightPrior(Processor.get());
  if(GetFromManager<bool>(Settings["JarlskogAnalysis"], true)) Processor->PerformJarlskogAnalysis();
}

ReweightPrior()

void ReweightPrior ( MCMCProcessor *  Processor )

inline

Definition at line 681 of file ProcessMCMC.cpp.

{
  YAML::Node card_yaml = M3OpenConfig(config.c_str());
  YAML::Node Settings = card_yaml["ProcessMCMC"];
 
  const auto& Prior = Settings["PriorReweighting"];
  std::vector<std::string> Names = Prior[0].as<std::vector<std::string>>();
  std::vector<double> NewCentral = Prior[1].as<std::vector<double>>();
  std::vector<double> NewError = Prior[2].as<std::vector<double>>();
 
  Processor->ReweightPrior(Names, NewCentral, NewError);
}

TMatrixIntoTH2D()

TH2D * TMatrixIntoTH2D ( TMatrixDSym *  Matrix, const std::string &  title  )

inline

KS: Convert TMatrix to TH2D, mostly useful for making fancy plots.

Definition at line 695 of file ProcessMCMC.cpp.

{
  TH2D* hMatrix = new TH2D(title.c_str(), title.c_str(), Matrix->GetNrows(), 0.0, Matrix->GetNrows(), Matrix->GetNcols(), 0.0, Matrix->GetNcols());
  for(int i = 0; i < Matrix->GetNrows(); i++)
  {
    for(int j = 0; j < Matrix->GetNcols(); j++)
    {
      //KS: +1 because there is offset in histogram relative to TMatrix
      hMatrix->SetBinContent(i+1,j+1, (*Matrix)(i,j));
    }
  }
  return hMatrix;
}

Variable Documentation

config

std::string config

Definition at line 29 of file ProcessMCMC.cpp.

FileNames

std::vector<std::string> FileNames

Definition at line 27 of file ProcessMCMC.cpp.

nFiles

int nFiles

Definition at line 26 of file ProcessMCMC.cpp.

TitleNames

std::vector<std::string> TitleNames

Definition at line 28 of file ProcessMCMC.cpp.