nuint09_qel2.C
Go to the documentation of this file.
1 //
2 // NuINT09 Conference, Benchmark Calculations (GENIE contribution)
3 //
4 // QEL.2:
5 // Cross section as a function of proton kinetic energy at E_nu= 0.5, 1.0 GeV.
6 //
7 // Inputs:
8 // - sample id: 0 (nu_mu+C12), 1 (nu_mu+O16), 2 (nu_mu+Fe56)
9 // - include_fsi: 0 (no), 1 (yes)
10 //
11 // Costas Andreopoulos, STFC / Rutherford Appleton Laboratory
12 //
13 #include <iomanip>
14 
15 //
16 // consts
17 //
18 const int kNSamples = 3;
19 const int kNWCur = 1;
20 const int kNEnergies = 2;
21 const int kNRunsPerCase = 5;
22 const int kNEvtPerRun = 100000;
23 
24 const char * kLabel[kNSamples] =
25 {
26  /* 0 */ "nu_mu_C12",
27  /* 1 */ "nu_mu_O16",
28  /* 2 */ "nu_mu_Fe56"
29 };
30 
32 {
33  /* indices : sample ; cc/nc ; energy */
34  {
35  /* 0,0,0 (nu_mu C12, CC, 0.5 GeV) */ { {200100, 200101, 200102, 200103, 200104},
36  /* 0,0,1 (nu_mu C12, CC, 1.0 GeV) */ {200200, 200201, 200202, 200203, 200204} }
37  },
38  {
39  /* 1,0,0 (nu_mu O16, CC, 0.5 GeV) */ { {210100, 210101, 210102, 210103, 210104},
40  /* 1,0,1 (nu_mu O16, CC, 1.0 GeV) */ {210200, 210201, 210202, 210203, 210204} }
41  },
42  {
43  /* 2,0,0 (nu_mu Fe56, CC, 0.5 GeV) */ { {220100, 220101, 220102, 220103, 220104},
44  /* 2,0,1 (nu_mu Fe56, CC, 1.0 GeV) */ {220200, 220201, 220202, 220203, 220204} }
45  }
46 };
47 int kZ[3] =
48 {
49  // Z for nuclear target at each sample
50  /* 0 */ 6,
51  /* 1 */ 8,
52  /* 2 */ 26
53 };
54 int kA[3] =
55 {
56  // A for nuclear target at each sample
57  /* 0 */ 12,
58  /* 1 */ 16,
59  /* 2 */ 56
60 };
61 
62 void nuint09_qel2(int isample, int include_fsi=0)
63 {
64  cout << " ***** running: QEL.2" << endl;
65  if(include_fsi==0) { cout << "-FSI" << endl; }
66  else
67  if(include_fsi==1) { cout << "+FSI" << endl; }
68  else
69  return;
70 
71  if(isample<0 || isample >= kNSamples) return;
72 
73  const char * label = kLabel[isample];
74 
75  int A = kA[isample];
76  int Z = kZ[isample];
77  int N = A-Z;
78 
79  // get cross section graphs
80 
81  TFile fsig("../sig/splines.root","read");
82  TDirectory * sig_dir = (TDirectory *) fsig.Get(label);
83 
84  TGraph * sig_graph_qelcc = (TGraph*) sig_dir->Get("qel_cc_n");
85 
86  // range & spacing
87 
88  const int nKEp = 60;
89  const double KEpmin = 0.01;
90  const double KEpmax = 1.50;
91 
92  // create output stream
93 
94  ostringstream out_filename;
95  if(include_fsi==0) { out_filename << label << ".qel_2.dsig_dKEp.data"; }
96  else { out_filename << label << ".qel_2.dsig_dKEp_withFSI.data"; }
97 
98  ofstream out_stream(out_filename.str().c_str(), ios::out);
99 
100  // write out txt file
101 
102  out_stream << "# [" << label << "]" << endl;
103  out_stream << "# " << endl;
104  out_stream << "# [QEL.2]:" << endl;
105  out_stream << "# Cross section as a function of f/s proton kinetic energy at E_nu= 0.5, 1.0 GeV." << endl;
106  out_stream << "# " << endl;
107  out_stream << "# Note:" << endl;
108  if(include_fsi==1)
109  {
110  out_stream << "# - INCLUDING FSI " << endl;
111  }
112  out_stream << "# - proton energies are _kinetic_ energies " << endl;
113  out_stream << "# - proton kinetic energy KE in GeV, between KEmin = " << KEpmin << " GeV, KEmax = " << KEpmax << " GeV " << endl;
114  out_stream << "# - cross sections in 1E-38 cm^2 / GeV" << endl;
115  out_stream << "# - quoted cross section is nuclear cross section per nucleon contributing in the scattering (eg only neutrons for nu_mu QELCC)" << endl;
116  out_stream << "# Columns:" << endl;
117  out_stream << "# | KE(p) | sig(QELCC; Ev=0.5GeV) | sig(QELCC; Ev=1.0GeV) | " << endl;
118 
119  out_stream << std::fixed << setprecision(6);
120 
121  //
122  // load event data
123  //
124 
125  TChain * chain = new TChain("gst");
126 
127  // loop over CC/NC cases
128  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
129  // loop over energies
130  for(int ie=0; ie<kNEnergies; ie++) {
131  // loop over runs for current case
132  for(int ir=0; ir<kNRunsPerCase; ir++) {
133  // build filename
134  ostringstream filename;
135  int run_number = kRunNu[isample][iwkcur][ie][ir];
136  filename << "../gst/gntp." << run_number << ".gst.root";
137  // add to chain
138  cout << "Adding " << filename.str() << " to event chain" << endl;
139  chain->Add(filename.str().c_str());
140  }
141  }
142  }
143 
144 
145  //
146  // get QELCC cross sections at given energies for normalization purposes
147  //
148  double sig_qelcc_0500MeV = sig_graph_qelcc->Eval(0.5) / N;
149  double sig_qelcc_1000MeV = sig_graph_qelcc->Eval(1.0) / N;
150 
151  //
152  // book histograms
153  //
154  TH1D * hst_dsig_dKEp_qelcc_0500MeV = new TH1D("hst_dsig_dKEp_qelcc_0500MeV","dsig/dKEp, nu_mu QEL CC, Enu=0.5 GeV", nKEp, KEpmin, KEpmax);
155  TH1D * hst_dsig_dKEp_qelcc_1000MeV = new TH1D("hst_dsig_dKEp_qelcc_1000MeV","dsig/dKEp, nu_mu QEL CC, Enu=1.0 GeV", nKEp, KEpmin, KEpmax);
156 
157  //
158  // fill histograms
159  //
160  if(include_fsi==0) {
161  chain->Draw("(Ei-0.938)>>hst_dsig_dKEp_qelcc_0500MeV","qel&&cc&&Ev>0.49&&Ev<0.51&&pdgi==2212","GOFF");
162  chain->Draw("(Ei-0.938)>>hst_dsig_dKEp_qelcc_1000MeV","qel&&cc&&Ev>0.99&&Ev<1.01&&pdgi==2212","GOFF");
163  } else {
164  chain->Draw("(Ef-0.938)>>hst_dsig_dKEp_qelcc_0500MeV","qel&&cc&&Ev>0.49&&Ev<0.51&&pdgf==2212","GOFF");
165  chain->Draw("(Ef-0.938)>>hst_dsig_dKEp_qelcc_1000MeV","qel&&cc&&Ev>0.99&&Ev<1.01&&pdgf==2212","GOFF");
166  }
167 
168  //
169  // normalize
170  //
171  double norm_qelcc_0500MeV = hst_dsig_dKEp_qelcc_0500MeV -> Integral("width") / sig_qelcc_0500MeV;
172  double norm_qelcc_1000MeV = hst_dsig_dKEp_qelcc_1000MeV -> Integral("width") / sig_qelcc_1000MeV;
173 
174  if (norm_qelcc_0500MeV > 0) hst_dsig_dKEp_qelcc_0500MeV -> Scale(1./norm_qelcc_0500MeV);
175  if (norm_qelcc_1000MeV > 0) hst_dsig_dKEp_qelcc_1000MeV -> Scale(1./norm_qelcc_1000MeV);
176 
177  for(int i = 1; i <= hst_dsig_dKEp_qelcc_1000MeV->GetNbinsX(); i++) {
178 
179  double KEp = hst_dsig_dKEp_qelcc_1000MeV->GetBinCenter(i);
180 
181  double dsig_dKEp_qelcc_0500MeV = hst_dsig_dKEp_qelcc_0500MeV -> GetBinContent(i);
182  double dsig_dKEp_qelcc_1000MeV = hst_dsig_dKEp_qelcc_1000MeV -> GetBinContent(i);
183 
184  dsig_dKEp_qelcc_0500MeV = TMath::Max(0., dsig_dKEp_qelcc_0500MeV);
185  dsig_dKEp_qelcc_1000MeV = TMath::Max(0., dsig_dKEp_qelcc_1000MeV);
186 
187  out_stream << setw(15) << KEp
188  << setw(15) << dsig_dKEp_qelcc_0500MeV
189  << setw(15) << dsig_dKEp_qelcc_1000MeV
190  << endl;
191  }
192 
193  out_stream.close();
194 
195  delete chain;
196 }
double Integral(const Spectrum &s, const double pot, int cvnregion)
string filename
Definition: shutoffs.py:106
const int kRunNu[kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
Definition: nuint09_qel2.C:31
const int kNEvtPerRun
Definition: nuint09_qel2.C:22
const char * label
Float_t Z
Definition: plot.C:38
const int kNEnergies
Definition: nuint09_qel2.C:20
const int kNRunsPerCase
Definition: nuint09_qel2.C:21
int kZ[3]
Definition: nuint09_qel2.C:47
chain
Check that an output directory exists.
const char * kLabel[kNSamples]
Definition: nuint09_qel2.C:24
OStream cout
Definition: OStream.cxx:6
static const double A
Definition: Units.h:82
const int kNSamples
Definition: nuint09_qel2.C:18
void nuint09_qel2(int isample, int include_fsi=0)
Definition: nuint09_qel2.C:62
simulatedPeak Scale(1/simulationNormalisationFactor)
const int kNWCur
Definition: nuint09_qel2.C:19
int kA[3]
Definition: nuint09_qel2.C:54