nuint09_1pi4.C
Go to the documentation of this file.
1 //
2 // NuINT09 Conference, Benchmark Calculations (GENIE contribution)
3 //
4 // 1PI.4:
5 // 1pi+ cross section at E_nu= 1.0 and 1.5 GeV as a function of the hadronic invarant mass.
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_1pi4(int isample, int single_pion_sources=0, int stage=1)
50 {
51  cout << " ***** running: 1PI.4" << 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 graph
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 nW = 60;
66  const double Wmin = 0.80;
67  const double Wmax = 2.50;
68 
69  // create output stream
70 
71  ostringstream out_filename;
72  out_filename << label;
73 
74  if (single_pion_sources==0) out_filename << ".1pi_4a.";
75  else if (single_pion_sources==1) out_filename << ".1pi_4b.";
76  else if (single_pion_sources==2) out_filename << ".1pi_4c.";
77 
78  if(stage==0) out_filename << "no_FSI.";
79 
80  out_filename << label << "dsig1pi_dW.data";
81 
82  ofstream out_stream(out_filename.str().c_str(), ios::out);
83 
84  // write out txt file
85 
86  out_stream << "# [" << label << "]" << endl;
87  out_stream << "# " << endl;
88  out_stream << "# [1PI.4]:" << endl;
89  out_stream << "# 1pi+ cross section at E_nu= 1.0 and 1.5 GeV as a function of the hadronic invariant mass W" << endl;
90  if(stage==0) {
91  out_stream << "# ***** NO FSI: The {X pi+} state is a primary hadronic state" << endl;
92  }
93  if(single_pion_sources==0) {
94  out_stream << "# 1pi sources: All" << endl;
95  }
96  else if(single_pion_sources==1) {
97  out_stream << "# 1pi sources: P33(1232) resonance only" << endl;
98  }
99  else if(single_pion_sources==2) {
100  out_stream << "# 1pi sources: All resonances only" << endl;
101  }
102  out_stream << "# " << endl;
103  out_stream << "# Note:" << endl;
104  out_stream << "# - W in GeV, linear spacing between Wmin = " << Wmin << " GeV, Wmax = " << Wmax << " GeV " << endl;
105  out_stream << "# - cross sections in 1E-38 cm^2 / GeV " << endl;
106  out_stream << "# - quoted cross section is nuclear cross section divided with number of nucleons A" << endl;
107  out_stream << "# Columns:" << endl;
108  out_stream << "# | W | dsig(numu A -> mu- 1pi+ X; Enu = 1.0 GeV) | dsig(numu A -> mu- 1pi+ X; Enu = 1.5 GeV) | " << endl;
109 
110  out_stream << std::fixed << setprecision(6);
111 
112  //
113  // load event data
114  //
115 
116  TChain * chain = new TChain("gst");
117 
118  // loop over CC/NC cases
119  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
120  // loop over energies
121  for(int ie=0; ie<kNEnergies; ie++) {
122  // loop over runs for current case
123  for(int ir=0; ir<kNRunsPerCase; ir++) {
124  // build filename
125  ostringstream filename;
126  int run_number = kRunNu1PI4[isample][iwkcur][ie][ir];
127  filename << "../gst/gntp." << run_number << ".gst.root";
128  // add to chain
129  cout << "Adding " << filename.str() << " to event chain" << endl;
130  chain->Add(filename.str().c_str());
131  }
132  }
133  }
134 
135  //
136  // get CC1pi+ cross sections at given energies for normalization purposes
137  //
138  double sig_cc1pip_1000MeV = sig_graph_cc1pip -> Eval (1.0);
139  double sig_cc1pip_1500MeV = sig_graph_cc1pip -> Eval (1.5);
140 
141  //
142  // book histograms
143  //
144  TH1D * hst_dsig_dW_1000MeV = new TH1D("hst_dsig_dW_1000MeV","dsig/dW, numu A -> mu- 1pi+ X, Enu=1.0 GeV", nW, Wmin, Wmax);
145  TH1D * hst_dsig_dW_1500MeV = new TH1D("hst_dsig_dW_1500MeV","dsig/dW, numu A -> mu- 1pi+ X, Enu=1.5 GeV", nW, Wmin, Wmax);
146 
147  //
148  // fill histograms
149  //
150  if(stage==1) {
151  if(single_pion_sources==0) {
152  // all sources
153  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
154  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
155  }
156  else if(single_pion_sources==1) {
157  // P33(1232) only
158  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&resid==0&&res&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
159  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&resid==0&&res&&Ev>1.49&&Ev<1.51&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
160  }
161  else if(single_pion_sources==2) {
162  // all resonances only
163  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&res&&Ev>0.99&&Ev<1.01&&pdgf==211&&nfpip==1&&nfpim==0&&nfpi0==0","GOFF");
164  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&res&&Ev>1.49&&Ev<1.51&&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("Ws>>hst_dsig_dW_1000MeV","cc&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
171  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
172  }
173  else if(single_pion_sources==1) {
174  // P33(1232) only
175  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&resid==0&&res&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
176  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&resid==0&&res&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
177  }
178  else if(single_pion_sources==2) {
179  // all resonances only
180  chain->Draw("Ws>>hst_dsig_dW_1000MeV","cc&&res&&Ev>0.99&&Ev<1.01&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
181  chain->Draw("Ws>>hst_dsig_dW_1500MeV","cc&&res&&Ev>1.49&&Ev<1.51&&pdgi==211&&nipip==1&&nipim==0&&nipi0==0","GOFF");
182  }
183  }
184 
185  //
186  // normalize
187  //
188  double norm_cc1pip_1000MeV = hst_dsig_dW_1000MeV -> Integral("width") / sig_cc1pip_1000MeV;
189  double norm_cc1pip_1500MeV = hst_dsig_dW_1500MeV -> Integral("width") / sig_cc1pip_1500MeV;
190 
191  if (norm_cc1pip_1000MeV > 0) hst_dsig_dW_1000MeV -> Scale(1./norm_cc1pip_1000MeV);
192  if (norm_cc1pip_1500MeV > 0) hst_dsig_dW_1500MeV -> Scale(1./norm_cc1pip_1500MeV);
193 
194  for(int i=1; i <= hst_dsig_dW_1000MeV->GetNbinsX(); i++) {
195 
196  double W = hst_dsig_dW_1000MeV -> GetBinCenter(i);
197 
198  double dsig_dW_1000MeV = hst_dsig_dW_1000MeV -> GetBinContent(i);
199  double dsig_dW_1500MeV = hst_dsig_dW_1500MeV -> GetBinContent(i);
200 
201  dsig_dW_1000MeV = TMath::Max(0., dsig_dW_1000MeV);
202  dsig_dW_1500MeV = TMath::Max(0., dsig_dW_1500MeV);
203 
204  out_stream << setw(15) << W
205  << setw(15) << dsig_dW_1000MeV
206  << setw(15) << dsig_dW_1500MeV
207  << endl;
208  }
209 
210  out_stream.close();
211 
212  // visual inspection
213  //TCanvas * c1 = new TCanvas("c1","",20,20,500,500);
214  //...
215  //c1->Update();
216 
217  delete chain;
218  fsig.Close();
219 }
void nuint09_1pi4(int isample, int single_pion_sources=0, int stage=1)
Definition: nuint09_1pi4.C:49
double Integral(const Spectrum &s, const double pot, int cvnregion)
double Wmax[kNWBins]
string filename
Definition: shutoffs.py:106
const int kNEvtPerRun
Definition: nuint09_1pi4.C:23
const int kRunNu1PI4[kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
Definition: nuint09_1pi4.C:32
const char * label
const int kNRunsPerCase
Definition: nuint09_1pi4.C:22
chain
Check that an output directory exists.
const int kNWCur
Definition: nuint09_1pi4.C:20
const char * kLabel[kNSamples]
Definition: nuint09_1pi4.C:25
OStream cout
Definition: OStream.cxx:6
const int kNSamples
Definition: nuint09_1pi4.C:19
double Wmin[kNWBins]
const int kNEnergies
Definition: nuint09_1pi4.C:21
simulatedPeak Scale(1/simulationNormalisationFactor)
#define W(x)