Make_NuMuCC_Inc_XS_v2.C
Go to the documentation of this file.
1 
2 #include "CAFAna/Vars/Vars.h"
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Cuts/Cuts.h"
5 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Core/Ratio.h"
8 
15 
17 
18 #include "TFile.h"
19 #include "TH1.h"
20 #include "TH2.h"
21 #include "TH3.h"
22 
23 #include <iostream>
24 #include <cmath>
25 #include <string>
26 #include <iomanip>
27 #include <map>
28 
29 #include "Utilities/func/MathUtil.h"
30 
31 // directory utilities
32 #include <boost/algorithm/string/classification.hpp> // Include boost::is_any_of
33 #include <boost/algorithm/string/split.hpp> // Include for boost::split
34 #include <boost/filesystem.hpp>
35 #include <boost/filesystem/convenience.hpp>
36 
37 using namespace ana;
38 using namespace boost::filesystem;
39 
40 //const double ntargs = 3.88998e+31; //pre-calculated with updated fiducial volume(1e6 interations)
41 const double ntargs = 3.89137e+31; //pre-calculated with updated fiducial volume(1e10 interations)
42 const double pot = 8.09e20; //data POT
43 
44 const std::string fake_xs = "ppfxwgt";
45 const std::string xs_mode0 = "xs_nominal";
46 const std::string xs_mode1 = "xs_bkgcons";
47 
48 std::string dirIn = "input_files";
49 
50 const int Nuniv_ppfx = 0;
51 const int Nuniv_genie = 0;
52 
53 std::map<std::string, std::string> coreFiles{
54  {"nominal", "syst_specs_cv.root"}, // cv dataset
55  {"fluxes", "fake_data_fluxes.root"}, // Flux dataset
56  {"fakedata", "fake_data_fluxes.root"}, // Fake Data
57  {"genie", ""},
58  {"ppfx", ""}
59 };
60 
61 const std::vector<std::string> unfoldVars3D{
62  "EAvail"
63 };
64 
65 const std::vector<std::string> unfoldVars1D{
66  "ENu",
67  "Q2"
68 };
69 
70 std::map<std::string, std::string> var2FD{
71  {"EAvail", "fakedata"},
72  {"ENu", "fakedata_enu"},
73  {"Q2", "fakedata_q2"}
74 };
75 
76 const std::map<std::string, std::string> systFiles = {
77  // {"lightup", "syst_specs_lightup.root"},
78  // {"lightdown", "syst_specs_lightdw.root"},
79  // {"calibpos", "syst_specs_calibpos.root"},
80  // {"calibneg", "syst_specs_calibneg.root"},
81  // {"MuES_onlyMuKEDw", "syst_specs_muonedw.root"},
82  // {"MuES_onlyMuKEUp", "syst_specs_muoneup.root"},
83  {"ckv", "syst_specs_ckv.root"},
84  {"angleup", "syst_specs_angleup.root"},
85  {"angledw", "syst_specs_angledw.root"},
86  {"neutron_Up", "syst_specs_neutron.root"},
87  {"neutron_Dw", "syst_specs_neutron.root"},
88  {"neutron_nom", "syst_specs_neutron_nom.root"}
89 };
90 
91 const std::map<std::string, std::string> beamSysts = {
92  // {"2kA", "beamsysts/syst_specs_2kA.root"},
93  // {"02mmBeamSpotSize", "beamsysts/syst_specs_02mmBeamSpotSize.root"},
94  // {"1mmBeamShiftX", "beamsysts/syst_specs_1mmBeamShiftX.root"},
95  // {"1mmBeamShiftY", "beamsysts/syst_specs_1mmBeamShiftY.root"},
96  // {"3mmHorn1X", "beamsysts/syst_specs_3mmHorn1X.root"},
97  // {"3mmHorn1Y", "beamsysts/syst_specs_3mmHorn1Y.root"},
98  // {"3mmHorn2X", "beamsysts/syst_specs_3mmHorn2X.root"},
99  // {"3mmHorn2Y", "beamsysts/syst_specs_3mmHorn2Y.root"},
100  // {"7mmTargetZ", "beamsysts/syst_specs_7mmTargetZ.root"},
101  // {"MagneticFieldinDecayPipe","beamsysts/syst_specs_MagneticFieldinDecayPipe.root"},
102  // {"1mmHornWater", "beamsysts/syst_specs_1mmHornWater.root"}
103 };
104 
105 const std::vector<std::string> updownlabels = {"_Up", "_Dw"};
106 
107 
108 void make_xs(TDirectory* &dirInMC,TDirectory* &dirInDt, std::string sdOut, TH1* hInFlux, TFile* fileOut, Ratio* inv_pur){
109 
110  std::cout<<"=> In make_xs() "<<sdOut<<std::endl;
111 
112  ///////////////////////////
113  ///Getting all:
114  //1. Efficiency:
115  TDirectory * nuTrueDir = dirInMC->GetDirectory("MCSigNuTrue");
116  if (!nuTrueDir || nuTrueDir == NULL || nuTrueDir == 0)
117  nuTrueDir = dirInMC->GetDirectory("MCSigNuTre"); // Try both titles to fix typo
118  Spectrum* fMCSig = ana::LoadFrom<Spectrum> (dirInMC->GetDirectory("MCSig")).release();
119  Spectrum* fMCSigNuTrue = ana::LoadFrom<Spectrum> (nuTrueDir).release();
120  TH3* h3Deff = ana::ToTH3(*fMCSig , pot, kPOT, angbinsCustom, mukebins, eavailbins);
121  TH3* h3Deff_den = ana::ToTH3(*fMCSigNuTrue, pot, kPOT, angbinsCustom, mukebins, eavailbins);
122  TH3* h3DSel = (TH3*) h3Deff->Clone("h3DSel");
123  TH3* h3DSig = (TH3*) h3Deff_den->Clone("h3DSig");
124  TH2* h2Deff = (TH2*) h3Deff->Project3D("yx");
125  TH2* h2Deff_den = (TH2*) h3Deff_den->Project3D("yx");
126  h3Deff->Divide(h3Deff_den);
127  h2Deff->Divide(h2Deff_den);
128  h3Deff_den->Delete();
129  h2Deff_den->Delete();
130  delete fMCSig;
131  delete fMCSigNuTrue;
132 
133  //2. Signal:
134  Spectrum* fData = ana::LoadFrom<Spectrum> (dirInDt->GetDirectory(fake_xs.c_str())).release();
135  IBkgdEstimator* fBkgEst = ana::LoadFrom<IBkgdEstimator>(dirInMC->GetDirectory("BkgdEst")).release();
136  Spectrum *fBkg = new Spectrum(fBkgEst->Background());
137  Spectrum *fSignalEst = new Spectrum(*fData - *fBkg);
138  delete fBkgEst;
139  delete fBkg;
140 
141  //Background constraint
142  // Ratio fSignalEst_Bkg_Constraint = fData/inv_pur;
143  Ratio * fSignalEst_Bkg_Constraint = new Ratio(*fData, Spectrum(inv_pur->ToTH1(), fData->POT(), fData->Livetime()));
144 
145  //3. Unfolding:
146  ReweightableSpectrum* fRecoTrue = ana::LoadFrom<ReweightableSpectrum>(dirInMC->GetDirectory("RecoTrue")).release();
147 
148  UnfoldIterative * unfold = new UnfoldIterative(*fRecoTrue,5);
149  Spectrum * fSigUnfolded = new Spectrum(unfold->Truth(*fSignalEst));
150  delete unfold;
151 
152  UnfoldIterative * unfold_bkgc = new UnfoldIterative(*fRecoTrue,5);
153  Spectrum * fSigUnfolded_bkgc = new Spectrum(unfold_bkgc->Truth(Spectrum(fSignalEst_Bkg_Constraint->ToTH1(), fData->POT(), fData->Livetime())));
154  delete unfold_bkgc;
155  delete fRecoTrue;
156 
157  //4. Calculation:
158  double integrated_flux = hInFlux->Integral();
159  TH3* h3SigEst[2];
160  h3SigEst[0] = ana::ToTH3(*fSignalEst, pot, kPOT, angbinsCustom, mukebins, eavailbins);
161  h3SigEst[1] = ana::ToTH3(*fSignalEst_Bkg_Constraint, angbinsCustom, mukebins, eavailbins);
162  TH3* h3SigUnfolded[2];
163  h3SigUnfolded[0] = ana::ToTH3(*fSigUnfolded , pot, kPOT, angbinsCustom, mukebins, eavailbins);
164  h3SigUnfolded[1] = ana::ToTH3(*fSigUnfolded_bkgc, pot, kPOT, angbinsCustom, mukebins, eavailbins);
165  delete fSignalEst;
166  delete fSignalEst_Bkg_Constraint;
167  delete fSigUnfolded;
168  delete fSigUnfolded_bkgc;
169 
170  TH2* h2Dsigest[2];
171  TH1* h1Dsigest[2];
172  TH2* h2Dxs[2];
173  TH1* h1Dxs[2];
174  for(int ii=0;ii<2;ii++){
175  h3SigUnfolded[ii]->Divide(h3Deff);
176  h2Dxs[ii] = (TH2*)h3SigUnfolded[ii]->Project3D("yx");
177  h1Dxs[ii] = (TH1*)h3SigUnfolded[ii]->Project3D("z");
178 
179  h2Dxs[ii]->Scale(1./integrated_flux);
180  for(int ibin=1;ibin<=h1Dxs[ii]->GetNbinsX();ibin++){
181  double cont_num = h1Dxs[ii]->GetBinContent(ibin);
182  double cont_den = hInFlux->GetBinContent(ibin);
183  if(cont_den>0)h1Dxs[ii]->SetBinContent(ibin,cont_num/cont_den);
184  else h1Dxs[ii]->SetBinContent(ibin,0);
185  h1Dxs[ii]->SetBinError(ibin,0);
186  }
187  h2Dxs[ii]->Scale(1./ntargs);
188  h1Dxs[ii]->Scale(1./ntargs);
189  }
190 
191  //Storing
192  //xs_nominal: calcualting the signal estomation by S-B
193  fileOut->cd();
194  h3DSel->Write(("h3DSel_" + sdOut).c_str());
195  h3DSig->Write(("h3DSig_" + sdOut).c_str());
196  h3Deff->Write(("h3DEff_" + sdOut).c_str());
197  h2Deff->Write(("h2DEff_" + sdOut).c_str());
198  fileOut->cd(xs_mode0.c_str());
199  h2Dxs[0]->Write(("hXS2Dsec_" + sdOut).c_str());
200  h1Dxs[0]->Write(("hXS1Dsec_" + sdOut).c_str());
201  h3SigUnfolded[0]->Write(("hXS3Dsec_" + sdOut).c_str());
202  h3SigEst[0]->Write(("h3SigEst_" + sdOut).c_str());
203  //xs_bkgcons: calcualting the signal estomation by P*S
204  fileOut->cd();
205  fileOut->cd(xs_mode1.c_str());
206  h2Dxs[1]->Write(("hXS2Dsec_" + sdOut).c_str());
207  h1Dxs[1]->Write(("hXS1Dsec_" + sdOut).c_str());
208  h3SigUnfolded[1]->Write(("hXS3Dsec_" + sdOut).c_str());
209  h3SigEst[1]->Write(("h3SigEst_" + sdOut).c_str());
210 
211  //Leftover cleanup
212  h3DSel->Delete();
213  h3DSig->Delete();
214  h3Deff->Delete();
215  h2Deff->Delete();
216 }
217 
218 void make_xs_1D(TDirectory* &dirInMC,TDirectory* &dirInDt, std::string sdOut, TH1* hInFlux, TFile* fileOut, Ratio* inv_pur){
219  std::cout<<"=> In make_xs_1D() "<< sdOut << std::endl;
220 
221  ///////////////////////////
222  ///Getting all:
223  //1. Efficiency:
224  TDirectory * nuTrueDir = dirInMC->GetDirectory("MCSigNuTrue");
225  if (!nuTrueDir || nuTrueDir == NULL || nuTrueDir == 0)
226  nuTrueDir = dirInMC->GetDirectory("MCSigNuTre"); // Try both titles to fix typo
227  Spectrum* fMCSig = ana::LoadFrom<Spectrum> (dirInMC->GetDirectory("MCSig")).release();
228  Spectrum* fMCSigNuTrue = ana::LoadFrom<Spectrum> (nuTrueDir).release();
229  TH1* h1Deff = fMCSig->ToTH1(pot);
230  TH1* h1Deff_den = fMCSigNuTrue->ToTH1(pot);
231  TH1* h1DSel = (TH1*) h1Deff->Clone("h1DSel");
232  TH1* h1DSig = (TH1*) h1Deff_den->Clone("h1DSig");
233  h1Deff->Divide(h1Deff_den);
234  h1Deff_den->Delete();
235  delete fMCSig;
236  delete fMCSigNuTrue;
237 
238  //2. Signal:
239  Spectrum* fData = ana::LoadFrom<Spectrum> (dirInDt->GetDirectory(fake_xs.c_str())).release();
240  IBkgdEstimator* fBkgEst = ana::LoadFrom<IBkgdEstimator>(dirInMC->GetDirectory("BkgdEst")).release();
241  Spectrum *fBkg = new Spectrum(fBkgEst->Background());
242  Spectrum *fSignalEst = new Spectrum(*fData - *fBkg);
243  delete fBkgEst;
244  delete fBkg;
245 
246  //Background constraint
247  // Ratio fSignalEst_Bkg_Constraint = fData/inv_pur;
248  Ratio * fSignalEst_Bkg_Constraint = new Ratio(*fData, Spectrum(inv_pur->ToTH1(), fData->POT(), fData->Livetime()));
249 
250  //3. Unfolding:
251  ReweightableSpectrum* fRecoTrue = ana::LoadFrom<ReweightableSpectrum>(dirInMC->GetDirectory("RecoTrue")).release();
252 
253  UnfoldIterative * unfold = new UnfoldIterative(*fRecoTrue,5);
254  Spectrum * fSigUnfolded = new Spectrum(unfold->Truth(*fSignalEst));
255  delete unfold;
256 
257  UnfoldIterative * unfold_bkgc = new UnfoldIterative(*fRecoTrue,5);
258  Spectrum * fSigUnfolded_bkgc = new Spectrum(unfold_bkgc->Truth(Spectrum(fSignalEst_Bkg_Constraint->ToTH1(), fData->POT(), fData->Livetime())));
259  delete unfold_bkgc;
260  delete fRecoTrue;
261 
262  //4. Calculation:
263  double integrated_flux = hInFlux->Integral();
264  TH1* h1SigEst[2];
265  h1SigEst[0] = fSignalEst->ToTH1(pot);
266  h1SigEst[1] = fSignalEst_Bkg_Constraint->ToTH1();
267  TH1* h1SigUnfolded[2];
268  h1SigUnfolded[0] = fSigUnfolded->ToTH1(pot);
269  h1SigUnfolded[1] = fSigUnfolded_bkgc->ToTH1(pot);
270  delete fSignalEst;
271  delete fSignalEst_Bkg_Constraint;
272  delete fSigUnfolded;
273  delete fSigUnfolded_bkgc;
274 
275  TH1* h1Dsigest[2];
276  TH1* h1Dxs[2];
277  for(int ii=0;ii<2;ii++){
278  h1SigUnfolded[ii]->Divide(h1Deff);
279  h1Dxs[ii] = h1SigUnfolded[ii];
280 
281  h1Dxs[ii]->Scale(1./integrated_flux);
282  // for(int ibin=1;ibin<=h1Dxs[ii]->GetNbinsX();ibin++){
283  // double cont_num = h1Dxs[ii]->GetBinContent(ibin);
284  // double cont_den = hInFlux->GetBinContent(ibin);
285  // if(cont_den>0)h1Dxs[ii]->SetBinContent(ibin,cont_num/cont_den);
286  // else h1Dxs[ii]->SetBinContent(ibin,0);
287  // h1Dxs[ii]->SetBinError(ibin,0);
288  // }
289  h1Dxs[ii]->Scale(1./ntargs);
290  }
291 
292  //Storing
293  //xs_nominal: calcualting the signal estomation by S-B
294  fileOut->cd();
295  h1DSel->Write(("h1DSel_" + sdOut).c_str());
296  h1DSig->Write(("h1DSig_" + sdOut).c_str());
297  h1Deff->Write(("h1DEff_" + sdOut).c_str());
298  fileOut->cd(xs_mode0.c_str());
299  h1Dxs[0]->Write(("hXS1Dsec_" + sdOut).c_str());
300  h1SigUnfolded[0]->Write(("hXS1Dsec_" + sdOut).c_str());
301  h1SigEst[0]->Write(("h1SigEst_" + sdOut).c_str());
302  //xs_bkgcons: calcualting the signal estomation by P*S
303  fileOut->cd();
304  fileOut->cd(xs_mode1.c_str());
305  h1Dxs[1]->Write(("hXS1Dsec_" + sdOut).c_str());
306  h1SigUnfolded[1]->Write(("hXS1Dsec_" + sdOut).c_str());
307  h1SigEst[1]->Write(("h1SigEst_" + sdOut).c_str());
308 
309  //Leftover cleanup
310  h1DSel->Delete();
311  h1DSig->Delete();
312  h1Deff->Delete();
313 }
314 
315 void Make_NuMuCC_Inc_XS_v2(const std::string inputDir, const std::string fOutname){
316  ::dirIn = inputDir;
317  std::cout<<"=> start "<<std::endl;
318  //this script does the cv, det, muon energy scale and focusing systematics
319 
320  // Load core files
321  TFile* fsys = new TFile((dirIn + "/" + coreFiles["nominal"]).c_str());
322  TFile* fluxFile = new TFile((dirIn + "/" + coreFiles["fluxes"]).c_str());
323  TFile* fakedataFile = new TFile((dirIn + "/" + coreFiles["fakedata"]).c_str());
324  TFile* fmc;
325  assert(fsys && fsys->IsOpen());
326  assert(fluxFile && fluxFile->IsOpen());
327  assert(fakedataFile && fakedataFile->IsOpen());
328 
329  //flux:
330  TH1* hFlux_cv = Spectrum::LoadFrom(fluxFile->GetDirectory("fluxes/CV"))->ToTH1(pot);
331  hFlux_cv->Scale(1e-4); // nu/m^2 to nu/cm^2
332 
333  // Read ppfx
334  TH1* hFlux_ppfx[Nuniv_ppfx];
335  for(int i=0;i<Nuniv_ppfx;i++){
336  hFlux_ppfx[i] = Spectrum::LoadFrom(fluxFile->GetDirectory(Form("fluxes/ppfx%d",i)))->ToTH1(pot);
337  hFlux_ppfx[i]->Scale(1e-4);
338  }
339 
340  //file Out
341  TFile* fOut = new TFile(fOutname.c_str(), "RECREATE");
342  fOut->mkdir(xs_mode0.c_str());
343  fOut->mkdir(xs_mode1.c_str());
344 
345  // Load all the unfolding spectrums we'll need. Load 3D, 1D separately to keep track.
346  std::map<std::string, TDirectory*> dir_data;
347  std::map<std::string, TDirectory*> nuTrueDirs;
348  std::map<std::string, Spectrum*> MCSelected;
349  std::map<std::string, Spectrum*> MCSignal;
350  std::map<std::string, Spectrum*> MCSignal_Reco;
351  std::map<std::string, Spectrum*> MCBkg;
352  std::map<std::string, Ratio*> MC_Inv_Pur;
353  for (std::string var : unfoldVars3D){
354  dir_data[var] = fakedataFile->GetDirectory(var2FD[var].c_str(), true);
355  assert(dir_data[var] && dir_data[var] != NULL);
356  MCSelected[var] = Spectrum::LoadFrom(fsys->GetDirectory(("cv/" + var + "/MCSel_Reco").c_str(), true)).release();
357  MCSignal_Reco[var] = Spectrum::LoadFrom(fsys->GetDirectory(("cv/" + var + "/MCSig_Reco").c_str(), true)).release();
358  nuTrueDirs[var] = fsys->GetDirectory(("cv/" + var + "/MCSigNuTrue").c_str(), true);
359  if (!nuTrueDirs[var] || nuTrueDirs[var] == NULL || nuTrueDirs[var] == 0)
360  nuTrueDirs[var] = fsys->GetDirectory(("cv/" + var + "/MCSigNuTre").c_str(), true); // Fix typo in dir name
361  assert(nuTrueDirs[var] && nuTrueDirs[var] != NULL);
362  MCSignal[var] = Spectrum::LoadFrom(nuTrueDirs[var]).release();
363  MC_Inv_Pur[var] = new Ratio(*MCSelected[var], *MCSignal_Reco[var]);
364  }
365  for (std::string var : unfoldVars1D){
366  dir_data[var] = fakedataFile->GetDirectory(var2FD[var].c_str(), true);
367  assert(dir_data[var] && dir_data[var] != NULL);
368  MCSelected[var] = Spectrum::LoadFrom(fsys->GetDirectory(("cv/" + var + "/MCSel_Reco").c_str(), true)).release();
369  MCSignal_Reco[var] = Spectrum::LoadFrom(fsys->GetDirectory(("cv/" + var + "/MCSig_Reco").c_str(), true)).release();
370  nuTrueDirs[var] = fsys->GetDirectory(("cv/" + var + "/MCSigNuTrue").c_str(), true);
371  if (!nuTrueDirs[var] || nuTrueDirs[var] == NULL || nuTrueDirs[var] == 0)
372  nuTrueDirs[var] = fsys->GetDirectory(("cv/" + var + "/MCSigNuTre").c_str(), true); // Fix typo in dir name
373  assert(nuTrueDirs[var] && nuTrueDirs[var] != NULL);
374  MCSignal[var] = Spectrum::LoadFrom(nuTrueDirs[var]).release();
375  MC_Inv_Pur[var] = new Ratio(*MCSelected[var], *MCSignal_Reco[var]);
376  }
377 
378  // Holder for the current var directory
379  TDirectory* dir_mc;
380  //CV:
381  for(std::string var : unfoldVars3D)
382  {
383  dir_mc = fsys->GetDirectory(("cv/" + var).c_str(), true);
384  assert(dir_mc && dir_mc != NULL);
385  make_xs(dir_mc, dir_data[var], "cv_" + var, hFlux_cv, fOut, MC_Inv_Pur[var]);
386  }
387 
388  for(std::string var : unfoldVars1D)
389  {
390  dir_mc = fsys->GetDirectory(("cv/" + var).c_str(), true);
391  assert(dir_mc && dir_mc != NULL);
392  make_xs_1D(dir_mc, dir_data[var], "cv_" + var, hFlux_cv, fOut, MC_Inv_Pur[var]);
393  }
394 
395  // Systematic files and Muon Energy
396  for (std::pair<std::string, std::string> syst : systFiles){
397  std::string systName = syst.first;
398  std::string systFilename = syst.second;
399  TFile * systFile = new TFile((dirIn + "/" + systFilename).c_str());
400  for(std::string var : unfoldVars3D)
401  {
402  dir_mc = systFile->GetDirectory((systName + "/" + var).c_str(), true);
403  assert(dir_mc && dir_mc != NULL);
404  make_xs(dir_mc, dir_data[var], (systName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
405  }
406 
407  for(std::string var : unfoldVars1D)
408  {
409  dir_mc = systFile->GetDirectory((systName + "/" + var).c_str(), true);
410  assert(dir_mc && dir_mc != NULL);
411  make_xs_1D(dir_mc, dir_data[var], (systName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
412  }
413  systFile->Close();
414  }
415 
416  // Beam Systematics
417  for (std::pair<std::string, std::string> beamsyst : beamSysts){
418  std::string beamsystName = beamsyst.first;
419  std::string beamFilename = beamsyst.second;
420  TFile * beamsystFile = new TFile((dirIn + "/" + beamFilename).c_str());
421 
422  std::string beamsystupName = beamsystName + updownlabels[0];
423  std::string beamsystdwName = beamsystName + updownlabels[1];
424  for(std::string var : unfoldVars3D)
425  {
426  dir_mc = beamsystFile->GetDirectory((beamsystupName + "/" + var).c_str(), true);
427  assert(dir_mc && dir_mc != NULL);
428  make_xs(dir_mc, dir_data[var], (beamsystupName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
429 
430  dir_mc = beamsystFile->GetDirectory((beamsystdwName + "/" + var).c_str(), true);
431  assert(dir_mc && dir_mc != NULL);
432  make_xs(dir_mc, dir_data[var], (beamsystdwName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
433  }
434 
435  for(std::string var : unfoldVars1D)
436  {
437  dir_mc = beamsystFile->GetDirectory((beamsystupName + "/" + var).c_str(), true);
438  assert(dir_mc && dir_mc != NULL);
439  make_xs_1D(dir_mc, dir_data[var], (beamsystupName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
440 
441  dir_mc = beamsystFile->GetDirectory((beamsystdwName + "/" + var).c_str(), true);
442  assert(dir_mc && dir_mc != NULL);
443  make_xs_1D(dir_mc, dir_data[var], (beamsystdwName + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
444  }
445 
446  beamsystFile->Close();
447  }
448 
449  //PPFX
450  TFile * ppfxFile;
451  for(int i=0;i<Nuniv_ppfx;i++){
452  // Only update file every 4 univs (Please can we do this in a smarter way?)
453  int idx_file = i%4;
454  if(idx_file==0)
455  {
456  // Close and delete file if it exists
457  if (ppfxFile && ppfxFile != NULL){
458  ppfxFile->Close();
459  delete ppfxFile;
460  ppfxFile = NULL;
461  }
462 
463  // Replace the file
464  std::string ppfxFilename = dirIn + "/ppfx_univ_" + to_string(i) + "_to_" + to_string(i+3) + ".root";
465  if(boost::filesystem::exists(ppfxFilename.c_str()))
466  ppfxFile = new TFile(ppfxFilename.c_str(), "read");
467 
468  }
469  if (!ppfxFile || ppfxFile == NULL || ppfxFile->IsZombie())
470  continue;
471 
472  // Perform unfolding for all 3D and 1D vars included for this PPFX
473  std::string ppfxLabel = "ppfx" + to_string(i);
474  for(std::string var : unfoldVars3D)
475  {
476  dir_mc = ppfxFile->GetDirectory((ppfxLabel + "/" + var).c_str(), true);
477  assert(dir_mc && dir_mc != NULL);
478  make_xs(dir_mc, dir_data[var], (ppfxLabel + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
479  }
480 
481  for(std::string var : unfoldVars1D)
482  {
483  dir_mc = ppfxFile->GetDirectory((ppfxLabel + "/" + var).c_str(), true);
484  assert(dir_mc && dir_mc != NULL);
485  make_xs_1D(dir_mc, dir_data[var], (ppfxLabel + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
486  }
487  }
488  if (ppfxFile && ppfxFile != NULL){
489  ppfxFile->Close();
490  delete ppfxFile;
491  ppfxFile = NULL;
492  }
493 
494  //GENIE
495  TFile * genieFile;
496  for(int i=0;i<Nuniv_genie;i++){
497  int idx_file = i%4;
498  if(idx_file==0)
499  {
500  // Close and delete file if it exists
501  if (genieFile && genieFile != NULL){
502  genieFile->Close();
503  delete genieFile;
504  genieFile = NULL;
505  }
506 
507  // Replace the file
508  std::string genieFilename = dirIn + "/genie_univ_" + to_string(i) + "_to_" + to_string(i+3) + ".root";
509  if(boost::filesystem::exists(genieFilename.c_str()))
510  genieFile = new TFile(genieFilename.c_str(), "read");
511 
512  }
513  if (!genieFile || genieFile == NULL || genieFile->IsZombie())
514  continue;
515 
516  std::string genieLabel = "genie" + to_string(i);
517  for(std::string var : unfoldVars3D)
518  {
519  dir_mc = genieFile->GetDirectory((genieLabel + "/" + var).c_str(), true);
520  assert(dir_mc && dir_mc != NULL);
521  make_xs(dir_mc, dir_data[var], (genieLabel + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
522  }
523 
524  for(std::string var : unfoldVars1D)
525  {
526  dir_mc = genieFile->GetDirectory((genieLabel + "/" + var).c_str(), true);
527  assert(dir_mc && dir_mc != NULL);
528  make_xs_1D(dir_mc, dir_data[var], (genieLabel + "_" + var), hFlux_cv, fOut, MC_Inv_Pur[var]);
529  }
530  }
531  if (genieFile && genieFile != NULL){
532  genieFile->Close();
533  delete genieFile;
534  genieFile = NULL;
535  }
536 
537  fOut->Close();
538 }
const int Nuniv_genie
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
const int Nuniv_ppfx
Spectrum with the value of a second variable, allowing for reweighting
const std::string fake_xs
TFile * fmc
Definition: bdt_com.C:14
const Binning eavailbins
const double ntargs
virtual Spectrum Background() const =0
const std::map< std::string, std::string > systFiles
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:237
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const std::string xs_mode0
const std::vector< std::string > unfoldVars1D
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const std::map< std::string, std::string > beamSysts
const std::string xs_mode1
std::string dirIn
const std::vector< std::string > updownlabels
std::vector< float > Spectrum
Definition: Constants.h:610
double POT() const
Definition: Spectrum.h:227
Represent the ratio between two spectra.
Definition: Ratio.h:8
OStream cout
Definition: OStream.cxx:6
virtual Spectrum Truth(const Spectrum &reco) override
void Make_NuMuCC_Inc_XS_v2(const std::string inputDir, const std::string fOutname)
const Binning mukebins
Definition: NumuCCIncBins.h:90
std::map< std::string, std::string > var2FD
const double pot
void make_xs_1D(TDirectory *&dirInMC, TDirectory *&dirInDt, std::string sdOut, TH1 *hInFlux, TFile *fileOut, Ratio *inv_pur)
assert(nhit_max >=nhit_nbins)
std::unique_ptr< IBkgdEstimator > LoadFrom< IBkgdEstimator >(TDirectory *dir, const std::string &label)
std::map< std::string, std::string > coreFiles
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Float_t e
Definition: plot.C:35
const Binning angbinsCustom
Definition: NumuCCIncBins.h:72
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Same as ToTH2, but with 3 dimensions.
Definition: UtilsExt.cxx:162
const std::vector< std::string > unfoldVars3D
void make_xs(TDirectory *&dirInMC, TDirectory *&dirInDt, std::string sdOut, TH1 *hInFlux, TFile *fileOut, Ratio *inv_pur)
enum BeamMode string