MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
ProcessMCMC.cpp
Go to the documentation of this file.
1//MaCh3 includes
3#include "Manager/Manager.h"
4
8
10inline void ProcessMCMC(const std::string& inputFile);
12inline void MultipleProcessMCMC();
13inline void CalcBayesFactor(MCMCProcessor* Processor);
14inline void CalcSavageDickey(MCMCProcessor* Processor);
15inline void CalcBipolarPlot(MCMCProcessor* Processor);
16inline void CalcParameterEvolution(MCMCProcessor* Processor);
17inline void GetTrianglePlot(MCMCProcessor* Processor);
18inline void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, const std::string& inputFile);
19inline void ReweightPrior(MCMCProcessor* Processor);
21inline TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, const std::string& title);
23inline void KolmogorovSmirnovTest(const std::vector<std::unique_ptr<MCMCProcessor>>& Processor,
24 const std::unique_ptr<TCanvas>& Posterior,
25 const TString& canvasname);
26
28std::vector <std::string> FileNames;
29std::vector <std::string> TitleNames;
30std::string config;
31
32int main(int argc, char *argv[])
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}
71
72void ProcessMCMC(const std::string& inputFile)
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}
166
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}
346
347// 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
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}
367
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}
385
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}
400
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}
413
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}
432
434void DiagnoseCovarianceMatrix(MCMCProcessor* Processor, const std::string& inputFile)
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}
608
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}
621
622//KS: Convert TMatrix to TH2D, mostly useful for making fancy plots
623TH2D* TMatrixIntoTH2D(TMatrixDSym* Matrix, const std::string& title)
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}
636
637//KS: Perform KS test to check if two posteriors for the same parameter came from the same distribution
638void KolmogorovSmirnovTest(const std::vector<std::unique_ptr<MCMCProcessor>>& Processor,
639 const std::unique_ptr<TCanvas>& Posterior,
640 const TString& canvasname)
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}
void RemoveFitter(TH1D *hist, const std::string &name)
KS: Remove fitted TF1 from hist to make comparison easier.
#define _UNDEF_
Definition: MCMCProcessor.h:4
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
void SetMaCh3LoggerFormat()
Set messaging format of the logger.
Definition: MaCh3Logger.h:30
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
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[])
Definition: ProcessMCMC.cpp:32
std::string config
Definition: ProcessMCMC.cpp:30
TH2D * TMatrixIntoTH2D(TMatrixDSym *Matrix, const std::string &title)
KS: Convert TMatrix to TH2D, mostly useful for making fancy plots.
std::vector< std::string > TitleNames
Definition: ProcessMCMC.cpp:29
void CalcBipolarPlot(MCMCProcessor *Processor)
void CalcParameterEvolution(MCMCProcessor *Processor)
void ProcessMCMC(const std::string &inputFile)
Main function processing MCMC and Producing plots.
Definition: ProcessMCMC.cpp:72
int nFiles
Definition: ProcessMCMC.cpp:27
void MultipleProcessMCMC()
Function producing comparison of posterior and more betwen a few MCMC chains.
std::vector< std::string > FileNames
Definition: ProcessMCMC.cpp:28
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.
Definition: YamlHelper.h:298
#define M3OpenConfig(filename)
Macro to simplify calling LoadYaml with file and line info.
Definition: YamlHelper.h:560
Class responsible for processing MCMC chains, performing diagnostics, generating plots,...
Definition: MCMCProcessor.h:65
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.