tutorials/nominal_plotter.C at main · NUISANCEMC/tutorials · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
double get_true_max(TH1* hist){
  double maxVal = 0;
  for (int x = 0; x < hist->GetNbinsX(); ++x){
    double thisVal = hist->GetBinContent(x+1)+hist->GetBinError(x+1);
    if (thisVal > maxVal) maxVal = thisVal;
  }
  return maxVal;
}

void make_basic_plot(std::string input_file_name, std::string measurement_name, std::string output_plot_name){

  gROOT->SetBatch(1);

  // Options to make the histogram prettier
  gStyle->SetLineWidth(3);
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gStyle->SetOptFit(0);
  TGaxis::SetMaxDigits(3);

  gStyle->SetTextSize(0.06);
  gStyle->SetLabelSize(0.05,"xyzt");
  gStyle->SetTitleSize(0.06,"xyzt");

  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);
  gStyle->SetNdivisions(505, "XY");

  // ROOT sometimes tries to ignore what you tell it...
  gROOT ->ForceStyle();

  // Make a canvas (where we will draw the histogram)
  TCanvas *can = new TCanvas("can", "can", 800, 800);

  // Open the file
  TFile *in_file = new TFile(input_file_name.c_str());

  // Get the generator prediction histogram
  TH1D *mc_hist = (TH1D*)in_file->Get((measurement_name+"_MC").c_str());

  // Get the data histogram
  TH1D *data_hist = (TH1D*)in_file->Get((measurement_name+"_data").c_str());

  // This is a bit of a NUISANCE hack, the chisq between data and MC is stored as the MC histogram title
  std:string chisq = mc_hist->GetTitle();

  // Presentation is important
  mc_hist  ->SetLineWidth(3);
  mc_hist  ->SetLineColor(kRed);
  data_hist->SetLineWidth(3);
  data_hist->SetLineColor(kBlack);

  // Make a legend
  TLegend *leg = new TLegend(0.35, 0.82, 1.00, 0.99, "", "NDC");

  // By default ROOT legends can be ugly... make it nicer
  leg ->SetShadowColor(0);
  leg ->SetFillColor(0);
  leg ->SetLineWidth(0);
  leg ->SetTextSize(0.04);
  leg ->SetLineColor(kWhite);
  leg ->SetNColumns(2);

  // The information to be conveyed
  leg ->AddEntry(data_hist, "Data", "l");
  leg ->AddEntry(mc_hist, ("MC (#chi^{2} = "+chisq+")").c_str(), "l");

  // Sort out the maximum value
  double maxVal = get_true_max(data_hist);
  if (mc_hist->GetMaximum() > maxVal){
    maxVal = mc_hist->GetMaximum();
  }

  // Note that all style options should be applied to the first histogram we draw
  // Add a small buffer to make it look better
  data_hist ->SetMaximum(maxVal*1.1);

  // Now draw everything on the canvas
  can ->cd();
  data_hist ->Draw();
  data_hist ->GetYaxis()->SetTitleOffset(0.9);
  mc_hist ->Draw("HIST SAME");
  leg ->Draw("SAME");

  // Make the borders sensible
  gPad ->SetRightMargin(0.03);
  gPad ->SetTopMargin(0.2);
  gPad ->SetBottomMargin(0.15);
  gPad ->SetLeftMargin(0.15);
  can ->Update();

  // This actually makes the plot
  can ->SaveAs(output_plot_name.c_str());
  return;
};


// ROOT CINT expects a function with the same name as the script
// A feature is that subsequent command line arguments are passed to this function
void nominal_plotter(){

  std::string input_file_name  = "nominal_GENIEv3_example.root";
  std::string samples[] = {"MINERvA_CC0pinp_STV_XSec_1Dpmu_nu",
    "MINERvA_CC1pip_XSec_1DTpi_nu_2017",
    "MINERvA_CC1pi0_XSec_1DTpi_nu"};

  for (const auto& sample : samples) {
    std::string output_plot_name = "plots/nominal_"+sample+"_GENIEv3.png";
    make_basic_plot(input_file_name, sample, output_plot_name);
  }

  return;
}