wrong_sign_uncertainty.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Analysis/Style.h"
12 #include "CAFAna/Systs/BeamSysts.h"
13 
14 
15 #include "TCanvas.h"
16 #include "TStyle.h"
17 #include "TLegend.h"
18 #include "TLegendEntry.h"
19 #include "TPaletteAxis.h"
20 #include "THStack.h"
21 
26 // temporary measure to use prod4 trained muonid
27 #include "NDAna/muonid/NDXSecMuonPID.h"
28 namespace nuebarccinc {
29  const ana::Cut kMuonIDProd4Cut = ana::kMuonIDProd4 < -0.2;
30 }
31 
33 {
34  const int ret = gSystem->Load("libNDAnamuonid");
35  if(ret != 0) exit(1);
36 }
37 
38 using namespace ana;
39 TH2 * ToTH2(const TH1 * h, const HistAxis & xaxis, const HistAxis & yaxis);
40 TH1 * MostConservative(TH1 * up, TH1 * down, TH1* nominal);
41 void PlotMultiverse(Multiverse & mv, TH1 * nominal, std::string title, std::string name, TH1 * error = 0);
42 void Plot1D(TH1* h1, std::string h1_label, std::map<std::string, TH1*> systs, std::vector<std::string> to_plot, std::string name, double xmax=-1);
43 void Plot2D(TH1* total, std::map<std::string, TH1*> systs, std::vector<std::string> to_plot, std::string basename);
44 void abs(TH1* hist);
45 void MovePalette(TH1 * hist);
46 
47 std::vector<std::string> all_systs = {"ppfx", "genie", "calib_shape", "calib_shift", "ll", "cherenkov", "beam_focusing", "stat"};
48 std::vector<std::string> beam_systs = {"ppfx", "beam_focusing"};
49 
50 std::map<std::string, nuebarccinc::Style> syst_styles = {{"genie" , {"", kCyan , 1}},
51  {"ppfx" , {"", kMagenta , 1}},
52  {"ll" , {"", kOrange , 1}},
53  {"calib_shift" , {"", kGreen , 1}},
54  {"cherenkov", {"", kBlue+1 , 1}},
55  {"ll_up" , {"", kOrange , 9}},
56  {"ll_dw" , {"", kOrange , 2}},
57  {"calib_shift_up" , {"", kGreen , 9}},
58  {"calib_shift_dw" , {"", kGreen , 2}},
59  {"calib_shape" , {"", kRed-4 , 1}},
60  {"stat" , {"", kBlack , 2}},
61 
62  {"2kA" , {"", kCyan, 1}},
63  {"02mmBeamSpotSize", {"", kCyan + 1, 1}},
64  {"1mmBeamShiftX", {"", kMagenta + 2, 1}},
65  {"1mmBeamShiftY", {"", kMagenta+2, 2}},
66  {"3mmHorn1X", {"", kViolet, 1}},
67  {"3mmHorn1Y", {"", kViolet, 2}},
68  {"3mmHorn2X", {"", kAzure, 1}},
69  {"3mmHorn2Y", {"", kAzure, 2}},
70  {"7mmTargetZ", {"", kBlue-6, 1}},
71  {"MagneticFieldinDecayPipe", {"", kGreen+2, 1}},
72  {"1mmHornWater", {"", kYellow-1, 1}},
73  {"beam_focusing", {"", kViolet+1}}
74 
75 };
76 
77 
78 
79 std::map<std::string, double> xmax = {{"electron_kin", -1.},
80  {"enu" , 6.},
81  {"q2" , 2.}};
82 
83 std::map<string, const ana::BeamSyst*> beamSystematics = {
84  {"2kA", &ana::kBeamHornCurrent},
85  {"02mmBeamSpotSize", &ana::kBeamSpotSize},
86  {"1mmBeamShiftX", &ana::kBeamPosX},
87  {"1mmBeamShiftY", &ana::kBeamPosY},
88  {"3mmHorn1X", &ana::kBeamH1PosX},
89  {"3mmHorn1Y", &ana::kBeamH1PosY},
90  {"3mmHorn2X", &ana::kBeamH2PosX},
91  {"3mmHorn2Y", &ana::kBeamH2PosY},
92  {"7mmTargetZ", &ana::kBeamTarget},
93  {"MagneticFieldinDecayPipe", &ana::kBeamMagField},
94  {"1mmHornWater", &ana::kBeamGeomWater}
95 };
96 
98  std::string input_file_name = "/nova/ana/users/ddoyle/NuebarCCInc/WrongSign/wrong_sign_uncertainty.root",
99  std::string plot_dump = "/nova/ana/users/ddoyle/NuebarCCInc/WrongSign/plots")
100 {
102  std::map<std::string, const HistAxis*> axes;
103  axes["electron_kin"] = new HistAxis(nuebarccinc::kTrueElectronEVsCosStandardAxis);
105  axes["q2" ] = new HistAxis(nuebarccinc::kTrueQ2StandardAxis );
106 
107  std::map<std::string, const Cut*> cuts;
108  cuts["nue" ] = new Cut(nuebarccinc::kNueCC );
109  cuts["nuebar" ] = new Cut(nuebarccinc::kNuebarCC);
110  cuts["nunubar"] = new Cut(nuebarccinc::kNuebarCC || nuebarccinc::kNueCC);
111 
112  std::map<std::string, SpectrumLoader*> loaders;
113  loaders["nominal" ] = 0;
114  loaders["calib_shift_up"] = 0;
115  loaders["calib_shift_dw"] = 0;
116  loaders["ll_up"] = 0;
117  loaders["ll_dw"] = 0;
118  loaders["cherenkov"] = 0;
119  loaders["calib_shape"] = 0;
120  if(!make_plots) {
121  loaders["nominal" ] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_NOMINAL);
122  loaders["calib_shift_up"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_UP);
123  loaders["calib_shift_dw"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_DOWN);
124  loaders["calib_shape"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_SHAPE);
125  loaders["ll_up"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_LIGHT_UP);
127  loaders["cherenkov"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CHERENKOV);
128  const Var kCVWeight = kPPFXFluxCVWgt * kXSecCVWgt2020;
129 
130  const Cut kSelection =
133 
134  std::map<std::string, std::map<std::string, std::map<std::string, Spectrum*> > > spectra;
135  for(auto load_it = loaders.begin(); load_it != loaders.end(); load_it++) {
136  load_it->second->SetSpillCut(kStandardSpillCuts);
137  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
138  for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
139  spectra[load_it->first][cut_it->first][ax_it->first] = new Spectrum(*load_it->second,
140  *ax_it->second,
141  kSelection && *cut_it->second,
142  kNoShift,
143  kCVWeight);
144  }
145  }
146  }
147 
148  // beam systs
149  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
150  for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
151  for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
152  auto name_up = beam_it->first + "_up";
153  auto name_dw = beam_it->first + "_dw";
154 
155  spectra[name_up][cut_it->first][ax_it->first] = new Spectrum(*loaders["nominal"],
156  *ax_it->second,
157  kSelection && *cut_it->second,
158  ana::SystShifts(beam_it->second, +1),
159  kCVWeight);
160 
161  spectra[name_dw][cut_it->first][ax_it->first] = new Spectrum(*loaders["nominal"],
162  *ax_it->second,
163  kSelection && *cut_it->second,
164  ana::SystShifts(beam_it->second, -1),
165  kCVWeight);
166  }
167  }
168  }
169 
170  const std::vector<Var> ppfx = GetkPPFXFluxUnivWgt();
171  GenieMultiverseParameters genie_multiverse(100, getAllXsecSysts_2020());
172  std::vector<SystShifts> genie_shifts = genie_multiverse.GetSystShifts();
173 
174  std::map<std::string, std::map<std::string, Multiverse *> > mv_ppfx;
175  std::map<std::string, std::map<std::string, Multiverse *> > mv_genie;
176  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
177  for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
178  mv_ppfx[cut_it->first][ax_it->first] = new Multiverse(*loaders["nominal"],
179  *ax_it->second,
180  kSelection && *cut_it->second,
181  kNoShift,
182  ppfx,
183  kXSecCVWgt2020);
184  mv_genie[cut_it->first][ax_it->first] = new Multiverse(*loaders["nominal"],
185  *ax_it->second,
186  kSelection && *cut_it->second,
187  genie_shifts,
189  }
190  }
191 
192  for(auto load_it = loaders.begin(); load_it != loaders.end(); load_it++)
193  load_it->second->Go();
194 
195 
196  TFile* output = new TFile("wrong_sign_uncertainty.root", "recreate");
197  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
198  for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
199  mv_ppfx .at(cut_it->first).at(ax_it->first)->SaveTo(output, ("ppfx_" + cut_it->first + "_" + ax_it->first).c_str());
200  mv_genie.at(cut_it->first).at(ax_it->first)->SaveTo(output, ("genie_" + cut_it->first + "_" + ax_it->first).c_str());
201  for(auto load_it = loaders.begin(); load_it != loaders.end(); load_it++) {
202  spectra.at(load_it->first).at(cut_it->first).at(ax_it->first)->SaveTo(output, TString::Format("%s_%s_%s",
203  load_it->first.c_str(),
204  cut_it ->first.c_str(),
205  ax_it ->first.c_str()).Data());
206  }
207  for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
208  spectra.at(beam_it->first + "_up").at(cut_it->first).at(ax_it->first)->SaveTo(output, TString::Format("%s_%s_%s",
209  cut_it ->first.c_str(),
210  ax_it ->first.c_str()).Data());
211  spectra.at(beam_it->first + "_dw").at(cut_it->first).at(ax_it->first)->SaveTo(output, TString::Format("%s_%s_%s",
212  (beam_it->first + "_dw").c_str(),
213  cut_it ->first.c_str(),
214  ax_it ->first.c_str()).Data());
215  }
216  }
217  }
218 
219 
220 
221  output->Close();
222  }
223  else {
224  TFile * input = TFile::Open(input_file_name.c_str());
225 
226  double MCPOT = ((TH1*) input->Get("nominal_nue_q2/pot"))->GetBinContent(1);
227  std::map<std::string, double> pots;
228  std::map<std::string, std::map<std::string, Multiverse *> > mv_ppfx;
229  std::map<std::string, std::map<std::string, Multiverse *> > mv_genie;
230  std::map<std::string, std::map<std::string, std::map<std::string, TH1*> > > hists;
231  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
232  for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
233  mv_ppfx [ax_it->first][cut_it->first] = Multiverse::LoadFrom(input, ("ppfx_" + cut_it->first + "_" + ax_it->first).c_str())).release(;
234  mv_genie[ax_it->first][cut_it->first] = Multiverse::LoadFrom(input, ("genie_" + cut_it->first + "_" + ax_it->first).c_str())).release(;
235 
236  mv_ppfx .at(ax_it->first).at(cut_it->first)->Scale(kAna2020RHCPOT / MCPOT);
237  mv_genie.at(ax_it->first).at(cut_it->first)->Scale(kAna2020RHCPOT / MCPOT);
238 
239  for(auto load_it = loaders.begin(); load_it != loaders.end(); load_it++) {
240  auto spec = Spectrum::LoadFrom(input, TString::Format("%s_%s_%s",
241  load_it->first.c_str(),
242  cut_it ->first.c_str(),
243  ax_it ->first.c_str()).Data());
244  pots[load_it->first] = spec->POT();
245  auto pot = kAna2020RHCPOT;//spec->POT();
246  hists[ax_it->first][cut_it->first][load_it->first] = spec->ToTH1(pot);
247  }
248 
249  for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
250  auto spec_up = Spectrum::LoadFrom(input, TString::Format("%s_%s_%s",
251  (beam_it->first + "_up").c_str(),
252  cut_it ->first .c_str(),
253  ax_it ->first .c_str()).Data());
254  auto spec_dw = Spectrum::LoadFrom(input, TString::Format("%s_%s_%s",
255  (beam_it->first + "_dw").c_str(),
256  cut_it ->first .c_str(),
257  ax_it ->first .c_str()).Data());
258 
259  hists[ax_it->first][cut_it->first][beam_it->first + "_up"] = spec_up->ToTH1(kAna2020RHCPOT);
260  hists[ax_it->first][cut_it->first][beam_it->first + "_dw"] = spec_dw->ToTH1(kAna2020RHCPOT);
261  }
262  }
263  }
264 
265  // make fractions
266  std::map<std::string, std::map<std::string, TH1*> > nuebar_fractions;
267  std::map<std::string, Multiverse*> mv_ppfx_nuebar_fractions;
268  std::map<std::string, Multiverse*> mv_genie_nuebar_fractions;
269  std::map<std::string, Multiverse*> mv_ppfx_nue_nuebar_ratios;
270  std::map<std::string, Multiverse*> mv_genie_nue_nuebar_ratios;
271  std::map<std::string, std::map<std::string, TH1*> > nue_nuebar_ratios;
272  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
273  mv_ppfx_nuebar_fractions [ax_it->first] = new Multiverse(*mv_ppfx .at(ax_it->first).at("nuebar"));
274  mv_genie_nuebar_fractions[ax_it->first] = new Multiverse(*mv_genie.at(ax_it->first).at("nuebar"));
275 
276  *mv_ppfx_nuebar_fractions.at(ax_it->first) /= *mv_ppfx .at(ax_it->first).at("nunubar");
277  *mv_genie_nuebar_fractions.at(ax_it->first) /= *mv_genie.at(ax_it->first).at("nunubar");
278 
279  mv_ppfx_nue_nuebar_ratios[ax_it->first] = new Multiverse(*mv_ppfx.at(ax_it->first).at("nue"));
280  *mv_ppfx_nue_nuebar_ratios.at(ax_it->first) /= *mv_ppfx.at(ax_it->first).at("nuebar");
281 
282  mv_genie_nue_nuebar_ratios[ax_it->first] = new Multiverse(*mv_genie.at(ax_it->first).at("nue"));
283  *mv_genie_nue_nuebar_ratios.at(ax_it->first) /= *mv_genie.at(ax_it->first).at("nuebar");
284 
285  for(auto loader_it = loaders.begin(); loader_it != loaders.end(); loader_it++) {
286  nuebar_fractions[ax_it->first][loader_it->first] = (TH1*) hists.at(ax_it->first).at("nuebar").at(loader_it->first)->Clone(UniqueName().c_str());
287  nuebar_fractions.at(ax_it->first).at(loader_it->first)->Divide(hists.at(ax_it->first).at("nunubar").at(loader_it->first));
288 
289  nue_nuebar_ratios[ax_it->first][loader_it->first] = (TH1*) hists.at(ax_it->first).at("nue").at(loader_it->first)->Clone(UniqueName().c_str());
290  nue_nuebar_ratios.at(ax_it->first).at(loader_it->first)->Divide(hists.at(ax_it->first).at("nuebar").at(loader_it->first));
291  }
292 
293  // two sided file systs
294  for(std::string syst : {"calib_shift", "ll"}) {
295  TH1 * tmp_up = (TH1*) hists.at(ax_it->first).at("nuebar").at(syst + "_up")->Clone(UniqueName().c_str());
296  TH1 * tmp_dw = (TH1*) hists.at(ax_it->first).at("nuebar").at(syst + "_dw")->Clone(UniqueName().c_str());
297  tmp_up->Divide(hists.at(ax_it->first).at("nunubar").at(syst + "_up"));
298  tmp_dw->Divide(hists.at(ax_it->first).at("nunubar").at(syst + "_dw"));
299 
300  nuebar_fractions[ax_it->first][syst] = MostConservative(tmp_up, tmp_dw, nuebar_fractions.at(ax_it->first).at("nominal"));
301 
302  delete tmp_up;
303  delete tmp_dw;
304 
305  tmp_up = (TH1*) hists.at(ax_it->first).at("nue").at(syst + "_up")->Clone(UniqueName().c_str());
306  tmp_dw = (TH1*) hists.at(ax_it->first).at("nue").at(syst + "_dw")->Clone(UniqueName().c_str());
307  tmp_up->Divide(hists.at(ax_it->first).at("nuebar").at(syst + "_up"));
308  tmp_dw->Divide(hists.at(ax_it->first).at("nuebar").at(syst + "_dw"));
309 
310  nue_nuebar_ratios[ax_it->first][syst] = MostConservative(tmp_up, tmp_dw, nue_nuebar_ratios.at(ax_it->first).at("nominal"));
311  }
312 
313  // treat all beam systs as two sided even though connor said there are some that are one sided.
314  for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
315  TH1 * tmp_up = (TH1*) hists.at(ax_it->first).at("nuebar").at(beam_it->first + "_up")->Clone(UniqueName().c_str());
316  TH1 * tmp_dw = (TH1*) hists.at(ax_it->first).at("nuebar").at(beam_it->first + "_dw")->Clone(UniqueName().c_str());
317  tmp_up->Divide(hists.at(ax_it->first).at("nunubar").at(beam_it->first + "_up"));
318  tmp_dw->Divide(hists.at(ax_it->first).at("nunubar").at(beam_it->first + "_dw"));
319 
320  nuebar_fractions[ax_it->first][beam_it->first] = MostConservative(tmp_up, tmp_dw, nuebar_fractions.at(ax_it->first).at("nominal"));
321 
322  delete tmp_up;
323  delete tmp_dw;
324 
325  tmp_up = (TH1*) hists.at(ax_it->first).at("nue") .at(beam_it->first + "_up")->Clone(UniqueName().c_str());
326  tmp_dw = (TH1*) hists.at(ax_it->first).at("nue") .at(beam_it->first + "_dw")->Clone(UniqueName().c_str());
327  tmp_up->Divide(hists.at(ax_it->first).at("nuebar").at(beam_it->first + "_up"));
328  tmp_dw->Divide(hists.at(ax_it->first).at("nuebar").at(beam_it->first + "_dw"));
329 
330  nue_nuebar_ratios[ax_it->first][beam_it->first] = MostConservative(tmp_up, tmp_dw, nue_nuebar_ratios.at(ax_it->first).at("nominal"));
331 
332  }
333 
334  // cherenkov
335  nuebar_fractions.at(ax_it->first)["cherenkov"] = (TH1*) hists.at(ax_it->first).at("nuebar").at("cherenkov")->Clone(UniqueName().c_str());
336  nuebar_fractions.at(ax_it->first).at("cherenkov")->Divide(hists.at(ax_it->first).at("nunubar").at("cherenkov"));
337 
338  nue_nuebar_ratios.at(ax_it->first)["cherenkov"] = (TH1*) hists.at(ax_it->first).at("nue").at("cherenkov")->Clone(UniqueName().c_str());
339  nue_nuebar_ratios.at(ax_it->first).at("cherenkov")->Divide(hists.at(ax_it->first).at("nuebar").at("cherenkov"));
340 
341  // calib shape
342  nuebar_fractions.at(ax_it->first)["calib_shape"] = (TH1*) hists.at(ax_it->first).at("nuebar").at("calib_shape")->Clone(UniqueName().c_str());
343  nuebar_fractions.at(ax_it->first).at("calib_shape")->Divide(hists.at(ax_it->first).at("nunubar").at("calib_shape"));
344 
345  nue_nuebar_ratios.at(ax_it->first)["calib_shape"] = (TH1*) hists.at(ax_it->first).at("nue").at("calib_shape")->Clone(UniqueName().c_str());
346  nue_nuebar_ratios.at(ax_it->first).at("calib_shape")->Divide(hists.at(ax_it->first).at("nuebar").at("calib_shape"));
347 
348 
349  // multiverses
350  TH1 * ppfx_up = mv_ppfx_nuebar_fractions.at(ax_it->first)->GetPlusOneSigmaShift (nuebar_fractions.at(ax_it->first).at("nominal"));
351  TH1 * ppfx_dw = mv_ppfx_nuebar_fractions.at(ax_it->first)->GetMinusOneSigmaShift(nuebar_fractions.at(ax_it->first).at("nominal"));
352 
353  TH1 * genie_up = mv_genie_nuebar_fractions.at(ax_it->first)->GetPlusOneSigmaShift (nuebar_fractions.at(ax_it->first).at("nominal"));
354  TH1 * genie_dw = mv_genie_nuebar_fractions.at(ax_it->first)->GetMinusOneSigmaShift(nuebar_fractions.at(ax_it->first).at("nominal"));
355 
356  nuebar_fractions[ax_it->first]["ppfx"] = MostConservative(ppfx_up , ppfx_dw , nuebar_fractions.at(ax_it->first).at("nominal"));
357  nuebar_fractions[ax_it->first]["genie"] = MostConservative(genie_up, genie_dw, nuebar_fractions.at(ax_it->first).at("nominal"));
358 
359 
360  delete ppfx_up;
361  delete ppfx_dw;
362  delete genie_up;
363  delete genie_dw;
364 
365  ppfx_up = mv_ppfx_nue_nuebar_ratios.at(ax_it->first)->GetPlusOneSigmaShift (nue_nuebar_ratios.at(ax_it->first).at("nominal"));
366  ppfx_dw = mv_ppfx_nue_nuebar_ratios.at(ax_it->first)->GetMinusOneSigmaShift(nue_nuebar_ratios.at(ax_it->first).at("nominal"));
367  nue_nuebar_ratios.at(ax_it->first)["ppfx"] = MostConservative(ppfx_up, ppfx_dw, nue_nuebar_ratios.at(ax_it->first).at("nominal"));
368 
369  genie_up = mv_genie_nue_nuebar_ratios.at(ax_it->first)->GetPlusOneSigmaShift (nue_nuebar_ratios.at(ax_it->first).at("nominal"));
370  genie_dw = mv_genie_nue_nuebar_ratios.at(ax_it->first)->GetMinusOneSigmaShift(nue_nuebar_ratios.at(ax_it->first).at("nominal"));
371  nue_nuebar_ratios.at(ax_it->first)["genie"] = MostConservative(genie_up, genie_dw, nue_nuebar_ratios.at(ax_it->first).at("nominal"));
372  }
373 
374  for(auto it1 = nuebar_fractions.begin(); it1 != nuebar_fractions.end(); it1++)
375  for(auto it2 = nuebar_fractions.at(it1->first).begin(); it2 != nuebar_fractions.at(it1->first).end(); it2++)
376  it2->second->GetYaxis()->SetTitle("#bar{#nu_{e}} Fraction");
377 
378  // uncertainties
379  std::map<std::string, std::map<std::string, TH1*> > nuebar_fractions_uncert_abs;
380  std::map<std::string, std::map<std::string, TH1*> > nuebar_fractions_uncert_frac;
381  std::map<std::string, std::map<std::string, TH1*> > nue_nuebar_ratios_uncert_abs;
382  std::map<std::string, std::map<std::string, TH1*> > nue_nuebar_ratios_uncert_frac;
383  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
384  for(std::string syst : {"ppfx", "genie", "calib_shift", "calib_shape", "ll", "cherenkov"} ) {
385  nuebar_fractions_uncert_abs[ax_it->first][syst] = (TH1*) nuebar_fractions.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
386  nuebar_fractions_uncert_abs[ax_it->first][syst]->Add(nuebar_fractions.at(ax_it->first).at("nominal"), -1);
387  abs(nuebar_fractions_uncert_abs[ax_it->first][syst]);
388 
389  nuebar_fractions_uncert_frac[ax_it->first][syst] = (TH1*) nuebar_fractions_uncert_abs.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
390  nuebar_fractions_uncert_frac.at(ax_it->first).at(syst)->Divide(nuebar_fractions.at(ax_it->first).at("nominal"));
391 
392 
393  // ratios
394  nue_nuebar_ratios_uncert_abs[ax_it->first][syst] = (TH1*) nue_nuebar_ratios.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
395  nue_nuebar_ratios_uncert_abs[ax_it->first][syst]->Add(nue_nuebar_ratios.at(ax_it->first).at("nominal"), -1);
396  abs(nue_nuebar_ratios_uncert_abs[ax_it->first][syst]);
397 
398  nue_nuebar_ratios_uncert_frac[ax_it->first][syst] = (TH1*) nue_nuebar_ratios_uncert_abs.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
399  nue_nuebar_ratios_uncert_frac.at(ax_it->first).at(syst)->Divide(nue_nuebar_ratios.at(ax_it->first).at("nominal"));
400  }
401 
402  // beam uncertainties
403  for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
404  std::string syst = beam_it->first;
405  nuebar_fractions_uncert_abs[ax_it->first][syst] = (TH1*) nuebar_fractions.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
406  nuebar_fractions_uncert_abs[ax_it->first][syst]->Add(nuebar_fractions.at(ax_it->first).at("nominal"), -1);
407  abs(nuebar_fractions_uncert_abs[ax_it->first][syst]);
408 
409  nuebar_fractions_uncert_frac[ax_it->first][syst] = (TH1*) nuebar_fractions_uncert_abs.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
410  nuebar_fractions_uncert_frac.at(ax_it->first).at(syst)->Divide(nuebar_fractions.at(ax_it->first).at("nominal"));
411 
412  // ratios
413  nue_nuebar_ratios_uncert_abs[ax_it->first][syst] = (TH1*) nue_nuebar_ratios.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
414  nue_nuebar_ratios_uncert_abs[ax_it->first][syst]->Add(nue_nuebar_ratios.at(ax_it->first).at("nominal"), -1);
415  abs(nue_nuebar_ratios_uncert_abs[ax_it->first][syst]);
416 
417  nue_nuebar_ratios_uncert_frac[ax_it->first][syst] = (TH1*) nue_nuebar_ratios_uncert_abs.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
418  nue_nuebar_ratios_uncert_frac.at(ax_it->first).at(syst)->Divide(nue_nuebar_ratios.at(ax_it->first).at("nominal"));
419  }
420 
421  // combine beam uncertainties into one
422  nuebar_fractions_uncert_frac.at(ax_it->first)["beam_focusing"] = (TH1*) nuebar_fractions_uncert_frac.at(ax_it->first).at("2kA")->Clone(UniqueName().c_str());
423  nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->Reset();
424 
425  nue_nuebar_ratios_uncert_frac.at(ax_it->first)["beam_focusing"] = (TH1*) nue_nuebar_ratios_uncert_frac.at(ax_it->first).at("2kA")->Clone(UniqueName().c_str());
426  nue_nuebar_ratios_uncert_frac.at(ax_it->first).at("beam_focusing")->Reset();
427 
428  auto NX = nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->GetNbinsX();
429  auto NY = nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->GetNbinsY();
430  auto NZ = nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->GetNbinsZ();
431  for(auto x = 1; x <= NX; x++) {
432  for(auto y = 1; y <= NY; y++) {
433  for(auto z = 1; z <= NZ; z++) {
434  double rat = 0;
435  double val = 0;
436  for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
437  rat += std::pow(nue_nuebar_ratios_uncert_frac.at(ax_it->first).at(beam_it->first)->GetBinContent(x, y, z), 2);
438  val += std::pow(nuebar_fractions_uncert_frac.at(ax_it->first).at(beam_it->first)->GetBinContent(x, y, z), 2);
439  }
440  nue_nuebar_ratios_uncert_frac.at(ax_it->first).at("beam_focusing")->SetBinContent(x, y, z, std::sqrt(rat));
441  nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->SetBinContent(x, y, z, std::sqrt(val));
442  }
443  }
444  }
445 
446  // MC statistics
447  nuebar_fractions_uncert_abs .at(ax_it->first)["stat"] = (TH1*) nuebar_fractions.at(ax_it->first).at("nominal")->Clone(UniqueName().c_str());
448  nuebar_fractions_uncert_frac.at(ax_it->first)["stat"] = (TH1*) nuebar_fractions.at(ax_it->first).at("nominal")->Clone(UniqueName().c_str());
449 
450  NX = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsX();
451  NY = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsY();
452  NZ = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsZ();
453  for(auto ix = 1; ix <= NX; ix++) {
454  for(auto iy = 1; iy <= NY; iy++) {
455  for(auto iz = 1; iz <= NZ; iz++) {
456  // rescale back to full MC statistics
457 
458  auto num = hists.at(ax_it->first).at("nuebar") .at("nominal")->GetBinContent(ix, iy, iz) * MCPOT / kAna2020RHCPOT;
459  auto den = hists.at(ax_it->first).at("nunubar").at("nominal")->GetBinContent(ix, iy, iz) * MCPOT / kAna2020RHCPOT;
460  auto err = std::sqrt(1/num + 1/den);
461  auto nom = num / den;
462  nuebar_fractions_uncert_abs .at(ax_it->first).at("stat")->SetBinContent(ix, iy, iz, nom * err);
463  nuebar_fractions_uncert_frac.at(ax_it->first).at("stat")->SetBinContent(ix, iy, iz, err );
464 
465  }
466  }
467  }
468  }
469 
470  // sum in quadrature
471  std::map<std::string, TH1*> nsel_uncert;
472  std::map<std::string, TH1*> mc_frac_uncert;
473  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
474  nsel_uncert [ax_it->first] = (TH1*) hists.at(ax_it->first).at("nuebar").at("nominal")->Clone(UniqueName().c_str());
475  mc_frac_uncert[ax_it->first] = (TH1*) nuebar_fractions.at(ax_it->first).at("nominal")->Clone(UniqueName().c_str());
476  nsel_uncert .at(ax_it->first)->Reset();
477  mc_frac_uncert.at(ax_it->first)->Reset();
478  mc_frac_uncert.at(ax_it->first)->SetTitle("#bar{#nu}_{e} Fraction Uncertainty");
479  nsel_uncert .at(ax_it->first)->SetTitle("Selected #bar{#nu}_{e} Uncertainty");
480 
481  if(ax_it->first == "electron_kin")
482  mc_frac_uncert.at(ax_it->first)->GetZaxis()->SetTitle("Fractional Uncertainty");
483  else
484  mc_frac_uncert.at(ax_it->first)->GetYaxis()->SetTitle("Fractional Uncertainty");
485 
486  auto NX = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsX();
487  auto NY = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsY();
488  auto NZ = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsZ();
489 
490  for(auto ix = 1; ix <= NX; ix++) {
491  for(auto iy = 1; iy <= NY; iy++) {
492  for(auto iz = 1; iz <= NZ; iz++) {
493  double val = 0;
494  for(std::string syst : {"ppfx", "genie", "calib_shape", "calib_shift", "ll", "cherenkov", "beam_focusing", "stat"}) {
495  val += std::pow(nuebar_fractions_uncert_frac.at(ax_it->first).at(syst)->GetBinContent(ix, iy, iz), 2);
496  }// syst
497  mc_frac_uncert.at(ax_it->first)->SetBinContent(ix, iy, iz, std::sqrt(val));
498 
499  auto n = hists.at(ax_it->first).at("nunubar").at("nominal")->GetBinContent(ix, iy, iz);
500  nsel_uncert.at(ax_it->first)->SetBinContent(ix, iy, iz, std::sqrt(1/n + val));
501 
502  } // iz
503  } // iy
504  } // ix
505  } // ax_it
506 
507  mc_frac_uncert.at("enu")->GetXaxis()->SetRangeUser(0, 6);
508  mc_frac_uncert.at("q2" )->GetXaxis()->SetRangeUser(0, 2);
509  nue_nuebar_ratios.at("enu").at("nominal")->GetXaxis()->SetRangeUser(0, 6);
510  nue_nuebar_ratios.at("q2").at("nominal") ->GetXaxis()->SetRangeUser(0, 2);
511 
512  Plot1D(mc_frac_uncert.at("enu"), "Total", nuebar_fractions_uncert_frac.at("enu"), all_systs, plot_dump + "/enu_fractional_uncertainties_all.pdf");
513  Plot1D(mc_frac_uncert.at("q2") , "Total", nuebar_fractions_uncert_frac.at("q2") , all_systs, plot_dump + "/q2_fractional_uncertainties_all.pdf" );
514  Plot1D(mc_frac_uncert.at("enu"), "Total", nuebar_fractions_uncert_frac.at("enu"), beam_systs, plot_dump + "/enu_fractional_uncertainties_beam.pdf");
515  Plot1D(mc_frac_uncert.at("q2") , "Total", nuebar_fractions_uncert_frac.at("q2") , beam_systs, plot_dump + "/q2_fractional_uncertainties_beam.pdf" );
516  for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
517  // Plot1D(nuebar_fractions .at(ax_it->first).at("nominal"), "Nominal", nuebar_fractions.at(ax_it->first) , all_systs, plot_dump + "/" + ax_it->first + "_nuebar_fraction_all.pdf" , xmax.at(ax_it->first));
518  // Plot1D(nue_nuebar_ratios.at(ax_it->first).at("nominal"), "Nominal", nue_nuebar_ratios.at(ax_it->first), all_systs, plot_dump + "/" + ax_it->first + "_nue_nuebar_ratios_all.pdf", xmax.at(ax_it->first));
519  //
520  // Plot1D(nuebar_fractions .at(ax_it->first).at("nominal"), "Nominal", nuebar_fractions.at(ax_it->first) , beam_systs, plot_dump + "/" + ax_it->first + "_nuebar_fraction_beam.pdf" , xmax.at(ax_it->first));
521  // Plot1D(nue_nuebar_ratios.at(ax_it->first).at("nominal"), "Nominal", nue_nuebar_ratios.at(ax_it->first), beam_systs, plot_dump + "/" + ax_it->first + "_nue_nuebar_ratios_beam.pdf", xmax.at(ax_it->first));
522 
523  // Plot1D(hists.at(ax_it->first).at("nuebar" ).at("nominal"), "Nominal", hists.at(ax_it->first).at("nuebar" ), plot_dump + "/" + ax_it->first + "_nuebar_events.pdf" , xmax.at(ax_it->first));
524  // Plot1D(hists.at(ax_it->first).at("nue" ).at("nominal"), "Nominal", hists.at(ax_it->first).at("nue" ), plot_dump + "/" + ax_it->first + "_nue_events.pdf" , xmax.at(ax_it->first));
525  // Plot1D(hists.at(ax_it->first).at("nunubar").at("nominal"), "Nominal", hists.at(ax_it->first).at("nunubar"), plot_dump + "/" + ax_it->first + "_nunubar_events.pdf", xmax.at(ax_it->first));
526 
527  PlotMultiverse(*mv_ppfx .at(ax_it->first).at("nuebar" ), hists.at(ax_it->first).at("nuebar" ).at("nominal"), "Selected #bar{#nu_{e}} PPFX" , plot_dump + "/" + ax_it->first + "_nuebar_events_ppfx.pdf" );
528  PlotMultiverse(*mv_genie.at(ax_it->first).at("nuebar" ), hists.at(ax_it->first).at("nuebar" ).at("nominal"), "Selected #bar{#nu_{e}} GENIE", plot_dump + "/" + ax_it->first + "_nuebar_events_genie.pdf");
529  PlotMultiverse(*mv_ppfx .at(ax_it->first).at("nue" ), hists.at(ax_it->first).at("nue" ).at("nominal"), "Selected #nu_{e} PPFX" , plot_dump + "/" + ax_it->first + "_nue_events_ppfx.pdf" );
530  PlotMultiverse(*mv_genie.at(ax_it->first).at("nue" ), hists.at(ax_it->first).at("nue" ).at("nominal"), "Selected #nu_{e} GENIE", plot_dump + "/" + ax_it->first + "_nue_events_genie.pdf");
531  PlotMultiverse(*mv_ppfx .at(ax_it->first).at("nunubar"), hists.at(ax_it->first).at("nunubar").at("nominal"), "Selected #bar{#nu_{e}} + #nu_{e} PPFX" , plot_dump + "/" + ax_it->first + "_nunubar_events_ppfx.pdf" );
532  PlotMultiverse(*mv_genie.at(ax_it->first).at("nunubar"), hists.at(ax_it->first).at("nunubar").at("nominal"), "Selected #bar{#nu_{e}} + #nu_{e} GENIE", plot_dump + "/" + ax_it->first + "_nunubar_events_genie.pdf");
533 
534  PlotMultiverse(*mv_ppfx_nuebar_fractions .at(ax_it->first), nuebar_fractions.at(ax_it->first).at("nominal"), "#bar{#nu_{e}} Fraction PPFX" , plot_dump + "/" + ax_it->first + "_nuebar_fraction_ppfx.pdf",
535  nuebar_fractions.at(ax_it->first).at("ppfx"));
536 
537  PlotMultiverse(*mv_genie_nuebar_fractions.at(ax_it->first), nuebar_fractions.at(ax_it->first).at("nominal"), "#bar{#nu_{e}} Fraction GENIE", plot_dump + "/" + ax_it->first + "_nuebar_fraction_genie.pdf",
538  nuebar_fractions.at(ax_it->first).at("genie"));
539 
540  PlotMultiverse(*mv_ppfx_nue_nuebar_ratios.at(ax_it->first), nue_nuebar_ratios.at(ax_it->first).at("nominal"), "#nu_{e} / #bar{#nu}_{e} PPFX", plot_dump + "/" + ax_it->first + "_nue_nuebar_ratio_ppfx.pdf",
541  nue_nuebar_ratios.at(ax_it->first).at("ppfx"));
542 
543  PlotMultiverse(*mv_genie_nue_nuebar_ratios.at(ax_it->first), nue_nuebar_ratios.at(ax_it->first).at("nominal"), "#nu_{e} / #bar{#nu}_{e} GENIE", plot_dump + "/" + ax_it->first + "_nue_nuebar_ratio_genie.pdf",
544  nue_nuebar_ratios.at(ax_it->first).at("genie"));
545 
546 
547  }
548  gStyle->SetPaintTextFormat(".2f");
549  Plot2D(mc_frac_uncert.at("electron_kin"), nuebar_fractions_uncert_frac.at("electron_kin"), all_systs, plot_dump + "/electron_kin_fractional_uncertainty");
550 
551  TCanvas * c = new TCanvas();
552  nsel_uncert.at("q2")->GetYaxis()->SetTitle("Fractional Uncertainty");
553  nsel_uncert.at("q2")->GetXaxis()->SetRangeUser(0, xmax.at("q2"));
554  nsel_uncert.at("q2")->Draw("hist");
555  c->Print(TString::Format("%s/%s_nsel_uncert.pdf", plot_dump.c_str(), "q2").Data());
556 
557  nsel_uncert.at("enu")->GetYaxis()->SetTitle("Fractional Uncertainty");
558  nsel_uncert.at("enu")->GetXaxis()->SetRangeUser(0, xmax.at("enu"));
559  nsel_uncert.at("enu")->Draw("hist");
560  c->Print(TString::Format("%s/%s_nsel_uncert.pdf", plot_dump.c_str(), "enu").Data());
561 
562 
563  c->SetRightMargin(0.2);
564  auto nsel_uncert_electron_kin_2d = ToTH2Helper(std::unique_ptr<TH1>((TH1*) nsel_uncert.at("electron_kin")->Clone(UniqueName().c_str())),
567  nsel_uncert_electron_kin_2d->GetZaxis()->SetTitle("Fractional Uncertainty");
568 
569  nsel_uncert_electron_kin_2d->GetXaxis()->SetRangeUser(0.85, 1);
570  nsel_uncert_electron_kin_2d->GetYaxis()->SetRangeUser(0, 6);
571  nsel_uncert_electron_kin_2d->GetZaxis()->SetRangeUser(0, 1);
572  nsel_uncert_electron_kin_2d->SetMarkerColor(kWhite);
573  nsel_uncert_electron_kin_2d->GetXaxis()->SetTitle("True cos#theta_{e}");
574  nsel_uncert_electron_kin_2d->GetYaxis()->SetTitle("True Electron Energy (GeV)");
575  MovePalette(nsel_uncert_electron_kin_2d);
576  nsel_uncert_electron_kin_2d->Draw("colz text");
577 
578  c->Print(TString::Format("%s/%s_nsel_uncert.pdf", plot_dump.c_str(), "electron_kin").Data());
579 
580  auto nnunubar_2d = ToTH2Helper(std::unique_ptr<TH1>((TH1*) hists.at("electron_kin").at("nunubar").at("nominal")->Clone(UniqueName().c_str())),
583 
584  nnunubar_2d->SetTitle("Selected N_{#bar{#nu}_{e}} + N_{#nu_{e}}");
585  nnunubar_2d->GetXaxis()->SetRangeUser(0.85, 1);
586  nnunubar_2d->GetYaxis()->SetRangeUser(0, 6);
587  nnunubar_2d->SetMarkerColor(kWhite);
588  nnunubar_2d->GetXaxis()->SetTitle("True cos#theta_{e}");
589  nnunubar_2d->GetYaxis()->SetTitle("True Electron Energy (GeV)");
590  MovePalette(nnunubar_2d);
591  nnunubar_2d->Draw("colz text");
592 
593  c->Print(TString::Format("%s/%s_nsel_nunubar.pdf", plot_dump.c_str(), "electron_kin").Data());
594 
595  auto nom_nuebar_frac_2d = ToTH2Helper(std::unique_ptr<TH1>((TH1*) nuebar_fractions.at("electron_kin").at("nominal")->Clone(UniqueName().c_str())),
598 
599  nom_nuebar_frac_2d->SetTitle("s_{#bar{#nu}_{e}}");
600  nom_nuebar_frac_2d->GetXaxis()->SetRangeUser(0.85, 1);
601  nom_nuebar_frac_2d->GetYaxis()->SetRangeUser(0, 6);
602  nom_nuebar_frac_2d->GetZaxis()->SetRangeUser(0, 1);
603  nom_nuebar_frac_2d->SetMarkerColor(kWhite);
604  nom_nuebar_frac_2d->GetXaxis()->SetTitle("True cos#theta_{e}");
605  nom_nuebar_frac_2d->GetYaxis()->SetTitle("True Electron Energy (GeV)");
606  MovePalette(nom_nuebar_frac_2d);
607  nom_nuebar_frac_2d->Draw("colz text");
608 
609  c->Print(TString::Format("%s/%s_nuebar_fraction_2d.pdf", plot_dump.c_str(), "electron_kin").Data());
610 
611 
612 
613  // calculate required MC POT that would make MC stat. uncert
614  // negligible -- ~1%
615 
616  std::map<std::string, TH1*> req_pot_1d;
617  std::map<std::string, TH1*> req_pot_2d;
618  std::map<std::string, TH1*> req_pot_wrt_prod5;
619  auto target_uncert = 0.01;
620  for(auto loader_it = loaders.begin(); loader_it != loaders.end(); loader_it++) {
621  req_pot_1d[loader_it->first] = (TH1*) nuebar_fractions.at("electron_kin").at(loader_it->first)->Clone(UniqueName().c_str());
622  req_pot_1d.at(loader_it->first)->Reset();
623 
624  for(auto i = 1; i <= req_pot_1d.at(loader_it->first)->GetNbinsX(); i++) {
625  double s = nuebar_fractions.at("electron_kin").at(loader_it->first)->GetBinContent(i);
626  double n = hists.at("electron_kin").at("nunubar").at(loader_it->first)->GetBinContent(i);
627  double pot = pots.at(loader_it->first);
628 
629  if(s != 0. && s != 1)
630  req_pot_1d.at(loader_it->first)->SetBinContent(i, (1+s)/s * std::pow(target_uncert, -2) * pot / n);
631  }
632 
633 
634  req_pot_2d[loader_it->first] = ToTH2Helper(std::unique_ptr<TH1>((TH1*) req_pot_1d.at(loader_it->first)->Clone(UniqueName().c_str())),
637  req_pot_wrt_prod5[loader_it->first] = (TH1*) req_pot_2d.at(loader_it->first)->Clone(UniqueName().c_str());
638  req_pot_wrt_prod5.at(loader_it->first)->Scale(1./pots.at(loader_it->first));
639 
640 
641  gStyle->SetPaintTextFormat("1.1g");
642  gStyle->SetPalette(kDeepSea);
643  c->SetLogz();
644  req_pot_2d.at(loader_it->first)->SetTitle(TString::Format("Required MC POT (%.1f%% Stat. Uncert)", target_uncert * 100));
645  req_pot_2d.at(loader_it->first)->GetXaxis()->SetRangeUser(0.85, 1);
646  req_pot_2d.at(loader_it->first)->GetYaxis()->SetRangeUser(0, 6);
647  req_pot_2d.at(loader_it->first)->SetMarkerColor(kWhite);
648  req_pot_2d.at(loader_it->first)->GetXaxis()->SetTitle("True cos#theta_{e}");
649  req_pot_2d.at(loader_it->first)->GetYaxis()->SetTitle("True Electron Energy (GeV)");
650  req_pot_2d.at(loader_it->first)->GetZaxis()->SetRangeUser(1e15, 1e25);
651  MovePalette(req_pot_2d.at(loader_it->first));
652  req_pot_2d.at(loader_it->first)->Draw("colz text");
653 
654 
655  c->Print(TString::Format("%s/req_pot_%s.pdf", plot_dump.c_str(), loader_it->first.c_str()).Data());
656 
657  req_pot_wrt_prod5.at(loader_it->first)->SetTitle("Required MC POT / Prod5 POT");
658  req_pot_wrt_prod5.at(loader_it->first)->GetXaxis()->SetRangeUser(0.85, 1);
659  req_pot_wrt_prod5.at(loader_it->first)->GetYaxis()->SetRangeUser(0, 6);
660  req_pot_wrt_prod5.at(loader_it->first)->SetMaximum(500);
661  req_pot_wrt_prod5.at(loader_it->first)->SetMarkerColor(kWhite);
662  req_pot_wrt_prod5.at(loader_it->first)->GetXaxis()->SetTitle("True cos#theta_{e}");
663  req_pot_wrt_prod5.at(loader_it->first)->GetYaxis()->SetTitle("True Electron Energy (GeV)");
664  MovePalette(req_pot_wrt_prod5.at(loader_it->first));
665  req_pot_wrt_prod5.at(loader_it->first)->Draw("colz text");
666 
667 
668  gStyle->SetPaintTextFormat(".1f");
669  c->SetLogz(false);
670  c->Print(TString::Format("%s/req_pot_wrt_prod5%s.pdf", plot_dump.c_str(), loader_it->first.c_str()).Data());
671  }
672 
673  delete c;
674  c = new TCanvas();
675  TLegend * leg = 0;
676  auto kNuebarColor = kGreen;
677  auto kNueColor = kMagenta;
678  // 1D stacked nnunubar plots
679  for(auto loader_it = loaders.begin(); loader_it != loaders.end(); loader_it++) {
680  for(std::string axis : {"enu", "q2"}) {
681  if(leg)
682  delete leg;
683  leg = new TLegend(0.6, 0.89, 0.8, 0.7);
684 
685  THStack * stack = new THStack(UniqueName().c_str(), "");
686  hists.at(axis).at("nue") .at(loader_it->first)->SetFillColor(kNueColor);
687  hists.at(axis).at("nuebar").at(loader_it->first)->SetFillColor(kNuebarColor);
688  hists.at(axis).at("nue") .at(loader_it->first)->SetLineColor(kBlack);
689  hists.at(axis).at("nuebar").at(loader_it->first)->SetLineColor(kBlack);
690  hists.at(axis).at("nue") .at(loader_it->first)->SetLineStyle(1);
691  hists.at(axis).at("nuebar").at(loader_it->first)->SetLineStyle(1);
692 
693  stack->Add(hists.at(axis).at("nue" ).at(loader_it->first));
694  stack->Add(hists.at(axis).at("nuebar").at(loader_it->first));
695  leg->AddEntry(hists.at(axis).at("nue" ).at(loader_it->first), "#nu_{e}" , "f");
696  leg->AddEntry(hists.at(axis).at("nuebar").at(loader_it->first), "#bar{#nu}_{e}", "f");
697 
698 
699  stack->Draw("hist");
700  stack->GetXaxis()->SetRangeUser(0, xmax.at(axis));
701  stack->GetXaxis()->SetTitle(hists.at(axis).at("nue").at(loader_it->first)->GetXaxis()->GetTitle());
702  stack->GetYaxis()->SetTitle("Events");
703  leg->Draw("same");
704  c->Print(TString::Format("%s/%s_nunubar_stacked_%s.pdf",
705  plot_dump.c_str(),
706  axis.c_str(),
707  loader_it->first.c_str()).Data());
708 
709  }
710  }
711 
712 
713  }
714 }
715 
716 void abs(TH1* hist)
717 {
718  for(auto x = 1; x <= hist->GetNbinsX(); x++) {
719  for(auto y = 1; y <= hist->GetNbinsY(); y++) {
720  for(auto z = 1; z <= hist->GetNbinsZ(); z++) {
721  auto val = hist->GetBinContent(x, y, z);
722  hist->SetBinContent(x, y, z, std::abs(val));
723  }
724  }
725  }
726 }
727 
728 void PlotMultiverse(Multiverse & mv, TH1 * nominal, std::string title, std::string name, TH1 * error)
729 {
730  auto hists = mv.GetUniversesHist(kCyan);
731 
732  double maxy = -1000;
733  for(auto ihist = 0u; ihist < hists.size(); ihist++)
734  maxy = std::max(maxy, hists[ihist]->GetMaximum());
735 
736  TCanvas * c = new TCanvas();
737  nominal->GetYaxis()->SetRangeUser(0, 1.2 * maxy);
738  nominal->SetLineColor(kBlue);
739  nominal->SetTitle(title.c_str());
740 
741  nominal->Draw("hist");
742 
743  for(auto ihist = 0u; ihist < hists.size(); ihist++)
744  hists[ihist]->Draw("hist same");
745 
746 
747  if(error) {
748  for(auto x = 1; x <= nominal->GetNbinsX(); x++) {
749  nominal->SetBinError(x, std::abs(error->GetBinContent(x) - nominal->GetBinContent(x)));
750  }
751  nominal->Draw("hist same e1");
752  }
753  else
754  nominal->Draw("hist same");
755 
756 
757  c->Print(name.c_str());
758 }
759 
760 void Plot2D(TH1* total, std::map<std::string, TH1*> systs, std::vector<std::string> to_plot, std::string basename)
761 {
762 
763  TCanvas * c = new TCanvas();
764  c->SetRightMargin(0.2);
765  std::map<std::string, TH1*> hists2d;
766  hists2d["total"] = ToTH2Helper(std::unique_ptr<TH1>((TH1*) total->Clone(UniqueName().c_str())),
769  for(auto syst : to_plot) {
770  hists2d[syst] = ToTH2Helper(std::unique_ptr<TH1>((TH1*) systs.at(syst)->Clone(UniqueName().c_str())),
773  }
774 
775  for(auto hist_it = hists2d.begin(); hist_it != hists2d.end(); hist_it++) {
776  hist_it->second->SetTitle("");
777  hist_it->second->GetXaxis()->SetTitle("True cos#theta_{e}");
778  hist_it->second->GetYaxis()->SetTitle("True Electron Energy (GeV)");
779  hist_it->second->GetXaxis()->SetRangeUser(0.85, 1);
780  hist_it->second->GetYaxis()->SetRangeUser(0, 6);
781  hist_it->second->SetMaximum(1);
782  hist_it->second->SetMarkerColor(kWhite);
783  MovePalette(hist_it->second);
784  hist_it->second->Draw("colz text");
785  c->Print((basename + "_" + hist_it->first + ".pdf").c_str());
786  }
787 
788 }
789 
790 void Plot1D(TH1* h1, std::string h1_label, std::map<std::string, TH1*> systs, std::vector<std::string> to_plot, std::string name, double xmax) {
791  TCanvas * c = new TCanvas();
792  TLegend * leg = new TLegend();
793 
794  // find max y
795  double maxy = -999;
796  for(auto syst : to_plot)
797  maxy = std::max(maxy, systs.at(syst)->GetMaximum());
798 
799  maxy = std::max(maxy, h1->GetMaximum());
800 
801 
802  if(xmax > 0)
803  h1->GetXaxis()->SetRangeUser(0, xmax);
804 
805  h1->GetYaxis()->SetRangeUser(0, maxy * 1.2);
806  h1->SetLineColor(kBlack);
807  h1->Draw("hist");
808  leg->AddEntry(h1, h1_label.c_str(), "l");
809 
810  for(auto syst : to_plot ) {
811  if(syst == "nominal") continue;
812  systs.at(syst)->SetLineColor(syst_styles.at(syst).color);
813  systs.at(syst)->SetLineStyle(syst_styles.at(syst).line_style);
814  leg->AddEntry(systs.at(syst), syst.c_str(), "l");
815  systs.at(syst)->Draw("hist same");
816  }
817  h1->Draw("hist same");
818  leg->Draw();
819  c->Print(name.c_str());
820 }
821 TH1 * MostConservative(TH1 * up, TH1 * down, TH1 * nominal)
822 {
823  auto ret = (TH1*) up->Clone(UniqueName().c_str());
824  auto NX = up->GetNbinsX();
825  auto NY = up->GetNbinsY();
826  auto NZ = up->GetNbinsZ();
827  for(auto x = 1; x <= NX; x++) {
828  for(auto y = 1; y <= NY; y++) {
829  for(auto z = 1; z <= NZ; z++) {
830  auto nom = nominal->GetBinContent(x, y, z);
831  auto shift_up = std::abs(nom - up->GetBinContent(x, y, z));
832  auto shift_down = std::abs(nom - down->GetBinContent(x, y, z));
833 
834  ret->SetBinContent(x, y, z, std::max(shift_up, shift_down) + nom);
835  }
836  }
837  }
838  return ret;
839 }
840 
841 void MovePalette(TH1 * hist)
842 {
843  hist->Draw("colz");
844  gPad->Update();
845  auto palette = (TPaletteAxis*) hist->GetListOfFunctions()->FindObject("palette");
846  palette->SetX1NDC(0.81);
847  palette->SetX2NDC(0.85);
848 }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
tree Draw("slc.nhit")
ofstream output
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
TH2 * ToTH2(const TH1 *h, const HistAxis &xaxis, const HistAxis &yaxis)
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const BeamSyst kBeamTarget((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"7mmTargetZ","+/- 7mm Target z Position")
Target z position shift +/-7mm.
Definition: BeamSysts.h:113
Spectrum * GetPlusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:443
const ana::Binning eelecbins
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const ana::Cut kNuebarCC
int pots
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
double maxy
std::vector< std::string > beam_systs
Float_t den
Definition: plot.C:36
const BeamSyst kBeamHornCurrent((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"2kA","2kAHornCurrent","+/- 2kA Horn Current")
Horn Current +/-2kA.
Definition: BeamSysts.h:97
const ana::Binning costhetabins
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
void Plot2D(TH1 *total, std::map< std::string, TH1 * > systs, std::vector< std::string > to_plot, std::string basename)
const std::string PROD5_MC_RHC_LIGHT_UP
void abs(TH1 *hist)
const std::string PROD5_MC_RHC_CALIB_DOWN
TString hists[nhists]
Definition: bdt_com.C:3
const ana::HistAxis kTrueElectronEVsCosStandardAxis
const std::string PROD5_MC_RHC_LIGHT_DOWN
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
Spectrum * GetMinusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:450
float abs(float number)
Definition: d0nt_math.hpp:39
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
const BeamSyst kBeamGeomWater((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmHornWater","+/- 1mm water on Horn 1")
Water layer on horn 1: +/- 1mm.
Definition: BeamSysts.h:122
const XML_Char * s
Definition: expat.h:262
std::vector< const ISyst * > getAllXsecSysts_2020()
void Cut(double x)
Definition: plot_outliers.C:1
const BeamSyst kBeamH1PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn1X","+/- 3mm Horn 1 X Position")
Horn 1 and 2 position +-3mm in X and Y separately.
Definition: BeamSysts.h:107
void MovePalette(TH1 *hist)
void load_libs_muonid()
void Plot1D(TH1 *h1, std::string h1_label, std::map< std::string, TH1 * > systs, std::vector< std::string > to_plot, std::string name, double xmax=-1)
const std::string PROD5_MC_RHC_NOMINAL
void wrong_sign_uncertainty(bool make_plots=true, std::string input_file_name="/nova/ana/users/ddoyle/NuebarCCInc/WrongSign/wrong_sign_uncertainty.root", std::string plot_dump="/nova/ana/users/ddoyle/NuebarCCInc/WrongSign/plots")
const std::string PROD5_MC_RHC_CALIB_UP
const std::string PROD5_MC_RHC_CALIB_SHAPE
#define pot
std::string xaxis
std::vector< float > Spectrum
Definition: Constants.h:527
z
Definition: test.py:28
std::map< std::string, nuebarccinc::Style > syst_styles
const ana::Cut kPreselectionLoose
const double kAna2020RHCPOT
Definition: Exposures.h:235
const std::string PROD5_MC_RHC_CHERENKOV
const SystShifts kNoShift
Definition: SystShifts.h:115
const BeamSyst kBeamPosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmBeamShiftY","Beam Position Y")
Definition: BeamSysts.h:104
const ana::HistAxis kTrueQ2StandardAxis
const Cut kSelection
Definition: SINCpi0_Cuts.h:328
TH1F * h1
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const BeamSyst kBeamH1PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn1Y","+/- 3mm Horn 1 Y Position")
Definition: BeamSysts.h:108
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const BeamSyst kBeamSpotSize((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"0.2mmBeamSpotSize","0p2mmBeamSpotSize"," 1.3 +/- 0.2 mm Spot Size")
Beam Spot Size 1.3 +/- 0.2 mm both XY.
Definition: BeamSysts.h:100
int num
Definition: f2_nu.C:119
const BeamSyst kBeamH2PosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn2X","+/- 3mm Horn 2 X Position")
Definition: BeamSysts.h:109
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< SystShifts > GetSystShifts()
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
const std::vector< TH1 * > GetUniversesHist()
Definition: Multiverse.cxx:120
exit(0)
std::string yaxis
TH2 * ToTH2Helper(std::unique_ptr< TH1 > h1, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Helper for ana::ToTH2.
Definition: Utilities.cxx:276
const BeamSyst kBeamPosX((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"1mmBeamShiftX","Beam Position X")
Beam position on target +-1 mm, X/Y separately.
Definition: BeamSysts.h:103
const ana::Cut kMuonIDProd4Cut
Definition: datamc.C:28
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const ana::Cut kNueCC
void PlotMultiverse(Multiverse &mv, TH1 *nominal, std::string title, std::string name, TH1 *error=0)
const BeamSyst kBeamH2PosY((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"3mmHorn2Y","+/- 3mm Horn 2 Y Position")
Definition: BeamSysts.h:110
const BeamSyst kBeamMagField((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"Magnetic Field in Decay Pipe","MagnFieldDecayPipe","Magnetic Field in Decay Pipe")
Constant magnetic field in decay pipe.
Definition: BeamSysts.h:116
const ana::HistAxis kTrueNeutrinoEnergyStandardAxis
TH1 * MostConservative(TH1 *up, TH1 *down, TH1 *nominal)
void make_plots(TFile *f, TH2 *hVsRun, TGraph *g, int run, std::string suffix, double xpos, TH2 *&cut)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
void Divide(Multiverse &, bool=false)
Definition: Multiverse.cxx:139
std::vector< std::string > all_systs
std::map< string, const ana::BeamSyst * > beamSystematics
const Var kXSecCVWgt2020
Definition: XsecTunes.h:106