MaCh3  2.6.0
Reference Guide
Functions | Variables
PredictivePlotting.cpp File Reference
#include "PlottingUtils/PlottingUtils.h"
#include "PlottingUtils/PlottingManager.h"
#include <numeric>
Include dependency graph for PredictivePlotting.cpp:

Go to the source code of this file.

Functions

std::vector< std::string > FindSamples (const std::string &File)
 
std::vector< int > FindDimensions (const std::string &File, const std::vector< std::string > &Samples)
 
std::vector< std::vector< std::string > > FindModes (const std::string &File, const std::vector< std::string > &SampleNames)
 
void PretifyHistogram (TH1 *Hist, const std::string &SampleName)
 
double GetPValue (const TH2D *hist)
 
void PrintPosteriorPValue (const YAML::Node &Settings, const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames)
 
void OverlayViolin (const YAML::Node &Settings, const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames, const std::vector< int > &SampleDimension, const std::unique_ptr< TCanvas > &canv)
 
void OverlayPredicitve (const YAML::Node &Settings, const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames, const std::vector< int > &SampleDimension, const std::unique_ptr< TCanvas > &canv)
 
void OverlayPredicitveByMode (const YAML::Node &Settings, const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames, const std::vector< int > &SampleDimension, const std::vector< std::vector< std::string >> &Modes, const std::unique_ptr< TCanvas > &canv)
 
void GetMeanError (TH1D *hist, double &Mean, double &Error)
 KS: Get mean and error from gaussian fit to event distribution. More...
 
void PrintPosteriorEventRates (const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames)
 KS Print event rates in Latex like table. More...
 
void PrintPosteriorFractionalUncertainties (const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames)
 KS: Print Fractional Uncertainties into Latex table format. More...
 
double GetLLH (TH1 *hist)
 
void PrintPredictiveLLH (const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames)
 KS Print Predictive LLH into Latex table format. More...
 
void PredictivePlotting (const std::string &ConfigName, const std::vector< std::string > &FileNames)
 
int main (int argc, char **argv)
 

Variables

M3::Plotting::PlottingManager * PlotMan = nullptr
 

Detailed Description

Author
Kamil Skwarczynski

Definition in file PredictivePlotting.cpp.

Function Documentation

◆ FindDimensions()

std::vector<int> FindDimensions ( const std::string &  File,
const std::vector< std::string > &  Samples 
)

Definition at line 49 of file PredictivePlotting.cpp.

50 {
51  TFile *file = M3::Open(File, "READ", __FILE__, __LINE__);
52  TDirectoryFile *PredicitveDir = file->Get<TDirectoryFile>("Predictive");
53 
54  std::vector<int> SampleDimension;
55  for (const auto& sample : Samples)
56  {
57  // Get directory for this sample
58  TDirectoryFile* SampleDir = PredicitveDir->Get<TDirectoryFile>(sample.c_str());
59 
60  int Dimension = 0;
61 
62  while (true)
63  {
64  // Construct name Tutorial_mc_dimX
65  std::string histName = fmt::format("{}_mc_dim{}", sample, Dimension);
66 
67  TH2D* hist = SampleDir->Get<TH2D>(histName.c_str());
68  if (!hist) break; // stop when next dimension does not exist
69 
70  Dimension++;
71  }
72 
73  MACH3LOG_DEBUG("Sample '{}' has dimension {}", sample, Dimension);
74  SampleDimension.push_back(Dimension);
75  }
76 
77  file->Close();
78  delete file;
79 
80  return SampleDimension;
81 }
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:34
TFile * Open(const std::string &Name, const std::string &Type, const std::string &File, const int Line)
Opens a ROOT file with the given name and mode.

◆ FindModes()

std::vector<std::vector<std::string> > FindModes ( const std::string &  File,
const std::vector< std::string > &  SampleNames 
)

Definition at line 84 of file PredictivePlotting.cpp.

86 {
87  TFile *file = M3::Open(File, "READ", __FILE__, __LINE__);
88  TDirectoryFile *PredictiveDir = file->Get<TDirectoryFile>("Predictive");
89 
90  std::vector<std::vector<std::string>> ModeNames(SampleNames.size());
91 
92  for(size_t iSample = 0; iSample < SampleNames.size(); iSample++)
93  {
94  TDirectoryFile* SampleDir = PredictiveDir->Get<TDirectoryFile>(SampleNames[iSample].c_str());
95  if(!SampleDir) continue;
96 
97  TDirectoryFile* ByModeDir = SampleDir->Get<TDirectoryFile>("ByMode");
98  if(!ByModeDir) continue;
99 
100  // Loop over all keys in ByModeDir
101  TIter next(ByModeDir->GetListOfKeys());
102  TKey* key;
103 
104  while ((key = static_cast<TKey*>(next())))
105  {
106  TObject* obj = key->ReadObj();
107  if (!obj->InheritsFrom("TH1")) continue;
108 
109  std::string histName = obj->GetName();
110 
111  // Example: sample_mode_dim0 → extract "mode"
112  std::string prefix = SampleNames[iSample] + "_";
113  std::string suffix = "_dim0";
114 
115  if (histName.find(prefix) == 0 &&
116  histName.rfind(suffix) == histName.size() - suffix.size())
117  {
118  std::string Mode = histName.substr(prefix.size(),
119  histName.size() - prefix.size() - suffix.size());
120  MACH3LOG_DEBUG("Found mode '{}' for sample {}", Mode, SampleNames[iSample]);
121  ModeNames[iSample].push_back(Mode);
122  }
123  }
124  }
125 
126  file->Close();
127  delete file;
128 
129  return ModeNames;
130 }

◆ FindSamples()

std::vector<std::string> FindSamples ( const std::string &  File)

Definition at line 17 of file PredictivePlotting.cpp.

18 {
19  TFile *file = M3::Open(File, "READ", __FILE__, __LINE__);
20  TDirectoryFile *PredicitveDir = file->Get<TDirectoryFile>("Predictive");
21 
22  std::vector<std::string> SampleNames;
23  //Get all entries in input file
24  TIter next(PredicitveDir->GetListOfKeys());
25  TKey *key = nullptr;
26 
27  // Loop through all entries
28  while ((key = static_cast<TKey*>(next()))) {
29  // get directory names, ignore flux
30  auto classname = std::string(key->GetClassName());
31  auto dirname = std::string(key->GetName());
32 
33  if (classname != "TDirectoryFile") continue;
34  dirname = std::string(key->GetName());
35 
36  if(dirname == "Total") continue;
37  if(dirname == "BetaParameters") continue;
38  if(dirname == "Correlations") continue;
39 
40  SampleNames.push_back(dirname);
41  MACH3LOG_DEBUG("Entering Sample {}", dirname);
42  }
43 
44  file->Close();
45  delete file;
46  return SampleNames;
47 }

◆ GetLLH()

double GetLLH ( TH1 *  hist)

Definition at line 672 of file PredictivePlotting.cpp.

673 {
674  std::string TempTile = hist->GetTitle();
675  std::string temp = "=";
676 
677  std::string::size_type SizeType = TempTile.find(temp);
678  TempTile.erase(0, SizeType+1);
679  double llh = atof(TempTile.c_str());
680  return llh;
681 }

◆ GetMeanError()

void GetMeanError ( TH1D *  hist,
double &  Mean,
double &  Error 
)

KS: Get mean and error from gaussian fit to event distribution.

Definition at line 564 of file PredictivePlotting.cpp.

564  {
565  TF1 *Gauss = hist->GetFunction("Fit"); //This name is hardcoded be careful
566  //KS: Get mean and error from Gauss
567  Mean = Gauss->GetParameter(1);
568  Error = Gauss->GetParameter(2);
569 
570  //KS: Get mean and error from HPD
571  //Mean = hist->GetMean();
572  //Error = hpost->GetRMS();
573 }
double ** Mean
Definition: RHat.cpp:63

◆ GetPValue()

double GetPValue ( const TH2D *  hist)

Definition at line 142 of file PredictivePlotting.cpp.

143 {
144  double pvalue = 0;
145  std::string TempTile = hist->GetTitle();
146  std::string temp = "=";
147 
148  std::string::size_type SizeType = TempTile.find(temp);
149  TempTile.erase(0, SizeType+1);
150  pvalue = atof(TempTile.c_str());
151  return pvalue;
152 }

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 789 of file PredictivePlotting.cpp.

790 {
792  if (argc < 3)
793  {
794  MACH3LOG_ERROR("Need at least two arguments, {} <Config.Yaml> <Prior/Post_PredOutput.root>", argv[0]);
795  throw MaCh3Exception(__FILE__, __LINE__);
796  }
797  std::string ConfigName = std::string(argv[1]);
798  // Collect all remaining arguments as file names
799  std::vector<std::string> FileNames;
800  for (int i = 2; i < argc; ++i) {
801  FileNames.emplace_back(argv[i]);
802  }
803 
804  PlotMan = new M3::Plotting::PlottingManager();
805  PlotMan->initialise();
806 
807  PredictivePlotting(ConfigName, FileNames);
808 
809  if(PlotMan) delete PlotMan;
810  return 0;
811 }
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:37
void SetMaCh3LoggerFormat()
Set messaging format of the logger.
Definition: MaCh3Logger.h:60
void PredictivePlotting(const std::string &ConfigName, const std::vector< std::string > &FileNames)
M3::Plotting::PlottingManager * PlotMan
std::vector< std::string > FileNames
Definition: ProcessMCMC.cpp:28
Custom exception class used throughout MaCh3.

◆ OverlayPredicitve()

void OverlayPredicitve ( const YAML::Node &  Settings,
const std::vector< TFile * > &  InputFiles,
const std::vector< std::string > &  SampleNames,
const std::vector< int > &  SampleDimension,
const std::unique_ptr< TCanvas > &  canv 
)

Definition at line 278 of file PredictivePlotting.cpp.

283 {
284  MACH3LOG_INFO("Starting {}", __func__);
285  canv->Clear();
286 
287  TPad* pad1 = new TPad("pad1","pad1",0,0.25,1,1);
288  pad1->AppendPad();
289  TPad* pad2 = new TPad("pad2","pad2",0,0,1,0.25);
290  pad2->AppendPad();
291 
292  pad1->SetGrid();
293  pad2->SetGrid();
294 
295  pad1->SetLeftMargin(canv->GetLeftMargin());
296  pad1->SetRightMargin(canv->GetRightMargin());
297  pad1->SetTopMargin(canv->GetTopMargin());
298  pad1->SetBottomMargin(0);
299 
300  pad2->SetLeftMargin(canv->GetLeftMargin());
301  pad2->SetRightMargin(canv->GetRightMargin());
302  pad2->SetTopMargin(0);
303  pad2->SetBottomMargin(0.28);
304 
305  auto PosteriorColor = Get<std::vector<Color_t >>(Settings["PosteriorColor"], __FILE__, __LINE__);
306  auto Titles = Get<std::vector<std::string>>(Settings["FileTitle"], __FILE__, __LINE__);
307 
308  if(Titles.size() < InputFiles.size() || PosteriorColor.size() < InputFiles.size()){
309  MACH3LOG_ERROR("Passed {} files, while only {} titles and {} colors", InputFiles.size(), Titles.size(), PosteriorColor.size());
310  throw MaCh3Exception(__FILE__, __LINE__);
311  }
312  for(size_t iSample = 0; iSample < SampleNames.size(); iSample++)
313  {
314  const int nFiles = static_cast<int>(InputFiles.size());
315  auto SampleName = SampleNames[iSample];
316  const int nDims = (SampleDimension[iSample] == 2) ? 2 : 1;
317  for(int iDim = 0; iDim < nDims; iDim++) {
318  std::string DataLocation = "";
319  if(nDims == 2) {
320  DataLocation = "Predictive/" + SampleName + "/Data_" + SampleName + "_Dim" + std::to_string(iDim);
321  } else {
322  DataLocation = "SampleFolder/data_" + SampleName;
323  }
324  TH1D* hist = InputFiles[0]->Get<TH1D>((DataLocation).c_str());
325 
326  auto BinWidthScale = PlotMan->style().getBinWidthScale(hist->GetXaxis()->GetTitle());
327  std::unique_ptr<TH1D> DataHist = M3::Clone(hist);
328  auto DataPoissonErrors = PoissonGraphScaled(DataHist.get(), BinWidthScale);
329  M3::ScaleHistogram(DataHist.get(), BinWidthScale);
330 
331  DataHist->SetLineColor(kBlack);
332  DataPoissonErrors->SetLineColor(kBlack);
333  //KS: +1 for data, we want to get integral before scaling of the histogram
334  std::vector<double> Integral(nFiles+1);
335  Integral[nFiles] = DataHist->Integral();
336  std::vector<std::unique_ptr<TH1D>> PredHist(nFiles);
337 
338  for(int iFile = 0; iFile < nFiles; iFile++)
339  {
340  InputFiles[iFile]->cd();
341  std::string HistLocation = "";
342  if(nDims == 2) {
343  HistLocation = "Predictive/" + SampleName + "/" + SampleName + "_mc_PostPred_dim" + std::to_string(iDim);
344  } else {
345  HistLocation = "Predictive/" + SampleName + "/" + SampleName + "_mc_PostPred";
346  }
347  PredHist[iFile] = M3::Clone(InputFiles[iFile]->Get<TH1D>((HistLocation).c_str()));
348  Integral[iFile] = PredHist[iFile]->Integral();
349  PredHist[iFile]->SetLineColor(PosteriorColor[iFile]);
350  PredHist[iFile]->SetMarkerColor(PosteriorColor[iFile]);
351  PredHist[iFile]->SetFillColorAlpha(PosteriorColor[iFile], 0.35);
352  PredHist[iFile]->SetFillStyle(1001);
353  PretifyHistogram(PredHist[iFile].get(), SampleName);
354  }
355  pad1->cd();
356 
357  PredHist[0]->Draw("p e2");
358  for(int iFile = 1; iFile < nFiles; iFile++) {
359  PredHist[iFile]->Draw("p e2 same");
360  }
361  DataPoissonErrors->Draw("p same");
362 
363  auto legend = std::make_unique<TLegend>(0.50,0.52,0.90,0.88);
364  legend->AddEntry(DataPoissonErrors.get(), Form("Data, #int=%.0f", Integral[nFiles]),"le");
365  for(int ig = 0; ig < nFiles; ig++ ) {
366  legend->AddEntry(PredHist[ig].get(), Form("%s, #int=%.2f", Titles[ig].c_str(), Integral[ig]), "lpf");
367  }
368  legend->SetLineStyle(0);
369  legend->SetTextSize(0.03);
370  legend->Draw();
371 
373  pad2->cd();
374 
375  auto line = std::make_unique<TLine>(PredHist[0]->GetXaxis()->GetBinLowEdge(PredHist[0]->GetXaxis()->GetFirst()), 1.0, PredHist[0]->GetXaxis()->GetBinUpEdge(PredHist[0]->GetXaxis()->GetLast()), 1.0);
376 
377  line->SetLineWidth(2);
378  line->SetLineColor(kBlack);
379  line->Draw("");
380 
381  std::unique_ptr<TH1D> RatioPlotData = M3::Clone(DataHist.get());
382  std::vector<std::unique_ptr<TH1D>> RatioPlot(nFiles);
383 
384  for(int ig = 0; ig < nFiles; ig++ )
385  {
386  RatioPlot[ig] = M3::Clone(DataHist.get());
387  RatioPlot[ig]->SetLineColor(PosteriorColor[ig]);
388  RatioPlot[ig]->SetMarkerColor(PosteriorColor[ig]);
389  RatioPlot[ig]->SetFillColorAlpha(PosteriorColor[ig], 0.35);
390  RatioPlot[ig]->SetFillStyle(1001);
391  RatioPlot[ig]->GetYaxis()->SetTitle("Data/MC");
392  auto PrettyX = PlotMan->style().prettifyKinematicName(PredHist[0]->GetXaxis()->GetTitle());
393  RatioPlot[ig]->GetXaxis()->SetTitle(PrettyX.c_str());
394  RatioPlot[ig]->SetBit(TH1D::kNoTitle);
395  RatioPlot[ig]->GetXaxis()->SetTitleSize(0.12);
396  RatioPlot[ig]->GetYaxis()->SetTitleOffset(0.4);
397  RatioPlot[ig]->GetYaxis()->SetTitleSize(0.10);
398 
399  RatioPlot[ig]->GetXaxis()->SetLabelSize(0.10);
400  RatioPlot[ig]->GetYaxis()->SetLabelSize(0.10);
401 
402  RatioPlot[ig]->Divide(PredHist[ig].get());
403  PassErrorToRatioPlot(RatioPlot[ig].get(), PredHist[ig].get(), DataHist.get());
404  }
405 
406  RatioPlotData->Divide(DataHist.get());
407  PassErrorToRatioPlot(RatioPlotData.get(), DataHist.get(), DataHist.get());
408 
409  M3::Plotting::SetSymmetricRatioRange(RatioPlot);
410 
411  RatioPlot[0]->Draw("p e2");
412  for(int ig = 1; ig < nFiles; ig++ ) {
413  RatioPlot[ig]->Draw("p e2 same");
414  }
415  RatioPlotData->Draw("he same");
416 
417  canv->Print("Overlay_Predictive.pdf", "pdf");
418  }
419  }
420 
421  delete pad1;
422  delete pad2;
423 }
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:35
void PretifyHistogram(TH1 *Hist, const std::string &SampleName)
int nFiles
Definition: ProcessMCMC.cpp:27
std::unique_ptr< TGraphAsymmErrors > PoissonGraphScaled(const TH1D *hist, double scale, double cl)
Create a TGraphAsymmErrors from a histogram using exact Poisson confidence intervals instead of symme...
void PassErrorToRatioPlot(TH1D *RatioHist, TH1D *Hist1, TH1D *DataHist)
Propagate numerator uncertainties to a ratio histogram.
std::unique_ptr< ObjectType > Clone(const ObjectType *obj, const std::string &name="")
KS: Creates a copy of a ROOT-like object and wraps it in a smart pointer.
void ScaleHistogram(TH1 *Sample_Hist, const double scale)
Scale histogram to get divided by bin width.

◆ OverlayPredicitveByMode()

void OverlayPredicitveByMode ( const YAML::Node &  Settings,
const std::vector< TFile * > &  InputFiles,
const std::vector< std::string > &  SampleNames,
const std::vector< int > &  SampleDimension,
const std::vector< std::vector< std::string >> &  Modes,
const std::unique_ptr< TCanvas > &  canv 
)

Definition at line 425 of file PredictivePlotting.cpp.

431 {
432  MACH3LOG_INFO("Starting {}", __func__);
433  canv->cd();
434  constexpr auto DefaultColor = kBlack;
435  auto Titles = Get<std::vector<std::string>>(Settings["FileTitle"], __FILE__, __LINE__);
436  auto RelevantModesName = Get<std::vector<std::string>>(Settings["RelevantModesName"], __FILE__, __LINE__);
437  auto RelevantColors = Get<std::vector<Color_t>>(Settings["RelevantModesColors"], __FILE__, __LINE__);
438  int nRelevantModes = static_cast<int>(RelevantModesName.size());
439  const int nFiles = static_cast<int>(InputFiles.size());
440  if(Titles.size() < InputFiles.size()){
441  MACH3LOG_ERROR("Passed {} files, while only {} titles", InputFiles.size(), Titles.size());
442  throw MaCh3Exception(__FILE__, __LINE__);
443  }
444  if(RelevantModesName.size() != RelevantColors.size()) {
445  MACH3LOG_ERROR("Colors ({}) doesn't match relevant modes {}", RelevantColors.size(), RelevantModesName.size());
446  throw MaCh3Exception(__FILE__, __LINE__);
447  }
448  for(int iFile = 0; iFile < nFiles; iFile++ )
449  {
450  for(size_t iSample = 0; iSample < SampleNames.size(); iSample++)
451  {
452  auto SampleName = SampleNames[iSample];
453  for(int iDim = 0; iDim < SampleDimension[iSample]; iDim++)
454  {
455  const int nDims = (SampleDimension[iSample] == 2) ? 2 : 1;
456  std::string HistLocation = "";
457  if(nDims == 2) {
458  HistLocation = "Predictive/" + SampleName + "/" + SampleName + "_mc_PostPred_dim" + std::to_string(iDim);
459  } else {
460  HistLocation = "Predictive/" + SampleName + "/" + SampleName + "_mc_PostPred";
461  }
462  std::unique_ptr<TH1D> Sample_MC_Full = M3::Clone(InputFiles[iFile]->Get<TH1D>((HistLocation).c_str()));
463  Sample_MC_Full->SetLineColor(kOrange);
464  Sample_MC_Full->SetLineWidth(2);
465  Sample_MC_Full->SetMarkerColor(kOrange);
466  PretifyHistogram(Sample_MC_Full.get(), SampleName);
467 
468  std::string DataLocation = "";
469  std::unique_ptr<TH1D> Sample_Data;
470  if(nDims == 2) {
471  DataLocation = "Predictive/" + SampleName + "/Data_" + SampleName + "_Dim" + std::to_string(iDim);
472  } else if(nDims == 1) {
473  DataLocation = "SampleFolder/data_" + SampleName;
474  }
475  if(DataLocation != "") {
476  Sample_Data = M3::Clone(InputFiles[iFile]->Get<TH1D>((DataLocation).c_str()));
477  Sample_Data->SetLineColor(kBlack);
478  Sample_Data->SetLineWidth(2);
479  Sample_Data->SetMarkerColor(kBlack);
480  PretifyHistogram(Sample_Data.get(), SampleName);
481  }
482  int nModes = static_cast<int>(Modes[iSample].size());
483  // Simple map to keep track which mode is relevant and which will be added to "Other"
484  std::vector<bool> isRelevantMode(nModes, false);
485  std::vector<Color_t > ColorMap(nModes, DefaultColor);
486  for(int iMode = 0; iMode < nModes; iMode++) {
487  for(int iRelevant = 0; iRelevant < nRelevantModes; iRelevant++) {
488  if(Modes[iSample][iMode] == RelevantModesName[iRelevant]) {
489  isRelevantMode[iMode] = true;
490  ColorMap[iMode] = RelevantColors[iRelevant];
491  }
492  }
493  }
494  auto Sample_Stack = std::make_unique<THStack>(SampleName.c_str(), SampleName.c_str());
495  // This will hold values for "Other" modes
496  std::unique_ptr<TH1D> Sample_MC_Other;
497  // KS: Store histogram for each mode
498  std::vector<std::unique_ptr<TH1D>> Sample_MC(nModes);
499  std::vector<double> Integrals(nModes, 0.);
500  for(int iMode = 0; iMode < nModes; iMode++)
501  {
502  std::string FileLocaction = "Predictive/" + SampleName + "/ByMode/" + SampleName
503  + "_" + Modes[iSample][iMode] + "_dim" + std::to_string(iDim);
504  auto SpectraByMode = InputFiles[iFile]->Get<TH2D>((FileLocaction).c_str());
505  if(SpectraByMode == nullptr){
506  MACH3LOG_ERROR("Something went wrong and didn't find histogram: {}", FileLocaction);
507  throw MaCh3Exception(__FILE__, __LINE__);
508  }
509  Sample_MC[iMode] = MakeSummaryFromSpectra(SpectraByMode, SpectraByMode->GetTitle());
510  Integrals[iMode] = Sample_MC[iMode]->Integral();
511  PretifyHistogram(Sample_MC[iMode].get(), SampleName);
512 
513  if(Sample_MC_Other == nullptr) {
514  Sample_MC_Other = M3::Clone(Sample_MC[iMode].get());
515  Sample_MC_Other->Reset();
516  Sample_MC_Other->SetFillColor(DefaultColor);
517  Sample_MC_Other->SetLineColor(DefaultColor);
518  }
519  if(!isRelevantMode[iMode]) {
520  Sample_MC_Other->Add(Sample_MC[iMode].get());
521  }
522  Sample_MC[iMode]->SetFillColor(ColorMap[iMode]);
523  Sample_MC[iMode]->SetLineColor(ColorMap[iMode]);
524  } // end loop over modes
525  Sample_Stack->Add(Sample_MC_Other.get());
526  // KS: We do this other way around as we want to have most relevant modes first
527  for(int iMode = nModes-1; iMode >= 0; iMode--) {
528  if(isRelevantMode[iMode]) Sample_Stack->Add( Sample_MC[iMode].get() );
529  }
530  Sample_Stack->Draw("hist");
531  Sample_MC_Full->Draw("SAME he");
532  if(Sample_Data) Sample_Data->Draw("SAME pe");
533  canv->cd();
534  Sample_Stack->GetXaxis();
535  Sample_Stack->SetTitle(Sample_MC_Other->GetTitle());
536  Sample_Stack->GetXaxis()->SetTitle(Sample_MC_Other->GetXaxis()->GetTitle());
537  Sample_Stack->GetYaxis()->SetTitle(Sample_MC_Other->GetYaxis()->GetTitle());
538 
539  double FullIntegral = std::accumulate(Integrals.begin(), Integrals.end(), 0.0);
540  double OtherIntegral = 0.;
541  TLegend legend(0.50,0.52,0.85,0.88);
542  if(Sample_Data) legend.AddEntry(Sample_Data.get(), "Data","ple");
543  legend.AddEntry(Sample_MC_Full.get(), Titles[iFile].c_str(),"fple");
544  for(int iMode = 0; iMode < nModes; iMode++) {
545  if(isRelevantMode[iMode]) {
546  std::string Label = Form("%s: %.1f%%", Modes[iSample][iMode].c_str(), Integrals[iMode]/FullIntegral*100);
547  legend.AddEntry(Sample_MC[iMode].get(), Label.c_str() ,"lf");
548  } else{
549  OtherIntegral += Integrals[iMode]/FullIntegral;
550  }
551  }
552  legend.AddEntry(Sample_MC_Other.get(), Form("Other: %.1f%%", OtherIntegral*100), "lf");
553  legend.SetTextSize(0.03);
554  legend.Draw();
555 
556  canv->Print("Overlay_Predictive.pdf", "pdf");
557  } // end loop over dimensions
558  } // end loop over samples
559  } // end loop over files
560 }
std::unique_ptr< TH1D > MakeSummaryFromSpectra(const TH2D *Spectra, const std::string &name)
Build a 1D posterior-predictive summary from a violin spectrum.

◆ OverlayViolin()

void OverlayViolin ( const YAML::Node &  Settings,
const std::vector< TFile * > &  InputFiles,
const std::vector< std::string > &  SampleNames,
const std::vector< int > &  SampleDimension,
const std::unique_ptr< TCanvas > &  canv 
)

Definition at line 220 of file PredictivePlotting.cpp.

225 {
226  MACH3LOG_INFO("Starting {}", __func__);
227  canv->Clear();
228 
229  canv->SetTopMargin(0.10);
230  canv->SetBottomMargin(0.12);
231  canv->SetRightMargin(0.075);
232  canv->SetLeftMargin(0.14);
233 
234  auto PosteriorColor = Get<std::vector<Color_t >>(Settings["PosteriorColor"], __FILE__, __LINE__);
235  auto Titles = Get<std::vector<std::string>>(Settings["FileTitle"], __FILE__, __LINE__);
236  const int nFiles = static_cast<int>(InputFiles.size());
237 
238  //KS: No idea why but ROOT changed treatment of violin in R6. If you have non uniform binning this will results in very hard to see violin plots.
239  TCandle::SetScaledViolin(false);
240  for(size_t iSample = 0; iSample < SampleNames.size(); iSample++)
241  {
242  for(int iDim = 0; iDim < SampleDimension[iSample]; iDim++)
243  {
244  std::vector<std::unique_ptr<TH2D>> ViolinHist(nFiles);
245  for(int iFile = 0; iFile < nFiles; iFile++)
246  {
247  InputFiles[iFile]->cd();
248  ViolinHist[iFile] = M3::Clone(InputFiles[iFile]->Get<TH2D>(("Predictive/" + SampleNames[iSample]
249  + "/" + SampleNames[iSample] + "_mc_dim" + iDim).Data()));
250  ViolinHist[iFile]->SetTitle(PlotMan->style().prettifySampleName(SampleNames[iSample]).c_str());
251  ViolinHist[iFile]->SetLineColor(PosteriorColor[iFile]);
252  ViolinHist[iFile]->SetMarkerColor(PosteriorColor[iFile]);
253  ViolinHist[iFile]->SetFillColorAlpha(PosteriorColor[iFile], 0.35);
254  ViolinHist[iFile]->SetFillStyle(1001);
255  ViolinHist[iFile]->GetXaxis()->SetTitle(PlotMan->style().prettifyKinematicName(
256  ViolinHist[iFile]->GetXaxis()->GetTitle()).c_str());
257  ViolinHist[iFile]->GetYaxis()->SetTitle("Events");
258  }
259 
260  ViolinHist[0]->Draw("violinX(03100300)");
261  for(int iFile = 1; iFile < nFiles; iFile++) {
262  ViolinHist[iFile]->Draw("violinX(03100300) same");
263  }
264 
265  TLegend legend(0.50, 0.52, 0.90, 0.88);
266  for(int ig = 0; ig < nFiles; ig++) {
267  legend.AddEntry(ViolinHist[ig].get(), Form("%s", Titles[ig].c_str()), "lpf");
268  }
269  legend.SetLineStyle(0);
270  legend.SetTextSize(0.03);
271  legend.Draw();
272 
273  canv->Print("Overlay_Predictive.pdf", "pdf");
274  }
275  }
276 }

◆ PredictivePlotting()

void PredictivePlotting ( const std::string &  ConfigName,
const std::vector< std::string > &  FileNames 
)

Definition at line 731 of file PredictivePlotting.cpp.

733 {
734  auto canvas = std::make_unique<TCanvas>("canv", "canv", 1080, 1080);
735  // set the paper & margin sizes
736  canvas->SetTopMargin(0.11);
737  canvas->SetBottomMargin(0.16);
738  canvas->SetRightMargin(0.075);
739  canvas->SetLeftMargin(0.12);
740  canvas->SetGrid();
741 
742  gStyle->SetOptStat(0); //Set 0 to disable statistic box
743  gStyle->SetPalette(51);
744  gStyle->SetLegendBorderSize(0); //This option disables legends borders
745  gStyle->SetFillStyle(0);
746 
747  //To avoid TCanvas::Print> messages
748  gErrorIgnoreLevel = kWarning;
749 
750  auto Samples = FindSamples(FileNames[0]);
751  auto Dimensions = FindDimensions(FileNames[0], Samples);
752  auto Modes = FindModes(FileNames[0], Samples);
753 
754  std::vector<TFile*> InputFiles(FileNames.size());
755  for(size_t i = 0; i < FileNames.size(); i++) {
756  InputFiles[i] = M3::Open(FileNames[i], "READ", __FILE__, __LINE__);
757  }
758 
759  // Load the YAML file
760  YAML::Node Config = M3OpenConfig(ConfigName);
761  // Access the "MatrixPlotter" section
762  YAML::Node settings = Config["PredictivePlotting"];
763  canvas->Print("Overlay_Predictive.pdf[", "pdf");
764 
765  // Make overlay of 1D hists
766  OverlayPredicitve(settings, InputFiles, Samples, Dimensions, canvas);
767  // Make overlay of violin plots
768  OverlayViolin(settings, InputFiles, Samples, Dimensions, canvas);
769  // Make By Mode post pred
770  if(Modes[0].size() != 0) OverlayPredicitveByMode(settings, InputFiles, Samples, Dimensions, Modes, canvas);
771  // Get PValue per sample
772  PrintPosteriorPValue(settings, InputFiles, Samples);
773  // KS: Print Fractional Uncertainties into Latex table format
774  PrintPosteriorEventRates(InputFiles, Samples);
775  // KS: Print Fractional Uncertainties into Latex table format
776  PrintPosteriorFractionalUncertainties(InputFiles, Samples);
777  // KS: Print Predictive LLH into Latex table format
778  PrintPredictiveLLH(InputFiles, Samples);
779  canvas->Print("Overlay_Predictive.pdf]", "pdf");
780 
781  for(size_t i = 0; i < FileNames.size(); i++)
782  {
783  InputFiles[i]->Close();
784  delete InputFiles[i];
785  }
786 }
void OverlayPredicitve(const YAML::Node &Settings, const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames, const std::vector< int > &SampleDimension, const std::unique_ptr< TCanvas > &canv)
void PrintPosteriorPValue(const YAML::Node &Settings, const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames)
std::vector< std::string > FindSamples(const std::string &File)
void OverlayPredicitveByMode(const YAML::Node &Settings, const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames, const std::vector< int > &SampleDimension, const std::vector< std::vector< std::string >> &Modes, const std::unique_ptr< TCanvas > &canv)
void PrintPredictiveLLH(const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames)
KS Print Predictive LLH into Latex table format.
void OverlayViolin(const YAML::Node &Settings, const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames, const std::vector< int > &SampleDimension, const std::unique_ptr< TCanvas > &canv)
std::vector< std::vector< std::string > > FindModes(const std::string &File, const std::vector< std::string > &SampleNames)
std::vector< int > FindDimensions(const std::string &File, const std::vector< std::string > &Samples)
void PrintPosteriorFractionalUncertainties(const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames)
KS: Print Fractional Uncertainties into Latex table format.
void PrintPosteriorEventRates(const std::vector< TFile * > &InputFiles, const std::vector< std::string > &SampleNames)
KS Print event rates in Latex like table.
#define M3OpenConfig(filename)
Macro to simplify calling LoadYaml with file and line info.
Definition: YamlHelper.h:589

◆ PretifyHistogram()

void PretifyHistogram ( TH1 *  Hist,
const std::string &  SampleName 
)

Definition at line 133 of file PredictivePlotting.cpp.

133  {
134  Hist->SetTitle(PlotMan->style().prettifySampleName(SampleName).c_str());
135  auto BinWidthScale = PlotMan->style().getBinWidthScale(Hist->GetXaxis()->GetTitle());
136  auto PrettyX = PlotMan->style().prettifyKinematicName(Hist->GetXaxis()->GetTitle());
137  Hist->GetXaxis()->SetTitle(PrettyX.c_str());
138  Hist->GetYaxis()->SetTitle(fmt::format("Events/{:.0f}", BinWidthScale).c_str());
139  M3::ScaleHistogram(Hist, BinWidthScale);
140 }

◆ PrintPosteriorEventRates()

void PrintPosteriorEventRates ( const std::vector< TFile * > &  InputFiles,
const std::vector< std::string > &  SampleNames 
)

KS Print event rates in Latex like table.

Definition at line 576 of file PredictivePlotting.cpp.

577  {
578  MACH3LOG_INFO("Starting {}", __func__);
579  MACH3LOG_INFO("");
580 
581  double mean, error;
582  //KS: We now prepare to make tables for TN etc.
583  std::cout<<"\\begin{table}[htb]"<<std::endl;
584  std::cout<<"\\centering"<<std::endl;
585  std::cout<<"\\begin{tabular}{ | l |";
586  for(unsigned int f = 0; f < InputFiles.size(); f++)
587  {
588  std::cout<<" c |";
589  }
590  std::cout<<"} \\hline"<<std::endl;
591  std::cout<<"Sample ";
592  for(unsigned int f = 0; f < InputFiles.size(); f++)
593  {
594  std::cout<<"& Event Rates ";
595  }
596  std::cout<<"\\\\ \\hline"<<std::endl;
597  for(unsigned int i = 0; i < SampleNames.size(); i++)
598  {
599  std::cout<<SampleNames[i];
600  std::string TempString = "Predictive/" + SampleNames[i]+"/"+SampleNames[i]+"_sum";
601  for(unsigned int f = 0; f < InputFiles.size(); f++)
602  {
603  TH1D *hist = static_cast<TH1D*>(InputFiles[f]->Get(TempString.c_str()));
604  GetMeanError(hist, mean, error);
605  std::cout<<" & "<<mean<<" $\\pm$ "<<error;
606  }
607  std::cout<<" \\\\"<<std::endl;
608  }
609  std::cout<<"Total";
610  for(unsigned int f = 0; f < InputFiles.size(); f++)
611  {
612  TH1D *histTot = static_cast<TH1D*>(InputFiles[f]->Get("Predictive/Total/Total_sum"));
613  GetMeanError(histTot, mean, error);
614  std::cout<<" & "<<mean<<" $\\pm$ "<<error;
615  }
616  std::cout<<" \\\\"<<std::endl;
617  std::cout<<"\\hline"<<std::endl;
618  std::cout<<"\\end{tabular}"<<std::endl;
619  std::cout<<"\\end{table}"<<std::endl;
620  MACH3LOG_INFO("");
621 }
void GetMeanError(TH1D *hist, double &Mean, double &Error)
KS: Get mean and error from gaussian fit to event distribution.

◆ PrintPosteriorFractionalUncertainties()

void PrintPosteriorFractionalUncertainties ( const std::vector< TFile * > &  InputFiles,
const std::vector< std::string > &  SampleNames 
)

KS: Print Fractional Uncertainties into Latex table format.

Definition at line 624 of file PredictivePlotting.cpp.

625  {
626  MACH3LOG_INFO("Starting {}", __func__);
627  MACH3LOG_INFO("");
628  double mean, error;
629 
630  //KS: Fractional uncertainties on the prior and posterior predictive event rates.
631  std::cout<<"\\begin{table}[htb]"<<std::endl;
632  std::cout<<"\\centering"<<std::endl;
633  std::cout<<"\\begin{tabular}{ | l |";
634  for(unsigned int f = 0; f < InputFiles.size(); f++)
635  {
636  std::cout<<" c |";
637  }
638  std::cout<<"} \\hline"<<std::endl;
639 
640  std::cout<<"Sample ";
641  for(unsigned int f = 0; f < InputFiles.size(); f++)
642  {
643  std::cout<<"& $\\delta N / N (\\%)$";
644  }
645  std::cout<<"\\\\ \\hline"<<std::endl;
646 
647  for(unsigned int i = 0; i < SampleNames.size(); i++)
648  {
649  std::cout<<SampleNames[i];
650  std::string TempString = "Predictive/" + SampleNames[i]+"/"+SampleNames[i]+"_sum";
651  for(unsigned int f = 0; f < InputFiles.size(); f++)
652  {
653  TH1D *hist = static_cast<TH1D*>(InputFiles[f]->Get(TempString.c_str()));
654  GetMeanError(hist, mean, error);
655  std::cout<<" & "<<error/mean*100;
656  }
657  std::cout<<" \\\\"<<std::endl;
658  }
659  std::cout<<"Total";
660  for(unsigned int f = 0; f < InputFiles.size(); f++)
661  {
662  TH1D *histTotal = static_cast<TH1D*>(InputFiles[f]->Get("Predictive/Total/Total_sum"));
663  GetMeanError(histTotal, mean, error);
664  std::cout<<" & "<<error/mean*100;
665  }
666  std::cout<<"\\\\ \\hline"<<std::endl;
667  std::cout<<"\\end{tabular}"<<std::endl;
668  std::cout<<"\\end{table}"<<std::endl;
669  MACH3LOG_INFO("");
670 }

◆ PrintPosteriorPValue()

void PrintPosteriorPValue ( const YAML::Node &  Settings,
const std::vector< TFile * > &  InputFiles,
const std::vector< std::string > &  SampleNames 
)

Definition at line 154 of file PredictivePlotting.cpp.

157 {
158  MACH3LOG_INFO("Starting {}", __func__);
159  auto Titles = Get<std::vector<std::string>>(Settings["FileTitle"], __FILE__, __LINE__);
160  std::vector<std::vector<double>> FlucDrawVec(InputFiles.size());
161  // KS: Alternatively try "_drawfluc_draw"
162  std::string FlucutationType = "_predfluc_draw";
163  //KS: P-values per each sample
164  std::cout<<"\\begin{table}[htb]"<<std::endl;
165  std::cout<<"\\centering"<<std::endl;
166  std::cout<<"\\begin{tabular}{ | l | ";
167 
168  for(unsigned int f = 0; f < InputFiles.size(); f++)
169  {
170  std::cout<<"c | ";
171  }
172 
173  std::cout<<"} \\hline"<<std::endl;
174  std::cout<<"Sample ";
175  for(unsigned int f = 0; f < InputFiles.size(); f++)
176  {
177  std::cout<<"& \\multicolumn{1}{| c |}{" + Titles[f] +" p-value} ";
178  }
179  std::cout<<"\\\\"<<std::endl;
180  for(unsigned int f = 0; f < InputFiles.size(); f++)
181  {
182  std::cout<<" & Fluctuation of Prediction ";
183  }
184  std::cout<<"\\\\ \\hline"<<std::endl;
185  for(unsigned int i = 0; i < SampleNames.size(); i++)
186  {
187  std::cout<<SampleNames[i];
188  for(unsigned int f = 0; f < InputFiles.size(); f++)
189  {
190  std::string TempString = "Predictive/" + SampleNames[i]+"/"+SampleNames[i] + FlucutationType;
191  TH2D *hist2D = InputFiles[f]->Get<TH2D>(TempString.c_str());
192  double FlucDraw = GetPValue(hist2D);
193  std::cout<<" & "<<FlucDraw;
194  FlucDrawVec[f].push_back(FlucDraw);
195  }
196  std::cout<<" \\\\"<<std::endl;
197  }
198  std::cout<<"Total ";
199  for(unsigned int f = 0; f < InputFiles.size(); f++)
200  {
201  TH2D *hFlucPred = InputFiles[f]->Get<TH2D>(("Predictive/Total/Total" + FlucutationType).c_str());
202  double FlucDraw = GetPValue(hFlucPred);
203  std::cout<<" & "<<FlucDraw;
204  }
205  std::cout<<" \\\\ \\hline"<<std::endl;
206  std::cout<<"\\hline"<<std::endl;
207  std::cout<<"\\end{tabular}"<<std::endl;
208  std::cout<<"\\end{table}"<<std::endl;
209 
210  auto Threshold = GetFromManager<double>(Settings["Significance"], 0.05, __FILE__ , __LINE__);
211  for(unsigned int f = 0; f < InputFiles.size(); f++)
212  {
213  MACH3LOG_INFO("Calculating Shape for file {}", Titles[f]);
214 
215  CheckBonferoniCorrectedpValue(SampleNames, FlucDrawVec[f], Threshold);
216  MACH3LOG_INFO("Combined pvalue following Fisher method: {:.4f}", FisherCombinedPValue(FlucDrawVec[f]));
217  }
218 }
double GetPValue(const TH2D *hist)
double FisherCombinedPValue(const std::vector< double > &pvalues)
KS: Combine p-values using Fisher's method.
void CheckBonferoniCorrectedpValue(const std::vector< std::string > &SampleNameVec, const std::vector< double > &PValVec, const double Threshold)
KS: For more see https://www.t2k.org/docs/technotes/429/TN429_v8#page=63.

◆ PrintPredictiveLLH()

void PrintPredictiveLLH ( const std::vector< TFile * > &  InputFiles,
const std::vector< std::string > &  SampleNames 
)

KS Print Predictive LLH into Latex table format.

Definition at line 684 of file PredictivePlotting.cpp.

685  {
686  MACH3LOG_INFO("Starting {}", __func__);
687  MACH3LOG_INFO("");
688 
689  std::vector<double> Total(InputFiles.size());
690  //KS: We now prepare to make tables for TN etc.
691  std::cout<<"\\begin{table}[htb]"<<std::endl;
692  std::cout<<"\\centering"<<std::endl;
693  std::cout<<"\\begin{tabular}{ | l |";
694  for(unsigned int f = 0; f < InputFiles.size(); f++)
695  {
696  Total[f] = 0.;
697  std::cout<<" c |";
698  }
699  std::cout<<"} \\hline"<<std::endl;
700  std::cout<<"Sample ";
701  for(unsigned int f = 0; f < InputFiles.size(); f++)
702  {
703  std::cout<<"& 2#log#mathcal{L}_{stat} ";
704  }
705  std::cout<<"\\\\ \\hline"<<std::endl;
706  for(unsigned int i = 0; i < SampleNames.size(); i++)
707  {
708  std::cout<<SampleNames[i];
709  std::string TempString = "Predictive/" + SampleNames[i]+"/"+SampleNames[i]+"_mc_PostPred";
710  for(unsigned int f = 0; f < InputFiles.size(); f++)
711  {
712  TH1 *hist = static_cast<TH1*>(InputFiles[f]->Get(TempString.c_str()));
713 
714  double llh = GetLLH(hist);
715  std::cout<<" & "<<llh;
716  Total[f] += llh;
717  }
718  std::cout<<" \\\\"<<std::endl;
719  }
720  std::cout<<"Total";
721  for(unsigned int f = 0; f < InputFiles.size(); f++) {
722  std::cout<<" & "<<Total[f];
723  }
724  std::cout<<" \\\\"<<std::endl;
725  std::cout<<"\\hline"<<std::endl;
726  std::cout<<"\\end{tabular}"<<std::endl;
727  std::cout<<"\\end{table}"<<std::endl;
728  std::cout<<" "<<std::endl;
729 }
double GetLLH(TH1 *hist)

Variable Documentation

◆ PlotMan

M3::Plotting::PlottingManager* PlotMan = nullptr
Warning
KS: keep raw pointer or ensure manual delete of PlotMan. If spdlog in automatically deleted before PlotMan then destructor has some spdlog and this could cause segfault

Definition at line 15 of file PredictivePlotting.cpp.