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-2017, 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 
39 void make_sk_xsec_table
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 
51  XSecSplineList * xspl = XSecSplineList::Instance();
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;
69  numu_H1. Configure (kPdgNuMu, pdg::IonPdgCodeToZ(kPdgTgtFreeP), pdg::IonPdgCodeToA(kPdgTgtFreeP));
70  numu_O16. Configure (kPdgNuMu, pdg::IonPdgCodeToZ(kPdgTgtO16), pdg::IonPdgCodeToA(kPdgTgtO16) );
71  numubar_H1. Configure (kPdgAntiNuMu, pdg::IonPdgCodeToZ(kPdgTgtFreeP), pdg::IonPdgCodeToA(kPdgTgtFreeP));
72  numubar_O16.Configure (kPdgAntiNuMu, pdg::IonPdgCodeToZ(kPdgTgtO16), pdg::IonPdgCodeToA(kPdgTgtO16) );
73  nue_H1. Configure (kPdgNuE, pdg::IonPdgCodeToZ(kPdgTgtFreeP), pdg::IonPdgCodeToA(kPdgTgtFreeP));
74  nue_O16. Configure (kPdgNuE, pdg::IonPdgCodeToZ(kPdgTgtO16), pdg::IonPdgCodeToA(kPdgTgtO16) );
75  nuebar_H1. Configure (kPdgAntiNuE, pdg::IonPdgCodeToZ(kPdgTgtFreeP), pdg::IonPdgCodeToA(kPdgTgtFreeP));
76  nuebar_O16. Configure (kPdgAntiNuE, pdg::IonPdgCodeToZ(kPdgTgtO16), pdg::IonPdgCodeToA(kPdgTgtO16) );
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 }
genie::GEVGDriver::CreateXSecSumSpline
void CreateXSecSumSpline(int nk, double Emin, double Emax, bool inlogE=true)
Definition: GEVGDriver.cxx:476
genie::kPdgNuE
const int kPdgNuE
Definition: PDGCodes.h:25
genie
include "Numerical/GSFunc.h"
Definition: AlgCmp.h:26
genie::kPdgAntiNuE
const int kPdgAntiNuE
Definition: PDGCodes.h:26
genie::kPdgNuMu
const int kPdgNuMu
Definition: PDGCodes.h:27
genie::pdg::IonPdgCodeToA
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:74
genie::XSecSplineList::Instance
static XSecSplineList * Instance()
Definition: XSecSplineList.cxx:103
genie::units::cm2
static const double cm2
Definition: Units.h:77
cout
the ParameterSet object passed in for the configuration of a destination should be the only source that can affect the behavior of that destination This is to eliminate the dependencies of configuring a destination from multiple mostly from the defaults It suppresses possible glitches about changing the configuration file somewhere outside of a destination segament might still affect the behavior of that destination In the previous configuration for a specific the value of a certain e may come from following and have been suppressed It the configuring ParameterSet object for each destination will be required to carry a parameter list as complete as possible If a parameter still cannot be found in the ParameterSet the configuration code will go look for a hardwired default directly The model is a great simplicity comparing with the previous especially when looking for default values Another great advantage is most of the parameters now have very limited places that allows to appear Usually they can only appear at one certain level in a configuration file For in the old configuring model or in a default ParameterSet object inside of a or in a category or in a severity object This layout of multiple sources for a single parameter brings great confusion in both writing a configuration and in processing the configuration file Under the new the only allowed place for the parameter limit to appear is inside of a category which is inside of a destination object Other improvements simplify the meaning of a destination name In the old a destination name has multiple folds of meanings the e cout and cerr have the special meaning of logging messages to standard output or standard error the name also serves as the output filename if the destination is a file these names are also references to look up for detailed configurations in configuring the MessageFacility The multi purpose of the destination name might cause some unwanted behavior in either writing or parsing the configuration file To amend in the new model the destination name is now merely a name for a which might represent the literal purpose of this or just an id All other meanings of the destinations names now go into the destination ParameterSet as individual such as the type parameter and filename parameter Following is the deatiled rule for the new configuring Everything that is related with MessageFacility configuration must be wrapped in a single ParameterSet object with the name MessageFacility The MessageFacility ParameterSet object contains a series of top level parameters These parameters can be chosen a vector of string listing the name of debug enabled models Or use *to enable debug messages in all modules a vector of string a vector of string a vector of string a ParameterSet object containing the list of all destinations The destinations ParameterSet object is a combination of ParameterSet objects for individual destinations There are two types of destinations that you can insert in the destinations ParameterSet ordinary including cout
Definition: MessageFacilityAsStandAlonePackage.txt:82
genie::kPdgTgtO16
const int kPdgTgtO16
Definition: PDGCodes.h:175
genie::XSecSplineList::LoadFromXml
XmlParserStatus_t LoadFromXml(string filename, bool keep=false)
Definition: XSecSplineList.cxx:293
E
const Double_t E[nknots]
Definition: testXsec.C:25
make_sk_xsec_table
void make_sk_xsec_table(string genie_cross_section_xml_file, double EvMin=0.025, double EvMax=15.000, double dEv=0.050)
Definition: make_sk_xsec_table.C:40
genie::kPdgTgtFreeP
const int kPdgTgtFreeP
Definition: PDGCodes.h:171
outf
TFile * outf
Definition: testXsec.C:51
genie::GEVGDriver
GENIE Event Generation Driver. A minimalist user interface object for generating neutrino interaction...
Definition: GEVGDriver.h:55
genie::kPdgAntiNuMu
const int kPdgAntiNuMu
Definition: PDGCodes.h:28
Configure
void Configure(string mesg)
Definition: gEvServ.cxx:196
genie::GEVGDriver::Configure
void Configure(int nu_pdgc, int Z, int A)
Definition: GEVGDriver.cxx:173
genie::pdg::IonPdgCodeToZ
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:66
genie::GEVGDriver::UseSplines
void UseSplines(void)
Definition: GEVGDriver.cxx:544
genie::XmlParserStatus_t
enum genie::EXmlParseStatus XmlParserStatus_t
genie::XSecSplineList
List of cross section vs energy splines.
Definition: XSecSplineList.h:46
Generated by   doxygen 1.8.11