specprod_numuccinc.C
Go to the documentation of this file.
1 #include "CAFAna/Vars/Vars.h"
2 #include "CAFAna/Core/Binning.h"
3 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Systs/Systs.h"
12 
14 // #include "CAFAna/Vars/XsecTunes.h"
16 
22 
23 #include "TFile.h"
24 #include "TH1.h"
25 #include "TH2.h"
26 
27 #include <iostream>
28 #include <cmath>
29 #include <sstream>
30 #include <string>
31 #include <fstream>
32 #include <iomanip>
33 
34 #include "Utilities/func/MathUtil.h"
35 
36 using namespace ana;
37 
38 vector<string> dataset_lists = {
39  "/development/NDAna/numucc_inc/nersc/defnames/nd_fhc_ndrespin_caf_0-3_of_5.txt",
40  "/development/NDAna/numucc_inc/nersc/defnames/nd_fhc_ndrespin_caf_4_of_5.txt",
41  "/development/NDAna/numucc_inc/nersc/defnames/prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1.txt",
42  "/development/NDAna/numucc_inc/nersc/defnames/prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1.txt",
43  "/development/NDAna/numucc_inc/nersc/defnames/prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1.txt",
44  "/development/NDAna/numucc_inc/nersc/defnames/prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1.txt",
45  "/development/NDAna/numucc_inc/nersc/defnames/prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1.txt",
46  "/development/NDAna/numucc_inc/nersc/defnames/prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1.txt",
47 };
48 
49 vector<string> FileListFromDataset(int ds_idx)
50 {
51  vector<string> all_filenames;
52  ifstream inf(dataset_lists[ds_idx]);
53  string fn;
54  vector<string> chosen_filenames;
55 
56  while(inf >> fn)
57  if(ds_idx == 0 || ds_idx == 1)
58  all_filenames.push_back("/input/nominal/"+fn);
59  else if(ds_idx == 2)
60  all_filenames.push_back("/input/calibup/"+fn);
61  else if(ds_idx == 3)
62  all_filenames.push_back("/input/calibdown/"+fn);
63  else if(ds_idx == 4)
64  all_filenames.push_back("/input/cherenkov/"+fn);
65  else if(ds_idx == 5)
66  all_filenames.push_back("/input/neg_offset/"+fn);
67  else if(ds_idx == 6)
68  all_filenames.push_back("/input/pos_offset/"+fn);
69  else if(ds_idx == 7)
70  all_filenames.push_back("/input/shape/"+fn);
71  else
72  all_filenames.push_back("/input/"+fn);
73 
74  std::cout << "Debug ds_idx all_filenames.size " << ds_idx << " " << all_filenames.size() << std::endl;
75 
76  return all_filenames;
77 }
78 
79 // start_file_num and end_file_num are arguments for NERSC running.
80 void specprod_numuccinc(string fOutname, int ifirst_genie_univ, int ilast_genie_univ, int ifirst_ppfx_univ, int ilast_ppfx_univ, int ido_allother, int start_file_num = -1, int end_file_num = -1){
81 
82  std::cout<<"=> numucc_inc_xs() start "<<std::endl;
83 
84  ////////// Information /////////////////////////////////////////////////////
85 
86  int Nuniv_genie = 0;
87  if(ifirst_genie_univ > -1 && ilast_genie_univ > -1){
88  Nuniv_genie = ilast_genie_univ - ifirst_genie_univ + 1;
89  if(Nuniv_genie<0)Nuniv_genie = 0;
90  }
91  int Nuniv_ppfx = 0;
92  if(ifirst_ppfx_univ > -1 && ilast_ppfx_univ > -1){
93  Nuniv_ppfx = ilast_ppfx_univ - ifirst_ppfx_univ + 1;
94  if(Nuniv_ppfx<0)Nuniv_ppfx = 0;
95  }
96  bool do_allother;
97  // 1 CV
98  // 6 det systs
99  // 4 Muon energy scale: 2 shifting Emu and Enu (for checking), and 2 only shifting Emu
100  // 22 Foc: 2 Horn Current, 2 Beam Spot Size, 4 Beam Shift XY, 8 Beam Horn Pos XY and H1 and 2,
101  // 2 Target Z, 2 Geo Mag Field, 2 IC Water Layer
102  // Tot 33
103  int Nshift_other = 0;
104  if(ido_allother == 0){
105  do_allother = false;
106  Nshift_other = 0;
107  }
108  else if(ido_allother == 1){
109  do_allother = true;
110  Nshift_other = 33;
111  }
112  else{
113  std::cout<<"Only 0 or 1 is allowed"<<std::endl;
114  exit (1);
115  }
116  int Ntot_spec = Nuniv_genie + Nuniv_ppfx + Nshift_other;
117 
118  ///////////////////////////////////////////////////////////////
119 
120  const int wanted_pdg = 14;
121  const int number_of_genie_universes = 1000;
122  const Binning FluxBinning = Binning::Simple(480, 0.0, 120.0);
123  const int Ntotal_univ_ppfx = 100;
124 
126  const Var wgt_xs = VarFromNuTruthVar(wgt_xsST, 1);
128  const Var wgt = VarFromNuTruthVar(wgtST, 1);
129  const std::vector<NuTruthVar> vppfxST = ana::GetkPPFXFluxUnivWgtST();
130  const std::vector<Var> vppfx = ana::GetkPPFXFluxUnivWgt();
131  GenieMultiverseParameters verse(number_of_genie_universes, getAllXsecSysts_2018_RPAFix());
132  std::vector<ana::SystShifts> genie_shifts = verse.GetSystShifts();
133 
134  //////// Focusing weights ///////////////////////////////////////////////////////
135  //Long definitions per focusing systematic weights. We use it as a weight because we also use it to derive the flux.
136  const char* path_srt = getenv("SRT_PUBLIC_CONTEXT");
137  std::string beamfoc_file = std::string(path_srt) + "/CAFAna/data/beam/TABeamSyst.root";
138  TFile fFoc(beamfoc_file.c_str(),"read");
139  const int Nbsys_type = 11;
140  const std::string bsys_type[Nbsys_type] = {"2kA","0.2mmBeamSpotSize","1mmBeamShiftX","1mmBeamShiftY","3mmHorn1X","3mmHorn1Y","3mmHorn2X","3mmHorn2Y",
141  "7mmTargetZ","Magnetic Field in Decay Pipe","1mmHornWater"};
142  const std::string better_bsys_type[Nbsys_type] = {"2kA","02mmBeamSpotSize","1mmBeamShiftX","1mmBeamShiftY","3mmHorn1X","3mmHorn1Y","3mmHorn2X","3mmHorn2Y",
143  "7mmTargetZ","MagneticFieldinDecayPipe","1mmHornWater"};
144  const int Nhel = 4;
145  std::string hel[Nhel] = {"numu","anumu","nue","anue"};
146  const int Nud = 2;
147  std::string ud[Nud] = {"plussigma","minussigma"};
148  std::string ud_name[Nud] = {"bfoc_up","bfoc_dw"};
149 
150  TH1D* hfoc_up[Nbsys_type][Nhel];
151  TH1D* hfoc_dw[Nbsys_type][Nhel];
152  const int Nfoc_array = Nbsys_type*Nhel;
153 
154  for(int i=0;i<Nbsys_type;i++){
155  for(int j=0;j<Nhel;j++){
156  hfoc_up[i][j] = (TH1D*)fFoc.Get(Form("enuMF_FHC_%s_ND_%s_%s",bsys_type[i].c_str(),hel[j].c_str(),ud[0].c_str()));
157  hfoc_dw[i][j] = (TH1D*)fFoc.Get(Form("enuMF_FHC_%s_ND_%s_%s",bsys_type[i].c_str(),hel[j].c_str(),ud[1].c_str()));
158 
159  }
160  }
161 
162  const TH1D* hfoc_up_array[Nfoc_array] = {hfoc_up[0][0] ,hfoc_up[0][1] ,hfoc_up[0][2] ,hfoc_up[0][3],
163  hfoc_up[1][0] ,hfoc_up[1][1] ,hfoc_up[1][2] ,hfoc_up[1][3],
164  hfoc_up[2][0] ,hfoc_up[2][1] ,hfoc_up[2][2] ,hfoc_up[2][3],
165  hfoc_up[3][0] ,hfoc_up[3][1] ,hfoc_up[3][2] ,hfoc_up[3][3],
166  hfoc_up[4][0] ,hfoc_up[4][1] ,hfoc_up[4][2] ,hfoc_up[4][3],
167  hfoc_up[5][0] ,hfoc_up[5][1] ,hfoc_up[5][2] ,hfoc_up[5][3],
168  hfoc_up[6][0] ,hfoc_up[6][1] ,hfoc_up[6][2] ,hfoc_up[6][3],
169  hfoc_up[7][0] ,hfoc_up[7][1] ,hfoc_up[7][2] ,hfoc_up[7][3],
170  hfoc_up[8][0] ,hfoc_up[8][1] ,hfoc_up[8][2] ,hfoc_up[8][3],
171  hfoc_up[9][0] ,hfoc_up[9][1] ,hfoc_up[9][2] ,hfoc_up[9][3],
172  hfoc_up[10][0],hfoc_up[10][1],hfoc_up[10][2],hfoc_up[10][3]};
173  const TH1D* hfoc_dw_array[Nfoc_array] = {hfoc_dw[0][0] ,hfoc_dw[0][1] ,hfoc_dw[0][2] ,hfoc_dw[0][3],
174  hfoc_dw[1][0] ,hfoc_dw[1][1] ,hfoc_dw[1][2] ,hfoc_dw[1][3],
175  hfoc_dw[2][0] ,hfoc_dw[2][1] ,hfoc_dw[2][2] ,hfoc_dw[2][3],
176  hfoc_dw[3][0] ,hfoc_dw[3][1] ,hfoc_dw[3][2] ,hfoc_dw[3][3],
177  hfoc_dw[4][0] ,hfoc_dw[4][1] ,hfoc_dw[4][2] ,hfoc_dw[4][3],
178  hfoc_dw[5][0] ,hfoc_dw[5][1] ,hfoc_dw[5][2] ,hfoc_dw[5][3],
179  hfoc_dw[6][0] ,hfoc_dw[6][1] ,hfoc_dw[6][2] ,hfoc_dw[6][3],
180  hfoc_dw[7][0] ,hfoc_dw[7][1] ,hfoc_dw[7][2] ,hfoc_dw[7][3],
181  hfoc_dw[8][0] ,hfoc_dw[8][1] ,hfoc_dw[8][2] ,hfoc_dw[8][3],
182  hfoc_dw[9][0] ,hfoc_dw[9][1] ,hfoc_dw[9][2] ,hfoc_dw[9][3],
183  hfoc_dw[10][0],hfoc_dw[10][1],hfoc_dw[10][2],hfoc_dw[10][3]};
184 
185  const NuTruthVar BeamFoc_2KA_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
186  int ihel = -1;
187  if(nu->pdg == 14)ihel=0;
188  if(nu->pdg ==-14)ihel=1;
189  if(nu->pdg == 12)ihel=2;
190  if(nu->pdg ==-12)ihel=3;
191  if(ihel < 0)return 1.;
192  double nuener = nu->E;
193  int nubin = hfoc_up_array[4*0+ihel]->GetXaxis()->FindBin(nuener);
194  double wgt = hfoc_up_array[4*0+ihel]->GetBinContent(nubin);
195  return wgt;
196  });
197  const NuTruthVar BeamFoc_2KA_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
198  int ihel = -1;
199  if(nu->pdg == 14)ihel=0;
200  if(nu->pdg ==-14)ihel=1;
201  if(nu->pdg == 12)ihel=2;
202  if(nu->pdg ==-12)ihel=3;
203  if(ihel < 0)return 1.;
204  double nuener = nu->E;
205  int nubin = hfoc_dw_array[4*0+ihel]->GetXaxis()->FindBin(nuener);
206  double wgt = hfoc_dw_array[4*0+ihel]->GetBinContent(nubin);
207  return wgt;
208  });
209  const NuTruthVar BeamFoc_02mmBeamSpotSize_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
210  int ihel = -1;
211  if(nu->pdg == 14)ihel=0;
212  if(nu->pdg ==-14)ihel=1;
213  if(nu->pdg == 12)ihel=2;
214  if(nu->pdg ==-12)ihel=3;
215  if(ihel < 0)return 1.;
216  double nuener = nu->E;
217  int nubin = hfoc_up_array[4*1+ihel]->GetXaxis()->FindBin(nuener);
218  double wgt = hfoc_up_array[4*1+ihel]->GetBinContent(nubin);
219  return wgt;
220  });
221  const NuTruthVar BeamFoc_02mmBeamSpotSize_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
222  int ihel = -1;
223  if(nu->pdg == 14)ihel=0;
224  if(nu->pdg ==-14)ihel=1;
225  if(nu->pdg == 12)ihel=2;
226  if(nu->pdg ==-12)ihel=3;
227  if(ihel < 0)return 1.;
228  double nuener = nu->E;
229  int nubin = hfoc_dw_array[4*1+ihel]->GetXaxis()->FindBin(nuener);
230  double wgt = hfoc_dw_array[4*1+ihel]->GetBinContent(nubin);
231  return wgt;
232  });
233  const NuTruthVar BeamFoc_1mmBeamShiftX_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
234  int ihel = -1;
235  if(nu->pdg == 14)ihel=0;
236  if(nu->pdg ==-14)ihel=1;
237  if(nu->pdg == 12)ihel=2;
238  if(nu->pdg ==-12)ihel=3;
239  if(ihel < 0)return 1.;
240  double nuener = nu->E;
241  int nubin = hfoc_up_array[4*2+ihel]->GetXaxis()->FindBin(nuener);
242  double wgt = hfoc_up_array[4*2+ihel]->GetBinContent(nubin);
243  return wgt;
244  });
245  const NuTruthVar BeamFoc_1mmBeamShiftX_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
246  int ihel = -1;
247  if(nu->pdg == 14)ihel=0;
248  if(nu->pdg ==-14)ihel=1;
249  if(nu->pdg == 12)ihel=2;
250  if(nu->pdg ==-12)ihel=3;
251  if(ihel < 0)return 1.;
252  double nuener = nu->E;
253  int nubin = hfoc_dw_array[4*2+ihel]->GetXaxis()->FindBin(nuener);
254  double wgt = hfoc_dw_array[4*2+ihel]->GetBinContent(nubin);
255  return wgt;
256  });
257  const NuTruthVar BeamFoc_1mmBeamShiftY_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
258  int ihel = -1;
259  if(nu->pdg == 14)ihel=0;
260  if(nu->pdg ==-14)ihel=1;
261  if(nu->pdg == 12)ihel=2;
262  if(nu->pdg ==-12)ihel=3;
263  if(ihel < 0)return 1.;
264  double nuener = nu->E;
265  int nubin = hfoc_up_array[4*3+ihel]->GetXaxis()->FindBin(nuener);
266  double wgt = hfoc_up_array[4*3+ihel]->GetBinContent(nubin);
267  return wgt;
268  });
269  const NuTruthVar BeamFoc_1mmBeamShiftY_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
270  int ihel = -1;
271  if(nu->pdg == 14)ihel=0;
272  if(nu->pdg ==-14)ihel=1;
273  if(nu->pdg == 12)ihel=2;
274  if(nu->pdg ==-12)ihel=3;
275  if(ihel < 0)return 1.;
276  double nuener = nu->E;
277  int nubin = hfoc_dw_array[4*3+ihel]->GetXaxis()->FindBin(nuener);
278  double wgt = hfoc_dw_array[4*3+ihel]->GetBinContent(nubin);
279  return wgt;
280  });
281  const NuTruthVar BeamFoc_3mmHorn1X_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
282  int ihel = -1;
283  if(nu->pdg == 14)ihel=0;
284  if(nu->pdg ==-14)ihel=1;
285  if(nu->pdg == 12)ihel=2;
286  if(nu->pdg ==-12)ihel=3;
287  if(ihel < 0)return 1.;
288  double nuener = nu->E;
289  int nubin = hfoc_up_array[4*4+ihel]->GetXaxis()->FindBin(nuener);
290  double wgt = hfoc_up_array[4*4+ihel]->GetBinContent(nubin);
291  return wgt;
292  });
293  const NuTruthVar BeamFoc_3mmHorn1X_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
294  int ihel = -1;
295  if(nu->pdg == 14)ihel=0;
296  if(nu->pdg ==-14)ihel=1;
297  if(nu->pdg == 12)ihel=2;
298  if(nu->pdg ==-12)ihel=3;
299  if(ihel < 0)return 1.;
300  double nuener = nu->E;
301  int nubin = hfoc_dw_array[4*4+ihel]->GetXaxis()->FindBin(nuener);
302  double wgt = hfoc_dw_array[4*4+ihel]->GetBinContent(nubin);
303  return wgt;
304  });
305  const NuTruthVar BeamFoc_3mmHorn1Y_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
306  int ihel = -1;
307  if(nu->pdg == 14)ihel=0;
308  if(nu->pdg ==-14)ihel=1;
309  if(nu->pdg == 12)ihel=2;
310  if(nu->pdg ==-12)ihel=3;
311  if(ihel < 0)return 1.;
312  double nuener = nu->E;
313  int nubin = hfoc_up_array[4*5+ihel]->GetXaxis()->FindBin(nuener);
314  double wgt = hfoc_up_array[4*5+ihel]->GetBinContent(nubin);
315  return wgt;
316  });
317  const NuTruthVar BeamFoc_3mmHorn1Y_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
318  int ihel = -1;
319  if(nu->pdg == 14)ihel=0;
320  if(nu->pdg ==-14)ihel=1;
321  if(nu->pdg == 12)ihel=2;
322  if(nu->pdg ==-12)ihel=3;
323  if(ihel < 0)return 1.;
324  double nuener = nu->E;
325  int nubin = hfoc_dw_array[4*5+ihel]->GetXaxis()->FindBin(nuener);
326  double wgt = hfoc_dw_array[4*5+ihel]->GetBinContent(nubin);
327  return wgt;
328  });
329  const NuTruthVar BeamFoc_3mmHorn2X_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
330  int ihel = -1;
331  if(nu->pdg == 14)ihel=0;
332  if(nu->pdg ==-14)ihel=1;
333  if(nu->pdg == 12)ihel=2;
334  if(nu->pdg ==-12)ihel=3;
335  if(ihel < 0)return 1.;
336  double nuener = nu->E;
337  int nubin = hfoc_up_array[4*6+ihel]->GetXaxis()->FindBin(nuener);
338  double wgt = hfoc_up_array[4*6+ihel]->GetBinContent(nubin);
339  return wgt;
340  });
341  const NuTruthVar BeamFoc_3mmHorn2X_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
342  int ihel = -1;
343  if(nu->pdg == 14)ihel=0;
344  if(nu->pdg ==-14)ihel=1;
345  if(nu->pdg == 12)ihel=2;
346  if(nu->pdg ==-12)ihel=3;
347  if(ihel < 0)return 1.;
348  double nuener = nu->E;
349  int nubin = hfoc_dw_array[4*6+ihel]->GetXaxis()->FindBin(nuener);
350  double wgt = hfoc_dw_array[4*6+ihel]->GetBinContent(nubin);
351  return wgt;
352  });
353  const NuTruthVar BeamFoc_3mmHorn2Y_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
354  int ihel = -1;
355  if(nu->pdg == 14)ihel=0;
356  if(nu->pdg ==-14)ihel=1;
357  if(nu->pdg == 12)ihel=2;
358  if(nu->pdg ==-12)ihel=3;
359  if(ihel < 0)return 1.;
360  double nuener = nu->E;
361  int nubin = hfoc_up_array[4*7+ihel]->GetXaxis()->FindBin(nuener);
362  double wgt = hfoc_up_array[4*7+ihel]->GetBinContent(nubin);
363  return wgt;
364  });
365  const NuTruthVar BeamFoc_3mmHorn2Y_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
366  int ihel = -1;
367  if(nu->pdg == 14)ihel=0;
368  if(nu->pdg ==-14)ihel=1;
369  if(nu->pdg == 12)ihel=2;
370  if(nu->pdg ==-12)ihel=3;
371  if(ihel < 0)return 1.;
372  double nuener = nu->E;
373  int nubin = hfoc_dw_array[4*7+ihel]->GetXaxis()->FindBin(nuener);
374  double wgt = hfoc_dw_array[4*7+ihel]->GetBinContent(nubin);
375  return wgt;
376  });
377  const NuTruthVar BeamFoc_7mmTargetZ_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
378  int ihel = -1;
379  if(nu->pdg == 14)ihel=0;
380  if(nu->pdg ==-14)ihel=1;
381  if(nu->pdg == 12)ihel=2;
382  if(nu->pdg ==-12)ihel=3;
383  if(ihel < 0)return 1.;
384  double nuener = nu->E;
385  int nubin = hfoc_up_array[4*8+ihel]->GetXaxis()->FindBin(nuener);
386  double wgt = hfoc_up_array[4*8+ihel]->GetBinContent(nubin);
387  return wgt;
388  });
389  const NuTruthVar BeamFoc_7mmTargetZ_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
390  int ihel = -1;
391  if(nu->pdg == 14)ihel=0;
392  if(nu->pdg ==-14)ihel=1;
393  if(nu->pdg == 12)ihel=2;
394  if(nu->pdg ==-12)ihel=3;
395  if(ihel < 0)return 1.;
396  double nuener = nu->E;
397  int nubin = hfoc_dw_array[4*8+ihel]->GetXaxis()->FindBin(nuener);
398  double wgt = hfoc_dw_array[4*8+ihel]->GetBinContent(nubin);
399  return wgt;
400  });
401  const NuTruthVar BeamFoc_MagneticFieldinDecayPipe_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
402  int ihel = -1;
403  if(nu->pdg == 14)ihel=0;
404  if(nu->pdg ==-14)ihel=1;
405  if(nu->pdg == 12)ihel=2;
406  if(nu->pdg ==-12)ihel=3;
407  if(ihel < 0)return 1.;
408  double nuener = nu->E;
409  int nubin = hfoc_up_array[4*9+ihel]->GetXaxis()->FindBin(nuener);
410  double wgt = hfoc_up_array[4*9+ihel]->GetBinContent(nubin);
411  return wgt;
412  });
413  const NuTruthVar BeamFoc_MagneticFieldinDecayPipe_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
414  int ihel = -1;
415  if(nu->pdg == 14)ihel=0;
416  if(nu->pdg ==-14)ihel=1;
417  if(nu->pdg == 12)ihel=2;
418  if(nu->pdg ==-12)ihel=3;
419  if(ihel < 0)return 1.;
420  double nuener = nu->E;
421  int nubin = hfoc_dw_array[4*9+ihel]->GetXaxis()->FindBin(nuener);
422  double wgt = hfoc_dw_array[4*9+ihel]->GetBinContent(nubin);
423  return wgt;
424  });
425  const NuTruthVar BeamFoc_1mmHornWater_WgtUpST([hfoc_up_array](const caf::SRNeutrinoProxy* nu){
426  int ihel = -1;
427  if(nu->pdg == 14)ihel=0;
428  if(nu->pdg ==-14)ihel=1;
429  if(nu->pdg == 12)ihel=2;
430  if(nu->pdg ==-12)ihel=3;
431  if(ihel < 0)return 1.;
432  double nuener = nu->E;
433  int nubin = hfoc_up_array[4*10+ihel]->GetXaxis()->FindBin(nuener);
434  double wgt = hfoc_up_array[4*10+ihel]->GetBinContent(nubin);
435  return wgt;
436  });
437  const NuTruthVar BeamFoc_1mmHornWater_WgtDwST([hfoc_dw_array](const caf::SRNeutrinoProxy* nu){
438  int ihel = -1;
439  if(nu->pdg == 14)ihel=0;
440  if(nu->pdg ==-14)ihel=1;
441  if(nu->pdg == 12)ihel=2;
442  if(nu->pdg ==-12)ihel=3;
443  if(ihel < 0)return 1.;
444  double nuener = nu->E;
445  int nubin = hfoc_dw_array[4*10+ihel]->GetXaxis()->FindBin(nuener);
446  double wgt = hfoc_dw_array[4*10+ihel]->GetBinContent(nubin);
447  return wgt;
448  });
449 
450  const Var BeamFoc_2KA_WgtUp = VarFromNuTruthVar(BeamFoc_2KA_WgtUpST, 1);
451  const Var BeamFoc_2KA_WgtDw = VarFromNuTruthVar(BeamFoc_2KA_WgtDwST, 1);
452  const Var BeamFoc_02mmBeamSpotSize_WgtUp = VarFromNuTruthVar(BeamFoc_02mmBeamSpotSize_WgtUpST, 1);
453  const Var BeamFoc_02mmBeamSpotSize_WgtDw = VarFromNuTruthVar(BeamFoc_02mmBeamSpotSize_WgtDwST, 1);
454  const Var BeamFoc_1mmBeamShiftX_WgtUp = VarFromNuTruthVar(BeamFoc_1mmBeamShiftX_WgtUpST, 1);
455  const Var BeamFoc_1mmBeamShiftX_WgtDw = VarFromNuTruthVar(BeamFoc_1mmBeamShiftX_WgtDwST, 1);
456  const Var BeamFoc_1mmBeamShiftY_WgtUp = VarFromNuTruthVar(BeamFoc_1mmBeamShiftY_WgtUpST, 1);
457  const Var BeamFoc_1mmBeamShiftY_WgtDw = VarFromNuTruthVar(BeamFoc_1mmBeamShiftY_WgtDwST, 1);
458  const Var BeamFoc_3mmHorn1X_WgtUp = VarFromNuTruthVar(BeamFoc_3mmHorn1X_WgtUpST, 1);
459  const Var BeamFoc_3mmHorn1X_WgtDw = VarFromNuTruthVar(BeamFoc_3mmHorn1X_WgtDwST, 1);
460  const Var BeamFoc_3mmHorn1Y_WgtUp = VarFromNuTruthVar(BeamFoc_3mmHorn1Y_WgtUpST, 1);
461  const Var BeamFoc_3mmHorn1Y_WgtDw = VarFromNuTruthVar(BeamFoc_3mmHorn1Y_WgtDwST, 1);
462  const Var BeamFoc_3mmHorn2X_WgtUp = VarFromNuTruthVar(BeamFoc_3mmHorn2X_WgtUpST, 1);
463  const Var BeamFoc_3mmHorn2X_WgtDw = VarFromNuTruthVar(BeamFoc_3mmHorn2X_WgtDwST, 1);
464  const Var BeamFoc_3mmHorn2Y_WgtUp = VarFromNuTruthVar(BeamFoc_3mmHorn2Y_WgtUpST, 1);
465  const Var BeamFoc_3mmHorn2Y_WgtDw = VarFromNuTruthVar(BeamFoc_3mmHorn2Y_WgtDwST, 1);
466  const Var BeamFoc_7mmTargetZ_WgtUp = VarFromNuTruthVar(BeamFoc_7mmTargetZ_WgtUpST, 1);
467  const Var BeamFoc_7mmTargetZ_WgtDw = VarFromNuTruthVar(BeamFoc_7mmTargetZ_WgtDwST, 1);
468  const Var BeamFoc_MagneticFieldinDecayPipe_WgtUp = VarFromNuTruthVar(BeamFoc_MagneticFieldinDecayPipe_WgtUpST, 1);
469  const Var BeamFoc_MagneticFieldinDecayPipe_WgtDw = VarFromNuTruthVar(BeamFoc_MagneticFieldinDecayPipe_WgtDwST, 1);
470  const Var BeamFoc_1mmHornWater_WgtUp = VarFromNuTruthVar(BeamFoc_1mmHornWater_WgtUpST, 1);
471  const Var BeamFoc_1mmHornWater_WgtDw = VarFromNuTruthVar(BeamFoc_1mmHornWater_WgtDwST, 1);
472 
473  const NuTruthVar BeamFoc_WgtST[Nbsys_type][Nud] = { {BeamFoc_2KA_WgtUpST,BeamFoc_2KA_WgtDwST},
474  {BeamFoc_02mmBeamSpotSize_WgtUpST, BeamFoc_02mmBeamSpotSize_WgtDwST},
475  {BeamFoc_1mmBeamShiftX_WgtUpST,BeamFoc_1mmBeamShiftX_WgtDwST},
476  {BeamFoc_1mmBeamShiftY_WgtUpST,BeamFoc_1mmBeamShiftY_WgtDwST},
477  {BeamFoc_3mmHorn1X_WgtUpST,BeamFoc_3mmHorn1X_WgtDwST},
478  {BeamFoc_3mmHorn1Y_WgtUpST,BeamFoc_3mmHorn1Y_WgtDwST},
479  {BeamFoc_3mmHorn2X_WgtUpST,BeamFoc_3mmHorn2X_WgtDwST},
480  {BeamFoc_3mmHorn2Y_WgtUpST,BeamFoc_3mmHorn2Y_WgtDwST},
481  {BeamFoc_7mmTargetZ_WgtUpST,BeamFoc_7mmTargetZ_WgtDwST},
482  {BeamFoc_MagneticFieldinDecayPipe_WgtUpST,BeamFoc_MagneticFieldinDecayPipe_WgtDwST},
483  {BeamFoc_1mmHornWater_WgtUpST,BeamFoc_1mmHornWater_WgtDwST}};
484  const Var BeamFoc_Wgt[Nbsys_type][Nud] = { {BeamFoc_2KA_WgtUp,BeamFoc_2KA_WgtDw},
485  {BeamFoc_02mmBeamSpotSize_WgtUp, BeamFoc_02mmBeamSpotSize_WgtDw},
486  {BeamFoc_1mmBeamShiftX_WgtUp,BeamFoc_1mmBeamShiftX_WgtDw},
487  {BeamFoc_1mmBeamShiftY_WgtUp,BeamFoc_1mmBeamShiftY_WgtDw},
488  {BeamFoc_3mmHorn1X_WgtUp,BeamFoc_3mmHorn1X_WgtDw},
489  {BeamFoc_3mmHorn1Y_WgtUp,BeamFoc_3mmHorn1Y_WgtDw},
490  {BeamFoc_3mmHorn2X_WgtUp,BeamFoc_3mmHorn2X_WgtDw},
491  {BeamFoc_3mmHorn2Y_WgtUp,BeamFoc_3mmHorn2Y_WgtDw},
492  {BeamFoc_7mmTargetZ_WgtUp,BeamFoc_7mmTargetZ_WgtDw},
493  {BeamFoc_MagneticFieldinDecayPipe_WgtUp,BeamFoc_MagneticFieldinDecayPipe_WgtDw},
494  {BeamFoc_1mmHornWater_WgtUp,BeamFoc_1mmHornWater_WgtDw}};
495 
496  //////Declarations for the muon energy scale systematics:
497  const int Nmues = 4;
498  const string mues[Nmues] = {"MuESUp","MuESDw","MuES_onlyMuKEUp","MuES_onlyMuKEDw"};
499  const Cut cut_mues[Nmues] = {kAllNumuCCCutsUp,kAllNumuCCCutsDw,kAllNumuCCCutsUp,kAllNumuCCCutsDw};
502 
503  //////// Datasets ///////////////////////////////////////////////////////
504  const int Ndatasets = 8;
505  const std::string datasets[Ndatasets] = {"defname: nd_fhc_ndrespin_caf_1_of_5 or defname: nd_fhc_ndrespin_caf_2_of_5 or defname: nd_fhc_ndrespin_caf_3_of_5 or defname: nd_fhc_ndrespin_caf_4_of_5",
506  "defname: nd_fhc_ndrespin_caf_0_of_5",
507  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1",
508  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1",
509  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1",
510  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1",
511  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1",
512  "dataset_def_name_newest_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1"};
513 
514  const std::string dataset_name[Ndatasets] = {"MC","mockdata","lightdown","lightup","ckv","calibneg","calibpos","calibshape"};
515 
516  std::vector<Cut> cuts = {!kIsTrueSig};
517 
518  //////// Declarations ///////////////////////////////////////////////////////
519 
520  //Flux
521  std::vector<Spectrum*> vFlux;
522  //Efficiency
523  std::vector< std::pair<Spectrum*,Spectrum*> > veff;
524  std::vector<Spectrum*> vMCSig;
525  std::vector<Spectrum*> vMCSigNuTree;
526  //Background
527  std::vector<TrivialBkgdEstimator*> vbkg;
528  std::vector<IBkgdEstimator*> vBkgdEst;
529  //Fake-data
530  std::vector<Spectrum*> vDt;
531  //Migration matrix
532  std::vector<ReweightableSpectrum*> vRecoTrue;
533 
534  veff.reserve(Ntot_spec);
535  vMCSig.reserve(Ntot_spec);
536  vMCSigNuTree.reserve(Ntot_spec);
537  vbkg.reserve(Ntot_spec);
538  vBkgdEst.reserve(Ntot_spec);
539  vRecoTrue.reserve(Ntot_spec);
540  if(do_allother){
541  vFlux.reserve(123); // 1 CV, 100 ppfx and 22 beam shift
542  vDt.reserve(2); // Mock-data: no weight and PPFX weight
543  }
544  std::vector<SpectrumLoaderBase*> vloader;
545 
546 
547  // bool NERSC = (start_file_num >= 0) && (end_file_num >= start_file_num);
548  bool NERSC = true;
549  for(int i=0;i<Ndatasets;i++){
550  if(!do_allother && i>0) continue;
551  if(!NERSC)
552  vloader.push_back( new SpectrumLoader(datasets[i]) );
553  else
554  {
555  vector<string> filelist = FileListFromDataset(i);
556  int nfiles = filelist.size();
557  if((start_file_num < nfiles) && (end_file_num < nfiles))
558  {
559  vector<string>::const_iterator idx1 = filelist.begin() + start_file_num;
560  vector<string>::const_iterator idx2 = filelist.begin() + end_file_num + 1;
561  vector<string> chosen_files(idx1, idx2);
562  vloader.push_back(new SpectrumLoader(chosen_files));
563  }
564  else
565  vloader.push_back(&kNullLoader);
566  }
567 
568  vloader[i]->SetSpillCut(kStandardSpillCuts);
569  }
570 
571  //////// Spectra ///////////////////////////////////////////////////////
572  //GENIE
573  for(int i=ifirst_genie_univ;i<=ilast_genie_univ && ifirst_genie_univ >= 0;i++){
574  veff.push_back( ana::Efficiency(*vloader[0],kTrueMuKEVsCosVsEnuStandardAxisST,kAllNumuCCCuts, kIsTrueSigST, genie_shifts[i], wgtST) );
575  vbkg.push_back( new TrivialBkgdEstimator(*vloader[0], kRecoMuKEVsCosVsEnuStandardAxis, kAllNumuCCCuts, cuts, genie_shifts[i], wgt) );
578  genie_shifts[i], wgt) );
579  }
580  //PPFX
581  for(int i=ifirst_ppfx_univ;i<=ilast_ppfx_univ && ifirst_ppfx_univ >= 0;i++){
582  veff.push_back( ana::Efficiency(*vloader[0],kTrueMuKEVsCosVsEnuStandardAxisST,kAllNumuCCCuts, kIsTrueSigST, kNoShift, wgt_xsST * vppfxST[i]) );
583  vbkg.push_back( new TrivialBkgdEstimator(*vloader[0], kRecoMuKEVsCosVsEnuStandardAxis, kAllNumuCCCuts, cuts, kNoShift, wgt_xs * vppfx[i]) );
586  kNoShift, wgt_xs * vppfx[i]) );
587  }
588  if(do_allother){
589  //CV:
591  vbkg.push_back( new TrivialBkgdEstimator(*vloader[0], kRecoMuKEVsCosVsEnuStandardAxis, kAllNumuCCCuts, cuts, kNoShift,wgt) );
594  kNoShift, wgt) );
595  //Det syst: datasets: 2-7
596  for(int i=2;i<Ndatasets;i++){
598  vbkg.push_back( new TrivialBkgdEstimator(*vloader[i], kRecoMuKEVsCosVsEnuStandardAxis, kAllNumuCCCuts, cuts, kNoShift,wgt) );
601  kNoShift, wgt) );
602  }
603  //Focusing systematics
604  // - Up:
605  for(int i=0;i<Nbsys_type;i++){
606  veff.push_back( ana::Efficiency(*vloader[0], kTrueMuKEVsCosVsEnuStandardAxisST,kAllNumuCCCuts, kIsTrueSigST, kNoShift, wgtST*BeamFoc_WgtST[i][0]) );
607  vbkg.push_back( new TrivialBkgdEstimator(*vloader[0], kRecoMuKEVsCosVsEnuStandardAxis, kAllNumuCCCuts, cuts, kNoShift,wgt*BeamFoc_Wgt[i][0]) );
610  kNoShift, wgt*BeamFoc_Wgt[i][0]) );
611  }
612  // - Down:
613  for(int i=0;i<Nbsys_type;i++){
614  veff.push_back( ana::Efficiency(*vloader[0], kTrueMuKEVsCosVsEnuStandardAxisST,kAllNumuCCCuts, kIsTrueSigST, kNoShift, wgtST*BeamFoc_WgtST[i][1]) );
615  vbkg.push_back( new TrivialBkgdEstimator(*vloader[0], kRecoMuKEVsCosVsEnuStandardAxis, kAllNumuCCCuts, cuts, kNoShift,wgt*BeamFoc_Wgt[i][1]) );
618  kNoShift, wgt*BeamFoc_Wgt[i][1]) );
619  }
620  //Muon energy scale systematics
621  for(int i=0;i<Nmues;i++){
622  veff.push_back( ana::Efficiency(*vloader[0],kTrueMuKEVsCosVsEnuStandardAxisST,cut_mues[i], kIsTrueSigST, kNoShift, wgtST) );
623  vbkg.push_back( new TrivialBkgdEstimator(*vloader[0], axis_mues[i], cut_mues[i], cuts, kNoShift,wgt) );
624  vRecoTrue.push_back( new ReweightableSpectrum(*vloader[0], axis_mues[i], kTrueMuKEVsCosVsEnuStandardAxis,
625  cut_mues[i] && CutFromNuTruthCut(kIsTrueSigST),
626  kNoShift, wgt) );
627  }
628  }
629  //To the final spectra:
630  for(unsigned i=0;i<veff.size();i++){
631  vMCSig.push_back( veff[i].first );
632  vMCSigNuTree.push_back( veff[i].second );
633  vBkgdEst.push_back( vbkg[i]);
634  }
635 
636  if(do_allother){
637  ////////////// Fluxes /////////////////////////////////////////////////
638  vFlux.push_back( DeriveFlux(*vloader[0], FluxBinning, wanted_pdg, &vtxmin, &vtxmax, kNoShift, kPPFXFluxCVWgtST) );
639  for(int i=0;i<Ntotal_univ_ppfx;i++){
640  vFlux.push_back( DeriveFlux(*vloader[0], FluxBinning, wanted_pdg, &vtxmin, &vtxmax, kNoShift, vppfxST[i]) );
641  }
642  for(int i=0;i<Nbsys_type;i++){
643  vFlux.push_back( DeriveFlux(*vloader[0], FluxBinning, wanted_pdg, &vtxmin, &vtxmax, kNoShift, kPPFXFluxCVWgtST * BeamFoc_WgtST[i][0]) );
644  vFlux.push_back( DeriveFlux(*vloader[0], FluxBinning, wanted_pdg, &vtxmin, &vtxmax, kNoShift, kPPFXFluxCVWgtST * BeamFoc_WgtST[i][1]) );
645  }
646 
647  ////////////// Fake data /////////////////////////////////////////////////
648  vDt.push_back( new Spectrum(*vloader[1], kRecoMuKEVsCosVsEnuStandardAxis, kAllNumuCCCuts, kNoShift, kUnweighted) );
649  vDt.push_back( new Spectrum(*vloader[1], kRecoMuKEVsCosVsEnuStandardAxis, kAllNumuCCCuts, kNoShift, kPPFXFluxCVWgt) );
650  }
651 
652  ////////////// Let's Go/////////////////////////////////////////////////
653  for(int i=0;i<Ndatasets;i++){
654  if(!do_allother && i>0)continue;
655  vloader[i]->Go();
656  }
657 
658  /////////////// Storing ////////////////////////////////////////////////
659 
660  if(fOutname.compare("genie") == 0)
661  fOutname = Form("/output/%s_univ_%d_to_%d_files_%d_to_%d.root",
662  fOutname.c_str(), ifirst_genie_univ, ilast_genie_univ,
663  start_file_num, end_file_num);
664  else if(fOutname.compare("ppfx") == 0)
665  fOutname = Form("/output/%s_univ_%d_to_%d_files_%d_to_%d.root",
666  fOutname.c_str(), ifirst_ppfx_univ, ilast_ppfx_univ,
667  start_file_num, end_file_num);
668  else if (fOutname.compare("systs") == 0)
669  fOutname = Form("/output/%s_files_%d_to_%d.root",
670  fOutname.c_str(), start_file_num, end_file_num);
671  else {
672  cerr << endl << "ERROR! Valid modes are \"genie\", \"ppfx\" or \"systs\". Exiting..." << endl;
673  exit(-1);
674  }
675 
676  TFile* fOut = new TFile(fOutname.c_str(), "RECREATE");
677  fOut->cd();
678  TDirectory* dir;
679  //Spectra:
680  int icurr = -1;
681  // - GENIE
682  for(int i=ifirst_genie_univ;i<=ilast_genie_univ && ifirst_genie_univ >= 0;i++){
683  icurr++;
684  dir = fOut->mkdir(Form("genie%d",i));
685  dir->cd();
686  vMCSig[icurr]->SaveTo(dir->mkdir("MCSig"));
687  vMCSigNuTree[icurr]->SaveTo(dir->mkdir("MCSigNuTre"));
688  vBkgdEst[icurr]->SaveTo(dir->mkdir("BkgdEst"));
689  vRecoTrue[icurr]->SaveTo(dir->mkdir("RecoTrue"));
690 
691  }
692  // - PPFX
693  for(int i=ifirst_ppfx_univ;i<=ilast_ppfx_univ && ifirst_ppfx_univ >= 0;i++){
694  icurr++;
695  dir = fOut->mkdir(Form("ppfx%d",i));
696  dir->cd();
697  vMCSig[icurr]->SaveTo(dir->mkdir("MCSig"));
698  vMCSigNuTree[icurr]->SaveTo(dir->mkdir("MCSigNuTre"));
699  vBkgdEst[icurr]->SaveTo(dir->mkdir("BkgdEst"));
700  vRecoTrue[icurr]->SaveTo(dir->mkdir("RecoTrue"));
701 
702  }
703  if(do_allother){
704  // - CV
705  {
706  icurr++;
707  dir = fOut->mkdir("cv");
708  dir->cd();
709  vMCSig[icurr]->SaveTo(dir->mkdir("MCSig"));
710  vMCSigNuTree[icurr]->SaveTo(dir->mkdir("MCSigNuTre"));
711  vBkgdEst[icurr]->SaveTo(dir->mkdir("BkgdEst"));
712  vRecoTrue[icurr]->SaveTo(dir->mkdir("RecoTrue"));
713  }
714  // - Det Syst
715  for(int i=2;i<Ndatasets;i++){
716  icurr++;
717  dir = fOut->mkdir(Form("%s",dataset_name[i].c_str()));
718  dir->cd();
719  vMCSig[icurr]->SaveTo(dir->mkdir("MCSig"));
720  vMCSigNuTree[icurr]->SaveTo(dir->mkdir("MCSigNuTre"));
721  vBkgdEst[icurr]->SaveTo(dir->mkdir("BkgdEst"));
722  vRecoTrue[icurr]->SaveTo(dir->mkdir("RecoTrue"));
723  }
724  // - Beam focusing Up
725  for(int i=0;i<Nbsys_type;i++){
726  icurr++;
727  dir = fOut->mkdir(Form("%s_Up",better_bsys_type[i].c_str()));
728  dir->cd();
729  vMCSig[icurr]->SaveTo(dir->mkdir("MCSig"));
730  vMCSigNuTree[icurr]->SaveTo(dir->mkdir("MCSigNuTre"));
731  vBkgdEst[icurr]->SaveTo(dir->mkdir("BkgdEst"));
732  vRecoTrue[icurr]->SaveTo(dir->mkdir("RecoTrue"));
733  }
734  // - Beam focusing Dw
735  for(int i=0;i<Nbsys_type;i++){
736  icurr++;
737  dir = fOut->mkdir(Form("%s_Dw",better_bsys_type[i].c_str()));
738  dir->cd();
739  vMCSig[icurr]->SaveTo(dir->mkdir("MCSig"));
740  vMCSigNuTree[icurr]->SaveTo(dir->mkdir("MCSigNuTre"));
741  vBkgdEst[icurr]->SaveTo(dir->mkdir("BkgdEst"));
742  vRecoTrue[icurr]->SaveTo(dir->mkdir("RecoTrue"));
743  }
744  // - Muon energy scale:
745  for(int i=0;i<Nmues;i++){
746  icurr++;
747  dir = fOut->mkdir(Form("%s",mues[i].c_str()));
748  dir->cd();
749  vMCSig[icurr]->SaveTo(dir->mkdir("MCSig"));
750  vMCSigNuTree[icurr]->SaveTo(dir->mkdir("MCSigNuTre"));
751  vBkgdEst[icurr]->SaveTo(dir->mkdir("BkgdEst"));
752  vRecoTrue[icurr]->SaveTo(dir->mkdir("RecoTrue"));
753  }
754  // - Flux:
755  {
756  dir = fOut->mkdir("fluxes");
757  dir->cd();
758  unsigned int iflux = 0;
759  vFlux[iflux]->SaveTo(dir->mkdir("CV"));
760  iflux++;
761  for(int i=0;i<Ntotal_univ_ppfx;i++){
762  vFlux[iflux]->SaveTo(dir->mkdir(Form("ppfx%d",i)));
763  iflux++;
764  }
765  for(int i=0;i<Nbsys_type;i++){
766  vFlux[iflux]->SaveTo(dir->mkdir(Form("%s_Up",better_bsys_type[i].c_str())));
767  vFlux[iflux+1]->SaveTo(dir->mkdir(Form("%s_Dw",better_bsys_type[i].c_str())));
768  iflux += 2;
769  }
770  }
771  //Fake data:
772  {
773  dir = fOut->mkdir("fakedata");
774  dir->cd();
775  vDt[0]->SaveTo(dir->mkdir("unwgt"));
776  vDt[1]->SaveTo(dir->mkdir("ppfxwgt"));
777  }
778  }
779 
780  std::cout<<"=> numucc_inc_xs() end "<<std::endl;
781 
782 }
const NuTruthHistAxis kTrueMuKEVsCosVsEnuStandardAxisST("True T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kTrueMuKEVsCosVsEnuST)
const HistAxis kRecoMuKEVsCosVsEnuStandardAxis_onlyMuKEDw("Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnu_onlyMuKEDw)
const int Nuniv_ppfx
const NuTruthCut kIsTrueSigST
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TFile * inf
Definition: drawXsec.C:1
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const TVector3 vtxmin(-130,-176, 225)
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
const ana::Var wgt_xs
Spectrum with the value of a second variable, allowing for reweighting
OStream cerr
Definition: OStream.cxx:7
caf::Proxy< float > E
Definition: SRProxy.h:520
const HistAxis kTrueMuKEVsCosVsEnuStandardAxis
const int Nuniv_genie
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
std::string getenv(std::string const &name)
const HistAxis kRecoMuKEVsCosVsEnuStandardAxisUp("Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnuUp)
const HistAxis kRecoMuKEVsCosVsEnuStandardAxisDw("Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnuDw)
const double j
Definition: BetheBloch.cxx:29
std::vector< NuTruthVar > GetkPPFXFluxUnivWgtST()
Definition: PPFXWeights.cxx:33
const NuTruthVar kXSecCVWgt2018_smallerDISScale_NT
Definition: XsecTunes.h:48
std::vector< float > Spectrum
Definition: Constants.h:610
const ana::Var wgt
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
void specprod_numuccinc(std::string outfileName, int firstGenie, int lastGenie, int firstPPFX=-1, int lastPPFX=-1)
std::pair< Spectrum *, Spectrum * > Efficiency(SpectrumLoaderBase &loader, const NuTruthHistAxis &axis, const Cut &seln, const NuTruthCut &truthSeln, const SystShifts &shift, const NuTruthVar &weight)
Definition: Efficiency.cxx:10
const HistAxis kRecoMuKEVsCosVsEnuStandardAxis_onlyMuKEUp("Reco. T_{#mu} vs cos #{theta} vs Neutrino Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnu_onlyMuKEUp)
const Cut kIsTrueSig
std::vector< const ISyst * > getAllXsecSysts_2018_RPAFix()
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:7
Spectrum * DeriveFlux(SpectrumLoaderBase &loader, const Binning &bins, int pdg, const TVector3 *min, const TVector3 *max, const SystShifts &shift, const NuTruthVar weight)
Definition: Flux.cxx:58
const HistAxis kRecoMuKEVsCosVsEnuStandardAxis("Reco. T_{#mu} vs cos #{theta} vs Available Energy (GeV)", angvsmukevsebins, kRecoMuKEVsCosVsEnu)
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
TDirectory * dir
Definition: macro.C:5
vector< string > dataset_lists
const ana::NuTruthVar wgtST
const Cut kAllNumuCCCutsUp
std::vector< SystShifts > GetSystShifts()
vector< string > FileListFromDataset(int ds_idx)
exit(0)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
const TVector3 vtxmax(160, 160, 1000)
Template for Var and SpillVar.
const Cut kAllNumuCCCutsDw
std::vector< Dataset > datasets
Definition: Periods.h:3
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
const NuTruthVar kPPFXFluxCVWgtST([](const caf::SRNeutrinoProxy *nu){ if(nu->rwgt.ppfx.cv!=nu->rwgt.ppfx.cv){return 1.f;}if(nu->rwgt.ppfx.cv >90){return 1.f;}return float(nu->rwgt.ppfx.cv);})
weight events with the flux PPFX Central value correction.
Definition: PPFXWeights.h:12
Just return the MC prediction for the background.
const Cut kAllNumuCCCuts
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
const ana::NuTruthVar wgt_xsST
enum BeamMode string