MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
Functions | Variables
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.

349{
350 YAML::Node card_yaml = M3OpenConfig(config.c_str());
351 YAML::Node Settings = card_yaml["ProcessMCMC"];
352
353 std::vector<std::string> ParNames;
354 std::vector<std::vector<double>> Model1Bounds;
355 std::vector<std::vector<double>> Model2Bounds;
356 std::vector<std::vector<std::string>> ModelNames;
357 for (const auto& dg : Settings["BayesFactor"])
358 {
359 ParNames.push_back(dg[0].as<std::string>());
360 ModelNames.push_back(dg[1].as<std::vector<std::string>>());
361 Model1Bounds.push_back(dg[2].as<std::vector<double>>());
362 Model2Bounds.push_back(dg[3].as<std::vector<double>>());
363 }
364
365 Processor->GetBayesFactor(ParNames, Model1Bounds, Model2Bounds, ModelNames);
366}
std::string config
Definition: ProcessMCMC.cpp:30
#define M3OpenConfig(filename)
Macro to simplify calling LoadYaml with file and line info.
Definition: YamlHelper.h:560
void GetBayesFactor(const std::vector< std::string > &ParName, const std::vector< std::vector< double > > &Model1Bounds, const std::vector< std::vector< double > > &Model2Bounds, const std::vector< std::vector< std::string > > &ModelNames)
Calculate Bayes factor for vector of params, and model boundaries.

◆ CalcBipolarPlot()

void CalcBipolarPlot ( MCMCProcessor Processor)
inline

Definition at line 401 of file ProcessMCMC.cpp.

402{
403 YAML::Node card_yaml = M3OpenConfig(config.c_str());
404 YAML::Node Settings = card_yaml["ProcessMCMC"];
405
406 std::vector<std::string> ParNames;
407 for (const auto& d : Settings["BipolarPlot"])
408 {
409 ParNames.push_back(d[0].as<std::string>());
410 }
411 Processor->GetPolarPlot(ParNames);
412}
void GetPolarPlot(const std::vector< std::string > &ParNames)
Make funny polar plot.

◆ CalcParameterEvolution()

void CalcParameterEvolution ( MCMCProcessor Processor)
inline

Definition at line 386 of file ProcessMCMC.cpp.

387{
388 YAML::Node card_yaml = M3OpenConfig(config.c_str());
389 YAML::Node Settings = card_yaml["ProcessMCMC"];
390
391 std::vector<std::string> ParNames;
392 std::vector<int> Intervals;
393 for (const auto& d : Settings["ParameterEvolution"])
394 {
395 ParNames.push_back(d[0].as<std::string>());
396 Intervals.push_back(d[1].as<int>());
397 }
398 Processor->ParameterEvolution(ParNames, Intervals);
399}
void ParameterEvolution(const std::vector< std::string > &Names, const std::vector< int > &NIntervals)
Make .gif of parameter evolution.

◆ CalcSavageDickey()

void CalcSavageDickey ( MCMCProcessor Processor)
inline

Definition at line 368 of file ProcessMCMC.cpp.

369{
370 YAML::Node card_yaml = M3OpenConfig(config.c_str());
371 YAML::Node Settings = card_yaml["ProcessMCMC"];
372
373 std::vector<std::string> ParNames;
374 std::vector<double> EvaluationPoint;
375 std::vector<std::vector<double>> Bounds;
376
377 for (const auto& d : Settings["SavageDickey"])
378 {
379 ParNames.push_back(d[0].as<std::string>());
380 EvaluationPoint.push_back(d[1].as<double>());
381 Bounds.push_back(d[2].as<std::vector<double>>());
382 }
383 Processor->GetSavageDickey(ParNames, EvaluationPoint, Bounds);
384}
void GetSavageDickey(const std::vector< std::string > &ParName, const std::vector< double > &EvaluationPoint, const std::vector< std::vector< double > > &Bounds)
Calculate Bayes factor for point like hypothesis using SavageDickey.

◆ 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.

435{
436 //Turn of plots from Processor
437 Processor->SetPrintToPDF(false);
438 // Open a TCanvas to write the posterior onto
439 auto Canvas = std::make_unique<TCanvas>("Canvas", "Canvas", 0, 0, 1024, 1024);
440 Canvas->SetGrid();
441 gStyle->SetOptStat(0);
442 gStyle->SetOptTitle(0);
443 Canvas->SetTickx();
444 Canvas->SetTicky();
445 Canvas->SetBottomMargin(0.1f);
446 Canvas->SetTopMargin(0.05f);
447 Canvas->SetRightMargin(0.15f);
448 Canvas->SetLeftMargin(0.10f);
449
450 //KS: Fancy colours
451 const int NRGBs = 10;
452 TColor::InitializeColors();
453 Double_t stops[NRGBs] = { 0.00, 0.10, 0.25, 0.35, 0.50, 0.60, 0.65, 0.75, 0.90, 1.00 };
454 Double_t red[NRGBs] = { 0.50, 1.00, 1.00, 0.25, 0.00, 0.10, 0.50, 1.00, 0.75, 0.55 };
455 Double_t green[NRGBs] = { 0.00, 0.25, 1.00, 0.25, 0.00, 0.60, 0.90, 1.00, 0.75, 0.75 };
456 Double_t blue[NRGBs] = { 0.00, 0.25, 1.00, 1.00, 0.50, 0.60, 0.90, 1.00, 0.05, 0.05 };
457 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, 255);
458 gStyle->SetNumberContours(255);
459
460 std::string OutName = inputFile;
461 OutName = OutName.substr(0, OutName.find(".root"));
462 Canvas->Print(Form("Correlation_%s.pdf[", OutName.c_str()), "pdf");
463 Canvas->Print(Form("Covariance_%s.pdf[", OutName.c_str()), "pdf");
464
465 YAML::Node card_yaml = M3OpenConfig(config.c_str());
466 YAML::Node Settings = card_yaml["ProcessMCMC"];
467
468 const int entries = int(Processor->GetnSteps());
469 const int NIntervals = GetFromManager<int>(Settings["NIntervals"], 5);
470 const int IntervalsSize = entries/NIntervals;
471 //We start with burn from 0 (no burn in at all)
472 int BurnIn = 0;
473 MACH3LOG_INFO("Diagnosing matrices with entries={}, NIntervals={} and IntervalsSize={}", entries, NIntervals, IntervalsSize);
474
475 TMatrixDSym *Covariance = nullptr;
476 TMatrixDSym *Correlation = nullptr;
477
478 TH2D *CovariancePreviousHist = nullptr;
479 TH2D *CorrelationPreviousHist = nullptr;
480
481 TH2D *CovarianceHist = nullptr;
482 TH2D *CorrelationHist = nullptr;
483
484 //KS: Get first covariances, we need two for comparison...
485 Processor->SetStepCut(BurnIn);
486 Processor->GetCovariance(Covariance, Correlation);
487
488 CovariancePreviousHist = TMatrixIntoTH2D(Covariance, "Covariance");
489 CorrelationPreviousHist = TMatrixIntoTH2D(Correlation, "Correlation");
490
491 delete Covariance;
492 Covariance = nullptr;
493 delete Correlation;
494 Correlation = nullptr;
495
496 //KS: Loop over all desired cuts
497 for(int k = 1; k < NIntervals; ++k)
498 {
499 BurnIn = k*IntervalsSize;
500 Processor->SetStepCut(BurnIn);
501 Processor->GetCovariance(Covariance, Correlation);
502 Processor->ResetHistograms();
503
504 CovarianceHist = TMatrixIntoTH2D(Covariance, "Covariance");
505 CorrelationHist = TMatrixIntoTH2D(Correlation, "Correlation");
506
507 TH2D *CovarianceDiff = static_cast<TH2D*>(CovarianceHist->Clone("Covariance_Ratio"));
508 TH2D *CorrelationDiff = static_cast<TH2D*>(CorrelationHist->Clone("Correlation_Ratio"));
509
510 //KS: Bit messy but quite often covariance is 0 is divided by 0 is problematic so
511 #ifdef MULTITHREAD
512 #pragma omp parallel for
513 #endif
514 for (int j = 1; j < CovarianceDiff->GetXaxis()->GetNbins()+1; ++j)
515 {
516 for (int i = 1; i < CovarianceDiff->GetYaxis()->GetNbins()+1; ++i)
517 {
518 if( std::fabs (CovarianceDiff->GetBinContent(j, i)) < 1.e-5 && std::fabs (CovariancePreviousHist->GetBinContent(j, i)) < 1.e-5)
519 {
520 CovarianceDiff->SetBinContent(j, i, _UNDEF_);
521 CovariancePreviousHist->SetBinContent(j, i, _UNDEF_);
522 }
523 if( std::fabs (CorrelationDiff->GetBinContent(j, i)) < 1.e-5 && std::fabs (CorrelationPreviousHist->GetBinContent(j, i)) < 1.e-5)
524 {
525 CorrelationDiff->SetBinContent(j, i, _UNDEF_);
526 CorrelationPreviousHist->SetBinContent(j, i, _UNDEF_);
527 }
528 }
529 }
530 //Divide matrices
531 CovarianceDiff->Divide(CovariancePreviousHist);
532 CorrelationDiff->Divide(CorrelationPreviousHist);
533
534 //Now it is time for fancy names etc.
535 for (int j = 0; j < CovarianceDiff->GetXaxis()->GetNbins(); ++j)
536 {
537 TString Title = "";
538 double Prior = 1.0;
539 double PriorError = 1.0;
540
541 Processor->GetNthParameter(j, Prior, PriorError, Title);
542
543 CovarianceDiff->GetXaxis()->SetBinLabel(j+1, Title);
544 CovarianceDiff->GetYaxis()->SetBinLabel(j+1, Title);
545 CorrelationDiff->GetXaxis()->SetBinLabel(j+1, Title);
546 CorrelationDiff->GetYaxis()->SetBinLabel(j+1, Title);
547 }
548 CovarianceDiff->GetXaxis()->SetLabelSize(0.015f);
549 CovarianceDiff->GetYaxis()->SetLabelSize(0.015f);
550 CorrelationDiff->GetXaxis()->SetLabelSize(0.015f);
551 CorrelationDiff->GetYaxis()->SetLabelSize(0.015f);
552
553 std::stringstream ss;
554 ss << "BCut_";
555 ss << BurnIn;
556 ss << "/";
557 ss << "BCut_";
558 ss << (k-1)*IntervalsSize;
559 std::string str = ss.str();
560
561 TString Title = "Cov " + str;
562 CovarianceDiff->GetZaxis()->SetTitle( Title );
563 Title = "Corr " + str;
564 CorrelationDiff->GetZaxis()->SetTitle(Title);
565
566 CovarianceDiff->SetMinimum(-2);
567 CovarianceDiff->SetMaximum(2);
568 CorrelationDiff->SetMinimum(-2);
569 CorrelationDiff->SetMaximum(2);
570
571 Canvas->cd();
572 CovarianceDiff->Draw("colz");
573 Canvas->Print(Form("Covariance_%s.pdf", OutName.c_str()), "pdf");
574
575 CorrelationDiff->Draw("colz");
576 Canvas->Print(Form("Correlation_%s.pdf", OutName.c_str()), "pdf");
577
578 //KS: Current hist become previous as we need it for further comparison
579 delete CovariancePreviousHist;
580 CovariancePreviousHist = static_cast<TH2D*>(CovarianceHist->Clone());
581 delete CorrelationPreviousHist;
582 CorrelationPreviousHist = static_cast<TH2D*>(CorrelationHist->Clone());
583
584 delete CovarianceHist;
585 CovarianceHist = nullptr;
586 delete CorrelationHist;
587 CorrelationHist = nullptr;
588
589 delete CovarianceDiff;
590 delete CorrelationDiff;
591 delete Covariance;
592 Covariance = nullptr;
593 delete Correlation;
594 Correlation = nullptr;
595 }
596 Canvas->cd();
597 Canvas->Print(Form("Covariance_%s.pdf]", OutName.c_str()), "pdf");
598 Canvas->Print(Form("Correlation_%s.pdf]", OutName.c_str()), "pdf");
599
600 Processor->SetPrintToPDF(true);
601 if(Covariance != nullptr) delete Covariance;
602 if(Correlation != nullptr) delete Correlation;
603 if(CovariancePreviousHist != nullptr) delete CovariancePreviousHist;
604 if(CorrelationPreviousHist != nullptr) delete CorrelationPreviousHist;
605 if(CovarianceHist != nullptr) delete CovarianceHist;
606 if(CorrelationHist != nullptr) delete CorrelationHist;
607}
#define _UNDEF_
Definition: MCMCProcessor.h:4
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
TH2D * TMatrixIntoTH2D(TMatrixDSym *Matrix, const std::string &title)
KS: Convert TMatrix to TH2D, mostly useful for making fancy plots.
void GetNthParameter(const int param, double &Prior, double &PriorError, TString &Title) const
Get properties of parameter by passing it number.
void ResetHistograms()
Reset 2D posteriors, in case we would like to calculate in again with different BurnInCut.
void SetPrintToPDF(const bool PlotOrNot)
Long64_t GetnSteps()
Get Number of Steps that Chain has, for merged chains will not be the same nEntries.
void GetCovariance(TMatrixDSym *&Cov, TMatrixDSym *&Corr)
Get the post-fit covariances and correlations.
void SetStepCut(const std::string &Cuts)
Set the step cutting by string.

◆ GetTrianglePlot()

void GetTrianglePlot ( MCMCProcessor Processor)
inline

Definition at line 414 of file ProcessMCMC.cpp.

414 {
415 YAML::Node card_yaml = M3OpenConfig(config.c_str());
416 YAML::Node Settings = card_yaml["ProcessMCMC"];
417
418 for (const auto& dg : Settings["TrianglePlot"])
419 {
420 std::string ParName = dg[0].as<std::string>();
421
422 std::vector<std::string> NameVec = dg[1].as<std::vector<std::string>>();
423 Processor->MakeTrianglePlot(NameVec,
424 GetFromManager<std::vector<double>>(Settings["CredibleIntervals"], {0.99, 0.90, 0.68}),
425 GetFromManager<std::vector<short int>>(Settings["CredibleIntervalsColours"], {436, 430, 422}),
426 GetFromManager<std::vector<double>>(Settings["CredibleRegions"], {0.99, 0.90, 0.68}),
427 GetFromManager<std::vector<short int>>(Settings["CredibleRegionStyle"], {2, 1, 3}),
428 GetFromManager<std::vector<short int>>(Settings["CredibleRegionColor"], {413, 406, 416}),
429 GetFromManager<bool>(Settings["CredibleInSigmas"], false));
430 }
431}
Type GetFromManager(const YAML::Node &node, Type defval, const std::string File="", const int Line=1)
Get content of config file if node is not found take default value specified.
Definition: YamlHelper.h:298
void MakeTrianglePlot(const std::vector< std::string > &ParNames, const std::vector< double > &CredibleIntervals={0.99, 0.90, 0.68 }, const std::vector< Color_t > &CredibleIntervalsColours={kCyan+4, kCyan-2, kCyan-10}, const std::vector< double > &CredibleRegions={0.99, 0.90, 0.68}, const std::vector< Style_t > &CredibleRegionStyle={kDashed, kSolid, kDotted}, const std::vector< Color_t > &CredibleRegionColor={kGreen-3, kGreen-10, kGreen}, const bool CredibleInSigmas=false)
Make fancy triangle plot for selected parameters.

◆ 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.

641{
642 const Color_t CumulativeColor[] = {kBlue-1, kRed, kGreen+2};
643 const Style_t CumulativeStyle[] = {kSolid, kDashed, kDotted};
644
645 for(int i = 0; i < Processor[0]->GetNParams(); ++i)
646 {
647 // This holds the posterior density
648 std::vector<TH1D*> hpost(nFiles);
649 std::vector<TH1D*> CumulativeDistribution(nFiles);
650
651 TString Title;
652 double Prior = 1.0;
653 double PriorError = 1.0;
654
655 Processor[0]->GetNthParameter(i, Prior, PriorError, Title);
656 bool Skip = false;
657 for (int ik = 0 ; ik < nFiles; ik++)
658 {
659 int Index = 0;
660 if(ik == 0 ) Index = i;
661 else
662 {
663 // KS: If somehow this chain doesn't given params we skip it
664 Index = Processor[ik]->GetParamIndexFromName(hpost[0]->GetTitle());
665 if(Index == _UNDEF_)
666 {
667 Skip = true;
668 break;
669 }
670 }
671 hpost[ik] = static_cast<TH1D*>(Processor[ik]->GetHpost(Index)->Clone());
672 CumulativeDistribution[ik] = static_cast<TH1D*>(Processor[ik]->GetHpost(Index)->Clone());
673 CumulativeDistribution[ik]->Fill(0., 0.);
674 CumulativeDistribution[ik]->Reset();
675 CumulativeDistribution[ik]->SetMaximum(1.);
676 TString TempTittle = Title+" Kolmogorov Smirnov";
677 CumulativeDistribution[ik]->SetTitle(TempTittle);
678
679 TempTittle = Title+" Value";
680 CumulativeDistribution[ik]->GetXaxis()->SetTitle(TempTittle);
681 CumulativeDistribution[ik]->GetYaxis()->SetTitle("Cumulative Probability");
682
683 CumulativeDistribution[ik]->SetLineWidth(2);
684 CumulativeDistribution[ik]->SetLineColor(CumulativeColor[ik]);
685 CumulativeDistribution[ik]->SetLineStyle(CumulativeStyle[ik]);
686 }
687
688 // Don't plot if this is a fixed histogram (i.e. the peak is the whole integral)
689 if(hpost[0]->GetMaximum() == hpost[0]->Integral()*1.5 || Skip)
690 {
691 for (int ik = 0; ik < nFiles; ik++)
692 {
693 delete hpost[ik];
694 delete CumulativeDistribution[ik];
695 }
696 continue;
697 }
698
699 for (int ik = 0 ; ik < nFiles; ik++)
700 {
701 const int NumberOfBins = hpost[ik]->GetXaxis()->GetNbins();
702 double Cumulative = 0;
703 const double Integral = hpost[ik]->Integral();
704 for (int j = 1; j < NumberOfBins+1; ++j)
705 {
706 Cumulative += hpost[ik]->GetBinContent(j)/Integral;
707 CumulativeDistribution[ik]->SetBinContent(j, Cumulative);
708 }
709 //KS: Set overflow to 1 just in case
710 CumulativeDistribution[ik]->SetBinContent(NumberOfBins+1, 1.);
711 }
712
713 std::vector<int> TestStatBin(nFiles, 0);
714 std::vector<double> TestStatD(nFiles, -999);
715 std::vector<std::unique_ptr<TLine>> LineD(nFiles);
716 //Find KS statistic
717 for (int ik = 1 ; ik < nFiles; ik++)
718 {
719 const int NumberOfBins = CumulativeDistribution[0]->GetXaxis()->GetNbins();
720 for (int j = 1; j < NumberOfBins+1; ++j)
721 {
722 const double BinValue = CumulativeDistribution[0]->GetBinCenter(j);
723 const int BinNumber = CumulativeDistribution[ik]->FindBin(BinValue);
724 //KS: Calculate D statistic for this bin, only save it if it's bigger than previously found value
725 double TempDstat = std::fabs(CumulativeDistribution[0]->GetBinContent(j) - CumulativeDistribution[ik]->GetBinContent(BinNumber));
726 if(TempDstat > TestStatD[ik])
727 {
728 TestStatD[ik] = TempDstat;
729 TestStatBin[ik] = j;
730 }
731 }
732 }
733
734 for (int ik = 0 ; ik < nFiles; ik++)
735 {
736 LineD[ik] = std::make_unique<TLine>(CumulativeDistribution[0]->GetBinCenter(TestStatBin[ik]), 0, CumulativeDistribution[0]->GetBinCenter(TestStatBin[ik]), CumulativeDistribution[0]->GetBinContent(TestStatBin[ik]));
737 LineD[ik]->SetLineColor(CumulativeColor[ik]);
738 LineD[ik]->SetLineWidth(2.0);
739 }
740 CumulativeDistribution[0]->Draw();
741 for (int ik = 0 ; ik < nFiles; ik++)
742 CumulativeDistribution[ik]->Draw("SAME");
743
744 auto leg = std::make_unique<TLegend>(0.15, 0.7, 0.5, 0.90);
745 leg->SetTextSize(0.04f);
746 for (int ik = 0; ik < nFiles; ik++)
747 leg->AddEntry(CumulativeDistribution[ik], TitleNames[ik].c_str(), "l");
748 for (int ik = 1; ik < nFiles; ik++)
749 leg->AddEntry(LineD[ik].get(), Form("#Delta D = %.4f", TestStatD[ik]), "l");
750
751 leg->SetLineColor(0);
752 leg->SetLineStyle(0);
753 leg->SetFillColor(0);
754 leg->SetFillStyle(0);
755 leg->Draw("SAME");
756
757 for (int ik = 1; ik < nFiles; ik++)
758 LineD[ik]->Draw("sam");
759
760 Posterior->cd();
761 Posterior->Print(canvasname);
762
763 for (int ik = 0; ik < nFiles; ik++)
764 {
765 delete hpost[ik];
766 delete CumulativeDistribution[ik];
767 }
768 } //End loop over parameter
769}
std::vector< std::string > TitleNames
Definition: ProcessMCMC.cpp:29
int nFiles
Definition: ProcessMCMC.cpp:27

◆ main()

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

Definition at line 32 of file ProcessMCMC.cpp.

33{
35 nFiles = 0;
36 if (argc != 3 && argc !=6 && argc != 8)
37 {
38 MACH3LOG_ERROR("How to use: {} <Config> <MCMM_ND_Output.root>", argv[0]);
39 throw MaCh3Exception(__FILE__ , __LINE__ );
40 }
41
42 if (argc == 3)
43 {
44 MACH3LOG_INFO("Producing single fit output");
45 config = argv[1];
46 std::string filename = argv[2];
47 ProcessMCMC(filename);
48 }
49 // If we want to compare two or more fits (e.g. binning changes or introducing new params/priors)
50 else if (argc == 6 || argc == 8)
51 {
52 MACH3LOG_INFO("Producing two fit comparison");
53 config = argv[1];
54
55 FileNames.push_back(argv[2]);
56 TitleNames.push_back(argv[3]);
57
58 FileNames.push_back(argv[4]);
59 TitleNames.push_back(argv[5]);
60 //KS: If there is third file add it
61 if(argc == 8)
62 {
63 FileNames.push_back(argv[6]);
64 TitleNames.push_back(argv[7]);
65 }
66
68 }
69 return 0;
70}
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
void SetMaCh3LoggerFormat()
Set messaging format of the logger.
Definition: MaCh3Logger.h:30
void ProcessMCMC(const std::string &inputFile)
Main function processing MCMC and Producing plots.
Definition: ProcessMCMC.cpp:72
void MultipleProcessMCMC()
Function producing comparison of posterior and more betwen a few MCMC chains.
std::vector< std::string > FileNames
Definition: ProcessMCMC.cpp:28
Custom exception class for MaCh3 errors.

◆ MultipleProcessMCMC()

void MultipleProcessMCMC ( )
inline

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

Definition at line 167 of file ProcessMCMC.cpp.

168{
169 YAML::Node card_yaml = M3OpenConfig(config.c_str());
170 YAML::Node Settings = card_yaml["ProcessMCMC"];
171
172 constexpr Color_t PosteriorColor[] = {kBlue-1, kRed, kGreen+2};
173 //constexpr Style_t PosteriorStyle[] = {kSolid, kDashed, kDotted};
174 nFiles = int(FileNames.size());
175 std::vector<std::unique_ptr<MCMCProcessor>> Processor(nFiles);
176
177 if(!Settings["BurnInSteps"])
178 {
179 MACH3LOG_WARN("BurnInSteps not set, defaulting to 20%");
180 }
181
182 for (int ik = 0; ik < nFiles; ik++)
183 {
184 MACH3LOG_INFO("File for study: {}", FileNames[ik]);
185 // Make the processor
186 Processor[ik] = std::make_unique<MCMCProcessor>(FileNames[ik]);
187 Processor[ik]->SetOutputSuffix(("_" + std::to_string(ik)).c_str());
188
189 Processor[ik]->SetExcludedTypes(GetFromManager<std::vector<std::string>>(Settings["ExcludedTypes"], {}));
190 Processor[ik]->SetExcludedNames(GetFromManager<std::vector<std::string>>(Settings["ExcludedNames"], {}));
191
192 //Apply additional cuts to 1D posterior
193 Processor[ik]->SetPosterior1DCut(GetFromManager<std::string>(Settings["Posterior1DCut"], ""));
194
195 Processor[ik]->SetPlotRelativeToPrior(GetFromManager<bool>(Settings["PlotRelativeToPrior"], false));
196 Processor[ik]->SetFancyNames(GetFromManager<bool>(Settings["FancyNames"], true));
197 Processor[ik]->Initialise();
198
199 if(Settings["BurnInSteps"])
200 {
201 Processor[ik]->SetStepCut(Settings["BurnInSteps"].as<int>());
202 }
203 else
204 {
205 Processor[ik]->SetStepCut(static_cast<int>(Processor[ik]->GetnSteps()/5));
206 }
207 }
208 //KS: Multithreading here is very tempting but there are some issues with root that need to be resovled :(
209 for (int ik = 0; ik < nFiles; ik++)
210 {
211 // Make the postfit
212 Processor[ik]->MakePostfit();
213 Processor[ik]->DrawPostfit();
214 }
215
216 // Open a TCanvas to write the posterior onto
217 auto Posterior = std::make_unique<TCanvas>("PosteriorMulti", "PosteriorMulti", 0, 0, 1024, 1024);
218 gStyle->SetOptStat(0);
219 gStyle->SetOptTitle(0);
220 Posterior->SetGrid();
221 Posterior->SetBottomMargin(0.1f);
222 Posterior->SetTopMargin(0.05f);
223 Posterior->SetRightMargin(0.03f);
224 Posterior->SetLeftMargin(0.10f);
225
226 FileNames[0] = FileNames[0].substr(0, FileNames[0].find(".root")-1);
227 TString canvasname = FileNames[0];
228 for (int ik = 1; ik < nFiles; ik++)
229 {
230 while (FileNames[ik].find("/") != std::string::npos)
231 {
232 FileNames[ik] = FileNames[ik].substr(FileNames[ik].find("/")+1, FileNames[ik].find(".root")-1);
233 }
234 canvasname = canvasname + "_"+FileNames[ik];
235 }
236
237 canvasname = canvasname +".pdf[";
238
239 Posterior->Print(canvasname);
240 // Once the pdf file is open no longer need to bracket
241 canvasname.ReplaceAll("[","");
242
243 for(int i = 0; i < Processor[0]->GetNParams(); ++i)
244 {
245 // This holds the posterior density
246 std::vector<TH1D*> hpost(nFiles);
247 std::vector<std::unique_ptr<TLine>> hpd(nFiles);
248 hpost[0] = static_cast<TH1D *>(Processor[0]->GetHpost(i)->Clone());
249
250 bool Skip = false;
251 for (int ik = 1 ; ik < nFiles; ik++)
252 {
253 // KS: If somehow this chain doesn't given params we skip it
254 const int Index = Processor[ik]->GetParamIndexFromName(hpost[0]->GetTitle());
255 if(Index == _UNDEF_)
256 {
257 Skip = true;
258 break;
259 }
260 hpost[ik] = static_cast<TH1D *>(Processor[ik]->GetHpost(Index)->Clone());
261 }
262
263 // Don't plot if this is a fixed histogram (i.e. the peak is the whole integral)
264 if(hpost[0]->GetMaximum() == hpost[0]->Integral()*1.5 || Skip)
265 {
266 for (int ik = 0; ik < nFiles; ik++)
267 delete hpost[ik];
268
269 continue;
270 }
271 for (int ik = 0; ik < nFiles; ik++)
272 {
273 RemoveFitter(hpost[ik], "Gauss");
274
275 // Set some nice colours
276 hpost[ik]->SetLineColor(PosteriorColor[ik]);
277 //hpost[ik]->SetLineStyle(PosteriorStyle[ik]);
278 hpost[ik]->SetLineWidth(2);
279
280 // Area normalise the distributions
281 hpost[ik]->Scale(1./hpost[ik]->Integral(), "width");
282 }
283 TString Title;
284 double Prior = 1.0;
285 double PriorError = 1.0;
286
287 Processor[0]->GetNthParameter(i, Prior, PriorError, Title);
288
289 // Now make the TLine for the Asimov
290 auto Asimov = std::make_unique<TLine>(Prior, hpost[0]->GetMinimum(), Prior, hpost[0]->GetMaximum());
291 Asimov->SetLineColor(kRed-3);
292 Asimov->SetLineWidth(2);
293 Asimov->SetLineStyle(kDashed);
294
295 // Make a nice little TLegend
296 auto leg = std::make_unique<TLegend>(0.12, 0.7, 0.6, 0.97);
297 leg->SetTextSize(0.03f);
298 leg->SetFillColor(0);
299 leg->SetFillStyle(0);
300 leg->SetLineColor(0);
301 leg->SetLineStyle(0);
302 TString asimovLeg = Form("#splitline{Prior}{x = %.2f , #sigma = %.2f}", Prior, PriorError);
303 leg->AddEntry(Asimov.get(), asimovLeg, "l");
304
305 for (int ik = 0; ik < nFiles; ik++)
306 {
307 TString rebinLeg = Form("#splitline{%s}{#mu = %.2f, #sigma = %.2f}", TitleNames[ik].c_str(), hpost[ik]->GetMean(), hpost[ik]->GetRMS());
308 leg->AddEntry(hpost[ik], rebinLeg, "l");
309
310 hpd[ik] = std::make_unique<TLine>(hpost[ik]->GetBinCenter(hpost[ik]->GetMaximumBin()), hpost[ik]->GetMinimum(),
311 hpost[ik]->GetBinCenter(hpost[ik]->GetMaximumBin()), hpost[ik]->GetMaximum());
312 hpd[ik]->SetLineColor(hpost[ik]->GetLineColor());
313 hpd[ik]->SetLineWidth(2);
314 hpd[ik]->SetLineStyle(kSolid);
315 }
316
317 // Find the maximum value to nicely resize hist
318 double maximum = 0;
319 for (int ik = 0; ik < nFiles; ik++) maximum = std::max(maximum, hpost[ik]->GetMaximum());
320 for (int ik = 0; ik < nFiles; ik++) hpost[ik]->SetMaximum(1.3*maximum);
321
322 hpost[0]->Draw("hist");
323 for (int ik = 1; ik < nFiles; ik++) hpost[ik]->Draw("hist same");
324 Asimov->Draw("same");
325 for (int ik = 0; ik < nFiles; ik++) hpd[ik]->Draw("same");
326 leg->Draw("same");
327 Posterior->cd();
328 Posterior->Print(canvasname);
329 for (int ik = 0; ik < nFiles; ik++) {
330 delete hpost[ik];
331 }
332 }//End loop over parameters
333
334 // Finally draw the parameter plot onto the PDF
335 // Close the .pdf file with all the posteriors
336 Posterior->cd();
337 Posterior->Clear();
338
339 if(GetFromManager<bool>(Settings["PerformKStest"], true)) KolmogorovSmirnovTest(Processor, Posterior, canvasname);
340
341 // Close the pdf file
342 MACH3LOG_INFO("Closing pdf {}", canvasname);
343 canvasname+="]";
344 Posterior->Print(canvasname);
345}
void RemoveFitter(TH1D *hist, const std::string &name)
KS: Remove fitted TF1 from hist to make comparison easier.
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
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...

◆ ProcessMCMC()

void ProcessMCMC ( const std::string &  inputFile)
inline

Main function processing MCMC and Producing plots.

Definition at line 72 of file ProcessMCMC.cpp.

73{
74 MACH3LOG_INFO("File for study: {} with config {}", inputFile, config);
75 // Make the processor)
76 auto Processor = std::make_unique<OscProcessor>(inputFile);
77
78 YAML::Node card_yaml = M3OpenConfig(config.c_str());
79 YAML::Node Settings = card_yaml["ProcessMCMC"];
80
81 const bool PlotCorr = GetFromManager<bool>(Settings["PlotCorr"], false);
82
83 Processor->SetExcludedTypes(GetFromManager<std::vector<std::string>>(Settings["ExcludedTypes"], {}));
84 Processor->SetExcludedNames(GetFromManager<std::vector<std::string>>(Settings["ExcludedNames"], {}));
85 //Apply additional cuts to 1D posterior
86 Processor->SetPosterior1DCut(GetFromManager<std::string>(Settings["Posterior1DCut"], ""));
87
88 if(PlotCorr) Processor->SetOutputSuffix("_drawCorr");
89 //KS:Turn off plotting detector and some other setting, should be via some config
90 Processor->SetPlotRelativeToPrior(GetFromManager<bool>(Settings["PlotRelativeToPrior"], false));
91 Processor->SetPrintToPDF(GetFromManager<bool>(Settings["PrintToPDF"], true));
92
93 //KS: Whether you want prior error bands for parameters with flat prior or not
94 Processor->SetPlotErrorForFlatPrior(GetFromManager<bool>(Settings["PlotErrorForFlatPrior"], true));
95 Processor->SetFancyNames(GetFromManager<bool>(Settings["FancyNames"], true));
96 Processor->SetPlotBinValue(GetFromManager<bool>(Settings["PlotBinValue"], false));
97 //KS: Plot only 2D posteriors with correlations greater than 0.2
98 Processor->SetPost2DPlotThreshold(GetFromManager<double>(Settings["Post2DPlotThreshold"], 0.2));
99
100 Processor->Initialise();
101
102 if(Settings["BurnInSteps"])
103 {
104 Processor->SetStepCut(Settings["BurnInSteps"].as<int>());
105 }
106 else
107 {
108 MACH3LOG_WARN("BurnInSteps not set, defaulting to 20%");
109 Processor->SetStepCut(static_cast<int>(Processor->GetnSteps()/5));
110 }
111
112 if(Settings["Thinning"])
113 {
114 if(Settings["Thinning"][0].as<bool>()){
115 Processor->ThinMCMC(Settings["Thinning"][1].as<int>());
116 }
117 }
118 // Make the postfit
119 Processor->MakePostfit();
120 Processor->DrawPostfit();
121 //KS: Should set via config whether you want below or not
122 if(GetFromManager<bool>(Settings["MakeCredibleIntervals"], true)) {
123 Processor->MakeCredibleIntervals(GetFromManager<std::vector<double>>(Settings["CredibleIntervals"], {0.99, 0.90, 0.68}),
124 GetFromManager<std::vector<short int>>(Settings["CredibleIntervalsColours"], {436, 430, 422}),
125 GetFromManager<bool>(Settings["CredibleInSigmas"], false));
126 }
127 if(GetFromManager<bool>(Settings["CalcBayesFactor"], true)) CalcBayesFactor(Processor.get());
128 if(GetFromManager<bool>(Settings["CalcSavageDickey"], true)) CalcSavageDickey(Processor.get());
129 if(GetFromManager<bool>(Settings["CalcBipolarPlot"], false)) CalcBipolarPlot(Processor.get());
130 if(GetFromManager<bool>(Settings["CalcParameterEvolution"], false)) CalcParameterEvolution(Processor.get());
131
132 if(PlotCorr)
133 {
134 Processor->SetSmoothing(GetFromManager<bool>(Settings["Smoothing"], true));
135 // Make the covariance matrix
136 //We have different treatment for multithread
137//#ifdef MULTITHREAD
138 Processor->CacheSteps();
139 //KS: Since we cached let's make fancy violins :)
140 if(GetFromManager<bool>(Settings["MakeViolin"], true)) Processor->MakeViolin();
141 Processor->MakeCovariance_MP();
142//#else
143 //Processor->MakeCovariance();
144//#endif
145 Processor->DrawCovariance();
146
147 auto const &MakeSubOptimality = Settings["MakeSubOptimality"];
148 if(MakeSubOptimality[0].as<bool>()) Processor->MakeSubOptimality(MakeSubOptimality[1].as<int>());
149
150 if(GetFromManager<bool>(Settings["MakeCredibleRegions"], false)) {
151 Processor->MakeCredibleRegions(GetFromManager<std::vector<double>>(Settings["CredibleRegions"], {0.99, 0.90, 0.68}),
152 GetFromManager<std::vector<short int>>(Settings["CredibleRegionStyle"], {2, 1, 3}),
153 GetFromManager<std::vector<short int>>(Settings["CredibleRegionColor"], {413, 406, 416}),
154 GetFromManager<bool>(Settings["CredibleInSigmas"], false),
155 GetFromManager<bool>(Settings["Draw2DPosterior"], true),
156 GetFromManager<bool>(Settings["DrawBestFit"], true));
157 }
158 if(GetFromManager<bool>(Settings["GetTrianglePlot"], true)) GetTrianglePlot(Processor.get());
159
160 //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
161 if(GetFromManager<bool>(Settings["DiagnoseCovarianceMatrix"], false)) DiagnoseCovarianceMatrix(Processor.get(), inputFile);
162 }
163 if(GetFromManager<bool>(Settings["ReweightPrior"], false)) ReweightPrior(Processor.get());
164 if(GetFromManager<bool>(Settings["JarlskogAnalysis"], true)) Processor->PerformJarlskogAnalysis();
165}
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...
void CalcBipolarPlot(MCMCProcessor *Processor)
void CalcParameterEvolution(MCMCProcessor *Processor)
void ReweightPrior(MCMCProcessor *Processor)
void CalcBayesFactor(MCMCProcessor *Processor)
void CalcSavageDickey(MCMCProcessor *Processor)

◆ ReweightPrior()

void ReweightPrior ( MCMCProcessor Processor)
inline

Definition at line 609 of file ProcessMCMC.cpp.

610{
611 YAML::Node card_yaml = M3OpenConfig(config.c_str());
612 YAML::Node Settings = card_yaml["ProcessMCMC"];
613
614 const auto& Prior = Settings["PriorReweighting"];
615 std::vector<std::string> Names = Prior[0].as<std::vector<std::string>>();
616 std::vector<double> NewCentral = Prior[1].as<std::vector<double>>();
617 std::vector<double> NewError = Prior[2].as<std::vector<double>>();
618
619 Processor->ReweightPrior(Names, NewCentral, NewError);
620}
void ReweightPrior(const std::vector< std::string > &Names, const std::vector< double > &NewCentral, const std::vector< double > &NewError)
Reweight Prior by giving new central value and new error.

◆ 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.

624{
625 TH2D* hMatrix = new TH2D(title.c_str(), title.c_str(), Matrix->GetNrows(), 0.0, Matrix->GetNrows(), Matrix->GetNcols(), 0.0, Matrix->GetNcols());
626 for(int i = 0; i < Matrix->GetNrows(); i++)
627 {
628 for(int j = 0; j < Matrix->GetNcols(); j++)
629 {
630 //KS: +1 because there is offset in histogram relative to TMatrix
631 hMatrix->SetBinContent(i+1,j+1, (*Matrix)(i,j));
632 }
633 }
634 return hMatrix;
635}

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.