nuint09_qel5.C
Go to the documentation of this file.
1 //
2 // NuINT09 Conference, Benchmark Calculations (GENIE contribution)
3 //
4 // QEL.5:
5 // dSigma / dOmega_mu dKE_mu at E_nu = 0.5 and 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_qel5(int isample)
62 {
63  cout << " ***** running: QEL.5" << 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  const int ncosth = 30;
87  const double costhmin = -1;
88  const double costhmax = +1;
89 
90  // create output stream
91 
92  ostringstream out_filename;
93  out_filename << label << ".qel_5.d2sig_dKEmudOmega.data";
94  ofstream out_stream(out_filename.str().c_str(), ios::out);
95 
96  // write out txt file
97 
98  out_stream << "# [" << label << "]" << endl;
99  out_stream << "# " << endl;
100  out_stream << "# [QEL.5]:" << endl;
101  out_stream << "# dSigma / dOmega_mu dKE_mu at E_nu = 0.5 and 1.0 GeV" << endl;
102  out_stream << "# " << endl;
103  out_stream << "# Note:" << endl;
104  out_stream << "# - muon energies are _kinetic_ energies " << endl;
105  out_stream << "# - muon kinetic energy KE in GeV, between KEmin = " << KEmumin << " GeV, KEmax = " << KEmumax << " GeV " << endl;
106  out_stream << "# - cross sections in 1E-38 cm^2 /GeV /sterad" << endl;
107  out_stream << "# - quoted cross section is nuclear cross section per nucleon contributing in the scattering (eg only neutrons for nu_mu QELCC)" << endl;
108  out_stream << "# Columns:" << endl;
109  out_stream << "# | KE(mu) | cos(theta_mu) | sig(QELCC; Ev=0.5GeV) | sig(QELCC; Ev=1.0GeV) | " << endl;
110 
111  out_stream << std::fixed << setprecision(6);
112 
113  //
114  // load event data
115  //
116 
117  TChain * chain = new TChain("gst");
118 
119  // loop over CC/NC cases
120  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
121  // loop over energies
122  for(int ie=0; ie<kNEnergies; ie++) {
123  // loop over runs for current case
124  for(int ir=0; ir<kNRunsPerCase; ir++) {
125  // build filename
126  ostringstream filename;
127  int run_number = kRunNu[isample][iwkcur][ie][ir];
128  filename << "../gst/gntp." << run_number << ".gst.root";
129  // add to chain
130  cout << "Adding " << filename.str() << " to event chain" << endl;
131  chain->Add(filename.str().c_str());
132  }
133  }
134  }
135 
136  //
137  // get QELCC cross sections at given energies for normalization purposes
138  //
139  double sig_qelcc_0500MeV = sig_graph_qelcc->Eval(0.5) / N;
140  double sig_qelcc_1000MeV = sig_graph_qelcc->Eval(1.0) / N;
141 
142  //
143  // book histograms
144  //
145  TH2D * hst_d2sig_dKEmudOmg_qelcc_0500MeV = new TH2D("hst_d2sig_dKEmudOmg_qelcc_0500MeV",
146  "dsig/d2KEmudOmega, nu_mu QEL CC, Enu=0.5 GeV", nKEmu, KEmumin, KEmumax, ncosth, costhmin, costhmax);
147  TH2D * hst_d2sig_dKEmudOmg_qelcc_1000MeV = new TH2D("hst_d2sig_dKEmudOmg_qelcc_1000MeV",
148  "dsig/d2KEmudOmega, nu_mu QEL CC, Enu=1.0 GeV", nKEmu, KEmumin, KEmumax, ncosth, costhmin, costhmax);
149 
150  //
151  // fill histograms
152  //
153  chain->Draw("pzl/sqrt(pzl*pzl+pyl*pyl+pxl*pxl):(El-0.105658)>>hst_d2sig_dKEmudOmg_qelcc_0500MeV","qel&&cc&&Ev>0.49&&Ev<0.51","GOFF");
154  chain->Draw("pzl/sqrt(pzl*pzl+pyl*pyl+pxl*pxl):(El-0.105658)>>hst_d2sig_dKEmudOmg_qelcc_1000MeV","qel&&cc&&Ev>0.99&&Ev<1.01","GOFF");
155 
156  //
157  // normalize
158  //
159  double norm_qelcc_0500MeV = hst_d2sig_dKEmudOmg_qelcc_0500MeV -> Integral("width") * 2*TMath::Pi() / sig_qelcc_0500MeV;
160  double norm_qelcc_1000MeV = hst_d2sig_dKEmudOmg_qelcc_1000MeV -> Integral("width") * 2*TMath::Pi() / sig_qelcc_1000MeV;
161 
162  if (norm_qelcc_0500MeV > 0) hst_d2sig_dKEmudOmg_qelcc_0500MeV -> Scale(1./norm_qelcc_0500MeV);
163  if (norm_qelcc_1000MeV > 0) hst_d2sig_dKEmudOmg_qelcc_1000MeV -> Scale(1./norm_qelcc_1000MeV);
164 
165  for(int i = 1; i <= hst_d2sig_dKEmudOmg_qelcc_1000MeV->GetNbinsX(); i++) {
166  for(int j = 1; j <= hst_d2sig_dKEmudOmg_qelcc_1000MeV->GetNbinsY(); j++) {
167 
168  double KEmu = hst_d2sig_dKEmudOmg_qelcc_1000MeV -> GetXaxis() -> GetBinCenter(i);
169  double costhmu = hst_d2sig_dKEmudOmg_qelcc_1000MeV -> GetYaxis() -> GetBinCenter(j);
170 
171  double d2sig_dKEmudOmg_qelcc_0500MeV = hst_d2sig_dKEmudOmg_qelcc_0500MeV -> GetBinContent(i,j);
172  double d2sig_dKEmudOmg_qelcc_1000MeV = hst_d2sig_dKEmudOmg_qelcc_1000MeV -> GetBinContent(i,j);
173 
174  d2sig_dKEmudOmg_qelcc_0500MeV = TMath::Max(0., d2sig_dKEmudOmg_qelcc_0500MeV);
175  d2sig_dKEmudOmg_qelcc_1000MeV = TMath::Max(0., d2sig_dKEmudOmg_qelcc_1000MeV);
176 
177  out_stream << setw(15) << KEmu
178  << setw(15) << costhmu
179  << setw(15) << d2sig_dKEmudOmg_qelcc_0500MeV
180  << setw(15) << d2sig_dKEmudOmg_qelcc_1000MeV
181  << endl;
182  }//costh
183  }//E
184 
185  out_stream.close();
186 
187  delete chain;
188 }
const int kNWCur
Definition: nuint09_qel5.C:18
correl_xv GetYaxis() -> SetDecimals()
int kA[3]
Definition: nuint09_qel5.C:53
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_qel5.C:30
const int kNEvtPerRun
Definition: nuint09_qel5.C:21
void nuint09_qel5(int isample)
Definition: nuint09_qel5.C:61
const char * label
Float_t Z
Definition: plot.C:38
int kZ[3]
Definition: nuint09_qel5.C:46
const int kNEnergies
Definition: nuint09_qel5.C:19
correl_xv GetXaxis() -> SetDecimals()
const char * kLabel[kNSamples]
Definition: nuint09_qel5.C:23
const int kNRunsPerCase
Definition: nuint09_qel5.C:20
chain
Check that an output directory exists.
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
static const double A
Definition: Units.h:82
simulatedPeak Scale(1/simulationNormalisationFactor)
const int kNSamples
Definition: nuint09_qel5.C:17