nuint09_1pi3.C
Go to the documentation of this file.
1 //
2 // NuINT09 Conference, Benchmark Calculations (GENIE contribution)
3 //
4 // 1PI.3:
5 // d^Sigma / dOmega_pi+ dKE_pi+ at E_nu= 1.0 and 1.5 GeV
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 kNEnergies = 2;
22 const int kNRunsPerCase = 5;
23 const int kNEvtPerRun = 100000;
24 
25 const char * kLabel[kNSamples] =
26 {
27  // available samples
28  /* 0 */ "nu_mu_C12",
29  /* 1 */ "nu_mu_O16",
30  /* 2 */ "nu_mu_Fe56"
31 };
33 {
34  /* indices : sample ; cc/nc ; energy */
35  {
36  /* 0,0,0 (nu_mu C12, CC, 1.0 GeV) */ { {900200, 900201, 900202, 900203, 900204},
37  /* 0,0,1 (nu_mu C12, CC, 1.5 GeV) */ {900300, 900301, 900302, 900303, 900304} }
38  },
39  {
40  /* 1,0,0 (nu_mu O16, CC, 1.0 GeV) */ { {910200, 910201, 910202, 910203, 910204},
41  /* 1,0,1 (nu_mu O16, CC, 1.5 GeV) */ {910300, 910301, 910302, 910303, 910304} }
42  },
43  {
44  /* 2,0,0 (nu_mu Fe56, CC, 1.0 GeV) */ { {920200, 920201, 920202, 920203, 920204},
45  /* 2,0,1 (nu_mu Fe56, CC, 1.5 GeV) */ {920300, 920301, 920302, 920303, 920304} }
46  }
47 };
48 
49 void nuint09_1pi3(int isample, int single_pion_sources=0, int stage=1)
50 {
51  cout << " ***** running: 1PI.3" << endl;
52 
53  if(isample<0 || isample >= kNSamples) return;
54  if(single_pion_sources<0 || single_pion_sources>2) return;
55 
56  const char * label = kLabel[isample];
57 
58  // get cross section graphs
59 
60  TFile fsig("./cc1pip_tmp.root","read"); // generated at [1PI.1]
61  TGraph * sig_graph_cc1pip = (TGraph*) fsig.Get("CC1pip");
62 
63  // range & spacing
64 
65  const int nKEpip = 60;
66  const double KEpipmin = 0.00;
67  const double KEpipmax = 1.50;
68 
69  const int ncosth = 30;
70  const double costhmin = -1;
71  const double costhmax = +1;
72 
73  // create output stream
74 
75  ostringstream out_filename;
76  out_filename << label;
77 
78  if (single_pion_sources==0) out_filename << ".1pi_3a.";
79  else if (single_pion_sources==1) out_filename << ".1pi_3b.";
80  else if (single_pion_sources==2) out_filename << ".1pi_3c.";
81 
82  if(stage==0) out_filename << "no_FSI.";
83 
84  out_filename << label << "d2sig1pi_dKEpidOmega.data";
85 
86  ofstream out_stream(out_filename.str().c_str(), ios::out);
87 
88  // write out txt file
89 
90  out_stream << "# [" << label << "]" << endl;
91  out_stream << "# " << endl;
92  out_stream << "# [1PI.3]:" << endl;
93  out_stream << "# d^Sigma / dOmega_pi+ dKE_pi+ at E_nu= 1.0 and 1.5 GeV" << endl;
94  if(stage==0) {
95  out_stream << "# ***** NO FSI: The {X pi+} state is a primary hadronic state" << endl;
96  }
97  if(single_pion_sources==0) {
98  out_stream << "# 1pi sources: All" << endl;
99  }
100  else if(single_pion_sources==1) {
101  out_stream << "# 1pi sources: P33(1232) resonance only" << endl;
102  }
103  else if(single_pion_sources==2) {
104  out_stream << "# 1pi sources: All resonances only" << endl;
105  }
106  out_stream << "# " << endl;
107  out_stream << "# Note:" << endl;
108  out_stream << "# - pi+ kinetic energy KE in GeV, linear spacing between KEmin = " << KEpipmin << " GeV, KEmax = " << KEpipmax << " GeV " << endl;
109  out_stream << "# - cross sections in 1E-38 cm^2 /GeV /sr" << endl;
110  out_stream << "# - quoted cross section is nuclear cross section divided with number of nucleons A" << endl;
111  out_stream << "# Columns:" << endl;
112  out_stream << "# | KE(pi+) | cos(theta_pi+) | dsig(numu A -> mu- 1pi+ X; Enu = 1.0 GeV) | dsig(numu A -> mu- 1pi+ X; Enu = 1.5 GeV) | " << endl;
113 
114  out_stream << std::fixed << setprecision(6);
115 
116  //
117  // load event data
118  //
119 
120  TChain * chain = new TChain("gst");
121 
122  // loop over CC/NC cases
123  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
124  // loop over energies
125  for(int ie=0; ie<kNEnergies; ie++) {
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 = kRunNu1PI3[isample][iwkcur][ie][ir];
131  filename << "../gst/gntp." << run_number << ".gst.root";
132  // add to chain
133  cout << "Adding " << filename.str() << " to event chain" << endl;
134  chain->Add(filename.str().c_str());
135  }
136  }
137  }
138 
139  //
140  // get CC1pi+ cross sections at given energies for normalization purposes
141  //
142  double sig_cc1pip_1000MeV = sig_graph_cc1pip -> Eval (1.0);
143  double sig_cc1pip_1500MeV = sig_graph_cc1pip -> Eval (1.5);
144 
145  //
146  // book histograms
147  //
148  TH2D * hst_d2sig_dKEpipdOmg_1000MeV = new TH2D("hst_d2sig_dKEpipdOmg_1000MeV",
149  "d2sig / dKEpi+ dOmega, numu A -> mu- 1pi+ X, Enu=1.0 GeV", nKEpip, KEpipmin, KEpipmax, ncosth, costhmin, costhmax);
150  TH2D * hst_d2sig_dKEpipdOmg_1500MeV = new TH2D("hst_d2sig_dKEpipdOmg_1500MeV",
151  "d2sig / dKEpi+ dOmega, numu A -> mu- 1pi+ X, Enu=1.5 GeV", nKEpip, KEpipmin, KEpipmax, ncosth, costhmin, costhmax);
152 
153  //
154  // fill histograms
155  //
156  if(stage==1) {
157  if(single_pion_sources==0) {
158  // all sources
159  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
160  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
161  }
162  else if(single_pion_sources==1) {
163  // P33(1232) only
164  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&resid==0&&res&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
165  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&resid==0&&res&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
166  }
167  else if(single_pion_sources==2) {
168  // all resonances only
169  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&res&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
170  chain->Draw("pzf/sqrt(pzf*pzf+pyf*pyf+pxf*pxf):(Ef-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&res&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
171  }
172  }
173  else if(stage==0) {
174  if(single_pion_sources==0) {
175  // all sources
176  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
177  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
178  }
179  else if(single_pion_sources==1) {
180  // P33(1232) only
181  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&resid==0&&res&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
182  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&resid==0&&res&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
183  }
184  else if(single_pion_sources==2) {
185  // all resonances only
186  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1000MeV","cc&&res&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
187  chain->Draw("pzi/sqrt(pzi*pzi+pyi*pyi+pxi*pxi):(Ei-0.139)>>hst_d2sig_dKEpipdOmg_1500MeV","cc&&res&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
188  }
189  }
190 
191  //
192  // normalize
193  //
194  double norm_cc1pip_1000MeV = hst_d2sig_dKEpipdOmg_1000MeV -> Integral("width") * 2*TMath::Pi() / sig_cc1pip_1000MeV;
195  double norm_cc1pip_1500MeV = hst_d2sig_dKEpipdOmg_1500MeV -> Integral("width") * 2*TMath::Pi() / sig_cc1pip_1500MeV;
196 
197  if (norm_cc1pip_1000MeV > 0) hst_d2sig_dKEpipdOmg_1000MeV -> Scale(1./norm_cc1pip_1000MeV);
198  if (norm_cc1pip_1500MeV > 0) hst_d2sig_dKEpipdOmg_1500MeV -> Scale(1./norm_cc1pip_1500MeV);
199 
200  for(int i = 1; i <= hst_d2sig_dKEpipdOmg_1500MeV->GetNbinsX(); i++) {
201  for(int j = 1; j <= hst_d2sig_dKEpipdOmg_1500MeV->GetNbinsY(); j++) {
202 
203  double KEpip = hst_d2sig_dKEpipdOmg_1500MeV -> GetXaxis() -> GetBinCenter(i);
204  double costhpip = hst_d2sig_dKEpipdOmg_1500MeV -> GetYaxis() -> GetBinCenter(j);
205 
206  double d2sig_dKEpipdOmg_1000MeV = hst_d2sig_dKEpipdOmg_1000MeV -> GetBinContent(i,j);
207  double d2sig_dKEpipdOmg_1500MeV = hst_d2sig_dKEpipdOmg_1500MeV -> GetBinContent(i,j);
208 
209  d2sig_dKEpipdOmg_1000MeV = TMath::Max(0., d2sig_dKEpipdOmg_1000MeV);
210  d2sig_dKEpipdOmg_1500MeV = TMath::Max(0., d2sig_dKEpipdOmg_1500MeV);
211 
212  out_stream << setw(15) << KEpip
213  << setw(15) << costhpip
214  << setw(15) << d2sig_dKEpipdOmg_1000MeV
215  << setw(15) << d2sig_dKEpipdOmg_1500MeV
216  << endl;
217  }
218  }
219 
220  out_stream.close();
221 
222  // visual inspection
223  //TCanvas * c1 = new TCanvas("c1","",20,20,500,500);
224  //...
225  //c1->Update();
226 
227  fsig.Close();
228  delete chain;
229 }
correl_xv GetYaxis() -> SetDecimals()
void nuint09_1pi3(int isample, int single_pion_sources=0, int stage=1)
Definition: nuint09_1pi3.C:49
double Integral(const Spectrum &s, const double pot, int cvnregion)
const int kNRunsPerCase
Definition: nuint09_1pi3.C:22
const int kRunNu1PI3[kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
Definition: nuint09_1pi3.C:32
string filename
Definition: shutoffs.py:106
const int kNWCur
Definition: nuint09_1pi3.C:20
const char * label
correl_xv GetXaxis() -> SetDecimals()
chain
Check that an output directory exists.
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
const int kNEvtPerRun
Definition: nuint09_1pi3.C:23
const int kNEnergies
Definition: nuint09_1pi3.C:21
simulatedPeak Scale(1/simulationNormalisationFactor)
const int kNSamples
Definition: nuint09_1pi3.C:19
const char * kLabel[kNSamples]
Definition: nuint09_1pi3.C:25