MaCh3  2.4.2
Reference Guide
GetPostfitParamPlots.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <sstream>
3 #include <iomanip>
4 #include <algorithm>
5 
6 #include "PlottingUtils/PlottingUtils.h"
7 #include "PlottingUtils/PlottingManager.h"
8 
10 #include "TROOT.h"
11 #include "TGaxis.h"
12 #include "TString.h"
13 #include "TStyle.h"
14 #include "TH1.h"
15 #include "TH2.h"
16 #include "TF1.h"
17 #include "TLegend.h"
18 #include "TPad.h"
19 #include "TCanvas.h"
20 #include "TTree.h"
21 #include "TFile.h"
22 #include "TVectorD.h"
23 #include "TCandle.h"
24 #include "TFrame.h"
25 #include "TGraphAsymmErrors.h"
27 
46 
47 //this file has lots of usage of the ROOT plotting interface that only takes floats, turn this warning off for this CU for now
48 #pragma GCC diagnostic ignored "-Wfloat-conversion"
49 #pragma GCC diagnostic ignored "-Wconversion"
50 
51 MaCh3Plotting::PlottingManager *PlotMan;
52 TH1D *Prefit;
53 
56 
57 std::vector<int> NDSamplesBins;
58 std::vector<std::string> NDSamplesNames;
59 int nBins;
60 TCanvas *canv;
61 
62 std::string SaveName;
63 
64 TPad *p3;
65 TPad *p4;
66 
67 // KS: Color for 0 - prefit, 1 postfit, 2 another postfit, 3 you know the drill
68 constexpr Color_t PlotColor[] = {kRed, kBlack, kBlue, kGreen};
69 std::string plotType;
70 
71 void copyParToBlockHist(const int localBin, const std::string& paramName, TH1D* blockHist,
72  const std::string& type, const int fileId, const bool setLabels = true){
73  // Set the values in the sub-histograms
74  MACH3LOG_DEBUG("copying data from at local bin {}: for parameter {}", localBin, paramName);
75  MACH3LOG_DEBUG(" Fitter specific name: {}", PlotMan->input().translateName(fileId, MaCh3Plotting::kPostFit, paramName));
76  MACH3LOG_DEBUG(" value: {:.4f}", PlotMan->input().getPostFitValue(fileId, paramName, type));
77  MACH3LOG_DEBUG(" error: {:.4f}", PlotMan->input().getPostFitError(fileId, paramName, type));
78 
79  blockHist->SetBinContent(localBin +1, PlotMan->input().getPostFitValue(fileId, paramName, type));
80  blockHist->SetBinError(localBin +1, PlotMan->input().getPostFitError(fileId, paramName, type));
81 
82  if(setLabels){
83  blockHist->GetXaxis()->SetBinLabel(localBin +1, paramName.c_str());
84  blockHist->GetXaxis()->LabelsOption("v");
85  }
86 }
87 
88 inline void InitializePads(TCanvas* canvas, TPad*& pad3, TPad*& pad4) {
89  // Initialize TPad p3
90  pad3 = new TPad("Top", "Top", 0.0, 0.4, 1.0, 1.0);
91  pad3->SetLeftMargin(canvas->GetLeftMargin());
92  pad3->SetRightMargin(canvas->GetRightMargin());
93  pad3->SetTopMargin(canvas->GetTopMargin());
94  pad3->SetBottomMargin(0);
95  pad3->SetGrid();
96 
97  // Initialize TPad p4
98  pad4 = new TPad("Bottom", "Bottom", 0.0, 0.0, 1.0, 0.4);
99  pad4->SetLeftMargin(canvas->GetLeftMargin());
100  pad4->SetRightMargin(canvas->GetRightMargin());
101  pad4->SetTopMargin(0);
102  pad4->SetBottomMargin(0.75);
103  pad4->SetGrid();
104 }
105 
106 void CopyViolinToBlock(TH2D* FullViolin, TH2D* ReducedViolin, const std::vector<std::string>& ParamNames) {
107  for(unsigned int i = 0; i < ParamNames.size(); i++)
108  {
109  int ParamBinId = M3::_BAD_INT_;
110  for (int ix = 0; ix < FullViolin->GetXaxis()->GetNbins(); ++ix) {
111  if(FullViolin->GetXaxis()->GetBinLabel(ix+1) == ParamNames[i])
112  {
113  ParamBinId = ix+1;
114  break;
115  }
116  }
117  if(ParamBinId == M3::_BAD_INT_) {
118  MACH3LOG_WARN("Didn't find param {}", ParamNames[i]);
119  continue;
120  }
121  //KS Fill content of reduced violin
122  for (int iy = 0; iy < FullViolin->GetYaxis()->GetNbins(); ++iy) {
123  ReducedViolin->SetBinContent(i+1, iy+1, FullViolin->GetBinContent(ParamBinId, iy+1));
124  ReducedViolin->GetXaxis()->SetBinLabel(i+1, ParamNames[i].c_str());
125  }
126  }
127  ReducedViolin->SetFillColor(FullViolin->GetFillColor());
128  ReducedViolin->SetFillColorAlpha(FullViolin->GetMarkerColor(), 0.35);
129  ReducedViolin->SetLineColor(FullViolin->GetMarkerColor());
130 
131  ReducedViolin->SetMarkerColor(FullViolin->GetMarkerColor());
132  ReducedViolin->SetMarkerStyle(FullViolin->GetMarkerStyle());
133  ReducedViolin->SetMarkerSize(FullViolin->GetMarkerSize());
134 
135  ReducedViolin->GetYaxis()->SetTitleOffset(FullViolin->GetTitleOffset());
136  ReducedViolin->GetYaxis()->SetTitle(FullViolin->GetYaxis()->GetTitle());
137  ReducedViolin->GetXaxis()->LabelsOption("v");
138 }
139 
140 void PrettifyTitles(TH1D *Hist) {
141  for (int i = 0; i < Hist->GetXaxis()->GetNbins(); ++i)
142  {
143  std::string title = Hist->GetXaxis()->GetBinLabel(i+1);
144  title = PlotMan->style().prettifyParamName(title);
145  Hist->GetXaxis()->SetBinLabel(i+1, title.c_str());
146  }
147 }
148 
149 void PrettifyTitles(TH2D *Hist) {
150  for (int i = 0; i < Hist->GetXaxis()->GetNbins(); ++i)
151  {
152  std::string title = Hist->GetXaxis()->GetBinLabel(i+1);
153 
154  title = PlotMan->style().prettifyParamName(title);
155  Hist->GetXaxis()->SetBinLabel(i+1, title.c_str());
156  }
157 }
158 
159 bool ReadSettings(const std::shared_ptr<TFile>& File1)
160 {
161  MACH3LOG_DEBUG("Reading settings for file {}", File1->GetName());
162  #ifdef DEBUG
163  File1->ls();
164  #endif
165  TTree *Settings = (File1->Get<TTree>("Settings"));
166 
167  // can't find settings tree :(
168  if (!Settings) return false;
169 
170  MACH3LOG_DEBUG("Got settings tree");
171  MaCh3Utils::Print(Settings);
172 
173  Settings->SetBranchAddress("NDParameters", &NDParameters);
174  Settings->SetBranchAddress("NDParametersStartingPos", &NDParametersStartingPos);
175 
176  std::vector<int> *NDSamples_Bins = 0;
177  std::vector<std::string> *NDSamples_Names = 0;
178  Settings->SetBranchAddress("NDSamplesNames", &NDSamples_Names);
179  Settings->SetBranchAddress("NDSamplesBins", &NDSamples_Bins);
180 
181  Settings->GetEntry(0);
182 
183  NDSamplesNames = *NDSamples_Names;
184  NDSamplesBins = *NDSamples_Bins;
185 
186  MACH3LOG_DEBUG("Read settings tree successfully");
187  return true;
188 }
189 
190 std::unique_ptr<TH1D> makeRatio(TH1D *PrefitCopy, TH1D *PostfitCopy, bool setAxes) {
191  // set up the ratio hist
192  std::unique_ptr<TH1D> Ratio = M3::Clone(PrefitCopy);
193  Ratio->GetYaxis()->SetTitle("(x_{Post}-#mu_{Prior})/#sigma_{Prior}");
194  Ratio->SetMinimum(-3.7);
195  Ratio->SetMaximum(3.7);
196 
197  for (int j = 0; j < Ratio->GetXaxis()->GetNbins(); ++j)
198  {
199  if ( PrefitCopy->GetBinError(j+1) > 1.e-5 )
200  {
201  Ratio->SetBinContent(j+1, (PostfitCopy->GetBinContent(j+1)-PrefitCopy->GetBinContent(j+1))/PrefitCopy->GetBinError(j+1));
202 
203  double up = (PostfitCopy->GetBinContent(j+1)+PostfitCopy->GetBinError(j+1)-PrefitCopy->GetBinContent(j+1))/PrefitCopy->GetBinError(j+1);
204  double down = (PostfitCopy->GetBinContent(j+1)-PostfitCopy->GetBinError(j+1)-PrefitCopy->GetBinContent(j+1))/PrefitCopy->GetBinError(j+1);
205 
206  double maximum = up-Ratio->GetBinContent(j+1);
207  double minimum = Ratio->GetBinContent(j+1)-down;
208 
209  Ratio->SetBinError(j+1, std::max(maximum, minimum));
210  }
211  //KS: Most likely flat prior
212  else {
213  Ratio->SetBinContent(j+1, (PostfitCopy->GetBinContent(j+1)-PrefitCopy->GetBinContent(j+1)));
214 
215  double up = (PostfitCopy->GetBinContent(j+1)+PostfitCopy->GetBinError(j+1)-PrefitCopy->GetBinContent(j+1));
216  double down = (PostfitCopy->GetBinContent(j+1)-PostfitCopy->GetBinError(j+1)-PrefitCopy->GetBinContent(j+1));
217 
218  double maximum = up-Ratio->GetBinContent(j+1);
219  double minimum = Ratio->GetBinContent(j+1)-down;
220 
221  Ratio->SetBinError(j+1, std::max(maximum, minimum));
222  }
223  } //end loop over parameters
224 
225  if(setAxes){
226  Ratio->SetFillStyle(0);
227  Ratio->SetFillColor(0);
228 
229  Ratio->SetLineColor(PostfitCopy->GetLineColor());
230  if (Ratio->GetLineColor() == 0) Ratio->SetLineColor(kBlack);
231  Ratio->SetMarkerColor(PostfitCopy->GetMarkerColor());
232 
233  Ratio->SetLineWidth(PlotMan->getOption<int>("plotLineWidth"));
234  Ratio->SetTitle("");
235 
236  Ratio->SetMarkerSize(2);
237  Ratio->SetMarkerStyle(20);
238 
239  Ratio->GetYaxis()->SetTitleSize(25);
240  Ratio->GetYaxis()->SetTitleFont(43);
241  Ratio->GetYaxis()->SetTitleOffset(2.0);
242  Ratio->GetYaxis()->SetLabelFont(43);
243  Ratio->GetYaxis()->SetLabelSize(25);
244  Ratio->GetYaxis()->CenterTitle();
245  Ratio->GetYaxis()->SetNdivisions(5,2,0);
246 
247  Ratio->GetXaxis()->SetTitleSize(25);
248  Ratio->GetXaxis()->SetTitleFont(43);
249  Ratio->GetXaxis()->SetTitleOffset(4.0);
250  Ratio->GetXaxis()->SetLabelOffset(0.025);
251  Ratio->GetXaxis()->SetLabelFont(43);
252  Ratio->GetXaxis()->SetLabelSize(25);
253  }
254 
255  return Ratio;
256 }
257 
258 void DrawPlots(TCanvas *plotCanv, TH1D* PrefitCopy, const std::vector<std::unique_ptr<TH1D>>& PostfitVec, TPad *mainPad, TPad *ratioPad) {
259  // Draw!
260  plotCanv->cd();
261  mainPad->Draw();
262  mainPad->cd();
263  PrefitCopy->GetYaxis()->SetTitle("Parameter Value");
264 
265  PrefitCopy->GetYaxis()->SetLabelSize(0.);
266  PrefitCopy->GetYaxis()->SetTitleSize(0.05);
267  PrefitCopy->GetYaxis()->SetTitleOffset(1.3);
268  PrefitCopy->Draw("e2");
269 
270  for (int fileId = 0; fileId < static_cast<int>(PostfitVec.size()); fileId++) {
271  TH1D *postFitHist = PostfitVec[fileId].get();
272 
273  postFitHist->SetMarkerColor(TColor::GetColorPalette(fileId));
274  postFitHist->SetLineColor(TColor::GetColorPalette(fileId));
275  postFitHist->SetMarkerStyle(7);
276  postFitHist->SetLineStyle(1+fileId);
277  postFitHist->SetLineWidth(PlotMan->getOption<int>("plotLineWidth"));
278 
279  postFitHist->Draw("e1, same");
280  }
281 
282  plotCanv->Update();
283  auto axis = std::make_unique<TGaxis>(PrefitCopy->GetXaxis()->GetBinLowEdge(PrefitCopy->GetXaxis()->GetFirst()), gPad->GetUymin()+0.01,
284  PrefitCopy->GetXaxis()->GetBinLowEdge(PrefitCopy->GetXaxis()->GetFirst()), gPad->GetUymax(),
285  gPad->GetUymin()+0.01, gPad->GetUymax(), 510, "");
286  axis->SetLabelFont(43);
287  axis->SetLabelSize(25);
288  axis->Draw();
289 
290  plotCanv->cd();
291  ratioPad->Draw();
292  ratioPad->cd();
293 
294  std::vector<std::unique_ptr<TH1D>> ratioHists;
295  // save pointers to these so we can delete them once we are done
296  ratioHists.push_back(makeRatio(PrefitCopy, PostfitVec[0].get(), true));
297 
298  ratioHists[0]->Draw("p");
299  for(int postFitIdx = 1; postFitIdx < static_cast<int>(PostfitVec.size()); postFitIdx++){
300  ratioHists.push_back(makeRatio(PrefitCopy, PostfitVec[postFitIdx].get(), true));
301 
302  ratioHists[postFitIdx]->SetMarkerColor(TColor::GetColorPalette(postFitIdx));
303  ratioHists[postFitIdx]->SetLineColor(TColor::GetColorPalette(postFitIdx));
304  ratioHists[postFitIdx]->SetMarkerStyle(7);
305  ratioHists[postFitIdx]->SetLineStyle(1+postFitIdx);
306  ratioHists[postFitIdx]->SetLineWidth(PlotMan->getOption<int>("plotLineWidth"));
307 
308  ratioHists[postFitIdx]->Draw("p same");
309  }
310 
311  // draw lines across the plot at +-1 and 0
312  TLine line(ratioHists[0]->GetXaxis()->GetBinLowEdge(ratioHists[0]->GetXaxis()->GetFirst()), 0.0, ratioHists[0]->GetXaxis()->GetBinLowEdge(ratioHists[0]->GetXaxis()->GetLast()+1), 0.0);
313  TLine line2(ratioHists[0]->GetXaxis()->GetBinLowEdge(ratioHists[0]->GetXaxis()->GetFirst()), 1.0, ratioHists[0]->GetXaxis()->GetBinLowEdge(ratioHists[0]->GetXaxis()->GetLast()+1), 1.0);
314  TLine line3(ratioHists[0]->GetXaxis()->GetBinLowEdge(ratioHists[0]->GetXaxis()->GetFirst()), -1.0, ratioHists[0]->GetXaxis()->GetBinLowEdge(ratioHists[0]->GetXaxis()->GetLast()+1), -1.0);
315 
316  line.SetLineColor(kRed);
317  line.SetLineStyle(kDashed);
318  line.SetLineWidth(PlotMan->getOption<int>("refLineWidth"));
319  line2.SetLineColor(kRed);
320  line2.SetLineStyle(kDashed);
321  line2.SetLineWidth(PlotMan->getOption<int>("refLineWidth"));
322  line3.SetLineColor(kRed);
323  line3.SetLineStyle(kDashed);
324  line3.SetLineWidth(PlotMan->getOption<int>("refLineWidth"));
325 
326  line.Draw("same");
327  line2.Draw("same");
328  line3.Draw("same");
329 
330  plotCanv->Print((SaveName).c_str());
331 }
332 
334 {
335  // get the names of the blocks of parameters to group together
336  std::vector<std::string> const blockNames = PlotMan->getOption<std::vector<std::string>>("paramGroups");
337  const int nPlots = static_cast<int>(blockNames.size());
338 
339  for (int i = 0; i < nPlots; i++)
340  {
341  // get the configuration for this parameter
342  std::string blockName = blockNames[i];
343  YAML::Node paramBlock = PlotMan->getOption(blockName);
344  auto blockTitle = paramBlock[0].as<std::string>();
345  auto blockLimits = paramBlock[1].as<std::vector<double>>();
346  auto blockContents = paramBlock[2].as<std::vector<std::string>>();
347 
348  // get num of params in the block
349  const int nParams = static_cast<int>(blockContents.size());
350 
351  // set some plot things
352  auto blockHist_prefit = std::make_unique<TH1D>(blockName.c_str(), blockTitle.c_str(),
353  nParams, 0.0, static_cast<double>(nParams));
354 
355  PlotMan->style().setTH1Style(blockHist_prefit.get(), PlotMan->getOption<std::string>("prefitHistStyle"));
356 
357  // set the errors for the prefit block hist
358  for(int localBin=0; localBin < nParams; localBin ++){
359  // the "local" bin is the params index within the group of parameters
360  std::string paramName = blockContents[localBin];
361  copyParToBlockHist(localBin, paramName, blockHist_prefit.get(), "Prior", 0);
362  }
363 
364  // now set for the postfit blocks for all files
365  std::vector <std::unique_ptr<TH1D>> blockHist_postfit_Vec;
366  for(unsigned int fileId = 0; fileId < PlotMan->getNFiles(); fileId++) {
367  auto blockHist_postfit = std::make_unique<TH1D>((blockName + PlotMan->getFileName(fileId)).c_str(),
368  blockTitle.c_str(), nParams, 0.0, static_cast<double>(nParams));
369 
370  // loop through all the parameters in this block and set the contents in the blocks TH1
371  for(int localBin=0; localBin < nParams; localBin ++){
372  // the "local" bin is the params index within the group of parameters
373  std::string paramName = blockContents[localBin];
374  copyParToBlockHist(localBin, paramName, blockHist_postfit.get(), "", fileId);
375  }
376  blockHist_postfit_Vec.push_back(std::move(blockHist_postfit));
377  }
378  // set the y axis limits we got from config
379  blockHist_prefit->GetYaxis()->SetRangeUser(blockLimits[0], blockLimits[1]);
380 
381  // Do some fancy replacements
382  PrettifyTitles(blockHist_prefit.get());
383 
384  DrawPlots(canv, blockHist_prefit.get(), blockHist_postfit_Vec, p3, p4);
385  }
386 }
387 
389 {
390  // these for non named params where we don't need as much space
391  auto p1 = std::make_unique<TPad>("p1", "p1", 0.0, 0.3, 1.0, 1.0);
392  auto p2 = std::make_unique<TPad>("p2", "p2", 0.0, 0.0, 1.0, 0.3);
393  p1->SetLeftMargin(canv->GetLeftMargin());
394  p1->SetRightMargin(canv->GetRightMargin());
395  p1->SetTopMargin(canv->GetTopMargin());
396  p1->SetBottomMargin(0);
397  p1->SetGrid();
398  p2->SetLeftMargin(canv->GetLeftMargin());
399  p2->SetRightMargin(canv->GetRightMargin());
400  p2->SetTopMargin(0);
401  p2->SetBottomMargin(0.25);
402  p2->SetGrid();
403 
404  p1->SetLogx(true);
405  p2->SetLogx(true);
406 
407  // get the names of the blocks of parameters to group together
408  std::vector<std::string> const fluxBlockNames = PlotMan->getOption<std::vector<std::string>>("fluxGroups");
409  auto const fluxBinningTable = PlotMan->getOption("FluxBinning");
410 
411  const int FluxPlots = static_cast<int>(fluxBlockNames.size());
412 
413  for (int i = 0; i < FluxPlots; i++)
414  {
415  // get the configuration for this block
416  std::string fluxBlockName = fluxBlockNames[i];
417  YAML::Node paramBlock = PlotMan->getOption(fluxBlockName);
418  std::string blockTitle = paramBlock[0].as<std::string>();
419  std::vector<double> blockLimits = paramBlock[1].as<std::vector<double>>();
420  std::string blockBinningName = paramBlock[2].as<std::string>();
421  std::vector<int> blockContents = paramBlock[3].as<std::vector<int>>();
422 
423  // get the binning for this block of flux params
424  std::vector<double> binning = fluxBinningTable[blockBinningName].as<std::vector<double>>();
425 
426  // get num of params in the block
427  int nParams = blockContents[1] - blockContents[0] +1;
428  // check for sanity
429  if(nParams <= 0 || blockContents.size() > 2){
430  MACH3LOG_CRITICAL("Invalid flux parameter block endpoints specified for {}", fluxBlockName);
431  MACH3LOG_CRITICAL(" Should have the form [<low index>, <up index>]");
432  throw MaCh3Exception(__FILE__ , __LINE__ );
433  }
434  if (nParams != static_cast<int>(binning.size()) - 1) {
435  MACH3LOG_CRITICAL("Binning provided for flux param block {} does not match the number of parameters specified for the block", fluxBlockName);
436  MACH3LOG_CRITICAL(" Provided {} parameters but {} bins", nParams, binning.size() -1);
437  throw MaCh3Exception(__FILE__ , __LINE__ );
438  }
439 
440  auto blockHist_prefit = std::make_unique<TH1D>(fluxBlockName.c_str(), blockTitle.c_str(), nParams, binning.data());
441  blockHist_prefit->GetYaxis()->SetTitle("Parameter Variation");
442  blockHist_prefit->GetXaxis()->SetTitle("E_{#nu} (GeV)");
443  blockHist_prefit->GetXaxis()->SetTitleOffset(blockHist_prefit->GetXaxis()->GetTitleOffset()*1.2);
444  PlotMan->style().setTH1Style(blockHist_prefit.get(), PlotMan->getOption<std::string>("prefitHistStyle"));
445  // set the errors for the prefit block hist
446  for(int fluxParId = blockContents[0]; fluxParId <= blockContents[1]; fluxParId++){
447  int localBin = fluxParId - blockContents[0];
448  std::string paramName = "b_" + std::to_string(fluxParId);
449  copyParToBlockHist(localBin, paramName, blockHist_prefit.get(), "Prior", 0, false);
450  }
451 
452  // now set for the postfit blocks for all files
453  std::vector <std::unique_ptr<TH1D>> blockHist_postfit_Vec;
454  for(unsigned int fileId = 0; fileId < PlotMan->getNFiles(); fileId++){
455  auto blockHist_postfit = std::make_unique<TH1D>(fluxBlockName.c_str(), blockTitle.c_str(), nParams, binning.data());
456 
457  for(int fluxParId = blockContents[0]; fluxParId <= blockContents[1]; fluxParId++){
458  int localBin = fluxParId - blockContents[0];
459  std::string paramName = "b_" + std::to_string(fluxParId);
460 
461  copyParToBlockHist(localBin, paramName, blockHist_postfit.get(), "", fileId, false);
462  }
463  blockHist_postfit_Vec.push_back(std::move(blockHist_postfit));
464  }
465  // set the y axis limits we got from config
466  blockHist_prefit->GetYaxis()->SetRangeUser(blockLimits[0], blockLimits[1]);
467  DrawPlots(canv, blockHist_prefit.get(), blockHist_postfit_Vec, p1.get(), p2.get());
468  }
469 
470  canv->cd();
471  canv->SetLogx(false);
472  canv->SetBottomMargin(canv->GetBottomMargin()*1.7);
473 }
474 
477 {
478  MACH3LOG_INFO("ND detector parameters: {}", NDParameters);
479  Prefit->GetYaxis()->SetTitleOffset(Prefit->GetYaxis()->GetTitleOffset()*1.2);
480 
481  TPad* pTop = nullptr;
482  TPad* pDown = nullptr;
483  InitializePads(canv, pTop, pDown);
484 
485  int NDbinCounter = NDParametersStartingPos;
486  int Start = NDbinCounter;
487 
488  MACH3LOG_INFO("Running on {} samples", NDSamplesNames.size());
489 
490  for (unsigned int i = 0; i < NDSamplesNames.size(); ++i)
491  {
492  MACH3LOG_DEBUG("--- On sample {}", NDSamplesNames[i]);
493  NDbinCounter += NDSamplesBins[i];
494 
495  std::vector<std::unique_ptr<TH1D>> PostfitNDDetHistVec(PlotMan->getNFiles());
496  TH1D *PreFitNDDetHist = PlotMan->input().getFile(0).file->Get<TH1D>(Form("param_%s_prefit", NDSamplesNames[i].c_str()));
497  PlotMan->style().setTH1Style(PreFitNDDetHist, PlotMan->getOption<std::string>("prefitHistStyle"));
498 
499  std::string temp = NDSamplesNames[i].c_str();
500  while (temp.find("_") != std::string::npos) {
501  temp.replace(temp.find("_"), 1, std::string(" "));
502  }
503  PreFitNDDetHist->SetTitle(temp.c_str());
504  PreFitNDDetHist->GetXaxis()->SetRangeUser(Start, NDbinCounter);
505 
506  MACH3LOG_DEBUG(" Start bin: {} :: End bin: {}", Start, NDbinCounter);
507  // set the x range for the postfits
508  for(unsigned int fileId = 0; fileId < PlotMan->getNFiles(); fileId++){
509  PostfitNDDetHistVec[fileId] = M3::Clone(PlotMan->input().getFile(fileId).file->Get<TH1D>(Form("param_%s_%s", NDSamplesNames[i].c_str(), plotType.c_str())));
510  }
511 
512  //KS: We dont' need name for every nd param
513  for(int j = 0; j < NDSamplesBins[i]; ++j)
514  {
515  bool ProductOfTen = false;
516  if(j % 10) ProductOfTen = true;
517  if(j != 0 && ProductOfTen) PreFitNDDetHist->GetXaxis()->SetBinLabel(Start+j+1, " ");
518  else{
519  PreFitNDDetHist->GetXaxis()->SetBinLabel(Start+j+1, Form("Det Variation Bin %i", Start+j));
520  }
521  }
522 
523  PreFitNDDetHist->GetYaxis()->SetRangeUser(PlotMan->getOption<double>("detParYRange_low"), PlotMan->getOption<double>("detParYRange_high"));
524 
525  Start += NDSamplesBins[i];
526 
527  DrawPlots(canv, PreFitNDDetHist, PostfitNDDetHistVec, pTop, pDown);
528  canv->Update();
529  }
530  delete pTop;
531  delete pDown;
532 }
533 
535 {
536  gStyle->SetPalette(51);
537 
538  auto blankCanv = std::make_unique<TCanvas>("blankCanv", "blankCanv", 2048, 2048);
539  blankCanv->SaveAs("RidgePlots.pdf[");
540 
541  // get the names of the blocks of parameters to group together
542  const auto blockNames = PlotMan->getOption<std::vector<std::string>>("paramGroups");
543  const int nPlots = static_cast<int>(blockNames.size());
544 
545  constexpr double padTopMargin = 0.9;
546  constexpr double padBottomMargin = 0.1;
547  constexpr double padOverlap = 0.9;
548  constexpr double ridgeLineWidth = 1.0;
549  for (int i = 0; i < nPlots; i++)
550  {
551  // get the configuration for this parameter
552  std::string blockName = blockNames[i];
553  auto const &paramBlock = PlotMan->getOption(blockName);
554  auto blockTitle = paramBlock[0].as<std::string>();
555  auto blockLimits = paramBlock[1].as<std::vector<double>>();
556  auto blockContents = paramBlock[2].as<std::vector<std::string>>();
557 
558  // the directory of histograms
559  TDirectoryFile *posteriorDir = PlotMan->input().getFile(0).file->Get<TDirectoryFile>("Post_1d_hists");
560 
561  // get num of params in the block
562  int nParams = static_cast<int>(blockContents.size());
563 
564  if (nParams == 1) {
565  MACH3LOG_WARN("{} doesn't work for single param", __func__);
566  continue;
567  }
568  auto ridgeCanv = std::make_unique<TCanvas>("RidgePlotCanv", "RidgePlotCanv", 2048, 2048);
569  ridgeCanv->Divide(1,1+nParams, 0.01, 0.0);
570 
571  auto title = std::make_unique<TLatex>();
572  title->SetTextAlign(21);
573  title->SetTextSize(0.03);
574  title->DrawLatex(0.5, 0.95, blockTitle.c_str());
575 
576  auto label = std::make_unique<TLatex>();
577  label->SetTextAlign(31);
578  label->SetTextSize(0.02);
579 
580  auto line = std::make_unique<TLine>();
581  line->SetLineColor(kBlack);
582  line->SetLineWidth(ridgeLineWidth);
583 
584  // use this to set the limits and also to plot the x axis and grid
585  auto axisPlot = std::make_unique<TH1D>("axis plot", "", 1, blockLimits[0], blockLimits[1]);
586 
587  std::vector<std::unique_ptr<TH1D>> axisPlot_holder(nParams);
588  std::vector<std::unique_ptr<TPad>> graph_holder(nParams);
589  for(int parId = 0; parId < nParams; parId++) {
590  std::string paramName = blockContents[parId];
591 
592  TH1D *posteriorDist = nullptr;
593  // get the list of objects in the directory
594  TIter next(posteriorDir->GetListOfKeys());
595  while (TKey* key = static_cast<TKey*>(next())) {
596  // check if the end of the param name matches with the MaCh3 name, do this so we exclude things like nds_ at the start of the name
597  std::string str(key->GetTitle());
598  std::string name = PlotMan->input().translateName(0, MaCh3Plotting::kPostFit, paramName);
599  uint pos = str.find(name);
600  bool foundPar = (pos == str.length() - name.length());
601 
602  MACH3LOG_TRACE("Looking for {} in {}", name, str);
603  if(foundPar){
604  MACH3LOG_TRACE("Found it");
605  posteriorDist = posteriorDir->Get<TH1D>(key->GetName());
606  }
607  }
608 
609  if(posteriorDist == nullptr){
610  MACH3LOG_WARN("Couldn't find parameter {} when making ridgeline plots", paramName);
611  MACH3LOG_WARN("It could be fixed param");
612  continue;
613  }
614 
615  // EM: do some funky scaling so that we always get evenly spaced pads in the range [bottomMargin, TopMargin] with the specified overlap
616  double padAnchor = padBottomMargin + (static_cast<double>(nParams - parId - 1) /
617  static_cast<double>(nParams - 1)) * (padTopMargin - padBottomMargin);
618  double padWidth = (padTopMargin - padBottomMargin) / static_cast<double>(nParams);
619  double norm = (padTopMargin - padBottomMargin);
620 
621  double padTop = padWidth * (1.0 + padOverlap) * (padTopMargin - padAnchor) / norm + padAnchor;
622  double padBottom = padAnchor - padWidth * (1.0 + padOverlap) * (padAnchor - padBottomMargin) / norm;
623 
624  auto pad = std::make_unique<TPad>(paramName.c_str(), "", 0.3, padBottom, 0.9, padTop, -1, 0, -1);
625  ridgeCanv->cd();
626 
627  pad->SetBottomMargin(0.0);
628  pad->SetTopMargin(0.0);
629  pad->SetLeftMargin(0.0);
630  pad->SetRightMargin(0.0);
631 
632  pad->Draw();
633  pad->cd();
634  pad->SetFillStyle(4000);
635 
636  gPad->SetFrameFillStyle(4000);
637  posteriorDist->GetFunction("Gauss")->SetBit(TF1::kNotDraw);
638  posteriorDist->SetTitle("");
639  posteriorDist->SetLineWidth(ridgeLineWidth);
640 
641  auto axisPlot_tmp = M3::Clone(axisPlot.get(), Form("AxisPlot_%s", paramName.c_str()));
642  axisPlot_tmp->Draw("A");
643  posteriorDist->Draw("H SAME");
644 
645  axisPlot_tmp->GetYaxis()->SetRangeUser(0.0, 0.7 *posteriorDist->GetMaximum());
646  posteriorDist->SetLineColor(kWhite);
647  posteriorDist->SetFillColorAlpha(TColor::GetColorPalette(floor(static_cast<float>(parId) *
648  TColor::GetNumberOfColors() / static_cast<float>(nParams))), 0.85);
649 
650  posteriorDist->GetXaxis()->SetRangeUser(blockLimits[0], blockLimits[1]);
651  posteriorDist->GetYaxis()->SetTitle(paramName.c_str());
652 
653  //EM: Have to get the frame drawn by the histogram and then set it to be transparent... it took me an hour to get rid of this one line on a plot
654  gPad->Modified(); gPad->Update();
655  TFrame *frame = gPad->GetFrame();
656  frame->SetLineColorAlpha(0, 0.0);
657 
658  ridgeCanv->cd();
659  label->DrawLatexNDC(0.29, padBottom + 0.005, PlotMan->style().prettifyParamName(paramName).c_str());
660  line->DrawLine(0.1, padBottom, 0.9, padBottom);
661 
662  axisPlot_holder[parId] = std::move(axisPlot_tmp);
663  graph_holder[parId] = std::move(pad);
664  }
665 
666  ridgeCanv->cd();
667  ridgeCanv->SetGrid(1,1);
668  auto axisPad = std::make_unique<TPad>("AxisPad", "", 0.3, 0.0, 0.9, 1.0, -1, 0, -1);
669  axisPad->SetLeftMargin(0.0);
670  axisPad->SetRightMargin(0.0);
671  axisPad->Draw();
672  axisPad->cd();
673  axisPad->SetGrid(1,1);
674  axisPad->SetFrameFillStyle(4000);
675 
676  axisPlot->GetXaxis()->SetTickSize(0.01);
677  axisPlot->GetXaxis()->SetTitle("Parameter Variation");
678  axisPlot->GetYaxis()->SetLabelOffset(9999);
679  axisPlot->GetYaxis()->SetLabelSize(0);
680  axisPlot->GetYaxis()->SetTickSize(0);
681  axisPlot->GetYaxis()->SetAxisColor(0,0.0);
682  axisPlot->Draw("AXIS");
683  axisPlot->Draw("AXIG SAME");
684 
685  axisPlot->SetFillStyle(4000);
686  axisPad->SetFillStyle(4000);
687 
688  axisPad->SetGrid(1,1);
689  gPad->Modified(); gPad->Update();
690  gPad->SetFrameFillStyle(4000);
691 
692  gPad->Modified(); gPad->Update();
693  TFrame *frame = gPad->GetFrame();
694  frame->SetLineColorAlpha(0, 0.0);
695 
696  ridgeCanv->SaveAs("RidgePlots.pdf");
697  }
698  blankCanv->SaveAs("RidgePlots.pdf]");
699 }
700 
702 {
703  SaveName = PlotMan->getOutputName();
704 
705  //KS: By default we take HPD values, but by changing "plotType" you can use for example Gauss
706  plotType = "HPD";
707  //plotType = "gaus";
708 
709  MACH3LOG_INFO("Plotting {} errors", plotType);
710 
711  // if we have one MaCh3 nd file then we can get settings from it
712  bool plotNDDet = false;
713  for (size_t fileId = 0; fileId < PlotMan->input().getNInputFiles(); fileId++) {
714  if(!ReadSettings(PlotMan->input().getFile(0).file)) {
715  MACH3LOG_INFO("at least one file provided does not have 'settings' tree indicating it is not MaCh3 ND file");
716  MACH3LOG_INFO(" sadly this means I cannot plot ND Det parameters as this is only supported for MaCh3 ND files for now... sorry :(");
717  plotNDDet = false;
718  }
719  }
720 
721  canv = new TCanvas("canv", "canv", 1024, 1024);
722  //gStyle->SetPalette(51);
723  gStyle->SetOptStat(0); //Set 0 to disable statistic box
724  canv->SetLeftMargin(0.12);
725  canv->SetBottomMargin(0.12);
726  canv->SetTopMargin(0.08);
727  canv->SetRightMargin(0.04);
728 
729  canv->Print((SaveName+"[").c_str());
730 
731  // these for named parameters where we need a nice big gap at the botto to fit the names
733 
734  // Make a Legend page
735  auto leg = std::make_unique<TLegend>(0.0, 0.0, 1.0, 1.0);
736  // make a dummy TH1 to set out legend
737  Prefit = new TH1D();
738  PlotMan->style().setTH1Style(Prefit, PlotMan->getOption<std::string>("prefitHistStyle"));
739  leg->AddEntry(Prefit, "Prior", "lpf");
740 
741  for(unsigned int fileId = 0; fileId < PlotMan->getNFiles(); fileId++){
742  TH1D *postFitHist_tmp = new TH1D();
743  postFitHist_tmp->SetBit(kCanDelete);
744 
745  postFitHist_tmp->SetMarkerColor(TColor::GetColorPalette(fileId));
746  postFitHist_tmp->SetLineColor(TColor::GetColorPalette(fileId));
747  postFitHist_tmp->SetMarkerStyle(7);
748  postFitHist_tmp->SetLineStyle(1+fileId);
749  postFitHist_tmp->SetLineWidth(PlotMan->getOption<int>("plotLineWidth"));
750  leg->AddEntry(postFitHist_tmp, PlotMan->getFileLabel(fileId).c_str(), "lpf");
751  }
752 
753  canv->cd();
754  canv->Clear();
755  leg->Draw();
756  canv->Print((SaveName).c_str());
757 
759 
760  MakeFluxPlots();
761 
762  //KS: By default we don't run ProcessMCMC with PlotDet as this take some time, in case we did let's make fancy plots
763  if(plotNDDet & (NDParameters > 0)) MakeNDDetPlots();
764 
765  canv->Print((SaveName+"]").c_str());
766 
767  MakeRidgePlots();
768 
769  delete canv;
770  delete Prefit;
771 }
772 
773 std::unique_ptr<TGraphAsymmErrors> MakeTGraphAsymmErrors(const std::shared_ptr<TFile>& File, std::vector<int> Index = {})
774 {
775  int GraphBins = Index.size() == 0 ? nBins : Index.size();
776  std::vector<double> x(GraphBins);
777  std::vector<double> y(GraphBins);
778  std::vector<double> exl(GraphBins);
779  std::vector<double> eyl(GraphBins);
780  std::vector<double> exh(GraphBins);
781  std::vector<double> eyh(GraphBins);
782 
783  TH1D* PostHist = static_cast<TH1D*>(File->Get( ("param_xsec_"+plotType).c_str() ));
784 
785  auto Errors_HPD_Positive = static_cast<TVectorD*>(File->Get( "Errors_HPD_Positive" ));
786  auto Errors_HPD_Negative = static_cast<TVectorD*>(File->Get( "Errors_HPD_Negative" ));
787  //KS: I am tempted to multithread this...
788  for(int i = 0; i < GraphBins; ++i)
789  {
790  int Counter = Index.size() == 0 ? i : Index[i];
791  //KS: We are extracting value from three object each having different numbering scheme, I have checked carefully so this is correct please don't change all these +1 +0.5 etc. it just work...
792  x[i] = i + 0.5;
793  y[i] = PostHist->GetBinContent(Counter+1);
794 
795  //KS: We don't want x axis errors as they are confusing in Violin plot
796  exh[i] = 0.00001;
797  exl[i] = 0.00001;
798  eyh[i] = (*Errors_HPD_Positive)(Counter);
799  eyl[i] = (*Errors_HPD_Negative)(Counter);
800  }
801  auto PostGraph = std::make_unique<TGraphAsymmErrors>(GraphBins, x.data(), y.data(), exl.data(), exh.data(), eyl.data(), eyh.data());
802  PostGraph->SetTitle("");
803 
804  return PostGraph;
805 }
806 
809 {
810  //KS: Should be in some config... either way it control whether you plot symmetric or asymmetric error bars
811  bool PlotAssym = true;
812 
813  //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.
814  TCandle::SetScaledViolin(false);
815 
816  std::string OutputName = "";
817  for(unsigned int fileId = 0; fileId < PlotMan->getNFiles(); fileId++){
818  MACH3LOG_INFO("File {}: {} ", fileId, PlotMan->getFileName(fileId));
819  OutputName += PlotMan->getFileName(fileId);
820  OutputName = OutputName.substr(0, OutputName.find(".root"));
821  }
822  MACH3LOG_INFO("Making Violin Plot");
823  OutputName += "_Violin";
824  if(PlotAssym) OutputName += "_Assym";
825 
826  auto canvas = std::make_unique<TCanvas>("canv", "canv", 1024, 1024);
827  canvas->SetGrid();
828  gStyle->SetOptStat(0);
829  //KS: Remove errors on X axis as they are confusing in violin type of plot
830  if(!PlotAssym) gStyle->SetErrorX(0.0001);
831  canvas->SetTickx();
832  canvas->SetTicky();
833  canvas->SetBottomMargin(0.25);
834  canvas->SetTopMargin(0.08);
835  canvas->SetRightMargin(0.03);
836  canvas->SetLeftMargin(0.10);
837  canvas->Print((OutputName+".pdf[").c_str());
838  canvas->SetGrid();
839 
840  if(PlotMan->input().getFile(0).file->Get<TH2D>( "param_violin_prior" ) == nullptr)
841  {
842  MACH3LOG_WARN("Couldn't find violin plot, make sure method from MCMCProcessor is being called");
843  return;
844  }
845  std::unique_ptr<TH2D> ViolinPre = M3::Clone(PlotMan->input().getFile(0).file->Get<TH2D>( "param_violin_prior" ));
846  // Do some fancy replacements
847  ViolinPre->SetFillColor(kRed);
848  ViolinPre->SetFillColorAlpha(kRed, 0.35);
849  ViolinPre->SetMarkerColor(kRed);
850  ViolinPre->SetMarkerStyle(20);
851  ViolinPre->SetMarkerSize(0.5);
852 
853  ViolinPre->GetYaxis()->SetTitleOffset(1.3);
854  ViolinPre->GetYaxis()->SetTitle("Parameter Value");
855  ViolinPre->GetXaxis()->LabelsOption("v");
856 
857  std::unique_ptr<TH1D> Postfit = M3::Clone(PlotMan->input().getFile(0).file->Get<TH1D>( ("param_xsec_"+plotType).c_str() ));
858  Postfit->SetMarkerColor(kRed);
859  Postfit->SetLineColor(kRed);
860  Postfit->SetMarkerStyle(7);
861 
862  std::vector<std::unique_ptr<TH2D>> Violin(PlotMan->getNFiles());
863  for(unsigned int fileId = 0; fileId < PlotMan->getNFiles(); fileId++) {
864  Violin[fileId] = M3::Clone(PlotMan->input().getFile(fileId).file->Get<TH2D>( "param_violin" ));
865  if(Violin[fileId] == nullptr)
866  {
867  MACH3LOG_ERROR("Couldn't find violin plot, make sure method from MCMCProcessor is being called");
868  return;
869  }
870  //KS: I know hardcoded but we can figure out later...
871  if(fileId == 0){
872  Violin[fileId]->SetFillColor(kBlue);
873  Violin[fileId]->SetFillColorAlpha(kBlue, 0.35);
874  Violin[fileId]->SetMarkerColor(kBlue);
875  Violin[fileId]->SetMarkerStyle(20);
876  Violin[fileId]->SetMarkerSize(0.5);
877  } else if (fileId == 1) {
878  Violin[fileId]->SetMarkerColor(kGreen);
879  Violin[fileId]->SetLineColor(kGreen);
880  Violin[fileId]->SetFillColor(kGreen);
881  Violin[fileId]->SetFillColorAlpha(kGreen, 0.35);
882  } else if (fileId == 2) {
883  Violin[fileId]->SetMarkerColor(kMagenta);
884  Violin[fileId]->SetLineColor(kMagenta);
885  Violin[fileId]->SetFillColor(kMagenta);
886  Violin[fileId]->SetFillColorAlpha(kMagenta, 0.35);
887  } else {
888  MACH3LOG_ERROR("Too many file, not implemented...");
889  throw MaCh3Exception(__FILE__ , __LINE__ );
890  }
891  }
892 
893  std::unique_ptr<TGraphAsymmErrors> PostGraphAll = MakeTGraphAsymmErrors(PlotMan->input().getFile(0).file);
894  // Make a Legend page
895  auto leg = std::make_unique<TLegend>(0.0, 0.0, 1.0, 1.0);
896  if (ViolinPre != nullptr) leg->AddEntry(ViolinPre.get(), "Prior", "lpf");
897  for(unsigned int fileId = 0; fileId < PlotMan->getNFiles(); fileId++) {
898  leg->AddEntry(Violin[fileId].get(), PlotMan->getFileLabel(fileId).c_str(), "lpf");
899  }
900  if(PlotAssym) leg->AddEntry(PostGraphAll.get(), "HPD Assym", "lp");
901  else leg->AddEntry(Postfit.get(), "HPD", "lpf");
902 
903  canvas->cd();
904  canvas->Clear();
905  leg->Draw();
906  canvas->Print((OutputName+".pdf").c_str());
907 
908  // get the names of the blocks of parameters to group together
909  std::vector<std::string> const blockNames = PlotMan->getOption<std::vector<std::string>>("paramGroups");
910  const int nPlots = static_cast<int>(blockNames.size());
911 
912  for (int i = 0; i < nPlots; i++)
913  {
914  // get the configuration for this parameter
915  std::string blockName = blockNames[i];
916  YAML::Node paramBlock = PlotMan->getOption(blockName);
917  std::string blockTitle = paramBlock[0].as<std::string>();
918  std::vector<double> blockLimits = paramBlock[1].as<std::vector<double>>();
919  std::vector<std::string> blockContents = paramBlock[2].as<std::vector<std::string>>();
920 
921  // get num of params in the block
922  const int nParams = static_cast<int>(blockContents.size());
923 
924  // set some plot things
925  auto blockHist_prefit = std::make_unique<TH2D>((blockName + "_Prefit").c_str(), blockTitle.c_str(), nParams, 0.0, static_cast<double>(nParams),
926  ViolinPre->GetYaxis()->GetNbins(), ViolinPre->GetYaxis()->GetXmin(), ViolinPre->GetYaxis()->GetXmax());
927  CopyViolinToBlock(ViolinPre.get(), blockHist_prefit.get(), blockContents);
928  // set the y axis limits we got from config
929  blockHist_prefit->GetYaxis()->SetRangeUser(blockLimits[0], blockLimits[1]);
930 
931  std::vector<std::unique_ptr<TH2D>> blockHist(PlotMan->getNFiles());
932  for(unsigned int fileId = 0; fileId < PlotMan->getNFiles(); fileId++) {
933  blockHist[fileId] = std::make_unique<TH2D>((blockTitle + "Violin" + fileId).Data(), (blockTitle + "Violin" + fileId).Data(),
934  nParams, 0.0, static_cast<double>(nParams), Violin[fileId]->GetYaxis()->GetNbins(),
935  Violin[fileId]->GetYaxis()->GetXmin(), Violin[fileId]->GetYaxis()->GetXmax());
936  CopyViolinToBlock(Violin[fileId].get(), blockHist[fileId].get(), blockContents);
937  blockHist[fileId]->GetYaxis()->SetRangeUser(blockLimits[fileId], blockLimits[1]);
938  }
939  // Do some fancy replacements
940  PrettifyTitles(blockHist_prefit.get());
941 
942  // set some plot things
943  auto blockHist_Best = std::make_unique<TH1D>(blockName.c_str(), blockTitle.c_str(),
944  nParams, 0.0, static_cast<double>(nParams));
945  // set the errors for the prefit block hist
946  for(int localBin=0; localBin < nParams; localBin ++){
947  // the "local" bin is the params index within the group of parameters
948  std::string paramName = blockContents[localBin];
949  copyParToBlockHist(localBin, paramName, blockHist_Best.get(), "", 0);
950  }
951 
952  std::vector<int> Index;
953  for(unsigned int is = 0; is < blockContents.size(); is++) {
954  int ParamBinId = M3::_BAD_INT_;
955  for (int ix = 0; ix < ViolinPre->GetXaxis()->GetNbins(); ++ix) {
956  if(ViolinPre->GetXaxis()->GetBinLabel(ix+1) == blockContents[is]) {
957  ParamBinId = ix;
958  break;
959  }
960  }
961  Index.push_back(ParamBinId);
962  }
963  std::unique_ptr<TGraphAsymmErrors> PostGraph = MakeTGraphAsymmErrors(PlotMan->input().getFile(0).file, Index);
964  PostGraph->SetMarkerColor(kBlack);
965  PostGraph->SetLineColor(kBlack);
966  PostGraph->SetMarkerStyle(7);
967  PostGraph->SetLineWidth(2);
968  PostGraph->SetLineStyle(kSolid);
969 
970  blockHist_prefit->Draw("violinX(03100300)");
971  for(unsigned int fileId = 0; fileId < PlotMan->getNFiles(); fileId++) {
972  blockHist[fileId]->Draw("violinX(03100300) SAME");
973  }
974 
975  if(PlotAssym) PostGraph->Draw("P SAME");
976  else Postfit->Draw("SAME");
977  canvas->Print((OutputName+".pdf").c_str());
978  }
979  canvas->Print((OutputName+".pdf]").c_str());
980 }
981 
983 void Get2DComparison(const std::string& FileName1, const std::string& FileName2)
984 {
985  auto canvas = std::make_unique<TCanvas>("canvas", "canvas", 0, 0, 1024, 1024);
986  canvas->SetBottomMargin(0.1f);
987  canvas->SetTopMargin(0.05f);
988  canvas->SetRightMargin(0.03f);
989  canvas->SetLeftMargin(0.15f);
990 
991  // Open the two ROOT files
992  TFile* File1 = M3::Open(FileName1, "READ", __FILE__, __LINE__);
993  TFile* File2 = M3::Open(FileName2, "READ", __FILE__, __LINE__);
994 
995  // Get the Post_2d_hists directory from both files
996  TDirectory* Dir1 = File1->Get<TDirectory>("Post_2d_hists");
997  TDirectory* Dir2 = File2->Get<TDirectory>("Post_2d_hists");
998 
999  if (!Dir1 || !Dir2) {
1000  MACH3LOG_WARN("Post_2d_hists directory not found in one or both files while running {}.", __func__);
1001  File1->Close();
1002  delete File1;
1003  File2->Close();
1004  delete File2;
1005  return;
1006  }
1007 
1008  // Get the list of keys in the first directory
1009  TIter next1(Dir1->GetListOfKeys());
1010  TKey* key1 = nullptr;
1011 
1012  // Prepare the output PDF filename
1013  std::string SaveName2D = "2DComparison_" + FileName1 + "_" + FileName2;
1014  SaveName2D = SaveName2D.substr(0, SaveName2D.find(".root"));
1015  SaveName2D = SaveName2D + ".pdf";
1016 
1017  canvas->Print((SaveName2D+"[").c_str());
1018  // Loop over keys in the first directory
1019  while ((key1 = static_cast<TKey*>(next1()))) {
1020  TString histName = key1->GetName();
1021 
1022  // Check if the key is a TH2D
1023  if (TString(key1->GetClassName()) == "TH2D") {
1024  TH2D* hist1 = static_cast<TH2D*>(key1->ReadObj());
1025 
1026  // Try to get the histogram with the same name from the second directory
1027  TH2D* hist2 = static_cast<TH2D*>(Dir2->Get(histName));
1028 
1029  if (hist2) {
1030  hist1->SetTitle("");
1031  hist1->SetTitle("");
1032 
1033  // Prettify axis titles
1034  std::string Xtitle = PlotMan->style().prettifyParamName(hist1->GetXaxis()->GetTitle());
1035  std::string Ytitle = PlotMan->style().prettifyParamName(hist1->GetYaxis()->GetTitle());
1036 
1037  // Adjust the axis ranges of hist1 to include both histograms
1038  double xmin = std::min(hist1->GetXaxis()->GetXmin(), hist2->GetXaxis()->GetXmin());
1039  double xmax = std::max(hist1->GetXaxis()->GetXmax(), hist2->GetXaxis()->GetXmax());
1040  double ymin = std::min(hist1->GetYaxis()->GetXmin(), hist2->GetYaxis()->GetXmin());
1041  double ymax = std::max(hist1->GetYaxis()->GetXmax(), hist2->GetYaxis()->GetXmax());
1042 
1043  hist1->GetXaxis()->SetRangeUser(xmin, xmax);
1044  hist1->GetYaxis()->SetRangeUser(ymin, ymax);
1045 
1046  hist1->GetXaxis()->SetTitle(Xtitle.c_str());
1047  hist1->GetYaxis()->SetTitle(Ytitle.c_str());
1048 
1049  hist1->SetLineColor(kBlue);
1050  hist1->SetLineStyle(kSolid);
1051  hist1->SetLineWidth(2);
1052 
1053  hist2->SetLineColor(kRed);
1054  hist2->SetLineStyle(kDashed);
1055  hist2->SetLineWidth(2);
1056 
1057  hist1->Draw("CONT3");
1058  hist2->Draw("CONT3 SAME");
1059 
1060  auto Legend = std::make_unique<TLegend>(0.20, 0.7, 0.4, 0.92);
1061  Legend->AddEntry(hist1, PlotMan->getFileLabel(0).c_str(), "l");
1062  Legend->AddEntry(hist2, PlotMan->getFileLabel(1).c_str(), "l");
1063  Legend->SetTextSize(0.03);
1064  Legend->SetLineColor(0);
1065  Legend->SetLineStyle(0);
1066  Legend->SetFillColor(0);
1067  Legend->SetFillStyle(0);
1068  Legend->SetBorderSize(0);
1069  Legend->Draw("SAME");
1070  canvas->Print((SaveName2D).c_str());
1071  }
1072  }
1073  }
1074  canvas->Print((SaveName2D+"]").c_str());
1075 
1076  File1->Close();
1077  delete File1;
1078  File2->Close();
1079  delete File2;
1080 }
1081 
1082 int main(int argc, char *argv[])
1083 {
1085  // Avoid Info in <TCanvas::Print>
1086  gErrorIgnoreLevel = kWarning;
1087 
1088  PlotMan = new MaCh3Plotting::PlottingManager();
1089  PlotMan->parseInputs(argc, argv);
1090  #ifdef DEBUG
1091  PlotMan->input().getFile(0).file->ls();
1092  #endif
1093  PlotMan->setExec("GetPostfitParamPlots");
1094 
1095  PlotMan->style().setPalette(PlotMan->getOption<std::string>("colorPalette"));
1096 
1098  GetViolinPlots();
1099 
1100  if (PlotMan->input().getNInputFiles() == 2)
1101  {
1102  std::string filename1 = PlotMan->getFileName(0);
1103  std::string filename2 = PlotMan->getFileName(1);
1104  Get2DComparison(filename1, filename2);
1105  }
1106 
1107  delete PlotMan;
1108  return 0;
1109 }
#define _MaCh3_Safe_Include_Start_
KS: Avoiding warning checking for headers.
Definition: Core.h:126
#define _MaCh3_Safe_Include_End_
KS: Restore warning checking after including external headers.
Definition: Core.h:140
void DrawPlots(TCanvas *plotCanv, TH1D *PrefitCopy, const std::vector< std::unique_ptr< TH1D >> &PostfitVec, TPad *mainPad, TPad *ratioPad)
int main(int argc, char *argv[])
std::vector< int > NDSamplesBins
TH1D * Prefit
void InitializePads(TCanvas *canvas, TPad *&pad3, TPad *&pad4)
void copyParToBlockHist(const int localBin, const std::string &paramName, TH1D *blockHist, const std::string &type, const int fileId, const bool setLabels=true)
MaCh3Plotting::PlottingManager * PlotMan
void MakeNDDetPlots()
std::unique_ptr< TH1D > makeRatio(TH1D *PrefitCopy, TH1D *PostfitCopy, bool setAxes)
bool ReadSettings(const std::shared_ptr< TFile > &File1)
constexpr Color_t PlotColor[]
std::string SaveName
void Get2DComparison(const std::string &FileName1, const std::string &FileName2)
KS: Make comparison of 2D Posteriors.
TPad * p3
void MakeRidgePlots()
void MakeFluxPlots()
void GetViolinPlots()
KS: Make fancy violin plots.
int NDParametersStartingPos
std::vector< std::string > NDSamplesNames
TPad * p4
void PrettifyTitles(TH1D *Hist)
std::unique_ptr< TGraphAsymmErrors > MakeTGraphAsymmErrors(const std::shared_ptr< TFile > &File, std::vector< int > Index={})
void GetPostfitParamPlots()
int NDParameters
void CopyViolinToBlock(TH2D *FullViolin, TH2D *ReducedViolin, const std::vector< std::string > &ParamNames)
TCanvas * canv
std::string plotType
void MakeParameterPlots()
#define MACH3LOG_CRITICAL
Definition: MaCh3Logger.h:38
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:34
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:37
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:35
void SetMaCh3LoggerFormat()
Set messaging format of the logger.
Definition: MaCh3Logger.h:61
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:36
#define MACH3LOG_TRACE
Definition: MaCh3Logger.h:33
Custom exception class used throughout MaCh3.
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.
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.
constexpr static const int _BAD_INT_
Default value used for int initialisation.
Definition: Core.h:55
void Print(const TTree *tree)
Definition: Monitor.cpp:325