nuint09_qel1.C
Go to the documentation of this file.
1 //
2 // NuINT09 Conference, Benchmark Calculations (GENIE contribution)
3 //
4 // QEL.1:
5 // QEL total cross section as a function of energy
6 //
7 // Inputs:
8 // - sample id: 0 (nu_mu+C12), 1 (nu_mu+O16), 2 (nu_mu+Fe56)
9 //
10 // Costas Andreopoulos, STFC / Rutherford Appleton Laboratory
11 //
12 #include <iomanip>
13 
14 const int kNSamples = 3;
15 
16 const char * kLabel[kNSamples] =
17 {
18  // available samples
19  /* 0 */ "nu_mu_C12",
20  /* 1 */ "nu_mu_O16",
21  /* 2 */ "nu_mu_Fe56"
22 };
23 int kZ[kNSamples] =
24 {
25  // Z for nuclear target at each sample
26  /* 0 */ 6,
27  /* 1 */ 8,
28  /* 2 */ 26
29 };
30 int kA[kNSamples] =
31 {
32  // A for nuclear target at each sample
33  /* 0 */ 12,
34  /* 1 */ 16,
35  /* 2 */ 56
36 };
37 
38 void nuint09_qel1(int isample)
39 {
40  cout << " ***** running: QEL.1" << endl;
41 
42  if(isample<0 || isample >= kNSamples) return;
43 
44  const char * label = kLabel[isample];
45 
46  int A = kA[isample];
47  int Z = kZ[isample];
48  int N = A-Z;
49 
50  // get cross section graphs
51 
52  TFile fsig("../sig/splines.root","read");
53  TDirectory * sig_dir = (TDirectory *) fsig.Get(label);
54 
55  TGraph * sig_graph_qelcc = (TGraph*) sig_dir->Get("qel_cc_n");
56 
57  // range & spacing
58 
59  const int ne = 60;
60  const double emin = 0.05;
61  const double emax = 30.00;
62 
63  const double dloge = (TMath::Log10(emax) - TMath::Log10(emin)) / (ne-1);
64 
65  // create output stream
66 
67  ostringstream out_filename;
68  out_filename << label << ".qel_1.sig_vs_Enu.data";
69  ofstream out_stream(out_filename.str().c_str(), ios::out);
70 
71  // write out txt file
72 
73  out_stream << "# [" << label << "]" << endl;
74  out_stream << "# " << endl;
75  out_stream << "# [QEL.1]:" << endl;
76  out_stream << "# QEL total cross section as a function of neutrino energy" << endl;
77  out_stream << "# " << endl;
78  out_stream << "# Note:" << endl;
79  out_stream << "# - neutrino energy E in GeV, log spacing between Emin = " << emin << " GeV, Emax = " << emax << " GeV " << endl;
80  out_stream << "# - cross sections in 1E-38 cm^2 " << endl;
81  out_stream << "# - quoted cross section is nuclear cross section per nucleon contributing in the scattering (eg only neutrons for nu_mu QELCC)" << endl;
82  out_stream << "# Columns:" << endl;
83  out_stream << "# | Energy | sig(nu_mu + n[bound] -> mu- + p) | " << endl;
84 
85  out_stream << std::fixed << setprecision(6);
86 
87  for(int i=0; i < ne; i++) {
88 
89  double e = TMath::Power(10., TMath::Log10(emin) + i * dloge);
90 
91  double sig_qelcc = sig_graph_qelcc->Eval(e) / N;
92  sig_qelcc = TMath::Max(0., sig_qelcc);
93 
94  out_stream << setw(15) << e << setw(15) << sig_qelcc << endl;
95  }
96 
97  out_stream.close();
98 
99  // visual inspection
100  TCanvas * c1 = new TCanvas("c1","",20,20,500,500);
101  sig_graph_qelcc->Draw("alp");
102  c1->Update();
103 }
const char * kLabel[kNSamples]
Definition: nuint09_qel1.C:16
const char * label
const double emin
Float_t Z
Definition: plot.C:38
int kZ[kNSamples]
Definition: nuint09_qel1.C:23
const double emax
const int kNSamples
Definition: nuint09_qel1.C:14
OStream cout
Definition: OStream.cxx:6
void nuint09_qel1(int isample)
Definition: nuint09_qel1.C:38
static const double A
Definition: Units.h:82
int kA[kNSamples]
Definition: nuint09_qel1.C:30
c1
Definition: demo5.py:24
Float_t e
Definition: plot.C:35