make_sk_xsec_table.C
Go to the documentation of this file.
1 //________________________________________________________________________________________
2 /*!
3 
4 \macro make_sk_xsec_table.C
5 
6 \brief Write-out the cross-section table needed by SuperK softw for MC normalization
7 
8 \usage % genie 'make_sk_xsec_table.C("/some/path/genie_t2k_splines.xml")'
9 
10 \author Costas Andreopoulos <costas.andreopoulos@stfc.ac.uk>
11  University of Liverpool & STFC Rutherford Appleton Lab
12 
13  Ryan Terri <r.terri \at qmul.ac.uk>>
14  Queen Mary, University of London
15 
16 \created Nov 24, 2008
17 
18 \cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
19  For the full text of the license visit http://copyright.genie-mc.org
20  or see $GENIE/LICENSE
21 */
22 //_________________________________________________________________________________________
23 
24 #include <fstream>
25 #include <iomanip>
26 #include <cassert>
27 
28 #include "Conventions/Units.h"
29 
30 using std::cout;
31 using std::endl;
32 using std::setw;
33 using std::setprecision;
34 using std::setfill;
35 using std::ios;
36 
37 using namespace genie;
38 
40 (
41  string genie_cross_section_xml_file,
42  double EvMin = 0.025, // minimum energy in output file, GeV
43  double EvMax = 15.000, // maximum energy in output file, GeV
44  double dEv = 0.050 // energy binning, GeV
45 )
46 {
47  //
48  // Load cross-section splines
49  //
50 
52  XmlParserStatus_t ist = xspl->LoadFromXml(genie_cross_section_xml_file);
53  assert(ist == kXmlOK);
54 
55  //
56  // Create GENIE & configure event generation drivers for all init states
57  //
58 
59  cout << "Initializing event generation drivers!" << endl;
60 
61  GEVGDriver numu_H1;
62  GEVGDriver numu_O16;
63  GEVGDriver numubar_H1;
64  GEVGDriver numubar_O16;
65  GEVGDriver nue_H1;
66  GEVGDriver nue_O16;
67  GEVGDriver nuebar_H1;
68  GEVGDriver nuebar_O16;
77  numu_H1. UseSplines();
78  numu_O16. UseSplines();
79  numubar_H1. UseSplines();
80  numubar_O16.UseSplines();
81  nue_H1. UseSplines();
82  nue_O16. UseSplines();
83  nuebar_H1. UseSplines();
84  nuebar_O16. UseSplines();
85 
86  //
87  // Instruct the drivers to sum-up the cross-section splines for all simulated processes
88  // (compute total cross-section)
89  //
90 
91  cout << "Calculating total interaction cross-sections" << endl;
92 
93  int splNumKnots = 300; // number of knots in total cross section spline
94  double splEvMin = 0.010; // min energy in the spline
95  double splEvMax = 50.000; // max energy in the spline
96  bool splInLogE = true; // spline knots distributed logarithmicaly in energy
97  numu_H1. CreateXSecSumSpline (splNumKnots, splEvMin, splEvMax, splInLogE);
98  numu_O16. CreateXSecSumSpline (splNumKnots, splEvMin, splEvMax, splInLogE);
99  numubar_H1. CreateXSecSumSpline (splNumKnots, splEvMin, splEvMax, splInLogE);
100  numubar_O16.CreateXSecSumSpline (splNumKnots, splEvMin, splEvMax, splInLogE);
101  nue_H1. CreateXSecSumSpline (splNumKnots, splEvMin, splEvMax, splInLogE);
102  nue_O16. CreateXSecSumSpline (splNumKnots, splEvMin, splEvMax, splInLogE);
103  nuebar_H1. CreateXSecSumSpline (splNumKnots, splEvMin, splEvMax, splInLogE);
104  nuebar_O16. CreateXSecSumSpline (splNumKnots, splEvMin, splEvMax, splInLogE);
105 
106  //
107  // Write out the cross-section table in the format required by SKDETSIM
108  //
109 
110  cout << "Writing out the GENIE cross-section table" << endl;
111 
112  ofstream outf("./genie_sk_xsec_table.dat",ios::out);
113 
114  double w_H1 = 2; // # of H in H2O
115  double w_O16 = 1; // # of O in H2O
116 
117  double Ev = EvMin;
118  while(1) {
119  double xsec_numu_H1 = numu_H1. XSecSumSpline() -> Evaluate(Ev) / (1E-38 * units::cm2);
120  double xsec_numu_O16 = numu_O16. XSecSumSpline() -> Evaluate(Ev) / (1E-38 * units::cm2);
121  double xsec_numubar_H1 = numubar_H1. XSecSumSpline() -> Evaluate(Ev) / (1E-38 * units::cm2);
122  double xsec_numubar_O16 = numubar_O16. XSecSumSpline() -> Evaluate(Ev) / (1E-38 * units::cm2);
123  double xsec_nue_H1 = nue_H1. XSecSumSpline() -> Evaluate(Ev) / (1E-38 * units::cm2);
124  double xsec_nue_O16 = nue_O16. XSecSumSpline() -> Evaluate(Ev) / (1E-38 * units::cm2);
125  double xsec_nuebar_H1 = nuebar_H1. XSecSumSpline() -> Evaluate(Ev) / (1E-38 * units::cm2);
126  double xsec_nuebar_O16 = nuebar_O16. XSecSumSpline() -> Evaluate(Ev) / (1E-38 * units::cm2);
127 
128  double xsec_numu_H20 = w_H1 * xsec_numu_H1 + w_O16 * xsec_numu_O16 ;
129  double xsec_numubar_H20 = w_H1 * xsec_numubar_H1 + w_O16 * xsec_numubar_O16 ;
130  double xsec_nue_H20 = w_H1 * xsec_nue_H1 + w_O16 * xsec_nue_O16 ;
131  double xsec_nuebar_H20 = w_H1 * xsec_nuebar_H1 + w_O16 * xsec_nuebar_O16 ;
132 
133  outf << setfill(' ') << setw(10) << std::fixed << setprecision(5) << Ev
134  << setfill(' ') << setw(12) << std::fixed << setprecision(5) << xsec_numu_H20
135  << setfill(' ') << setw(12) << std::fixed << setprecision(5) << xsec_numubar_H20
136  << setfill(' ') << setw(12) << std::fixed << setprecision(5) << xsec_nue_H20
137  << setfill(' ') << setw(12) << std::fixed << setprecision(5) << xsec_nuebar_H20
138  << endl;
139 
140  Ev += dEv;
141  if(Ev > EvMax) break;
142  }
143 
144  outf.close();
145 
146  cout << "Done!" << endl;
147 }
void CreateXSecSumSpline(int nk, double Emin, double Emax, bool inlogE=true)
Definition: GEVGDriver.cxx:441
const int kPdgNuE
Definition: PDGCodes.h:28
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
const int kPdgAntiNuE
Definition: PDGCodes.h:29
const int kPdgNuMu
Definition: PDGCodes.h:30
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:61
static XSecSplineList * Instance()
static const double cm2
Definition: Units.h:77
const int kPdgTgtO16
Definition: PDGCodes.h:180
void make_sk_xsec_table(string genie_cross_section_xml_file, double EvMin=0.025, double EvMax=15.000, double dEv=0.050)
Float_t E
Definition: plot.C:20
const int kPdgTgtFreeP
Definition: PDGCodes.h:176
TFile * outf
Definition: testXsec.C:51
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Definition: GEVGDriver.h:55
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
OStream cout
Definition: OStream.cxx:6
void Configure(string mesg)
Definition: gEvServ.cxx:196
void Configure(int nu_pdgc, int Z, int A)
Definition: GEVGDriver.cxx:138
assert(nhit_max >=nhit_nbins)
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:53
void UseSplines(void)
Definition: GEVGDriver.cxx:509
enum genie::EXmlParseStatus XmlParserStatus_t
List of cross section vs energy splines.
XmlParserStatus_t LoadFromXml(const string &filename, bool keep=false)