HighLAND
FluxWeighting.cxx
1 #include "FluxWeighting.hxx"
2 #include "ND280AnalysisUtils.hxx"
3 
4 //********************************************************************
5 FluxWeighting::FluxWeighting(std::string const& fluxfolder,
6  std::string const& version,
7  std::string const& tuning,
8  std::string const& planeident) {
9  //********************************************************************
10 
11  bool load_run5 = false;
12  bool load_run6 = false;
13  if (fluxfolder.find("13a") != std::string::npos) {
14  load_run5 = true;
15  if (fluxfolder.find("av1.1") != std::string::npos) load_run6 = true;
16  }
17 
18  std::cout << " ------------------------------------------------------ "
19  << std::endl;
20  std::cout << " ------------------------------------------------------ "
21  << std::endl;
22  std::cout << " Flux reweight is enabled " << std::endl;
23  std::cout << " Folder : " << fluxfolder << std::endl;
24  std::cout << " File : "
25  << fluxfolder + "/runXXX/" + planeident + "_" + version +
26  "_runXXX.root"
27  << std::endl;
28  std::cout << " Tuning : " << tuning << std::endl;
29  std::cout << " ------------------------------------------------------ "
30  << std::endl;
31  std::cout << " ------------------------------------------------------ "
32  << std::endl;
33 
34  std::string file1 = "/run1/" + planeident + "_" + version + "_run1.root";
35  std::string file2 = "/run2/" + planeident + "_" + version + "_run2.root";
36  std::string file3b = "/run3b/" + planeident + "_" + version + "_run3b.root";
37  std::string file3c = "/run3c/" + planeident + "_" + version + "_run3c.root";
38  std::string file4 = "/run4/" + planeident + "_" + version + "_run4.root";
39  std::string file5a = "/run5a/" + planeident + "_" + version + "_run5a.root";
40  std::string file5b = "/run5b/" + planeident + "_" + version + "_run5b.root";
41  std::string file5c =
42  "/run5c/" + planeident + "_" + version + "_run5c_antinumode.root";
43  std::string file6b =
44  "/run6b/" + planeident + "_" + version + "_run6b_antinumode.root";
45  std::string file6c =
46  "/run6c/" + planeident + "_" + version + "_run6c_antinumode.root";
47  std::string file6d =
48  "/run6d/" + planeident + "_" + version + "_run6d_antinumode.root";
49  std::string file6e =
50  "/run6e/" + planeident + "_" + version + "_run6e_antinumode.root";
51 
52  run1_file = new TFile(file1.insert(0, fluxfolder).c_str(), "READ");
53  run2_file = new TFile(file2.insert(0, fluxfolder).c_str(), "READ");
54  run3b_file = new TFile(file3b.insert(0, fluxfolder).c_str(), "READ");
55  run3c_file = new TFile(file3c.insert(0, fluxfolder).c_str(), "READ");
56  run4_file = new TFile(file4.insert(0, fluxfolder).c_str(), "READ");
57 
58  if (load_run5) {
59  run5a_file = new TFile(file5a.insert(0, fluxfolder).c_str(), "READ");
60  run5b_file = new TFile(file5b.insert(0, fluxfolder).c_str(), "READ");
61  run5c_file = new TFile(file5c.insert(0, fluxfolder).c_str(), "READ");
62  } else {
63  run5a_file = NULL;
64  run5b_file = NULL;
65  run5c_file = NULL;
66  }
67 
68  if (load_run6) {
69  run6b_file = new TFile(file6b.insert(0, fluxfolder).c_str(), "READ");
70  run6c_file = new TFile(file6c.insert(0, fluxfolder).c_str(), "READ");
71  run6d_file = new TFile(file6d.insert(0, fluxfolder).c_str(), "READ");
72  run6e_file = new TFile(file6e.insert(0, fluxfolder).c_str(), "READ");
73  } else {
74  run6b_file = NULL;
75  run6c_file = NULL;
76  run6d_file = NULL;
77  run6e_file = NULL;
78  }
79 
80  // Get the histograms
81  run1_weight_numu = dynamic_cast<TH1D*>(run1_file->Get(
82  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
83  if (run1_weight_numu) {
84  run1_weight_numu = static_cast<TH1D*>(run1_weight_numu->Clone());
85  run1_weight_numu->SetDirectory(NULL);
86  }
87  run1_weight_numubar = dynamic_cast<TH1D*>(run1_file->Get(
88  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
89  if (run1_weight_numubar) {
90  run1_weight_numubar = static_cast<TH1D*>(run1_weight_numubar->Clone());
91  run1_weight_numubar->SetDirectory(NULL);
92  }
93  run1_weight_nue = dynamic_cast<TH1D*>(run1_file->Get(
94  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
95  if (run1_weight_nue) {
96  run1_weight_nue = static_cast<TH1D*>(run1_weight_nue->Clone());
97  run1_weight_nue->SetDirectory(NULL);
98  }
99  run1_weight_nuebar = dynamic_cast<TH1D*>(run1_file->Get(
100  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
101  if (run1_weight_nuebar) {
102  run1_weight_nuebar = static_cast<TH1D*>(run1_weight_nuebar->Clone());
103  run1_weight_nuebar->SetDirectory(NULL);
104  }
105 
106  run2_weight_numu = dynamic_cast<TH1D*>(run2_file->Get(
107  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
108  if (run2_weight_numu) {
109  run2_weight_numu = static_cast<TH1D*>(run2_weight_numu->Clone());
110  run2_weight_numu->SetDirectory(NULL);
111  }
112  run2_weight_numubar = dynamic_cast<TH1D*>(run2_file->Get(
113  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
114  if (run2_weight_numubar) {
115  run2_weight_numubar = static_cast<TH1D*>(run2_weight_numubar->Clone());
116  run2_weight_numubar->SetDirectory(NULL);
117  }
118  run2_weight_nue = dynamic_cast<TH1D*>(run2_file->Get(
119  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
120  if (run2_weight_nue) {
121  run2_weight_nue = static_cast<TH1D*>(run2_weight_nue->Clone());
122  run2_weight_nue->SetDirectory(NULL);
123  }
124  run2_weight_nuebar = dynamic_cast<TH1D*>(run2_file->Get(
125  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
126  if (run2_weight_nuebar) {
127  run2_weight_nuebar = static_cast<TH1D*>(run2_weight_nuebar->Clone());
128  run2_weight_nuebar->SetDirectory(NULL);
129  }
130 
131  run3b_weight_numu = dynamic_cast<TH1D*>(run3b_file->Get(
132  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
133  if (run3b_weight_numu) {
134  run3b_weight_numu = static_cast<TH1D*>(run3b_weight_numu->Clone());
135  run3b_weight_numu->SetDirectory(NULL);
136  }
137  run3b_weight_numubar = dynamic_cast<TH1D*>(run3b_file->Get(
138  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
139  if (run3b_weight_numubar) {
140  run3b_weight_numubar = static_cast<TH1D*>(run3b_weight_numubar->Clone());
141  run3b_weight_numubar->SetDirectory(NULL);
142  }
143  run3b_weight_nue = dynamic_cast<TH1D*>(run3b_file->Get(
144  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
145  if (run3b_weight_nue) {
146  run3b_weight_nue = static_cast<TH1D*>(run3b_weight_nue->Clone());
147  run3b_weight_nue->SetDirectory(NULL);
148  }
149  run3b_weight_nuebar = dynamic_cast<TH1D*>(run3b_file->Get(
150  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
151  if (run3b_weight_nuebar) {
152  run3b_weight_nuebar = static_cast<TH1D*>(run3b_weight_nuebar->Clone());
153  run3b_weight_nuebar->SetDirectory(NULL);
154  }
155 
156  run3c_weight_numu = dynamic_cast<TH1D*>(run3c_file->Get(
157  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
158  if (run3c_weight_numu) {
159  run3c_weight_numu = static_cast<TH1D*>(run3c_weight_numu->Clone());
160  run3c_weight_numu->SetDirectory(NULL);
161  }
162  run3c_weight_numubar = dynamic_cast<TH1D*>(run3c_file->Get(
163  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
164  if (run3c_weight_numubar) {
165  run3c_weight_numubar = static_cast<TH1D*>(run3c_weight_numubar->Clone());
166  run3c_weight_numubar->SetDirectory(NULL);
167  }
168  run3c_weight_nue = dynamic_cast<TH1D*>(run3c_file->Get(
169  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
170  if (run3c_weight_nue) {
171  run3c_weight_nue = static_cast<TH1D*>(run3c_weight_nue->Clone());
172  run3c_weight_nue->SetDirectory(NULL);
173  }
174  run3c_weight_nuebar = dynamic_cast<TH1D*>(run3c_file->Get(
175  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
176  if (run3c_weight_nuebar) {
177  run3c_weight_nuebar = static_cast<TH1D*>(run3c_weight_nuebar->Clone());
178  run3c_weight_nuebar->SetDirectory(NULL);
179  }
180 
181  run4_weight_numu = dynamic_cast<TH1D*>(run4_file->Get(
182  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
183  if (run4_weight_numu) {
184  run4_weight_numu = static_cast<TH1D*>(run4_weight_numu->Clone());
185  run4_weight_numu->SetDirectory(NULL);
186  }
187  run4_weight_numubar = dynamic_cast<TH1D*>(run4_file->Get(
188  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
189  if (run4_weight_numubar) {
190  run4_weight_numubar = static_cast<TH1D*>(run4_weight_numubar->Clone());
191  run4_weight_numubar->SetDirectory(NULL);
192  }
193  run4_weight_nue = dynamic_cast<TH1D*>(run4_file->Get(
194  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
195  if (run4_weight_nue) {
196  run4_weight_nue = static_cast<TH1D*>(run4_weight_nue->Clone());
197  run4_weight_nue->SetDirectory(NULL);
198  }
199  run4_weight_nuebar = dynamic_cast<TH1D*>(run4_file->Get(
200  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
201  if (run4_weight_nuebar) {
202  run4_weight_nuebar = static_cast<TH1D*>(run4_weight_nuebar->Clone());
203  run4_weight_nuebar->SetDirectory(NULL);
204  }
205 
206  if (load_run5) {
207  run5a_weight_numu = dynamic_cast<TH1D*>(run5a_file->Get(
208  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
209  if (run5a_weight_numu) {
210  run5a_weight_numu = static_cast<TH1D*>(run5a_weight_numu->Clone());
211  run5a_weight_numu->SetDirectory(NULL);
212  }
213  run5a_weight_numubar = dynamic_cast<TH1D*>(run5a_file->Get(
214  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
215  if (run5a_weight_numubar) {
216  run5a_weight_numubar = static_cast<TH1D*>(run5a_weight_numubar->Clone());
217  run5a_weight_numubar->SetDirectory(NULL);
218  }
219  run5a_weight_nue = dynamic_cast<TH1D*>(run5a_file->Get(
220  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
221  if (run5a_weight_nue) {
222  run5a_weight_nue = static_cast<TH1D*>(run5a_weight_nue->Clone());
223  run5a_weight_nue->SetDirectory(NULL);
224  }
225  run5a_weight_nuebar = dynamic_cast<TH1D*>(run5a_file->Get(
226  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
227  if (run5a_weight_nuebar) {
228  run5a_weight_nuebar = static_cast<TH1D*>(run5a_weight_nuebar->Clone());
229  run5a_weight_nuebar->SetDirectory(NULL);
230  }
231 
232  run5b_weight_numu = dynamic_cast<TH1D*>(run5b_file->Get(
233  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
234  if (run5b_weight_numu) {
235  run5b_weight_numu = static_cast<TH1D*>(run5b_weight_numu->Clone());
236  run5b_weight_numu->SetDirectory(NULL);
237  }
238  run5b_weight_numubar = dynamic_cast<TH1D*>(run5b_file->Get(
239  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
240  if (run5b_weight_numubar) {
241  run5b_weight_numubar = static_cast<TH1D*>(run5b_weight_numubar->Clone());
242  run5b_weight_numubar->SetDirectory(NULL);
243  }
244  run5b_weight_nue = dynamic_cast<TH1D*>(run5b_file->Get(
245  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
246  if (run5b_weight_nue) {
247  run5b_weight_nue = static_cast<TH1D*>(run5b_weight_nue->Clone());
248  run5b_weight_nue->SetDirectory(NULL);
249  }
250  run5b_weight_nuebar = dynamic_cast<TH1D*>(run5b_file->Get(
251  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
252  if (run5b_weight_nuebar) {
253  run5b_weight_nuebar = static_cast<TH1D*>(run5b_weight_nuebar->Clone());
254  run5b_weight_nuebar->SetDirectory(NULL);
255  }
256 
257  run5c_weight_numu = dynamic_cast<TH1D*>(run5c_file->Get(
258  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
259  if (run5c_weight_numu) {
260  run5c_weight_numu = static_cast<TH1D*>(run5c_weight_numu->Clone());
261  run5c_weight_numu->SetDirectory(NULL);
262  }
263  run5c_weight_numubar = dynamic_cast<TH1D*>(run5c_file->Get(
264  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
265  if (run5c_weight_numubar) {
266  run5c_weight_numubar = static_cast<TH1D*>(run5c_weight_numubar->Clone());
267  run5c_weight_numubar->SetDirectory(NULL);
268  }
269  run5c_weight_nue = dynamic_cast<TH1D*>(run5c_file->Get(
270  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
271  if (run5c_weight_nue) {
272  run5c_weight_nue = static_cast<TH1D*>(run5c_weight_nue->Clone());
273  run5c_weight_nue->SetDirectory(NULL);
274  }
275  run5c_weight_nuebar = dynamic_cast<TH1D*>(run5c_file->Get(
276  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
277  if (run5c_weight_nuebar) {
278  run5c_weight_nuebar = static_cast<TH1D*>(run5c_weight_nuebar->Clone());
279  run5c_weight_nuebar->SetDirectory(NULL);
280  }
281  } else {
282  run5a_weight_numu = NULL;
283  run5a_weight_numubar = NULL;
284  run5a_weight_nue = NULL;
285  run5a_weight_nuebar = NULL;
286 
287  run5b_weight_numu = NULL;
288  run5b_weight_numubar = NULL;
289  run5b_weight_nue = NULL;
290  run5b_weight_nuebar = NULL;
291 
292  run5c_weight_numu = NULL;
293  run5c_weight_numubar = NULL;
294  run5c_weight_nue = NULL;
295  run5c_weight_nuebar = NULL;
296  }
297  if (load_run6) {
298  run6b_weight_numu = dynamic_cast<TH1D*>(run6b_file->Get(
299  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
300  if (run6b_weight_numu) {
301  run6b_weight_numu = static_cast<TH1D*>(run6b_weight_numu->Clone());
302  run6b_weight_numu->SetDirectory(NULL);
303  }
304  run6b_weight_numubar = dynamic_cast<TH1D*>(run6b_file->Get(
305  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
306  if (run6b_weight_numubar) {
307  run6b_weight_numubar = static_cast<TH1D*>(run6b_weight_numubar->Clone());
308  run6b_weight_numubar->SetDirectory(NULL);
309  }
310  run6b_weight_nue = dynamic_cast<TH1D*>(run6b_file->Get(
311  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
312  if (run6b_weight_nue) {
313  run6b_weight_nue = static_cast<TH1D*>(run6b_weight_nue->Clone());
314  run6b_weight_nue->SetDirectory(NULL);
315  }
316  run6b_weight_nuebar = dynamic_cast<TH1D*>(run6b_file->Get(
317  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
318  if (run6b_weight_nuebar) {
319  run6b_weight_nuebar = static_cast<TH1D*>(run6b_weight_nuebar->Clone());
320  run6b_weight_nuebar->SetDirectory(NULL);
321  }
322  run6c_weight_numu = dynamic_cast<TH1D*>(run6c_file->Get(
323  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
324  if (run6c_weight_numu) {
325  run6c_weight_numu = static_cast<TH1D*>(run6c_weight_numu->Clone());
326  run6c_weight_numu->SetDirectory(NULL);
327  }
328  run6c_weight_numubar = dynamic_cast<TH1D*>(run6c_file->Get(
329  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
330  if (run6c_weight_numubar) {
331  run6c_weight_numubar = static_cast<TH1D*>(run6c_weight_numubar->Clone());
332  run6c_weight_numubar->SetDirectory(NULL);
333  }
334  run6c_weight_nue = dynamic_cast<TH1D*>(run6c_file->Get(
335  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
336  if (run6c_weight_nue) {
337  run6c_weight_nue = static_cast<TH1D*>(run6c_weight_nue->Clone());
338  run6c_weight_nue->SetDirectory(NULL);
339  }
340  run6c_weight_nuebar = dynamic_cast<TH1D*>(run6c_file->Get(
341  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
342  if (run6c_weight_nuebar) {
343  run6c_weight_nuebar = static_cast<TH1D*>(run6c_weight_nuebar->Clone());
344  run6c_weight_nuebar->SetDirectory(NULL);
345  }
346  run6d_weight_numu = dynamic_cast<TH1D*>(run6d_file->Get(
347  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
348  if (run6d_weight_numu) {
349  run6d_weight_numu = static_cast<TH1D*>(run6d_weight_numu->Clone());
350  run6d_weight_numu->SetDirectory(NULL);
351  }
352  run6d_weight_numubar = dynamic_cast<TH1D*>(run6d_file->Get(
353  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
354  if (run6d_weight_numubar) {
355  run6d_weight_numubar = static_cast<TH1D*>(run6d_weight_numubar->Clone());
356  run6d_weight_numubar->SetDirectory(NULL);
357  }
358  run6d_weight_nue = dynamic_cast<TH1D*>(run6d_file->Get(
359  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
360  if (run6d_weight_nue) {
361  run6d_weight_nue = static_cast<TH1D*>(run6d_weight_nue->Clone());
362  run6d_weight_nue->SetDirectory(NULL);
363  }
364  run6d_weight_nuebar = dynamic_cast<TH1D*>(run6d_file->Get(
365  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
366  if (run6d_weight_nuebar) {
367  run6d_weight_nuebar = static_cast<TH1D*>(run6d_weight_nuebar->Clone());
368  run6d_weight_nuebar->SetDirectory(NULL);
369  }
370  run6e_weight_numu = dynamic_cast<TH1D*>(run6e_file->Get(
371  ("enu_" + planeident + "_" + tuning + "_numu_ratio").c_str()));
372  if (run6e_weight_numu) {
373  run6e_weight_numu = static_cast<TH1D*>(run6e_weight_numu->Clone());
374  run6e_weight_numu->SetDirectory(NULL);
375  }
376  run6e_weight_numubar = dynamic_cast<TH1D*>(run6e_file->Get(
377  ("enu_" + planeident + "_" + tuning + "_numub_ratio").c_str()));
378  if (run6e_weight_numubar) {
379  run6e_weight_numubar = static_cast<TH1D*>(run6e_weight_numubar->Clone());
380  run6e_weight_numubar->SetDirectory(NULL);
381  }
382  run6e_weight_nue = dynamic_cast<TH1D*>(run6e_file->Get(
383  ("enu_" + planeident + "_" + tuning + "_nue_ratio").c_str()));
384  if (run6e_weight_nue) {
385  run6e_weight_nue = static_cast<TH1D*>(run6e_weight_nue->Clone());
386  run6e_weight_nue->SetDirectory(NULL);
387  }
388  run6e_weight_nuebar = dynamic_cast<TH1D*>(run6e_file->Get(
389  ("enu_" + planeident + "_" + tuning + "_nueb_ratio").c_str()));
390  if (run6e_weight_nuebar) {
391  run6e_weight_nuebar = static_cast<TH1D*>(run6e_weight_nuebar->Clone());
392  run6e_weight_nuebar->SetDirectory(NULL);
393  }
394 
395  } else {
396  run6b_weight_numu = NULL;
397  run6b_weight_numubar = NULL;
398  run6b_weight_nue = NULL;
399  run6b_weight_nuebar = NULL;
400  run6c_weight_numu = NULL;
401  run6c_weight_numubar = NULL;
402  run6c_weight_nue = NULL;
403  run6c_weight_nuebar = NULL;
404  run6d_weight_numu = NULL;
405  run6d_weight_numubar = NULL;
406  run6d_weight_nue = NULL;
407  run6d_weight_nuebar = NULL;
408  run6e_weight_numu = NULL;
409  run6e_weight_numubar = NULL;
410  run6e_weight_nue = NULL;
411  run6e_weight_nuebar = NULL;
412  }
413 
414  if (run1_file) {
415  if (run1_file) {
416  run1_file->Close();
417  }
418  delete run1_file;
419  if (run2_file) {
420  run2_file->Close();
421  }
422  delete run2_file;
423  if (run3b_file) {
424  run3b_file->Close();
425  }
426  delete run3b_file;
427  if (run3c_file) {
428  run3c_file->Close();
429  }
430  delete run3c_file;
431  if (run4_file) {
432  run4_file->Close();
433  }
434  delete run4_file;
435  if (run5a_file) {
436  if (run5a_file) {
437  run5a_file->Close();
438  }
439  delete run5a_file;
440  if (run5b_file) {
441  run5b_file->Close();
442  }
443  delete run5b_file;
444  if (run5c_file) {
445  run5c_file->Close();
446  }
447  delete run5c_file;
448  if (run6b_file) {
449  if (run6b_file) {
450  run6b_file->Close();
451  }
452  delete run6b_file;
453  if (run6c_file) {
454  run6c_file->Close();
455  }
456  delete run6c_file;
457  if (run6d_file) {
458  run6d_file->Close();
459  }
460  delete run6d_file;
461  if (run6e_file) {
462  run6e_file->Close();
463  }
464  delete run6e_file;
465  }
466  }
467  }
468 
469  if (!run1_weight_numu || !run1_weight_numubar || !run1_weight_nue ||
470  !run1_weight_nuebar) {
471  std::cerr
472  << "ERROR: Flux weighting was requested, but could not be initialised."
473  << std::endl;
474  std::cerr
475  << " Specify a valid flux file and tuning in your parameters file, or"
476  << std::endl;
477  std::cerr
478  << " disable the flux weighting. Flux files can be downloaded from"
479  << std::endl;
480  std::cerr << " http://www.t2k.org/beam/NuFlux/FluxRelease" << std::endl;
481  exit(EXIT_FAILURE);
482  }
483 }
484 
485 //********************************************************************
486 FluxWeighting::~FluxWeighting() {
487  //********************************************************************
488  if (run1_weight_numu) {
489  delete run1_weight_numu;
490  delete run2_weight_numu;
491  delete run3b_weight_numu;
492  delete run3c_weight_numu;
493  delete run4_weight_numu;
494  if (run5a_weight_numu) {
495  delete run5a_weight_numu;
496  delete run5b_weight_numu;
497  delete run5c_weight_numu;
498  if (run6b_weight_numu) {
499  delete run6b_weight_numu;
500  delete run6c_weight_numu;
501  delete run6d_weight_numu;
502  delete run6e_weight_numu;
503  }
504  }
505  }
506  if (run1_weight_numubar) {
507  delete run1_weight_numubar;
508  delete run2_weight_numubar;
509  delete run3b_weight_numubar;
510  delete run3c_weight_numubar;
511  delete run4_weight_numubar;
512  if (run5a_weight_numubar) {
513  delete run5a_weight_numubar;
514  delete run5b_weight_numubar;
515  delete run5c_weight_numubar;
516  if (run6b_weight_numubar) {
517  delete run6b_weight_numubar;
518  delete run6c_weight_numubar;
519  delete run6d_weight_numubar;
520  delete run6e_weight_numubar;
521  }
522  }
523  }
524  if (run1_weight_nue) {
525  delete run1_weight_nue;
526  delete run2_weight_nue;
527  delete run3b_weight_nue;
528  delete run3c_weight_nue;
529  delete run4_weight_nue;
530  if (run5a_weight_nue) {
531  delete run5a_weight_nue;
532  delete run5b_weight_nue;
533  delete run5c_weight_nue;
534  if (run6b_weight_nue) {
535  delete run6b_weight_nue;
536  delete run6c_weight_nue;
537  delete run6d_weight_nue;
538  delete run6e_weight_nue;
539  }
540  }
541  }
542  if (run1_weight_nuebar) {
543  delete run1_weight_nuebar;
544  delete run2_weight_nuebar;
545  delete run3b_weight_nuebar;
546  delete run3c_weight_nuebar;
547  delete run4_weight_nuebar;
548  if (run5a_weight_nuebar) {
549  delete run5a_weight_nuebar;
550  delete run5b_weight_nuebar;
551  delete run5c_weight_nuebar;
552  if (run6b_weight_nuebar) {
553  delete run6b_weight_nuebar;
554  delete run6c_weight_nuebar;
555  delete run6d_weight_nuebar;
556  delete run6e_weight_nuebar;
557  }
558  }
559  }
560 }
561 
562 //********************************************************************
563 void FluxWeighting::UpdateBunchWeight(AnaBunchB& bunch, AnaTrueVertexB* vertex,
564  int RunPeriod) {
565  //********************************************************************
566  Float_t weight = GetWeight(vertex, RunPeriod);
567  bunch.Weight *= weight;
568 }
569 
570 //********************************************************************
571 void FluxWeighting::UpdateEventWeight(AnaEventB& event,
572  AnaTrueVertexB* vertex) {
573  //********************************************************************
574  Float_t weight =
575  GetWeight(vertex, anaUtils::GetRunPeriod(event.EventInfo.Run));
576  event.Weight *= weight;
577 }
578 
579 //********************************************************************
580 Float_t FluxWeighting::GetWeight(AnaTrueVertexB* vertex, int RunPeriod) {
581  //********************************************************************
582  Float_t weight = 1.;
583 
584  if (vertex) {
585  weight = GetWeight(*vertex, RunPeriod);
586  }
587 
588  return weight;
589 }
590 
591 //********************************************************************
592 Float_t FluxWeighting::GetWeight(const AnaTrueVertexB& vertex, int RunPeriod) {
593  //********************************************************************
594 
595  Float_t weight = 1.;
596  TH1D* hist = NULL;
597 
598  switch (vertex.NuPDG) {
599  case 14:
600  switch (RunPeriod) {
601  case 0:
602  hist = run1_weight_numu;
603  break;
604  case 1:
605  hist = run2_weight_numu;
606  break;
607  case 2:
608  hist = run2_weight_numu;
609  break;
610  case 3:
611  hist = run3b_weight_numu;
612  break;
613  case 4:
614  hist = run3c_weight_numu;
615  break;
616  case 5:
617  hist = run4_weight_numu;
618  break;
619  case 6:
620  hist = run4_weight_numu;
621  break;
622  case 7:
623  hist = run5a_weight_numu;
624  break;
625  case 8:
626  hist = run5c_weight_numu;
627  break;
628  case 9:
629  hist = run6b_weight_numu;
630  break;
631  case 10:
632  hist = run6c_weight_numu;
633  break;
634  case 11:
635  hist = run6d_weight_numu;
636  break;
637  case 12:
638  hist = run6e_weight_numu;
639  break;
640  default:
641  std::cerr << "Unknown run number " << RunPeriod
642  << ": setting event weight to 1." << std::endl;
643  break;
644  }
645  break;
646  case -14:
647  switch (RunPeriod) {
648  case 0:
649  hist = run1_weight_numubar;
650  break;
651  case 1:
652  hist = run2_weight_numubar;
653  break;
654  case 2:
655  hist = run2_weight_numubar;
656  break;
657  case 3:
658  hist = run3b_weight_numubar;
659  break;
660  case 4:
661  hist = run3c_weight_numubar;
662  break;
663  case 5:
664  hist = run4_weight_numubar;
665  break;
666  case 6:
667  hist = run4_weight_numubar;
668  break;
669  case 7:
670  hist = run5a_weight_numubar;
671  break;
672  case 8:
673  hist = run5c_weight_numubar;
674  break;
675  case 9:
676  hist = run6b_weight_numubar;
677  break;
678  case 10:
679  hist = run6c_weight_numubar;
680  break;
681  case 11:
682  hist = run6d_weight_numubar;
683  break;
684  case 12:
685  hist = run6e_weight_numubar;
686  break;
687  default:
688  std::cerr << "Unknown run number " << RunPeriod
689  << ": setting event weight to 1." << std::endl;
690  break;
691  }
692  break;
693  case 12:
694  switch (RunPeriod) {
695  case 0:
696  hist = run1_weight_nue;
697  break;
698  case 1:
699  hist = run2_weight_nue;
700  break;
701  case 2:
702  hist = run2_weight_nue;
703  break;
704  case 3:
705  hist = run3b_weight_nue;
706  break;
707  case 4:
708  hist = run3c_weight_nue;
709  break;
710  case 5:
711  hist = run4_weight_nue;
712  break;
713  case 6:
714  hist = run4_weight_nue;
715  break;
716  case 7:
717  hist = run5a_weight_nue;
718  break;
719  case 8:
720  hist = run5c_weight_nue;
721  break;
722  case 9:
723  hist = run6b_weight_nue;
724  break;
725  case 10:
726  hist = run6c_weight_nue;
727  break;
728  case 11:
729  hist = run6d_weight_nue;
730  break;
731  case 12:
732  hist = run6e_weight_nue;
733  break;
734  default:
735  std::cerr << "Unknown run number " << RunPeriod
736  << ": setting event weight to 1." << std::endl;
737  break;
738  }
739  break;
740  case -12:
741  switch (RunPeriod) {
742  case 0:
743  hist = run1_weight_nuebar;
744  break;
745  case 1:
746  hist = run2_weight_nuebar;
747  break;
748  case 2:
749  hist = run2_weight_nuebar;
750  break;
751  case 3:
752  hist = run3b_weight_nuebar;
753  break;
754  case 4:
755  hist = run3c_weight_nuebar;
756  break;
757  case 5:
758  hist = run4_weight_nuebar;
759  break;
760  case 6:
761  hist = run4_weight_nuebar;
762  break;
763  case 7:
764  hist = run5a_weight_nuebar;
765  break;
766  case 8:
767  hist = run5c_weight_nuebar;
768  break;
769  case 9:
770  hist = run6b_weight_nuebar;
771  break;
772  case 10:
773  hist = run6c_weight_nuebar;
774  break;
775  case 11:
776  hist = run6d_weight_nuebar;
777  break;
778  case 12:
779  hist = run6e_weight_nuebar;
780  break;
781  default:
782  std::cerr << "Unknown run number " << RunPeriod
783  << ": setting event weight to 1." << std::endl;
784  break;
785  }
786  break;
787  default:
788  std::cerr << "Unknown neutrino flavour " << vertex.NuPDG
789  << ": setting event weight to 1." << std::endl;
790  break;
791  }
792 
793  if (hist) {
794  int bin = hist->FindBin(vertex.NuEnergy / 1000.);
795  weight = hist->GetBinContent(bin);
796  } else {
797  std::cerr << "Couldn't get correct flux weighting histogram - check you "
798  "have the right tuning folder specified in the parameters file"
799  << std::endl;
800  }
801 
802  return weight;
803 }
AnaTrueVertexB::NuPDG
Int_t NuPDG
The PDG code of the incoming neutrino.
Definition: BaseDataClasses.hxx:119
AnaTrueVertexB
Representation of a true Monte Carlo vertex.
Definition: BaseDataClasses.hxx:97
AnaEventB::EventInfo
AnaEventInfoB EventInfo
Run, sunrun, event, time stamp, etc.
Definition: BaseDataClasses.hxx:1123
FluxWeighting::FluxWeighting
FluxWeighting(std::string const &fluxfolder, std::string const &version, std::string const &tuning, std::string const &planeident="nd5")
Definition: FluxWeighting.cxx:5
AnaTrueVertexB::NuEnergy
Float_t NuEnergy
The true energy of the incoming neutrino.
Definition: BaseDataClasses.hxx:122
anaUtils::GetRunPeriod
int GetRunPeriod(int run, int subrun=-1)
Returns the run period (sequentially: 0,1,2,3,4,5 ...)
Definition: ND280AnalysisUtils.cxx:17
FluxWeighting::UpdateBunchWeight
void UpdateBunchWeight(AnaBunchB &bunch, AnaTrueVertexB *vertex, int RunPeriod)
Definition: FluxWeighting.cxx:563
FluxWeighting::GetWeight
Float_t GetWeight(AnaTrueVertexB *vertex, int RunPeriod)
Definition: FluxWeighting.cxx:580
AnaBunchC::Weight
Float_t Weight
The weight to apply to this bunch (nominally 1). An example is the beam flux weight.
Definition: CoreDataClasses.hxx:183
FluxWeighting::UpdateEventWeight
void UpdateEventWeight(AnaEventB &event, AnaTrueVertexB *vertex)
Definition: FluxWeighting.cxx:571
AnaEventInfoB::Run
Int_t Run
The ND280 run number.
Definition: BaseDataClasses.hxx:942
Generated by   doxygen 1.8.13