nuint09_1pi1.C
Go to the documentation of this file.
1 //
2 // NuINT09 Conference, Benchmark Calculations (GENIE contribution)
3 //
4 // 1PI.1:
5 // Total cross section for 1 pion (pi+ only) production as a function of neutrino energy
6 //
7 // Inputs:
8 // - sample id: 0 (nu_mu+C12), 1 (nu_mu+O16), 2 (nu_mu+Fe56)
9 // - single pion source: 0 (all), 1 (P33(1232) resonance only), 2 (resonances only)
10 // - stage: 0 -> primary Xpi+, 1 -> final state Xpi+
11 //
12 // Costas Andreopoulos, STFC / Rutherford Appleton Laboratory
13 //
14 #include <iomanip>
15 
16 //
17 // consts
18 //
19 const int kNSamples = 3;
20 const int kNWCur = 1;
21 const int kNRunsPerCase = 10;
22 const int kNEvtPerRun = 100000;
23 
24 const char * kLabel[kNSamples] =
25 {
26  // available samples
27  /* 0 */ "nu_mu_C12",
28  /* 1 */ "nu_mu_O16",
29  /* 2 */ "nu_mu_Fe56"
30 };
32 {
33  /* indices : sample ; cc/nc */
34  {
35  /* 0,0 (nu_mu C12, CC) */ {909900, 909901, 909902, 909903, 909904, 909905, 909906, 909907, 909908, 909909}
36  },
37  {
38  /* 1,0 (nu_mu O16, CC) */ {919900, 919901, 919902, 919903, 919904, 919905, 919906, 919907, 919908, 919909}
39  },
40  {
41  /* 2,0 (nu_mu Fe56, CC) */ {929900, 929901, 929902, 929903, 929904, 929905, 929906, 929907, 929908, 929909}
42  }
43 };
44 int kA[kNSamples] =
45 {
46  // A for nuclear target at each sample
47  /* 0 */ 12,
48  /* 1 */ 16,
49  /* 2 */ 56
50 };
51 
52 void nuint09_1pi1(int isample, int single_pion_sources=0, int stage=1)
53 {
54  cout << " ***** running: 1PI.1" << endl;
55 
56  if(isample<0 || isample >= kNSamples) return;
57  if(single_pion_sources<0 || single_pion_sources>2) return;
58 
59  const char * label = kLabel[isample];
60 
61  int A = kA[isample];
62 
63  // get cross section graphs
64 
65  TFile fsig("../sig/splines.root","read");
66  TDirectory * sig_dir = (TDirectory *) fsig.Get(label);
67 
68  TGraph * sig_graph_totcc = (TGraph*) sig_dir->Get("tot_cc");
69 
70  // range & spacing
71 
72  const int nEnu = 60;
73  const double Enumin = 0.05;
74  const double Enumax = 30.00;
75 
76  // create output stream
77 
78  ostringstream out_filename;
79  out_filename << label;
80 
81  if (single_pion_sources==0) out_filename << ".1pi_1a.";
82  else if (single_pion_sources==1) out_filename << ".1pi_1b.";
83  else if (single_pion_sources==2) out_filename << ".1pi_1c.";
84 
85  if(stage==0) out_filename << "no_FSI.";
86 
87  out_filename << "sig1pi_vs_Enu.data";
88  ofstream out_stream(out_filename.str().c_str(), ios::out);
89 
90  // write out txt file
91 
92  out_stream << "# [" << label << "]" << endl;
93  out_stream << "# " << endl;
94  out_stream << "# [1PI.1]:" << endl;
95  out_stream << "# Total cross section for 1 pion (pi+ only) production as a function of energy" << endl;
96  if(stage==0) {
97  out_stream << "# ***** NO FSI: The {X pi+} state is a primary hadronic state" << endl;
98  }
99  if(single_pion_sources==0) {
100  out_stream << "# 1pi sources: All" << endl;
101  }
102  else if(single_pion_sources==1) {
103  out_stream << "# 1pi sources: P33(1232) resonance only" << endl;
104  }
105  else if(single_pion_sources==2) {
106  out_stream << "# 1pi sources: All resonances only" << endl;
107  }
108  out_stream << "# " << endl;
109  out_stream << "# Note:" << endl;
110  out_stream << "# - neutrino energy E in GeV, linear spacing between Emin = " << Enumin << " GeV, Emax = " << Enumax << " GeV " << endl;
111  out_stream << "# - cross sections in 1E-38 cm^2 " << endl;
112  out_stream << "# - quoted cross section is nuclear cross section divided with number of nucleons A" << endl;
113  out_stream << "# Columns:" << endl;
114  out_stream << "# | Energy | sig(nu_mu + A -> mu- 1pi+ X) | " << endl;
115 
116  out_stream << std::fixed << setprecision(6);
117 
118  //
119  // load event data
120  //
121 
122  TChain * chain = new TChain("gst");
123 
124  // loop over CC/NC cases
125  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
126  // loop over runs for current case
127  for(int ir=0; ir<kNRunsPerCase; ir++) {
128  // build filename
129  ostringstream filename;
130  int run_number = kRunNu1PI1[isample][iwkcur][ir];
131 
132  cout << "isample = " << isample << ", iwkcur = " << iwkcur << ", ir = " << ir << ", run = " << run_number << endl;
133 
134  filename << "../gst/gntp." << run_number << ".gst.root";
135  // add to chain
136  cout << "Adding " << filename.str() << " to event chain" << endl;
137  chain->Add(filename.str().c_str());
138  }
139  }
140 
141  //
142  // book histograms
143  //
144  TH1D * hst_cc1pip = new TH1D("hst_cc1pip","CC1pi+ events, Enu spectrum", nEnu, Enumin, Enumax);
145  TH1D * hst_allcc = new TH1D("hst_allcc", "all CC events, Enu spectrum", nEnu, Enumin, Enumax);
146 
147  //
148  // fill histograms
149  //
150 
151  chain->Draw("Ev>>hst_allcc", "cc","GOFF");
152 
153  if(stage==1) {
154  if(single_pion_sources==0) {
155  //all sources
156  chain->Draw("Ev>>hst_cc1pip","cc&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
157  }
158  else if(single_pion_sources==1) {
159  //P33(1232) only
160  chain->Draw("Ev>>hst_cc1pip","cc&&resid==0&&res&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
161  }
162  else if(single_pion_sources==2) {
163  //all resonances only
164  chain->Draw("Ev>>hst_cc1pip","cc&&res&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
165  }
166  }
167  else if(stage==0) {
168  if(single_pion_sources==0) {
169  //all sources
170  chain->Draw("Ev>>hst_cc1pip","cc&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
171  }
172  else if(single_pion_sources==1) {
173  //P33(1232) only
174  chain->Draw("Ev>>hst_cc1pip","cc&&resid==0&&res&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
175  }
176  else if(single_pion_sources==2) {
177  //all resonances only
178  chain->Draw("Ev>>hst_cc1pip","cc&&res&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
179  }
180  }
181 
182  cout << "CC1pi+ nevents: " << hst_cc1pip->GetEntries() << endl;
183  cout << "ALLCC nevents: " << hst_allcc->GetEntries() << endl;
184 
185  //
186  // compute sig(CC1pi+) and write out in txt file expected by NuINT organizers
187  // Also write out root graph in temp file to be re-used at later benchmarking calc
188  //
189 
190  double e[nEnu], sig[nEnu];
191 
192  for(int i=1; i <= hst_cc1pip->GetNbinsX(); i++) {
193 
194  double Enu = hst_cc1pip->GetBinCenter(i);
195 
196  double Ncc1pip = hst_cc1pip -> GetBinContent(i);
197  double Nallcc = hst_allcc -> GetBinContent(i);
198 
199  double sig_cc1pip=0;
200  if(Nallcc>0) { sig_cc1pip = (Ncc1pip/Nallcc) * sig_graph_totcc->Eval(Enu) / A; }
201 
202  out_stream << setw(15) << Enu << setw(15) << sig_cc1pip << endl;
203 
204  e [i-1] = Enu;
205  sig[i-1] = sig_cc1pip;
206  }
207 
208  out_stream.close();
209 
210  TFile ftmp("./cc1pip_tmp.root","recreate");
211  TGraph * grCC1pip = new TGraph(nEnu,e,sig);
212  grCC1pip->Write("CC1pip");
213  hst_allcc->Write();
214  hst_cc1pip->Write();
215 
216  // visual inspection
217  TCanvas * c1 = new TCanvas("c1","",20,20,500,500);
218  grCC1pip->Draw("alp");
219  c1->Update();
220 
221  ftmp.Close();
222 
223  delete chain;
224 }
const int kNRunsPerCase
Definition: nuint09_1pi1.C:21
string filename
Definition: shutoffs.py:106
const char * label
const int kNWCur
Definition: nuint09_1pi1.C:20
int kA[kNSamples]
Definition: nuint09_1pi1.C:44
chain
Check that an output directory exists.
const int kRunNu1PI1[kNSamples][kNWCur][kNRunsPerCase]
Definition: nuint09_1pi1.C:31
const int kNEvtPerRun
Definition: nuint09_1pi1.C:22
OStream cout
Definition: OStream.cxx:6
static const double A
Definition: Units.h:82
const int kNSamples
Definition: nuint09_1pi1.C:19
c1
Definition: demo5.py:24
Float_t e
Definition: plot.C:35
const char * kLabel[kNSamples]
Definition: nuint09_1pi1.C:24
void nuint09_1pi1(int isample, int single_pion_sources=0, int stage=1)
Definition: nuint09_1pi1.C:52