10inline void ProcessMCMC(
const std::string& inputFile);
21inline TH2D*
TMatrixIntoTH2D(TMatrixDSym* Matrix,
const std::string& title);
24 const std::unique_ptr<TCanvas>& Posterior,
25 const TString& canvasname);
32int main(
int argc,
char *argv[])
36 if (argc != 3 && argc !=6 && argc != 8)
38 MACH3LOG_ERROR(
"How to use: {} <Config> <MCMM_ND_Output.root>", argv[0]);
46 std::string filename = argv[2];
50 else if (argc == 6 || argc == 8)
76 auto Processor = std::make_unique<OscProcessor>(inputFile);
79 YAML::Node Settings = card_yaml[
"ProcessMCMC"];
81 const bool PlotCorr = GetFromManager<bool>(Settings[
"PlotCorr"],
false);
83 Processor->SetExcludedTypes(
GetFromManager<std::vector<std::string>>(Settings[
"ExcludedTypes"], {}));
84 Processor->SetExcludedNames(
GetFromManager<std::vector<std::string>>(Settings[
"ExcludedNames"], {}));
86 Processor->SetPosterior1DCut(GetFromManager<std::string>(Settings[
"Posterior1DCut"],
""));
88 if(PlotCorr) Processor->SetOutputSuffix(
"_drawCorr");
90 Processor->SetPlotRelativeToPrior(GetFromManager<bool>(Settings[
"PlotRelativeToPrior"],
false));
91 Processor->SetPrintToPDF(GetFromManager<bool>(Settings[
"PrintToPDF"],
true));
94 Processor->SetPlotErrorForFlatPrior(GetFromManager<bool>(Settings[
"PlotErrorForFlatPrior"],
true));
95 Processor->SetFancyNames(GetFromManager<bool>(Settings[
"FancyNames"],
true));
96 Processor->SetPlotBinValue(GetFromManager<bool>(Settings[
"PlotBinValue"],
false));
98 Processor->SetPost2DPlotThreshold(GetFromManager<double>(Settings[
"Post2DPlotThreshold"], 0.2));
100 Processor->Initialise();
102 if(Settings[
"BurnInSteps"])
104 Processor->SetStepCut(Settings[
"BurnInSteps"].as<int>());
109 Processor->SetStepCut(
static_cast<int>(Processor->GetnSteps()/5));
112 if(Settings[
"Thinning"])
114 if(Settings[
"Thinning"][0].as<bool>()){
115 Processor->ThinMCMC(Settings[
"Thinning"][1].as<int>());
119 Processor->MakePostfit();
120 Processor->DrawPostfit();
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));
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());
134 Processor->SetSmoothing(GetFromManager<bool>(Settings[
"Smoothing"],
true));
138 Processor->CacheSteps();
140 if(GetFromManager<bool>(Settings[
"MakeViolin"],
true)) Processor->MakeViolin();
141 Processor->MakeCovariance_MP();
145 Processor->DrawCovariance();
147 auto const &MakeSubOptimality = Settings[
"MakeSubOptimality"];
148 if(MakeSubOptimality[0].as<bool>()) Processor->MakeSubOptimality(MakeSubOptimality[1].as<int>());
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));
158 if(GetFromManager<bool>(Settings[
"GetTrianglePlot"],
true))
GetTrianglePlot(Processor.get());
161 if(GetFromManager<bool>(Settings[
"DiagnoseCovarianceMatrix"],
false))
DiagnoseCovarianceMatrix(Processor.get(), inputFile);
163 if(GetFromManager<bool>(Settings[
"ReweightPrior"],
false))
ReweightPrior(Processor.get());
164 if(GetFromManager<bool>(Settings[
"JarlskogAnalysis"],
true)) Processor->PerformJarlskogAnalysis();
170 YAML::Node Settings = card_yaml[
"ProcessMCMC"];
172 constexpr Color_t PosteriorColor[] = {kBlue-1, kRed, kGreen+2};
175 std::vector<std::unique_ptr<MCMCProcessor>> Processor(
nFiles);
177 if(!Settings[
"BurnInSteps"])
182 for (
int ik = 0; ik <
nFiles; ik++)
186 Processor[ik] = std::make_unique<MCMCProcessor>(
FileNames[ik]);
187 Processor[ik]->SetOutputSuffix((
"_" + std::to_string(ik)).c_str());
189 Processor[ik]->SetExcludedTypes(
GetFromManager<std::vector<std::string>>(Settings[
"ExcludedTypes"], {}));
190 Processor[ik]->SetExcludedNames(
GetFromManager<std::vector<std::string>>(Settings[
"ExcludedNames"], {}));
193 Processor[ik]->SetPosterior1DCut(GetFromManager<std::string>(Settings[
"Posterior1DCut"],
""));
195 Processor[ik]->SetPlotRelativeToPrior(GetFromManager<bool>(Settings[
"PlotRelativeToPrior"],
false));
196 Processor[ik]->SetFancyNames(GetFromManager<bool>(Settings[
"FancyNames"],
true));
197 Processor[ik]->Initialise();
199 if(Settings[
"BurnInSteps"])
201 Processor[ik]->SetStepCut(Settings[
"BurnInSteps"].as<int>());
205 Processor[ik]->SetStepCut(
static_cast<int>(Processor[ik]->GetnSteps()/5));
209 for (
int ik = 0; ik <
nFiles; ik++)
212 Processor[ik]->MakePostfit();
213 Processor[ik]->DrawPostfit();
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);
228 for (
int ik = 1; ik <
nFiles; ik++)
230 while (
FileNames[ik].find(
"/") != std::string::npos)
234 canvasname = canvasname +
"_"+
FileNames[ik];
237 canvasname = canvasname +
".pdf[";
239 Posterior->Print(canvasname);
241 canvasname.ReplaceAll(
"[",
"");
243 for(
int i = 0; i < Processor[0]->GetNParams(); ++i)
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());
251 for (
int ik = 1 ; ik <
nFiles; ik++)
254 const int Index = Processor[ik]->GetParamIndexFromName(hpost[0]->GetTitle());
260 hpost[ik] =
static_cast<TH1D *
>(Processor[ik]->GetHpost(Index)->Clone());
264 if(hpost[0]->GetMaximum() == hpost[0]->Integral()*1.5 || Skip)
266 for (
int ik = 0; ik <
nFiles; ik++)
271 for (
int ik = 0; ik <
nFiles; ik++)
276 hpost[ik]->SetLineColor(PosteriorColor[ik]);
278 hpost[ik]->SetLineWidth(2);
281 hpost[ik]->Scale(1./hpost[ik]->Integral(),
"width");
285 double PriorError = 1.0;
287 Processor[0]->GetNthParameter(i, Prior, PriorError, Title);
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);
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");
305 for (
int ik = 0; ik <
nFiles; ik++)
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");
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);
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);
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");
328 Posterior->Print(canvasname);
329 for (
int ik = 0; ik <
nFiles; ik++) {
339 if(GetFromManager<bool>(Settings[
"PerformKStest"],
true))
KolmogorovSmirnovTest(Processor, Posterior, canvasname);
344 Posterior->Print(canvasname);
351 YAML::Node Settings = card_yaml[
"ProcessMCMC"];
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"])
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>>());
365 Processor->
GetBayesFactor(ParNames, Model1Bounds, Model2Bounds, ModelNames);
371 YAML::Node Settings = card_yaml[
"ProcessMCMC"];
373 std::vector<std::string> ParNames;
374 std::vector<double> EvaluationPoint;
375 std::vector<std::vector<double>> Bounds;
377 for (
const auto& d : Settings[
"SavageDickey"])
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>>());
389 YAML::Node Settings = card_yaml[
"ProcessMCMC"];
391 std::vector<std::string> ParNames;
392 std::vector<int> Intervals;
393 for (
const auto& d : Settings[
"ParameterEvolution"])
395 ParNames.push_back(d[0].as<std::string>());
396 Intervals.push_back(d[1].as<int>());
404 YAML::Node Settings = card_yaml[
"ProcessMCMC"];
406 std::vector<std::string> ParNames;
407 for (
const auto& d : Settings[
"BipolarPlot"])
409 ParNames.push_back(d[0].as<std::string>());
416 YAML::Node Settings = card_yaml[
"ProcessMCMC"];
418 for (
const auto& dg : Settings[
"TrianglePlot"])
420 std::string ParName = dg[0].as<std::string>();
422 std::vector<std::string> NameVec = dg[1].as<std::vector<std::string>>();
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));
439 auto Canvas = std::make_unique<TCanvas>(
"Canvas",
"Canvas", 0, 0, 1024, 1024);
441 gStyle->SetOptStat(0);
442 gStyle->SetOptTitle(0);
445 Canvas->SetBottomMargin(0.1f);
446 Canvas->SetTopMargin(0.05f);
447 Canvas->SetRightMargin(0.15f);
448 Canvas->SetLeftMargin(0.10f);
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);
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");
466 YAML::Node Settings = card_yaml[
"ProcessMCMC"];
468 const int entries = int(Processor->
GetnSteps());
469 const int NIntervals = GetFromManager<int>(Settings[
"NIntervals"], 5);
470 const int IntervalsSize = entries/NIntervals;
473 MACH3LOG_INFO(
"Diagnosing matrices with entries={}, NIntervals={} and IntervalsSize={}", entries, NIntervals, IntervalsSize);
475 TMatrixDSym *Covariance =
nullptr;
476 TMatrixDSym *Correlation =
nullptr;
478 TH2D *CovariancePreviousHist =
nullptr;
479 TH2D *CorrelationPreviousHist =
nullptr;
481 TH2D *CovarianceHist =
nullptr;
482 TH2D *CorrelationHist =
nullptr;
492 Covariance =
nullptr;
494 Correlation =
nullptr;
497 for(
int k = 1; k < NIntervals; ++k)
499 BurnIn = k*IntervalsSize;
507 TH2D *CovarianceDiff =
static_cast<TH2D*
>(CovarianceHist->Clone(
"Covariance_Ratio"));
508 TH2D *CorrelationDiff =
static_cast<TH2D*
>(CorrelationHist->Clone(
"Correlation_Ratio"));
512 #pragma omp parallel for
514 for (
int j = 1; j < CovarianceDiff->GetXaxis()->GetNbins()+1; ++j)
516 for (
int i = 1; i < CovarianceDiff->GetYaxis()->GetNbins()+1; ++i)
518 if( std::fabs (CovarianceDiff->GetBinContent(j, i)) < 1.e-5 && std::fabs (CovariancePreviousHist->GetBinContent(j, i)) < 1.e-5)
520 CovarianceDiff->SetBinContent(j, i,
_UNDEF_);
521 CovariancePreviousHist->SetBinContent(j, i,
_UNDEF_);
523 if( std::fabs (CorrelationDiff->GetBinContent(j, i)) < 1.e-5 && std::fabs (CorrelationPreviousHist->GetBinContent(j, i)) < 1.e-5)
525 CorrelationDiff->SetBinContent(j, i,
_UNDEF_);
526 CorrelationPreviousHist->SetBinContent(j, i,
_UNDEF_);
531 CovarianceDiff->Divide(CovariancePreviousHist);
532 CorrelationDiff->Divide(CorrelationPreviousHist);
535 for (
int j = 0; j < CovarianceDiff->GetXaxis()->GetNbins(); ++j)
539 double PriorError = 1.0;
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);
548 CovarianceDiff->GetXaxis()->SetLabelSize(0.015f);
549 CovarianceDiff->GetYaxis()->SetLabelSize(0.015f);
550 CorrelationDiff->GetXaxis()->SetLabelSize(0.015f);
551 CorrelationDiff->GetYaxis()->SetLabelSize(0.015f);
553 std::stringstream ss;
558 ss << (k-1)*IntervalsSize;
559 std::string str = ss.str();
561 TString Title =
"Cov " + str;
562 CovarianceDiff->GetZaxis()->SetTitle( Title );
563 Title =
"Corr " + str;
564 CorrelationDiff->GetZaxis()->SetTitle(Title);
566 CovarianceDiff->SetMinimum(-2);
567 CovarianceDiff->SetMaximum(2);
568 CorrelationDiff->SetMinimum(-2);
569 CorrelationDiff->SetMaximum(2);
572 CovarianceDiff->Draw(
"colz");
573 Canvas->Print(Form(
"Covariance_%s.pdf", OutName.c_str()),
"pdf");
575 CorrelationDiff->Draw(
"colz");
576 Canvas->Print(Form(
"Correlation_%s.pdf", OutName.c_str()),
"pdf");
579 delete CovariancePreviousHist;
580 CovariancePreviousHist =
static_cast<TH2D*
>(CovarianceHist->Clone());
581 delete CorrelationPreviousHist;
582 CorrelationPreviousHist =
static_cast<TH2D*
>(CorrelationHist->Clone());
584 delete CovarianceHist;
585 CovarianceHist =
nullptr;
586 delete CorrelationHist;
587 CorrelationHist =
nullptr;
589 delete CovarianceDiff;
590 delete CorrelationDiff;
592 Covariance =
nullptr;
594 Correlation =
nullptr;
597 Canvas->Print(Form(
"Covariance_%s.pdf]", OutName.c_str()),
"pdf");
598 Canvas->Print(Form(
"Correlation_%s.pdf]", OutName.c_str()),
"pdf");
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;
612 YAML::Node Settings = card_yaml[
"ProcessMCMC"];
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>>();
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++)
628 for(
int j = 0; j < Matrix->GetNcols(); j++)
631 hMatrix->SetBinContent(i+1,j+1, (*Matrix)(i,j));
639 const std::unique_ptr<TCanvas>& Posterior,
640 const TString& canvasname)
642 const Color_t CumulativeColor[] = {kBlue-1, kRed, kGreen+2};
643 const Style_t CumulativeStyle[] = {kSolid, kDashed, kDotted};
645 for(
int i = 0; i < Processor[0]->GetNParams(); ++i)
648 std::vector<TH1D*> hpost(
nFiles);
649 std::vector<TH1D*> CumulativeDistribution(
nFiles);
653 double PriorError = 1.0;
655 Processor[0]->GetNthParameter(i, Prior, PriorError, Title);
657 for (
int ik = 0 ; ik <
nFiles; ik++)
660 if(ik == 0 ) Index = i;
664 Index = Processor[ik]->GetParamIndexFromName(hpost[0]->GetTitle());
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);
679 TempTittle = Title+
" Value";
680 CumulativeDistribution[ik]->GetXaxis()->SetTitle(TempTittle);
681 CumulativeDistribution[ik]->GetYaxis()->SetTitle(
"Cumulative Probability");
683 CumulativeDistribution[ik]->SetLineWidth(2);
684 CumulativeDistribution[ik]->SetLineColor(CumulativeColor[ik]);
685 CumulativeDistribution[ik]->SetLineStyle(CumulativeStyle[ik]);
689 if(hpost[0]->GetMaximum() == hpost[0]->Integral()*1.5 || Skip)
691 for (
int ik = 0; ik <
nFiles; ik++)
694 delete CumulativeDistribution[ik];
699 for (
int ik = 0 ; ik <
nFiles; ik++)
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)
706 Cumulative += hpost[ik]->GetBinContent(j)/Integral;
707 CumulativeDistribution[ik]->SetBinContent(j, Cumulative);
710 CumulativeDistribution[ik]->SetBinContent(NumberOfBins+1, 1.);
713 std::vector<int> TestStatBin(
nFiles, 0);
714 std::vector<double> TestStatD(
nFiles, -999);
715 std::vector<std::unique_ptr<TLine>> LineD(
nFiles);
717 for (
int ik = 1 ; ik <
nFiles; ik++)
719 const int NumberOfBins = CumulativeDistribution[0]->GetXaxis()->GetNbins();
720 for (
int j = 1; j < NumberOfBins+1; ++j)
722 const double BinValue = CumulativeDistribution[0]->GetBinCenter(j);
723 const int BinNumber = CumulativeDistribution[ik]->FindBin(BinValue);
725 double TempDstat = std::fabs(CumulativeDistribution[0]->GetBinContent(j) - CumulativeDistribution[ik]->GetBinContent(BinNumber));
726 if(TempDstat > TestStatD[ik])
728 TestStatD[ik] = TempDstat;
734 for (
int ik = 0 ; ik <
nFiles; ik++)
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);
740 CumulativeDistribution[0]->Draw();
741 for (
int ik = 0 ; ik <
nFiles; ik++)
742 CumulativeDistribution[ik]->Draw(
"SAME");
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");
751 leg->SetLineColor(0);
752 leg->SetLineStyle(0);
753 leg->SetFillColor(0);
754 leg->SetFillStyle(0);
757 for (
int ik = 1; ik <
nFiles; ik++)
758 LineD[ik]->Draw(
"sam");
761 Posterior->Print(canvasname);
763 for (
int ik = 0; ik <
nFiles; ik++)
766 delete CumulativeDistribution[ik];
void RemoveFitter(TH1D *hist, const std::string &name)
KS: Remove fitted TF1 from hist to make comparison easier.
void SetMaCh3LoggerFormat()
Set messaging format of the logger.
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...
int main(int argc, char *argv[])
TH2D * TMatrixIntoTH2D(TMatrixDSym *Matrix, const std::string &title)
KS: Convert TMatrix to TH2D, mostly useful for making fancy plots.
std::vector< std::string > TitleNames
void CalcBipolarPlot(MCMCProcessor *Processor)
void CalcParameterEvolution(MCMCProcessor *Processor)
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.
std::vector< std::string > FileNames
void ReweightPrior(MCMCProcessor *Processor)
void CalcBayesFactor(MCMCProcessor *Processor)
void CalcSavageDickey(MCMCProcessor *Processor)
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...
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.
#define M3OpenConfig(filename)
Macro to simplify calling LoadYaml with file and line info.
Class responsible for processing MCMC chains, performing diagnostics, generating plots,...
void GetNthParameter(const int param, double &Prior, double &PriorError, TString &Title) const
Get properties of parameter by passing it number.
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.
void ResetHistograms()
Reset 2D posteriors, in case we would like to calculate in again with different BurnInCut.
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.
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 GetPolarPlot(const std::vector< std::string > &ParNames)
Make funny polar plot.
void ParameterEvolution(const std::vector< std::string > &Names, const std::vector< int > &NIntervals)
Make .gif of parameter evolution.
void GetCovariance(TMatrixDSym *&Cov, TMatrixDSym *&Corr)
Get the post-fit covariances and correlations.
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.
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.
void SetStepCut(const std::string &Cuts)
Set the step cutting by string.
Custom exception class for MaCh3 errors.