NusSystsMaker.cxx
Go to the documentation of this file.
2 
7 
9 
10 #include "TCanvas.h"
11 #include "TLegend.h"
12 #include "TLegendEntry.h"
13 #include "TObjString.h"
14 #include "TKey.h"
15 #include "TVectorD.h"
16 
17 namespace ana {
18 
19  // -------------------------------------------------------------------------------------
21  const HistAxis* axis, const HistAxis* numu_axis,
22  const Cut* fd_cut, const Cut* nd_cut, const Cut* numu_cut,
23  const SystShifts* shift_data, const SystShifts* shift_mc,
24  const Var* weight,
25  Loaders* loaders, Loaders* loaders_nd) {
26 
27  IPrediction* pred = nullptr;
28 
29  if (pred_type == NusSystsPredType::kExtrap) {
30  SterileGenerator genExtrap(*axis, *numu_axis,
31  *fd_cut, *nd_cut, *numu_cut,
32  *shift_data, *weight);
33  if (loaders_nd)
34  pred = genExtrap.Generate(*loaders_nd, *loaders, *shift_mc).release();
35  else
36  pred = genExtrap.Generate(*loaders, *shift_mc).release();
37  }
38 
39  else if (pred_type == NusSystsPredType::kFD) {
40  FDPredictionGenerator genFD(*axis, *fd_cut, *shift_data, *weight);
41  pred = genFD.Generate(*loaders, *shift_mc).release();
42  }
43 
44  else if (pred_type == NusSystsPredType::kND) {
45  NDPredictionGenerator genND(*axis, *nd_cut, *shift_data, *weight);
46  pred = genND.Generate(*loaders, *shift_mc).release();
47  }
48 
49  return pred;
50 
51  } // function GetPrediction
52 
53  // --------------------------------------------------------------------------
54  TH1* GetSpectrum(IPrediction* pred, bool signal, double pot) {
56  calc->SetDm(4, 0);
57  return pred->PredictComponent(calc, Flavors::kAll,
58  signal ? Current::kNC : Current::kCC,
59  Sign::kBoth).ToTH1(pot);
60  }
61 
62  // --------------------------------------------------------------------------
64 
65  TString mystring(instring.c_str());
66 
67  std::vector<TString> in = {"0","1","2","3","4","5",
68  "6","7","8","9","_","#"," ",".","-","(",")"};
69  std::vector<TString> out = {"zero","one","two","three","four","five",
70  "six","seven","eight","nine","","","","","XX","",""};
71  for(size_t i = 0; i < in.size(); ++i)
72  mystring.ReplaceAll(in[i],out[i]);
73 
74  return std::string(mystring.Data());
75  }
76 
77  // -------------------------------------------------------------------------------------
78  /// NusSystShiftSigma constructor
80  Name(name),
81  Sigma(sigma) {
82  }
83 
84  // -------------------------------------------------------------------------------------
85  /// Print the sigmas for a given shift
86  void NusSystShiftSigma::PrintSigma(std::ostream& os) {
87  os << " " << std::showpos << Sigma << ": "
88  << Shift.first << " (signal), "
89  << Shift.second << " (background)"
90  << std::endl;
91 
92  }
93 
94  // -------------------------------------------------------------------------------------
95  void NusSystShiftSigma::SaveTo(TDirectory* dir) const {
96 
97  TDirectory* tmp = gDirectory;
98 
99  dir->cd();
100  TObjString("NusSystShiftSigma").Write("type");
101 
102  TObjString(Name.c_str()).Write("name");
103  TVectorD sigma(1);
104  sigma[0] = Sigma;
105  sigma.Write("sigma");
106 
107  Prediction->SaveTo(dir->mkdir("pred"));
108  TH1* shiftsig = ShiftedSpectrum.first;
109  TH1* shiftbg = ShiftedSpectrum.second;
110  if (shiftsig)
111  shiftsig->Write("shiftsig");
112  if (shiftbg)
113  shiftbg->Write("shiftbg");
114  TH1* shiftratsig = ShiftedRatio.first;
115  TH1* shiftratbg = ShiftedRatio.second;
116  if (shiftratsig)
117  shiftratsig->Write("shiftratsig");
118  if (shiftratbg)
119  shiftratbg->Write("shiftratbg");
120 
121  TVectorD v(2);
122  v[0] = Shift.first;
123  v[1] = Shift.second;
124  v.Write("shift");
125 
126  tmp->cd();
127 
128  }
129 
130  // -------------------------------------------------------------------------------------
132 
133  DontAddDirectory guard;
134 
135  TObjString* tag = (TObjString*)dir->Get("type");
136  assert(tag);
137  assert(tag->GetString() == "NusSystShiftSigma");
138 
139  TObjString* name = (TObjString*)dir->Get("name");;
140  const TVectorD sigma = *(TVectorD*)dir->Get("sigma");
141  NusSystShiftSigma* nusSystShiftSigma = new NusSystShiftSigma(name->GetString().Data(),
142  sigma[0]);
143 
144  nusSystShiftSigma->Prediction = LoadFromFile<IPrediction>(static_cast<TFile*>(dir), "pred").release();
145  TH1* shiftsig = (TH1*)dir->Get("shiftsig");
146  TH1* shiftbg = (TH1*)dir->Get("shiftbg");
147  if (shiftsig and shiftbg)
148  nusSystShiftSigma->ShiftedSpectrum = std::make_pair(shiftsig, shiftbg);
149  TH1* shiftratsig = (TH1*)dir->Get("shiftratsig");
150  TH1* shiftratbg = (TH1*)dir->Get("shiftratbg");
151  if (shiftratsig and shiftratbg)
152  nusSystShiftSigma->ShiftedRatio = std::make_pair(shiftratsig, shiftratbg);
153 
154  const TVectorD v = *(TVectorD*)dir->Get("shift");
155  nusSystShiftSigma->Shift = std::make_pair(v[0], v[1]);
156 
157  return nusSystShiftSigma;
158 
159  }
160 
161  // -------------------------------------------------------------------------------------
163  fName(name),
164  fLabel(label)
165  {}
166 
167  // -------------------------------------------------------------------------------------
169  {}
170 
171  // -------------------------------------------------------------------------------------
172  std::vector<int> NusSystShift::GetSigmas() {
173  std::vector<int> sigmas;
174  for (const auto& sigma : fSigmas)
175  sigmas.push_back(sigma.first);
176  return sigmas;
177  }
178 
179  // -------------------------------------------------------------------------------------
180  void NusSystShift::PrintShift(std::ostream& os) {
181  os << std::endl << " " << fLabel << " shifts:" << std::endl;
182  for (const auto& sigma : fSigmas)
183  sigma.second->PrintSigma(os);
184  return;
185  }
186 
187  // -------------------------------------------------------------------------------------
189  // Require -1, +1 shifts
190  bool good = true;
191  if (!fSigmas.count(-1) or !fSigmas.count(+1))
192  good = false;
193  return good;
194  }
195 
196  // -------------------------------------------------------------------------------------
197  TCanvas* NusSystShift::DrawSig(TH1* nominal) {
198 
199  // Drawing tools ------------
200  TCanvas* canv = new TCanvas("canv", "", 800, 800);
201  TPad* pad_spec_sig_total = new TPad("pad_spec_sig_total", "", 0., 0.375, 1., 1.);
202  TPad* pad_ratio_sig_total = new TPad("pad_ratio_sig_total", "", 0., 0., 1., 0.375);
203  TLegend* leg = new TLegend(0.5, 0.5, 0.88, 0.88);
204  std::map<int, TLegendEntry*> leg_entries;
205 
206  // Draw nominal ------------
207  TH1* h_nom = (TH1*)nominal->Clone("h_nom");
208  for (size_t bin = 1; bin <= size_t(h_nom->GetNbinsX()); ++bin)
209  h_nom->SetBinContent(bin, h_nom->GetBinContent(bin)/h_nom->GetBinWidth(bin));
210  h_nom->GetYaxis()->SetTitle("Events / GeV");
211  h_nom->SetLineStyle(kSolid);
212  h_nom->SetLineColor(kBlack);
213  h_nom->SetLineWidth(2);
214  pad_spec_sig_total->cd();
215  h_nom->Draw("hist");
216 
217  // ratio
218  TH1* nom_ratio = (TH1*)nominal->Clone("nom_ratio");
219  nom_ratio->GetYaxis()->SetTitle("Shifted / Nominal");
220  for (Int_t bin = 0; bin <= nom_ratio->GetNbinsX(); ++bin)
221  nom_ratio->SetBinContent(bin, 1.);
222  pad_ratio_sig_total->cd();
223  nom_ratio->Draw("hist");
224  leg_entries[0] = new TLegendEntry(nominal, "Nominal", "l");
225  double ymax_ratio = 0;
226 
227  // Draw shifts ------------
228  int it = 0;
229  for (std::map<int, NusSystShiftSigma*>::const_iterator sigma = fSigmas.begin(); sigma != fSigmas.end(); ++sigma, ++it) {
230  if (sigma->second->Sigma == 0)
231  continue;
232  // Total ------------
233  TH1* sigma_sig = (TH1*)sigma->second->ShiftedSpectrum.first->Clone();
234  for (size_t bin = 1; bin <= size_t(sigma_sig->GetNbinsX()); ++bin)
235  sigma_sig->SetBinContent(bin, sigma_sig->GetBinContent(bin)/sigma_sig->GetBinWidth(bin));
236  sigma_sig->GetYaxis()->SetTitle("Events / GeV");
237  sigma_sig->SetLineStyle(kDashed);
238  sigma_sig->SetLineColor(abs(sigma->second->Sigma)+1);
239  sigma_sig->SetLineWidth(2);
240  pad_spec_sig_total->cd();
241  sigma_sig->DrawClone("hist same");
242  TH1* sigma_ratio_sig = (TH1*)sigma_sig->Clone("sigma_ratio_sig");
243  sigma_ratio_sig->Divide(h_nom);
244  sigma_ratio_sig->SetBinContent(0, 1);
245  for (size_t i = 1; i <= (size_t)sigma_ratio_sig->GetNbinsX(); ++i) {
246  double val = sigma_ratio_sig->GetBinContent(i);
247  if (val == 0) {
248  sigma_ratio_sig->SetBinContent(i, 1);
249  val = 1;
250  }
251  if (fabs(val-1) > ymax_ratio)
252  ymax_ratio = fabs(val-1);
253  }
254  pad_ratio_sig_total->cd();
255  sigma_ratio_sig->DrawClone("hist same");
256  leg_entries[abs(sigma->second->Sigma)]
257  = new TLegendEntry(sigma_sig->Clone(), Form("+/-%d#sigma", sigma->second->Sigma), "l");
258  }
259 
260  // Draw legend ------------
261  for (const auto& leg_entry : leg_entries)
262  leg->AddEntry(leg_entry.second->GetObject(),
263  leg_entry.second->GetLabel(),
264  leg_entry.second->GetOption());
265  pad_spec_sig_total->cd();
266  leg->Draw();
267 
268  nom_ratio->SetMinimum(1-(1.1*ymax_ratio));
269  nom_ratio->SetMaximum(1+(1.1*ymax_ratio));
270 
271  // Draw all together ------------
272  canv->cd();
273  canv->Clear();
274  pad_spec_sig_total->Draw();
275  pad_ratio_sig_total->Draw();
276 
277  return canv;
278 
279  } // function NusSystShift::DrawSig
280 
281  // -------------------------------------------------------------------------------------
282  TCanvas* NusSystShift::DrawBkg(TH1* nominal) {
283 
284  // Drawing tools ------------
285  TCanvas* canv = new TCanvas("canv", "", 800, 800);
286  TPad* pad_spec_bkg_total = new TPad("pad_spec_bkg_total", "", 0., 0.375, 1., 1.);
287  TPad* pad_ratio_bkg_total = new TPad("pad_ratio_bkg_total", "", 0., 0., 1., 0.375);
288  TLegend* leg = new TLegend(0.5, 0.5, 0.88, 0.88);
289  std::map<int, TLegendEntry*> leg_entries;
290 
291  // Draw nominal ------------
292  TH1* h_nom = (TH1*)nominal->Clone("h_nom");
293  for (size_t bin = 1; bin <= size_t(h_nom->GetNbinsX()); ++bin)
294  h_nom->SetBinContent(bin, h_nom->GetBinContent(bin)/h_nom->GetBinWidth(bin));
295  h_nom->GetYaxis()->SetTitle("Events / GeV");
296  nominal->SetLineStyle(kSolid);
297  nominal->SetLineColor(kBlack);
298  nominal->SetLineWidth(2);
299  pad_spec_bkg_total->cd();
300  nominal->Draw("hist");
301 
302  // ratio
303  TH1* nom_ratio = (TH1*)nominal->Clone("nom_ratio");
304  nom_ratio->GetYaxis()->SetTitle("Shifted / Nominal");
305  for (Int_t bin = 0; bin <= nom_ratio->GetNbinsX(); ++bin)
306  nom_ratio->SetBinContent(bin, 1.);
307  pad_ratio_bkg_total->cd();
308  nom_ratio->Draw("hist");
309  leg_entries[0] = new TLegendEntry(nominal, "Nominal", "l");
310  double ymax_ratio = 0;
311 
312  // Draw shifts ------------
313  int it = 0;
314  for (std::map<int, NusSystShiftSigma*>::const_iterator sigma = fSigmas.begin();
315  sigma != fSigmas.end(); ++sigma, ++it) {
316  if (sigma->second->Sigma == 0)
317  continue;
318  // Total ------------
319  TH1* sigma_bkg = (TH1*)sigma->second->ShiftedSpectrum.second->Clone();
320  for (size_t bin = 1; bin <= size_t(sigma_bkg->GetNbinsX()); ++bin)
321  sigma_bkg->SetBinContent(bin, sigma_bkg->GetBinContent(bin)/sigma_bkg->GetBinWidth(bin));
322  sigma_bkg->GetYaxis()->SetTitle("Events / GeV");
323  sigma_bkg->SetLineStyle(kDashed);
324  sigma_bkg->SetLineColor(abs(sigma->second->Sigma)+1);
325  sigma_bkg->SetLineWidth(2);
326  pad_spec_bkg_total->cd();
327  sigma_bkg->DrawClone("hist same");
328  TH1* sigma_ratio_bkg = (TH1*)sigma_bkg->Clone("sigma_ratio_bkg");
329  sigma_ratio_bkg->Divide(h_nom);
330  sigma_ratio_bkg->SetBinContent(0, 1);
331  for (size_t i = 1; i <= (size_t)sigma_ratio_bkg->GetNbinsX(); ++i) {
332  double val = sigma_ratio_bkg->GetBinContent(i);
333  if (val == 0) {
334  sigma_ratio_bkg->SetBinContent(i, 1);
335  val = 1;
336  }
337  if (fabs(val-1) > ymax_ratio)
338  ymax_ratio = fabs(val-1);
339  }
340  pad_ratio_bkg_total->cd();
341  sigma_ratio_bkg->DrawClone("hist same");
342  leg_entries[abs(sigma->second->Sigma)]
343  = new TLegendEntry(sigma_bkg->Clone(), Form("+/-%d#sigma",
344  sigma->second->Sigma), "l");
345  }
346 
347  // Draw legend ------------
348  for (const auto& leg_entry : leg_entries)
349  leg->AddEntry(leg_entry.second->GetObject(),
350  leg_entry.second->GetLabel(),
351  leg_entry.second->GetOption());
352  pad_spec_bkg_total->cd();
353  leg->Draw();
354 
355  nom_ratio->SetMinimum(1-(1.1*ymax_ratio));
356  nom_ratio->SetMaximum(1+(1.1*ymax_ratio));
357 
358  // Draw all together ------------
359  canv->cd();
360  canv->Clear();
361  pad_spec_bkg_total->Draw();
362  pad_ratio_bkg_total->Draw();
363 
364  return canv;
365 
366  } // function NusSystShift::DrawBkg
367 
368  // -------------------------------------------------------------------------------------
369  void NusSystShift::DrawShift(TDirectory* outFile, std::pair<TH1*, TH1*>& nominal) {
370 
371  // Out directories ------------
372  std::string dir_sig_shift = Form("Signal/Shifts/%s", fName.c_str());
373  std::string dir_bg_shift = Form("Background/Shifts/%s", fName.c_str());
374  std::string dir_sig_shift_sigmas = dir_sig_shift+"/Sigmas/";
375  std::string dir_bg_shift_sigmas = dir_bg_shift+"/Sigmas/";
376  outFile->mkdir(dir_sig_shift.c_str());
377  outFile->mkdir(dir_bg_shift.c_str());
378  outFile->mkdir(dir_sig_shift_sigmas.c_str());
379  outFile->mkdir(dir_bg_shift_sigmas.c_str());
380 
381  // Drawing tools ------------
382  TCanvas* canv = new TCanvas("canv", "", 800, 800);
383  TPad* pad_spec_sig = new TPad("pad_spec_sig", "", 0., 0.375, 1., 1.);
384  TPad* pad_ratio_sig = new TPad("pad_ratio_sig", "", 0., 0., 1., 0.375);
385  TPad* pad_spec_bg = new TPad("pad_spec_bg", "", 0., 0.375, 1., 1.);
386  TPad* pad_ratio_bg = new TPad("pad_ratio_bg", "", 0., 0., 1., 0.375);
387  TPad* pad_spec_sig_total = new TPad("pad_spec_sig_total", "", 0., 0.375, 1., 1.);
388  TPad* pad_ratio_sig_total = new TPad("pad_ratio_sig_total", "", 0., 0., 1., 0.375);
389  TPad* pad_spec_bg_total = new TPad("pad_spec_bg_total", "", 0., 0.375, 1., 1.);
390  TPad* pad_ratio_bg_total = new TPad("pad_ratio_bg_total", "", 0., 0., 1., 0.375);
391  TLegend* leg = new TLegend(0.5, 0.5, 0.88, 0.88);
392  std::map<int, TLegendEntry*> leg_entries;
393 
394  // Draw nominal ------------
395  // signal
396  TH1* nom_sig = (TH1*)nominal.first->Clone("nom_sig");
397  for (size_t bin = 1; bin <= size_t(nom_sig->GetNbinsX()); ++bin)
398  nom_sig->SetBinContent(bin, nom_sig->GetBinContent(bin)/nom_sig->GetBinWidth(bin));
399  nom_sig->GetYaxis()->SetTitle("Events / GeV");
400  nom_sig->SetLineStyle(kSolid);
401  nom_sig->SetLineColor(kBlack);
402  nom_sig->SetLineWidth(2);
403  pad_spec_sig_total->cd();
404  nom_sig->Draw("hist");
405  // background
406  TH1* nom_bg = (TH1*)nominal.second->Clone("nom_bg");
407  for (size_t bin = 1; bin <= size_t(nom_bg->GetNbinsX()); ++bin)
408  nom_bg->SetBinContent(bin, nom_bg->GetBinContent(bin)/nom_bg->GetBinWidth(bin));
409  nom_bg->GetYaxis()->SetTitle("Events / GeV");
410  nom_bg->SetLineStyle(kSolid);
411  nom_bg->SetLineColor(kBlack);
412  nom_bg->SetLineWidth(2);
413  pad_spec_bg_total->cd();
414  nom_bg->Draw("hist");
415  // ratio
416  TH1* nom_ratio = (TH1*)nom_sig->Clone("nom_ratio");
417  nom_ratio->GetYaxis()->SetTitle("Shifted / Nominal");
418  for (Int_t bin = 1; bin <= nom_ratio->GetNbinsX(); ++bin)
419  nom_ratio->SetBinContent(bin, 1.);
420  pad_ratio_sig_total->cd();
421  nom_ratio->Draw("hist");
422  pad_ratio_bg_total->cd();
423  nom_ratio->Draw("hist");
424  leg_entries[0] = new TLegendEntry(nom_sig, "Nominal", "l");
425 
426  // Draw shifts ------------
427  int it = 0;
428  for (std::map<int, NusSystShiftSigma*>::const_iterator sigma = fSigmas.begin(); sigma != fSigmas.end(); ++sigma, ++it) {
429  if (sigma->second->Sigma == 0)
430  continue;
431  // Total ------------
432  // signal
433  TH1* sigma_sig = (TH1*)sigma->second->ShiftedSpectrum.first->Clone();
434  for (size_t bin = 1; bin <= size_t(sigma_sig->GetNbinsX()); ++bin)
435  sigma_sig->SetBinContent(bin, sigma_sig->GetBinContent(bin)/sigma_sig->GetBinWidth(bin));
436  sigma_sig->GetYaxis()->SetTitle("Events / GeV");
437  sigma_sig->SetLineStyle(kDashed);
438  sigma_sig->SetLineColor(abs(sigma->second->Sigma)+1);
439  sigma_sig->SetLineWidth(2);
440  pad_spec_sig_total->cd();
441  sigma_sig->DrawClone("hist same");
442  TH1* sigma_ratio_sig = (TH1*)sigma_sig->Clone("sigma_ratio_sig");
443  sigma_ratio_sig->Divide(nom_sig);
444  pad_ratio_sig_total->cd();
445  sigma_ratio_sig->DrawClone("hist same");
446  // background
447  TH1* sigma_bg = (TH1*)sigma->second->ShiftedSpectrum.second->Clone();
448  for (size_t bin = 1; bin <= size_t(sigma_bg->GetNbinsX()); ++bin)
449  sigma_bg->SetBinContent(bin, sigma_bg->GetBinContent(bin)/sigma_bg->GetBinWidth(bin));
450  sigma_bg->GetYaxis()->SetTitle("Events / GeV");
451  sigma_bg->SetLineStyle(kDashed);
452  sigma_bg->SetLineColor(abs(sigma->second->Sigma)+1);
453  sigma_bg->SetLineWidth(2);
454  pad_spec_bg_total->cd();
455  sigma_bg->DrawClone("hist same");
456  TH1* sigma_ratio_bg = (TH1*)sigma_bg->Clone("sigma_ratio_bg");
457  sigma_ratio_bg->Divide(nom_bg);
458  pad_ratio_bg_total->cd();
459  sigma_ratio_bg->DrawClone("hist same");
460  leg_entries[abs(sigma->second->Sigma)]
461  = new TLegendEntry(sigma_sig->Clone(), Form("+/-%d#sigma", sigma->second->Sigma), "l");
462  // Individual ------------
463  // signal
464  pad_spec_sig->cd();
465  pad_spec_sig->Clear();
466  nom_sig->Draw("hist");
467  sigma_sig->SetLineColor(kRed);
468  sigma_sig->Draw("hist same");
469  pad_ratio_sig->cd();
470  pad_ratio_sig->Clear();
471  nom_ratio->Draw("hist");
472  sigma_ratio_sig->SetLineColor(kRed);
473  sigma_ratio_sig->Draw("hist same");
474  canv->cd();
475  canv->Clear();
476  pad_spec_sig->DrawClone();
477  pad_ratio_sig->DrawClone();
478  outFile->cd(dir_sig_shift_sigmas.c_str());
479  canv->Write(Form("%s%dSigma", fName.c_str(), sigma->second->Sigma));
480  // background
481  pad_spec_bg->cd();
482  pad_spec_bg->Clear();
483  nom_bg->Draw("hist");
484  sigma_bg->SetLineColor(kRed);
485  sigma_bg->Draw("hist same");
486  pad_ratio_bg->cd();
487  pad_ratio_bg->Clear();
488  nom_ratio->Draw("hist");
489  sigma_ratio_bg->SetLineColor(kRed);
490  sigma_ratio_bg->Draw("hist same");
491  canv->cd();
492  canv->Clear();
493  pad_spec_bg->DrawClone();
494  pad_ratio_bg->DrawClone();
495  outFile->cd(dir_bg_shift_sigmas.c_str());
496  canv->Write(Form("%s%dSigma", fName.c_str(), sigma->second->Sigma));
497 
498  }
499 
500  // Draw legend ------------
501  for (const auto& leg_entry : leg_entries)
502  leg->AddEntry(leg_entry.second->GetObject(),
503  leg_entry.second->GetLabel(),
504  leg_entry.second->GetOption());
505  pad_spec_sig_total->cd();
506  leg->Draw();
507  pad_spec_bg_total->cd();
508  leg->Draw();
509 
510  // Draw all together ------------
511  // signal
512  canv->cd();
513  canv->Clear();
514  pad_spec_sig_total->Draw();
515  pad_ratio_sig_total->Draw();
516  outFile->cd(dir_sig_shift.c_str());
517  canv->Write(Form("%sSigAllSigmas", fName.c_str()));
518  // background
519  canv->cd();
520  canv->Clear();
521  pad_spec_bg_total->Draw();
522  pad_ratio_bg_total->Draw();
523  outFile->cd(dir_bg_shift.c_str());
524  canv->Write(Form("%sBgAllSigmas", fName.c_str()));
525 
526  // clean up
527  // need to figure out what needs to be deleted
528  delete canv;
529 
530  } // function NusSystShift::DrawShift
531 
532  // --------------------------------------------------------------------------
533  std::string NusSystShift::WriteTable(std::pair<TH1*, TH1*> nominal,
534  std::string name) {
535 
536  std::ostringstream table;
537 
538  table << "\\newcommand{\\" << name << "}\n";
539  table << "{\\centerline{\n";
540  table << "\\begin{tabular}{ l | r r r | r r }\n";
541  table << " & Nom.&(+)&(-)&\\%(+)&\\%(-) \\\\ \\hline\n";
542 
543  // We write +1 and -1 sigmas to a table. First check they both exist!
544 
545  std::vector<int> sigmas = this->GetSigmas();
546  if (std::find(sigmas.begin(), sigmas.end(), 1) == sigmas.end())
547  assert(false and "+1 shift required to call WriteTable()!");
548  bool use_down = std::find(sigmas.begin(), sigmas.end(), -1) != sigmas.end();
549 
550  // Get nominal event counts for signal and background
551 
552  double sig_nom = 0, bkg_nom = 0;
553  size_t nbins = nominal.first->GetNbinsX();
554  for (size_t it = 1; it <= nbins; ++it) {
555  sig_nom += nominal.first->GetBinContent(it);
556  bkg_nom += nominal.second->GetBinContent(it);
557  }
558 
559  // Now calculate the up and down totals
560 
561  double sig_up = sig_nom, bkg_up = bkg_nom, sig_down = 0, bkg_down = 0;
562  if (use_down) {
563  sig_down = sig_nom;
564  bkg_down = bkg_nom;
565  }
566 
567  for (size_t it = 1; it <= nbins; ++it) {
568  sig_up += this->GetShiftedSpectrum(1).first->GetBinContent(it)
569  - nominal.first->GetBinContent(it);
570  sig_down += this->GetShiftedSpectrum(-1).first->GetBinContent(it)
571  - nominal.first->GetBinContent(it);
572  bkg_up += this->GetShiftedSpectrum(1).second->GetBinContent(it)
573  - nominal.second->GetBinContent(it);
574  if (use_down) {
575  sig_down += this->GetShiftedSpectrum(-1).first->GetBinContent(it)
576  - nominal.first->GetBinContent(it);
577  bkg_down += this->GetShiftedSpectrum(-1).second->GetBinContent(it)
578  - nominal.second->GetBinContent(it);
579  }
580  }
581 
582  // Write them to the table and return it
583 
584  table << "Signal (NC) & " << sig_nom << " & " << sig_up << " & ";
585  if (use_down) table << sig_down;
586  else table << "N/A";
587  table << " & " << 100*((sig_up/sig_nom)-1) << " & ";
588  if (use_down) table << 100*((sig_down/sig_nom)-1);
589  else table << "N/A";
590  table << " \\\\\n";
591  table << "Background (CC) & " << bkg_nom << " & " << bkg_up << " & ";
592  if (use_down) table << bkg_down;
593  else table << "N/A";
594  table << " & " << 100*((bkg_up/bkg_nom)-1) << " & ";
595  if (use_down) table << 100*((bkg_down/bkg_nom)-1);
596  else table << "N/A";
597  table << " \\\\\n\\end{tabular}\n}}\n";
598 
599  return table.str();
600 
601  } // function NusSystShift::WriteTable
602 
603  // -------------------------------------------------------------------------------------
604  std::pair<TH1*, TH1*> NusSystShift::GetShiftedSpectrum(int sigma) {
605  if (!fSigmas.count(sigma)) {
606  std::cout << "Error: attempting to retrieve " << sigma << " sigma shift from "
607  << fName << "; does not exist. Aborting..." << std::endl;
608  abort();
609  }
610  return fSigmas[sigma]->ShiftedSpectrum;
611  }
612 
613  // -------------------------------------------------------------------------------------
614  void NusSystShift::ProcessShift(std::pair<TH1*, TH1*>& nominal_spectrum,
615  double pot) {
616  for (auto& sigma : fSigmas) {
617  TH1* hShiftSig = GetSpectrum(sigma.second->Prediction, true, pot);
618  TH1* hShiftBg = GetSpectrum(sigma.second->Prediction, false, pot);
619  TH1* hShiftRatioSig = (TH1*)hShiftSig->Clone("shiftratsig");
620  hShiftRatioSig->Divide(nominal_spectrum.first);
621  TH1* hShiftRatioBg = (TH1*)hShiftBg->Clone("shiftratbg");
622  hShiftRatioBg->Divide(nominal_spectrum.second);
623  double shiftSig = std::abs(hShiftSig->Integral() - nominal_spectrum.first->Integral()) /
624  nominal_spectrum.first->Integral();
625  double shiftBg = std::abs(hShiftBg->Integral() - nominal_spectrum.second->Integral()) /
626  nominal_spectrum.second->Integral();
627  sigma.second->ShiftedSpectrum = std::make_pair(hShiftSig, hShiftBg);
628  sigma.second->ShiftedRatio = std::make_pair(hShiftRatioSig, hShiftRatioBg);
629  sigma.second->Shift = std::make_pair(shiftSig, shiftBg);
630  }
631  return;
632  }
633 
634  // -------------------------------------------------------------------------------------
635  /// Base class implementation of MakePredictions, which just throws an error
637  HistAxis* nus_axis, HistAxis* numu_axis,
638  Cut* fd_cut, Cut* nd_cut, Cut* numu_cut,
639  const Var* weight, Loaders* loaders) {
640 
641  // This function has to be defined in the base class for <reasons> but the
642  // user should never actually call it so it just throws an error
643 
644  assert(false and "Base class version of MakePredictions should never be called!");
645 
646  } // function NusSystShift::MakePredictions
647 
648  // -------------------------------------------------------------------------------------
649  /// SaveTo implementation for NusSystShift
650  void NusSystShift::SaveTo(TDirectory* dir) const {
651 
652  TDirectory* tmp = gDirectory;
653 
654  dir->cd();
655  TObjString("NusSystShift").Write("type");
656 
657  TObjString(fName.c_str()).Write("name");
658  TObjString(fLabel.c_str()).Write("label");
659 
660  const char* sigma_dir = "Sigmas";
661  dir->mkdir(sigma_dir);
662  for (const auto& sigma : fSigmas) {
663  const char* sigma_subdir = Form("%s/Sigma%d", sigma_dir, sigma.second->Sigma);
664  dir->mkdir(sigma_subdir);
665  sigma.second->SaveTo(dir->GetDirectory(sigma_subdir));
666  }
667 
668  tmp->cd();
669 
670  } // function NusSystShift::SaveTo
671 
672  // -------------------------------------------------------------------------------------
673  /// LoadFrom implementation for NusSystShift
674  std::unique_ptr<NusSystShift> NusSystShift::LoadFrom(TDirectory* dir) {
675 
676  DontAddDirectory guard;
677 
678  TObjString* tag = (TObjString*)dir->Get("type");
679  assert(tag);
680  assert(tag->GetString() == "NusSystShift");
681 
682  TObjString* name = (TObjString*)dir->Get("name");;
683  TObjString* label = (TObjString*)dir->Get("label");;
684  std::unique_ptr<NusSystShift> nusSystShift
685  = std::make_unique<NusSystShift>(name->GetString().Data(), label->GetString().Data());
686 
687  std::map<int, NusSystShiftSigma*> sigmas;
688  TIter next(dir->GetDirectory("Sigmas")->GetListOfKeys());
689  TKey* key;
690  while ((key = (TKey*)next())) {
691  NusSystShiftSigma* sigma
692  = NusSystShiftSigma::LoadFrom(dir->GetDirectory(Form("Sigmas/%s", key->GetName())));
693  sigmas[sigma->Sigma] = sigma;
694  }
695  nusSystShift->fSigmas = sigmas;
696 
697  return nusSystShift;
698 
699  } // function NusSystShift::LoadFrom
700 
701  // -------------------------------------------------------------------------------------
702  /// NusSystLoaderShift constructor
704  NusSystShift(name, label) {
705  }
706 
707  // -------------------------------------------------------------------------------------
708  /// Add sigma for loader shift with the same ND and FD loader
710  int sigma) {
711 
712  this->AddSigma(loaders, loaders, sigma);
713 
714  } // function NusSystLoaderShift::AddSigma
715 
716  // -------------------------------------------------------------------------------------
717  /// Add sigma for loader shift with separate ND and FD loaders
719  Loaders* loaders_fd,
720  int sigma) {
721 
722  NusSystShiftSigma* shift_sigma = new NusSystShiftSigma(fName, sigma);
723  fSigmas[sigma] = shift_sigma;
724  fLoaders[sigma] = std::make_pair(loaders_nd, loaders_fd);
725 
726  } // function NusSystLoaderShift::AddSigma
727 
728  // -------------------------------------------------------------------------------------
729  /// Add on/off for loader shift
731 
732  this->AddOnOff(loaders, loaders);
733 
734  } // function NusSystLoaderShift::AddOnOff
735 
736  // -------------------------------------------------------------------------------------
737  /// Add on/off for loader shift with separate ND and FD loaders
739  Loaders* loaders_fd) {
740 
741  for (const int sigma : {-1,1}) {
742  NusSystShiftSigma* shift_sigma = new NusSystShiftSigma(fName, sigma);
743  fSigmas[sigma] = shift_sigma;
744  fLoaders[sigma] = std::pair<Loaders*, Loaders*>(loaders_nd, loaders_fd);
745  }
746 
747  // std::cout << "Calling CheckLoaders in AddOnOff" << std::endl;
748  // CheckLoaders();
749 
750  } // function NusSystLoaderShift::AddOnOff
751 
752  // -------------------------------------------------------------------------------------
753  /// Make shifted prediction for systematic shift derived from alternate MC loader
755  HistAxis* nus_axis, HistAxis* numu_axis,
756  Cut* fd_cut, Cut* nd_cut, Cut* numu_cut,
757  const Var* weight, Loaders* loaders) {
758 
759  // CheckLoaders();
760 
761  for (auto& sigma : fSigmas) {
762 
763  sigma.second->Prediction = GetPrediction(pred_type,
764  nus_axis, numu_axis,
765  fd_cut, nd_cut, numu_cut,
766  &kNoShift, &kNoShift,
767  weight,
768  fLoaders[sigma.second->Sigma].second,
769  fLoaders[sigma.second->Sigma].first);
770  }
771 
772  } // function NusSystLoaderShift::MakePredictions
773 
774  // -------------------------------------------------------------------------------------
775  /// NusSystWeightShift constructor
777  NusSystShift(name, label) {
778  }
779 
780  // -------------------------------------------------------------------------------------
781  /// NusSystWeightShift constructor
783  for (auto& weight : fWeights)
784  delete weight.second;
785  }
786 
787  // -------------------------------------------------------------------------------------
788  /// Add sigma for weight shift
789  void NusSystWeightShift::AddSigma(const Var* weight, int sigma) {
790 
791  NusSystShiftSigma* shift_sigma = new NusSystShiftSigma(fName, sigma);
792  fSigmas[sigma] = shift_sigma;
793  fWeights[sigma] = weight;
794 
795  } // function NusSystWeightShift::AddSigma
796 
797  // -------------------------------------------------------------------------------------
798  /// Add on/off for weight shift
800 
801  for (const auto& sigma : {-1,1}) {
802  NusSystShiftSigma* shift_sigma = new NusSystShiftSigma(fName, sigma);
803  fSigmas[sigma] = shift_sigma;
804  fWeights[sigma] = weight;
805  }
806 
807  } // function NusSystWeightShift::AddOnOff
808 
809  // -------------------------------------------------------------------------------------
810  /// Make shifted prediction for systematic shift derived from CAFAna weight
812  HistAxis* nus_axis, HistAxis* numu_axis,
813  Cut* fd_cut, Cut* nd_cut, Cut* numu_cut,
814  const Var* weight, Loaders* loaders) {
815 
816  for (auto& sigma : fSigmas)
817  sigma.second->Prediction = GetPrediction(pred_type,
818  nus_axis, numu_axis,
819  fd_cut, nd_cut, numu_cut,
820  &kNoShift, &kNoShift,
821  fWeights[sigma.second->Sigma],
822  loaders);
823 
824  } // function NusSystWeightShift::MakePrediction
825 
826  // -------------------------------------------------------------------------------------
827  /// NusSystSystShift constructor
829  NusSystShift(name, label) {
830  }
831 
832  // -------------------------------------------------------------------------------------
833  /// Add sigma for syst shift
835  NusSystShiftSigma* shift_sigma = new NusSystShiftSigma(fName, sigma);
836  fSigmas[sigma] = shift_sigma;
837  fSysts[sigma] = syst;
838  return;
839  }
840 
841  // -------------------------------------------------------------------------------------
842  /// Make shifted prediction for systematic shift derived from SystShift
844  HistAxis* nus_axis, HistAxis* numu_axis,
845  Cut* fd_cut, Cut* nd_cut, Cut* numu_cut,
846  const Var* weight, Loaders* loaders) {
847 
848  for (auto& sigma : fSigmas)
849  sigma.second->Prediction = GetPrediction(pred_type,
850  nus_axis, numu_axis,
851  fd_cut, nd_cut, numu_cut,
852  &kNoShift, fSysts[sigma.second->Sigma],
853  weight, loaders);
854 
855  } // function NusSystSystShift::MakePrediction
856 
857  // -------------------------------------------------------------------------------------
858  /// NusSystMaker constructor
860  fName(name),
861  fLabel(label),
862  fSampleName(sample_name),
863  fNusAxis(nullptr),
864  fNuMuAxis(nullptr),
865  fNDCut(nullptr),
866  fFDCut(nullptr),
867  fNuMuCut(nullptr),
868  fLoaders(nullptr),
869  fWeight(nullptr),
870  fNominalPred(nullptr) {
871 
872  } // NusSystMaker constructor
873 
874  // -------------------------------------------------------------------------------------
875  /// Get the shift for a given name
877  if (!fShifts.count(shift_name)) {
878  std::cout << "Error: requested shift " << shift_name << " does not exist. Aborting..." << std::endl;
879  abort();
880  }
881  return fShifts[shift_name];
882 
883  } // function NusSystMaker::GetShift
884 
885  // -------------------------------------------------------------------------------------
886  /// Get the total shifts for a given sigma (sig+bkg)
887  std::pair<double, double> NusSystMaker::GetTotalShift(int sigma) {
888  if (!fTotalShifts.count(sigma)) {
889  std::cout << "Error: sigma " << sigma << " does not exist. Aborting..." << std::endl;
890  abort();
891  }
892  return fTotalShifts[sigma];
893 
894  } // function NusSystShift::GetTotalShift
895 
896  // -------------------------------------------------------------------------------------
897  /// Get the total shift for a given sigma and channel selection (sig/bkg)
898  double NusSystMaker::GetTotalShift(int sigma, bool signal) {
899 
900  std::pair<double, double> totalShift = this->GetTotalShift(sigma);
901  return signal ? totalShift.first : totalShift.second;
902 
903  } // function NusSystShift::GetTotalShift
904 
905  // -------------------------------------------------------------------------------------
906  /// Get the sigmas for a given systematic shift
907  std::vector<int> NusSystMaker::GetSigmas() {
908 
909  std::map<int, std::vector<std::string> > shifts_sigmas;
910  for (const auto& shift : fShifts) {
911  std::vector<int> shift_sigmas = shift.second->GetSigmas();
912  for (const auto& shift_sigma : shift_sigmas)
913  shifts_sigmas[shift_sigma].push_back(shift.second->GetName());
914  }
915  std::vector<int> sigmas;
916  for (const auto& shift_sigmas : shifts_sigmas)
917  if (shift_sigmas.second.size() == fShifts.size())
918  sigmas.push_back(shift_sigmas.first);
919 
920  return sigmas;
921 
922  } // function NusSystMaker::GetSigmas
923 
924  // -------------------------------------------------------------------------------------
925  /// Print the total shifts for each shift owned by this syst maker
926  void NusSystMaker::PrintSyst(std::ostream& os) {
927 
928  os << std::endl << std::endl << std::showpos
929  << "Systematic " << fLabel << " has total shift:" << std::endl
930  << " -1: " << fTotalShifts[-1].first << " (signal), " << fTotalShifts[-1].second << " (background)" << std::endl
931  << " +1: " << fTotalShifts[+1].first << " (signal), " << fTotalShifts[+1].second << " (background)" << std::endl;
932  os << std::endl << "Individual shifts:" << std::endl;
933  for (const auto& shift : fShifts)
934  shift.second->PrintShift(os);
935 
936  } // function NusSystMaker::PrintSyst
937 
938  // -------------------------------------------------------------------------------------
939  /// Add NusSystShift to this syst maker
941  // bool good = shift->CheckShifts();
942  // if (!good) {
943  // std::cout << "Error: attempting to add shift " << shift->GetName()
944  // << " to systematic " << fName << " without required +1/-1 sigma shifts. Aborting..."
945  // << std::endl;
946  // abort();
947  // }
948  fShifts[shift->GetName()] = shift;
949 
950  } // function NusSystMaker::AddShift
951 
952  // -------------------------------------------------------------------------------------
953  /// Make plots of all the shifted spectra and ratios for this systmaker
955 
956  // Nominal
957 
958  TH1* nom_sig = (TH1*)fSpectra[0].first->Clone();
959  for (size_t bin = 1; bin <= size_t(nom_sig->GetNbinsX()); ++bin)
960  nom_sig->SetBinContent(bin, nom_sig->GetBinContent(bin)/nom_sig->GetBinWidth(bin));
961  nom_sig->GetYaxis()->SetTitle("Events / GeV");
962  nom_sig->SetLineColor(kBlack);
963  nom_sig->SetLineStyle(kSolid);
964  nom_sig->SetLineWidth(2);
965  TH1* nom_ratio = (TH1*)nom_sig->Clone();
966  nom_ratio->GetYaxis()->SetTitle("Shifted / Nominal");
967  for (size_t bin = 0; bin <= (size_t)nom_ratio->GetNbinsX(); ++bin)
968  nom_ratio->SetBinContent(bin, 1.);
969 
970  // Draw total shifts -- all sigmas
971 
972  TCanvas* canv = new TCanvas("canv_sig", "", 800, 800);
973  TPad* pad_spec = new TPad("pad_spec_sig", "", 0., 0.375, 1., 1. );
974  TPad* pad_ratio = new TPad("pad_ratio_sig", "", 0., 0., 1., 0.375);
975  std::map<int, TLegendEntry*> leg_entry_sigma;
976  pad_spec->cd();
977  nom_sig->Draw("hist");
978  pad_ratio->cd();
979  nom_ratio->Draw("hist");
980  leg_entry_sigma[0] = new TLegendEntry(nom_sig, "Nominal", "l");
981  double ymax_ratio = 0;
982 
983  // Loop over each sigma
984 
985  for (const auto& sigma : this->GetSigmas()) {
986  if (sigma == 0)
987  continue;
988 
989  // Draw the spectrum
990 
991  TH1* sigma_sig = (TH1*)fSpectra[sigma].first->Clone();
992  for (size_t bin = 1; bin <= size_t(sigma_sig->GetNbinsX()); ++bin)
993  sigma_sig->SetBinContent(bin, sigma_sig->GetBinContent(bin)/sigma_sig->GetBinWidth(bin));
994  sigma_sig->GetYaxis()->SetTitle("Events / GeV");
995  double hmax = nom_sig->GetMaximum() > sigma_sig->GetMaximum() ?
996  nom_sig->GetMaximum() : sigma_sig->GetMaximum();
997  nom_sig->SetMaximum(1.1*hmax);
998  sigma_sig->SetLineColor(abs(sigma)+1);
999  sigma_sig->SetLineStyle(kDashed);
1000  sigma_sig->SetLineWidth(2);
1001  pad_spec->cd();
1002  sigma_sig->Draw("hist same");
1003  if (!leg_entry_sigma.count(abs(sigma)))
1004  leg_entry_sigma[abs(sigma)] = new TLegendEntry(sigma_sig, Form("+/-%d #sigma", abs(sigma)), "l");
1005 
1006  // Draw the ratio
1007 
1008  TH1* ratio_sig = (TH1*)sigma_sig->Clone();
1009  ratio_sig->Divide(nom_sig);
1010  ratio_sig->SetBinContent(0, 1);
1011  for (size_t i = 1; i <= (size_t)ratio_sig->GetNbinsX(); ++i) {
1012  double val = ratio_sig->GetBinContent(i);
1013  if (val == 0) {
1014  ratio_sig->SetBinContent(i, 1);
1015  val = 1;
1016  }
1017  if (fabs(1-val) > ymax_ratio)
1018  ymax_ratio = fabs(1-val);
1019  }
1020  pad_ratio->cd();
1021  ratio_sig->Draw("hist same");
1022  }
1023  TLegend* leg_sigma = new TLegend(0.5, 0.5, 0.88, 0.88);
1024  for (const auto& leg_entry : leg_entry_sigma)
1025  leg_sigma->AddEntry(leg_entry.second->GetObject(),
1026  leg_entry.second->GetLabel(),
1027  leg_entry.second->GetOption());
1028  pad_spec->cd();
1029  leg_sigma->Draw("same");
1030 
1031  nom_ratio->SetMinimum(1-(1.1*ymax_ratio));
1032  nom_ratio->SetMaximum(1+(1.1*ymax_ratio));
1033 
1034  canv->cd();
1035  canv->Clear();
1036  pad_spec->Draw();
1037  pad_ratio->Draw();
1038 
1039  return canv;
1040 
1041  } // function NusSystMaker::DrawSig
1042 
1043  // -------------------------------------------------------------------------------------
1044  /// Make plots of all the shifted spectra and ratios for this systmaker
1046 
1047  // Nominal
1048 
1049  TH1* nom_bkg = (TH1*)fSpectra[0].second->Clone();
1050  for (size_t bin = 1; bin <= size_t(nom_bkg->GetNbinsX()); ++bin)
1051  nom_bkg->SetBinContent(bin, nom_bkg->GetBinContent(bin)/nom_bkg->GetBinWidth(bin));
1052  nom_bkg->GetYaxis()->SetTitle("Events / GeV");
1053  nom_bkg->SetLineColor(kBlack);
1054  nom_bkg->SetLineStyle(kSolid);
1055  nom_bkg->SetLineWidth(2);
1056  TH1* nom_ratio = (TH1*)nom_bkg->Clone();
1057  nom_ratio->GetYaxis()->SetTitle("Shifted / Nominal");
1058  for (size_t bin = 0; bin <= (size_t)nom_ratio->GetNbinsX(); ++bin)
1059  nom_ratio->SetBinContent(bin, 1.);
1060 
1061  // Draw total shifts -- all sigmas
1062 
1063  TCanvas* canv = new TCanvas("canv_bkg", "", 800, 800);
1064  TPad* pad_spec = new TPad("pad_spec_bkg", "", 0., 0.375, 1., 1. );
1065  TPad* pad_ratio = new TPad("pad_ratio_bkg", "", 0., 0., 1., 0.375);
1066  std::map<int, TLegendEntry*> leg_entry_sigma;
1067  pad_spec->cd();
1068  nom_bkg->Draw("hist");
1069  pad_ratio->cd();
1070  nom_ratio->Draw("hist");
1071  leg_entry_sigma[0] = new TLegendEntry(nom_bkg, "Nominal", "l");
1072  double ymax_ratio = 0;
1073 
1074  // Loop over each sigma
1075  for (const auto& sigma : this->GetSigmas()) {
1076  if (sigma == 0)
1077  continue;
1078 
1079  // Draw the spectrum
1080  TH1* sigma_bkg = (TH1*)fSpectra[sigma].second->Clone();
1081  for (size_t bin = 1; bin <= size_t(sigma_bkg->GetNbinsX()); ++bin)
1082  sigma_bkg->SetBinContent(bin, sigma_bkg->GetBinContent(bin)/sigma_bkg->GetBinWidth(bin));
1083  sigma_bkg->GetYaxis()->SetTitle("Events / GeV");
1084  double hmax = nom_bkg->GetMaximum() > sigma_bkg->GetMaximum() ?
1085  nom_bkg->GetMaximum() : sigma_bkg->GetMaximum();
1086  nom_bkg->SetMaximum(1.1*hmax);
1087  sigma_bkg->SetLineColor(abs(sigma)+1);
1088  sigma_bkg->SetLineStyle(kDashed);
1089  sigma_bkg->SetLineWidth(2);
1090  pad_spec->cd();
1091  sigma_bkg->Draw("hist same");
1092  if (!leg_entry_sigma.count(abs(sigma)))
1093  leg_entry_sigma[abs(sigma)] = new TLegendEntry(sigma_bkg, Form("+/-%d #sigma", abs(sigma)), "l");
1094 
1095  // Draw the ratio
1096  TH1* ratio_bkg = (TH1*)sigma_bkg->Clone();
1097  ratio_bkg->Divide(nom_bkg);
1098  ratio_bkg->SetBinContent(0, 1);
1099  for (size_t i = 1; i <= (size_t)ratio_bkg->GetNbinsX(); ++i) {
1100  double val = ratio_bkg->GetBinContent(i);
1101  if (val == 0) {
1102  ratio_bkg->SetBinContent(i, 1);
1103  val = 1;
1104  }
1105  if (fabs(1-val) > ymax_ratio)
1106  ymax_ratio = fabs(1-val);
1107  }
1108  pad_ratio->cd();
1109  ratio_bkg->Draw("hist same");
1110  }
1111  TLegend* leg_sigma = new TLegend(0.5, 0.5, 0.88, 0.88);
1112  for (const auto& leg_entry : leg_entry_sigma)
1113  leg_sigma->AddEntry(leg_entry.second->GetObject(),
1114  leg_entry.second->GetLabel(),
1115  leg_entry.second->GetOption());
1116  pad_spec->cd();
1117  leg_sigma->Draw("same");
1118 
1119  nom_ratio->SetMinimum(1-(1.1*ymax_ratio));
1120  nom_ratio->SetMaximum(1+(1.1*ymax_ratio));
1121 
1122  canv->cd();
1123  canv->Clear();
1124  pad_spec->Draw();
1125  pad_ratio->Draw();
1126 
1127  return canv;
1128 
1129  } // function NusSystMaker::DrawBkg
1130 
1131  // -------------------------------------------------------------------------------------
1132  /// Make plots of all the shifted spectra and ratios for this systmaker
1134 
1135  std::string dir_top = fName;
1136  std::string dir_signal = dir_top+"/Signal/";
1137  std::string dir_background = dir_top+"/Background/";
1138  std::string dir_sig_shifts = dir_signal+"/Shifts/";
1139  std::string dir_bg_shifts = dir_background+"/Shifts/";
1140  outFile->mkdir(dir_top.c_str());
1141  outFile->mkdir(dir_signal.c_str());
1142  outFile->mkdir(dir_background.c_str());
1143  outFile->mkdir(dir_sig_shifts.c_str());
1144  outFile->mkdir(dir_bg_shifts.c_str());
1145 
1146  TCanvas* canv_sig = DrawSig();
1147  TCanvas* canv_bkg = DrawBkg();
1148  outFile->cd(dir_signal.c_str());
1149  canv_sig->Write("TotalShiftSig");
1150  outFile->cd(dir_background.c_str());
1151  canv_bkg->Write("TotalShiftBg");
1152 
1153  // Draw total shifts -- +1/-1 sigma
1154 
1155  // Draw individual shifts
1156  for (const auto& shift : fShifts)
1157  shift.second->DrawShift(outFile->GetDirectory(dir_top.c_str()), fSpectra[0]);
1158 
1159  // clean up
1160  delete canv_sig;
1161  delete canv_bkg;
1162 
1163  } // function NusSystMaker::DrawSyst
1164 
1165  // --------------------------------------------------------------------------
1166  /// Write a table with up & down shifts for this systematic
1168 
1169  std::ostringstream table;
1170 
1171  table << "\\newcommand{\\" << name << "}\n";
1172  table << "{\\centerline{\n";
1173  table << "\\begin{tabular}{ l | r r r | r r }\n";
1174  table << " & Nom.&(+)&(-)&\\%(+)&\\%(-) \\\\ \\hline\n";
1175 
1176  // We write +1 and -1 sigmas to a table. First check they both exist!
1177  std::vector<int> sigmas = this->GetSigmas();
1178  if (std::find(sigmas.begin(), sigmas.end(), 1) == sigmas.end())
1179  assert(false and "+1 sigma shift required to call WriteTable()!");
1180  bool use_down = std::find(sigmas.begin(), sigmas.end(), -1) != sigmas.end();
1181 
1182  // Get nominal event counts for signal and background
1183  double sig_nom = 0, bkg_nom = 0;
1184  size_t nbins = fSpectra[0].first->GetNbinsX();
1185  for (size_t it = 1; it <= nbins; ++it) {
1186  sig_nom += fSpectra[0].first->GetBinContent(it);
1187  bkg_nom += fSpectra[0].second->GetBinContent(it);
1188  }
1189 
1190  // Now calculate the up and down totals
1191  double sig_up = sig_nom, bkg_up = bkg_nom, sig_down = 0, bkg_down = 0;
1192  if (use_down) {
1193  sig_down = sig_nom;
1194  bkg_down = bkg_nom;
1195  }
1196 
1197  for (auto& shift : fShifts) {
1198  for (size_t it = 1; it <= nbins; ++it) {
1199  sig_up += shift.second->GetShiftedSpectrum(1).first->GetBinContent(it)
1200  - fSpectra[0].first->GetBinContent(it);
1201  bkg_up += shift.second->GetShiftedSpectrum(1).second->GetBinContent(it)
1202  - fSpectra[0].second->GetBinContent(it);
1203  if (use_down) {
1204  sig_down += shift.second->GetShiftedSpectrum(-1).first->GetBinContent(it)
1205  - fSpectra[0].first->GetBinContent(it);
1206  bkg_down += shift.second->GetShiftedSpectrum(-1).second->GetBinContent(it)
1207  - fSpectra[0].second->GetBinContent(it);
1208  }
1209  }
1210  }
1211 
1212  // Write them to the table and return it
1213  table << "Signal (NC) & " << sig_nom << " & " << sig_up << " & ";
1214  if (use_down) table << sig_down;
1215  else table << "N/A";
1216  table << " & " << 100*((sig_up/sig_nom)-1) << " & ";
1217  if (use_down) table << 100*((sig_down/sig_nom)-1);
1218  else table << "N/A";
1219  table << " \\\\\n";
1220  table << "Background (CC) & " << bkg_nom << " & " << bkg_up << " & ";
1221  if (use_down) table << bkg_down;
1222  else table << "N/A";
1223  table << " & " << 100*((bkg_up/bkg_nom)-1) << " & ";
1224  if (use_down) table << 100*((bkg_down/bkg_nom)-1);
1225  else table << "N/A";
1226  table << " \\\\\n\\end{tabular}\n}}\n";
1227 
1228  return table.str();
1229 
1230  } // function NusSystMaker::WriteTable
1231 
1232  // -------------------------------------------------------------------------------------
1233  /// Make predictions for all the NusSystShift objects owned by this syst maker
1235 
1236  fNominalPred = GetPrediction(pred_type,
1239  &kNoShift, &kNoShift,
1240  fWeight, fLoaders);
1241 
1242  for (const auto& shift : fShifts) {
1243  shift.second->MakePredictions(pred_type,
1246  fWeight, fLoaders);
1247  }
1248  } // function NusSystMaker::MakePredictions
1249 
1250  // -------------------------------------------------------------------------------------
1252 
1253  // Get nominal spectrum
1255  GetSpectrum(fNominalPred, false, pot));
1256 
1257  // Process all shifts
1258  for (const auto& shift : fShifts) {
1259  shift.second->ProcessShift(fSpectra[0], pot);
1260  NusSystShift* systShift = static_cast<NusSystShift*>(shift.second);
1261  fShifts[shift.first] = systShift;
1262  }
1263 
1264  // Get total +1/-1 sigma systematic shift -- required
1265  // Also get total sigma shift for all sigmas in this syst
1266  std::vector<int> sigmas = this->GetSigmas();
1267  Int_t nbins = fSpectra[0].first->GetNbinsX();
1268  for (const auto& sigma : sigmas) {
1269  std::vector<double> evtDiffSig(nbins, 0.), evtDiffBg(nbins, 0.);
1270  TH1* shift_sig = (TH1*)fSpectra[0].first->Clone("shift_sig");
1271  TH1* shift_bg = (TH1*)fSpectra[0].second->Clone("shift_bg");
1272  double tot_shift_sig = 0., tot_shift_bg = 0.;
1273  for (Int_t bin = 1; bin <= nbins; ++bin) {
1274  double val_sig = 0., val_bg = 0.;
1275  // If there's more than one shift, we should add them in quadrature
1276  if (fShifts.size() > 1) {
1277  for (auto& shift : fShifts) {
1278  double bin_shift_sig = shift.second->GetShiftedSpectrum(sigma).first->GetBinContent(bin) -
1279  fSpectra[0].first->GetBinContent(bin);
1280  double bin_shift_bg = shift.second->GetShiftedSpectrum(sigma).second->GetBinContent(bin) -
1281  fSpectra[0].second->GetBinContent(bin);
1282  val_sig += pow(bin_shift_sig,2);
1283  val_bg += pow(bin_shift_bg,2);
1284  }
1285  double sign = sigma > 0 ? 1.0 : -1.0;
1286  shift_sig->SetBinContent(bin, fSpectra[0].first->GetBinContent(bin) + (sign * sqrt(val_sig)));
1287  shift_bg->SetBinContent(bin, fSpectra[0].second->GetBinContent(bin) + (sign * sqrt(val_bg)));
1288  }
1289  // If there's only one shift, just use the raw value
1290  else {
1291  val_sig = fShifts.begin()->second->GetShiftedSpectrum(sigma).first->GetBinContent(bin)
1292  - fSpectra[0].first->GetBinContent(bin);
1293  val_bg = fShifts.begin()->second->GetShiftedSpectrum(sigma).second->GetBinContent(bin)
1294  - fSpectra[0].second->GetBinContent(bin);
1295  shift_sig->SetBinContent(bin, fShifts.begin()->second->GetShiftedSpectrum(sigma).first->GetBinContent(bin));
1296  shift_bg->SetBinContent(bin, fShifts.begin()->second->GetShiftedSpectrum(sigma).second->GetBinContent(bin));
1297  }
1298  tot_shift_sig += val_sig;
1299  tot_shift_bg += val_bg;
1300  }
1301  fTotalShifts[sigma] = std::make_pair(tot_shift_sig, tot_shift_bg);
1302  fSpectra[sigma] = std::make_pair(shift_sig, shift_bg);
1303  }
1304 
1305  } // function NusSystMaker::ProcessSyst
1306 
1307  // -------------------------------------------------------------------------------------
1309  fNusAxis = axis;
1310  fNuMuAxis = axis;
1311  }
1312 
1313  // -------------------------------------------------------------------------------------
1314  void NusSystMaker::SetCuts(Cut* fd_cut, Cut* nd_cut, Cut* numu_cut) {
1315  fNDCut = nd_cut;
1316  fFDCut = fd_cut;
1317  fNuMuCut = numu_cut;
1318  }
1319 
1320  // -------------------------------------------------------------------------------------
1322  fLoaders = loaders;
1323  }
1324 
1325  // -------------------------------------------------------------------------------------
1327  fWeight = weight;
1328  }
1329 
1330  // -------------------------------------------------------------------------------------
1331  /// Function to make ISyst for this NusSystMaker
1333 
1334  if (fSpectra.size() == 0) {
1335  std::cout << "There are no spectra owned by this NusSystMaker!"
1336  << " Did you forget to call Go()?" << std::endl;
1337  return nullptr;
1338  }
1339 
1340  // Make ratios from shifted spectra
1341  std::map<int, std::pair<TH1*, TH1*>> shifts;
1342  for (const auto& sigma : fSpectra) {
1343  TH1* sig_ratio = (TH1*)sigma.second.first->Clone(Form("SigRatioSigma%d", sigma.first));
1344  TH1* bg_ratio = (TH1*)sigma.second.second->Clone(Form("BgRatioSigma%d", sigma.first));
1345  sig_ratio->Divide(fSpectra.at(0).first);
1346  bg_ratio->Divide(fSpectra.at(0).second);
1347  shifts[sigma.first] = std::make_pair(sig_ratio, bg_ratio);
1348  }
1349 
1350  // Return new ISyst
1351  if (fSampleName.empty())
1352  return new NusISyst(fName, fLabel, fSampleName, shifts);
1353  else
1354  return new NusISyst(fSampleName + "_" + fName, fLabel, fSampleName, shifts);
1355 
1356 
1357  } // function NusSystMaker::MakeISyst
1358 
1359  // -------------------------------------------------------------------------------------
1360  /// SaveTo implementation for NusSystMaker
1361  void NusSystMaker::SaveTo(TDirectory* dir) const {
1362 
1363  TDirectory* tmp = gDirectory;
1364 
1365  dir->cd();
1366  TObjString("NusSystMaker").Write("type");
1367 
1368  TObjString(fName.c_str()).Write("name");
1369  TObjString(fLabel.c_str()).Write("label");
1370  TObjString(fSampleName.c_str()).Write("sample_name");
1371 
1372  fNominalPred->SaveTo(dir->mkdir("nompred"));
1373 
1374  const char* shifts_dir = "Shifts";
1375  dir->mkdir("Shifts");
1376  for (const auto& shift : fShifts) {
1377  const char* shift_subdir = Form("%s/%s", shifts_dir, shift.first.c_str());
1378  dir->mkdir(shift_subdir);
1379  shift.second->SaveTo(dir->GetDirectory(shift_subdir));
1380  }
1381 
1382  tmp->cd();
1383 
1384  } // function NusSystMaker::SaveTo
1385 
1386  // -------------------------------------------------------------------------------------
1387  /// LoadFrom implementation for NusSystMaker
1388  std::unique_ptr<NusSystMaker> NusSystMaker::LoadFrom(TDirectory* dir) {
1389 
1390  DontAddDirectory guard;
1391 
1392  TObjString* tag = (TObjString*)dir->Get("type");
1393  assert(tag);
1394  assert(tag->GetString() == "NusSystMaker");
1395 
1396  TObjString* name = (TObjString*)dir->Get("name");
1397  assert(name);
1398 
1399  TObjString* label = (TObjString*)dir->Get("label");
1400  assert(label);
1401 
1402  TObjString* sample_name = (TObjString*)dir->Get("sample_name");
1403  assert(sample_name);
1404 
1405  std::unique_ptr<NusSystMaker> nusSystMaker
1406  = std::make_unique<NusSystMaker>(name->GetString().Data(),
1407  label->GetString().Data(), sample_name->GetString().Data());
1408 
1409  // Load prediction
1410 
1411  nusSystMaker->fNominalPred = LoadFromFile<IPrediction>(static_cast<TFile*>(dir), "nompred").release();
1412 
1413  // Load shift makers
1414 
1415  std::map<std::string, NusSystShift*> shifts;
1416 
1417  TIter shift_next(dir->GetDirectory("Shifts")->GetListOfKeys());
1418  TKey* shift_key;
1419  while ((shift_key = (TKey*)shift_next())) {
1421  = NusSystShift::LoadFrom(dir->GetDirectory(Form("Shifts/%s", shift_key->GetName()))).release();
1422  shifts[shift->GetName()] = shift;
1423  }
1424  nusSystMaker->fShifts = shifts;
1425 
1426  return std::move(nusSystMaker);
1427 
1428  } // function NusSystMaker::LoadFrom
1429 
1430  // -------------------------------------------------------------------------------------
1431  /// NusSystematicsMaker constructor
1433  fName(name),
1434  fPredType(pred_type),
1435  fNusAxis(nullptr),
1436  fNuMuAxis(nullptr),
1437  fNDCut(nullptr),
1438  fFDCut(nullptr),
1439  fNuMuCut(nullptr),
1440  fLoaders(nullptr),
1441  fLightLoaders(nullptr),
1442  fWeight(nullptr)
1443  {
1444  if (std::find(kSupportedPredTypes.begin(), kSupportedPredTypes.end(), pred_type) == kSupportedPredTypes.end()) {
1445  std::cout << "NusSystematicsMaker error: Unsupported prediction type. Aborting..." << std::endl;
1446  abort();
1447  }
1448  } // NusSystematicsMaker constructor
1449 
1450  // -------------------------------------------------------------------------------------
1451  /// Get the names of all attached systematics
1452  std::vector<std::string> NusSystematicsMaker::GetSystematicsNames() {
1453 
1454  std::vector<std::string> systNames;
1455  for (const auto& syst : fSystMakers)
1456  systNames.push_back(syst.first);
1457  return systNames;
1458 
1459  } // function NusSystematicsMaker::GetSystematicsNames
1460 
1461  // -------------------------------------------------------------------------------------
1462  /// Make ISyst objects for each attached systematic and return them
1463  std::vector<NusISyst*> NusSystematicsMaker::MakeISysts() {
1464 
1465  std::vector<NusISyst*> systs;
1466  for (const auto& maker : fSystMakers)
1467  systs.push_back(maker.second->MakeISyst());
1468 
1469  return systs;
1470 
1471  } // function NusSystematicsMaker::MakeISysts
1472 
1473  // -------------------------------------------------------------------------------------
1475 
1476  syst->SetAxes(fNusAxis, fNuMuAxis);
1477  syst->SetCuts(fFDCut, fNDCut, fNuMuCut);
1478  syst->SetNominalWeight(fWeight);
1479  if (syst->GetName().find("LightLevel") != std::string::npos)
1481  else
1482  syst->SetNominalLoaders(fLoaders);
1483  syst->MakePredictions(fPredType);
1484  fSystMakers[syst->GetName()] = syst;
1485  }
1486 
1487  // -------------------------------------------------------------------------------------
1488  void NusSystematicsMaker::PrintSysts(std::ostream& os) {
1489  for (const auto& syst : fSystMakers)
1490  syst.second->PrintSyst(os);
1491  return;
1492  }
1493 
1494  // -------------------------------------------------------------------------------------
1496  for (const auto& syst : fSystMakers)
1497  syst.second->DrawSyst(outFile);
1498  return;
1499  }
1500 
1501  // -------------------------------------------------------------------------------------
1503  for (const auto& syst : fSystMakers) {
1504  syst.second->ProcessSyst(pot);
1505  fSystMakers[syst.first] = syst.second;
1506  }
1507  }
1508 
1509  // -------------------------------------------------------------------------------------
1510  void NusSystematicsMaker::SetAxes(HistAxis* nus_axis, HistAxis* numu_axis) {
1511  fNusAxis = nus_axis;
1512  fNuMuAxis = numu_axis;
1513  }
1514 
1515  // -------------------------------------------------------------------------------------
1516  void NusSystematicsMaker::SetCuts(Cut* fd_cut, Cut* nd_cut, Cut* numu_cut) {
1517  fNDCut = nd_cut;
1518  fFDCut = fd_cut;
1519  fNuMuCut = numu_cut;
1520  }
1521 
1522  // -------------------------------------------------------------------------------------
1524  fLoaders = loaders;
1525  fLightLoaders = loaders_light;
1526  }
1527 
1528  // -------------------------------------------------------------------------------------
1530  fWeight = weight;
1531  }
1532 
1533  // -------------------------------------------------------------------------------------
1534  /// function to add another NusSystematicsMaker to this one
1536 
1537  std::map<std::string, NusSystMaker*> syst_makers;
1538  for (const auto& syst_maker : syst_makers) {
1539  if (fSystMakers.count(syst_maker.first) and !overwrite)
1540  continue;
1541  fSystMakers[syst_maker.first] = syst_maker.second;
1542  }
1543  } // function NusSystematicsMaker::Add
1544 
1545  // -------------------------------------------------------------------------------------
1546  void NusSystematicsMaker::SaveTo(TDirectory* dir, bool separate) const {
1547 
1548  TDirectory* tmp = gDirectory;
1549 
1550  dir->cd();
1551  TObjString("NusSystematicsMaker").Write("type");
1552 
1553  TObjString(fName.c_str()).Write("name");
1554 
1555  int pred_type = -1;
1556  switch (fPredType) {
1558  pred_type = 0;
1559  break;
1560  case NusSystsPredType::kFD:
1561  pred_type = 1;
1562  break;
1563  case NusSystsPredType::kND:
1564  pred_type = 2;
1565  break;
1566  default:
1567  break;
1568  }
1569  TVectorD predtype(1);
1570  predtype[0] = pred_type;
1571  predtype.Write("predtype");
1572 
1573  TObjArray syst_names;
1574  dir->mkdir("SystMakers");
1575  for (const auto& syst_maker : fSystMakers) {
1576  // std::cout << "Saving SystMaker " << syst_maker.second->GetName() << std::endl;
1577  syst_names.Add(new TObjString(syst_maker.first.c_str()));
1578  const char* syst_maker_subdir = Form("SystMakers/%s", syst_maker.first.c_str());
1579  TDirectory* syst_maker_dir;
1580  TFile* syst_maker_file = nullptr;
1581  if (separate) {
1582  syst_maker_file = new TFile(Form("%s_%s.root", dir->GetName(), syst_maker.first.c_str()), "RECREATE");
1583  syst_maker_file->mkdir(dir->GetName());
1584  syst_maker_file->mkdir(Form("%s/%s", dir->GetName(), "SystMakers"));
1585  syst_maker_file->mkdir(Form("%s/%s", dir->GetName(), syst_maker_subdir));
1586  syst_maker_dir = syst_maker_file->GetDirectory(Form("%s/%s", dir->GetName(), syst_maker_subdir));
1587  }
1588  else {
1589  dir->mkdir(syst_maker_subdir);
1590  syst_maker_dir = dir->GetDirectory(syst_maker_subdir);
1591  }
1592  syst_maker.second->SaveTo(syst_maker_dir);
1593  if (syst_maker_file) delete syst_maker_file;
1594  }
1595  dir->cd();
1596  syst_names.Write("systnames", TObject::kSingleKey);
1597 
1598  tmp->cd();
1599 
1600  } // function NusSystematicsMaker::SaveTo
1601 
1602  //----------------------------------------------------------------------
1603  std::unique_ptr<NusSystematicsMaker> NusSystematicsMaker::LoadFrom(
1604  TDirectory* dir, std::string path) {
1605 
1606  DontAddDirectory guard;
1607 
1608  TObjString* tag = (TObjString*)dir->Get("type");
1609  assert(tag);
1610  assert(tag->GetString() == "NusSystematicsMaker");
1611 
1612  TObjString* ts_name = (TObjString*)dir->Get("name");
1613  assert(ts_name);
1614  std::string name(ts_name->GetString().Data());
1615 
1616  const TVectorD predtype = *(TVectorD*)dir->Get("predtype");
1618  switch (static_cast<int>(std::round(predtype[0]))) {
1619  case 0:
1620  pred_type = NusSystsPredType::kExtrap;
1621  break;
1622  case 1:
1623  pred_type = NusSystsPredType::kFD;
1624  break;
1625  case 2:
1626  pred_type = NusSystsPredType::kND;
1627  break;
1628  }
1629  TObjArray* systnames = (TObjArray*)dir->Get("systnames");
1630 
1631  std::unique_ptr<NusSystematicsMaker> nus_systs_maker = std::make_unique<NusSystematicsMaker>(name, pred_type);
1632 
1633  std::map<std::string, NusSystMaker*> syst_makers;
1634  TIter next(dir->GetDirectory("SystMakers")->GetListOfKeys());
1635  TKey* key;
1636  while ((key = (TKey*)next())) {
1637  NusSystMaker* syst_maker
1638  = NusSystMaker::LoadFrom(dir->GetDirectory(Form("SystMakers/%s", key->GetName()))).release();
1639  syst_makers[syst_maker->GetName()] = syst_maker;
1640  }
1641 
1642  if (syst_makers.size() == 0) {
1643  for (int syst = 0; syst < systnames->GetEntries(); ++syst) {
1644  const char* systname = ((TObjString*)systnames->At(syst))->GetString().Data();
1645  TFile* systfile = new TFile(Form("%s/%s_%s.root", path.c_str(), dir->GetName(), systname), "READ");
1646  NusSystMaker* syst_maker
1647  = NusSystMaker::LoadFrom(systfile->GetDirectory(Form("%s/SystMakers/%s", dir->GetName(), systname))).release();
1648  syst_makers[syst_maker->GetName()] = syst_maker;
1649  delete systfile;
1650  }
1651  }
1652 
1653  nus_systs_maker->fSystMakers = syst_makers;
1654 
1655  return std::move(nus_systs_maker);
1656 
1657  } // function NusSystematicsMaker::LoadFrom
1658 
1659  // -------------------------------------------------------------------------------------
1661  std::map<int, std::pair<TH1*, TH1*>> shifts) :
1662  ISyst(name, label),
1663  fSampleName(sample_name),
1664  fShifts(shifts)
1665  {}
1666 
1667  // -------------------------------------------------------------------------------------
1669 
1670  // Clean up the weight histograms
1671 
1672  for (auto shift : fShifts) {
1673  if (shift.second.first) delete shift.second.first;
1674  if (shift.second.second) delete shift.second.second;
1675  }
1676 
1677  } // NusISyst destructor
1678 
1679  // -------------------------------------------------------------------------------------
1680  /// Implementation of standard ISyst::Shift function for NusISyst
1681  void NusISyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const {
1682 
1683  // Find the energy estimate for this slice
1684  double calE = kNus18Energy(sr);
1685 
1686  // Find the signal or background oscillation channel
1688 
1689  weight *= WeightFor(chan, sigma, calE);
1690 
1691  } // function NusISyst::Shift
1692 
1693  // -------------------------------------------------------------------------------------
1694  /// SaveTo implementation for NusISyst
1695  void NusISyst::SaveTo(TDirectory* dir) const {
1696 
1697  TDirectory* tmp = gDirectory;
1698 
1699  // Go into directory
1700  dir->cd();
1701  TObjString("NusISyst").Write("type");
1702 
1703  // Write names
1704  TObjString(ShortName().c_str()).Write("name");
1705  TObjString(LatexName().c_str()).Write("label");
1706  TObjString(fSampleName.c_str()).Write("sample_name");
1707 
1708  TDirectory* sigmas_dir = dir->mkdir("Sigmas");
1709 
1710  // Save each shift
1711  for (auto const& sigma : fShifts) {
1712 
1713  TDirectory* sigma_dir = sigmas_dir->mkdir(Form("Sigma%d", sigma.first));
1714 
1715  // Save sigma value
1716  TVectorD v_sigma(1);
1717  v_sigma[0] = sigma.first;
1718  sigma_dir->WriteTObject(&v_sigma, "sigma");
1719 
1720  // Save shift histograms
1721  sigma_dir->WriteTObject(sigma.second.first, "sig");
1722  sigma_dir->WriteTObject(sigma.second.second, "bkg");
1723 
1724  }
1725 
1726  tmp->cd();
1727 
1728  } // function NusISyst::SaveTo
1729 
1730  // -------------------------------------------------------------------------------------
1731  /// LoadFrom implementation for NusISyst
1732  std::unique_ptr<NusISyst> NusISyst::LoadFrom(TDirectory* dir) {
1733 
1734  DontAddDirectory guard;
1735 
1736  // Check this is a NusISyst
1737  TObjString* tag = (TObjString*)dir->Get("type");
1738  assert(tag);
1739  assert(tag->GetString() == "NusISyst");
1740 
1741  // Load names
1742  TObjString* name = (TObjString*)dir->Get("name");
1743  assert(name);
1744  TObjString* label = (TObjString*)dir->Get("label");
1745  assert(label);
1746  TObjString* sample_name = (TObjString*)dir->Get("sample_name");
1747  assert(sample_name);
1748 
1749  std::map<int, std::pair<TH1*, TH1*>> shifts;
1750 
1751  // Load all shifts
1752  TIter next(dir->GetDirectory("Sigmas")->GetListOfKeys());
1753  TKey* key;
1754  while ((key = (TKey*)next())) {
1755  TDirectory* sigma_dir = (TDirectory*)key->ReadObj();
1756  const TVectorD v_sigma = *(TVectorD*)sigma_dir->Get("sigma");
1757  TH1* h_signal = (TH1*)sigma_dir->Get("sig");
1758  TH1* h_background = (TH1*)sigma_dir->Get("bkg");
1759  shifts[v_sigma[0]] = std::make_pair(h_signal, h_background);
1760  }
1761 
1762  std::unique_ptr<NusISyst> ret = std::make_unique<NusISyst>(name->GetString().Data(),
1763  label->GetString().Data(),sample_name->GetString().Data(), shifts);
1764 
1765  return std::move(ret);
1766 
1767  } // function NusISyst::LoadFrom
1768 
1769  // -------------------------------------------------------------------------------------
1770  /// Get channel (signal or background) from this event's standard record
1772 
1774  if (sr->mc.nnu == 0)
1775  channel = NusSystsChannel::kNC;
1776  else if (!sr->mc.nu[0].iscc)
1777  channel = NusSystsChannel::kNC;
1778  else if (sr->mc.nu[0].iscc)
1779  channel = NusSystsChannel::kBG;
1780  if (channel == NusSystsChannel::kUnknown) {
1781  std::cout << "Error: Unknown oscillation channel." << std::endl;
1782  abort();
1783  }
1784 
1785  return channel;
1786 
1787  } // function NusISyst::GetNusChannel
1788 
1789  // -------------------------------------------------------------------------------------
1790  /// Get weight for a given energy, channel and sigma shift
1791  double NusISyst::WeightFor(NusSystsChannel channel, double sigma, double calE) const {
1792 
1793  // If it's a one-sided shift (ie. no -1 sigma weights) then we want to return
1794  // 1 for a negative sigma
1795  if (sigma < 0 and not fShifts.count(-1)) return 1;
1796  const int bin = fShifts.at(0).first->FindBin(calE);
1797 
1798  auto low = fShifts.begin();
1799  auto end = fShifts.end();
1800  std::advance(end, -2);
1801  if (sigma < fShifts.begin()->first)
1802  low = fShifts.begin();
1803  else if (sigma >= end->first)
1804  low = end;
1805  else {
1806  for (auto it = fShifts.begin(); it != end; ++it) {
1807  if (sigma >= it->first) {
1808  low = it;
1809  break;
1810  }
1811  }
1812  }
1813 
1814  auto high = low;
1815  ++high;
1816 
1817  // // Why would we have templates differing by more than 1 sigma?
1818  // // fracpart below assumes this
1819  assert(high->first - low->first == 1);
1820 
1821  TH1* h_low = nullptr;
1822  TH1* h_high = nullptr;
1823 
1824  if (channel == NusSystsChannel::kNC) {
1825  h_low = low->second.first;
1826  h_high = high->second.first;
1827  }
1828  else {
1829  h_low = low->second.second;
1830  h_high = high->second.second;
1831  }
1832 
1833  const double fracpart = sigma - low->first;
1834  const double ret = (fracpart*h_high->GetBinContent(bin)) +
1835  ((1-fracpart)*h_low->GetBinContent(bin));
1836 
1837  return std::max(0., ret); // Keep the LL from blowing up
1838 
1839  } // function NusISyst::WeightFor
1840 
1841 } // namespace ana
virtual void MakePredictions(NusSystsPredType pred_type, HistAxis *nus_axis, HistAxis *numu_axis, Cut *fd_cut, Cut *nd_cut, Cut *numu_cut, const Var *weight, Loaders *loaders) override
Make shifted prediction for systematic shift derived from SystShift.
NusSystShift * GetShift(std::string shift_name)
Get the shift for a given name.
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
std::pair< double, double > Shift
Definition: NusSystsMaker.h:93
std::string MakeLatexCommandName(const std::string &instring)
std::string fLabel
void Add(NusSystematicsMaker *systMaker, bool overwrite)
function to add another NusSystematicsMaker to this one
TCanvas * DrawBkg(TH1 *nominal)
osc::OscCalculatorDumb calc
static std::unique_ptr< NusISyst > LoadFrom(TDirectory *dir)
LoadFrom implementation for NusISyst.
Oscillation analysis framework, runs over CAF files outside of ART.
const Var kNus18Energy([](const caf::SRProxy *sr){bool h_FHC=sr->spill.isFHC;bool h_RHC=sr->spill.isRHC;double cale=sr->slc.calE;double recoE=0.;if(h_FHC &&sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE18 *cale;else if(h_FHC &&sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE18 *cale;else if(h_RHC &&sr->hdr.det==caf::kFARDET) recoE=FDscaleCalE18RHC *cale;else if(h_RHC &&sr->hdr.det==caf::kNEARDET) recoE=NDscaleCalE18RHC *cale;return recoE;})
Definition: NusVars.h:17
set< int >::iterator it
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
std::map< int, std::pair< Loaders *, Loaders * > > fLoaders
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
~NusSystWeightShift()
NusSystWeightShift constructor.
Adapt the PMNS_Sterile calculator to standard interface.
void DrawSysts(TFile *outFile)
void DrawShift(TDirectory *outFile, std::pair< TH1 *, TH1 * > &nominal)
VectorProxy< SRNeutrinoProxy > nu
Definition: SRProxy.h:1890
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:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
std::map< int, std::pair< TH1 *, TH1 * > > fShifts
IPrediction * GetPrediction(NusSystsPredType pred_type, const HistAxis *axis, const HistAxis *numu_axis, const Cut *fd_cut, const Cut *nd_cut, const Cut *numu_cut, const SystShifts *shift_data, const SystShifts *shift_mc, const Var *weight, Loaders *loaders, Loaders *loaders_nd)
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
virtual void SaveTo(TDirectory *dir) const
TCanvas * DrawBkg()
Make plots of all the shifted spectra and ratios for this systmaker.
const Var weight
std::string fName
std::map< int, std::pair< double, double > > fTotalShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void ProcessShift(std::pair< TH1 *, TH1 * > &nominal_spectrum, double pot)
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:17
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Var * fWeight
constexpr T pow(T x)
Definition: pow.h:75
NusSystShift(std::string name, std::string label)
std::string fSampleName
static std::unique_ptr< NusSystematicsMaker > LoadFrom(TDirectory *dir, std::string path="")
TCanvas * canv
virtual void MakePredictions(NusSystsPredType pred_type, HistAxis *nus_axis, HistAxis *numu_axis, Cut *fd_cut, Cut *nd_cut, Cut *numu_cut, const Var *weight, Loaders *loaders) override
Make shifted prediction for systematic shift derived from CAFAna weight.
void AddShift(NusSystShift *shift)
Add NusSystShift to this syst maker.
Proxy for StandardRecord.
Definition: SRProxy.h:2237
Float_t tmp
Definition: plot.C:36
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:335
void AddSigma(const Var *weight, int sigma)
Add sigma for weight shift.
TCanvas * DrawSig(TH1 *nominal)
Generates extrapolated NC predictions using ProportionalDecomp.
fvar< T > round(const fvar< T > &x)
Definition: round.hpp:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
float abs(float number)
Definition: d0nt_math.hpp:39
std::string fSampleName
std::string WriteTable(std::pair< TH1 *, TH1 * > nominal, std::string name)
std::map< int, NusSystShiftSigma * > fSigmas
void SaveTo(TDirectory *dir) const
SaveTo implementation for NusISyst.
HistAxis * fNuMuAxis
IPrediction * fNominalPred
std::string WriteTable(std::string name)
Write a table with up & down shifts for this systematic.
const char * label
static std::unique_ptr< NusSystMaker > LoadFrom(TDirectory *dir)
LoadFrom implementation for NusSystMaker.
void AddSystematic(NusSystMaker *syst)
virtual void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Implementation of standard ISyst::Shift function for NusISyst.
const int nbins
Definition: cellShifts.C:15
Charged-current interactions.
Definition: IPrediction.h:39
Definition: Shift.h:6
Proxy< short int > nnu
Definition: SRProxy.h:1899
const std::vector< NusSystsPredType > kSupportedPredTypes
Definition: NusSystsMaker.h:65
TFile * outFile
Definition: PlotXSec.C:135
std::pair< TH1 *, TH1 * > GetShiftedSpectrum(int sigma)
std::string fLabel
void AddOnOff(const Var *weight)
Add on/off for weight shift.
IPrediction * Prediction
Definition: NusSystsMaker.h:90
void SetNominalWeight(const Var *weight)
void SetNominalLoaders(Loaders *loaders)
void SetCuts(Cut *fd_cut, Cut *nd_cut, Cut *numu_cut=nullptr)
void SaveTo(TDirectory *dir) const
SaveTo implementation for NusSystMaker.
#define pot
std::vector< std::string > GetSystematicsNames()
Get the names of all attached systematics.
void SetAxes(HistAxis *axis, HistAxis *numu_axis=nullptr)
const std::string path
Definition: plot_BEN.C:43
std::map< int, std::pair< double, double > > GetTotalShift()
void ProcessSyst(double pot)
std::vector< int > GetSigmas()
NusSystsPredType fPredType
std::map< std::string, NusSystShift * > fShifts
NusSystsPredType
Definition: NusSystsMaker.h:51
std::string fName
void SetNominalLoaders(Loaders *loader, Loaders *loader_light=nullptr)
void DrawSyst(TFile *outFile)
Make plots of all the shifted spectra and ratios for this systmaker.
float bin[41]
Definition: plottest35.C:14
static std::unique_ptr< NusSystShift > LoadFrom(TDirectory *dir)
LoadFrom implementation for NusSystShift.
void PrintSyst(std::ostream &os=std::cout)
Print the total shifts for each shift owned by this syst maker.
def sign(x)
Definition: canMan.py:204
const SystShifts kNoShift
Definition: SystShifts.h:112
OStream cout
Definition: OStream.cxx:6
std::map< int, std::pair< TH1 *, TH1 * > > fSpectra
void PrintSysts(std::ostream &os=std::cout)
HistAxis * fNusAxis
virtual const std::string & LatexName() const final
The name used on plots (ROOT&#39;s TLatex syntax)
Definition: ISyst.h:30
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
NusSystWeightShift(std::string name, std::string label)
NusSystWeightShift constructor.
void SaveTo(TDirectory *dir) const
NusSystsChannel GetNusChannel(caf::SRProxy *sr) const
Get channel (signal or background) from this event&#39;s standard record.
NusSystsChannel
Definition: NusSystsMaker.h:39
Template for Cut and SpillCut.
Definition: Cut.h:15
ifstream in
Definition: comparison.C:7
TCanvas * DrawSig()
Make plots of all the shifted spectra and ratios for this systmaker.
void SetNominalWeight(const Var *weight)
string GetString(xmlDocPtr xml_doc, string node_path)
std::vector< NusISyst * > MakeISysts()
Make ISyst objects for each attached systematic and return them.
TDirectory * dir
Definition: macro.C:5
NusISyst(std::string name, std::string label, std::string sample_name="", std::map< int, std::pair< TH1 *, TH1 * > > shifts={})
void AddOnOff(Loaders *loaders)
Add on/off for loader shift.
void PrintSigma(std::ostream &os=std::cout)
Print the sigmas for a given shift.
std::map< std::string, std::string > systNames
Definition: PlotUnfolding.C:32
std::vector< Loaders * > loaders
Definition: syst_header.h:385
string syst
Definition: plotSysts.py:176
NusSystMaker(std::string name, std::string label, std::string sample_name="")
NusSystMaker constructor.
TH1 * GetSpectrum(IPrediction *pred, bool signal, double pot)
void PrintShift(std::ostream &os=std::cout)
virtual void MakePredictions(NusSystsPredType pred_type, HistAxis *nus_axis, HistAxis *numu_axis, Cut *fd_cut, Cut *nd_cut, Cut *numu_cut, const Var *weight, Loaders *loaders) override
Make shifted prediction for systematic shift derived from alternate MC loader.
double WeightFor(NusSystsChannel channel, double sigma, double calE) const
Get weight for a given energy, channel and sigma shift.
std::pair< TH1 *, TH1 * > ShiftedRatio
Definition: NusSystsMaker.h:92
SRTruthBranchProxy mc
Definition: SRProxy.h:2253
Neutral-current interactions.
Definition: IPrediction.h:40
std::string GetName()
std::vector< int > GetSigmas()
Get the sigmas for a given systematic shift.
static NusSystShiftSigma * LoadFrom(TDirectory *dir)
void SaveTo(TDirectory *dir) const
SaveTo implementation for NusSystShift.
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
std::pair< TH1 *, TH1 * > ShiftedSpectrum
Definition: NusSystsMaker.h:91
void AddSigma(SystShifts *syst, int sigma)
Add sigma for syst shift.
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void AddSigma(Loaders *loaders, int sigma)
Add sigma for loader shift with the same ND and FD loader.
void SetAxes(HistAxis *axis, HistAxis *numu_axis=nullptr)
NusISyst * MakeISyst()
Function to make ISyst for this NusSystMaker.
std::unique_ptr< IPrediction > Generate(Loaders &loaders, const SystShifts &shiftMC=kNoShift) const override
void SetDm(int i, double dm)
void MakePredictions(NusSystsPredType pred_type)
Make predictions for all the NusSystShift objects owned by this syst maker.
osc::OscCalculatorSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
All neutrinos, any flavor.
Definition: IPrediction.h:26
Prevent histograms being added to the current directory.
Definition: Utilities.h:66
Generates Near Detector predictions.
std::map< int, const Var * > fWeights
void next()
Definition: show_event.C:84
virtual void MakePredictions(NusSystsPredType pred_type, HistAxis *nus_axis, HistAxis *numu_axis, Cut *fd_cut, Cut *nd_cut, Cut *numu_cut, const Var *weight, Loaders *loaders)
Base class implementation of MakePredictions, which just throws an error.
std::string GetName()
void SaveTo(TDirectory *dir, bool separate=false) const
NusSystLoaderShift(std::string name, std::string label)
NusSystLoaderShift constructor.
static constexpr Double_t sr
Definition: Munits.h:164
std::map< int, SystShifts * > fSysts
std::map< std::string, NusSystMaker * > fSystMakers
NusSystSystShift(std::string name, std::string label)
NusSystSystShift constructor.
void SetCuts(Cut *fd_cut, Cut *nd_cut, Cut *numu_cut=nullptr)
Template for Var and SpillVar.
Definition: Var.h:16
gm Write()
NusSystematicsMaker(std::string name, NusSystsPredType pred_type)
NusSystematicsMaker constructor.