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