5#include "TGraphAsymmErrors.h"
22 double MinVale = 1234567890.;
23 TAxis *xaxis = hist->GetXaxis();
24 for (
int x = 1; x <= hist->GetNbinsX(); x++)
26 if (xaxis->GetBinLowEdge(x) > minRange && xaxis->GetBinUpEdge(x) < maxRange)
28 if (MinVale > hist->GetBinContent(x))
29 MinVale = hist->GetBinContent(x);
39 int failure_count = 0;
41 for (
int bin = 2; bin <= hist->GetNbinsX(); ++bin)
43 if (fabs(hist->GetBinContent(bin) - 1.0) > tolerance)
45 if (++failure_count > max_failures)
54void MakePlot(TString fname1, TString flabel1, TString fname2, TString flabel2, TString fname3, TString flabel3, TString fname4, TString flabel4)
56 TCanvas *c1 =
new TCanvas(
"c1",
" ", 0, 0, 800, 630);
57 gStyle->SetOptStat(0);
59 gErrorIgnoreLevel = kWarning;
62 TFile *infile = TFile::Open(fname1.Data());
64 bool add_legend =
false;
66 TFile *infile2 = NULL;
69 infile2 = TFile::Open(fname2.Data());
72 TFile *infile3 = NULL;
75 infile3 = TFile::Open(fname3.Data());
78 TFile *infile4 = NULL;
81 infile4 = TFile::Open(fname4.Data());
85 TIter next(infile->GetListOfKeys());
86 while ((key =
static_cast<TKey *
>(next())))
88 std::string dirname = std::string(key->GetName());
89 if (std::string(key->GetClassName()) !=
"TDirectoryFile")
92 if ((dirname ==
"LogL") || (dirname ==
"Batched_means"))
97 infile->cd(dirname.c_str());
98 TIter nextsub(gDirectory->GetListOfKeys());
99 c1->Print(Form(
"%s.pdf[", dirname.c_str()),
"pdf");
101 while ((subkey =
static_cast<TKey *
>(nextsub())))
104 TLegend *leg =
new TLegend(0.7, 0.7, 0.9, 0.9);
105 leg->SetFillColorAlpha(kWhite,
float(0.7));
106 std::string name = std::string(subkey->GetName());
107 name = dirname +
"/" + name;
109 if (std::string(subkey->GetClassName()) !=
"TH1D")
119 MACH3LOG_INFO(
"Looking for {} from file {}", name.c_str(), fname1.Data());
120 blarb[0] =
static_cast<TH1D *
>(infile->Get(name.c_str())->Clone());
122 if (TMath::IsNaN(blarb[0]->GetBinContent(1)))
126 blarb[0]->SetLineStyle(kSolid);
127 blarb[0]->SetLineColor(kRed);
130 leg->AddEntry(blarb[0], flabel1.Data(),
"l");
132 if (dirname ==
"AccProb")
133 blarb[0]->GetYaxis()->SetRangeUser(0, 1.0);
134 if (name ==
"AccProb/AcceptanceProbability")
138 blarb[1] =
static_cast<TH1D *
>(infile2->Get(name.c_str())->Clone());
140 blarb[1]->SetLineStyle(kDashed);
141 blarb[1]->SetLineColor(kBlue);
142 blarb[1]->Draw(
"same");
143 leg->AddEntry(blarb[1], flabel2.Data(),
"l");
147 blarb[2] =
static_cast<TH1D *
>(infile3->Get(name.c_str())->Clone());
149 blarb[2]->SetLineStyle(kDotted);
150 blarb[2]->SetLineColor(kGreen);
151 blarb[2]->Draw(
"same");
152 leg->AddEntry(blarb[2], flabel3.Data(),
"l");
156 blarb[3] =
static_cast<TH1D *
>(infile4->Get(name.c_str())->Clone());
158 blarb[3]->SetLineStyle(kDashDotted);
159 blarb[3]->SetLineColor(kOrange);
160 blarb[3]->Draw(
"same");
161 leg->AddEntry(blarb[3], flabel4.Data(),
"l");
168 c1->Print(Form(
"%s.pdf", dirname.c_str()),
"pdf");
171 gDirectory->cd(
"..");
172 c1->Print(Form(
"%s.pdf]", dirname.c_str()),
"pdf");
184void PlotAutoCorr(TString fname1, TString flabel1, TString fname2, TString flabel2, TString fname3, TString flabel3, TString fname4, TString flabel4)
192 std::vector<TString> flabel = {flabel1, flabel2, flabel3, flabel4};
195 infile[0] = TFile::Open(fname[0].Data());
201 infile[1] = TFile::Open(fname[1].Data());
206 infile[2] = TFile::Open(fname[2].Data());
211 infile[3] = TFile::Open(fname[3].Data());
215 TCanvas *c1 =
new TCanvas(
"c1",
" ", 0, 0, 800, 630);
216 gStyle->SetOptStat(0);
218 gErrorIgnoreLevel = kWarning;
220 c1->Print(
"Auto_Corr_PerFile.pdf[",
"pdf");
221 for (
int ik = 0; ik < Nfiles; ik++)
223 TIter next(infile[ik]->GetListOfKeys());
226 TLegend *leg =
new TLegend(0.7, 0.7, 0.9, 0.9);
228 while ((key =
static_cast<TKey *
>(next())))
230 std::string dirname = std::string(key->GetName());
233 if ((dirname !=
"Auto_corr"))
236 infile[ik]->cd(dirname.c_str());
237 TIter nextsub(gDirectory->GetListOfKeys());
240 bool FirstTime =
true;
242 while ((subkey =
static_cast<TKey *
>(nextsub())))
244 std::string name = std::string(subkey->GetName());
245 name = dirname +
"/" + name;
247 if (std::string(subkey->GetClassName()) !=
"TH1D")
250 TH1D *blarb =
static_cast<TH1D *
>(infile[ik]->Get(name.c_str())->Clone());
252 if (TMath::IsNaN(blarb->GetBinContent(1)))
256 delete blarb->GetListOfFunctions()->FindObject(
"Fitter");
260 if (MinValue >= 0.80)
261 blarb->SetLineColor(kRed);
262 else if (MinValue >= 0.40 && MinValue < 0.80)
263 blarb->SetLineColor(kOrange);
264 else if (MinValue > 0.20 && MinValue < 0.40)
265 blarb->SetLineColor(kYellow);
266 else if (MinValue <= 0.20)
267 blarb->SetLineColor(kGreen);
268 blarb->GetXaxis()->UnZoom();
271 blarb->SetTitle(Form(
"Auto_Corr_%s.pdf", fname[ik].Data()));
273 blarb->SetLineStyle(kDashed);
278 TString integral = Form(
"Integral: %.2f", blarb->Integral());
284 gDirectory->cd(
"..");
291 c1->Print(
"Auto_Corr_PerFile.pdf",
"pdf");
294 c1->Print(
"Auto_Corr_PerFile.pdf]",
"pdf");
302 int nBins = hist->GetNbinsX();
305 for (
int i = 0; i <
nBins; ++i)
307 x[i] = hist->GetBinCenter(i + 1);
308 ymin[i] = hist->GetBinContent(i + 1);
309 ymax[i] = hist->GetBinContent(i + 1);
312 TGraph *minGraph =
new TGraph(
nBins, x.data(), ymin.data());
313 TGraph *maxGraph =
new TGraph(
nBins, x.data(), ymax.data());
315 minGraph->SetLineColor(color);
316 maxGraph->SetLineColor(color);
317 minGraph->SetFillColorAlpha(color,
float(0.3));
318 maxGraph->SetFillColorAlpha(color,
float(0.3));
320 return {minGraph, maxGraph};
325 if (histograms.empty())
327 throw std::invalid_argument(
"Empty histogram vector provided");
330 int nBins = histograms[0]->GetNbinsX();
333 for (
int bin = 0; bin <
nBins; bin++)
335 x[bin] = histograms[0]->GetBinCenter(bin + 1);
337 for (
const auto &hist : histograms)
339 double content = hist->GetBinContent(bin + 1);
340 ymin[bin] = std::min(ymin[bin], content);
341 ymax[bin] = std::max(ymax[bin], content);
346 TGraphAsymmErrors *band =
new TGraphAsymmErrors(
nBins);
347 for (
int i = 0; i <
nBins; ++i)
349 band->SetPoint(i, x[i], (ymax[i] + ymin[i]) / 2.0);
350 band->SetPointError(i, 0, 0, (ymax[i] - ymin[i]) / 2.0, (ymax[i] - ymin[i]) / 2.0);
353 band->SetFillColorAlpha(color,
float(0.2));
354 band->SetLineWidth(1);
361 if (histograms.empty())
363 throw std::invalid_argument(
"Empty histogram vector provided");
366 TH1D *min_hist =
static_cast<TH1D *
>(histograms[0]->Clone());
367 TH1D *max_hist =
static_cast<TH1D *
>(histograms[0]->Clone());
369 for (
const auto &hist : histograms)
371 for (
int bin = 1; bin <= min_hist->GetNbinsX(); ++bin)
373 double current_min = min_hist->GetBinContent(bin);
374 double current_max = max_hist->GetBinContent(bin);
375 double bin_content = hist->GetBinContent(bin);
377 min_hist->SetBinContent(bin, std::min(current_min, bin_content));
378 max_hist->SetBinContent(bin, std::max(current_max, bin_content));
382 return {min_hist, max_hist};
388 int ¶meter_count,
389 std::vector<TH1D *> &histograms)
391 TIter next(autocor_dir->GetListOfKeys());
394 while ((key =
dynamic_cast<TKey *
>(next())))
396 TH1D *current_hist =
nullptr;
397 autocor_dir->GetObject(key->GetName(), current_hist);
400 current_hist->GetMaximum() <= 0 ||
406 current_hist->SetDirectory(
nullptr);
407 histograms.push_back(current_hist);
411 average_hist =
static_cast<TH1D *
>(current_hist->Clone());
412 average_hist->SetDirectory(
nullptr);
416 average_hist->Add(current_hist);
424 int ¶meter_count,
425 std::vector<TH1D *> &histograms)
427 std::unique_ptr<TFile> input_file(TFile::Open(file_path));
428 if (!input_file || input_file->IsZombie())
430 throw std::runtime_error(
"Could not open file: " + std::string(file_path.Data()));
433 TDirectoryFile *autocor_dir =
nullptr;
434 input_file->GetObject(
"Auto_corr", autocor_dir);
439 "Auto_corr directory not found in file: " + std::string(file_path.Data()));
448 TH1D *average_hist =
nullptr;
449 int parameter_count = 0;
455 catch (
const std::exception &e)
457 MACH3LOG_ERROR(
"Error processing file {} : {}", input_file, e.what());
460 if (average_hist && parameter_count > 0)
462 MACH3LOG_INFO(
"Processed {} parameters from {} files", parameter_count, input_file);
463 average_hist->Scale(1.0 / parameter_count);
470 const std::vector<TH1D *> &averages,
471 const std::vector<TString> &hist_labels,
472 const TString &output_name,
473 bool draw_min_max =
true,
474 bool draw_all =
false,
475 bool draw_errors =
true)
477 TCanvas *canvas =
new TCanvas(
"AverageAC",
"Average Auto Correlation", 800, 600);
482 Double_t stops[8] = {0.0, 0.25, 0.5, 0.75};
483 Double_t red[8] = {0.83203125, 0.796875, 0.0, 0.9375};
484 Double_t green[8] = {0.3671875, 0.47265625, 0.4453125, 0.890625};
485 Double_t blue[8] = {0.0, 0.65234375, 0.6953125, 0.2578125};
486 TColor::CreateGradientColorTable(8, stops, red, green, blue, nb);
488 TLegend *leg =
new TLegend(0.5, 0.7, 0.9, 0.9);
489 leg->SetFillColorAlpha(kWhite,
float(0.7));
493 for (
size_t i = 0; i < histograms.size(); ++i)
495 auto colour =
static_cast<Color_t
>(TColor::GetColorPalette(
static_cast<Int_t
>(i * nb / histograms.size())));
499 band->SetLineWidth(0);
500 band->SetTitle(
"Average Auto Correlation");
501 band->GetXaxis()->SetTitle(
"lag");
502 band->GetYaxis()->SetTitle(
"Autocorrelation Function");
509 band->Draw(
"3 SAME");
515 for (
size_t i = 0; i < averages.size(); ++i)
517 auto colour =
static_cast<Color_t
>(TColor::GetColorPalette(
static_cast<Int_t
>(i * nb / histograms.size())));
521 auto error_hist =
static_cast<TH1D *
>(averages[i]->Clone());
522 error_hist->SetFillColorAlpha(colour,
float(0.3));
523 error_hist->SetLineWidth(0);
524 error_hist->SetTitle(
"Average Auto Correlation");
525 error_hist->GetXaxis()->SetTitle(
"lag");
526 error_hist->GetYaxis()->SetTitle(
"Autocorrelation Function");
527 error_hist->Draw(
"E3 SAME");
531 for (
const auto &hist : histograms[i])
533 hist->SetLineColorAlpha(colour,
float(0.05));
534 hist->SetLineWidth(1);
535 hist->Draw(
"HIST SAME");
539 averages[i]->SetLineColor(colour);
540 averages[i]->SetLineWidth(2);
541 averages[i]->SetTitle(
"Average Auto Correlation");
542 averages[i]->GetXaxis()->SetTitle(
"lag");
543 averages[i]->GetYaxis()->SetTitle(
"Autocorrelation Function");
544 averages[i]->Draw(
"HIST SAME");
546 TString int_str = TString::Format(
"Average integrated autocorrelation %f", averages[i]->Integral());
548 leg->AddEntry(averages[i], hist_labels[i] + int_str,
"l");
551 if (averages.size() > 1)
556 canvas->SaveAs(output_name +
".pdf");
561 std::vector<TString> hist_labels,
562 const TString &output_name,
563 bool draw_min_max =
true)
568 std::vector<std::vector<TH1D *>> histograms;
569 std::vector<TH1D *> averages;
571 for (
int i = 0; i < static_cast<int>(input_files.size()); i++)
573 TString folder = input_files[i];
576 if (folder.IsNull() || folder ==
DUMMYFILE)
578 MACH3LOG_WARN(
"Skipping empty or dummy folder: {}", folder.Data());
582 std::vector<TH1D *> histograms_i;
584 histograms.push_back(histograms_i);
587 CompareAverageAC(histograms, averages, hist_labels, output_name, draw_min_max);
590int main(
int argc,
char *argv[])
593 if (argc != 3 && argc != 5 && argc != 7 && argc != 9)
596 MACH3LOG_ERROR(
"How to use: {} DiagMCMC_Output.root Plot Name", argv[0]);
611 PlotAverageACMult({argv[1], argv[3]}, {argv[2], argv[4]},
"Average_Auto_Corr",
true);
617 PlotAverageACMult({argv[1], argv[3], argv[5]}, {argv[2], argv[4], argv[6]},
"Average_Auto_Corr",
true);
621 MakePlot(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7], argv[8]);
622 PlotAutoCorr(argv[1], argv[2], argv[3], argv[4], argv[5], argv[6], argv[7], argv[8]);
623 PlotAverageACMult({argv[1], argv[3], argv[5], argv[7]}, {argv[2], argv[4], argv[6], argv[8]},
"Average_Auto_Corr",
true);
void RemoveFitter(TH1D *hist, const std::string &name)
KS: Remove fitted TF1 from hist to make comparison easier.
void SetMaCh3LoggerFormat()
Set messaging format of the logger.
int main(int argc, char *argv[])
bool IsHistogramAllOnes(TH1D *hist, double tolerance=0.001, int max_failures=100)
HW: Check if histogram is flat within a given tolerance.
void ProcessAutoCorrelationDirectory(TDirectoryFile *autocor_dir, TH1D *&average_hist, int ¶meter_count, std::vector< TH1D * > &histograms)
double GetMinimumInRange(TH1D *hist, double minRange, double maxRange)
KS: function which looks for minimum in given range.
void MakePlot(TString fname1, TString flabel1, TString fname2, TString flabel2, TString fname3, TString flabel3, TString fname4, TString flabel4)
std::pair< TH1D *, TH1D * > CalculateMinMaxHistograms(const std::vector< TH1D * > &histograms)
void PlotAverageACMult(std::vector< TString > input_files, std::vector< TString > hist_labels, const TString &output_name, bool draw_min_max=true)
void CompareAverageAC(const std::vector< std::vector< TH1D * > > &histograms, const std::vector< TH1D * > &averages, const std::vector< TString > &hist_labels, const TString &output_name, bool draw_min_max=true, bool draw_all=false, bool draw_errors=true)
std::pair< TGraph *, TGraph * > CreateMinMaxBand(TH1D *hist, Color_t color)
HW: Create a band of minimum and maximum values from a histogram.
TGraph * CalculateMinMaxBand(const std::vector< TH1D * > &histograms, Color_t color)
void PlotAutoCorr(TString fname1, TString flabel1, TString fname2, TString flabel2, TString fname3, TString flabel3, TString fname4, TString flabel4)
TH1D * AutocorrProcessInputs(const TString &input_file, std::vector< TH1D * > &histograms)
void ProcessDiagnosticFile(const TString &file_path, TH1D *&average_hist, int ¶meter_count, std::vector< TH1D * > &histograms)
Custom exception class for MaCh3 errors.