sk_sample_norm_abs.C
Go to the documentation of this file.
1 //________________________________________________________________________________________
2 /*!
3 
4 \macro sk_sample_norm_abs.C
5 
6 \brief Calculate absolute normalization of SuperK event samples.
7 
8 \details The program calculates the number of events of each flavour per 1E+21 POT
9  per 22.5 kton water fiducial.
10 
11  N = Integral{
12  (d3Flux / dE dS dI) * sig(E) * (Na/A) * rho*L * dE*dS*dI}
13  = Integral{
14  (d3Flux / dE dS dI) * sig(E) * (Na/A) * dE*dM*dI} (1)
15 
16  where
17  - (d3Flux / dE dS dI): numu flux per energy bin, per unit area, per POT
18  - sigma(E): total numu cross section on water
19  - Na: Avogadro's number
20  - A: mass number for water
21  - rho: water density
22  - L: path length
23  - M: mass
24  - I: beam intensity (POT)
25 
26  The SuperK flux is given in #neutrinos per 0.05 GeV per cm2 per 1E+21 POTs.
27  Input water cross sections are given in 1E-38 cm2.
28  Equation (1) becomes
29 
30  N = 6.023E+23 x 1E-38 x (Mfv/A) x Ebinsz x Sum_{i} { F_{i} * sig_{i} }
31 
32  where
33  - N : expected number of events for NF x 1E+21 POT
34  - Mfv: fiducial volume mass
35  - NF : number of 1E+21 POT worth of flux files chained together to produce
36  the flux histograms
37  - Ebinsz: energy bin size (0.05 GeV)
38  - F_{i): flux in bin i (in number of flux neutrinos / Ebinsz / (NF x 1E+21 POT))
39  - sig_{i}: cross section evaluated at centre of bin i (1E-38 cm^2)
40 
41  Notes:
42  Ptot = P(H1) + P(O16) = sigma(H1)*w(H1)*rho/A(H1) + sigma(O16)*w(O16)*rho/A(O16)
43  rho: water density and w: mass contribution, w(H1)=2/18, w(O16)=16/18
44  So: Ptot ~ sigma(H1) * (2/18) + sigma(O16) * (16/18) / 16 =>
45  Ptot ~ (2*sigma(H1)+sigma(O16))/18 =>
46  Ptot ~ sigma(H20)/A(H20)
47 
48 \inputs
49  - xsecfile:
50  Neutrino - water cross section file.
51  The output of $GENIE/src/contrib/t2k/make_sk_xsec_table.C
52  - skfluxfile :
53  Input neutrino flux file.
54  The output of $GENIE/src/scripts/production/misc/generate_sk_flux_histograms.C
55  - IF :
56  Number of POTs per flux simulation file used for filling-in the flux histograms
57  (typically 1E+21 POT)
58  - NF :
59  Number of flux files used filling-in the flux histograms
60 
61 
62 \author Costas Andreopoulos <costas.andreopoulos@stfc.ac.uk>
63  University of Liverpool & STFC Rutherford Appleton Lab
64 
65 \created Nov 24, 2008
66 
67 \cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
68  For the full text of the license visit http://copyright.genie-mc.org
69  or see $GENIE/LICENSE
70 */
71 //_________________________________________________________________________________________
72 
73 
75 (
76  const char * xsecfile, // Neutrino - water cross section file
77  const char * skfluxfile, // Input neutrino flux file.
78  double IF = 1E+21, // Number of POTs per flux file
79  int NF = 500 // Number of flux files used
80 )
81 
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 }
const int Na
double dE
Loaders::FluxType flux
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.
Float_t E
Definition: plot.C:20
OStream cout
Definition: OStream.cxx:6
static const double A
Definition: Units.h:82
TTree xsec_water
Definition: FillPIDs.h:16
double Emax
int nE