Functions
sk_sample_norm_abs.C File Reference

Go to the source code of this file.

Functions

void sk_sample_norm_abs (const char *xsecfile, const char *skfluxfile, double IF=1E+21, int NF=500)
 Calculate absolute normalization of SuperK event samples. More...
 

Function Documentation

void sk_sample_norm_abs ( const char *  xsecfile,
const char *  skfluxfile,
double  IF = 1E+21,
int  NF = 500 
)

Calculate absolute normalization of SuperK event samples.

sk_sample_norm_abs.C

The program calculates the number of events of each flavour per 1E+21 POT per 22.5 kton water fiducial.

N = Integral{ (d3Flux / dE dS dI) * sig(E) * (Na/A) * rho*L * dE*dS*dI} = Integral{ (d3Flux / dE dS dI) * sig(E) * (Na/A) * dE*dM*dI} (1)

where

  • (d3Flux / dE dS dI): numu flux per energy bin, per unit area, per POT
  • sigma(E): total numu cross section on water
  • Na: Avogadro's number
  • A: mass number for water
  • rho: water density
  • L: path length
  • M: mass
  • I: beam intensity (POT)

The SuperK flux is given in #neutrinos per 0.05 GeV per cm2 per 1E+21 POTs. Input water cross sections are given in 1E-38 cm2. Equation (1) becomes

N = 6.023E+23 x 1E-38 x (Mfv/A) x Ebinsz x Sum_{i} { F_{i} * sig_{i} }

where

  • N : expected number of events for NF x 1E+21 POT
  • Mfv: fiducial volume mass
  • NF : number of 1E+21 POT worth of flux files chained together to produce the flux histograms
  • Ebinsz: energy bin size (0.05 GeV)
  • F_{i): flux in bin i (in number of flux neutrinos / Ebinsz / (NF x 1E+21 POT))
  • sig_{i}: cross section evaluated at centre of bin i (1E-38 cm^2)

Notes: Ptot = P(H1) + P(O16) = sigma(H1)*w(H1)*rho/A(H1) + sigma(O16)*w(O16)*rho/A(O16) rho: water density and w: mass contribution, w(H1)=2/18, w(O16)=16/18 So: Ptot ~ sigma(H1) * (2/18) + sigma(O16) * (16/18) / 16 => Ptot ~ (2*sigma(H1)+sigma(O16))/18 => Ptot ~ sigma(H20)/A(H20)

  • xsecfile: Neutrino - water cross section file. The output of $GENIE/src/contrib/t2k/make_sk_xsec_table.C
  • skfluxfile : Input neutrino flux file. The output of $GENIE/src/scripts/production/misc/generate_sk_flux_histograms.C
  • IF : Number of POTs per flux simulation file used for filling-in the flux histograms (typically 1E+21 POT)
  • NF : Number of flux files used filling-in the flux histograms
Author
Costas Andreopoulos costa.nosp@m.s.an.nosp@m.dreop.nosp@m.oulo.nosp@m.s@stf.nosp@m.c.ac.nosp@m..uk University of Liverpool & STFC Rutherford Appleton Lab

Nov 24, 2008

Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 75 of file sk_sample_norm_abs.C.

82 {
83  // consts
84  double Mfv = 2.25E+10; // want final event numbers shown for this ficucial mass (gr)
85  double I0 = 1E+21; // want final event numbers shown for this POT exposure
86  double Na = 6.023E+23;
87  double A = 18; // gr
88 
89  // binning in xsection and flux files
90 
91  double Emin = 0.0;
92  double Emax = 15.0;
93  double dE = 0.050;
94  int nE = (Emax-Emin)/dE;
95 
96  // load genie cross sections for water
97 
98  TTree xsec_water;
99  xsec_water.ReadFile(xsecfile, "E/D:xsec_numu/D:xsec_numubar/D:xsec_nue/D:xsec_nuebar/D");
100  TH1D * xnumu = new TH1D("xnumu", "", nE, Emin, Emax);
101  TH1D * xnumubar = new TH1D("xnumubar", "", nE, Emin, Emax);
102  TH1D * xnue = new TH1D("xnue", "", nE, Emin, Emax);
103  TH1D * xnuebar = new TH1D("xnuebar", "", nE, Emin, Emax);
104  xsec_water.Draw("E>>xnumu", "xsec_numu", "goff");
105  xsec_water.Draw("E>>xnumubar", "xsec_numubar", "goff");
106  xsec_water.Draw("E>>xnue", "xsec_nue", "goff");
107  xsec_water.Draw("E>>xnuebar", "xsec_nuebar", "goff");
108 
109  // load SuperK flux histograms
110 
111  TFile * flux = new TFile(skfluxfile, "read");
112  TH1D * fnumu = (TH1D*) flux->Get("numu_flux");
113  TH1D * fnumubar = (TH1D*) flux->Get("numubar_flux");
114  TH1D * fnue = (TH1D*) flux->Get("nue_flux"); // instrinsic beam nue
115  TH1D * fnuebar = (TH1D*) flux->Get("nuebar_flux"); // ...
116 
117  // integrate flux x cross section and apply exposure / fiducial mass and dimensional factors
118 
119  double Nnumu = 0;
120  double Nnumubar = 0;
121  double Nnue = 0;
122  double Nnuebar = 0;
123  double Nnuesig = 0; // numu->nue
124 
125  for(int i=1; i<=fnumu->GetNbinsX(); i++) {
126  double E = fnumu->GetBinCenter(i);
127 
128  double xsec_numu = xnumu -> GetBinContent (xnumu -> FindBin(E));
129  double flux_numu = fnumu -> GetBinContent (fnumu -> FindBin(E));
130  double xsec_numubar = xnumubar -> GetBinContent (xnumubar -> FindBin(E));
131  double flux_numubar = fnumubar -> GetBinContent (fnumubar -> FindBin(E));
132  double xsec_nue = xnue -> GetBinContent (xnue -> FindBin(E));
133  double flux_nue = fnue -> GetBinContent (fnue -> FindBin(E));
134  double xsec_nuebar = xnuebar -> GetBinContent (xnuebar -> FindBin(E));
135  double flux_nuebar = fnuebar -> GetBinContent (fnuebar -> FindBin(E));
136 
137  cout << "E = " << E << " GeV";
138  cout << " - numu : sigma(H20) = " << xsec_numu << " x1E-38 cm2, flux(@SK) = " << flux_numu
139  << " /" << dE << " GeV /" << (NF*IF) << " POT /cm2" << endl;
140  cout << " - numubar: sigma(H20) = " << xsec_numubar << " x1E-38 cm2, flux(@SK) = " << flux_numubar
141  << " /" << dE << " GeV /" << (NF*IF) << " POT /cm2" << endl;
142  cout << " - nue : sigma(H20) = " << xsec_nue << " x1E-38 cm2, flux(@SK) = " << flux_nue
143  << " /" << dE << " GeV /" << (NF*IF) << " POT /cm2" << endl;
144  cout << " - nuebar : sigma(H20) = " << xsec_nuebar << " x1E-38 cm2, flux(@SK) = " << flux_nuebar
145  << " /" << dE << " GeV /" << (NF*IF) << " POT /cm2" << endl;
146 
147  Nnumu += ( flux_numu * xsec_numu );
148  Nnumubar += ( flux_numubar * xsec_numubar );
149  Nnue += ( flux_nue * xsec_nue );
150  Nnuebar += ( flux_nuebar * xsec_nuebar );
151  Nnuesig += ( flux_numu * xsec_nue ); // 100% numu->nue
152  }
153 
154  double f = Na * 1E-38 * (Mfv/A) * I0 / (NF*IF);
155 
156  Nnumu *= f;
157  Nnumubar *= f;
158  Nnue *= f;
159  Nnuebar *= f;
160  Nnuesig *= f;
161 
162  // print-out results
163 
164  cout << endl;
165  cout << endl;
166  cout << "---------------------------------------------------------------------------" << endl;
167  cout << "species | # of SK events per "
168  << I0 << " POT per " << Mfv/1E+9 << " kton fiducial"<< endl;
169  cout << "---------------------------------------------------------------------------" << endl;
170  cout << "numu | " << Nnumu << endl;
171  cout << "numubar | " << Nnumubar << endl;
172  cout << "nue(bkg) | " << Nnue << endl;
173  cout << "nuebar(bkg) | " << Nnuebar << endl;
174  cout << "nue(sig,100% numu->nue) | " << Nnuesig << endl;
175  cout << "---------------------------------------------------------------------------" << endl;
176  cout << endl;
177 
178 }
i
size_t i(0)
f
FILE * f
Definition: loadlibs.C:26
dE
double dE
Definition: write_out_h20_xsec.C:9
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::units::A
static const double A
Definition: Units.h:82
xsec_water
TTree xsec_water
Definition: write_out_h20_xsec.C:19
Emax
double Emax
Definition: write_out_h20_xsec.C:8
nE
int nE
Definition: write_out_h20_xsec.C:10
Generated by   doxygen 1.8.11