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.
void	MultipleProcessMCMC ()
	Function producing comparison of posterior and more betwen a few MCMC chains.
void	CalcBayesFactor (MCMCProcessor *Processor)
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.
void	ReweightPrior (MCMCProcessor *Processor)
TH2D *	TMatrixIntoTH2D (TMatrixDSym *Matrix, const std::string &title)
	KS: Convert TMatrix to TH2D, mostly useful for making fancy plots.
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.
int	main (int argc, char *argv[])

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

Definition at line 348 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 401 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 386 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 368 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 434 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, _UNDEF_);
          CovariancePreviousHist->SetBinContent(j, i, _UNDEF_);
        }
        if( std::fabs (CorrelationDiff->GetBinContent(j, i)) < 1.e-5 && std::fabs (CorrelationPreviousHist->GetBinContent(j, i)) < 1.e-5)
        {
          CorrelationDiff->SetBinContent(j, i, _UNDEF_);
          CorrelationPreviousHist->SetBinContent(j, i, _UNDEF_);
        }
      }
    }
    //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;
}

GetTrianglePlot()

void GetTrianglePlot ( MCMCProcessor *  Processor )

inline

Definition at line 414 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 638 of file ProcessMCMC.cpp.

{
  const Color_t CumulativeColor[] = {kBlue-1, kRed, kGreen+2};
  const Style_t CumulativeStyle[] = {kSolid, kDashed, kDotted};
 
  for(int i = 0; i < Processor[0]->GetNParams(); ++i) 
  {
    // This holds the posterior density
    std::vector<TH1D*> hpost(nFiles);
    std::vector<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 == _UNDEF_)
        {
          Skip = true;
          break;
        }
      }
      hpost[ik] = static_cast<TH1D*>(Processor[ik]->GetHpost(Index)->Clone());
      CumulativeDistribution[ik] = static_cast<TH1D*>(Processor[ik]->GetHpost(Index)->Clone());
      CumulativeDistribution[ik]->Fill(0., 0.);
      CumulativeDistribution[ik]->Reset();
      CumulativeDistribution[ik]->SetMaximum(1.);
      TString TempTittle = Title+" Kolmogorov Smirnov";
      CumulativeDistribution[ik]->SetTitle(TempTittle);
      
      TempTittle = Title+" Value";
      CumulativeDistribution[ik]->GetXaxis()->SetTitle(TempTittle);
      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)
    {
      for (int ik = 0; ik < nFiles;  ik++)
      {
        delete hpost[ik];
        delete CumulativeDistribution[ik];
      }
      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], 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);
    
    for (int ik = 0; ik < nFiles;  ik++)
    {
      delete hpost[ik];
      delete CumulativeDistribution[ik];
    }
  } //End loop over parameter
}

main()

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

Definition at line 32 of file ProcessMCMC.cpp.

{
  SetMaCh3LoggerFormat();
  nFiles = 0;
  if (argc != 3 && argc !=6 && argc != 8)
  {
    MACH3LOG_ERROR("How to use: {} <Config> <MCMM_ND_Output.root>", argv[0]);
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  if (argc == 3)
  {
    MACH3LOG_INFO("Producing single fit output");
    config = argv[1];
    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");
    config = argv[1];
 
    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 167 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"], {}));
 
    //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));
    }
  }
  //KS: Multithreading here is very tempting but there are some issues with root that need to be resovled :(
  for (int ik = 0; ik < nFiles;  ik++)
  {
    // Make the postfit
    Processor[ik]->MakePostfit();
    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.10f);
 
  FileNames[0] = FileNames[0].substr(0, FileNames[0].find(".root")-1);
  TString canvasname = FileNames[0];
  for (int ik = 1; ik < nFiles;  ik++)
  {
    while (FileNames[ik].find("/") != std::string::npos)
    {
      FileNames[ik] = FileNames[ik].substr(FileNames[ik].find("/")+1, FileNames[ik].find(".root")-1);
    }
    canvasname = canvasname + "_"+FileNames[ik];
  }
 
  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<TH1D*> hpost(nFiles);
    std::vector<std::unique_ptr<TLine>> hpd(nFiles);
    hpost[0] = static_cast<TH1D *>(Processor[0]->GetHpost(i)->Clone());
 
    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 == _UNDEF_)
      {
        Skip = true;
        break;
      }
      hpost[ik] = static_cast<TH1D *>(Processor[ik]->GetHpost(Index)->Clone());
    }
 
    // 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)
    {
      for (int ik = 0; ik < nFiles;  ik++)
        delete hpost[ik];
 
      continue;
    }
    for (int ik = 0; ik < nFiles;  ik++)
    {
      RemoveFitter(hpost[ik], "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(), "width");
    }
    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.12, 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],  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);
    for (int ik = 0; ik < nFiles;  ik++) {
      delete hpost[ik];
    }
  }//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 72 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"], {}));
  //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["Thinning"])
  {
    if(Settings["Thinning"][0].as<bool>()){
      Processor->ThinMCMC(Settings["Thinning"][1].as<int>());
    }
  }
  // Make the postfit
  Processor->MakePostfit();
  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
//#ifdef MULTITHREAD
    Processor->CacheSteps();
    //KS: Since we cached let's make fancy violins :)
    if(GetFromManager<bool>(Settings["MakeViolin"], true)) Processor->MakeViolin();
    Processor->MakeCovariance_MP();
//#else
    //Processor->MakeCovariance();
//#endif
    Processor->DrawCovariance();
 
    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 609 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 623 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 30 of file ProcessMCMC.cpp.

FileNames

std::vector<std::string> FileNames

Definition at line 28 of file ProcessMCMC.cpp.

nFiles

int nFiles

Definition at line 27 of file ProcessMCMC.cpp.

TitleNames

std::vector<std::string> TitleNames

Definition at line 29 of file ProcessMCMC.cpp.