nuint09_coh2.C
Go to the documentation of this file.
1 //
2 // NuINT09 Conference, Benchmark Calculations (GENIE contribution)
3 //
4 // COH.2:
5 // Cross section as a function of pi0 or pi+ energy at E_nu=0.5, 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 //
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 = 2;
19 const int kNEnergies = 3;
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) */ { {100100, 100101, 100102, 100103, 100104},
35  /* 0,0,1 (nu_mu C12, CC, 1.0 GeV) */ {100200, 100201, 100202, 100203, 100204},
36  /* 0,0,2 (nu_mu C12, CC, 1.5 GeV) */ {100300, 100301, 100302, 100303, 100304} },
37  /* 0,1,0 (nu_mu C12, NC, 0.5 GeV) */ { {101100, 101101, 101102, 101103, 101104},
38  /* 0,1,1 (nu_mu C12, NC, 1.0 GeV) */ {101200, 101201, 101202, 101203, 101204},
39  /* 0,1,2 (nu_mu C12, NC, 1.5 GeV) */ {101300, 101301, 101302, 101303, 101304} }
40  },
41  {
42  /* 1,0,0 (nu_mu O16, CC, 0.5 GeV) */ { {110100, 110101, 110102, 110103, 110104},
43  /* 1,0,1 (nu_mu O16, CC, 1.0 GeV) */ {110200, 110201, 110202, 110203, 110204},
44  /* 1,0,2 (nu_mu O16, CC, 1.5 GeV) */ {110300, 110301, 110302, 110303, 110304} },
45  /* 1,1,0 (nu_mu O16, NC, 0.5 GeV) */ { {111100, 111101, 111102, 111103, 111104},
46  /* 1,1,1 (nu_mu O16, NC, 1.0 GeV) */ {111200, 111201, 111202, 111203, 111204},
47  /* 1,1,2 (nu_mu O16, NC, 1.5 GeV) */ {111300, 111301, 111302, 111303, 111304} }
48  },
49  {
50  /* 2,0,0 (nu_mu Fe56, CC, 0.5 GeV) */ { {120100, 120101, 120102, 120103, 120104},
51  /* 2,0,1 (nu_mu Fe56, CC, 1.0 GeV) */ {120200, 120201, 120202, 120203, 120204},
52  /* 2,0,2 (nu_mu Fe56, CC, 1.5 GeV) */ {120300, 120301, 120302, 120303, 120304} },
53  /* 2,1,0 (nu_mu Fe56, NC, 0.5 GeV) */ { {121100, 121101, 121102, 121103, 121104},
54  /* 2,1,1 (nu_mu Fe56, NC, 1.0 GeV) */ {121200, 121201, 121202, 121203, 121204},
55  /* 2,1,2 (nu_mu Fe56, NC, 1.5 GeV) */ {121300, 121301, 121302, 121303, 121304} }
56  }
57 };
58 
59 
60 void nuint09_coh2(int isample)
61 {
62  cout << " ***** running: COH.2" << endl;
63 
64  if(isample<0 || isample >= kNSamples) return;
65 
66  const char * label = kLabel[isample];
67 
68  // get cross section graphs
69 
70  TFile fsig("../sig/splines.root","read");
71  TDirectory * sig_dir = (TDirectory *) fsig.Get(label);
72 
73  TGraph * sig_graph_cohcc = (TGraph*) sig_dir->Get("coh_cc");
74  TGraph * sig_graph_cohnc = (TGraph*) sig_dir->Get("coh_nc");
75 
76  // range & spacing
77 
78  const int nKEpi = 60;
79  const double KEpimin = 0.01;
80  const double KEpimax = 1.50;
81 
82  // create output stream
83 
84  ostringstream out_filename;
85  out_filename << label << ".coh_2.dsig_dKEpi.data";
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 << "# [COH.2]:" << endl;
93  out_stream << "# Cross section as a function of pi0 or pi+ energy at E_nu=0.5, 1.0 and 1.5 GeV" << endl;
94  out_stream << "# " << endl;
95  out_stream << "# Note:" << endl;
96  out_stream << "# - pion energies are _kinetic_ energies " << endl;
97  out_stream << "# - pion kinetic energy KE in GeV, between KEmin = " << KEpimin << " GeV, KEmax = " << KEpimax << " GeV " << endl;
98  out_stream << "# - cross sections in 1E-38 cm^2 / GeV" << endl;
99  out_stream << "# - for coherent scattering we quote _nuclear_ cross section " << endl;
100  out_stream << "# Columns:" << endl;
101  out_stream << "# | KE(pi) | sig(coh; CC; Ev=0.5GeV) | sig(coh; NC; Ev=0.5GeV) | sig(coh; CC; Ev=1.0GeV) | sig(coh; NC; Ev=1.0GeV) | sig(coh; CC; Ev=1.5GeV) | sig(coh; NC; Ev=1.5GeV) " << endl;
102 
103  out_stream << std::fixed << setprecision(6);
104 
105  //
106  // load event data
107  //
108 
109  TChain * chain = new TChain("gst");
110 
111  // loop over CC/NC cases
112  for(int iwkcur=0; iwkcur<kNWCur; iwkcur++) {
113  // loop over energies
114  for(int ie=0; ie<kNEnergies; ie++) {
115  // loop over runs for current case
116  for(int ir=0; ir<kNRunsPerCase; ir++) {
117  // build filename
118  ostringstream filename;
119  int run_number = kRunNu[isample][iwkcur][ie][ir];
120  filename << "../gst/gntp." << run_number << ".gst.root";
121  // add to chain
122  cout << "Adding " << filename.str() << " to event chain" << endl;
123  chain->Add(filename.str().c_str());
124  }
125  }
126  }
127 
128 
129  //
130  // get COH pi cross sections at given energies for normalization purposes
131  //
132  double sig_cohcc_0500MeV = sig_graph_cohcc->Eval(0.5);
133  double sig_cohnc_0500MeV = sig_graph_cohnc->Eval(0.5);
134  double sig_cohcc_1000MeV = sig_graph_cohcc->Eval(1.0);
135  double sig_cohnc_1000MeV = sig_graph_cohnc->Eval(1.0);
136  double sig_cohcc_1500MeV = sig_graph_cohcc->Eval(1.5);
137  double sig_cohnc_1500MeV = sig_graph_cohnc->Eval(1.5);
138 
139  //
140  // book histograms
141  //
142  TH1D * hst_dsig_dKEpi_cohcc_0500MeV = new TH1D("hst_dsig_dKEpi_cohcc_0500MeV","dsig/dKEpi, COH CC, Enu=0.5 GeV", nKEpi, KEpimin, KEpimax);
143  TH1D * hst_dsig_dKEpi_cohnc_0500MeV = new TH1D("hst_dsig_dKEpi_cohnc_0500MeV","dsig/dKEpi, COH NC, Enu=0.5 GeV", nKEpi, KEpimin, KEpimax);
144  TH1D * hst_dsig_dKEpi_cohcc_1000MeV = new TH1D("hst_dsig_dKEpi_cohcc_1000MeV","dsig/dKEpi, COH CC, Enu=1.0 GeV", nKEpi, KEpimin, KEpimax);
145  TH1D * hst_dsig_dKEpi_cohnc_1000MeV = new TH1D("hst_dsig_dKEpi_cohnc_1000MeV","dsig/dKEpi, COH NC, Enu=1.0 GeV", nKEpi, KEpimin, KEpimax);
146  TH1D * hst_dsig_dKEpi_cohcc_1500MeV = new TH1D("hst_dsig_dKEpi_cohcc_1500MeV","dsig/dKEpi, COH CC, Enu=1.5 GeV", nKEpi, KEpimin, KEpimax);
147  TH1D * hst_dsig_dKEpi_cohnc_1500MeV = new TH1D("hst_dsig_dKEpi_cohnc_1500MeV","dsig/dKEpi, COH NC, Enu=1.5 GeV", nKEpi, KEpimin, KEpimax);
148 
149  //
150  // fill histograms
151  //
152  chain->Draw("(Ef-0.13957)>>hst_dsig_dKEpi_cohcc_0500MeV","coh&&cc&&Ev>0.49&&Ev<0.51&&pdgf==211","GOFF");
153  chain->Draw("(Ef-0.13498)>>hst_dsig_dKEpi_cohnc_0500MeV","coh&&nc&&Ev>0.49&&Ev<0.51&&pdgf==111","GOFF");
154  chain->Draw("(Ef-0.13957)>>hst_dsig_dKEpi_cohcc_1000MeV","coh&&cc&&Ev>0.99&&Ev<1.01&&pdgf==211","GOFF");
155  chain->Draw("(Ef-0.13498)>>hst_dsig_dKEpi_cohnc_1000MeV","coh&&nc&&Ev>0.99&&Ev<1.01&&pdgf==111","GOFF");
156  chain->Draw("(Ef-0.13957)>>hst_dsig_dKEpi_cohcc_1500MeV","coh&&cc&&Ev>1.49&&Ev<1.51&&pdgf==211","GOFF");
157  chain->Draw("(Ef-0.13498)>>hst_dsig_dKEpi_cohnc_1500MeV","coh&&nc&&Ev>1.49&&Ev<1.51&&pdgf==111","GOFF");
158 
159  //
160  // normalize
161  //
162  double norm_cohcc_0500MeV = hst_dsig_dKEpi_cohcc_0500MeV -> Integral("width") / sig_cohcc_0500MeV;
163  double norm_cohnc_0500MeV = hst_dsig_dKEpi_cohnc_0500MeV -> Integral("width") / sig_cohnc_0500MeV;
164  double norm_cohcc_1000MeV = hst_dsig_dKEpi_cohcc_1000MeV -> Integral("width") / sig_cohcc_1000MeV;
165  double norm_cohnc_1000MeV = hst_dsig_dKEpi_cohnc_1000MeV -> Integral("width") / sig_cohnc_1000MeV;
166  double norm_cohcc_1500MeV = hst_dsig_dKEpi_cohcc_1500MeV -> Integral("width") / sig_cohcc_1500MeV;
167  double norm_cohnc_1500MeV = hst_dsig_dKEpi_cohnc_1500MeV -> Integral("width") / sig_cohnc_1500MeV;
168 
169  if (norm_cohcc_0500MeV > 0) hst_dsig_dKEpi_cohcc_0500MeV -> Scale(1./norm_cohcc_0500MeV);
170  if (norm_cohnc_0500MeV > 0) hst_dsig_dKEpi_cohnc_0500MeV -> Scale(1./norm_cohnc_0500MeV);
171  if (norm_cohcc_1000MeV > 0) hst_dsig_dKEpi_cohcc_1000MeV -> Scale(1./norm_cohcc_1000MeV);
172  if (norm_cohnc_1000MeV > 0) hst_dsig_dKEpi_cohnc_1000MeV -> Scale(1./norm_cohnc_1000MeV);
173  if (norm_cohcc_1500MeV > 0) hst_dsig_dKEpi_cohcc_1500MeV -> Scale(1./norm_cohcc_1500MeV);
174  if (norm_cohnc_1500MeV > 0) hst_dsig_dKEpi_cohnc_1500MeV -> Scale(1./norm_cohnc_1500MeV);
175 
176  for(int i = 1; i <= hst_dsig_dKEpi_cohnc_1500MeV->GetNbinsX(); i++) {
177 
178  double Epi = hst_dsig_dKEpi_cohnc_1500MeV->GetBinCenter(i);
179 
180  double dsig_dKEpi_cohcc_0500MeV = hst_dsig_dKEpi_cohcc_0500MeV -> GetBinContent(i);
181  double dsig_dKEpi_cohnc_0500MeV = hst_dsig_dKEpi_cohnc_0500MeV -> GetBinContent(i);
182  double dsig_dKEpi_cohcc_1000MeV = hst_dsig_dKEpi_cohcc_1000MeV -> GetBinContent(i);
183  double dsig_dKEpi_cohnc_1000MeV = hst_dsig_dKEpi_cohnc_1000MeV -> GetBinContent(i);
184  double dsig_dKEpi_cohcc_1500MeV = hst_dsig_dKEpi_cohcc_1500MeV -> GetBinContent(i);
185  double dsig_dKEpi_cohnc_1500MeV = hst_dsig_dKEpi_cohnc_1500MeV -> GetBinContent(i);
186 
187  dsig_dKEpi_cohcc_0500MeV = TMath::Max(0., dsig_dKEpi_cohcc_0500MeV);
188  dsig_dKEpi_cohnc_0500MeV = TMath::Max(0., dsig_dKEpi_cohnc_0500MeV);
189  dsig_dKEpi_cohcc_1000MeV = TMath::Max(0., dsig_dKEpi_cohcc_1000MeV);
190  dsig_dKEpi_cohnc_1000MeV = TMath::Max(0., dsig_dKEpi_cohnc_1000MeV);
191  dsig_dKEpi_cohcc_1500MeV = TMath::Max(0., dsig_dKEpi_cohcc_1500MeV);
192  dsig_dKEpi_cohnc_1500MeV = TMath::Max(0., dsig_dKEpi_cohnc_1500MeV);
193 
194  out_stream << setw(15) << Epi
195  << setw(15) << dsig_dKEpi_cohcc_0500MeV
196  << setw(15) << dsig_dKEpi_cohnc_0500MeV
197  << setw(15) << dsig_dKEpi_cohcc_1000MeV
198  << setw(15) << dsig_dKEpi_cohnc_1000MeV
199  << setw(15) << dsig_dKEpi_cohcc_1500MeV
200  << setw(15) << dsig_dKEpi_cohnc_1500MeV
201  << endl;
202  }
203 
204  out_stream.close();
205 
206  // visual inspection
207 
208  TCanvas * c1 = new TCanvas("c1","",20,20,500,500);
209  hst_dsig_dKEpi_cohcc_1500MeV->Draw();
210  hst_dsig_dKEpi_cohnc_1500MeV->Draw("same");
211  c1->Update();
212 
213  delete chain;
214 }
const int kNRunsPerCase
Definition: nuint09_coh2.C:20
const int kNWCur
Definition: nuint09_coh2.C:18
double Integral(const Spectrum &s, const double pot, int cvnregion)
string filename
Definition: shutoffs.py:106
const char * label
const char * kLabel[kNSamples]
Definition: nuint09_coh2.C:23
chain
Check that an output directory exists.
OStream cout
Definition: OStream.cxx:6
const int kRunNu[kNSamples][kNWCur][kNEnergies][kNRunsPerCase]
Definition: nuint09_coh2.C:30
const int kNSamples
Definition: nuint09_coh2.C:17
const int kNEvtPerRun
Definition: nuint09_coh2.C:21
c1
Definition: demo5.py:24
const int kNEnergies
Definition: nuint09_coh2.C:19
simulatedPeak Scale(1/simulationNormalisationFactor)
void nuint09_coh2(int isample)
Definition: nuint09_coh2.C:60