GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Validation.cpp
Go to the documentation of this file.
1 #include <TCanvas.h>
2 #include <TF2.h>
3 #include <TGraph.h>
4 #include <TH1D.h>
5 
6 #include <cmath>
7 #include <cstring>
8 #include <fstream>
9 #include <iostream>
10 #include <sstream>
11 #include <string>
12 #include <vector>
13 
14 #include "gst.h"
15 
16 using namespace std;
17 
18 namespace {
19 
20 
21 const double q2_error = 0.05;
22 
23 void fill_delnu_plot_genie(const gst& gtree, const double q2, TH1D* delnu_histo) {
24  const double nucleon_mass = 0.939;
25  if (fabs(q2 - gtree.Q2) < q2_error) {
26  const double nu = gtree.Ev - gtree.El;
27  const double delnu = nu - gtree.Q2 / (2 * nucleon_mass);
28  delnu_histo->Fill(delnu);
29  }
30 }
31 
32 const char* get_nucleus_name(const int nucZ) {
33  switch(nucZ) {
34  case 1:
35  return "H";
36  case 2:
37  return "He";
38  case 6:
39  return "C";
40  case 10:
41  return "Ne";
42  case 13:
43  return "Al";
44  case 18:
45  return "Ar";
46  case 26:
47  return "Fe";
48  case 82:
49  return "Pb";
50  }
51  return "ERROR";
52 }
53 
54 TH1D* get_delnu_plot_genie(const int nucZ, const int nucA, const double q2) {
55  char histo_title[500];
56  sprintf(histo_title, "^{%d}%s, Q^{2} = (%0.1f #pm %0.2f) GeV",
57  nucA, get_nucleus_name(nucZ), q2, q2_error);
58  TH1D* delnu = new TH1D("genie", histo_title, 100, -1.0, 2.0);
59  delnu->GetXaxis()->SetTitle("#nu - Q^{2} / (2 * M_{p}) (GeV)");
60  delnu->GetYaxis()->SetTitle("(1 / #sigma) d#sigma / d#nu (GeV^{-1})");
61  delnu->SetLineColor(kRed);
62  const int start_run = 500000;
63  const int end_run = 500010;
64  for (int run = start_run; run <= end_run; run++) {
65  char data_file[500];
66  sprintf(data_file, "./data/EFF/nu_14/%d_%d/gntp.%d.gst.root", nucZ, nucA, run);
67  gst gtree(data_file);
68  const int num_events_per_file = 50000;
69  for (int event = 0; event < num_events_per_file; event++) {
70  gtree.GetEntry(event);
71  fill_delnu_plot_genie(gtree, q2, delnu);
72  }
73  }
74  delnu->Scale(1.0 / delnu->Integral("width"));
75  return delnu;
76 }
77 
78 bool fill_superscaling_vectors(const int nucZ, const int nucA, const double q2,
79  vector<float>* delnu_vals, vector<float>* relative_xs) {
80  char superscaling_data_file[500];
81  sprintf(superscaling_data_file, "superscaling/%d_%d_%0.1f.txt", nucZ, nucA, q2);
82  ifstream infile(superscaling_data_file);
83  if (!infile.is_open()) {
84  return false;
85  }
86  string line;
87  while (getline(infile, line)) {
88  float delnu, xs, garbage;
89  if (!(istringstream(line) >> delnu >> garbage >> xs)) {
90  cerr << "Couldn't process line: " << line << endl;
91  return false;
92  }
93  delnu_vals->push_back(delnu);
94  relative_xs->push_back(xs);
95  }
96  return true;
97 }
98 
99 TH1D* get_delnu_plot_superscaling(const int nucZ, const int nucA, const double q2) {
100  vector<float> delnu_vals;
101  vector<float> relative_xs;
102  if(!fill_superscaling_vectors(nucZ, nucA, q2, &delnu_vals, &relative_xs)) {
103  return new TH1D("", "", 100, -1, 1);
104  }
105  TH1D* delnu = new TH1D("SuSc", "SuSc", delnu_vals.size() - 1, &delnu_vals[0]);
106  for (int i = 0; i < (int)delnu_vals.size(); i++) {
107  delnu->Fill(delnu_vals[i], relative_xs[i]);
108  }
109  delnu->SetLineColor(kBlack);
110  delnu->Scale(1.0 / delnu->Integral("width"));
111  return delnu;
112 }
113 
114 void validation_plot_delnu(const int nucZ, const int nucA, const double q2) {
115  TCanvas canvas;
116  TH1D* genie_delnu = get_delnu_plot_genie(nucZ, nucA, q2);
117  TH1D* susc_delnu = get_delnu_plot_superscaling(nucZ, nucA, q2);
118  genie_delnu->Draw();
119  susc_delnu->Draw("SAME");
120  char outfile[500];
121  sprintf(outfile, "./plots/delnu_%d_%d_%0.1f.pdf", nucZ, nucA, q2);
122  for (int i = 1; i < (int)strlen(outfile) - 4; i++) {
123  outfile[i] = outfile[i] == '.' ? 'p' : outfile[i];
124  }
125  canvas.Print(outfile);
126  delete genie_delnu;
127  delete susc_delnu;
128 }
129 
130 void validation_plot_energy_xs_nu_nubar(bool is_bar) {
131  TCanvas canvas;
132  TFile def_file("./splines/DEF_QE_splines.root", "READ");
133  TFile eff_file("./splines/TE_QE_splines.root", "READ");
134  char hist_name[500];
135  sprintf(hist_name, "nu_mu_%sC12/qel_cc_%s",
136  (is_bar ? "bar_" : ""), (is_bar ? "p" : "n"));
137  TGraph* def_energy_xs = (TGraph*)def_file.Get(hist_name);
138  def_energy_xs->SetLineColor(kOrange);
139  TGraph* eff_energy_xs = (TGraph*)eff_file.Get(hist_name);
140  eff_energy_xs->SetLineColor(kRed);
141  def_energy_xs->Apply(&TF2("name,", "y / 6.0"));
142  eff_energy_xs->Apply(&TF2("name", "y / 6.0"));
143  char title[5000];
144  sprintf(title, "%sNeutrino Cross Section - %s + %s #rightarrow %s + #mu^%s",
145  (is_bar ? "Anti-" : ""), (is_bar ? "#bar{#nu}" : "#nu"),
146  (is_bar ? "p" : "n"), (is_bar ? "n" : "p"), (is_bar ? "{+}" : "{-}"));
147  TH1F* frame = canvas.DrawFrame(1e-1, 0, 1e2, (is_bar ? 1.2 : 1.6));
148  frame->SetTitle(title);
149  frame->GetYaxis()->SetTitle("#sigma (10^{-38} cm^{2})");
150  frame->GetXaxis()->SetTitle("E_{#nu} (GeV)");
151  frame->Draw();
152  eff_energy_xs->Draw("SAME");
153  def_energy_xs->Draw("SAME");
154  canvas.SetLogx();
155  char file_name[500];
156  sprintf(file_name, "./plots/nu_%sxs_energy.pdf", (is_bar ? "bar_" : ""));
157  canvas.Print(file_name);
158 }
159 
160 void validation_plot_energy_xs() {
161  validation_plot_energy_xs_nu_nubar(true);
162  validation_plot_energy_xs_nu_nubar(false);
163 }
164 
165 void validation_plot_energy_ratio_nu_nubar(bool is_bar) {
166  TCanvas canvas;
167  TFile def_file("./splines/DEF_QE_splines.root", "READ");
168  TFile eff_file("./splines/TE_QE_splines.root", "READ");
169  char hist_name[500];
170  sprintf(hist_name, "nu_mu_%sC12/qel_cc_%s",
171  (is_bar ? "bar_" : ""), (is_bar ? "p" : "n"));
172  TGraph* def_energy_xs = (TGraph*)def_file.Get(hist_name);
173  TGraph* eff_energy_xs = (TGraph*)eff_file.Get(hist_name);
174  TGraph ratio;
175  ratio.Set(def_energy_xs->GetN());
176  for (int i = 0; i < def_energy_xs->GetN(); i++) {
177  double x_val, def_point, eff_point;
178  def_energy_xs->GetPoint(i, x_val, def_point);
179  eff_energy_xs->GetPoint(i, x_val, eff_point);
180  if (def_point > 0) {
181  ratio.SetPoint(i, x_val, eff_point / def_point);
182  } else {
183  ratio.SetPoint(i, x_val, 0);
184  }
185  }
186  ratio.SetLineColor(kRed);
187  ratio.Draw();
188  char title[5000];
189  sprintf(title, "%s + %s #rightarrow %s + #mu^%s, (Trans Enh in G^{#nu}_{M}, M_{A} = 1.014) / (M_{A} = 1.014)",
190  (is_bar ? "#bar{#nu}" : "#nu"), (is_bar ? "p" : "n"), (is_bar ? "n" : "p"), (is_bar ? "{+}" : "{-}"));
191  TH1F* frame = canvas.DrawFrame(0, (is_bar ? 0.85 : 0.9), 10, (is_bar ? 1.4 : 1.5));
192  frame->SetTitle(title);
193  frame->GetYaxis()->SetTitle("ratio #sigma");
194  frame->GetXaxis()->SetTitle("E^{#nu} (GeV)");
195  frame->Draw();
196  ratio.Draw("SAME");
197  char file_name[500];
198  sprintf(file_name, "./plots/nu_%sratio_energy.pdf", (is_bar ? "bar_" : ""));
199  canvas.Print(file_name);
200 }
201 
202 void validation_plot_energy_ratio() {
203  validation_plot_energy_ratio_nu_nubar(true);
204  validation_plot_energy_ratio_nu_nubar(false);
205 }
206 
207 double get_qe_increase_factor(double is_bar, double energy = 3.0) {
208  TFile def_file("./splines/DEF_QE_splines.root", "READ");
209  TFile eff_file("./splines/TE_QE_splines.root", "READ");
210  char hist_name[500];
211  sprintf(hist_name, "nu_mu_%sC12/qel_cc_%s",
212  (is_bar ? "bar_" : ""), (is_bar ? "p" : "n"));
213  TGraph* def_energy_xs = (TGraph*)def_file.Get(hist_name);
214  TGraph* eff_energy_xs = (TGraph*)eff_file.Get(hist_name);
215  for (int i = 0; i < def_energy_xs->GetN(); i++) {
216  double point_energy, def_xs, eff_xs;
217  def_energy_xs->GetPoint(i, point_energy, def_xs);
218  if (point_energy > energy) {
219  eff_energy_xs->GetPoint(i, point_energy, eff_xs);
220  return eff_xs / def_xs;
221  }
222  }
223  return -1;
224 }
225 
226 void populate_q2_plot(bool is_bar, const char* tag, TH1D* q2) {
227  const int start_run = 500000;
228  const int end_run = 500010;
229  for (int run = start_run; run <= end_run; run++) {
230  char data_file[500];
231  sprintf(data_file, "./data/%s/nu_%d/6_12/gntp.%d.gst.root", tag,
232  (is_bar ? -14 : 14), run);
233  gst gtree(data_file);
234  const int num_events_per_file = 50000;
235  for (int event = 0; event < num_events_per_file; event++) {
236  gtree.GetEntry(event);
237  q2->Fill(gtree.Q2);
238  }
239  }
240 }
241 
242 
243 void validation_plot_q2_ratio_nu_nubar(bool is_bar) {
244  TCanvas canvas;
245  const double increase_factor = get_qe_increase_factor(is_bar);
246  char title[500];
247  sprintf(title, "E_{#nu} = 3 GeV, %s + %s #rightarrow %s + #mu^%s, (Trans Enh in G^{#nu}_{M}, M_{A} = 1.014) / (M_{A} = 1.014)",
248  (is_bar ? "#bar{#nu}" : "#nu"), (is_bar ? "p" : "n"), (is_bar ? "n" : "p"), (is_bar ? "{+}" : "{-}"));
249  TH1D te_q2("", title, 50, 0, 2.5);
250  te_q2.GetYaxis()->SetTitle("ratio d#sigma / dQ^{2}");
251  te_q2.GetXaxis()->SetTitle("Q^{2} (GeV / c)^{2}");
252  TH1D def_q2("", "",50, 0, 2.5);
253  populate_q2_plot(is_bar, "TE", &te_q2);
254  populate_q2_plot(is_bar, "DEF", &def_q2);
255  te_q2.Scale(increase_factor);
256  te_q2.Divide(&def_q2);
257  te_q2.SetLineColor(kRed);
258  te_q2.Draw();
259  char filename[500];
260  sprintf(filename, "./plots/nu%s_ratio_q2.pdf", (is_bar ? "_bar" : ""));
261  canvas.Print(filename);
262 }
263 
264 
265 void validation_plot_q2_ratio() {
266  validation_plot_q2_ratio_nu_nubar(true);
267  validation_plot_q2_ratio_nu_nubar(false);
268 }
269 
270 } // namespace
271 
272 void Validation() {
273  const int nucZs[] = {1, 2, 6, 10, 13, 18, 26, 82 };
274  const int nucAs[] = {2, 4, 12, 20, 27, 40, 56, 208};
275  const int num_nucs = 8;
276 
277  const double q2s[] = {0.1, 0.3, 0.5, 0.7, 1.0, 1.2, 1.5, 2.0};
278  const int num_q2s = 8;
279  validation_plot_energy_xs();
280  validation_plot_energy_ratio();
281  validation_plot_q2_ratio();
282 
283  for (int i = 0; i < num_nucs; i++) {
284  for (int j = 0; j < num_q2s; j++) {
285  validation_plot_delnu(nucZs[i], nucAs[i], q2s[j]);
286  }
287  }
288 }
289 
290 int main() {
291  Validation();
292  return 0;
293 }
void Validation()
Definition: Validation.cpp:272
Double_t Ev
Definition: gst.h:47
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
const double e
string infile
Double_t El
Definition: gst.h:55
Definition: gst.h:8
Double_t Q2
Definition: gst.h:45
virtual Int_t GetEntry(Long64_t entry)
Definition: gst.h:206