eventRatesL.cc
Go to the documentation of this file.
1 #define eventRates_cxx
2 #include "eventRates.h"
3 #include "OscLib/OscCalculator.cxx"
4 #include <TH2.h>
5 #include <TH1D.h>
6 #include <TStyle.h>
7 #include <TCanvas.h>
8 #include <iomanip>
9 #include <iostream>
10 #include <fstream>
11 #include <sstream>
12 #include <stdio.h>
13 #include <math.h>
14 #include <assert.h>
15 #include <TGaxis.h>
16 #include <TParameter.h>
17 
18 
19 #ifndef __CINT__ // hide header stuff from CINT, assume load_dk2nu.C run first
20 #include "dk2nu/tree/dk2nu.h"
21 #include "dk2nu/tree/dkmeta.h"
22 #endif // ifndef __CINT__
23 
24 
25 using namespace std;
26 
27 int iread = 0;
28 
29 void eventRates::Loop()
30 {
31 
32  // In a ROOT session, you can do:
33  // Root > .L eventRates.C++ // The ++ is important because the code runs MUCH faster when compiled
34  // Root > eventRates t
35  // Root > t.Loop(); // Loop on all entries
36  //
37 
38  //
39  // Define the Reference pot
40  //
41 
42  double refpot = 1;
43  std::string potref_str = eventRates::GetPOTAsString(refpot);
44 
45  //
46  // Set Histogram Binning
47  //
48 
49  // simple binning for some alignment plots
50  int nbins = 40;
51  double xmin = 0.0;
52  double xmax = 20.0;
53 
54  // more complex binning for passing to fastMC
55  std::vector< Double_t > fastmc_bins;
56  // 0.125 GeV bins up to 8 GeV
57  for(int i = 0; i<8/0.125; i++)
58  fastmc_bins.push_back(i*.125);
59  // 0.5 GeV bins up to 20 GeV
60  for(int i = 0; i<(20-8)/0.5; i++)
61  fastmc_bins.push_back(8.0+i*.5);
62  // 2.0 GeV bins up to 120 GeV
63  for(int i = 0; i<(120-20)/2.0; i++)
64  fastmc_bins.push_back(20.0+i*2.0);
65 
66  //
67  // Declare all the histograms
68  //
69 
70  TH1D *fhNuMuFlux = new TH1D("numu_flux_forplots",
71  "numu_flux_forplots", nbins,xmin,xmax);
72  TH1D *fhNuMuBarFlux = new TH1D("numubar_flux_forplots",
73  "numubar_flux_forplots", nbins/2,xmin,xmax);
74  TH1D *fhNuEFlux = new TH1D("nue_flux_forplots",
75  "nue_flux_forplots", nbins/2,xmin,xmax);
76  TH1D *fhNuEBarFlux = new TH1D("nuebar_flux_forplots",
77  "nuebar_flux_forplots", nbins/2,xmin,xmax);
78 
79  TH1D *fhNuTauFlux = new TH1D("nutau_flux_forplots",
80  "nutau_flux_forplots", nbins/2,xmin,xmax);
81  TH1D *fhNuTauBarFlux = new TH1D("nutaubar_flux_forplots",
82  "nutaubar_flux_forplots", nbins/2,xmin,xmax);
83 
84  TH1D *fhNuMuCCEventRate = new TH1D("numu_cceventrate_forplots",
85  "numu_cceventrate_forplots", nbins,xmin,xmax);
86  TH1D *fhNuMuBarCCEventRate = new TH1D("numubar_cceventrate_forplots",
87  "numubar_cceventrate_forplots", nbins/2,xmin,xmax);
88  TH1D *fhNuECCEventRate = new TH1D("nue_cceventrate_forplots",
89  "nue_cceventrate_forplots", nbins/2,xmin,xmax);
90  TH1D *fhNuEBarCCEventRate = new TH1D("nuebar_cceventrate_forplots",
91  "nuebar_cceventrate_forplots", nbins/2,xmin,xmax);
92 
93  TH1D *fhNuTauCCEventRate = new TH1D("nutau_cceventrate_forplots",
94  "nutau_cceventrate_forplots", nbins/2,xmin,xmax);
95  TH1D *fhNuTauBarCCEventRate = new TH1D("nutaubar_cceventrate_forplots",
96  "nutaubar_cceventrate_forplots", nbins/2,xmin,xmax);
97 
98  TH1D *fhNuMuNCEventRate = new TH1D("numu_nceventrate_forplots",
99  "numu_nceventrate_forplots", nbins,xmin,xmax);
100  TH1D *fhNuMuBarNCEventRate = new TH1D("numubar_nceventrate_forplots",
101  "numubar_nceventrate_forplots", nbins/2,xmin,xmax);
102  TH1D *fhNuENCEventRate = new TH1D("nue_nceventrate_forplots",
103  "nue_nceventrate_forplots", nbins/2,xmin,xmax);
104  TH1D *fhNuEBarNCEventRate = new TH1D("nuebar_nceventrate_forplots",
105  "nuebar_nceventrate_forplots", nbins/2,xmin,xmax);
106 
107  TH1D *fhNuTauNCEventRate = new TH1D("nutau_nceventrate_forplots",
108  "nutau_nceventrate_forplots", nbins/2,xmin,xmax);
109  TH1D *fhNuTauBarNCEventRate = new TH1D("nutaubar_nceventrate_forplots",
110  "nutaubar_nceventrate_forplots", nbins/2,xmin,xmax);
111 
112 
113  TH1D *fhNuMuFlux_FastMC = new TH1D("numu_flux",
114  "numu_flux",fastmc_bins.size()-1,&fastmc_bins[0]);
115  TH1D *fhNuMuBarFlux_FastMC = new TH1D("numubar_flux",
116  "numubar_flux", fastmc_bins.size()-1,&fastmc_bins[0]);
117  TH1D *fhNuEFlux_FastMC = new TH1D("nue_flux",
118  "nue_flux", fastmc_bins.size()-1,&fastmc_bins[0]);
119  TH1D *fhNuEBarFlux_FastMC = new TH1D("nuebar_flux",
120  "nuebar_flux", fastmc_bins.size()-1,&fastmc_bins[0]);
121  TH1D *fhNuTauFlux_FastMC = new TH1D("nutau_flux",
122  "nutau_flux", fastmc_bins.size()-1,&fastmc_bins[0]);
123  TH1D *fhNuTauBarFlux_FastMC = new TH1D("nutaubar_flux",
124  "nutaubar_flux", fastmc_bins.size()-1,&fastmc_bins[0]);
125 
126  TH1D *fhNuMuCCEventRate_FastMC = new TH1D("numu_cceventrate",
127  "numu_cceventrate",fastmc_bins.size()-1,&fastmc_bins[0]);
128  TH1D *fhNuMuBarCCEventRate_FastMC = new TH1D("numubar_cceventrate",
129  "numubar_cceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
130  TH1D *fhNuECCEventRate_FastMC = new TH1D("nue_cceventrate",
131  "nue_cceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
132  TH1D *fhNuEBarCCEventRate_FastMC = new TH1D("nuebar_cceventrate",
133  "nuebar_cceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
134  TH1D *fhNuTauCCEventRate_FastMC = new TH1D("nutau_cceventrate",
135  "nutau_cceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
136  TH1D *fhNuTauBarCCEventRate_FastMC = new TH1D("nutaubar_cceventrate",
137  "nutaubar_cceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
138 
139  TH1D *fhNuMuNCEventRate_FastMC = new TH1D("numu_nceventrate",
140  "numu_nceventrate",fastmc_bins.size()-1,&fastmc_bins[0]);
141  TH1D *fhNuMuBarNCEventRate_FastMC = new TH1D("numubar_nceventrate",
142  "numubar_nceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
143  TH1D *fhNuENCEventRate_FastMC = new TH1D("nue_nceventrate",
144  "nue_nceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
145  TH1D *fhNuEBarNCEventRate_FastMC = new TH1D("nuebar_nceventrate",
146  "nuebar_nceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
147  TH1D *fhNuTauNCEventRate_FastMC = new TH1D("nutau_nceventrate",
148  "nutau_nceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
149  TH1D *fhNuTauBarNCEventRate_FastMC = new TH1D("nutaubar_nceventrate",
150  "nutaubar_nceventrate", fastmc_bins.size()-1,&fastmc_bins[0]);
151 
152  double globes_low = 0;
153  double globes_high = 125.25;
154  double globes_binwidth = 0.25;
155  int globes_nbins = ( globes_high - globes_low ) / globes_binwidth;
156 
157  TH1D *fhNuMuFlux_Globes = new TH1D("numu_flux_globes",
158  "numu_flux_globes", globes_nbins,globes_low,globes_high);
159  TH1D *fhNuMuBarFlux_Globes = new TH1D("numubar_flux_globes",
160  "numubar_flux_globes", globes_nbins, globes_low, globes_high);
161  TH1D *fhNuEFlux_Globes = new TH1D("nue_flux_globes",
162  "nue_flux_globes", globes_nbins, globes_low, globes_high);
163  TH1D *fhNuEBarFlux_Globes = new TH1D("nuebar_flux_globes",
164  "nuebar_flux_globes", globes_nbins, globes_low, globes_high);
165 
166  TH1D *fhNuTauFlux_Globes = new TH1D("nutau_flux_globes",
167  "nutau_flux_globes", globes_nbins, globes_low, globes_high);
168  TH1D *fhNuTauBarFlux_Globes = new TH1D("nutaubar_flux_globes",
169  "nutau_flux_globes", globes_nbins, globes_low, globes_high);
170 
171  TH1D *fhNuMuFluxOsc = new TH1D("numu_fluxosc_forplots",
172  "numu_fluxosc_forplots", nbins,xmin,xmax);
173  TH1D *fhNuMuBarFluxOsc = new TH1D("numubar_fluxosc_forplots",
174  "numubar_fluxosc_forplots", nbins/2,xmin,xmax);
175  TH1D *fhNuEFluxOsc = new TH1D("nue_fluxosc_forplots",
176  "nue_fluxosc_forplots", nbins/2,xmin,xmax);
177  TH1D *fhNuEBarFluxOsc = new TH1D("nuebar_fluxosc_forplots",
178  "nuebar_fluxosc_forplots", nbins/2,xmin,xmax);
179 
180  TH1D *fhNuTauFluxOsc = new TH1D("nutau_fluxosc_forplots",
181  "nutau_fluxosc_forplots", nbins/2,xmin,xmax);
182  TH1D *fhNuTauBarFluxOsc = new TH1D("nutaubar_fluxosc_forplots",
183  "nutaubar_fluxosc_forplots", nbins/2,xmin,xmax);
184 
185  TH1D *fhNuMuCCEventRateOsc = new TH1D("numu_cceventrateosc_forplots",
186  "numu_cceventrateosc_forplots", nbins,xmin,xmax);
187  TH1D *fhNuMuBarCCEventRateOsc = new TH1D("numubar_cceventrateosc_forplots",
188  "numubar_cceventrateosc_forplots", nbins/2,xmin,xmax);
189  TH1D *fhNuECCEventRateOsc = new TH1D("nue_cceventrateosc_forplots",
190  "nue_cceventrateosc_forplots", nbins/2,xmin,xmax);
191  TH1D *fhNuEBarCCEventRateOsc = new TH1D("nuebar_cceventrateosc_forplots",
192  "nuebar_cceventrateosc_forplots", nbins/2,xmin,xmax);
193 
194  TH1D *fhNuTauCCEventRateOsc = new TH1D("nutau_cceventrateosc_forplots",
195  "nutau_cceventrateosc_forplots", nbins/2,xmin,xmax);
196  TH1D *fhNuTauBarCCEventRateOsc = new TH1D("nutaubar_cceventrateosc_forplots",
197  "nutaubar_cceventrateosc_forplots", nbins/2,xmin,xmax);
198 
199  TH1D *fhNuMuNCEventRateOsc = new TH1D("numu_nceventrateosc_forplots",
200  "numu_nceventrateosc_forplots", nbins,xmin,xmax);
201  TH1D *fhNuMuBarNCEventRateOsc = new TH1D("numubar_nceventrateosc_forplots",
202  "numubar_nceventrateosc_forplots", nbins/2,xmin,xmax);
203  TH1D *fhNuENCEventRateOsc = new TH1D("nue_nceventrateosc_forplots",
204  "nue_nceventrateosc_forplots", nbins/2,xmin,xmax);
205  TH1D *fhNuEBarNCEventRateOsc = new TH1D("nuebar_nceventrateosc_forplots",
206  "nuebar_nceventrateosc_forplots", nbins/2,xmin,xmax);
207 
208  TH1D *fhNuTauNCEventRateOsc = new TH1D("nutau_nceventrateosc_forplots",
209  "nutau_nceventrateosc_forplots", nbins/2,xmin,xmax);
210  TH1D *fhNuTauBarNCEventRateOsc = new TH1D("nutaubar_nceventrateosc_forplots",
211  "nutaubar_nceventrateosc_forplots", nbins/2,xmin,xmax);
212 
213  TH1D *fhNuMuFluxOsc_FastMC = new TH1D("numu_fluxosc",
214  "numu_fluxosc",fastmc_bins.size()-1,&fastmc_bins[0]);
215  TH1D *fhNuMuBarFluxOsc_FastMC = new TH1D("numubar_fluxosc",
216  "numubar_fluxosc", fastmc_bins.size()-1,&fastmc_bins[0]);
217  TH1D *fhNuEFluxOsc_FastMC = new TH1D("nue_fluxosc",
218  "nue_fluxosc", fastmc_bins.size()-1,&fastmc_bins[0]);
219  TH1D *fhNuEBarFluxOsc_FastMC = new TH1D("nuebar_fluxosc",
220  "nuebar_fluxosc", fastmc_bins.size()-1,&fastmc_bins[0]);
221  TH1D *fhNuTauFluxOsc_FastMC = new TH1D("nutau_fluxosc",
222  "nutau_fluxosc", fastmc_bins.size()-1,&fastmc_bins[0]);
223  TH1D *fhNuTauBarFluxOsc_FastMC = new TH1D("nutaubar_fluxosc",
224  "nutaubar_fluxosc", fastmc_bins.size()-1,&fastmc_bins[0]);
225 
226  TH1D *fhNuMuCCEventRateOsc_FastMC = new TH1D("numu_cceventrateosc",
227  "numu_cceventrateosc",fastmc_bins.size()-1,&fastmc_bins[0]);
228  TH1D *fhNuMuBarCCEventRateOsc_FastMC = new TH1D("numubar_cceventrateosc",
229  "numubar_cceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
230  TH1D *fhNuECCEventRateOsc_FastMC = new TH1D("nue_cceventrateosc",
231  "nue_cceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
232  TH1D *fhNuEBarCCEventRateOsc_FastMC = new TH1D("nuebar_cceventrateosc",
233  "nuebar_cceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
234  TH1D *fhNuTauCCEventRateOsc_FastMC = new TH1D("nutau_cceventrateosc",
235  "nutau_cceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
236  TH1D *fhNuTauBarCCEventRateOsc_FastMC = new TH1D("nutaubar_cceventrateosc",
237  "nutaubar_cceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
238 
239  TH1D *fhNuMuNCEventRateOsc_FastMC = new TH1D("numu_nceventrateosc",
240  "numu_nceventrateosc",fastmc_bins.size()-1,&fastmc_bins[0]);
241  TH1D *fhNuMuBarNCEventRateOsc_FastMC = new TH1D("numubar_nceventrateosc",
242  "numubar_nceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
243  TH1D *fhNuENCEventRateOsc_FastMC = new TH1D("nue_nceventrateosc",
244  "nue_nceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
245  TH1D *fhNuEBarNCEventRateOsc_FastMC = new TH1D("nuebar_nceventrateosc",
246  "nuebar_nceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
247  TH1D *fhNuTauNCEventRateOsc_FastMC = new TH1D("nutau_nceventrateosc",
248  "nutau_nceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
249  TH1D *fhNuTauBarNCEventRateOsc_FastMC = new TH1D("nutaubar_nceventrateosc",
250  "nutaubar_nceventrateosc", fastmc_bins.size()-1,&fastmc_bins[0]);
251 
252  //
253  // Call Sumw2 to make sure histogram errors are propagated correctly
254  //
255 
256  fhNuMuFlux->Sumw2();
257  fhNuMuBarFlux->Sumw2();
258  fhNuEFlux->Sumw2();
259  fhNuEBarFlux->Sumw2();
260  fhNuTauFlux->Sumw2();
261  fhNuTauBarFlux->Sumw2();
262 
263  fhNuMuFlux_FastMC->Sumw2();
264  fhNuMuBarFlux_FastMC->Sumw2();
265  fhNuEFlux_FastMC->Sumw2();
266  fhNuEBarFlux_FastMC->Sumw2();
267  fhNuTauFlux_FastMC->Sumw2();
268  fhNuTauBarFlux_FastMC->Sumw2();
269 
270  fhNuMuFlux_Globes->Sumw2();
271  fhNuMuBarFlux_Globes->Sumw2();
272  fhNuEFlux_Globes->Sumw2();
273  fhNuEBarFlux_Globes->Sumw2();
274  fhNuTauFlux_Globes->Sumw2();
275  fhNuTauBarFlux_Globes->Sumw2();
276 
277  fhNuMuFluxOsc->Sumw2();
278  fhNuMuBarFluxOsc->Sumw2();
279  fhNuEFluxOsc->Sumw2();
280  fhNuEBarFluxOsc->Sumw2();
281  fhNuTauFluxOsc->Sumw2();
282  fhNuTauBarFluxOsc->Sumw2();
283 
284  fhNuMuFluxOsc_FastMC->Sumw2();
285  fhNuMuBarFluxOsc_FastMC->Sumw2();
286  fhNuEFluxOsc_FastMC->Sumw2();
287  fhNuEBarFluxOsc_FastMC->Sumw2();
288  fhNuTauFluxOsc_FastMC->Sumw2();
289  fhNuTauBarFluxOsc_FastMC->Sumw2();
290 
291  fhNuMuCCEventRate->Sumw2();
292  fhNuMuBarCCEventRate->Sumw2();
293  fhNuECCEventRate->Sumw2();
294  fhNuEBarCCEventRate->Sumw2();
295  fhNuTauCCEventRate->Sumw2();
296  fhNuTauBarCCEventRate->Sumw2();
297  fhNuMuNCEventRate->Sumw2();
298  fhNuMuBarNCEventRate->Sumw2();
299  fhNuENCEventRate->Sumw2();
300  fhNuEBarNCEventRate->Sumw2();
301  fhNuTauNCEventRate->Sumw2();
302  fhNuTauBarNCEventRate->Sumw2();
303 
304  fhNuMuCCEventRate_FastMC->Sumw2();
305  fhNuMuBarCCEventRate_FastMC->Sumw2();
306  fhNuECCEventRate_FastMC->Sumw2();
307  fhNuEBarCCEventRate_FastMC->Sumw2();
308  fhNuTauCCEventRate_FastMC->Sumw2();
309  fhNuTauBarCCEventRate_FastMC->Sumw2();
310  fhNuMuNCEventRate_FastMC->Sumw2();
311  fhNuMuBarNCEventRate_FastMC->Sumw2();
312  fhNuENCEventRate_FastMC->Sumw2();
313  fhNuEBarNCEventRate_FastMC->Sumw2();
314  fhNuTauNCEventRate_FastMC->Sumw2();
315  fhNuTauBarNCEventRate_FastMC->Sumw2();
316 
317  fhNuMuCCEventRateOsc->Sumw2();
318  fhNuMuBarCCEventRateOsc->Sumw2();
319  fhNuECCEventRateOsc->Sumw2();
320  fhNuEBarCCEventRateOsc->Sumw2();
321  fhNuTauCCEventRateOsc->Sumw2();
322  fhNuTauBarCCEventRateOsc->Sumw2();
323  fhNuMuNCEventRateOsc->Sumw2();
324  fhNuMuBarNCEventRateOsc->Sumw2();
325  fhNuENCEventRateOsc->Sumw2();
326  fhNuEBarNCEventRateOsc->Sumw2();
327  fhNuTauNCEventRateOsc->Sumw2();
328  fhNuTauBarNCEventRateOsc->Sumw2();
329 
330  fhNuMuCCEventRateOsc_FastMC->Sumw2();
331  fhNuMuBarCCEventRateOsc_FastMC->Sumw2();
332  fhNuECCEventRateOsc_FastMC->Sumw2();
333  fhNuEBarCCEventRateOsc_FastMC->Sumw2();
334  fhNuTauCCEventRateOsc_FastMC->Sumw2();
335  fhNuTauBarCCEventRateOsc_FastMC->Sumw2();
336  fhNuMuNCEventRateOsc_FastMC->Sumw2();
337  fhNuMuBarNCEventRateOsc_FastMC->Sumw2();
338  fhNuENCEventRateOsc_FastMC->Sumw2();
339  fhNuEBarNCEventRateOsc_FastMC->Sumw2();
340  fhNuTauNCEventRateOsc_FastMC->Sumw2();
341  fhNuTauBarNCEventRateOsc_FastMC->Sumw2();
342 
343 
344  //
345  // Set histogram titles
346  //
347 
348  std::string fluxtitle = "Neutrinos / GeV / m^{2} / POT";
349  std::string oscfluxtitle = "Oscillated Neutrinos / GeV / m^{2} / POT";
350  std::string cceventratetitle = "CC Events / POT";
351  std::string nceventratetitle = "CC Events / POT";
352 
353  SetTitles(fhNuMuFlux, "#nu_{#mu} Energy (GeV) ", "Unosc #nu_{#mu}s / GeV / m^{2} / POT");
354  SetTitles(fhNuMuBarFlux, "#bar{#nu}_{#mu} Energy (GeV)", "Unosc #bar{#nu}_{#mu}s / GeV / m^{2} / POT");
355  SetTitles(fhNuEFlux, "#nu_{e} Energy (GeV)", "Unosc #nu_{e}s / GeV / m^{2} / POT");
356  SetTitles(fhNuEBarFlux, "#bar{#nu}_{e} Energy (GeV)", "Unosc #bar{#nu}_{e}s / GeV / m^{2} / POT");
357  SetTitles(fhNuTauFlux, "#nu_{#tau} Energy (GeV)", "Unosc #nu_{#tau}s / GeV / m^{2} / POT");
358  SetTitles(fhNuTauBarFlux, "#bar{#nu}_{#tau} Energy (GeV)", "Unosc #bar{#nu}_{#tau}s / GeV / m^{2} / POT");
359 
360  SetTitles(fhNuMuFlux_FastMC, "#nu_{#mu} Energy (GeV) ", "Unosc #nu_{#mu}s / m^{2} / POT");
361  SetTitles(fhNuMuBarFlux_FastMC, "#bar{#nu}_{#mu} Energy (GeV)", "Unosc #bar{#nu}_{#mu}s / m^{2} / POT");
362  SetTitles(fhNuEFlux_FastMC, "#nu_{e} Energy (GeV)", "Unosc #nu_{e}s / m^{2} / POT");
363  SetTitles(fhNuEBarFlux_FastMC, "#bar{#nu}_{e} Energy (GeV)", "Unosc #bar{#nu}_{e}s / m^{2} / POT");
364  SetTitles(fhNuTauFlux_FastMC, "#nu_{#tau} Energy (GeV)", "Unosc #nu_{#tau}s / m^{2} / POT");
365  SetTitles(fhNuTauBarFlux_FastMC, "#bar{#nu}_{#tau} Energy (GeV)", "Unosc #bar{#nu}_{#tau}s / m^{2} / POT");
366 
367  SetTitles(fhNuMuFlux_Globes, "#nu_{#mu} Energy (GeV)", "Unosc #nu_{#mu}s / GeV / m^{2} / POT");
368  SetTitles(fhNuMuBarFlux_Globes, "#bar{#nu}_{#mu} Energy (GeV)", "Unosc #bar{#nu}_{#mu}s / GeV / m^{2} / POT");
369  SetTitles(fhNuEFlux_Globes, "#nu_{e} Energy (GeV)", "Unosc #nu_{e}s / GeV / m^{2} / POT");
370  SetTitles(fhNuEBarFlux_Globes, "#bar{#nu}_{e} Energy (GeV)", "Unosc #bar{#nu}_{e}s / GeV / m^{2} / POT");
371  SetTitles(fhNuTauFlux_Globes, "#nu_{#tau} Energy (GeV)", "Unosc #nu_{#tau}s / GeV / m^{2} / POT");
372  SetTitles(fhNuTauBarFlux_Globes, "#bar{#nu}_{#tau} Energy (GeV)", "Unosc #bar{#nu}_{#tau}s / GeV / m^{2} / POT");
373 
374  SetTitles(fhNuMuFluxOsc, "Energy (GeV)", "Oscillated #nu_{#mu}s / GeV / m^{2} / POT");
375  SetTitles(fhNuMuBarFluxOsc, "Energy (GeV)", "Oscillated #bar{#nu}_{#mu}s / GeV / m^{2} / POT");
376  SetTitles(fhNuEFluxOsc, "Energy (GeV)", "Oscillated #nu_{e}s / GeV / m^{2} / POT");
377  SetTitles(fhNuEBarFluxOsc, "Energy (GeV)", "Oscillated #bar{#nu}_{e}s / GeV / m^{2} / POT");
378  SetTitles(fhNuTauFluxOsc, "Energy (GeV)", "Oscillated #nu_{#tau}s / GeV / m^{2} / POT");
379  SetTitles(fhNuTauBarFluxOsc, "Energy (GeV)", "Oscillated #bar{#nu}_{#tau}s / GeV / m^{2} / POT");
380 
381  SetTitles(fhNuMuFluxOsc_FastMC, "Energy (GeV)", "Oscillated #nu_{#mu}s / m^{2} / POT");
382  SetTitles(fhNuMuBarFluxOsc_FastMC, "Energy (GeV)", "Oscillated #bar{#nu}_{#mu}s / m^{2} / POT");
383  SetTitles(fhNuEFluxOsc_FastMC, "Energy (GeV)", "Oscillated #nu_{e}s / m^{2} / POT");
384  SetTitles(fhNuEBarFluxOsc_FastMC, "Energy (GeV)", "Oscillated #bar{#nu}_{e}s / m^{2} / POT");
385  SetTitles(fhNuTauFluxOsc_FastMC, "Energy (GeV)", "Oscillated #nu_{#tau}s / GeV / m^{2} / POT");
386  SetTitles(fhNuTauBarFluxOsc_FastMC, "Energy (GeV)", "Oscillated #bar{#nu}_{#tau}s / GeV / m^{2} / POT");
387 
388 
389  SetTitles(fhNuMuCCEventRate, "Energy (GeV)", "#nu_{#mu} CC Events / GeV / kTon / POT");
390  SetTitles(fhNuMuBarCCEventRate, "Energy (GeV)", "#bar{#nu}_{#mu} CC Events / kTon / POT");
391  SetTitles(fhNuECCEventRate, "Energy (GeV)", "#nu_{e} CC Events / GeV / kTon / POT");
392  SetTitles(fhNuEBarCCEventRate, "Energy (GeV)", "#bar{#nu}_{e} CC Events / kTon / POT");
393  SetTitles(fhNuTauCCEventRate, "Energy (GeV)", "#nu_{#tau} CC Events / kTon / POT");
394  SetTitles(fhNuTauBarCCEventRate, "Energy (GeV)", "#bar{#nu}_{#tau} CC Events / kTon / POT");
395 
396  SetTitles(fhNuMuCCEventRate_FastMC, "Energy (GeV)", "#nu_{#mu} CC Events / kTon / POT");
397  SetTitles(fhNuMuBarCCEventRate_FastMC, "Energy (GeV)", "#bar{#nu}_{#mu} CC Events / kTon / POT");
398  SetTitles(fhNuECCEventRate_FastMC, "Energy (GeV)", "#nu_{e} CC Events / kTon / POT");
399  SetTitles(fhNuEBarCCEventRate_FastMC, "Energy (GeV)", "#bar{#nu}_{e} CC Events / kTon / POT");
400  SetTitles(fhNuTauCCEventRate_FastMC, "Energy (GeV)", "#nu_{#tau} CC Events / kTon / POT");
401  SetTitles(fhNuTauBarCCEventRate_FastMC, "Energy (GeV)", "#bar{#nu}_{#tau} CC Events / kTon / POT");
402 
403  SetTitles(fhNuMuCCEventRateOsc, "Energy (GeV)", "Oscillated #nu_{#mu} CC Events / GeV / kTon / POT");
404  SetTitles(fhNuMuBarCCEventRateOsc, "Energy (GeV)", "Oscillated #bar{#nu}_{#mu} CC Events / GeV / kTon / POT");
405  SetTitles(fhNuECCEventRateOsc, "Energy (GeV)", "Oscillated #nu_{e} CC Events / GeV / kTon / POT");
406  SetTitles(fhNuEBarCCEventRateOsc, "Energy (GeV)", "Oscillated #bar{#nu}_{e} CC Events / GeV / kTon / POT");
407  SetTitles(fhNuTauCCEventRateOsc, "Energy (GeV)", "Oscillated #nu_{#tau} CC Events / GeV / kTon / POT");
408  SetTitles(fhNuTauBarCCEventRateOsc, "Energy (GeV)", "Oscillated #bar{#nu}_{#tau} CC Events / GeV / kTon / POT");
409 
410  SetTitles(fhNuMuCCEventRateOsc_FastMC, "Energy (GeV)", "Oscillated #nu_{#mu} CC Events / kTon / POT");
411  SetTitles(fhNuMuBarCCEventRateOsc_FastMC, "Energy (GeV)", "Oscillated #bar{#nu}_{#mu} CC Events / kTon / POT");
412  SetTitles(fhNuECCEventRateOsc_FastMC, "Energy (GeV)", "Oscillated #nu_{e} CC Events / kTon / POT");
413  SetTitles(fhNuEBarCCEventRateOsc_FastMC, "Energy (GeV)", "Oscillated #bar{#nu}_{e} CC Events / kTon / POT");
414  SetTitles(fhNuTauCCEventRateOsc_FastMC, "Energy (GeV)", "Oscillated #nu_{#tau} CC Events / kTon / POT");
415  SetTitles(fhNuTauBarCCEventRateOsc_FastMC, "Energy (GeV)", "Oscillated #bar{#nu}_{#tau} CC Events / kTon / POT");
416 
417 
418  SetTitles(fhNuMuNCEventRate, "Energy (GeV)", "#nu_{#mu} NC Events / GeV / kTon / POT");
419  SetTitles(fhNuMuBarNCEventRate, "Energy (GeV)","#bar{#nu}_{#mu} NC Events / GeV / kTon / POT");
420  SetTitles(fhNuENCEventRate, "Energy (GeV)", "#nu_{e} NC Events / GeV / kTon / POT");
421  SetTitles(fhNuEBarNCEventRate, "Energy (GeV)", "#bar{#nu}_{e} NC Events / GeV / kTon / POT");
422  SetTitles(fhNuTauNCEventRate, "Energy (GeV)", "#nu_{#tau} NC Events / GeV / kTon / POT");
423  SetTitles(fhNuTauBarNCEventRate, "Energy (GeV)", "#bar{#nu}_{#tau} NC Events / GeV / kTon / POT");
424 
425  SetTitles(fhNuMuNCEventRate_FastMC, "Energy (GeV)", "#nu_{#mu} NC Events / kTon / POT");
426  SetTitles(fhNuMuBarNCEventRate_FastMC, "#bar{#nu}_{#mu} Energy (GeV)","#bar{#nu}_{#mu} NC Events / kTon / POT");
427  SetTitles(fhNuENCEventRate_FastMC, "#nu_{e} Energy (GeV)", "#nu_{e} NC Events / kTon / POT");
428  SetTitles(fhNuEBarNCEventRate_FastMC, "#bar{#nu}_{e} Energy (GeV)", "#bar{#nu}_{e} NC Events / kTon / POT");
429  SetTitles(fhNuTauNCEventRate_FastMC, "#nu_{#tau} Energy (GeV)", "#nu_{#tau} NC Events / kTon / POT");
430  SetTitles(fhNuTauBarNCEventRate_FastMC, "#bar{#nu}_{#tau} Energy (GeV)", "#bar{#nu}_{#tau} NC Events / kTon / POT");
431 
432  SetTitles(fhNuMuNCEventRateOsc, "#nu_{#mu} Energy (GeV)", "Oscillated #nu_{#mu} NC Events / GeV / kTon / POT");
433  SetTitles(fhNuMuBarNCEventRateOsc, "#bar{#nu}_{#mu} Energy (GeV)","Oscillated #bar{#nu}_{#mu} NC Events / GeV / kTon / POT");
434  SetTitles(fhNuENCEventRateOsc, "#nu_{e} Energy (GeV)", "Oscillated #nu_{e} NC Events / GeV / kTon / POT");
435  SetTitles(fhNuEBarNCEventRateOsc, "#bar{#nu}_{e} Energy (GeV)", "Oscillated #bar{#nu}_{e} NC Events / GeV / kTon / POT");
436  SetTitles(fhNuTauNCEventRateOsc, "#nu_{#tau} Energy (GeV)", "Oscillated #nu_{#tau} NC Events / GeV / kTon / POT");
437  SetTitles(fhNuTauBarNCEventRateOsc, "#bar{#nu}_{#tau} Energy (GeV)", "Oscillated #bar{#nu}_{#tau} NC Events / GeV / kTon / POT");
438 
439  SetTitles(fhNuMuNCEventRateOsc_FastMC, "#nu_{#mu} Energy (GeV)", "Oscillated #nu_{#mu} NC Events / kTon / POT");
440  SetTitles(fhNuMuBarNCEventRateOsc_FastMC, "#bar{#nu}_{#mu} Energy (GeV)","Oscillated #bar{#nu}_{#mu} NC Events / kTon / POT");
441  SetTitles(fhNuENCEventRateOsc_FastMC, "#nu_{e} Energy (GeV)", "Oscillated #nu_{e} NC Events / kTon / POT");
442  SetTitles(fhNuEBarNCEventRateOsc_FastMC, "#bar{#nu}_{e} Energy (GeV)", "Oscillated #bar{#nu}_{e} NC Events / kTon / POT");
443  SetTitles(fhNuTauNCEventRateOsc_FastMC, "#nu_{#tau} Energy (GeV)", "Oscillated #nu_{#tau} NC Events / kTon / POT");
444  SetTitles(fhNuTauBarNCEventRateOsc_FastMC, "#bar{#nu}_{#tau} Energy (GeV)", "Oscillated #bar{#nu}_{#tau} NC Events / kTon / POT");
445 
446  //
447  //start loop over entries in ntuple
448  //
449 
450  dk2nu = new bsim::Dk2Nu;
451  if(isDk2nu)
452  fChain->SetBranchAddress("dk2nu",&dk2nu);
453 
454 
455  Long64_t nentries = fChain->GetEntries();
456  std::cout << "Total number of Entries = " << nentries << std::endl;
457 
458  for (Long64_t jentry=0; jentry<nentries;jentry++)
459  {
460  if(isDk2nu) {
461  //Long64_t ientry = LoadTree(jentry);
462  fChain->GetEntry(jentry);
463  //if (ientry < 0) break;
464  }
465  else {
466  Long64_t ientry = LoadTree(jentry);
467  fChain->GetEntry(jentry);
468  if (ientry < 0) break;
469  }
470  ++iread;
471 
472  if(iread % 10000 == 0)
473  {
474  std::cout << "Reading Entry " << iread << std::endl;
475  }
476 
477  //
478  // Compute the detector location weight
479  // This calculation needs to know the coordinates
480  // of the detector where you want to plot the flux
481  // which is set by the eventRates constructor
482  //
483 
484  double nuenergyatsomedet = -999.0;
485  double detectorwghtatsomedet = -999.0;
486  std::vector<double> detvec;
487  detvec.push_back(detx);
488  detvec.push_back(dety);
489  detvec.push_back(detz);
490 
491  //this function computes the location weight and neutrino energy at detx,y,z
492  eventRates::GetWeight(detvec, detectorwghtatsomedet, nuenergyatsomedet);
493 
494  // Calculate the total weight, including location weight, importance weight
495  // and POT scale factor
496 
497  if(isDk2nu)
498  Nimpwt = dk2nu->decay.nimpwt;
499  double fluxwghtsomedet = (detectorwghtatsomedet*Nimpwt/3.1415)*(refpot/fTotalPOT);
500  // Multiply flux weights by cross sections to get event rate weights
501 
502  std::string current_string = "NC";
503  double cceventratewghtsomedet = fluxwghtsomedet * GetXSec((double)Ntype,nuenergyatsomedet,current_string);
504 
505  current_string = "NC";
506  double nceventratewghtsomedet = fluxwghtsomedet * GetXSec(Ntype,nuenergyatsomedet,current_string);
507 
508  // Get Oscillated Flavor
509  int Ntype_osc = GetOscillatedNeutrinoType(Ntype,nuenergyatsomedet);
510 
511  current_string = "CC";
512  double cceventratewghtsomedet_osc = fluxwghtsomedet * GetXSec((double)Ntype_osc,nuenergyatsomedet,current_string);
513  current_string = "NC";
514  double nceventratewghtsomedet_osc = fluxwghtsomedet * GetXSec(Ntype_osc,nuenergyatsomedet,current_string);
515  ///
516 
517  //Now fill the histograms
518  //by neutrino type
519  //
520  if(isDk2nu)
521  Ntype = dk2nu->decay.ntype;
522  if(Ntype == 56 || Ntype == 14)
523  {
524  fhNuMuFlux -> Fill(nuenergyatsomedet, fluxwghtsomedet);
525  fhNuMuCCEventRate -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
526  fhNuMuNCEventRate -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
527 
528  fhNuMuFlux_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
529  fhNuMuCCEventRate_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
530  fhNuMuNCEventRate_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
531  fhNuMuFlux_Globes-> Fill(nuenergyatsomedet, fluxwghtsomedet);
532  }
533  if(Ntype == 55 || Ntype==-14)
534  {
535  fhNuMuBarFlux -> Fill(nuenergyatsomedet, fluxwghtsomedet);
536  fhNuMuBarCCEventRate -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
537  fhNuMuBarNCEventRate -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
538 
539  fhNuMuBarFlux_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
540  fhNuMuBarCCEventRate_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
541  fhNuMuBarNCEventRate_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
542  fhNuMuBarFlux_Globes-> Fill(nuenergyatsomedet, fluxwghtsomedet);
543  }
544  if(Ntype == 53 || Ntype==12)
545  {
546  fhNuEFlux -> Fill(nuenergyatsomedet, fluxwghtsomedet);
547  fhNuECCEventRate -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
548  fhNuENCEventRate -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
549 
550  fhNuEFlux_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
551  fhNuECCEventRate_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
552  fhNuENCEventRate_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
553  fhNuEFlux_Globes-> Fill(nuenergyatsomedet, fluxwghtsomedet);
554  }
555  if(Ntype == 52 || Ntype==-12)
556  {
557  fhNuEBarFlux -> Fill(nuenergyatsomedet, fluxwghtsomedet);
558  fhNuEBarCCEventRate -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
559  fhNuEBarNCEventRate -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
560 
561  fhNuEBarFlux_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
562  fhNuEBarCCEventRate_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
563  fhNuEBarNCEventRate_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
564  fhNuEBarFlux_Globes-> Fill(nuenergyatsomedet, fluxwghtsomedet);
565  }
566 
567  if(Ntype == 59 || Ntype==16)
568  {
569  fhNuTauFlux -> Fill(nuenergyatsomedet, fluxwghtsomedet);
570  fhNuTauCCEventRate -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
571  fhNuTauNCEventRate -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
572 
573  fhNuTauFlux_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
574  fhNuTauCCEventRate_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
575  fhNuTauNCEventRate_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
576  fhNuTauFlux_Globes-> Fill(nuenergyatsomedet, fluxwghtsomedet);
577  }
578  if(Ntype == 59 || Ntype==-16)
579  {
580  fhNuTauBarFlux -> Fill(nuenergyatsomedet, fluxwghtsomedet);
581  fhNuTauBarCCEventRate -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
582  fhNuTauBarNCEventRate -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
583 
584  fhNuTauBarFlux_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
585  fhNuTauBarCCEventRate_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet);
586  fhNuTauBarNCEventRate_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet);
587  fhNuTauBarFlux_Globes-> Fill(nuenergyatsomedet, fluxwghtsomedet);
588  }
589 
590  if(Ntype_osc == 56 || Ntype_osc==14)
591  {
592  fhNuMuFluxOsc -> Fill(nuenergyatsomedet, fluxwghtsomedet);
593  fhNuMuCCEventRateOsc -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
594  fhNuMuNCEventRateOsc -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
595  fhNuMuFluxOsc_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
596  fhNuMuCCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
597  fhNuMuNCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
598  }
599  if(Ntype_osc == 55 || Ntype_osc==-14)
600  {
601  fhNuMuBarFluxOsc -> Fill(nuenergyatsomedet, fluxwghtsomedet);
602  fhNuMuBarCCEventRateOsc -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
603  fhNuMuBarNCEventRateOsc -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
604  fhNuMuBarFluxOsc_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
605  fhNuMuBarCCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
606  fhNuMuBarNCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
607 
608  }
609  if(Ntype_osc == 53 || Ntype_osc==12)
610  {
611  fhNuEFluxOsc -> Fill(nuenergyatsomedet, fluxwghtsomedet);
612  fhNuECCEventRateOsc -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
613  fhNuENCEventRateOsc -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
614  fhNuEFluxOsc_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
615  fhNuECCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
616  fhNuENCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
617  }
618  if(Ntype_osc == 52 || Ntype_osc==-12)
619  {
620  fhNuEBarFluxOsc -> Fill(nuenergyatsomedet, fluxwghtsomedet);
621  fhNuEBarCCEventRateOsc -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
622  fhNuEBarNCEventRateOsc -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
623  fhNuEBarFluxOsc_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
624  fhNuEBarCCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
625  fhNuEBarNCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
626 
627  }
628 
629  if(Ntype_osc == 59 || Ntype_osc==16)
630  {
631  fhNuTauFluxOsc -> Fill(nuenergyatsomedet, fluxwghtsomedet);
632  fhNuTauCCEventRateOsc -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
633  fhNuTauNCEventRateOsc -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
634  fhNuTauFluxOsc_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
635  fhNuTauCCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
636  fhNuTauNCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
637 
638  }
639  if(Ntype_osc == 58 || Ntype_osc==-16)
640  {
641  fhNuTauBarFluxOsc -> Fill(nuenergyatsomedet, fluxwghtsomedet);
642  fhNuTauBarCCEventRateOsc -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
643  fhNuTauBarNCEventRateOsc -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
644  fhNuTauBarFluxOsc_FastMC -> Fill(nuenergyatsomedet, fluxwghtsomedet);
645  fhNuTauBarCCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, cceventratewghtsomedet_osc);
646  fhNuTauBarNCEventRateOsc_FastMC -> Fill(nuenergyatsomedet, nceventratewghtsomedet_osc);
647  }
648  }//end loop over entries
649 
650  // normalize by bin width for plots but not for fastmc
651  fhNuMuFlux->Scale(1.0,"width");
652  fhNuMuCCEventRate->Scale(1.0,"width");
653  fhNuMuNCEventRate->Scale(1.0,"width");
654 
655  fhNuMuBarFlux->Scale(1.0,"width");
656  fhNuMuBarCCEventRate->Scale(1.0,"width");
657  fhNuMuBarNCEventRate->Scale(1.0,"width");
658 
659  fhNuEFlux->Scale(1.0,"width");
660  fhNuECCEventRate->Scale(1.0,"width");
661  fhNuENCEventRate->Scale(1.0,"width");
662 
663  fhNuEBarFlux->Scale(1.0,"width");
664  fhNuEBarCCEventRate->Scale(1.0,"width");
665  fhNuEBarNCEventRate->Scale(1.0,"width");
666 
667  fhNuTauFlux->Scale(1.0,"width");
668  fhNuTauCCEventRate->Scale(1.0,"width");
669  fhNuTauNCEventRate->Scale(1.0,"width");
670 
671  fhNuTauBarFlux->Scale(1.0,"width");
672  fhNuTauBarCCEventRate->Scale(1.0,"width");
673  fhNuTauBarNCEventRate->Scale(1.0,"width");
674 
675  fhNuMuFluxOsc->Scale(1.0,"width");
676  fhNuMuCCEventRateOsc->Scale(1.0,"width");
677  fhNuMuNCEventRateOsc->Scale(1.0,"width");
678 
679  fhNuMuBarFluxOsc->Scale(1.0,"width");
680  fhNuMuBarCCEventRateOsc->Scale(1.0,"width");
681  fhNuMuBarNCEventRateOsc->Scale(1.0,"width");
682 
683  fhNuEFluxOsc->Scale(1.0,"width");
684  fhNuECCEventRateOsc->Scale(1.0,"width");
685  fhNuENCEventRateOsc->Scale(1.0,"width");
686 
687  fhNuEBarFluxOsc->Scale(1.0,"width");
688  fhNuEBarCCEventRateOsc->Scale(1.0,"width");
689  fhNuEBarNCEventRateOsc->Scale(1.0,"width");
690 
691  fhNuTauFluxOsc->Scale(1.0,"width");
692  fhNuTauCCEventRateOsc->Scale(1.0,"width");
693  fhNuTauNCEventRateOsc->Scale(1.0,"width");
694 
695  fhNuTauBarFluxOsc->Scale(1.0,"width");
696  fhNuTauBarCCEventRateOsc->Scale(1.0,"width");
697  fhNuTauBarNCEventRateOsc->Scale(1.0,"width");
698 
699  fhNuMuFlux_Globes->Scale(1.0,"width");
700  fhNuMuBarFlux_Globes->Scale(1.0,"width");
701  fhNuEFlux_Globes->Scale(1.0,"width");
702  fhNuEBarFlux_Globes->Scale(1.0,"width");
703  fhNuTauFlux_Globes->Scale(1.0,"width");
704  fhNuTauBarFlux_Globes->Scale(1.0,"width");
705 
706  //
707  // Style histograms
708  //
709 
710  TGaxis::SetMaxDigits(2);
711  fhNuMuFlux_FastMC->GetYaxis()->SetTitleOffset(1.4);
712  fhNuMuBarFlux_FastMC->GetYaxis()->SetTitleOffset(1.4);
713  fhNuEFlux_FastMC->GetYaxis()->SetTitleOffset(1.4);
714  fhNuEBarFlux_FastMC->GetYaxis()->SetTitleOffset(1.4);
715  fhNuTauFlux_FastMC->GetYaxis()->SetTitleOffset(1.4);
716  fhNuTauBarFlux_FastMC->GetYaxis()->SetTitleOffset(1.4);
717 
718  fhNuMuCCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
719  fhNuMuBarCCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
720  fhNuECCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
721  fhNuEBarCCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
722  fhNuTauCCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
723  fhNuTauBarCCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
724 
725  fhNuMuNCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
726  fhNuMuBarNCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
727  fhNuENCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
728  fhNuEBarNCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
729  fhNuTauNCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
730  fhNuTauBarNCEventRate_FastMC->GetYaxis()->SetTitleOffset(1.4);
731 
732  // save histograms to file
733  ffilename = ffilename + "_" + detectorname;
734 
735 
736  // put the location of the detector in the file
737  TParameter<double> det_x("det_x",detx);
738  TParameter<double> det_y("det_y",dety);
739  TParameter<double> det_z("det_z",detz);
740 
741  // Save histograms to a root file
742  std::cout<<"writing"+ffilename+".root"<<std::endl;
743  TFile f((ffilename+".root").c_str(),"recreate");
744  det_x.Write();
745  det_y.Write();
746  det_z.Write();
747  fhNuMuFlux->Write();
748  fhNuMuBarFlux->Write();
749  fhNuEFlux->Write();
750  fhNuEBarFlux->Write();
751  fhNuTauFlux->Write();
752  fhNuTauBarFlux->Write();
753  fhNuMuCCEventRate->Write();
754  fhNuMuBarCCEventRate->Write();
755  fhNuECCEventRate->Write();
756  fhNuEBarCCEventRate->Write();
757  fhNuTauCCEventRate->Write();
758  fhNuTauBarCCEventRate->Write();
759  fhNuMuNCEventRate->Write();
760  fhNuMuBarNCEventRate->Write();
761  fhNuENCEventRate->Write();
762  fhNuEBarNCEventRate->Write();
763  fhNuTauNCEventRate->Write();
764  fhNuTauBarNCEventRate->Write();
765  fhNuMuFluxOsc->Write();
766  fhNuMuBarFluxOsc->Write();
767  fhNuEFluxOsc->Write();
768  fhNuEBarFluxOsc->Write();
769  fhNuTauFluxOsc->Write();
770  fhNuTauBarFluxOsc->Write();
771  fhNuMuCCEventRateOsc->Write();
772  fhNuMuBarCCEventRateOsc->Write();
773  fhNuECCEventRateOsc->Write();
774  fhNuEBarCCEventRateOsc->Write();
775  fhNuTauCCEventRateOsc->Write();
776  fhNuTauBarCCEventRateOsc->Write();
777  fhNuMuNCEventRateOsc->Write();
778  fhNuMuBarNCEventRateOsc->Write();
779  fhNuENCEventRateOsc->Write();
780  fhNuEBarNCEventRateOsc->Write();
781  fhNuTauNCEventRateOsc->Write();
782  fhNuTauBarNCEventRateOsc->Write();
783  f.Close();
784 
785  TFile g((ffilename+"_fastmc.root").c_str(),"recreate");
786  det_x.Write();
787  det_y.Write();
788  det_z.Write();
789  fhNuMuFlux_FastMC->Write();
790  fhNuMuBarFlux_FastMC->Write();
791  fhNuEFlux_FastMC->Write();
792  fhNuEBarFlux_FastMC->Write();
793  fhNuTauFlux_FastMC->Write();
794  fhNuTauBarFlux_FastMC->Write();
795  fhNuMuCCEventRate_FastMC->Write();
796  fhNuMuBarCCEventRate_FastMC->Write();
797  fhNuECCEventRate_FastMC->Write();
798  fhNuEBarCCEventRate_FastMC->Write();
799  fhNuTauCCEventRate_FastMC->Write();
800  fhNuTauBarCCEventRate_FastMC->Write();
801  fhNuMuNCEventRate_FastMC->Write();
802  fhNuMuBarNCEventRate_FastMC->Write();
803  fhNuENCEventRate_FastMC->Write();
804  fhNuEBarNCEventRate_FastMC->Write();
805  fhNuTauNCEventRate_FastMC->Write();
806  fhNuTauBarNCEventRate_FastMC->Write();
807 
808  fhNuMuFluxOsc_FastMC->Write();
809  fhNuMuBarFluxOsc_FastMC->Write();
810  fhNuEFluxOsc_FastMC->Write();
811  fhNuEBarFluxOsc_FastMC->Write();
812  fhNuTauFluxOsc_FastMC->Write();
813  fhNuTauBarFluxOsc_FastMC->Write();
814  fhNuMuCCEventRateOsc_FastMC->Write();
815  fhNuMuBarCCEventRateOsc_FastMC->Write();
816  fhNuECCEventRateOsc_FastMC->Write();
817  fhNuEBarCCEventRateOsc_FastMC->Write();
818  fhNuTauCCEventRateOsc_FastMC->Write();
819  fhNuTauBarCCEventRateOsc_FastMC->Write();
820  fhNuMuNCEventRateOsc_FastMC->Write();
821  fhNuMuBarNCEventRateOsc_FastMC->Write();
822  fhNuENCEventRateOsc_FastMC->Write();
823  fhNuEBarNCEventRateOsc_FastMC->Write();
824  fhNuTauNCEventRateOsc_FastMC->Write();
825  fhNuTauBarNCEventRateOsc_FastMC->Write();
826  g.Close();
827 
828  // Write Globes files
829  ofstream myfile;
830  myfile.open((ffilename+"_globes_flux.txt").c_str());
831  for(int i = 0; i<globes_nbins; i++) {
832  myfile<<fhNuMuFlux_Globes->GetBinCenter(i+1)<<" "<<
833  fhNuEFlux_Globes->GetBinContent(i+1)<<" "<<
834  fhNuMuFlux_Globes->GetBinContent(i+1)<<" "<<
835  fhNuTauFlux_Globes->GetBinContent(i+1)<<" "<<
836  fhNuEBarFlux_Globes->GetBinContent(i+1)<<" "<<
837  fhNuMuBarFlux_Globes->GetBinContent(i+1)<<" "<<
838  fhNuTauBarFlux_Globes->GetBinContent(i+1)<<" "<<std::endl;
839  }
840 
841 }
842 
843 //-------------------------------------------------------------------------------------
844 std::string eventRates::GetPOTAsString(const double dpot)
845 {
846 
847  std::stringstream potstrm;
848  potstrm << scientific << dpot;
849 
850  string potstr = potstrm.str();
851 
852  //
853  //get base
854  //
855  size_t baselength;
856  if(potstr.find("e",0) != string::npos)
857  {
858  baselength = potstr.find("e",0);
859  }
860  else if(potstr.find("E",0) != string::npos)
861  {
862  baselength = potstr.find("E",0);
863  }
864  else
865  {
866  cout << "eventRates::GetPOTAsString - PROBLEM: pot is not in scientific notation" << endl;
867  return "Problem";
868  }
869 
870  string base = potstr.substr(0, baselength);
871 
872  //
873  //get exp
874  //
875  size_t exppos;
876  if(potstr.find("+",baselength) != string::npos)
877  {
878  exppos = potstr.find("+",baselength);
879  }
880  else if(potstr.find("-",baselength) != string::npos)
881  {
882  exppos = potstr.find("-",baselength);
883  }
884  else
885  {
886  cout << "eventRates::GetPOTAsString - PROBLEM: pot is not in scientific notation" << endl;
887  return "Problem";
888  }
889 
890  string exp = potstr.substr(exppos);
891 
892 
893  //
894  //modify base string if needed
895  //
896  string baseNumber = base;
897 
898  size_t baseDecimalpos = base.find(".",0);
899  if(baseDecimalpos != string::npos)
900  {
901  size_t baseNotZeropos = base.find_last_not_of("0", string::npos);
902  if(baseNotZeropos != string::npos)
903  {
904  if(baseNotZeropos > baseDecimalpos)
905  {
906  baseNumber = base.substr(0,baseNotZeropos+1);
907  }
908  else
909  {
910  baseNumber = base.substr(0,baseDecimalpos+2);
911  }
912 
913  }
914  }
915  else
916  {
917  baseNumber = baseNumber + ".0";
918  }
919 
920  //
921  //modify exp string if needed
922  //
923 
924  string expSign = exp.substr(0, 1);
925  string expNumber = exp.substr(1, string::npos);
926 
927  size_t expNotZeropos = expNumber.find_first_not_of("0",0);
928  if(expNotZeropos != string::npos)
929  {
930  expNumber = expNumber.substr(expNotZeropos, string::npos);
931  }
932 
933 
934 
935  string potfinalstr;
936 
937  //
938  //put base and exp together
939  //
940 
941  if(baseNumber.empty() && expNumber.empty())
942  {
943  cout << "eventRates::GetPOTAsString - PROBLEM: base number and exp number are both empty" << endl;
944  return "Problem";
945  }
946 
947  if(baseNumber == "1.0")
948  {
949  if(expSign == "-")
950  potfinalstr = "10^{" + expSign + expNumber + "}";
951  else
952  potfinalstr = "10^{" + expNumber + "}";
953  }
954  else
955  {
956  if(expSign == "-")
957  potfinalstr = baseNumber + "#times10^{" + expSign + expNumber + "}";
958  else
959  potfinalstr = baseNumber + "#times10^{" + expNumber + "}";
960  }
961 
962 
963  /*
964  cout << "pot str = " << potstr << endl
965  << " base = " << base << endl
966  << " baseNumber = " << baseNumber << endl
967  << " expSign = " << expSign << endl
968  << " expNumber = " << expNumber << endl
969  << " potfinalstr = " << potfinalstr << endl;
970  */
971 
972 
973  return potfinalstr;
974 }
975 
976 
977 //---------------------------------------------------------------------------------------------
978 void eventRates::SetTitles(TH1* h, const std::string &xtitle, const std::string &ytitle)
979 {
980  if(!ytitle.empty())
981  {
982  h -> GetYaxis() -> SetTitle(ytitle.c_str());
983  h -> GetYaxis() -> CenterTitle();
984  }
985  if(!xtitle.empty())
986  {
987  h -> GetXaxis() -> SetTitle(xtitle.c_str());
988  h -> GetXaxis() -> CenterTitle();
989  }
990 }
991 
992 //---------------------------------------------------------------------------------------------
993 double eventRates::GetWeight(const std::vector<double> xdet,
994  double& nu_wght,
995  double& nu_energy)
996 {
997 
998 // if(iread > 60000)
999 // std::cout << "start iread = " << iread;
1000 
1001  //assumes units are GeV and cm
1002 
1003  const double rdet = 100.0; //in cm
1004  const double pimass = 0.13957; //in GeV
1005  const double kmass = 0.49368;
1006  const double k0mass = 0.49767;
1007  const double mumass = 0.105658389;
1008  const double taumass = 1.77682;
1009 
1010  //these are geant codes not PDG
1011  int nue = 53;
1012  int nuebar = 52;
1013  int numu = 56;
1014  int numubar = 55;
1015  int nutau = 59; // is this right?
1016  int nutaubar = 58; // is this right?
1017  int muplus = 5;
1018  int muminus = 6;
1019  if(isDk2nu) {
1020  nue = 12;
1021  nuebar = -12;
1022  numu = 14;
1023  numubar = -14;
1024  nutau = 16;
1025  nutaubar = -16;
1026  muplus = -13;
1027  muminus = 13;
1028  }
1029 
1030  // variables needed from the ntuple
1031  int m_ptype = -9999;
1032  double m_pdPx = -9999;
1033  double m_pdPy = -9999;
1034  double m_pdPz = -9999;
1035  double m_Vx = -9999;
1036  double m_Vy = -9999;
1037  double m_Vz = -9999;
1038  double m_ppdxdz = -9999;
1039  double m_ppdydz = -9999;
1040  double m_ppenergy = -9999;;
1041  double m_pppz = -9999;
1042  double m_muparpx = -9999;
1043  double m_muparpy = -9999;
1044  double m_muparpz = -9999;
1045  double m_mupare = -9999;
1046  int m_Ntype = -9999;
1047  double m_Necm = -9999;
1048 
1049  double m_parent_mass = -9999;
1050  if(isDk2nu) {
1051  m_ptype = dk2nu->decay.ptype;
1052  m_pdPx = dk2nu->decay.pdpx;
1053  m_pdPy = dk2nu->decay.pdpy;
1054  m_pdPz = dk2nu->decay.pdpz;
1055  m_Vx = dk2nu->decay.vx;
1056  m_Vy = dk2nu->decay.vy;
1057  m_Vz = dk2nu->decay.vz;
1058  m_ppdxdz = dk2nu->decay.ppdxdz;
1059  m_ppdydz = dk2nu->decay.ppdydz;
1060  m_ppenergy = dk2nu->decay.ppenergy;
1061  m_pppz = dk2nu->decay.pppz;
1062  m_muparpx = dk2nu->decay.muparpx;
1063  m_muparpy = dk2nu->decay.muparpy;
1064  m_muparpz = dk2nu->decay.muparpz;
1065  m_mupare = dk2nu->decay.mupare;
1066  m_Ntype = dk2nu->decay.ntype;
1067  m_Necm = dk2nu->decay.necm;
1068 
1069  if (m_ptype == 211 || m_ptype == -211) m_parent_mass = pimass;
1070  else if (m_ptype == 321 || m_ptype == -321) m_parent_mass = kmass;
1071  else if (m_ptype == 311 || m_ptype == 130 || m_ptype==310) m_parent_mass = k0mass;
1072  else if (m_ptype == 13 || m_ptype == -13) m_parent_mass = mumass;
1073  else {
1074  cout <<"eventRates::GetWeight - Wrong parent type!! "<< m_ptype << " = "
1075  << m_ptype << " Decay code = " << Ndecay <<endl;
1076 
1077  return -999;
1078  }
1079  }
1080  else{
1081  m_ptype = ptype;
1082  if (m_ptype == 8 || m_ptype == 9) {
1083  m_parent_mass = pimass;
1084  }
1085  else if (m_ptype == 11 || m_ptype == 12) m_parent_mass = kmass;
1086  else if (m_ptype == 10) m_parent_mass = k0mass;
1087  else if (m_ptype == 5 || m_ptype == 6) m_parent_mass = mumass;
1088  else {
1089  cout <<"eventRates::GetWeight - Wrong parent type!! "<< m_ptype << " = "
1090  << m_ptype << " Decay code = " << Ndecay <<endl;
1091  return -999;
1092  }
1093  m_pdPx = pdPx;
1094  m_pdPy = pdPy;
1095  m_pdPz = pdPz;
1096  m_Vx = Vx;
1097  m_Vy = Vy;
1098  m_Vz = Vz;
1099  m_ppdxdz = ppdxdz;
1100  m_ppdydz = ppdydz;
1101  m_ppenergy = ppenergy;
1102  m_pppz = pppz;
1103  m_muparpx = muparpx;
1104  m_muparpy = muparpy;
1105  m_muparpz = muparpz;
1106  m_mupare = mupare;
1107  m_Ntype = Ntype;
1108  m_Necm = Necm;
1109 
1110  }
1111 
1112  /*
1113  std::cout<<"ptype "<<m_ptype<<std::endl;
1114  std::cout<<"pdPx "<<m_pdPx<<std::endl;
1115  std::cout<<"pdPy "<<m_pdPy<<std::endl;
1116  std::cout<<"pdPz "<<m_pdPz<<std::endl;
1117  std::cout<<"Vx "<<m_Vx<<std::endl;
1118  std::cout<<"Vy "<<m_Vy<<std::endl;
1119  std::cout<<"Vz "<<m_Vz<<std::endl;
1120  std::cout<<"m_ppdxdz "<<m_ppdxdz<<std::endl;
1121  std::cout<<"m_ppdydz "<<m_ppdydz<<std::endl;
1122  std::cout<<"m_ppenergy "<<m_ppenergy<<std::endl;
1123  std::cout<<"m_pppz "<<m_pppz<<std::endl;
1124  std::cout<<"m_muparpx "<<m_muparpx<<std::endl;
1125  std::cout<<"m_muparpy "<<m_muparpy<<std::endl;
1126  std::cout<<"m_muparpz "<<m_muparpz<<std::endl;
1127  std::cout<<"m_mupare "<<m_mupare<<std::endl;
1128  std::cout<<"m_Ntype "<<m_Ntype<<std::endl;
1129  std::cout<<"parent mass "<<m_parent_mass<<std::endl;
1130  */
1131 
1132  double parent_energy = sqrt(m_pdPx*m_pdPx +
1133  m_pdPy*m_pdPy +
1134  m_pdPz*m_pdPz +
1135  m_parent_mass*m_parent_mass);
1136  double gamma = parent_energy / m_parent_mass;
1137  double gamma_sqr = gamma*gamma;
1138  double beta_mag = sqrt((gamma_sqr-1.)/gamma_sqr);
1139 
1140  double enuzr = m_Necm;
1141 
1142  double rad = sqrt((xdet[0]-m_Vx)*(xdet[0]-m_Vx) +
1143  (xdet[1]-m_Vy)*(xdet[1]-m_Vy) +
1144  (xdet[2]-m_Vz)*(xdet[2]-m_Vz));
1145 
1146  double parentp = sqrt((m_pdPx*m_pdPx)+
1147  (m_pdPy*m_pdPy)+
1148  (m_pdPz*m_pdPz));
1149  double costh_pardet = (m_pdPx*(xdet[0]-m_Vx) +
1150  m_pdPy*(xdet[1]-m_Vy) +
1151  m_pdPz*(xdet[2]-m_Vz))/(parentp*rad);
1152 
1153  if (costh_pardet>1.) costh_pardet = 1.;
1154  else if (costh_pardet<-1.) costh_pardet = -1.;
1155  double theta_pardet = acos(costh_pardet);
1156 
1157  double emrat = 1./(gamma*(1. - beta_mag * cos(theta_pardet)));
1158 
1159  nu_energy = emrat*enuzr;
1160 
1161  double sangdet = (rdet*rdet /(rad*rad)/ 4.);
1162 
1163  nu_wght = sangdet * emrat * emrat;
1164 
1165  //done for all except polarized muon
1166  // in which case need to modify weight
1167  if (m_ptype==muplus || m_ptype==muminus)
1168  {
1169  //boost new neutrino to mu decay cm
1170  double beta[3];
1171  double p_nu[3]; //nu momentum
1172  beta[0]=m_pdPx / parent_energy;
1173  beta[1]=m_pdPy / parent_energy;
1174  beta[2]=m_pdPz / parent_energy;
1175 
1176  p_nu[0] = (xdet[0]- m_Vx) * nu_energy / rad;
1177  p_nu[1] = (xdet[1]- m_Vy) * nu_energy / rad;
1178  p_nu[2] = (xdet[2]- m_Vz) * nu_energy / rad;
1179 
1180  double partial = gamma*(beta[0]*p_nu[0]+
1181  beta[1]*p_nu[1]+
1182  beta[2]*p_nu[2]);
1183  partial = nu_energy-partial / (gamma+1.);
1184  double p_dcm_nu[4];
1185  for (int i=0;i<3;i++) p_dcm_nu[i]=p_nu[i]-beta[i]*gamma*partial;
1186  p_dcm_nu[3]=0.;
1187  for (int i=0;i<3;i++) p_dcm_nu[3]+=p_dcm_nu[i]*p_dcm_nu[i];
1188  p_dcm_nu[3]=sqrt(p_dcm_nu[3]);
1189 
1190  //boost parent of mu to mu production cm
1191  gamma=m_ppenergy / m_parent_mass;
1192  beta[0] = m_ppdxdz * m_pppz / m_ppenergy;
1193  beta[1] = m_ppdydz * m_pppz / m_ppenergy;
1194  beta[2] = m_pppz / m_ppenergy;
1195  partial = gamma*(beta[0]*m_muparpx+
1196  beta[1]*m_muparpy+
1197  beta[2]*m_muparpz);
1198  partial = m_mupare - partial / (gamma+1.);
1199  double p_pcm_mp[4];
1200  p_pcm_mp[0]=m_muparpx-beta[0]*gamma*partial;
1201  p_pcm_mp[1]=m_muparpy-beta[1]*gamma*partial;
1202  p_pcm_mp[2]=m_muparpz-beta[2]*gamma*partial;
1203  p_pcm_mp[3]=0.;
1204  for (int i=0;i<3;i++) p_pcm_mp[3]+=p_pcm_mp[i]*p_pcm_mp[i];
1205  p_pcm_mp[3]=sqrt(p_pcm_mp[3]);
1206 
1207  double wt_ratio = 1.;
1208  //have to check p_pcm_mp
1209  //it can be 0 if mupar..=0. (I guess muons created in target??)
1210  if (p_pcm_mp[3] != 0. ) {
1211  //calc new decay angle w.r.t. (anti)spin direction
1212  double costh = (p_dcm_nu[0]*p_pcm_mp[0]+
1213  p_dcm_nu[1]*p_pcm_mp[1]+
1214  p_dcm_nu[2]*p_pcm_mp[2])/(p_dcm_nu[3]*p_pcm_mp[3]);
1215 
1216  if (costh>1.) costh = 1.;
1217  else if (costh<-1.) costh = -1.;
1218 
1219  //calc relative weight due to angle difference
1220  if (m_Ntype == nue || m_Ntype == nuebar)
1221  {
1222  wt_ratio = 1.-costh;
1223  }
1224  else if (m_Ntype == numu || m_Ntype == numubar)
1225  {
1226  double xnu = 2.* enuzr / mumass;
1227  wt_ratio = ( (3.-2.*xnu) - (1.-2.*xnu)*costh ) / (3.-2.*xnu);
1228  }
1229  else if (m_Ntype == nutau || m_Ntype == nutaubar)
1230  {
1231  double xnu = 2.* enuzr / taumass;
1232  wt_ratio = ( (3.-2.*xnu) - (1.-2.*xnu)*costh ) / (3.-2.*xnu);
1233  std::cout << "calculating weight for tau neutrino; this may not be correct" << std::endl;
1234  }
1235  else
1236  {
1237  std::cout << "eventRates:: Bad neutrino type = " << m_Ntype << std::endl;
1238  }
1239  }
1240  nu_wght *= wt_ratio;
1241  }
1242 
1243 // if(iread > 60000)
1244 // std::cout << " end iread = " << iread << std::endl;
1245 
1246  return nu_wght;
1247 }
1248 
1249 double eventRates::GetXSec( double nu_type,
1250  double nu_energy,
1251  std::string current)
1252 
1253 {
1254 
1255  if( current != "NC" && current != "CC") {
1256  cout <<" eventRates::GetXSec: Current other than NC or CC specified... I don't know what to do." << endl;
1257  return -999;
1258  }
1259 
1260  int file_index=0;
1261  if (nu_type == 52) file_index = 0; // nuebar
1262  if (nu_type == 53) file_index = 1; // nue
1263  if (nu_type == 55) file_index = 2; // numubar
1264  if (nu_type == 56) file_index = 3; // numu
1265  if (nu_type == 58) file_index = 4; // nutaubar
1266  if (nu_type == 59) file_index = 5; // nutau
1267 
1268  int current_index = 0;
1269  if( current == "NC")
1270  current_index = 1;
1271 
1272  // calculate cross section
1273  double thexsec = 0.;
1274 
1275  // scale factor
1276  double scale_factor = 6.026e-10;
1277  // xseccc = xseccc * 4.09e9
1278  // 10**-38 cm2 * 10**6 kg/kton * 3.8e20 POT/year *1iron/56/1.66e-27kg
1279  // * 10**-4 m2/cm2 = 4.09e9 conversion factor
1280  // since table was 26*sigma(nu-proton) + 30*sigma(nu-neutron)
1281  //
1282  // 10**-38 cm2 * 10**6 kg/kton *1iron/56/1.66e-27kg
1283  // * 10**-4 m2/cm2 = 1.076e-11 conversion factor
1284  //
1285  // 10**-38 cm2 * 10**6 kg/kton * 1 nucleon / 1.66e-27 kg
1286  // * 10**-4 m2/cm2 = 6.026e-10
1287 
1288  // if energy is higher than any available xsection point
1289  // return xsection of highest available energy point
1290  if ( nu_energy > f_e_arr[fnlines-1][file_index][current_index] ) {
1291  thexsec = f_xsec_arr[fnlines-1][file_index][current_index]*f_e_arr[fnlines-1][file_index][current_index]*scale_factor;
1292  }
1293  // if energy is lower than any available xsection point
1294  // return xsection of lowest available energy poitn
1295  else if ( nu_energy < f_e_arr[0][file_index][current_index] ) {
1296  thexsec = f_xsec_arr[0][file_index][current_index]*f_e_arr[0][file_index][current_index]*scale_factor;
1297  }
1298  else {
1299  // if not, find the xsections for energy values immendiately above and
1300  // below the requested energy
1301  int energy_index = 0;
1302  for(int i = 0; i< fnlines-1; i++) {
1303  if( nu_energy > f_e_arr[i][file_index][current_index] &&
1304  nu_energy < f_e_arr[i+1][file_index][current_index]) {
1305  energy_index = i;
1306  break;
1307  }
1308  }
1309 
1310  double sig1 = f_xsec_arr[energy_index][file_index][current_index];
1311  double sig2 = f_xsec_arr[energy_index+1][file_index][current_index];
1312 
1313 
1314  thexsec = sig1 + ((sig2 - sig1)/(f_e_arr[energy_index+1][file_index][current_index]-f_e_arr[energy_index][file_index][current_index]))*(nu_energy - f_e_arr[energy_index][file_index][current_index]);
1315  thexsec = thexsec * nu_energy;
1316 
1317  thexsec = thexsec * scale_factor;
1318  }
1319 
1320  return thexsec;
1321 }
1322 
1323 
1324 void eventRates::ReadXSecsFromFiles() {
1325 
1326  fnbins = 1500;
1327 
1328  std::string base = getenv("G4LBNEWORKDIR");
1329  base = base + "/ProductionScripts/data/argon_genie2.8.4/";
1330 
1331  for(int current = 1; current <= 2; current++) {
1332 
1333  string charge = "cc";
1334 
1335  if ( current == 2 ) {
1336  charge = "nc";
1337  }
1338 
1339  const int narr=6;
1340  string suffix[narr];
1341  if (current == 1) {
1342  suffix[0] = "_nuebar.dat";
1343  suffix[1] = "_nue.dat";
1344  suffix[2] = "_numubar.dat";
1345  suffix[3] = "_numu.dat";
1346  suffix[4] = "_nutaubar.dat";
1347  suffix[5] = "_nutau.dat";
1348  } else if (current==2) {
1349  suffix[0] = "_nuebar.dat";
1350  suffix[1] = "_nue.dat";
1351  suffix[2] = "_numubar.dat";
1352  suffix[3] = "_numu.dat";
1353  suffix[4] = "_nutaubar.dat";
1354  suffix[5] = "_nutau.dat";
1355  }
1356 
1357  string filename[narr];
1358  for (int flavor=0; flavor<narr; flavor++) {
1359  filename[flavor] = "xsec_"+charge+suffix[flavor];
1360  string tmpfilename = base + filename[flavor];
1361 
1362  fdat_file[flavor].open(tmpfilename.c_str());
1363  if (fdat_file[flavor].fail()) {
1364  cout << " File not found: " << filename[flavor] << endl;
1365  //assert(0);
1366  }
1367  else {
1368  cout << " Opened "<<filename[flavor] << "."<<endl;
1369  }
1370  double row[2];
1371  fnlines = 0;
1372 
1373  while ( fdat_file[flavor] >> row[0] >> row[1] ) {
1374 
1375  if (fnlines >= fnbins) {
1376  cout << " length of data file exceed array size. Fix me. " << filename << endl;
1377  assert(0);
1378  }
1379 
1380  f_e_arr[fnlines][flavor][current-1] = row[0];
1381  f_xsec_arr[fnlines][flavor][current-1] = row[1];
1382  fnlines++;
1383  }
1384  fdat_file[flavor].close();
1385  }
1386  }
1387  std::cout<<"finished reading cross section files"<<std::endl;
1388 }
1389 
1390 double eventRates::GetOscillatedNeutrinoType(double ntype, double E) {
1391 
1392 
1393  int flavbefore = 0;
1394  const int nue = 53;
1395  const int nuebar = 52;
1396  const int numu = 56;
1397  const int numubar = 55;
1398  const int nutau = 59;
1399  const int nutaubar = 58;
1400  if(isDk2nu)
1401  Ntype = dk2nu->decay.ntype;
1402  else if (ntype == nuebar) flavbefore = -12;
1403  else if (ntype == numubar) flavbefore = -14;
1404  else if (ntype == nutaubar) flavbefore = -16;
1405  else if (ntype == nue) flavbefore = 12;
1406  else if (ntype == numu) flavbefore = 14;
1407  else if (ntype == nutau) flavbefore = 16;
1408  else {
1409  std::cout<<"WARNING: unrecognized neutrino type: "<<ntype<<std::endl;
1410  }
1411 
1412  osc::OscCalculator calculator;
1413  calculator.SetRho(2.8); // from LBNE Snowmass Document
1414  calculator.SetL(sqrt(detx*detx+dety*dety+detz*detz)/100/1000); // cm to km
1415  calculator.SetDmsq21(7.50e-5); // PDG 2013
1416  calculator.SetDmsq32(0.00232); // PDG 2013
1417  calculator.SetTh12(0.591); // PDG 2013
1418  calculator.SetTh13(.16); // PDG 2013
1419  calculator.SetTh23(0.785); // Maximal
1420  calculator.SetdCP(0);
1421 
1422  double P_e = 0;
1423  double P_t = 0;
1424  double P_m = 0;
1425 
1426  if(flavbefore<0) {
1427  P_e = calculator.P(flavbefore, -12, E);
1428  P_m = calculator.P(flavbefore, -14, E);
1429  P_t = calculator.P(flavbefore, -16, E);
1430  }
1431  else {
1432  P_e = calculator.P(flavbefore, 12, E);
1433  P_m = calculator.P(flavbefore, 14, E);
1434  P_t = calculator.P(flavbefore, 16, E);
1435  }
1436 
1437  // throw a random number to determine which neutrino is produced
1438  double random_number = rand3->Rndm();
1439 
1440  int ntype_oscillated;
1441  if(flavbefore <0) {
1442  if(random_number < P_e)
1443  ntype_oscillated = nuebar;
1444  else if(random_number < P_e + P_t)
1445  ntype_oscillated = nutaubar;
1446  else
1447  ntype_oscillated = numubar;
1448  }
1449  else {
1450  if(random_number < P_e)
1451  ntype_oscillated = nue;
1452  else if(random_number < P_e + P_t)
1453  ntype_oscillated = nutau;
1454  else
1455  ntype_oscillated = numu;
1456  }
1457 
1458  return ntype_oscillated;
1459 
1460 }
1461 
eventRates::GetXSec
double GetXSec(double nu_type, double nu_energy, std::string current)
Definition: eventRatesL.cc:1249
osc::OscCalculator::SetTh12
void SetTh12(double th12)
Definition: OscCalculator.h:43
keras_to_tensorflow.f
f
Definition: keras_to_tensorflow.py:162
string
std::string string
Definition: nybbler.cc:12
osc::OscCalculator::SetTh23
void SetTh23(double th23)
Definition: OscCalculator.h:45
std
STL namespace.
rdet
const double rdet
Definition: LBNENuWeight.cc:13
train.filename
string filename
Definition: train.py:213
sensitivity_comps.i
int i
Definition: sensitivity_comps.py:15
osc::OscCalculator::P
double P(int flavBefore, int flavAfter, double E)
Definition: OscCalculator.cxx:42
eventRates::Loop
virtual void Loop()
Definition: eventRatesL.cc:29
osc::OscCalculator::SetRho
void SetRho(double rho)
Definition: OscCalculator.h:40
eventRates::ReadXSecsFromFiles
void ReadXSecsFromFiles()
Definition: eventRatesL.cc:1324
fail
const HLTPathStatus fail
Definition: EventSelExc_t.cpp:91
cet::getenv
std::string getenv(std::string const &name)
Definition: getenv.cc:15
osc::OscCalculator::SetDmsq21
void SetDmsq21(double dmsq21)
Definition: OscCalculator.h:41
eventRates::SetTitles
void SetTitles(TH1 *h, const std::string &xtitle="", const std::string &ytitle="")
Definition: eventRatesL.cc:978
current
static Entry * current
Definition: commentscan.cpp:3516
iread
int iread
Definition: eventRatesL.cc:27
eventRates::GetOscillatedNeutrinoType
double GetOscillatedNeutrinoType(double ntype, double E)
Definition: eventRatesL.cc:1390
eventRates::GetPOTAsString
std::string GetPOTAsString(const double dpot)
Definition: eventRatesL.cc:844
eventRates::GetWeight
double GetWeight(const std::vector< double > xdet, double &nu_wght, double &nu_energy)
Definition: eventRatesL.cc:993
osc::OscCalculator::SetdCP
void SetdCP(double dCP)
Definition: OscCalculator.h:46
osc::OscCalculator::SetDmsq32
void SetDmsq32(double dmsq32)
Definition: OscCalculator.h:42
osc::OscCalculator::SetTh13
void SetTh13(double th13)
Definition: OscCalculator.h:44
osc::OscCalculator::SetL
void SetL(double L)
Definition: OscCalculator.h:39
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
train_cnn.h
h
training ###############################
Definition: train_cnn.py:186
Generated by   doxygen 1.8.11